Comparative Evaluation of Covalent Docking Tools

six covalent docking tools, AutoDock4, CovDock, FITTED, GOLD, ICM-Pro and ... FITTED and ICM-Pro were found to be less sensitive up to 35 heavy atoms...
4 downloads 0 Views 11MB Size
Article Cite This: J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

pubs.acs.org/jcim

Comparative Evaluation of Covalent Docking Tools Andrea Scarpino, György G. Ferenczy, and György M. Keserű * Medicinal Chemistry Research Group, Research Centre for Natural Sciences, Hungarian Academy of Sciences, Magyar tudósok körútja 2, Budapest 1117, Hungary

Downloaded via TU MUENCHEN on July 7, 2018 at 12:19:56 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Increased interest in covalent drug discovery led to the development of computer programs predicting binding mode and affinity of covalent inhibitors. Here we compare the performance of six covalent docking tools, AutoDock4, CovDock, FITTED, GOLD, ICM-Pro, and MOE, for reproducing experimental binding modes in an unprecedently large and diverse set of covalent complexes. It was found that 40−60% of the top scoring ligand poses are within 2.0 Å RMSD from the experimental binding mode. This rate showed program dependent increase and achieved 50−90% when the best RMSD among the top ten scoring poses was considered. This performance is comparable to that of noncovalent docking tools and therefore suggests that anchoring the ligand does not necessarily improve the accuracy of the prediction. The effect of various ligand and protein features on the docking performance was investigated. At the level of warhead chemistry, higher success rate was found for Michael additions, nucleophilic additions and nucleophilic substitutions than for ring opening reactions and disulfide formation. Increasing ligand size and flexibility generally affects pose predictions unfavorably, although AutoDock4, FITTED, and ICM-Pro were found to be less sensitive up to 35 heavy atoms. Increasing the accessibility of the target cysteine tends to result in improved binding mode predictions. Docking programs show protein dependent performance suggesting a targetdependent choice of the optimal docking tool. It was found that noncovalent docking into Cys/Ala mutated proteins by ICMPro and Glide reproduced experimental binding modes with only slightly lower performance and at a significantly lower computational expense than covalent docking did. Overall, our results highlight the key factors influencing the docking performance of the investigated tools and they give guidelines for selecting the optimal combination of warheads, ligands, and tools for the system investigated. Results also identify the most important aspects to be considered for developing improved protocols for docking and virtual screening of covalent ligands.



INTRODUCTION Covalent inhibition of therapeutic protein targets has gained increasing interest in the last two decades.1,2 Covalent inhibitors are typically characterized by the formation of a chemical bond between an electrophilic moiety (warhead) on the ligand side and a nucleophilic residue, mostly cysteine, on the protein target. Following the most accepted model, covalent ligands act through a two-step mechanism: the initial formation of a noncovalent complex is essential to place the electrophilic warhead in the right position and orientation with respect to the targeted nucleophilic residue that forms the covalent bond in the second step. Depending on the nature of the electrophile, covalent binding can either lead to reversible or irreversible inhibition, resulting in a remarkably different profile as compared to the noncovalent counterpart. In terms © XXXX American Chemical Society

of binding thermodynamics and kinetics, covalent inhibitors can provide high potency and practically infinite residence time, which give the potential for lower drug dosage and less frequent administration. Nevertheless, the concept received firm scepticism for many years in the past, mainly due to risks of off-target modifications, generation of highly reactive metabolites and potential idiosyncratic reactions arising from the use of hyper-reactive and/or nonselective covalent binders.3−5 In addition to the recognized market share of clinically validated covalent inhibitors,6,7 their recent resurgence has been facilitated by the currently available design principles of safe and effective targeted covalent inhibitors Received: April 15, 2018 Published: June 11, 2018 A

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

Figure 1. Protein targets included in the data set, grouped by most representative target families.

(TCIs).8,9 Suggested strategies still heavily rely on the optimization of noncovalent interactions as primary selectivity keys, with additional focus on tailoring the electrophile reactivity to the desired target nucleophilicity.10,11 Together with the guidelines, several chemical,12,13 biochemical,14 and computational15−18 tools have been introduced for the prospective design and characterization of novel TCIs. A number of computational methods have been integrated to the already wide toolbox available for modeling canonical noncovalent binders. Indeed, since standard molecular docking protocols do not reflect the structural changes occurring upon covalent ligand binding, many software providers implemented new algorithms to handle covalent ligand interactions.19 To our knowledge, comparative analyses of covalent docking tools have been performed only on limited sets of protein−ligand complexes.20−22 The increasing number of structures with covalent ligands available in the Protein Data Bank23 (PDB) prompted us to evaluate these methods on a much wider data set of 207 complexes representing 54 targets. In addition, covalent complexes collected in previous comparative studies only included Michael acceptor and β-lactam type warheads that limited the analysis only to these types of mechanisms. Considering the much broader landscape of warheads and reaction mechanisms available for TCIs,24 we aimed to consider cysteine targeting warheads used in the present set of PDB complexes. We focused on cysteine bound covalent complexes due to the abundance of cysteine targeting covalent inhibitors including FDA approved drugs.2 Cysteines have several features that make them highly suitable for covalent targeting. Their low occurrence (2.3%) in the human proteome25,26 alleviates selectivity issues. Furthermore, the poor conservation of noncatalytic cysteines particularly in protein kinases27 is advantageous for selectivity reasons and to avoid the development of resistance. The cysteine thiol is highly reactive due to its high electron density and polarizability. Therefore, it can be targeted with less reactive ligands allowing to minimize potential side effects.28,29 Cysteines play a significant role in a variety of functions, and they are found

on diverse proteins such as proteases, oxidoreductases, and kinases. Therefore, cysteine-targeted electrophiles can be utilized to affect the function of a wide range of proteins.30 This is supported by proteomic reactivity profiles showing that a broad range of cysteine residues could act as reactive sites.31 Here we describe an unbiased evaluation of six different covalent docking algorithms including the flexible side-chain method developed in AutoDock4, CovDock, FITTED, GOLD, ICM-Pro, and MOE’s built-in covalent docking module. We evaluated the sampling and scoring efficiency of different docking tools in the prediction of experimental binding modes. The conceptually distinct docking procedures face the same set of challenges that consist of predicting traditional nonbonding interactions and modeling optimal geometries of the electrophilic warheads at the reaction site as key requirements for the generation of accurate binding conformations. Furthermore, conventional scoring functions do not allow the evaluation of the contribution of covalent bond formation that typically requires time-consuming quantum-chemical calculations. While developers managed to work around this limitation by relying on molecular mechanics, differences in reaction mechanisms prevent the direct comparison of ligands equipped with different warheads. Therefore, virtual screening applications are limited to warhead categories with similar reactivities. This limitation also applies to self-docking experiments, since the bonding energy is assumed to be nearly constant across a set of ligand conformers. Given the documented success of noncovalent docking in many applications,32−34 we were especially interested in the impact of these limitations on the performance of covalent docking tools and whether they represent a potential to effectively support covalent inhibitor programs.



METHODS Data Collection. Covalent protein−ligand complexes were collected from the PDB. Following the guidelines used for compiling noncovalent data sets35 used for the evaluation of noncovalent docking programs, we restricted our analysis to B

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

Figure 2. Distribution of warhead chemistries included in the docking data set.

binding site and the electron density around the ligand in problematic cases. As a result of this analysis we identified 207 high-quality covalent protein−ligand complexes. Ligand binding modes in this data set are fully supported by ligand electron densities and ligands are in contact with protein atoms from the given crystallographic subunit. Our data set includes 207 unique ligands and 54 targets distributed in different protein classes mainly spanning over proteases (60%), kinases (20%), and GTPases (7%) (Figure 1). Average values and standard deviations of relevant properties, such as molecular weight and lipophilicity and the number of H-bond acceptors and donors as obtained by the calculator plugins of ChemAxon,37 were MW = 421.4 ± 114.2, logP = 2.61 ± 1.63, HBA = 4.7 ± 1.6, and HBD = 2.0 ± 1.1, respectively. The formation of the covalent bond is modeled by matching the reactive atom on the ligand and modifying its structure upon covalent linkage with the targeted cysteine in all of the evaluated docking algorithms. Since the warhead chemistries applied in the covalent ligands involve different reaction energetics and product geometries we analyzed their composition in the data set. Covalent complexes were therefore annotated based on the seven types of reaction mechanisms occurring between the ligand and the protein: Michael addition, nucleophilic addition to nitrile, aldehyde and ketone, nucleophilic substitution, three-membered ring opening, and disulfide bond formation (Figure 2). The algorithms evaluated in this study present a variety of approaches to handle the binding mode prediction of covalent ligands. In order to ensure a fair head-to-head challenge we used predefined parameters in all cases. A maximum of ten best scoring ligand conformations were saved in each docking simulation to investigate both the scoring and sampling efficacy of the methodologies. Root-mean-squared deviation (RMSD) values were calculated by the Python script developed at Schrödinger (rmsd.py) providing pairwise symmetry-corrected RMSD between the cocrystallized ligand and the docking pose.

