Predicting absolute ligand binding free energies to a simple model site

A central challenge in structure-based ligand design is the accurate prediction of binding free energies. Here we apply alchemical free energy calculations in explicit solvent to predict ligand binding in a model cavity in T4 lysozyme. Even in this simple site, there are challenges. We made systematic improvements, beginning with single poses from docking, then including multiple poses, additional protein conformational changes, and using an improved charge model. Computed absolute binding free energies had an RMS error of 1.9 kcal/mol relative to previously determined experimental values. In blind prospective tests, the methods correctly discriminated between several true ligands and decoys in a set of putative binders identified by docking. In these prospective tests, the RMS error in predicted binding free energies relative to those subsequently determined experimentally was only 0.6 kcal/mol. X-ray crystal structures of the new ligands bound in the cavity corresponded closely to predictions from the free energy calculations, but sometimes differed from those predicted by docking. Finally, we examined the impact of holding the protein rigid, as in docking, with a view to learning how approximations made in docking affect accuracy and how they may be improved.     

Flexible docking using tabu search and an empirical estimate of binding affinity

This article describes the implementation of a new docking approach. The method uses a Tabu search methodology to dock flexibly ligand molecules into rigid receptor structures. It uses an empirical objective function with a small number of physically based terms derived from fitting experimental binding affinities for crystallographic complexes. This means that docking energies produced by the searching algorithm provide direct estimates of the binding affinities of the ligands. The method has been tested on 50 ligand-receptor complexes for which the experimental binding affinity and binding geometry are known. All water molecules are removed from the structures and ligand molecules are minimized in vacuo before docking. The lowest energy geometry produced by the docking protocol is within 1.5 Å root-mean square of the experimental binding mode for 86% of the complexes. The lowest energies produced by the docking are in fair agreement with the known free energies of binding for the ligands. Proteins 33:367–382, 1998. © 1998 Wiley-Liss, Inc.     

Incorporating replacement free energy of binding‐site waters in molecular docking

Binding‐site water molecules play a crucial role in protein‐ligand recognition, either being displaced upon ligand binding or forming water bridges to stabilize the complex. However, rigorously treating explicit binding‐site waters is challenging in molecular docking, which requires to fully sample ensembles of waters and to consider the free energy cost of replacing waters. Here, we describe a method to incorporate structural and energetic properties of binding‐site waters into molecular docking. We first developed a solvent property analysis (SPA) program to compute the replacement free energies of binding‐site water molecules by post‐processing molecular dynamics trajectories obtained from ligand‐free protein structure simulation in explicit water. Next, we implemented a distance‐dependent scoring term into DOCK scoring function to take account of the water replacement free energy cost upon ligand binding. We assessed this approach in protein targets containing important binding‐site waters, and we demonstrated that our approach is reliable in reproducing the crystal binding geometries of protein‐ligand‐water complexes, as well as moderately improving the ligand docking enrichment performance. In addition, SPA program (free available to academic users upon request) may be applied in identifying hot‐spot binding‐site residues and structure‐based lead optimization. Proteins 2014; 82:1765–1776. © 2014 Wiley Periodicals, Inc.     

SCORE: A New Empirical Method for Estimating the Binding Affinity of a Protein-Ligand Complex

A new method is presented to estimate the binding affinity of a protein-ligand complex with known three-dimensional structure. The method, SCORE, uses an empirical scoring function to describe the binding free energy, which includes terms to account for van der Waals contact, metal-ligand bonding, hydrogen bonding, desolvation effect, and deformation penalty upon the binding process. The coefficients of each term are obtained by multivariate regressional analysis of a diverse training set of 170 protein-ligand complexes. The final scoring function reproduces the binding free energies of the whole training set with a cross-validated deviation of 6.3 kJ/mol. The predictive ability of the function is further tested by a set of 11 endothiapepsin complexes and the internal consistency of the function is demonstrated in a stepwise procedure named Evolutionary Test. A major innovation of this method is the introduction of an atomic binding score which allows the researcher to inspect and optimize the lead compound rationally in a structure-based drug design scheme.     

