A combined treatment of hydration and dynamical effects for the modeling of host–guest binding thermodynamics: the SAMPL5 blinded challenge

As part of the SAMPL5 blinded experiment, we computed the absolute binding free energies of 22 host–guest complexes employing a novel approach based on the BEDAM single-decoupling alchemical free energy protocol with parallel replica exchange conformational sampling and the AGBNP2 implicit solvation model specifically customized to treat the effect of water displacement as modeled by the Hydration Site Analysis method with explicit solvation. Initial predictions were affected by the lack of treatment of ionic charge screening, which is very significant for these highly charged hosts, and resulted in poor relative ranking of negatively versus positively charged guests. Binding free energies obtained with Debye–Hückel treatment of salt effects were in good agreement with experimental measurements. Water displacement effects contributed favorably and very significantly to the observed binding affinities; without it, the modeling predictions would have grossly underestimated binding. The work validates the implicit/explicit solvation approach employed here and it shows that comprehensive physical models can be effective at predicting binding affinities of molecular complexes requiring accurate treatment of conformational dynamics and hydration.     

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.     

Exhaustive docking and solvated interaction energy scoring: lessons learned from the SAMPL4 challenge

We continued prospective assessments of the Wilma–solvated interaction energy (SIE) platform for pose prediction, binding affinity prediction, and virtual screening on the challenging SAMPL4 data sets including the HIV-integrase inhibitor and two host–guest systems. New features of the docking algorithm and scoring function are tested here prospectively for the first time. Wilma–SIE provides good correlations with actual binding affinities over a wide range of binding affinities that includes strong binders as in the case of SAMPL4 host–guest systems. Absolute binding affinities are also reproduced with appropriate training of the scoring function on available data sets or from comparative estimation of the change in target’s vibrational entropy. Even when binding modes are known, SIE predictions lack correlation with experimental affinities within dynamic ranges below 2 kcal/mol as in the case of HIV-integrase ligands, but they correctly signaled the narrowness of the dynamic range. Using a common protein structure for all ligands can reduce the noise, while incorporating a more sophisticated solvation treatment improves absolute predictions. The HIV-integrase virtual screening data set consists of promiscuous weak binders with relatively high flexibility and thus it falls outside of the applicability domain of the Wilma–SIE docking platform. Despite these difficulties, unbiased docking around three known binding sites of the enzyme resulted in over a third of ligands being docked within 2 Å from their actual poses and over half of the ligands docked in the correct site, leading to better-than-random virtual screening results.     

Using physics-based pose predictions and free energy perturbation calculations to predict binding poses and relative binding affinities for FXR ligands in the D3R Grand Challenge 2

Computer-aided drug design has become an integral part of drug discovery and development in the pharmaceutical and biotechnology industry, and is nowadays extensively used in the lead identification and lead optimization phases. The drug design data resource (D3R) organizes challenges against blinded experimental data to prospectively test computational methodologies as an opportunity for improved methods and algorithms to emerge. We participated in Grand Challenge 2 to predict the crystallographic poses of 36 Farnesoid X Receptor (FXR)-bound ligands and the relative binding affinities for two designated subsets of 18 and 15 FXR-bound ligands. Here, we present our methodology for pose and affinity predictions and its evaluation after the release of the experimental data. For predicting the crystallographic poses, we used docking and physics-based pose prediction methods guided by the binding poses of native ligands. For FXR ligands with known chemotypes in the PDB, we accurately predicted their binding modes, while for those with unknown chemotypes the predictions were more challenging. Our group ranked #1st (based on the median RMSD) out of 46 groups, which submitted complete entries for the binding pose prediction challenge. For the relative binding affinity prediction challenge, we performed free energy perturbation (FEP) calculations coupled with molecular dynamics (MD) simulations. FEP/MD calculations displayed a high success rate in identifying compounds with better or worse binding affinity than the reference (parent) compound. Our studies suggest that when ligands with chemical precedent are available in the literature, binding pose predictions using docking and physics-based methods are reliable; however, predictions are challenging for ligands with completely unknown chemotypes. We also show that FEP/MD calculations hold predictive value and can nowadays be used in a high throughput mode in a lead optimization project provided that crystal structures of sufficiently high quality are available.     

D3R Grand Challenge 4: ligand similarity and MM-GBSA-based pose prediction and affinity ranking for BACE-1 inhibitors