high-quality protein structures. Therefore, we selected PDB entries deposited after January 1, 1995, having a resolution of 2.5 Å or better with structure factors available. Limitation on the date reflects the technology improvements in the measurement, structure solving, and analysis of protein−ligand complexes observed in the last decades. Criteria in resolution gave further confidence on the overall quality of protein structures while structure factors allow the reproduction of electron density maps found to be critical in the evaluation of experimental poses. Next, we identified entries with ligands covalently bound to their protein target via cysteine analyzing the connectivity information available in the PDB files. Connectivity information was inspected in the advanced search interface of PDB by the “Link records” option in the “Structure Features” menu. Cysteine targeting covalent ligands were identified using “CYS” as “Component” and “SG” as “Atom label” that retrieved complexes in which a cysteine residue has been modified. High quality covalent protein−ligand complexes were filtered further using ligand information. Our evaluation of covalent docking tools is focused to synthetic ligands relevant for covalent drug discovery. Therefore, we analyzed the ligand properties using chemical pattern and physicochemical property filters. Analyzing the content of HET groups in the PDB files, we selected ligands that contain oxygen and nitrogen and restricted the suitable atom types to H, C, N, O, F, P, S, Cl, and Br. Crystallization agents, buffer components, cofactors, saccharides, tetrapeptides, or bigger oligopeptides were excluded. The remaining set of ligands were then subjected to physicochemical filters that identified ligands with heavy atoms between 8 and 50 fulfilling the extended Rule-of-Five criteria (molecular weight ≤ 700, 0 ≤ cLogP ≤ 7.5, number of H-bond donors ≤ 5, polar surface area ≤ 200 Å2, number of rotatable bonds ≤ 20).36 Finally, the resulting set of complexes were further analyzed by visual inspection. During this procedure we investigated a large number of complexes and considered the types of warheads together with their observed binding geometry, the protein contacts of the ligands within the C

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

generate PDBQT files, to process flexible and rigid components, to calculate grid maps, and to write out docking parameter files. The ligand was treated as a fully flexible residue by sampling it as part of the receptor in order to optimize its interactions within the binding pocket. Docking simulations were performed using default Lamarckian Genetic Algorithm settings. Ligand docking conformations were ranked by a semiempirical force field-based scoring function. CovDock. CovDock is the covalent docking workflow integrated in Schrödinger’s Small-Molecule Drug Discovery Suite. It is based on an automated multistep process devised to capture the essential features of the binding mechanism of covalent ligands. CovDock relies on a SMARTS definition of the reaction type by performing the docking simulation on all of the reaction sites matching the pattern. The method has been developed in two versions: a virtual screening mode (CovDock-VS),43 tailored to address throughput needs, and a “Pose Prediction” mode,20 enabling better accuracies for higher computational costs. For the purpose of a self-docking experiment, we used CovDock in the pose prediction mode. Its workflow starts with a Glide44−46 noncovalent docking simulation using positional constraints for pose selection (5 Å distance cutoff between the two atoms involved in the bond formation). In detail, the default protocol uses ConfGen47 to sample ligand conformations prior to docking, where the prereaction ligand is docked to the binding site mutating the reactive residue to alanine. This allows to bring the warhead close to the targeted cysteine by avoiding unfavorable clashes. The mutation is then reversed, and rotamer states of the nucleophilic cysteine residue are sampled to form a covalent bond with different top scoring noncovalent ligand poses based on geometric criteria. In the default pose prediction mode evaluated in this comparative study, covalent docking poses are scored and ranked by the Prime VSGB2.0 energy model, used to restore optimal physical geometries following covalent bond formation.48−50 The Affinity Score, an additional scoring function provided by Schrödinger to compare different compounds with intrinsically similar reactivities, combines both parts of the process, as it is calculated as the average of the initial prereaction and the postreaction in-place GlideScore. Note that Prime VSGB2.0 energy model was used throughout the study while Affinity Score was only used in a comparative analysis with the default scoring function. FITTED. FITTED51−53 is a docking tool that has recently been integrated in the FORECASTER Platform. Its latest implementation supports covalent docking and screening libraries with covalent and noncovalent ligands. This unique feature is based on an automated identification of the ligand warhead atoms, thus allowing to discriminate between the two binding mechanisms. Moreover, a covalent bond is formed during the docking simulation only if the warhead is predicted to be in close proximity to the targeted residue. Otherwise the ligand is considered to be a noncovalent binder regardless the presence of the warhead in its structure. However, FITTED does not allow to add custom warheads to the default set due to a detailed mechanistic parametrization of the reaction mechanisms. Therefore, about 15% of our data set (i.e., ligands reacting via nucleophilic substitution, ring opening and disulfide bridge formation) could not be evaluated at the time of analysis due to their missing implementation. Ligands were prepared through the SMART program, while protein structures were prepared and processed by means of the PREPARE and PROCESS modules, respectively. A genetic

The full list of PDB accession codes and information is available in Supporting Information Table S1. Ligand and Protein Preparation. In order to avoid any bias deriving from the binding conformation observed in the crystal structure, ligand structures were generated from 2D and prepared by using the corresponding preparation protocol of the given docking program, followed by an unconstrained energy minimization before docking. In case of prochiral ligands, CovDock, FITTED, ICM-Pro, and MOE automatically evaluated all stereoisomers possibly generated upon cysteine attachment, while AutoDock4 and GOLD were manually instructed to evaluate all stereoisomers. Protein structures were prepared according to the specific protein preparation protocol associated with the different docking tools. In all cases, protein structures were obtained by removing protein domains that were irrelevant for ligand docking. In addition, mimicking a real life virtual screening situation, cocrystallized covalent ligands, water molecules, metals, cofactors, and crystallization agents were removed from target structures. Binding site residues were defined by centering the box on the targeted cysteine. In particular, AutoDockTools was used to prepare flexible and rigid protein components for AutoDock4,38 by adding hydrogens and Gasteiger charges and then generate the corresponding structural files in the PDBQT format. Protein Preparation Wizard provided in the Schrödinger Suite39,40 was applied for the preparation of target structures for CovDock simulations. It was used to add hydrogens, to determine protonation states, to refine the protein’s H-bond network, to correct side chains orientations, and to perform a restrained minimization of hydrogen atoms. Protein structures for FITTED docking were prepared by the PREPARE module in the FORECASTER platform.41 It was used to add hydrogens, for protonation state generation and for structural optimizations. Default settings were used in GOLD’s preparation wizard in GOLD to add hydrogen atoms and assign protonation states. In ICM-Pro’s environment, protein structures were converted into ICM objects and processed to optimize H-bond network and to further optimize the orientation of other residues (His/Pro/ Asn/Gln,Cys). For MOE docking, protein structures were prepared by using MOE QuickPrep functionality with default settings to add hydrogens and partial charges, correct Asn/ Gln/His orientations, optimize the H-bond network (Protonate 3D), and perform an energy minimization. AutoDock4. AutoDock4 implemented two algorithms to model covalent protein−ligand systems including the two-point attractor and the flexible side-chain methods.42 Based on the comparative evaluation of 20 covalent complexes equipped with 5 warhead chemistries, the two-point attractor method was clearly outperformed by the latter. Therefore, we used the flexible side-chain method in AutoDock4 to assess the binding mode prediction on our set of covalent complexes. It requires the initial modification of the ligand structure by connecting the two terminal atoms of the reactive side-chain, Cβ-SG in the case of a cysteine, at the site of alkylation. Hence, the two additional atoms were used for the alignment on the nucleophilic cysteine by specifying a SMARTS pattern and the indices of the overlapping atoms through the code available at http://autodock.scripps.edu/resources/covalentdocking/ (accessed March 12, 2018). In order to avoid automated and undesirable bond formation between other ligand-protein atoms particularly close attention must be paid to the geometry of the initial modification. AutoDockTools was used to D

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

3O1G; 3O6T; 3OF8; 3OVX; 3RHY; 3SN8; 3SVV; 3SZ9; 3SZB; 3T9T; 3V4J; 3W2P; 3ZV9; 3ZVG; 4BPV; 4BQV; 4BS5; 4BS6; 4BSQ; 4CDC; 4CDD; 4CDE; 4CDF; 4D9U; 4DCD; 4DMX; 4DMY; 4GMX; 4GS6; 4HCU; 4HCV; 4I24; 4L1O; 4LUC; 4LV6; 4LYF; 4LYH; 4M1O; 4M1S; 4M1T; 4M1W; 4M1Y; 4M22; 4MZ4; 4MZO; 4MZS; 4PE0; 4PE1; 4PE4; 4PIQ; 4PIS; 4PNC; 4QPS; 4THI; 4TWY; 4US2; 4WVF; 4WX4; 4WX7; 4X6J; 4XCU; 4YHF; 4YRT; 4ZZM; 4ZZO; 5C1U; 5C1X; 5C1Y; 5C91; 5DAF; 5DGJ; 5DP4; 5DP5; 5DP7; 5DP8; 5DP9; 5DPA; 5E7R; 5EZS; 5F02; 5F2E; 5FX6; 5GIT; 5GNK; 5HG5; 5HG7; 5HG8; 5HG9; 5HZE; 5IYT; 5J7S; 5J87; 5J9Z; 5JH6; 5JK3; 5L6P; 5LCK; 5LWN; 5MAJ; 5MJB; 5MQY; 5N0Y; 5O8U; 5O8V; 5P9J; 5P9K; 5P9M; 5TDI; 5TG2; 5TOZ; 5TTS; 5TTU; 5TTV; 5UG8; 5UGC; 5V6S; 5V6V; 5V88.