Flexible relaxation of rigid-body docking solutions

Molecular Dynamics (MD) simulations have been performed on a set of rigid-body docking poses, carried out over 25 protein-protein complexes. The results show that fully flexible relaxation increases the fraction of native contacts (NC) by up to 70% for certain docking poses. The largest increase in the fraction of NC is observed for docking poses where anchor residues are able to sample their bound conformation. For each MD simulation, structural snap-shots were clustered and the centre of each cluster used as the MD-relaxed docking pose. A comparison between two energy-based scoring schemes, the first calculated for the MD-relaxed poses, the second for energy minimized poses, shows that the former are better in ranking complexes with large hydrophobic interfaces. Furthermore, complexes with large interfaces are generally ranked well, regardless of the type of relaxation method chosen, whereas complexes with small hydrophobic interfaces remain difficult to rank. In general, the results indicate that current force-fields are able to correctly describe direct intermolecular interactions between receptor and ligand molecules. However, these force-fields still fail in cases where protein-protein complexes are stabilized by subtle energy contributions.     

Estimating protein-ligand binding free energy: atomic solvation parameters for partition coefficient and solvation free energy calculation

Solvation energy calculation is one of the main difficulties for the estimation of protein-ligand binding free energy and the correct scoring in docking studies. We have developed a new solvation energy estimation method for protein-ligand binding based on atomic solvation parameter (ASP), which has been shown to improve the power of protein-ligand binding free energy predictions. The ASP set, designed to handle both proteins and organic compounds and derived from experimental n-octanol/water partition coefficient (log P) data, contains 100 atom types (united model that treats hydrogen atoms implicitly) or 119 atom types (all-atom model that treats hydrogen atoms explicitly). By using this unified ASP set, an algorithm was developed for solvation energy calculation and was further integrated into a score function for predicting protein-ligand binding affinity. The score function reproduced the absolute binding free energies of a test set of 50 protein-ligand complexes with a standard error of 8.31 kJ/mol. As a byproduct, a conformation-dependent log P calculation algorithm named ASPLOGP was also implemented. The predictive results of ASPLOGP for a test set of 138 compounds were r = 0.968, s = 0.344 for the all-atom model and r = 0.962, s = 0.367 for the united model, which were better than previous conformation-dependent approaches and comparable to fragmental and atom-based methods. ASPLOGP also gave good predictive results for small peptides. The score function based on the ASP model can be applied widely in protein-ligand interaction studies and structure-based drug design.     

Development and validation of an empirical free energy function for calculating protein–protein binding free energy surfaces

In a previous paper, we described a novel empirical free energy function that was used to accurately predict experimental binding free energies for a diverse test set of 31 protein–protein complexes to within ≈ 1.0 kcal. Here, we extend that work and show that an updated version of the function can be used to (1) accurately predict native binding free energies and (2) rank crystallographic, native-like and non-native binding modes in a physically realistic manner. The modified function includes terms designed to capture some of the unfavorable interactions that characterize non-native interfaces. The function was used to calculate one-dimensional binding free energy surfaces for 21 protein complexes. In roughly 90% of the cases tested, the function was used to place native-like and crystallographic binding modes in global free energy minima. Our analysis further suggests that buried hydrogen bonds might provide the key to distinguishing native from non-native interactions. To the best of our knowledge our function is the only one of its kind, a single expression that can be used to accurately calculate native and non-native binding free energies for a large number of proteins. Given the encouraging results presented in this paper, future work will focus on improving the function and applying it to the protein–protein docking problem.     

Discrete molecular dynamics distinguishes nativelike binding poses from decoys in difficult targets

