SAMPL4 & DOCK3.7: lessons for automated docking procedures

The Impact of Supervised Learning Methods in Ultralarge High-Throughput Docking

Structure-based virtual screening methods are, nowadays, one of the key pillars of computational drug discovery. In recent years, a series of studies have reported docking-based virtual screening campaigns of large databases ranging from hundreds to thousands of millions compounds, further identifying novel hits after experimental validation. As these larg-scale efforts are not generally accessible, machine learning-based protocols have emerged to accelerate the identification of virtual hits within an ultralarge chemical space, reaching impressive reductions in computational time. Herein, we illustrate the motivation and the problem behind the screening of large databases, providing an overview of key concepts and essential applications of machine learning-accelerated protocols, specifically concerning supervised learning methods. We also discuss where the field stands with these novel developments, highlighting possible insights for future studies.

Binding free energy predictions in host-guest systems using Autodock4. A retrospective analysis on SAMPL6, SAMPL7 and SAMPL8 challenges

We systematically tested the Autodock4 docking program for absolute binding free energy predictions using the host-guest systems from the recent SAMPL6, SAMPL7 and SAMPL8 challenges. We found that Autodock4 behaves surprisingly well, outperforming in many instances expensive molecular dynamics or quantum chemistry techniques, with an extremely favorable benefit-cost ratio. Some interesting features of Autodock4 predictions are revealed, yielding valuable hints on the overall reliability of docking screening campaigns in drug discovery projects.

Computation of host-guest binding free energies with a new quantum mechanics based mining minima algorithm

A new method called QM-VM2 is presented that efficiently combines statistical mechanics with quantum mechanical (QM) energy potentials in order to calculate noncovalent binding free energies of host-guest systems. QM-VM2 efficiently couples the use of semi-empirical QM (SEQM) energies and geometry optimizations with an underlying molecular mechanics (MM) based conformational search, to find low SEQM energy minima, and allows for processing of these minima at higher levels of ab initio QM theory. A progressive geometry optimization scheme is introduced as a means to increase conformational sampling efficiency. The newly implemented QM-VM2 is used to compute the binding free energies of the host molecule cucurbit[7]uril and a set of 15 guest molecules. The results are presented along with comparisons to experimentally determined binding affinities. For the full set of 15 host-guest complexes, which have a range of formal charges from +1 to +3, SEQM-VM2 based binding free energies show poor correlation with experiment, whereas for the ten +1 complexes only, a significant correlation (R2 = 0.8) is achieved. SEQM-VM2 generation of conformers followed by single-point ab initio QM calculations at the dispersion corrected restricted Hartree-Fock-D3(BJ) and TPSS-D3(BJ) levels of theory, as post-processing corrections, yields a reasonable correlation with experiment for the full set of host-guest complexes (R2 = 0.6 and R2 = 0.7, respectively) and an excellent correlation for the +1 formal charge set (R2 = 1.0 and R2 = 0.9, respectively), as long as a sufficiently large basis set (triple-zeta quality) is employed. The importance of the inclusion of configurational entropy, even at the MM level, for the achievement of good correlation with experiment was demonstrated by comparing the calculated ΔE values with experiment and finding a considerably poorer correlation with experiment than for the calculated free energy ΔE - TΔS. For the complete set of host-guest systems with the range of formal charges, it was observed that the deviation of the predicted binding free energy from experiment correlates somewhat with the net charge of the systems. This observation leads to a simple empirical interpolation scheme to improve the linear regression of the full set.

Artificial intelligence in the early stages of drug discovery

Although the use of computational methods within the pharmaceutical industry is well established, there is an urgent need for new approaches that can improve and optimize the pipeline of drug discovery and development. In spite of the fact that there is no unique solution for this need for innovation, there has recently been a strong interest in the use of Artificial Intelligence for this purpose. As a matter of fact, not only there have been major contributions from the scientific community in this respect, but there has also been a growing partnership between the pharmaceutical industry and Artificial Intelligence companies. Beyond these contributions and efforts there is an underlying question, which we intend to discuss in this review: can the intrinsic difficulties within the drug discovery process be overcome with the implementation of Artificial Intelligence? While this is an open question, in this work we will focus on the advantages that these algorithms provide over the traditional methods in the context of early drug discovery.

AutoPH4: An Automated Method for Generating Pharmacophore Models from Protein Binding Pockets

