Docking of small molecules to farnesoid X receptors using AutoDock Vina with the Convex-PL potential: lessons learned from D3R Grand Challenge 2

The 2016 D3R Grand Challenge 2 provided an opportunity to test multiple protein–ligand docking protocols on a set of ligands bound to farnesoid X receptor that has many available experimental structures. We participated in the Stage 1 of the Challenge devoted to the docking pose predictions, with the mean RMSD value of our submission poses of 2.9 Å. Here we present a thorough analysis of our docking predictions made with AutoDock Vina and the Convex-PL rescoring potential by reproducing our submission protocol and running a series of additional molecular docking experiments. We conclude that a correct receptor structure, or more precisely, the structure of the binding pocket, plays the crucial role in the success of our docking studies. We have also noticed the important role of a local ligand geometry, which seems to be not well discussed in literature. We succeed to improve our results up to the mean RMSD value of 2.15–2.33 Å  dependent on the models of the ligands, if docking these to all available homologous receptors. Overall, for docking of ligands of diverse chemical series we suggest to perform docking of each of the ligands to a set of multiple receptors that are homologous to the target.     

Binding affinities of the farnesoid X receptor in the D3R Grand Challenge 2 estimated by free-energy perturbation and docking

We have studied the binding of 102 ligands to the farnesoid X receptor within the D3R Grand Challenge 2016 blind-prediction competition. First, we employed docking with five different docking software and scoring functions. The selected docked poses gave an average root-mean-squared deviation of 4.2 Å. Consensus scoring gave decent results with a Kendall’s τ of 0.26?±?0.06 and a Spearman’s ρ of 0.41?±?0.08. For a subset of 33 ligands, we calculated relative binding free energies with free-energy perturbation. Five transformations between the ligands involved a change of the net charge and we implemented and benchmarked a semi-analytic correction (Rocklin et al., J Chem Phys 139:184103, 2013) for artifacts caused by the periodic boundary conditions and Ewald summation. The results gave a mean absolute deviation of 7.5 kJ/mol compared to the experimental estimates and a correlation coefficient of R ²?=?0.1. These results were among the four best in this competition out of 22 submissions. The charge corrections were significant (7–8 kJ/mol) and always improved the results. By employing 23 intermediate states in the free-energy perturbation, there was a proper overlap between all states and the precision was 0.1–0.7 kJ/mol. However, thermodynamic cycles indicate that the sampling was insufficient in some of the perturbations.     

Workflows and performances in the ranking prediction of 2016 D3R Grand Challenge 2: lessons learned from a collaborative effort

The 2016 D3R Grand Challenge 2 includes both pose and affinity or ranking predictions. This article is focused exclusively on affinity predictions submitted to the D3R challenge from a collaborative effort of the modeling and informatics group. Our submissions include ranking of 102 ligands covering 4 different chemotypes against the FXR ligand binding domain structure, and the relative binding affinity predictions of the two designated free energy subsets of 15 and 18 compounds. Using all the complex structures prepared in the same way allowed us to cover many types of workflows and compare their performances effectively. We evaluated typical workflows used in our daily structure-based design modeling support, which include docking scores, force field-based scores, QM/MM, MMGBSA, MD-MMGBSA, and MacroModel interaction energy estimations. The best performing methods for the two free energy subsets are discussed. Our results suggest that affinity ranking still remains very challenging; that the knowledge of more structural information does not necessarily yield more accurate predictions; and that visual inspection and human intervention are considerably important for ranking. Knowledge of the mode of action and protein flexibility along with visualization tools that depict polar and hydrophobic maps are very useful for visual inspection. QM/MM-based workflows were found to be powerful in affinity ranking and are encouraged to be applied more often. The standardized input and output enable systematic analysis and support methodology development and improvement for high level blinded predictions.     

Binding pose and affinity prediction in the 2016 D3R Grand Challenge 2 using the Wilma-SIE method