Virtual screening is one of the major tools used in computer-aided drug discovery. In structure-based virtual screening, the scoring function is critical to identifying the correct docking pose and accurately predicting the binding affinities of compounds. However, the performance of existing scoring functions has been shown to be uneven for different targets, and some important drug targets have proven especially challenging. In these targets, scoring functions cannot accurately identify the native or near-native binding pose of the ligand from among decoy poses, which affects both the accuracy of the binding affinity prediction and the ability of virtual screening to identify true binders in chemical libraries. Here, we present an approach to discriminating native poses from decoys in difficult targets for which several scoring functions failed to correctly identify the native pose. Our approach employs Discrete Molecular Dynamics simulations to incorporate protein-ligand dynamics and the entropic effects of binding. We analyze a collection of poses generated by docking and find that the residence time of the ligand in the native and nativelike binding poses is distinctly longer than that in decoy poses. This finding suggests that molecular simulations offer a unique approach to distinguishing the native (or nativelike) binding pose from decoy poses that cannot be distinguished using scoring functions that evaluate static structures. The success of our method emphasizes the importance of protein-ligand dynamics in the accurate determination of the binding pose, an aspect that is not addressed in typical docking and scoring protocols.     

Exploring the primary electron acceptor (QA)-site of the bacterial reaction center from Rhodobacter sphaeroides. Binding mode of vitamin K derivatives.

The functional replacement of the primary ubiquinone (QA) in the photosynthetic reaction center (RC) from Rhodobacter sphaeroides with synthetic vitamin K derivatives has provided a powerful tool to investigate the electron transfer mechanism. To investigate the binding mode of these quinones to the QA binding site we have determined the binding free energy and charge recombination rate from QA(-) to D+ (kAD) of 29 different 1,4-naphthoquinone derivatives with systematically altered structures. The most striking result was that none of the eight tested compounds carrying methyl groups in both positions 5 and 8 of the aromatic ring exhibited functional binding. To understand the binding properties of these quinones on a molecular level, the structures of the reaction center-naphthoquinone complexes were predicted with ligand docking calculations. All protein--ligand structures show hydrogen bonds between the carbonyl oxygens of the quinone and AlaM260 and HisM219 as found for the native ubiquinone-10 in the X-ray structure. The center-to-center distance between the naphthoquinones at QA and the native ubiquinone-10 at QB (the secondary electron acceptor) is essentially the same, compared to the native structure. A detailed analysis of the docking calculations reveals that 5,8-disubstitution prohibits binding due to steric clashes of the 5-methyl group with the backbone atoms of AlaM260 and AlaM249. The experimentally determined binding free energies were reproduced with an rmsd of approximately 4 kJ x mol(-1) in most cases providing a valuable tool for the design of new artificial electron acceptors and inhibitors.     

A knowledge-based forcefield for protein-protein interface design

A distance-dependent knowledge-based potential for protein-protein interactions is derived and tested for application in protein design. Information on residue type specific C(alpha) and C(beta) pair distances is extracted from complex crystal structures in the Protein Data Bank and used in the form of radial distribution functions. The use of only backbone and C(beta) position information allows generation of relative protein-protein orientation poses with minimal sidechain information. Further coarse-graining can be done simply in the same theoretical framework to give potentials for residues of known type interacting with unknown type, as in a one-sided interface design problem. Both interface design via pose generation followed by sidechain repacking and localized protein-protein docking tests are performed on 39 nonredundant antibody-antigen complexes for which crystal structures are available. As reference, Lennard-Jones potentials, unspecific for residue type and biasing toward varying degrees of residue pair separation are used as controls. For interface design, the knowledge-based potentials give the best combination of consistently designable poses, low RMSD to the known structure, and more tightly bound interfaces with no added computational cost. 77% of the poses could be designed to give complexes with negative free energies of binding. Generally, larger interface separation promotes designability, but weakens the binding of the resulting designs. A localized docking test shows that the knowledge-based nature of the potentials improves performance and compares respectably with more sophisticated all-atoms potentials.     

Knowledge-based scoring function to predict protein-ligand interactions