algorithm is used also by FITTED to drive the conformational search at the ligand side. Ligand conformations are ranked by energy. Therefore, the lowest energy is used to identify the best docking pose in the simulation. GOLD. GOLD54−56 is a genetic algorithm based automated docking program developed at the Cambridge Crystallographic Data Centre (CCDC). Its covalent implementation uses postreaction ligand conformations by relying on the superposition of a “link-atom” in both the ligand and protein structures. Input files should be modified manually to attach the reactive sulfur atom of the targeted cysteine to the electrophilic atom of the ligand. The link-atom of the ligand is forced to fit onto the link-atom of the protein and a bending energy term is calculated from the corresponding parameter file. Eventually, this energy term is added to the fitness docking score in order to ensure a correct geometry for the covalent linkage. Other contributions to the final score are calculated from the conventional noncovalent interactions generated by sampling the torsional degrees of freedom of the ligand within the binding pocket. Ligand docking conformations are ranked by the GOLD Fitness Score. The best docking pose in the simulation is assigned the highest (most positive) score. ICM-Pro. Molsoft LLC developed the docking code ICMPro57 that includes a unique protocol for covalent docking.58 Based on reaction type, ICM-Pro converts automatically the prereaction ligand into a “pseudo-ligand” by covalently attaching the side-chain of the cysteine to the electrophilic warhead in all the stereoisomers produced upon addition. Sidechain atom coordinates on the protein are first used to constrain the position of the pseudoligand and then removed to avoid unfavorable clashes during docking. Ligand conformations are sampled through Monte Carlo simulations in a set of grid maps calculated for the pocket.59 Finally, docking solutions are scored and ranked with a modified version of the full ICM scoring function excluding pairwise interactions among the atoms that are directly neighbor the covalent link. MOE. The covalent docking module of MOE (Molecular Operating Environment)60 relies on a reaction/transformation placement methodology to match the reactive groups on the ligand and the cysteine residue followed by the formation of the covalent bond between them. Protein and ligand structures were prepared by using default parameters. Covalent ligands were modeled and minimized in the prereaction form generated with the Builder panel. MarvinSketch was used to draw chemical reactions37 that were not already described in the predefined set. The default postprocessing protocol using minimization by the Rigid Receptor refinement scheme and rescoring by the GBVI/WSA dG scoring function61 was applied to estimate the free energy of ligand binding in kilocalories per mole. Ligand poses were ranked by the estimated binding free energy. The following is a list of the PDB codes used in this work: 1AIM; 1AU3; 1AU4; 1BGO; 1BMQ; 1CVZ; 1EWL; 1EWM; 1EWO; 1EWP; 1F29; 1F2A; 1F2B; 1F2C; 1F4C; 1M6D; 1ME3; 1MEM; 1MS6; 1NLJ; 1NPZ; 1NQC; 1Q6K; 1RE1; 1RHJ; 1RHM; 1RWN; 1RWP; 1RWX; 1THE; 1TU6; 1U9V; 1U9W; 1YK7; 1YT7; 2AIM; 2ALV; 2AMD; 2AUX; 2AUZ; 2AWZ; 2AX0; 2AX1; 2BDL; 2CNO; 2DCC; 2DW5; 2F7D; 2FQ9; 2FRA; 2FRQ; 2FT2; 2FYE; 2G6D; 2G8E; 2GH5; 2GX4; 2HWO; 2HXZ; 2IJN; 2OP3; 2OP9; 2R6N; 2R9F; 2XM7; 2XU1; 2XU4; 2XU5; 2XYG; 2XYP; 2YJ2; 2YJ8; 2YJ9; 2ZU5; 3B0R; 3B1U; 3BLT; 3BLU; 3C9W; 3H0E; 3HD3; 3HHA; 3HWN; 3I06; 3KW9; 3KWB; 3KWZ; 3N3G; 3N4C;



RESULTS AND DISCUSSION Overall Success Rates. Covalent docking methods were first evaluated to predict the experimental binding mode found in the 207 covalent complexes. The prediction accuracy was expressed in terms of the ligand heavy atoms RMSD values from the conformation cocrystallized in the PDB entry. Docking results were analyzed by calculating the RMSD of the best scoring solution (Top1), and the best RMSD among top ten scoring poses (Top10). Differences between Top1 and Top10 results indicate scoring errors, since the high quality binding pose is generated but not scored as the best. Low quality Top10 results might suggest an inappropriate sampling, since the experimental pose is not recovered, at least not in the top ten scored solutions. In addition, our evaluation on the best scoring pose (Top1) might be a surrogate performance metric in virtual screening situations considering only this pose for ligand ranking. As shown in Figure 3, when taking a RMSD cutoff of 2.0 Å as the criterion for successful predictions, ICM-Pro was able to retrieve the experimental conformation in Top1 in 62% of the

Figure 3. Overall success rates in the prediction of experimental binding modes. Rates are calculated as the percentage of complexes for which a ligand pose was predicted within the specific RMSD cutoff. Bright colors refer to success rates in the best scoring (Top1) docking conformation; lighter colors to rates within Top10 poses. The gray dashed area stacked on FITTED represents the percentage of complexes for which a covalent mechanism is not implemented yet (≈15%). E

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

sampling. GOLD and MOE have shown a large gap between the two sources of error, with 35% and 49% of failed predictions caused by a wrong sampling, respectively. The difference between the potential error sources was less pronounced for AutoDock4. Inaccurate sampling is obviously most remarkable for MOE because of the higher number of complexes that were erroneously predicted by its covalent docking module (≈63%). It should also be highlighted that MOE could provide ten docking solutions only in 56% of the cases, which further supports the influence of its sampling flaws on the performance. Overall, the performance of GOLD was more affected by sampling errors as compared to other genetic algorithm based docking tools, thus questioning its default conformational search. Finally, 26% of the failed Top1 predictions generated by CovDock were caused by sampling errors. Based on its workflow, it samples ligand conformations only on the prereaction form, i.e., before the covalent bond formation. A more detailed analysis will be provided later by inspecting the noncovalent step (Table 2). Cases with Scoring Errors. The complex 4us262 represents an example of scoring error, since four out of six covalent docking methods could retrieve a near-native conformation in Top10, but not in Top1 (Figure 4). In particular, the fragment-sized maleimide derivative covalently bound to Cys118 in the H-Ras:SOS complex does not make significant noncovalent interactions. This obviously contribute to the observed binding conformation, except for a hydrogen bond with Asn85 characterized by a rather suboptimal geometry. In addition, the overall structural change of the electrophile upon binding makes the prediction of its binding mode especially challenging. As a consequence, none of the programs could correctly predict its geometry in the best scoring pose suggesting that this error was likely due to an incorrect modeling of the geometry of the electrophile at the attachment point. As for GOLD, whose deviation was only slightly larger than the accepted threshold, this case can be rather considered as an exception, since it proved to be the best in predicting fragment-sized covalent inhibitors in Top1. In the case of the 2yj263 complex instead, the scoring error observed with CovDock, ICM-Pro, and MOE resulted from an inaccurate evaluation of the noncovalent interactions driving the binding mode of the nitrile-based covalent inhibitor of Cathepsin L1 (Figure 5). In particular, CovDock and MOE were not able to score Top1 ligand conformation with a halogen bond between the aryl-bromide and the CO group of Gly61 in the S3 pocket. In the case of ICM-Pro, the docking pose wrongly predicted in Top1 completely lacked the H-bond framework formed by the amide group in the vicinity of the warhead, therefore producing a large deviation from the experimental pose (12.6 Å). The same applies to the best scoring conformation observed with GOLD, while in the case of FITTED the Top1 pose missing the halogen bond also provided the lowest RMSD within ten solutions. In this example, the flexible side-chain method in AutoDock4 was the

complexes, followed by CovDock (59%), AutoDock4 (55%), and GOLD (53%). Since FITTED deals with a limited set of covalent chemistries, its success rate (48% in Top1) was normalized on the whole set to ensure a fair comparison among the tools (56% and 81% for Top1 and Top10 on actual set, respectively; Figure S1). A much lower performance was instead observed for MOE, which could make an accurate prediction only in 37% of the cases by using default settings. By changing the scenario to Top10 docking solutions, we could systematically retrieve more near-native conformations. In particular, ICM-Pro performed the best retrieving 88% of the ligands within 2.0 Å from the corresponding experimental pose, followed by AutoDock4 and CovDock (≈75%). To examine if there was a statistically significant difference between the various tools, we applied a Mann−Whitney U-test for the RMSD sets obtained. The analyses were performed pairwise both within the Top1 and the Top10 sets (Supplementary Table S2). The majority of the probability values (p) were less than 0.05 indicating statistically significant differences at the p < 0.05 level. p > 0.05 was obtained when AutoDock4, GOLD, and FITTED, the three tools using genetic algorithm to explore the conformational space, were pairwise compared. A lower level of differences is seen for ICM-Pro, CovDock, and FITTED, and this is in line with their similar (and top) performance. The differences between the Top1 and Top10 RMSD sets of each individual docking tool were also statistically analyzed, and they were all found to be different at the p < 0.05 significance level (Supplementary Table S3). A more detailed analysis of the results aimed to identify possible scoring and sampling errors of the docking tools. We considered as a scoring error when a near-native ligand conformation was retrieved within ten docking solutions but not scored as the best. A sampling error, instead, refers to docking that could not provide the correct ligand conformation within the top ten docking poses. The summary provided in Table 1 shows how the majority of unsuccessful results in Table 1. Distribution of Scoring and Sampling Errors Across Failed Top1 Covalent Docking Predictions When Considering 2.0 Å as the RMSD Cutoff covalent docking method