The Farnesoid X receptor (FXR) exhibits significant backbone movement in response to the binding of various ligands and can be a challenge for pose prediction algorithms. As part of the D3R Grand Challenge 2, we tested Wilma-SIE, a rigid-protein docking method, on a set of 36 FXR ligands for which the crystal structures had originally been blinded. These ligands covered several classes of compounds. To overcome the rigid protein limitations of the method, we used an ensemble of publicly available structures for FXR from the PDB. The use of the ensemble allowed Wilma-SIE to predict poses with average and median RMSDs of 2.3 and 1.4 Å, respectively. It was quite clear, however, that had we used a single structure for the receptor the success rate would have been much lower. The most successful predictions were obtained on chemical classes for which one or more crystal structures of the receptor bound to a molecule of the same class was available. In the absence of a crystal structure for the class, observing a consensus binding mode for the ligands of the class using one or more receptor structures of other classes seemed to be indicative of a reasonable pose prediction. Affinity prediction proved to be more challenging with generally poor correlation with experimental IC50s (Kendall tau?~?0.3). Even when the 36 crystal structures were used the accuracy of the predicted affinities was not appreciably improved. A possible cause of difficulty is the internal energy strain arising from conformational differences in the receptor across complexes, which may need to be properly estimated and incorporated into the SIE scoring function.     

Binding free energy predictions of farnesoid X receptor (FXR) agonists using a linear interaction energy (LIE) approach with reliability estimation: application to the D3R Grand Challenge 2

Computational protein binding affinity prediction can play an important role in drug research but performing efficient and accurate binding free energy calculations is still challenging. In the context of phase 2 of the Drug Design Data Resource (D3R) Grand Challenge 2 we used our automated eTOX ALLIES approach to apply the (iterative) linear interaction energy (LIE) method and we evaluated its performance in predicting binding affinities for farnesoid X receptor (FXR) agonists. Efficiency was obtained by our pre-calibrated LIE models and molecular dynamics (MD) simulations at the nanosecond scale, while predictive accuracy was obtained for a small subset of compounds. Using our recently introduced reliability estimation metrics, we could classify predictions with higher confidence by featuring an applicability domain (AD) analysis in combination with protein–ligand interaction profiling. The outcomes of and agreement between our AD and interaction-profile analyses to distinguish and rationalize the performance of our predictions highlighted the relevance of sufficiently exploring protein–ligand interactions during training and it demonstrated the possibility to quantitatively and efficiently evaluate if this is achieved by using simulation data only.     

Performance of multiple docking and refinement methods in the pose prediction D3R prospective Grand Challenge 2016

We describe the performance of multiple pose prediction methods for the D3R 2016 Grand Challenge. The pose prediction challenge includes 36 ligands, which represent 4 chemotypes and some miscellaneous structures against the FXR ligand binding domain. In this study we use a mix of fully automated methods as well as human-guided methods with considerations of both the challenge data and publicly available data. The methods include ensemble docking, colony entropy pose prediction, target selection by molecular similarity, molecular dynamics guided pose refinement, and pose selection by visual inspection. We evaluated the success of our predictions by method, chemotype, and relevance of publicly available data. For the overall data set, ensemble docking, visual inspection, and molecular dynamics guided pose prediction performed the best with overall mean RMSDs of 2.4, 2.2, and 2.2 Å respectively. For several individual challenge molecules, the best performing method is evaluated in light of that particular ligand. We also describe the protein, ligand, and public information data preparations that are typical of our binding mode prediction workflow.     

Docking-undocking combination applied to the D3R Grand Challenge 2015