Pharmacophore models are widely used in computational drug discovery (e.g., in the virtual screening of drug molecules) to capture essential information about interactions between ligands and a target protein. Generating pharmacophore models from protein structures is typically a manual process, but there has been growing interest in automated pharmacophore generation methods. Automation makes feasible the processing of large numbers of protein conformations, such as those generated by molecular dynamics (MD) simulations, and thus may help achieve the longstanding goal of incorporating protein flexibility into virtual screening workflows. Here, we present AutoPH4, a new automated method for generating pharmacophore models based on protein structures; we show that a virtual screening workflow incorporating AutoPH4 ranks compounds more accurately than any other pharmacophore-based virtual screening workflow for which results on a public benchmark have been reported. The strong performance of the virtual screening workflow indicates that the AutoPH4 component of the workflow generates high-quality pharmacophores, making AutoPH4 promising for use in future virtual screening workflows as well, such as ones that use conformations generated by MD simulations.

The SAMPL6 SAMPLing challenge: assessing the reliability and efficiency of binding free energy calculations

Approaches for computing small molecule binding free energies based on molecular simulations are now regularly being employed by academic and industry practitioners to study receptor-ligand systems and prioritize the synthesis of small molecules for ligand design. Given the variety of methods and implementations available, it is natural to ask how the convergence rates and final predictions of these methods compare. In this study, we describe the concept and results for the SAMPL6 SAMPLing challenge, the first challenge from the SAMPL series focusing on the assessment of convergence properties and reproducibility of binding free energy methodologies. We provided parameter files, partial charges, and multiple initial geometries for two octa-acid (OA) and one cucurbit[8]uril (CB8) host-guest systems. Participants submitted binding free energy predictions as a function of the number of force and energy evaluations for seven different alchemical and physical-pathway (i.e., potential of mean force and weighted ensemble of trajectories) methodologies implemented with the GROMACS, AMBER, NAMD, or OpenMM simulation engines. To rank the methods, we developed an efficiency statistic based on bias and variance of the free energy estimates. For the two small OA binders, the free energy estimates computed with alchemical and potential of mean force approaches show relatively similar variance and bias as a function of the number of energy/force evaluations, with the attach-pull-release (APR), GROMACS expanded ensemble, and NAMD double decoupling submissions obtaining the greatest efficiency. The differences between the methods increase when analyzing the CB8-quinine system, where both the guest size and correlation times for system dynamics are greater. For this system, nonequilibrium switching (GROMACS/NS-DS/SB) obtained the overall highest efficiency. Surprisingly, the results suggest that specifying force field parameters and partial charges is insufficient to generally ensure reproducibility, and we observe differences between seemingly converged predictions ranging approximately from 0.3 to 1.0 kcal/mol, even with almost identical simulations parameters and system setup (e.g., Lennard-Jones cutoff, ionic composition). Further work will be required to completely identify the exact source of these discrepancies. Among the conclusions emerging from the data, we found that Hamiltonian replica exchange-while displaying very small variance-can be affected by a slowly-decaying bias that depends on the initial population of the replicas, that bidirectional estimators are significantly more efficient than unidirectional estimators for nonequilibrium free energy calculations for systems considered, and that the Berendsen barostat introduces non-negligible artifacts in expanded ensemble simulations.

Structure-inspired design of β-arrestin-biased ligands for aminergic GPCRs

Development of biased ligands targeting G protein-coupled receptors (GPCRs) is a promising approach for current drug discovery. Although structure-based drug design of biased agonists remains challenging even with an abundance of GPCR crystal structures, we present an approach for translating GPCR structural data into β-arrestin-biased ligands for aminergic GPCRs. We identified specific amino acid-ligand contacts at transmembrane helix 5 (TM5) and extracellular loop 2 (EL2) responsible for Gi/o and β-arrestin signaling, respectively, and targeted those residues to develop biased ligands. For these ligands, we found that bias is conserved at other aminergic GPCRs that retain similar residues at TM5 and EL2. Our approach provides a template for generating arrestin-biased ligands by modifying predicted ligand interactions that block TM5 interactions and promote EL2 interactions. This strategy could facilitate the structure-guided design of arrestin-biased ligands at other GPCRs, including polypharmacological biased ligands.

Free Energy of Binding and Mechanism of Interaction for the MEEVD-TPR2A Peptide-Protein Complex