data set size

% error

% scoring error

% sampling error

AutoDock4 CovDock FITTED GOLD ICM-Pro MOE

207 207 175 207 207 207

45% 41% 44% 47% 38% 63%

20% 15% 25% 12% 26% 13%

25% 26% 19% 35% 12% 49%

FITTED and ICM-Pro (≈25%) are due to an imperfect scoring of a ligand binding conformation close to the cocrystallized structure, since a different docking pose was scored the best. For the remaining methods, a higher number of wrong predictions arises from an inaccurate conformational

Table 2. Noncovalent Docking Results on Failed Covalent Docking Predictions failed covalent predictions in Top1

failed covalent predictions in Top10

noncovalent docking method

successful predictions in Top1

successful predictions in Top10

successful predictions in Top100

successful predictions in Top1

successful predictions in Top10

successful predictions in Top100

Glide SP ICM-Pro

37% 32%

57% 71%

65% 86%

34% 12%

49% 56%

59% 68%

F

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

Figure 4. Scoring error arising from docking results on the PDB entry 4us2, a maleimide bound to H-Ras:SOS complex. The cocrystallized ligand is shown in purple sticks, the best docking pose (Top1), in orange, and the best pose in Top10, in green.

hinge is basically responsible for the sampling error detected, although misplacement of the solvent exposed tail of the ligand also contributed to the docking failure in certain cases. It is likely that such sampling errors could be effectively reduced by appropriate pharmacophore constraints as pre- or postdocking filters.43 Indeed, the well-known binding mechanism of kinase inhibitors would suggest the application of H-bond filters in proximity to the hinge region in order to ensure that only ATP-competitive conformations would be sampled in the docking simulation. Without these filters, MOE was exceptionally good in making an accurate prediction already in the best scoring pose, while other tools found it difficult to sample this rather flexible ligand. Factors Influencing the Docking Performance. Next we analyzed the impact of ligand and target properties on the performance of covalent docking tools. From the ligand side we considered the warhead chemistry and the physicochemical properties. Targets were evaluated by considering the accessibility of targeted cysteines and the volume of the binding site that accommodates them. Finally, the effect of the protein and the impact of noncovalent interactions were investigated. Warhead Chemistry. The type of the electrophilic warhead applied in the covalent binder influences the geometry of the bound ligand (Figure 8). Considering only the top scoring solution, ICM-Pro has shown the best performance (71%) in modeling Michael additions, the most represented

only one that properly identified the experimental binding mode with all the noncovalent interactions already in Top1. Cases with Sampling Errors. A specific example of sampling error caused by the improper treatment of structural changes occurring in the electrophile geometry upon covalent bond formation was found for HCV NS5b RNA polymerase structures (2ax0, 2ax1, and 2awz).64 The mechanism of inhibition in these cases is characterized by a large conformational change upon ligand binding. When considering its prereaction form, the rhodanine-based warhead is indeed located at the center of a highly conjugated planar system that gets broken upon covalent linkage to Cys366, leading to a structural rearrangement into a more packed conformation. As shown in Figure 6, most of the docking poses were predicted in an extended conformation that was clearly derived from the pool sampled from the initial prereaction form. Not surprisingly, the few accurate predictions were only made by protocols (i.e., AutoDock4 and GOLD) using a postreaction ligand structure, where conformational sampling takes place on the ligand showing no π-conjugation in the warhead. A distinct type of sampling error instead emerges for 5p9k,65 an acrylamide inhibitor bound to Bruton’s tyrosine kinase (BTK) (Figure 7). Indeed, the wrong predictions observed in this case were mainly the consequence of an inaccurate sampling of the 5-fluoro-pyrimidine scaffold involved in a noncovalent interaction with the backbone NH of Met477 in the hinge binding region. The lack of the hydrogen bond at the G

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

Figure 5. Scoring error arising from docking results on the PDB entry 2yj2, a nitrile-based covalent inhibitor bound to Cathepsin L1. The cocrystallized ligand is shown in purple sticks, the best docking pose (Top1), in orange, and the best pose in Top10, in green.

Ligand Descriptors. We were interested in exploring to what extent ligand properties influenced the performance of covalent docking methods. Our data set is well suited for this type of analysis since it contains ligands with a wide range of size and flexibility. By investigating the impact of the ligand size (Figure 9), a distinct descending trend in the success rate was identified for GOLD and MOE, while a steady performance (around 60%) was observed for AutoDock4 up to 35 heavy atoms. Within the same limit, ICM-Pro and FITTED have shown some improvement by increasing ligands size, indicating that they can properly capture the native binding mode of bulky structures. Fragment-sized covalent binders, however, were best retrieved by GOLD in Top1. Focusing to the difference between trend lines in Top1 and Top10 we could generally see that scoring errors are typically larger for fragments. This part of the set is particularly interesting for fragment-based ligand discovery (FBLD) applications aimed to the identification of covalent drug binders. FBLD is indeed considered a well-established approach providing several drugs and candidates in the last decades,66,67 and its role in the design of covalent drugs has recently been enabled by new methods for screening fragment size reactive probes directly in proteomes and cells.68 Analyzing the impact of ligand flexibility, a clear drop in the prediction rates was generally observed by the increasing number of rotatable bonds (Figure 10). ICM-Pro showed much less dependence on the number of rotatable bonds. In

type of reaction in the library. Together with CovDock (68%), it hence represents the method of choice for a virtual screening of compound libraries containing Michael acceptors. When considering nucleophilic additions involving aldehyde-, ketone-, and nitrile-based warheads, all of the methods successfully retrieved at least half of the experimental poses (≈60% for ICM-Pro and FITTED), except for MOE, whose rate virtually matched the overall one (39%). Although the reaction mechanism of these warheads is identical, different parametrization and product geometry required to analyze them separately. In this sense, ICM-Pro and FITTED have shown the highest success rate on the aldehyde subset (67%), GOLD performed better with regard to ketones (58%), while AutoDock4 proved to be most reliable on nitriles (62% in common with ICM-Pro), slightly more than FITTED (60%) and CovDock (58%). These findings suggest that, when screening libraries of these electrophiles, the use of the specific tools could provide better hit rates. All of the methods made a successful prediction for nucleophilic substitutions in nearly half of the complexes. On the contrary, poor performances were registered for both the scarcely represented threemembered ring opening reactions and disulfide bond formation, with GOLD being the best exception, as it was able to predict a reliable docking pose for 4 out of 5 ligands reacting through a thiol oxidation. As mentioned, these warheads were out of the scope of FITTED. H

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

Figure 6. Sampling error detected for 2awz, a rhodanine derivative bound to HCV NS5b. Large conformational changes occurring upon protein− ligand reaction can represent a major challenge for covalent docking protocols. The cocrystallized ligand is shown in purple sticks, the best docking pose (Top1), in orange, and the best pose in Top10, in green.

particular, we found that the increased flexibility of the ligand raised the relevance of sampling errors for CovDock, GOLD, and MOE. In AutoDock4, the widening range between Top1 and Top10 indicates that scoring errors contribute almost equally to the failed predictions (c.f. data shown in Table 1). It should be noted that the normalized trend lines of FITTED are related to a smaller set and therefore could not be compared directly to that of the other tools. As a complementary descriptor, we inspected the effect of the pharmacophore content on the covalent docking performance to estimate the impact of noncovalent interactions on the prediction accuracy (Figure 11). Using e-Pharmacophores69,70 of Schrö dinger, we analyzed the energetically favorable pharmacophore features available in the protein−ligand complexes. By increasing the number of pharmacophores we observed a drastic improvement in the success rates of FITTED and ICM-Pro. This finding indicates that their performances become more accurate when the binding mode is stabilized by multiple noncovalent interactions within the pocket. As for the remaining methods, we could not identify any distinctive trend other than a consistency of the prediction. Cysteine Accessibility and Binding Site Volume. On the protein side, docking results were analyzed by means of both the solvent accessible surface area (SASA) of the reactive cysteine and the binding site volume. The protein surroundings at the reaction site is indeed a key factor to consider not only for the electronic effect it exerts on the warhead reactivity

and cysteine nucleophilicity, but also for the steric influence on the accessibility of the nucleophilic cysteine. Moreover, treatment of nonbonded interactions of atoms at the reaction center affects the scoring accuracy. By crossing docking results with cysteines’ SASA, the rate of successful Top1 predictions was generally improved by increasing the surface accessibility (Figure 12a). This is in line with the expectation that a more spacious reactive site allows the ligand warhead to be placed in the correct position in the cysteine surroundings. Steric clashes occurring in more densely packed subpockets represent a challenge for the accuracy of the binding mode prediction. Only FITTED has deviated from this trend, since the majority of its errors were generated with SASA higher than 30 Å2. Interestingly, by looking at the distribution of warheads across the different groups of SASA, we found that the largest part of electrophiles used to target less accessible cysteines react through electrophilic addition (75%). By contrast, Michael acceptors were mostly used to target residues with a larger surface accessibility (70% for SASA higher than 30 Å2). As previously seen for warhead chemistries, the prediction accuracy of covalent docking programs has generally improved by increasing the occurrence of Michael additions in the subset, demonstrating that both factors (warhead geometry and cysteine accessibility) influence the performance of docking. As a complementary analysis on the noncovalent counterpart, we inspected the impact of the binding site volume on the I

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