Novel methods for drug discovery are constantly under development and independent exercises to test and validate them for different goals are extremely useful. The drug discovery data resource (D3R) Grand Challenge 2015 offers an excellent opportunity as an external assessment and validation experiment for Computer-Aided Drug Discovery methods. The challenge comprises two protein targets and prediction tests: binding mode and ligand ranking. We have faced both of them with the same strategy: pharmacophore-guided docking followed by dynamic undocking (a new method tested experimentally here) and, where possible, critical assessment of the results based on pre-existing information. In spite of using methods that are qualitative in nature, our results for binding mode and ligand ranking were amongst the best on Hsp90. Results for MAP4K4 were less positive and we track the different performance across systems to the level of previous knowledge about accessible conformational states. We conclude that docking is quite effective if supplemented by dynamic undocking and empirical information (e.g. binding hot spots, productive protein conformations). This setup is well suited for virtual screening, a frequent application that was not explicitly tested in this edition of the D3R Grand Challenge 2015. Protein flexibility remains as the main cause for hard failures.     

Lessons learned in induced fit docking and metadynamics in the Drug Design Data Resource Grand Challenge 2

Two of the major ongoing challenges in computational drug discovery are predicting the binding pose and affinity of a compound to a protein. The Drug Design Data Resource Grand Challenge 2 was developed to address these problems and to drive development of new methods. The challenge provided the 2D structures of compounds for which the organizers help blinded data in the form of 35 X-ray crystal structures and 102 binding affinity measurements and challenged participants to predict the binding pose and affinity of the compounds. We tested a number of pose prediction methods as part of the challenge; we found that docking methods that incorporate protein flexibility (Induced Fit Docking) outperformed methods that treated the protein as rigid. We also found that using binding pose metadynamics, a molecular dynamics based method, to score docked poses provided the best predictions of our methods with an average RMSD of 2.01 Å. We tested both structure-based (e.g. docking) and ligand-based methods (e.g. QSAR) in the affinity prediction portion of the competition. We found that our structure-based methods based on docking with Smina (Spearman ρ?=?0.614), performed slightly better than our ligand-based methods (ρ?=?0.543), and had equivalent performance with the other top methods in the competition. Despite the overall good performance of our methods in comparison to other participants in the challenge, there exists significant room for improvement especially in cases such as these where protein flexibility plays such a large role.     

Predicting the affinity of Farnesoid X Receptor ligands through a hierarchical ranking protocol: a D3R Grand Challenge 2 case study

The Drug Design Data Resource (D3R) Grand Challenges are blind contests organized to assess the state-of-the-art methods accuracy in predicting binding modes and relative binding free energies of experimentally validated ligands for a given target. The second stage of the D3R Grand Challenge 2 (GC2) was focused on ranking 102 compounds according to their predicted affinity for Farnesoid X Receptor. In this task, our workflow was ranked 5th out of the 77 submissions in the structure-based category. Our strategy consisted in (1) a combination of molecular docking using AutoDock 4.2 and manual edition of available structures for binding poses generation using SeeSAR, (2) the use of HYDE scoring for pose selection, and (3) a hierarchical ranking using HYDE and MM/GBSA. In this report, we detail our pose generation and ligands ranking protocols and provide guidelines to be used in a prospective computer aided drug design program.     

CDOCKER and $$\lambda$$-dynamics for prospective prediction in D3R Grand Challenge 2

The opportunity to prospectively predict ligand bound poses and free energies of binding to the Farnesoid X Receptor in the D3R Grand Challenge 2 provided a useful exercise to evaluate CHARMM based docking (CDOCKER) and \(\lambda\)-dynamics methodologies for use in “real-world” applications in computer aided drug design. In addition to measuring their current performance, several recent methodological developments have been analyzed retrospectively to highlight best procedural practices in future applications. For pose prediction with CDOCKER, when the protein structure used for rigid receptor docking was close to the crystallographic holo structure, reliable poses were obtained. Benzimidazoles, with a known holo receptor structure, were successfully docked with an average RMSD of 0.97 \(\AA\). Other non-benzimidazole ligands displayed less accuracy largely because the receptor structures we chose for docking were too different from the experimental holo structures. However, retrospective analysis has shown that when these ligands were re-docked into their holo structures, the average RMSD dropped to 1.18 \(\AA\) for all ligands. When sulfonamides and spiros were docked with the apo structure, which agrees more with their holo structure than the structures we chose, five out of six ligands were correctly docked. These docking results emphasize the need for flexible receptor docking approaches. For \(\lambda\)-dynamics techniques, including multisite \(\lambda\)-dynamics (MS\(\lambda\)D), reasonable agreement with experiment was observed for the 33 ligands investigated; root mean square errors of 2.08 and 1.67 kcal/mol were obtained for free energy sets 1 and 2, respectively. Retrospectively, soft-core potentials, adaptive landscape flattening, and biasing potential replica exchange (BP-REX) algorithms were critical to model large substituent perturbations with sufficient precision and within restrictive timeframes, such as was required with participation in Grand Challenge 2. These developments, their associated benefits, and proposed procedures for their use in future applications are discussed.     