The development and validation of a new knowledge-based scoring function (DrugScore) to describe the binding geometry of ligands in proteins is presented. It discriminates efficiently between well-docked ligand binding modes (root-mean-square deviation <2.0 A with respect to a crystallographically determined reference complex) and those largely deviating from the native structure, e.g. generated by computer docking programs. Structural information is extracted from crystallographically determined protein-ligand complexes using ReLiBase and converted into distance-dependent pair-preferences and solvent-accessible surface (SAS) dependent singlet preferences for protein and ligand atoms. Definition of an appropriate reference state and accounting for inaccuracies inherently present in experimental data is required to achieve good predictive power. The sum of the pair preferences and the singlet preferences is calculated based on the 3D structure of protein-ligand binding modes generated by docking tools. For two test sets of 91 and 68 protein-ligand complexes, taken from the Protein Data Bank (PDB), the calculated score recognizes poses generated by FlexX deviating <2 A from the crystal structure on rank 1 in three quarters of all possible cases. Compared to FlexX, this is a substantial improvement. For ligand geometries generated by DOCK, DrugScore is superior to the "chemical scoring" implemented into this tool, while comparable results are obtained using the "energy scoring" in DOCK. None of the presently known scoring functions achieves comparable power to extract binding modes in agreement with experiment. It is fast to compute, regards implicitly solvation and entropy contributions and produces correctly the geometry of directional interactions. Small deviations in the 3D structure are tolerated and, since only contacts to non-hydrogen atoms are regarded, it is independent from assumptions of protonation states.     

Improved Ligand-Protein Binding Affinity Predictions Using Multiple Binding Modes

Accurate ligand-protein binding affinity prediction, for a set of similar binders, is a major challenge in the lead optimization stage in drug development. In general, docking and scoring functions perform unsatisfactorily in this application. Docking calculations, followed by molecular dynamics simulations and free energy calculations can be applied to improve the predictions. However, for targets with large, flexible binding sites, with no experimentally determined binding modes for a set of ligands, insufficient sampling can decrease the accuracy of the free energy calculations. Cytochrome P450s, a protein family of major importance for drug metabolism, is an example of a challenging target for binding affinity predictions. As a result, the choice of starting structure from the docking solutions becomes crucial. In this study, an iterative scheme is introduced that includes multiple independent molecular dynamics simulations to obtain weighted ensemble averages to be used in the linear interaction energy method. The proposed scheme makes the initial pose selection less crucial for further simulation, as it automatically calculates the relative weights of the various poses. It also properly takes into account the possibility that multiple binding modes contribute similarly to the overall affinity, or of similar compounds occupying very different poses. The method was applied to a set of 12 compounds binding to cytochrome P450 2C9 and it displayed a root mean-square error of 2.9 kJ/mol.     

Lectin-sugar interaction. Calculated versus experimental binding energies.

Although a steadily increasing number of protein--ligand docking experiments have been performed successfully, there are only few studies concerning protein--sugar interactions. In this study, we investigate the interaction of wheat germ agglutinin (WGA) with N-acetylglucosamine and a number of its derivatives and predict the binding free energies using flexible docking techniques. To assess the quality of our predictions, we also determined those binding free energies experimentally in cell-binding studies. The predicted binding site, ligand orientation, and details of the binding mode are in perfect agreement with the known crystal structure of WGA with a sialoglycopeptide. Furthermore, we obtained an excellent linear correlation of our predicted binding free energies with both our own data and experimental data from the literature [Monsigny, M., Roche, A.C., Sene, C., Maget Dana, R. & Delmotte, F. (1980) Eur. J. Biochem. 104, 147-153.]. In both cases, predicted energies were within 1.0 kJ x mol(-1) of the experimental value. These results illustrate the usefulness of docking-based methods for the qualitative and quantitative prediction of protein--carbohydrate interactions. The insights gained from such theoretical studies may be used to complement the results from the still scarce crystal structures.     

Molecular docking study and development of an empirical binding free energy model for phosphodiesterase 4 inhibitors