The Drug Design Data Resource (D3R) Grand Challenges present an opportunity to assess, in the context of a blind predictive challenge, the accuracy and the limits of tools and methodologies designed to help guide pharmaceutical drug discovery projects. Here, we report the results of our participation in the D3R Grand Challenge 4 (GC4), which focused on predicting the binding poses and affinity ranking for compounds targeting the $$\beta$$-amyloid precursor protein (BACE-1). Our ligand similarity-based protocol using HYBRID (OpenEye Scientific Software) successfully identified poses close to the native binding mode for most of the ligands with less than 2 Å RMSD accuracy. Furthermore, we compared the performance of our HYBRID-based approach to that of AutoDock Vina and DOCK 6 and found that using a reference ligand to guide the docking process is a better strategy for pose prediction and helped HYBRID to perform better here. We also conducted end-point free energy estimates on molecules dynamics based ensembles of protein-ligand complexes using molecular mechanics combined with generalized Born surface area method (MM-GBSA). We found that the binding affinity ranking based on MM-GBSA scores have poor correlation with the experimental values. Finally, the main lessons from our participation in D3R GC4 are: (i) the generation of the macrocyclic conformers is a key step for successful pose prediction, (ii) the protonation states of the BACE-1 binding site should be treated carefully, (iii) the MM-GBSA method could not discriminate well between different predicted binding poses, and (iv) the MM-GBSA method does not perform well at predicting protein–ligand binding affinities here.     

LogP prediction performance with the SMD solvation model and the M06 density functional family for SAMPL6 blind prediction challenge molecules

This work presents a quantum mechanical model for predicting octanol-water partition coefficients of small protein-kinase inhibitor fragments as part of the SAMPL6 LogP Prediction Challenge. The model calculates solvation free energy differences using the M06-2X functional with SMD implicit solvation and the def2-SVP basis set. This model was identified as dqxk4 in the SAMPL6 Challenge and was the third highest performing model in the physical methods category with 0.49 log Root Mean Squared Error (RMSE) for predicting the 11 compounds in SAMPL6 blind prediction set. We also collaboratively investigated the use of empirical models to address model deficiencies for halogenated compounds at minimal additional computational cost. A mixed model consisting of the dqxk4 physical and hdpuj empirical models found improved performance at 0.34 log RMSE on the SAMPL6 dataset. This collaborative mixed model approach shows how empirical models can be leveraged to expediently improve performance in chemical spaces that are difficult for ab initio methods to simulate.     

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.     

Predicting cyclohexane/water distribution coefficients for the SAMPL5 challenge using MOSCED and the SMD solvation model

We present blind predictions using the solubility parameter based method MOSCED submitted for the SAMPL5 challenge on calculating cyclohexane/water distribution coefficients at 298 K. Reference data to parameterize MOSCED was generated with knowledge only of chemical structure by performing solvation free energy calculations using electronic structure calculations in the SMD continuum solvent. To maintain simplicity and use only a single method, we approximate the distribution coefficient with the partition coefficient of the neutral species. Over the final SAMPL5 set of 53 compounds, we achieved an average unsigned error of \(2.2\pm 0.2\) log units (ranking 15 out of 62 entries), the correlation coefficient (R) was \(0.6\pm 0.1\) (ranking 35), and \(72\pm 6\,\%\) of the predictions had the correct sign (ranking 30). While used here to predict cyclohexane/water distribution coefficients at 298 K, MOSCED is broadly applicable, allowing one to predict temperature dependent infinite dilution activity coefficients in any solvent for which parameters exist, and provides a means by which an excess Gibbs free energy model may be parameterized to predict composition dependent phase-equilibrium.     

SAMPL4 & DOCK3.7: lessons for automated docking procedures

The SAMPL4 challenges were used to test current automated methods for solvation energy, virtual screening, pose and affinity prediction of the molecular docking pipeline DOCK 3.7. Additionally, first-order models of binding affinity were proposed as milestones for any method predicting binding affinity. Several important discoveries about the molecular docking software were made during the challenge: (1) Solvation energies of ligands were five-fold worse than any other method used in SAMPL4, including methods that were similarly fast, (2) HIV Integrase is a challenging target, but automated docking on the correct allosteric site performed well in terms of virtual screening and pose prediction (compared to other methods) but affinity prediction, as expected, was very poor, (3) Molecular docking grid sizes can be very important, serious errors were discovered with default settings that have been adjusted for all future work. Overall, lessons from SAMPL4 suggest many changes to molecular docking tools, not just DOCK 3.7, that could improve the state of the art. Future difficulties and projects will be discussed.     