Combining self- and cross-docking as benchmark tools: the performance of DockBench in the D3R Grand Challenge 2

Molecular docking is a powerful tool in the field of computer-aided molecular design. In particular, it is the technique of choice for the prediction of a ligand pose within its target binding site. A multitude of docking methods is available nowadays, whose performance may vary depending on the data set. Therefore, some non-trivial choices should be made before starting a docking simulation. In the same framework, the selection of the target structure to use could be challenging, since the number of available experimental structures is increasing. Both issues have been explored within this work. The pose prediction of a pool of 36 compounds provided by D3R Grand Challenge 2 organizers was preceded by a pipeline to choose the best protein/docking-method couple for each blind ligand. An integrated benchmark approach including ligand shape comparison and cross-docking evaluations was implemented inside our DockBench software. The results are encouraging and show that bringing attention to the choice of the docking simulation fundamental components improves the results of the binding mode predictions.     

Molecular docking performance evaluated on the D3R Grand Challenge 2015 drug-like ligand datasets

The D3R Grand Challenge 2015 was focused on two protein targets: Heat Shock Protein 90 (HSP90) and Mitogen-Activated Protein Kinase Kinase Kinase Kinase 4 (MAP4K4). We used a protocol involving a preliminary analysis of the available data in PDB and PubChem BioAssay, and then a docking/scoring step using more computationally demanding parameters that were required to provide more reliable predictions. We could evidence that different docking software and scoring functions can behave differently on individual ligand datasets, and that the flexibility of specific binding site residues is a crucial element to provide good predictions.     

Sampling and refinement protocols for template-based macrocycle docking: 2018 D3R Grand Challenge 4

We describe a new template-based method for docking flexible ligands such as macrocycles to proteins. It combines Monte-Carlo energy minimization on the manifold, a fast manifold search method, with BRIKARD for complex flexible ligand searching, and with the MELD accelerator of Replica-Exchange Molecular Dynamics simulations for atomistic degrees of freedom. Here we test the method in the Drug Design Data Resource blind Grand Challenge competition. This method was among the best performers in the competition, giving sub-angstrom prediction quality for the majority of the targets.     

Deep neural network affinity model for BACE inhibitors in D3R Grand Challenge 4

Drug Design Data Resource (D3R) Grand Challenge 4 (GC4) offered a unique opportunity for designing and testing novel methodology for accurate docking and affinity prediction of ligands in an open and blinded manner. We participated in the beta-secretase 1 (BACE) Subchallenge which is comprised of cross-docking and redocking of 20 macrocyclic ligands to BACE and predicting binding affinity for 154 macrocyclic ligands. For this challenge, we developed machine learning models trained specifically on BACE. We developed a deep neural network (DNN) model that used a combination of both structure and ligand-based features that outperformed simpler machine learning models. According to the results released by D3R, we achieved a Spearman's rank correlation coefficient of 0.43(7) for predicting the affinity of 154 ligands. We describe the formulation of our machine learning strategy in detail. We compared the performance of DNN with linear regression, random forest, and support vector machines using ligand-based, structure-based, and combining both ligand and structure-based features. We compared different structures for our DNN and found that performance was highly dependent on fine optimization of the L2 regularization hyperparameter, alpha. We also developed a novel metric of ligand three-dimensional similarity inspired by crystallographic difference density maps to match ligands without crystal structures to similar ligands with known crystal structures. This report demonstrates that detailed parameterization, careful data training and implementation, and extensive feature analysis are necessary to obtain strong performance with more complex machine learning methods. Post hoc analysis shows that scoring functions based only on ligand features are competitive with those also using structural features. Our DNN approach tied for fifth in predicting BACE-ligand binding affinities.

     