Figure 7. Docking results for the PDB entry 5p9k, an acrylamide inhibitor bound to the Bruton’s tyrosine kinase. The cocrystallized ligand is shown in purple sticks, the best docking pose (Top1), in orange, and the best pose in Top10, in green.

and CatS (70%), while GOLD achieved above 65% and 75% on CatK and K-Ras, respectively. FITTED proved to be reliable with CatL1 (65%), EGFRs (77%), and A9XG43 (88%); ICM-Pro, instead, performed best with K-Ras (85%) and could exhaustively predict accurate conformations on EGFRs and A9XG43. MOE, on the other hand, has shown a success rate above 75% only for K-Ras. Among less represented targets with 3−9 structures, GOLD showed high performances against BTK (80%), CatC (100%), PADI4 (100%), and ITK (100%, shared with FITTED). On the other hand, ICM-Pro correctly predicted all of the experimental conformations of ligands bound to JAK3, TAK1, Caspase-1, ERK2, and Adenain. By focusing on protein families, proteases-bound ligands were overall retrieved with a success rate around 50%, except for MOE which did not exceed a 30% rate. Kinases were better predicted by ICM-Pro (90%), followed by AutoDock4, FITTED, and CovDock (≈70%). Finally, CovDock and ICM-Pro showed a success rate around 80% against GTPases, followed by GOLD and MOE (≈70%). Comparison to Noncovalent Docking. In an attempt to identify the source of covalent docking failures, we investigated those cases in which ICM-Pro and CovDock, the two methods having the best rate of predictions within 2.0 Å on the whole set, could not accurately reproduce the experimental conformation of the reference ligand. In particular, the multistep protocol used in CovDock suggested to analyze whether the first noncovalent docking stage would result in

performance as it is expected that a higher number of interactions would be permitted in larger pockets. In line with the trend observed when considering the effect of the pharmacophore content, ICM-Pro has shown a linear improvement in the success rate by increasing the size of the pocket (Figure 12b). The other tools do not appear to be affected to the same extent, although for all of them the best performance was registered for volumes larger than 300 Å3. In particular, CovDock and MOE seem to be less dependent on the binding site volume considering Top1 solutions. It is worth noticing that methods using genetic algorithms for conformational sampling (AutoDock4, FITTED, and GOLD) shared a highly similar behavior by changing the binding site volume. Impact of the Protein Target. Next we investigated the impact of protein targets on the prediction accuracy of covalent docking tools (Figure 13). Considering only targets having at least ten structures (i.e., Cathepsin L1, K-Ras, Cruzain, Cathepsin S, and Cathepsin K), the best average success rate of docking predictions within 2.0 Å was shown by CovDock in Top1 (64%) and ICM-Pro in Top10 (93%). Analyzing Top1 results for the same set of proteins revealed the target specific performance of the programs and might help to find a suitable tool for virtual screening against the target of interest. For example, AutoDock4 retrieved a near-native pose for all of the Enterovirus A71 genome polyprotein structures (UniProt ID: A9XG43), CovDock showed the highest success rates in covalent docking against K-Ras (85%), Cruzain (77%), J

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

Figure 8. Docking performance (% success rate with 2.0 Å as RMSD cutoff) obtained with different warhead chemistries indicating the number of complexes per reaction type in brackets. The gray area in FITTED refers to reaction mechanisms for which covalent docking is not supported yet (≈15%).

better binding mode predictions as compared to that modeling the covalent bond formation. Following the default workflow in the pose prediction mode, we evaluated the maximum of 100 poses obtained by docking the prereaction ligand structure against the protein target having its reactive cysteine mutated to alanine. As shown in Table 2, Glide SP was able to retrieve 65% of the failed covalent docking predictions, with more than one-third already as the best scoring conformation (37%). Considering only complexes with unsuccessful predictions from Top10 poses gave highly similar rates (59% in Top100, 34% in the best scoring pose). These results show that the Prime refinement of the covalent docked complex might be responsible for the loss of accurate docking poses identified in the previous noncovalent step with Glide. Next we analyzed the effect of different scoring functions (Affinity Score and Prime Energy) on ranking covalent docking poses. While the Prime Energy is used for ranking different docking poses of the same ligand (also in this study), the Affinity Score is calculated as the average of prereaction and postreaction GlideScore and is aimed to rank different compounds having similar intrinsic reactivities. However, the high percentage of native conformations retrieved with Glide prereaction docking for cases failed by the Prime Energy function prompted us to investigate the impact of scoring functions in depth. We found that the scoring functions gave highly comparable success rates both in Top1 and in Top10 (Figure 14). Mann−Whitney U-test was applied to the RMSD sets to analyze if there is a statistically significant difference between these scoring functions. The

probability values (p) obtained for the Top1 and Top10 sets are p = 0.72 and 0.76, respectively, indicating that there is no significant difference between the functions (Supplementary Table S4). We performed a similar analysis on ICM-Pro Top1 covalent failures using its conventional noncovalent docking protocol against Cys-to-Ala mutated proteins (Table 2). Also in this case, remarkable rates of accurate predictions were observed, with a total of 86% near-native conformations retrieved by noncovalent docking. More than two-thirds of these were among Top10 docking solutions (71%) and one-third already as the Top1 pose (32%). By focusing on the failed covalent docking outcomes in Top10 (25 complexes in total), ICM-Pro could make a correct prediction in the Top1 pose only in three complexes (12%) and showed a significantly improved performance (56%) in the Top10 poses. These results show that the complexes concerned are subject to sampling error in covalent docking and to scoring error in noncovalent docking. Finally, these results show that despite of the excluded direct contribution of the reaction center, the constrained pseudoligand approach might provide unfavorable electrophile geometries, thus affecting the overall binding conformation. Next we investigated the performance of noncovalent docking algorithms on the whole set of covalent complexes (Supplementary Table S5). In a previous analysis, noncovalent docking with the SCAR method has been performed for 75 covalent complexes.74 It was found that noncovalent docking into cysteine to glycine mutated proteins gave mean RMSD and percent of docked poses within 1 and 2 Å limits K

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

Figure 9. Impact of ligand size on the covalent docking prediction, expressed as success rate (2.0 Å as RMSD cutoff) against number of heavy atoms. Black and red lines refer to rates in Top1 and Top10, respectively. Ligands distribution across the different groups of heavy atoms is shown as a light blue area. Success rates and ligands distribution for FITTED are normalized to the actual set of complexes for which a covalent docking could be performed.

Figure 10. Impact of ligand flexibility on the covalent docking prediction, expressed as success rate (2.0 Å as RMSD cutoff) against number of rotatable bonds. Black and red lines refer to rates in Top1 and Top10, respectively. Ligands distribution across the different groups of rotatable bonds is shown as a light blue area. Success rates and ligands distribution for FITTED are normalized to the actual set of complexes for which a covalent docking could be performed.

comparable to those of covalent docking tools. We found that both Glide SP and ICM-Pro in Cys-to-Ala mutated proteins were able to make an accurate best scoring (Top1) prediction

in around 53% of the covalent complexes (Table 3), approximately 6% and 9% less than the respective covalent docking modules. Considering the Top10 solutions, the gap to L

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

Figure 11. Effect of the pharmacophore content on the covalent docking prediction, expressed as success rate (2.0 Å as RMSD cutoff) against the number of pharmacophore features in the complex. Structure-based pharmacophore models were generated with the e-Pharmacophores method in Schrödinger. Black and red lines refer to rates in Top1 and Top10, respectively. Ligands distribution across the different groups of features is shown as a light blue area. Success rates and ligands distribution for FITTED are normalized to the actual set of complexes for which a covalent docking was performed.

Figure 12. Impact of protein descriptors on the covalent docking performance: (a) solvent accessible surface area (SASA) of the reactive cysteine against success rate in Top1. SASA calculations were performed on unbound protein structures with POPS.71 Warhead distributions on the whole set (thus not relevant for FITTED) across the different groups of CYS SASA are shown on top of the respective fields. (b) Binding site volume against success rate in Top1. Volumes were calculated with SiteMap72,73 in the Schrödinger Suite. Success rates are based on a RMSD tolerance of 2.0 Å. Success rates in FITTED are normalized to the actual set of complexes for which a covalent docking was performed.

covalent counterparts was further reduced (2% and 6% difference for Glide and ICM-Pro, respectively). Overall, these findings suggest that dedicated covalent docking tools can improve the performance observed by using noncovalent methods neglecting the structural arrangements at the reaction center. However, the modest improve-

ment observed might question the use of covalent docking tools that sometimes require higher computational costs. This question might be especially important in virtual screening of large electrophile libraries. In contrast to the abundance of reported virtual screening applications aimed to the discovery of noncovalent ligands,75 there are only few cases of M

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

Figure 13. Performance distribution across proteins having at least three representative structures. Success rates are calculated as the percentage of docking predictions with a RMSD within 2.0 Å from the reference ligand. Targets are sorted by members count in ascending order, with Cathepsin K being the most represented group in the data set. “FITTED (not supp.)” complements the adjacent column, as it refers to the percentage of complexes for which a covalent docking was not supported.