Performance of HADDOCK and a simple contact-based protein–ligand binding affinity predictor in the D3R Grand Challenge 2

We present the performance of HADDOCK, our information-driven docking software, in the second edition of the D3R Grand Challenge. In this blind experiment, participants were requested to predict the structures and binding affinities of complexes between the Farnesoid X nuclear receptor and 102 different ligands. The models obtained in Stage1 with HADDOCK and ligand-specific protocol show an average ligand RMSD of 5.1 Å from the crystal structure. Only 6/35 targets were within 2.5 Å RMSD from the reference, which prompted us to investigate the limiting factors and revise our protocol for Stage2. The choice of the receptor conformation appeared to have the strongest influence on the results. Our Stage2 models were of higher quality (13 out of 35 were within 2.5 Å), with an average RMSD of 4.1 Å. The docking protocol was applied to all 102 ligands to generate poses for binding affinity prediction. We developed a modified version of our contact-based binding affinity predictor PRODIGY, using the number of interatomic contacts classified by their type and the intermolecular electrostatic energy. This simple structure-based binding affinity predictor shows a Kendall’s Tau correlation of 0.37 in ranking the ligands (7th best out of 77 methods, 5th/25 groups). Those results were obtained from the average prediction over the top10 poses, irrespective of their similarity/correctness, underscoring the robustness of our simple predictor. This results in an enrichment factor of 2.5 compared to a random predictor for ranking ligands within the top 25%, making it a promising approach to identify lead compounds in virtual screening.     

Testing the semi-explicit assembly model of aqueous solvation in the SAMPL4 challenge

Here, we test a method, called semi-explicit assembly (SEA), that computes the solvation free energies of molecules in water in the SAMPL4 blind test challenge. SEA was developed with the intention of being as accurate as explicit-solvent models, but much faster to compute. It is accurate because it uses pre-simulations of simple spheres in explicit solvent to obtain structural and thermodynamic quantities, and it is fast because it parses solute free energies into regionally additive quantities. SAMPL4 provided us the opportunity to make new tests of SEA. Our tests here lead us to the following conclusions: (1) The newest version, called Field-SEA, which gives improved predictions for highly charged ions, is shown here to perform as well as the earlier versions (dipolar and quadrupolar SEA) on this broad blind SAMPL4 test set. (2) We find that both the past and present SEA models give solvation free energies that are as accurate as TIP3P. (3) Using a new approach for force field parameter optimization, we developed improved hydroxyl parameters that ensure consistency with neat-solvent dielectric constants, and found that they led to improved solvation free energies for hydroxyl-containing compounds in SAMPL4. We also learned that these hydroxyl parameters are not just fixing solvent exposed oxygens in a general sense, and therefore do not improve predictions for carbonyl or carboxylic-acid groups. Other such functional groups will need their own independent optimizations for potential improvements. Overall, these tests in SAMPL4 indicate that SEA is an accurate, general and fast new approach to computing solvation free energies.     

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.     

D3R Grand Challenge 2: blind prediction of protein–ligand poses,affinity rankings,and relative binding free energies