Interaction with specific HSP90 residues as a scoring function: validation in the D3R Grand Challenge 2015

Here is reported the development of a novel scoring function that performs remarkably well at identifying the native binding pose of a subset of HSP90 inhibitors containing aminopyrimidine or resorcinol based scaffolds. This scoring function is called PocketScore, and consists of the interaction energy between a ligand and three residues in the binding pocket: Asp93, Thr184 and a water molecule. We integrated PocketScore into a molecular docking workflow, and used it to participate in the Drug Design Data Resource (D3R) Grand Challenge 2015 (GC2015). PocketScore was able to rank 180 molecules of the GC2015 according to their binding affinity with satisfactory performance. These results indicate that the specific residues considered by PocketScore are determinant to properly model the interaction between HSP90 and its subset of inhibitors containing aminopyrimidine or resorcinol based scaffolds. Moreover, the development of PocketScore aimed at improving docking power while neglecting the prediction of binding affinities, suggesting that accurate identification of native binding poses is a determinant factor for the performance of virtual screens.     

Optimal strategies for virtual screening of induced-fit and flexible target in the 2015 D3R Grand Challenge

Induced fit or protein flexibility can make a given structure less useful for docking and/or scoring. The 2015 Drug Design Data Resource (D3R) Grand Challenge provided a unique opportunity to prospectively test optimal strategies for virtual screening in these type of targets: heat shock protein 90 (HSP90), a protein with multiple ligand-induced binding modes; and mitogen-activated protein kinase kinase kinase kinase 4 (MAP4K4), a kinase with a large flexible pocket. Using previously known co-crystal structures, we tested predictions from methods that keep the receptor structure fixed and used (a) multiple receptor/ligand co-crystals as binding templates for minimization or docking (“close”), (b) methods that align or dock to a single receptor (“cross”), and (c) a hybrid approach that chose from multiple bound ligands as initial templates for minimization to a single receptor (“min-cross”). Pose prediction using our “close” models resulted in average ligand RMSDs of 0.32 and 1.6 Å for HSP90 and MAP4K4, respectively, the most accurate models of the community-wide challenge. On the other hand, affinity ranking using our “cross” methods performed well overall despite the fact that a fixed receptor cannot model ligand-induced structural changes,. In addition, “close” methods that leverage the co-crystals of the different binding modes of HSP90 also predicted the best affinity ranking. Our studies suggest that analysis of changes on the receptor structure upon ligand binding can help select an optimal virtual screening strategy.     

Docking rigid macrocycles using Convex-PL,AutoDock Vina,and RDKit in the D3R Grand Challenge 4

The D3R Grand Challenge 4 provided a brilliant opportunity to test macrocyclic docking protocols on a diverse high-quality experimental data. We participated in both pose and affinity prediction exercises. Overall, we aimed to use an automated structure-based docking pipeline built around a set of tools developed in our team. This exercise again demonstrated a crucial importance of the correct local ligand geometry for the overall success of docking. Starting from the second part of the pose prediction stage, we developed a stable pipeline for sampling macrocycle conformers. This resulted in the subangstrom average precision of our pose predictions. In the affinity prediction exercise we obtained average results. However, we could improve these when using docking poses submitted by the best predictors. Our docking tools including the Convex-PL scoring function are available at https://team.inria.fr/nano-d/software/.