In the present work, several computational methodologies were combined to develop a model for the prediction of PDE4B inhibitors' activity. The adequacy of applying the ligand docking approach, keeping the enzyme rigid, to the study of a series of PDE4 inhibitors was confirmed by a previous molecular dynamics analysis of the complete enzyme. An exhaustive docking procedure was performed to identify the most probable binding modes of the ligands to the enzyme, including the active site metal ions and the surrounding structural water molecules. The enzyme-inhibitor interaction enthalpies, refined by using the semiempirical molecular orbital approach, were combined with calculated solvation free energies and entropy considerations in an empirical free energy model that enabled the calculation of binding free energies that correlated very well with experimentally derived binding free energies. Our results indicate that both the inclusion of the structural water molecules close to the ions in the binding site and the use of a free energy model with a quadratic dependency on the ligand free energy of solvation are important aspects to be considered for molecular docking investigations involving the PDE4 enzyme family.     

A detailed comparison of current docking and scoring methods on systems of pharmaceutical relevance

A thorough evaluation of some of the most advanced docking and scoring methods currently available is described, and guidelines for the choice of an appropriate protocol for docking and virtual screening are defined. The generation of a large and highly curated test set of pharmaceutically relevant protein-ligand complexes with known binding affinities is described, and three highly regarded docking programs (Glide, GOLD, and ICM) are evaluated on the same set with respect to their ability to reproduce crystallographic binding orientations. Glide correctly identified the crystallographic pose within 2.0 A in 61% of the cases, versus 48% for GOLD and 45% for ICM. In general Glide appears to perform most consistently with respect to diversity of binding sites and ligand flexibility, while the performance of ICM and GOLD is more binding site-dependent and it is significantly poorer when binding is predominantly driven by hydrophobic interactions. The results also show that energy minimization and reranking of the top N poses can be an effective means to overcome some of the limitations of a given docking function. The same docking programs are evaluated in conjunction with three different scoring functions for their ability to discriminate actives from inactives in virtual screening. The evaluation, performed on three different systems (HIV-1 protease, IMPDH, and p38 MAP kinase), confirms that the relative performance of different docking and scoring methods is to some extent binding site-dependent. GlideScore appears to be an effective scoring function for database screening, with consistent performance across several types of binding sites, while ChemScore appears to be most useful in sterically demanding sites since it is more forgiving of repulsive interactions. Energy minimization of docked poses can significantly improve the enrichments in systems with sterically demanding binding sites. Overall Glide appears to be a safe general choice for docking, while the choice of the best scoring tool remains to a larger extent system-dependent and should be evaluated on a case-by-case basis.     

Validation of an automated procedure for the prediction of relative free energies of binding on a set of aldose reductase inhibitors

Among the available methods for predicting free energies of binding of ligands to a protein, the molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) and molecular mechanics generalized Born surface area (MM-GBSA) approaches have been validated for a relatively limited number of targets and compounds in the training set. Here, we report the results of an extensive study on a series of 28 inhibitors of aldose reductase with experimentally determined crystal structures and inhibitory activities, in which we evaluate the ability of MM-PBSA and MM-GBSA methods in predicting binding free energies using a number of different simulation conditions. While none of the methods proved able to predict absolute free energies of binding in quantitative agreement with the experimental values, calculated and experimental free energies of binding were significantly correlated. Comparing the predicted and experimental ΔG of binding, MM-PBSA proved to perform better than MM-GBSA, and within the MM-PBSA methods, the PBSA of Amber performed similarly to Delphi. In particular, significant relationships between experimental and computed free energies of binding were obtained using Amber PBSA and structures minimized with a distance-dependent dielectric function. Importantly, while free energy predictions are usually made on large collections of equilibrated structures sampled during molecular dynamics in water, we have found that a single minimized structure is a reasonable approximation if relative free energies of binding are to be calculated. This finding is particularly relevant, considering that the generation of equilibrated MD ensembles and the subsequent free energy analysis on multiple snapshots is computationally intensive, while the generation and analysis of a single minimized structure of a protein–ligand complex is relatively fast, and therefore suited for high-throughput virtual screening studies. At this aim, we have developed an automated workflow that integrates all the necessary steps required to generate structures and calculate free energies of binding. The procedure is relatively fast and able to screen automatically and iteratively molecules contained in databases and libraries of compounds. Taken altogether, our results suggest that the workflow can be a valuable tool for ligand identification and optimization, being able to automatically and efficiently refine docking poses, which sometimes may not be accurate, and rank the compounds based on more accurate scoring functions.     