The association between the MEEVD C-terminal peptide from the heat shock protein 90 (Hsp90) and tetratricopeptide repeat A (TPR2A) domain of the heat shock organizing protein (Hop) is a useful prototype to study the fundamental molecular details about the Hop-Hsp90 interaction. We study here the mechanism of binding/unbinding and compute the standard binding free energy and potential of mean force for the association of the MEEVD peptide to the TPR2A domain using the Adaptive Biasing Force (ABF) methodology. We observe conformational changes of the peptide and the protein receptor induced by binding. We measure the binding free energy of -8.4 kcal/mol, which is consistent with experimental estimates. The simulations achieve multiple unbinding and rebinding events along a consistent pathway connecting the binding site to solvent. The MEEVD peptide slowly dissociates disrupting the hydrogen bonds first, then tilting on the side while preserving the interaction with the side chain of residue Asp 5 of the peptide. After this initial displacement, the peptide completely dissociates and moves into the solvent. Rebinding of the MEEVD peptide from the solvent to the receptor binding site occurs slowly through the portal of entry. Unbinding and rebinding go through intermediate states characterized by the peptide interacting with a lateral helix, helix A1, of the receptor with mainly Asp 5, Val 4, and Glu 3 of the peptide. This newly discovered intermediate structure is characterized by numerous contacts with the receptor which lead to complete formation of the bound complex. The structure of the bound complex obtained after rebinding is structurally very similar to the crystal structure of the complex (0.48 Å root-mean square deviation). The residues Asp 5, Val 4, and Glu 3 adopt conformations and intermolecular contacts with excellent structural similarity with the native ones. Finally, the dissociation and reassociation of MEEVD induce hydration/dehydration transitions, which provide insights on the role of desolvation and solvation processes in protein-peptide binding.

Towards a critical evaluation of an empirical and volume-based solvation function for ligand docking

Molecular docking is an important tool for the discovery of new biologically active molecules given that the receptor structure is known. An excellent environment for the development of new methods and improvement of the current methods is being provided by the rapid growth in the number of proteins with known structure. The evaluation of the solvation energies outstands among the challenges for the modeling of the receptor-ligand interactions, especially in the context of molecular docking where a fast, though accurate, evaluation is ought to be achieved. Here we evaluated a variation of the desolvation energy model proposed by Stouten (Stouten P.F.W. et al, Molecular Simulation, 1993, 10: 97–120), or SV model. The SV model showed a linear correlation with experimentally determined solvation energies, as available in the database FreeSolv. However, when used in retrospective docking simulations using the benchmarks DUD, charged-matched DUD and DUD-Enhanced, the SV model resulted in poorer enrichments when compared to a pure force field model with no correction for solvation effects. The data provided here is consistent with other empirical solvation models employed in the context of molecular docking and indicates that a good model to account for solvent effects is still a goal to achieve. On the other hand, despite the inability to improve the enrichment of retrospective simulations, the SV solvation model showed an interesting ability to reduce the number of molecules with net charge -2 and -3 e among the top-scored molecules in a prospective test.

Bayesian Model Averaging for Ensemble-Based Estimates of Solvation-Free Energies

This paper applies the Bayesian Model Averaging statistical ensemble technique to estimate small molecule solvation free energies. There is a wide range of methods available for predicting solvation free energies, ranging from empirical statistical models to ab initio quantum mechanical approaches. Each of these methods is based on a set of conceptual assumptions that can affect predictive accuracy and transferability. Using an iterative statistical process, we have selected and combined solvation energy estimates using an ensemble of 17 diverse methods from the fourth Statistical Assessment of Modeling of Proteins and Ligands (SAMPL) blind prediction study to form a single, aggregated solvation energy estimate. Methods that possess minimal or redundant information are pruned from the ensemble and the evaluation process repeats until aggregate predictive performance can no longer be improved. We show that this process results in a final aggregate estimate that outperforms all individual methods by reducing estimate errors by as much as 91% to 1.2 kcal mol^-1 accuracy. This work provides a new approach for accurate solvation free energy prediction and lays the foundation for future work on aggregate models that can balance computational cost with prediction accuracy.

Binding-affinity predictions of HSP90 in the D3R Grand Challenge 2015 with docking, MM/GBSA, QM/MM, and free-energy simulations

We have estimated the binding affinity of three sets of ligands of the heat-shock protein 90 in the D3R grand challenge blind test competition. We have employed four different methods, based on five different crystal structures: first, we docked the ligands to the proteins with induced-fit docking with the Glide software and calculated binding affinities with three energy functions. Second, the docked structures were minimised in a continuum solvent and binding affinities were calculated with the MM/GBSA method (molecular mechanics combined with generalised Born and solvent-accessible surface area solvation). Third, the docked structures were re-optimised by combined quantum mechanics and molecular mechanics (QM/MM) calculations. Then, interaction energies were calculated with quantum mechanical calculations employing 970-1160 atoms in a continuum solvent, combined with energy corrections for dispersion, zero-point energy and entropy, ligand distortion, ligand solvation, and an increase of the basis set to quadruple-zeta quality. Fourth, relative binding affinities were estimated by free-energy simulations, using the multi-state Bennett acceptance-ratio approach. Unfortunately, the results were varying and rather poor, with only one calculation giving a correlation to the experimental affinities larger than 0.7, and with no consistent difference in the quality of the predictions from the various methods. For one set of ligands, the results could be strongly improved (after experimental data were revealed) if it was recognised that one of the ligands displaced one or two water molecules. For the other two sets, the problem is probably that the ligands bind in different modes than in the crystal structures employed or that the conformation of the ligand-binding site or the whole protein changes.