The Drug Design Data Resource (D3R) ran Grand Challenge 2 (GC2) from September 2016 through February 2017. This challenge was based on a dataset of structures and affinities for the nuclear receptor farnesoid X receptor (FXR), contributed by F. Hoffmann-La Roche. The dataset contained 102 IC50 values, spanning six orders of magnitude, and 36 high-resolution co-crystal structures with representatives of four major ligand classes. Strong global participation was evident, with 49 participants submitting 262 prediction submission packages in total. Procedurally, GC2 mimicked Grand Challenge 2015 (GC2015), with a Stage 1 subchallenge testing ligand pose prediction methods and ranking and scoring methods, and a Stage 2 subchallenge testing only ligand ranking and scoring methods after the release of all blinded co-crystal structures. Two smaller curated sets of 18 and 15 ligands were developed to test alchemical free energy methods. This overview summarizes all aspects of GC2, including the dataset details, challenge procedures, and participant results. We also consider implications for progress in the field, while highlighting methodological areas that merit continued development. Similar to GC2015, the outcome of GC2 underscores the pressing need for methods development in pose prediction, particularly for ligand scaffolds not currently represented in the Protein Data Bank (http://www.pdb.org), and in affinity ranking and scoring of bound ligands.     

Thermodynamic integration to predict host-guest binding affinities

An alchemical free energy method with explicit solvent molecular dynamics simulations was applied as part of the blind prediction contest SAMPL3 to calculate binding free energies for seven guests to an acyclic cucurbit-[n]uril host. The predictions included determination of protonation states for both host and guests, docking pose generation, and binding free energy calculations using thermodynamic integration. We found a root mean square error (RMSE) of 3.6 kcal mol(-1) from the reference experimental results, with an R(2) correlation of 0.51. The agreement with experiment for the largest contributor to this error, guest 6, is improved by 1.7 kcal mol(-1) when a periodicity-induced free energy correction is applied. The corrections for the other ligands were significantly smaller, and altogether the RMSE was reduced by 0.4 kcal mol(-1). We link properties of the host-guest systems during simulation to errors in the computed free energies. Overall, we show that charged host-guest systems studied here, initialized in unconfirmed docking poses, present a challenge to accurate alchemical simulation methods.     

Grand canonical Monte Carlo simulation of ligand-protein binding

A new application of the grand canonical thermodynamics ensemble to compute ligand-protein binding is described. The described method is sufficiently rapid that it is practical to compute ligand-protein binding free energies for a large number of poses over the entire protein surface, thus identifying multiple putative ligand binding sites. In addition, the method computes binding free energies for a large number of poses. The method is demonstrated by the simulation of two protein-ligand systems, thermolysin and T4 lysozyme, for which there is extensive thermodynamic and crystallographic data for the binding of small, rigid ligands. These low-molecular-weight ligands correspond to the molecular fragments used in computational fragment-based drug design. The simulations correctly identified the experimental binding poses and rank ordered the affinities of ligands in each of these systems.     

Binding affinities in the SAMPL3 trypsin and host–guest blind tests estimated with the MM/PBSA and LIE methods

We have estimated affinities for the binding of 34 ligands to trypsin and nine guest molecules to three different hosts in the SAMPL3 blind challenge, using the MM/PBSA, MM/GBSA, LIE, continuum LIE, and Glide score methods. For the trypsin challenge, none of the methods were able to accurately predict the experimental results. For the MM/GB(PB)SA and LIE methods, the rankings were essentially random and the mean absolute deviations were much worse than a null hypothesis giving the same affinity to all ligand. Glide scoring gave a Kendall's τ index better than random, but the ranking is still only mediocre, τ = 0.2. However, the range of affinities is small and most of the pairs of ligands have an experimental affinity difference that is not statistically significant. Removing those pairs improves the ranking metric to 0.4-1.0 for all methods except CLIE. Half of the trypsin ligands were non-binders according to the binding assay. The LIE methods could not separate the inactive ligands from the active ones better than a random guess, whereas MM/GBSA and MM/PBSA were slightly better than random (area under the receiver-operating-characteristic curve, AUC = 0.65-0.68), and Glide scoring was even better (AUC = 0.79). For the first host, MM/GBSA and MM/PBSA reproduce the experimental ranking fairly good, with τ = 0.6 and 0.5, respectively, whereas the Glide scoring was considerably worse, with a τ = 0.4, highlighting that the success of the methods is system-dependent.     

Rapid prediction of solvation free energy. 3. Application to the SAMPL2 challenge

The SAMPL2 hydration free energy blind prediction challenge consisted of a data set of 41 molecules divided into three subsets: explanatory, obscure and investigatory, where experimental hydration free energies were given for the explanatory, withheld for the obscure, and not known for the investigatory molecules. We employed two solvation models for this challenge, a linear interaction energy (LIE) model based on explicit-water molecular dynamics simulations, and the first-shell hydration (FiSH) continuum model previously calibrated to mimic LIE data. On the 23 compounds from the obscure (blind) dataset, the prospectively submitted LIE and FiSH models provided predictions highly correlated with experimental hydration free energy data, with mean-unsigned-errors of 1.69 and 1.71 kcal/mol, respectively. We investigated several parameters that may affect the performance of these models, namely, the solute flexibility for the LIE explicit-solvent model, the solute partial charging method, and the incorporation of the difference in intramolecular energy between gas and solution phases for both models. We extended this analysis to the various chemical classes that can be formed within the SAMPL2 dataset. Our results strengthen previous findings on the excellent accuracy and transferability of the LIE explicit-solvent approach to predict transfer free energies across a wide spectrum of functional classes. Further, the current results on the SAMPL2 test dataset provide additional support for the FiSH continuum model as a fast yet accurate alternative to the LIE explicit-solvent model. Overall, both the LIE explicit-solvent model and the FiSH continuum solvation model show considerable improvement on the SAMPL2 data set over our previous continuum electrostatics-dispersion solvation model used in the SAMPL1 blind challenge.     

Virtual screening of integrase inhibitors by large scale binding free energy calculations: the SAMPL4 challenge

As part of the SAMPL4 blind challenge, filtered AutoDock Vina ligand docking predictions and large scale binding energy distribution analysis method binding free energy calculations have been applied to the virtual screening of a focused library of candidate binders to the LEDGF site of the HIV integrase protein. The computational protocol leveraged docking and high level atomistic models to improve enrichment. The enrichment factor of our blind predictions ranked best among all of the computational submissions, and second best overall. This work represents to our knowledge the first example of the application of an all-atom physics-based binding free energy model to large scale virtual screening. A total of 285 parallel Hamiltonian replica exchange molecular dynamics absolute protein-ligand binding free energy simulations were conducted starting from docked poses. The setup of the simulations was fully automated, calculations were distributed on multiple computing resources and were completed in a 6-weeks period. The accuracy of the docked poses and the inclusion of intramolecular strain and entropic losses in the binding free energy estimates were the major factors behind the success of the method. Lack of sufficient time and computing resources to investigate additional protonation states of the ligands was a major cause of mispredictions. The experiment demonstrated the applicability of binding free energy modeling to improve hit rates in challenging virtual screening of focused ligand libraries during lead optimization.     

Exhaustive search and solvated interaction energy (SIE) for virtual screening and affinity prediction

We carried out a prospective evaluation of the utility of the SIE (solvation interaction energy) scoring function for virtual screening and binding affinity prediction. Since experimental structures of the complexes were not provided, this was an exercise in virtual docking as well. We used our exhaustive docking program, Wilma, to provide high-quality poses that were rescored using SIE to provide binding affinity predictions. We also tested the combination of SIE with our latest solvation model, first shell of hydration (FiSH), which captures some of the discrete properties of water within a continuum model. We achieved good enrichment in virtual screening of fragments against trypsin, with an area under the curve of about 0.7 for the receiver operating characteristic curve. Moreover, the early enrichment performance was quite good with 50% of true actives recovered with a 15% false positive rate in a prospective calculation and with a 3% false positive rate in a retrospective application of SIE with FiSH. Binding affinity predictions for both trypsin and host-guest complexes were generally within 2 kcal/mol of the experimental values. However, the rank ordering of affinities differing by 2 kcal/mol or less was not well predicted. On the other hand, it was encouraging that the incorporation of a more sophisticated solvation model into SIE resulted in better discrimination of true binders from binders. This suggests that the inclusion of proper Physics in our models is a fruitful strategy for improving the reliability of our binding affinity predictions.     

D3R grand challenge 2015: Evaluation of protein–ligand pose and affinity predictions

The Drug Design Data Resource (D3R) ran Grand Challenge 2015 between September 2015 and February 2016. Two targets served as the framework to test community docking and scoring methods: (1) HSP90, donated by AbbVie and the Community Structure Activity Resource (CSAR), and (2) MAP4K4, donated by Genentech. The challenges for both target datasets were conducted in two stages, with the first stage testing pose predictions and the capacity to rank compounds by affinity with minimal structural data; and the second stage testing methods for ranking compounds with knowledge of at least a subset of the ligand–protein poses. An additional sub-challenge provided small groups of chemically similar HSP90 compounds amenable to alchemical calculations of relative binding free energy. Unlike previous blinded Challenges, we did not provide cognate receptors or receptors prepared with hydrogens and likewise did not require a specified crystal structure to be used for pose or affinity prediction in Stage 1. Given the freedom to select from over 200 crystal structures of HSP90 in the PDB, participants employed workflows that tested not only core docking and scoring technologies, but also methods for addressing water-mediated ligand–protein interactions, binding pocket flexibility, and the optimal selection of protein structures for use in docking calculations. Nearly 40 participating groups submitted over 350 prediction sets for Grand Challenge 2015. This overview describes the datasets and the organization of the challenge components, summarizes the results across all submitted predictions, and considers broad conclusions that may be drawn from this collaborative community endeavor.