Docking study and binding free energy calculation of poly (ADP-ribose) polymerase inhibitors

Recently, the massively parallel computation of absolute binding free energy with a well-equilibrated system (MP-CAFEE) has been developed. The present study aimed to determine whether the MP-CAFEE method is useful for drug discovery research. In the drug discovery process, it is important for computational chemists to predict the binding affinity accurately without detailed structural information for protein / ligand complex. We investigated the absolute binding free energies for Poly (ADP-ribose) polymerase-1 (PARP-1) / inhibitor complexes, using the MP-CAFEE method. Although each docking model was used as an input structure, it was found that the absolute binding free energies calculated by MP-CAFEE are well consistent with the experimental ones. The accuracy of this method is much higher than that using molecular mechanics Poisson-Boltzmann / surface area (MM / PBSA). Although the simulation time is quite extensive, the reliable predictor of binding free energies would be a useful tool for drug discovery projects.     

Structure-activity relationship for extracellular block of K+ channels by tetraalkylammonium ions

External tetraalkylammonium ion binding to potassium channels is studied using microscopic molecular modelling methods and the experimental structure of the KcsA channel. Relative binding free energies of the KcsA complexes with Me4N+, Et4N+, and n-Pr4N+ are calculated with the molecular dynamics free energy perturbation approach together with automated ligand docking. The four-fold symmetry of the entrance cavity formed by the Tyr82 residues is found to provide stronger binding for the D2d than for the S4 conformation of the ligands. In agreement with experiment the Et4N+ blocker shows several kcal/mol better binding than the other tetraalkylammonium ions.     

Docking enzyme-inhibitor complexes using a preference-based free-energy surface

We present a docking scheme that utilizes both a surface complementarity screen as well as an energetic criterion based on surface area burial. Twenty rigid enzyme/inhibitor complexes with known coordinate sets are arbitrarily separated and reassembled to an average all-atom rms (root mean square) deviation of 1.0 Å from the native complexes. Docking is accomplished by a hierarchical search of geometrically compatible triplets of surface normals on each molecule. A pruned tree of possible bound configurations is built up using successive consideration of larger and larger triplets. The best scoring configurations are then passed through a free-energy screen where the lowest energy member is selected as the predicted native state. The free energy approximation is derived from observations of surface burial by atom pairs across the interface of known enzyme/inhibitor complexes. The occurrence of specific atom-atom surface burial, for a set of complexes with well-defined secondary structure both in the bound and unbound states, is parameterized to mimic the free energy of binding. The docking procedure guides the inhibitor into its native state using orientation and distance-dependent functions that reproduce the ideal model of free energies with an average rms deviation of 0.9 kcal/mol. For all systems studied, this docking procedure identifies a single, unique minimum energy configuration that is highly compatible with the native state. © 1996 Wiley-Liss, Inc.     

Structural study of biologically significant ligands with major birch pollen allergen Betv1 by docking and molecular dynamics simulation

The major birch pollen allergen, Betv1 of Betula verrucosa is the main causative agent of birch pollen allergy in humans. Betv1 is capable of binding several physiological ligands including fatty acids, flavones, cytokinins and sterols. Until now, no structural information from crystallography or NMR is available regarding binding mode of any of these ligands into the binding pocket of Betv1. In the present study thirteen ligands have been successfully docked into the hydrophobic cavity of Betv1 and binding free energies of the complexes have been calculated using AutoDock 3.0.5. A linear relationship with correlation coefficient (R2) of 0.6 is obtained between ΔG(b)s values plotted against their corresponding IC50 values. The complex formed between Betv1 and the best docking pose for each ligand has been optimized by molecular dynamics simulation. Here, we describe the ligand binding of Betv1, which provides insight into the biological function of this protein. This knowledge is required for structural alteration or inhibition of some of these ligands in order to modify the allergenic properties of this protein.