BEDAM binding free energy predictions for the SAMPL4 octa-acid host challenge

The binding energy distribution analysis method (BEDAM) protocol has been employed as part of the SAMPL4 blind challenge to predict the binding free energies of a set of octa-acid host-guest complexes. The resulting predictions were consistently judged as some of the most accurate predictions in this category of the SAMPL4 challenge in terms of quantitative accuracy and statistical correlation relative to the experimental values, which were not known at the time the predictions were made. The work has been conducted as part of a hands-on graduate class laboratory session. Collectively the students, aided by automated setup and analysis tools, performed the bulk of the calculations and the numerical and structural analysis. The success of the experiment confirms the reliability of the BEDAM methodology and it shows that physics-based atomistic binding free energy estimation models, when properly streamlined and automated, can be successfully employed by non-specialists.

The SAMPL4 host-guest blind prediction challenge: an overview

Prospective validation of methods for computing binding affinities can help assess their predictive power and thus set reasonable expectations for their performance in drug design applications. Supramolecular host-guest systems are excellent model systems for testing such affinity prediction methods, because their small size and limited conformational flexibility, relative to proteins, allows higher throughput and better numerical convergence. The SAMPL4 prediction challenge therefore included a series of host-guest systems, based on two hosts, cucurbit[7]uril and octa-acid. Binding affinities in aqueous solution were measured experimentally for a total of 23 guest molecules. Participants submitted 35 sets of computational predictions for these host-guest systems, based on methods ranging from simple docking, to extensive free energy simulations, to quantum mechanical calculations. Over half of the predictions provided better correlations with experiment than two simple null models, but most methods underperformed the null models in terms of root mean squared error and linear regression slope. Interestingly, the overall performance across all SAMPL4 submissions was similar to that for the prior SAMPL3 host-guest challenge, although the experimentalists took steps to simplify the current challenge. While some methods performed fairly consistently across both hosts, no single approach emerged as consistent top performer, and the nonsystematic nature of the various submissions made it impossible to draw definitive conclusions regarding the best choices of energy models or sampling algorithms. Salt effects emerged as an issue in the calculation of absolute binding affinities of cucurbit[7]uril-guest systems, but were not expected to affect the relative affinities significantly. Useful directions for future rounds of the challenge might involve encouraging participants to carry out some calculations that replicate each others' studies, and to systematically explore parameter options.

Blind prediction of HIV integrase binding from the SAMPL4 challenge

Here, we give an overview of the protein-ligand binding portion of the Statistical Assessment of Modeling of Proteins and Ligands 4 (SAMPL4) challenge, which focused on predicting binding of HIV integrase inhibitors in the catalytic core domain. The challenge encompassed three components--a small "virtual screening" challenge, a binding mode prediction component, and a small affinity prediction component. Here, we give summary results and statistics concerning the performance of all submissions at each of these challenges. Virtual screening was particularly challenging here in part because, in contrast to more typical virtual screening test sets, the inactive compounds were tested because they were thought to be likely binders, so only the very top predictions performed significantly better than random. Pose prediction was also quite challenging, in part because inhibitors in the set bind to three different sites, so even identifying the correct binding site was challenging. Still, the best methods managed low root mean squared deviation predictions in many cases. Here, we give an overview of results, highlight some features of methods which worked particularly well, and refer the interested reader to papers in this issue which describe specific submissions for additional details.

Combining in silico and in cerebro approaches for virtual screening and pose prediction in SAMPL4

The SAMPL challenges provide an ideal opportunity for unbiased evaluation and comparison of different approaches used in computational drug design. During the fourth round of this SAMPL challenge, we participated in the virtual screening and binding pose prediction on inhibitors targeting the HIV-1 integrase enzyme. For virtual screening, we used well known and widely used in silico methods combined with personal in cerebro insights and experience. Regular docking only performed slightly better than random selection, but the performance was significantly improved upon incorporation of additional filters based on pharmacophore queries and electrostatic similarities. The best performance was achieved when logical selection was added. For the pose prediction, we utilized a similar consensus approach that amalgamated the results of the Glide-XP docking with structural knowledge and rescoring. The pose prediction results revealed that docking displayed reasonable performance in predicting the binding poses. However, prediction performance can be improved utilizing scientific experience and rescoring approaches. In both the virtual screening and pose prediction challenges, the top performance was achieved by our approaches. Here we describe the methods and strategies used in our approaches and discuss the rationale of their performances.