poxvirus proteinase I7L inhibitors.58 Screening a virtual library of about 230 000 ketones and aldehydes against a homology model of the I7L ligand binding site, the authors selected 456 virtual hits by their docking score and chemical properties. Biochemical testing confirmed 97 inhibitors at 200 μM concentration, thus yielding a hit rate of around 21%. In another prospective application, a structure-based virtual screening aimed to the identification of Cathepsin K (CatK) covalent inhibitors was performed by means of the covalent docking module developed in GOLD.76 A focused library of 1247 compounds was screened against a CatK X-ray structure previously validated for binding mode prediction. In vitro testing of 39 predicted hits led to the identification of 21 actives (hit rate higher than 50%), 3 of which in the nanomolar range. Also FITTED was used in a prospective screening for the discovery of prolyl oligopeptidase (POP) inhibitors covalently bound to the active site Ser554.53 Starting from a focused data set of 2798 molecules equipped with either an aldehyde or a nitrile group as warhead, the authors inspected the rank-ordered library and identified one hit for further structure-based optimizations, finally leading to novel selective POP covalent inhibitors. We analyzed the information available on the compounds selected from the prospective virtual screening applications of ICM-Pro and GOLD in light of our results. Interestingly, we found that the property averages of hits selected using GOLD and ICM-Pro fall in

Figure 14. Docking success rates observed in CovDock using Prime Energy and Affinity Score functions. Bright colors refer to success rates in the best scoring docking conformation; lighter colors to rates in Top10 poses.

prospective screens based on the use of covalent docking algorithms. By focusing on the methods evaluated here, only three prospective applications have been reported in the literature to our knowledge (Supplementary Table S6). ICMPro has been used for virtual screening against ubiquitin-like

Table 3. Noncovalent Docking Analysis of Covalent Complexes by Glide and ICM-Pro covalent docking

noncovalent docking on failed covalent

noncovalent docking on whole set

noncovalent docking method

top

successful

failed

successful

failed

successful

failed

Glide SP

1 10 1 10

59% 74% 62% 88%

41% 26% 38% 12%

15% 9% 12% 2%

26% 17% 26% 11%

53% 72% 53% 82%

47% 28% 47% 18%

ICM-Pro

N

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

the other hand, a higher number of noncovalent pharmacophore features improved the performance of FITTED and ICM-Pro only. Our results suggest that covalent docking programs performed differently for various proteins, hence suggesting a target-dependent choice of the docking tool to be preferred. Furthermore, we found that increasing the accessibility of the target cysteine can lead to improved binding mode prediction accuracies. Finally, it was shown for the covalent complexes with failed predictions by the two best performing covalent tools CovDock and ICM-Pro, that noncovalent docking into the Cys/Ala mutated proteins resulted in correct binding conformations in many cases. This observation suggests that sampling and scoring are adversely affected by the presence of the covalent bond in the modeled complex, therefore leaving the usefulness of covalent docking in virtual screening applications up for debate. In summary, the range of performances on this challenging set of complexes suggests that solid success rates of pose predictions can be obtained. Our results on the influencing factors might help the users selecting the optimal docking tool for the particular system under investigation and provide guidance for further improvement of covalent docking methods. These considerations should include the warhead chemistry applied, the size, flexibility and pharmacophore content of the ligand, the accessibility of the targeted cysteine, and the volume of the available binding site. In general, we noticed that neglecting the energetic contributions of the covalent region could lead to accurate pose predictions. This approach puts greater emphasis on the noncovalent interactions and it works for those systems where they represent the main components of the binding. Conversely, in complexes where the covalent bond formation is the key driver (e.g., fragment-sized covalent binders), modeling the optimal geometry of the cysteine-bound electrophile represents the priority for the generation of reasonable binding poses. Another important aspect to consider concerns the conformational changes of the ligand occurring upon covalent linkage. Indeed, ligands characterized by such large structural modifications can be better predicted by tools sampling postreaction conformations in the docking simulation. Altogether, our study not only highlights the key factors influencing the docking performance of the tools investigated but more importantly identifies the most important aspects that should be considered for developing successful covalent docking protocols.

ranges of favorable pose prediction accuracy according to our study. The average number of heavy atoms is 24.0 and 25.0 for validated hits by GOLD and ICM-Pro (Table S6), respectively, both corresponding to ∼60% correct Top1 pose predictions with 2 Å RMSD criterion (Figure 9). The average number of rotatable bonds is 5.5 and 5.3 for validated hits by GOLD and ICM-Pro (Table S6), respectively, both being near to the optimal Top1 pose prediction accuracy with around 70% correct predictions (Figure 10). Furthermore, the remarkable hit rate reported for the virtual screening on CatK with GOLD is in line with the high rate of correct pose predictions we observed against CatK in Top1 with this program (Figure 13). Concerning retrospective virtual screening studies, CovDock-VS has shown high enrichment results43 when used with protein−ligand interaction filtering in retrospective structurebased virtual screening on four different targets such as HCV NS3 protease (EF1% = 52), Cathepsin K (EF1% = 9), EGFR (EF1% = 46), and XPO1 (EF1% = 33). Interestingly, however, without the knowledge-based H-bond filters the enrichment was reduced significantly (almost 2-fold in EF1% and 3-fold in EF10%). Furthermore, considering the binding mode of the 21 cocrystallized known actives CovDock-VS found accurate poses (RMSD below 2.0 Å) in 43% of the cases, as compared to 71% observed in the pose prediction mode. Although with fluctuating hit rates, also other docking-based platforms for virtual screening applications have recently emerged (DOCKovalent,77 DOCKTITE22), highlighting that approaches modeling the covalent bond formation through diverse procedures can lead to the identification of potent and selective covalent compounds.



CONCLUSIONS The recent rise of interest in covalent drugs has propelled the development of new algorithms aimed to predict the binding mode and binding energetics of covalent inhibitors. In this evaluation, we analyzed the performance of six covalent docking tools to reproduce the experimental binding mode of covalent ligands. The analysis was performed on a large and diverse set of 207 covalent complexes collected from the Protein Data Bank representing 54 targets. It was found that around 40−60% of the top scoring ligand poses were predicted to have a RMSD within 2.0 Å from the experimental binding mode. This rate increased to almost 90% for ICM-Pro, to the range of 65−75% for AutoDock4, CovDock, FITTED, and GOLD, and to 50% for MOE when top ten scoring poses are considered. This performance is comparable to that obtained for noncovalent docking tools27 and therefore suggests that anchoring the ligand does not necessarily improve the accuracy of the prediction. It has to be emphasized that for fair comparison these results were obtained using the default parameters of the tools and we omitted any a priori knowledge on the targets. However, it was also illustrated that pharmacophoric or conformational constraints derived from analogous complexes may be able to improve the accuracy of the prediction. We showed that the nature of the electrophilic warhead affects the success rate. Indeed, experimental poses were better reproduced for Michael additions, electrophilic additions and nucleophilic substitutions rather than for ring opening reactions and disulfide formation. The analysis of ligand descriptors has revealed that increasing size and flexibility generally affected the pose prediction unfavorably, although the performances of AutoDock4, FITTED, and ICMPro were found to be less sensitive up to 35 heavy atoms. On



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.8b00228. PDB accession codes, information and calculated RMSD for complexes included in the data set; FITTED success rates on actual set of complexes for which a covalent docking is supported; statistical analysis of RMSD sets; comparison between results obtained via covalent and noncovalent docking approaches; information on prospective covalent docking-based virtual screening applications (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (G.M.K.). O

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling ORCID

(16) Mulliner, D.; Wondrousch, D.; Schuurmann, G. Predicting Michael-Acceptor Reactivity and Toxicity Through Quantum Chemical Transition-State Calculations. Org. Biomol. Chem. 2011, 9, 8400−8412. (17) Krishnan, S.; Miller, R. M.; Tian, B.; Mullins, R. D.; Jacobson, M. P.; Taunton, J. Design of Reversible, Cysteine-Targeted Michael Acceptors Guided by Kinetic and Computational Analysis. J. Am. Chem. Soc. 2014, 136, 12624−12630. (18) Allgauer, D. S.; Jangra, H.; Asahara, H.; Li, Z.; Chen, Q.; Zipse, H.; Ofial, A. R.; Mayr, H. Quantification and Theoretical Analysis of the Electrophilicities of Michael Acceptors. J. Am. Chem. Soc. 2017, 139, 13318−13329. (19) Kumalo, H. M.; Bhakat, S.; Soliman, M. E. S. Theory and Applications of Covalent Docking in Drug Discovery: Merits and Pitfalls. Molecules 2015, 20, 1984−2000. (20) Zhu, K.; Borrelli, K. W.; Greenwood, J. R.; Day, T.; Abel, R.; Farid, R. S.; Harder, E. Docking Covalent Inhibitors: A Parameter Free Approach to Pose Prediction and Scoring. J. Chem. Inf. Model. 2014, 54, 1932−1940. (21) Ouyang, X.; Zhou, S.; Su, C. T. T.; Ge, Z.; Li, R.; Kwoh, C. K. CovalentDock: Automated Covalent Docking with Parameterized Covalent Linkage Energy Estimation and Molecular Geometry Constraints. J. Comput. Chem. 2013, 34, 326−336. (22) Scholz, C.; Knorr, S.; Hamacher, K.; Schmidt, B. DOCKTITEA Highly Versatile Step-by-Step Workflow for Covalent Docking and Virtual Screening in the Molecular Operating Environment. J. Chem. Inf. Model. 2015, 55, 398−406. (23) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. The Protein Data Bank. Nucl. Ac. Res. 2000, 28, 235−242. (24) De Cesco, S.; Kurian, J.; Dufresne, C.; Mittermaier, A.; Moitessier, N. Covalent Inhibitors Design and Discovery. Eur. J. Med. Chem. 2017, 138, 96−114. (25) Miseta, A.; Csutora, P. Relationship Between the Occurrence of Cysteine in Proteins and the Complexity of Organisms. Mol. Biol. Evol. 2000, 17, 1232−1239. (26) Hallenbeck, K. K.; Turner, D. M.; Renslo, A. R.; Arkin, M. R. Targeting Non-Catalytic Cysteine Residues Through StructureGuided Drug Discovery. Curr. Top. Med. Chem. 2016, 17, 4−15. (27) Zhao, Z.; Bourne, P. E. Progress with Covalent Small-Molecule Kinase Inhibitors. Drug Discovery Today 2018, 23, 727−735. (28) Shannon, D. A.; Weerapana, E. Covalent Protein Modification: The Current Landscape of Residue-Specific Electrophiles. Curr. Opin. Chem. Biol. 2015, 24, 18−26. (29) Lagoutte, R.; Patouret, R.; Winssinger, N. Covalent Inhibitors: An Opportunity for Rational Target Selectivity. Curr. Opin. Chem. Biol. 2017, 39, 54−63. (30) Pace, N. J.; Weerapana, E. Diverse Functional Roles of Reactive Cysteines. ACS Chem. Biol. 2013, 8, 283−296. (31) Weerapana, E.; Wang, C.; Simon, G. M.; Richter, F.; Khare, S.; Dillon, M. B.; Bachovchin, D. A.; Mowen, K.; Baker, D.; Cravatt, B. F. Quantitative Reactivity Profiling Predicts Functional Cysteines in Proteomes. Nature 2010, 468, 790−795. (32) Kellenberger, E.; Rodrigo, J.; Muller, P.; Rognan, D. Comparative Evaluation of Eight Docking Tools for Docking and Virtual Screening Accuracy. Proteins: Struct., Funct., Genet. 2004, 57, 225−242. (33) 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.; Tedesco, G.; Wall, I. D.; Woolven, J. M.; Peishoff, C. E.; Head, M. S. A Critical Assessment of Docking Programs and Scoring Functions. J. Med. Chem. 2006, 49, 5912−5931. (34) Wang, Z.; Sun, H.; Yao, X.; Li, D.; Xu, L.; Li, Y.; Tian, S.; Hou, T. Comprehensive Evaluation of Ten Docking Programs on a Diverse Set of Protein. Phys. Chem. Chem. Phys. 2016, 18, 12964−12975. (35) Hartshorn, M. J.; Verdonk, M. L.; Chessari, G.; Brewerton, S. C.; Mooij, W. T. M.; Mortenson, P. N.; Murray, C. W. Diverse, HighQuality Test Set for the Validation of Protein−Ligand Docking Performance. J. Med. Chem. 2007, 50, 726−741.

György M. Keserű : 0000-0003-1039-7809 Funding

This work has been supported by the Marie Sklodowska Curie Action (MSCA) Innovative Training Network grant FRAGNET. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors are grateful to Chris de Graaf, Márton Vass, and Albert Kooistra (Vrije Universiteit Amsterdam) for their assistance in the evaluation of GOLD and MOE software.



REFERENCES

(1) Smith, A. J. T.; Zhang, X.; Leach, A. G.; Houk, K. N. Beyond Picomolar Affinities: Quantitative Aspects of Noncovalent and Covalent Binding of Drugs to Proteins. J. Med. Chem. 2009, 52, 225−233. (2) Singh, J.; Petter, R. C.; Baillie, T. A.; Whitty, A. The Resurgence of Covalent Drugs. Nat. Rev. Drug Discovery 2011, 10, 307−317. (3) Uetrecht, J.; Naisbitt, D. J. Idiosyncratic Adverse Drug Reactions: Current Concepts. Pharmacol. Rev. 2013, 65, 779−808. (4) Lipinski, C.; Hopkins, A. Navigating Chemical Space for Biology and Medicine. Nature 2004, 432, 855−861. (5) Erve, J. C. Chemical Toxicology: Reactive Intermediates and Their Role in Pharmacology and Toxicology. Expert Opin. Drug Metab. Toxicol. 2006, 2, 923−946. (6) Vane, J.; Botting, R. The Mechanism of Action of Aspirin. Thromb. Res. 2003, 110, 255−258. (7) Bauer, R. A. Covalent Inhibitors in Drug Discovery: From Accidental Discoveries to Avoided Liabilities and Designed Therapies. Drug Discov. Drug Discovery Today 2015, 20, 1061−1073. (8) Johnson, D. S.; Weerapana, E.; Cravatt, B. F. Strategies for Discovering and Derisking Covalent, Irreversible Enzyme Inhibitors. Future Med. Chem. 2010, 2, 949−964. (9) Baillie, T. A. Targeted Covalent Inhibitors for Drug Design. Angew. Chem., Int. Ed. 2016, 55, 13408−13421. (10) Lonsdale, R.; Burgess, J.; Colclough, N.; Davies, N. L.; Lenz, E. M.; Orton, A. L.; Ward, R. A. Expanding the Armory: Predicting and Tuning Covalent Warhead Reactivity. J. Chem. Inf. Model. 2017, 57, 3124−3137. (11) Jöst, C.; Nitsche, C.; Scholz, T.; Roux, L.; Klein, C. D. Promiscuity and Selectivity in Covalent Enzyme Inhibition: A Systematic Study of Electrophilic Fragments. J. Med. Chem. 2014, 57, 7590−7599. (12) MacFaul, P. A.; Morley, A. D.; Crawford, J. J. A Simple in Vitro Assay for Assessing the Reactivity of Nitrile Containing Compounds. Bioorg. Med. Chem. Lett. 2009, 19, 1136−1138. (13) Flanagan, M. E.; Abramite, J. A.; Anderson, D. P.; Aulabaugh, A.; Dahal, U. P.; Gilbert, A. M.; Li, C.; Montgomery, J.; Oppenheimer, S. R.; Ryder, T.; Schuff, B. P.; Uccello, D. P.; Walker, G. S.; Wu, Y.; Brown, M. F.; Chen, J. M.; Hayward, M. M.; Noe, M. C.; Obach, R. S.; Philippe, L.; Shanmugasundaram, V.; Shapiro, M. J.; Starr, J.; Stroh, J.; Che, Y. Chemical and Computational Methods for the Characterization of Covalent Reactive Groups for the Prospective Design of Irreversible Inhibitors. J. Med. Chem. 2014, 57, 10072−10079. (14) Cravatt, B. F.; Wright, A. T.; Kozarich, J. W. Activity-Based Protein Profiling: From Enzyme Chemistry to Proteomic Chemistry. Annu. Rev. Biochem. 2008, 77, 383−414. (15) Krenske, E. H.; Petter, R. C.; Zhu, Z.; Houk, K. N. Transition States and Energetics of Nucleophilic Additions of Thiols to Substituted α,β-Unsaturated Ketones: Substituent Effects Involve Enone Stabilization, Product Branching, and Solvation. J. Org. Chem. 2011, 76, 5074−5081. P

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling (36) Doak, B. C.; Over, B.; Giordanetto, F.; Kihlberg, J. Oral Druggable Space Beyond the Rule of 5: Insights from Drugs and Clinical Candidates. Chem. Biol. 2014, 21, 1115−1142. (37) MarvinSketch, Marvin 17.21.0; ChemAxon (https://www. chemaxon.com). (38) Morris, G. M.; Huey, R.; Lindstrom, W.; Sanner, M. F.; Belew, R. K.; Goodsell, D. S.; Olson, A. J. AutoDock4 and AutoDockTools4: Automated Docking with Selective Receptor Flexibility. J. Comput. Chem. 2009, 30, 2785−2791. (39) Schrödinger Suite 2017-4 Protein Preparation Wizard; Epik, Schrödinger, LLC, New York, NY, 2018; Impact, Schrödinger, LLC, New York, NY, 2018; Prime, Schrödinger, LLC, New York, NY, 2017. (40) 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, 221−234. (41) Therrien, E.; Englebienne, P.; Arrowsmith, A. G.; MendozaSanchez, R.; Corbeil, C. R.; Weill, N.; Campagna-Slater, V.; Moitessier, N. Integrating Medicinal Chemistry, Organic/Combinatorial Chemistry, and Computational Chemistry for the Discovery of Selective Estrogen Receptor Modulators with FORECASTER, a Novel Platform for Drug Discovery. J. Chem. Inf. Model. 2012, 52, 210−224. (42) Bianco, G.; Forli, S.; Goodsell, D. S.; Olson, A. J. Covalent Docking Using Autodock: Two-point Attractor and Flexible Side Chain Methods. Protein Sci. 2016, 25, 295−301. (43) Toledo Warshaviak, D.; Golan, G.; Borrelli, K. W.; Zhu, K.; Kalid, O. Structure-Based Virtual Screening Approach for Discovery of Covalently Bound Ligands. J. Chem. Inf. Model. 2014, 54, 1941− 1950. (44) 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.; Shaw, D. E.; Francis, P.; Shenkin, P. S. Glide: a New Approach for Rapid, Accurate Docking and Scoring. 1. Method and Assessment of Docking Accuracy. J. Med. Chem. 2004, 47, 1739− 1749. (45) 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. (46) 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. (47) Watts, K. S.; Dalal, P.; Murphy, R. B.; Sherman, W.; Friesner, R. A.; Shelley, J. C. ConfGen: A Conformational Search Method for Efficient Generation of Bioactive Conformers. J. Chem. Inf. Model. 2010, 50, 534−546. (48) Jacobson, M. P.; Pincus, D. L.; Rapp, C. S.; Day, T. J. F.; Honig, B.; Shaw, D. E.; Friesner, R. A. A Hierarchical Approach to All-Atom Protein Loop Prediction. Proteins: Struct., Funct., Genet. 2004, 55, 351−367. (49) Jacobson, M. P.; Friesner, R. A.; Xiang, Z.; Honig, B. On the Role of the Crystal Environment in Determining Protein Side-Chain Conformations. J. Mol. Biol. 2002, 320, 597−608. (50) Li, J.; Abel, R.; Zhu, K.; Cao, Y.; Zhao, S.; Friesner, R. A. The VSGB 2.0 Model: A Next Generation Energy Model for High Resolution Protein Structure Modeling. Proteins: Struct., Funct., Genet. 2011, 79, 2794−2812. (51) Moitessier, N.; Pottel, J.; Therrien, E.; Englebienne, P.; Liu, Z.; Tomberg, A.; Corbeil, C. R. Medicinal Chemistry Projects Requiring Imaginative Structure-Based Drug Design Methods. Acc. Chem. Res. 2016, 49, 1646−1657. (52) Corbeil, C. R.; Englebienne, P.; Moitessier, N. Docking Ligands into Flexible and Solvated Macromolecules. 1. Development and Validation of FITTED 1.0. J. Chem. Inf. Model. 2007, 47, 435−449. (53) De Cesco, S.; Deslandes, S.; Therrien, E.; Levan, D.; Cueto, M.; Schmidt, R.; Cantin, L. D.; Mittermaier, A.; Juillerat-Jeanneret, L.;

Moitessier, N. Virtual Screening and Computational Optimization for the Discovery of Covalent Prolyl Oligopeptidase Inhibitors with Activity in Human Cells. J. Med. Chem. 2012, 55, 6306−6315. (54) Jones, G.; Willett, P.; Glen, R. C.; Leach, A. R.; Taylor, R. Development and Validation of a Genetic Algorithm for Flexible Docking. J. Mol. Biol. 1997, 267, 727−748. (55) Verdonk, M. L.; Cole, J. C.; Hartshorn, M. J.; Murray, C. W.; Taylor, R. D. Improved Protein-Ligand Docking Using GOLD. Proteins: Struct., Funct., Genet. 2003, 52, 609−623. (56) Cole, J. C.; Nissink, J. W. M.; Taylor, R. Protein-Ligand Docking and Virtual Screening with GOLD. In Virtual Screening in Drug Discovery; Shoichet, B., Alvarez, J., Eds.; Taylor & Francis CRC Press: Boca Raton, Florida, USA, 2005; pp 379−416. (57) Abagyan, R. A.; Totrov, M. M.; Kuznetsov, D. A. ICM: A New Method for Protein Modeling and Design: Applications to Docking and Structure Prediction from the Distorted Native Conformation. J. Comput. Chem. 1994, 15, 488−506. (58) Katritch, V.; Byrd, C. M.; Tseitin, V.; Dai, D.; Raush, E.; Totrov, M.; Abagyan, R.; Jordan, R.; Hruby, D. E. Discovery of Small Molecule Inhibitors of Ubiquitin-like Poxvirus Proteinase I7L Using Homology Modeling and Covalent Docking Approaches. J. Comput.Aided Mol. Des. 2007, 21, 549−558. (59) Abagyan, R. A.; Totrov, M. M. Biased Probability Monte Carlo Conformational Searches and Electrostatic Calculations for Peptides and Proteins. J. Mol. Biol. 1994, 235, 983−1002. (60) Molecular Operating Environment (MOE), 2013.08; Chemical Computing Group ULC, Montreal, QC, Canada, 2018. (61) Labute, P. The Generalized Born/Volume Integral Implicit Solvent Model: Estimation of the Free Energy of Hydration Using London Dispersion Instead of Atomic Surface Area. J. Comput. Chem. 2008, 29, 1693−1698. (62) Winter, J. J. G.; Anderson, M.; Blades, K.; Brassington, C.; Breeze, A. L.; Chresta, C.; Embrey, K.; Fairley, G.; Faulder, P.; Finlay, M. R.; Kettle, J. G.; Nowak, T.; Overman, R.; Patel, S. J.; Perkins, P.; Spadola, L.; Tart, J.; Tucker, J. A.; Wrigley, G. Small Molecule Binding Sites on the Ras:SOS Complex Can Be Exploited for Inhibition of Ras Activation. J. Med. Chem. 2015, 58, 2265−2274. (63) Hardegger, L. A.; Kuhn, B.; Spinnler, B.; Anselm, L.; Ecabert, R.; Stihle, M.; Gsell, B.; Thoma, R.; Diez, J.; Benz, J.; Plancher, J.; Hartmann, G.; Isshiki, Y.; Morikami, K.; Shimma, N.; Haap, W.; Banner, D. W.; Diederich, F. Halogen Bonding at the Active Sites of Human Cathepsin L and MEK1 Kinase: Efficient Interactions in Different Environments. ChemMedChem 2011, 6, 2048−2054. (64) Powers, J. P.; Piper, D. E.; Li, Y.; Mayorga, V.; Anzola, J.; Chen, J. M.; Jaen, J. C.; Lee, G.; Liu, J.; Peterson, G.; Tonn, G. R.; Ye, Q.; Walker, N. P. C.; Wang, Z. SAR and Mode of Action of Novel NonNucleoside Inhibitors of Hepatitis C NS5b RNA Polymerase. J. Med. Chem. 2006, 49, 1034−1046. (65) Bender, A. T.; Gardberg, A.; Pereira, A.; Johnson, T.; Wu, Y.; Grenningloh, R.; Head, J.; Morandi, F.; Haselmayer, P.; Liu-Bujalski, L. Bruton’s Tyrosine Kinase Inhibitors to Sequester Y551 and Prevent Phosphorylation Determines Potency for Inhibition of Fc Receptor but not B-Cell Receptor Signaling. Mol. Pharmacol. 2017, 91, 208− 219. (66) Hajduk, P. J.; Greer, J. A Decade of Fragment-Based Drug Design: Strategic Advances and Lessons Learned. Nat. Rev. Drug Discovery 2007, 6, 211−219. (67) Scott, D. E.; Coyne, A. G.; Hudson, S. A.; Abell, C. FragmentBased Approaches in Drug Discovery and Chemical Biology. Biochemistry 2012, 51, 4990−5003. (68) Backus, K. M.; Correia, B. E.; Lum, K. M.; Forli, S.; Horning, B. D.; González-Paez, G. E.; Chatterjee, S.; Lanning, B. R.; Teijaro, J. R.; Olson, A. J.; Wolan, D. W.; Cravatt, B. F. Proteome-Wide Covalent Ligand Discovery in Native Biological Systems. Nature 2016, 534, 570−574. (69) Salam, N. K.; Nuti, R.; Sherman, W. Novel Method for Generating Structure-Based Pharmacophores Using Energetic Analysis. J. Chem. Inf. Model. 2009, 49, 2356−2368. Q

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling (70) Loving, K.; Salam, N. K.; Sherman, W. Energetic Analysis of Fragment Docking and Application to Structure-Based Pharmacophore Hypothesis Generation. J. Comput.-Aided Mol. Des. 2009, 23, 541−554. (71) Cavallo, L.; Kleinjung, J.; Fraternali, F. POPS: a Fast Algorithm for Solvent Accessible Surface Areas at Atomic and Residue Level. Nucl. Ac. Res. 2003, 31, 3364−3366. (72) Schrödinger Release 2017-4: SiteMap, Schrödinger, LLC, New York, NY, 2018. (73) Halgren, T. Identifying and Characterizing Binding Sites and Assessing Druggability. J. Chem. Inf. Model. 2009, 49, 377−389. (74) Ai, Y.; Yu, L.; Tan, X.; Chai, Y.; Liu, S. Discovery of Covalent Ligands via Noncovalent Docking by Dissecting Covalent Docking Based on a Steric-Clashes Alleviating Receptor (SCAR) Strategy. J. Chem. Inf. Model. 2016, 56, 1563−1575. (75) 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. (76) Schröder, J.; Klinger, A.; Oellien, F.; Marhöfer, R. J.; Duszenko, M.; Selzer, P. M. Docking-Based Virtual Screening of Covalently Binding Ligands: An Orthogonal Lead Discovery Approach. J. Med. Chem. 2013, 56, 1478−1490. (77) London, N.; Miller, R. M.; Krishnan, S.; Uchida, K.; Irwin, J. J.; Eidam, O.; Gibold, L.; Cimermančič, P.; Bonnet, R.; Shoichet, B. K.; Taunton, J. Covalent Docking of Large Libraries for the Discovery of Chemical Probes. Nat. Chem. Biol. 2014, 10, 1066−1072.

R

DOI: 10.1021/acs.jcim.8b00228 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX