Density functional theory calculations on entire proteins for free energies of binding: Application to a model polar binding site

In drug optimization calculations, the molecular mechanics Poisson‐Boltzmann surface area (MM‐PBSA) method can be used to compute free energies of binding of ligands to proteins. The method involves the evaluation of the energy of configurations in an implicit solvent model. One source of errors is the force field used, which can potentially lead to large errors due to the restrictions in accuracy imposed by its empirical nature. To assess the effect of the force field on the calculation of binding energies, in this article we use large‐scale density functional theory (DFT) calculations as an alternative method to evaluate the energies of the configurations in a “QM‐PBSA” approach. Our DFT calculations are performed with a near‐complete basis set and a minimal parameter implicit solvent model, within the self‐consistent calculation, using the ONETEP program on protein–ligand complexes containing more than 2600 atoms. We apply this approach to the T4‐lysozyme double mutant L99A/M102Q protein, which is a well‐studied model of a polar binding site, using a set of eight small aromatic ligands. We observe that there is very good correlation between the MM and QM binding energies in vacuum but less so in the solvent. The relative binding free energies from DFT are more accurate than the ones from the MM calculations, and give markedly better agreement with experiment for six of the eight ligands. Furthermore, in contrast to MM‐PBSA, QM‐PBSA is able to correctly predict a nonbinder. Proteins 2014; 82:3335–3346. © 2014 Wiley Periodicals, Inc.     

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.     

PSI-DOCK: towards highly efficient and accurate flexible ligand docking

We have developed a new docking method, Pose-Sensitive Inclined (PSI)-DOCK, for flexible ligand docking. An improved SCORE function has been developed and used in PSI-DOCK for binding free energy evaluation. The improved SCORE function was able to reproduce the absolute binding free energies of a training set of 200 protein-ligand complexes with a correlation coefficient of 0.788 and a standard error of 8.13 kJ/mol. For ligand binding pose exploration, a unique searching strategy was designed in PSI-DOCK. In the first step, a tabu-enhanced genetic algorithm with a rapid shape-complementary scoring function is used to roughly explore and store potential binding poses of the ligand. Then, these predicted binding poses are optimized and compete against each other by using a genetic algorithm with the accurate SCORE function to determine the binding pose with the lowest docking energy. The PSI-DOCK 1.0 program is highly efficient in identifying the experimental binding pose. For a test dataset of 194 complexes, PSI-DOCK 1.0 achieved a 67% success rate (RMSD < 2.0 A) for only one run and a 74% success rate for 10 runs. PSI-DOCK can also predict the docking binding free energy with high accuracy. For a test set of 64 complexes, the correlation between the experimentally observed binding free energies and the docking binding free energies for 64 complexes is r = 0.777 with a standard deviation of 7.96 kJ/mol. Moreover, compared with other docking methods, PSI-DOCK 1.0 is extremely easy to use and requires minimum docking preparations. There is no requirement for the users to add hydrogen atoms to proteins because all protein hydrogen atoms and the flexibility of the terminal protein atoms are intrinsically taken into account in PSI-DOCK. There is also no requirement for the users to calculate partial atomic charges because PSI-DOCK does not calculate an electrostatic energy term. These features are not only convenient for the users but also help to avoid the influence of different preparation methods.     

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.     

The effect of multiple binding modes on empirical modeling of ligand docking to proteins.

A popular approach to the computational modeling of ligand/receptor interactions is to use an empirical free energy like model with adjustable parameters. Parameters are learned from one set of complexes, then used to predict another set. To improve these empirical methods requires an independent way to study their inherent errors. We introduce a toy model of ligand/receptor binding as a workbench for testing such errors. We study the errors incurred from the two state binding assumption--the assumption that a ligand is either bound in one orientation, or unbound. We find that the two state assumption can cause large errors in free energy predictions, but it does not affect rank order predictions significantly. We show that fitting parameters using data from high affinity ligands can reduce two state errors; so can using more physical models that do not use the two state assumption. We also find that when using two state models to predict free energies, errors are more severe on high affinity ligands than low affinity ligands. And we show that two state errors can be diagnosed by systematically adding new binding modes when predicting free energies: if predictions worsen as the modes are added, then the two state assumption in the fitting step may be at fault.     

Theoretical and Molecular Docking Study of Ketoconazole on Heptakis(2,3,6‐tri‐O‐methyl)‐β‐cyclodextrin as Chiral Selector

A molecular docking study, using molecular mechanics calculations with AutoDock and semi‐empirical PM3 calculations, was used to predict the enantiodiscrimination of heptakis(2,3,6‐tri‐O‐methyl)‐β‐cyclodextrin (TMβCD) and ketoconazole (KTZ) enantiomers. A Density Functional Theory (DFT) single‐point calculation at the level of B3LYP/6‐311G (d,p) was performed for the PM3‐optimized complexes to obtain more accurate binding energy and the electronic structures of the complexes. The difference in energies of the inclusion complexes between the KTZ enantiomers and TMβCD is probably a measure of chiral discrimination, which results in the separation of the enantiomers as observed in the experimental studies. Chirality 28:209–214, 2016. © 2015 Wiley Periodicals, Inc.     

Comparison of end-point continuum-solvation methods for the calculation of protein-ligand binding free energies

We have compared the predictions of ligand‐binding affinities from several methods based on end‐point molecular dynamics simulations and continuum solvation, that is, methods related to MM/PBSA (molecular mechanics combined with Poisson–Boltzmann and surface area solvation). Two continuum‐solvation models were considered, viz., the Poisson–Boltzmann (PB) and generalised Born (GB) approaches. The nonelectrostatic energies were also obtained in two different ways, viz., either from the sum of the bonded, van der Waals, nonpolar solvation energies, and entropy terms (as in MM/PBSA), or from the scaled protein–ligand van der Waals interaction energy (as in the linear interaction energy approach, LIE). Three different approaches to calculate electrostatic energies were tested, viz., the sum of electrostatic interaction energies and polar solvation energies, obtained either from a single simulation of the complex or from three independent simulations of the complex, the free protein, and the free ligand, or the linear‐response approximation (LRA). Moreover, we investigated the effect of scaling the electrostatic interactions by an effective internal dielectric constant of the protein (?_int). All these methods were tested on the binding of seven biotin analogues to avidin and nine 3‐amidinobenzyl‐1H‐indole‐2‐carboxamide inhibitors to factor Xa. For avidin, the best results were obtained with a combination of the LIE nonelectrostatic energies with the MM+GB electrostatic energies from a single simulation, using ?_int = 4. For fXa, standard MM/GBSA, based on one simulation and using ?_int = 4–10 gave the best result. The optimum internal dielectric constant seems to be slightly higher with PB than with GB solvation. © Proteins 2012; © 2012 Wiley Periodicals, Inc.     

Molecular recognition mechanism of FK506 binding protein: an all-electron fragment molecular orbital study

The fragment molecular orbital (FMO) method has enabled electronic structure calculations and geometry optimizations of very large molecules with ab initio quality. We applied the method to four FK506 binding protein (FKBP) complexes (denoted by their PDB codes 1fkb, 1fkf, 1fkg, and 1fki) containing rapamycin, FK506, and two synthetic ligands. The geometries of reduced complex models were optimized at the restricted Hartree–Fock (FMO‐RHF) level using the 3‐21G basis set, and then for a better estimate of binding, the energetics were refined at a higher level of theory (2nd order Møller–Plesset perturbation theory FMO‐MP2 with the 6‐31G* basis set). Thus, obtained binding energies were ?103.9 (?82.0), ?102.2 (?69.2), ?70.1 (?57.7), and ?71.3 (?55.3) kcal/mol for 1fkb, 1fkf, 1fkg, and 1fki, respectively, where the correlation contribution is given in parentheses. The results show that the electron correlation contribution to binding is extremely important, and it accounts for 70–80% of the binding energy. The molecular recognition mechanism of FKBP was analyzed in detail based on the FMO‐pair interactions between protein residues and the ligands. Solvation effects on the protein–ligand binding were estimated using the Poisson–Boltzmann/surface area model. Proteins 2007. © 2007 Wiley‐Liss, Inc.     

A simple and efficient dispersion correction to the Hartree–Fock theory

One of the most challenging problems in computational chemistry and in drug discovery is the accurate prediction of the binding energy between a ligand and a protein receptor. It is well known that the binding energy calculated with the Hartree–Fock molecular orbital theory (HF) lacks the dispersion interaction energy that significantly affects the accuracy of the total binding energy of a large molecular system. We propose a simple and efficient dispersion energy correction to the HF theory (HF-D_tq). The performance of HF-D_tq was compared with those of several recently proposed dispersion corrected density functional theory methods (DFT-Ds) as to the binding energies of 68 small non-covalent complexes. The overall performance of HF-D_tq was found to be nearly equivalent to that of more sophisticated B3LYP-D₃. HF-D_tq will thus be a useful and powerful method for accurately predicting the binding energy between a ligand and a protein, albeit it is a simple correction procedure based on HF.     

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.     

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.     

Different combinations of atomic interactions predict protein‐small molecule and protein‐DNA/RNA affinities with similar accuracy

Interactions between proteins and other molecules play essential roles in all biological processes. Although it is widely held that a protein's ligand specificity is determined primarily by its three‐dimensional structure, the general principles by which structure determines ligand binding remain poorly understood. Here we use statistical analyses of a large number of protein?ligand complexes with associated binding‐affinity measurements to quantitatively characterize how combinations of atomic interactions contribute to ligand affinity. We find that there are significant differences in how atomic interactions determine ligand affinity for proteins that bind small chemical ligands, those that bind DNA/RNA and those that interact with other proteins. Although protein‐small molecule and protein‐DNA/RNA binding affinities can be accurately predicted from structural data, models predicting one type of interaction perform poorly on the others. Additionally, the particular combinations of atomic interactions required to predict binding affinity differed between small‐molecule and DNA/RNA data sets, consistent with the conclusion that the structural bases determining ligand affinity differ among interaction types. In contrast to what we observed for small‐molecule and DNA/RNA interactions, no statistical models were capable of predicting protein?protein affinity with >60% correlation. We demonstrate the potential usefulness of protein‐DNA/RNA binding prediction as a possible tool for high‐throughput virtual screening to guide laboratory investigations, suggesting that quantitative characterization of diverse molecular interactions may have practical applications as well as fundamentally advancing our understanding of how molecular structure translates into function. Proteins 2015; 83:2100–2114. © 2015 The Authors. Proteins: Structure, Function, and Bioinformatics Published by Wiley Periodicals, Inc.     

Parameterization of an effective potential for protein–ligand binding from host–guest affinity data

Force field accuracy is still one of the “stalemates” in biomolecular modeling. Model systems with high quality experimental data are valuable instruments for the validation and improvement of effective potentials. With respect to protein–ligand binding, organic host–guest complexes have long served as models for both experimental and computational studies because of the abundance of binding affinity data available for such systems. Binding affinity data collected for cyclodextrin (CD) inclusion complexes, a popular model for molecular recognition, is potentially a more reliable resource for tuning energy parameters than hydration free energy measurements. Convergence of binding free energy calculations on CD host–guest systems can also be obtained rapidly, thus offering the opportunity to assess the robustness of these parameters. In this work, we demonstrate how implicit solvent parameters can be developed using binding affinity experimental data and the binding energy distribution analysis method (BEDAM) and validated using the Grid Inhomogeneous Solvation Theory analysis. These new solvation parameters were used to study protein–ligand binding in two drug targets against the HIV‐1 virus and improved the agreement between the calculated and the experimental binding affinities. This work illustrates how benchmark sets of high quality experimental binding affinity data and physics‐based binding free energy models can be used to evaluate and optimize force fields for protein–ligand systems. Copyright © 2015 John Wiley & Sons, Ltd.     

Feature selection and classification of protein–protein complexes based on their binding affinities using machine learning approaches

Protein–protein interactions are intrinsic to virtually every cellular process. Predicting the binding affinity of protein–protein complexes is one of the challenging problems in computational and molecular biology. In this work, we related sequence features of protein–protein complexes with their binding affinities using machine learning approaches. We set up a database of 185 protein–protein complexes for which the interacting pairs are heterodimers and their experimental binding affinities are available. On the other hand, we have developed a set of 610 features from the sequences of protein complexes and utilized Ranker search method, which is the combination of Attribute evaluator and Ranker method for selecting specific features. We have analyzed several machine learning algorithms to discriminate protein‐protein complexes into high and low affinity groups based on their K_d values. Our results showed a 10‐fold cross‐validation accuracy of 76.1% with the combination of nine features using support vector machines. Further, we observed accuracy of 83.3% on an independent test set of 30 complexes. We suggest that our method would serve as an effective tool for identifying the interacting partners in protein–protein interaction networks and human–pathogen interactions based on the strength of interactions. Proteins 2014; 82:2088–2096. © 2014 Wiley Periodicals, Inc.     

Novel statistical-thermodynamic methods to predict protein-ligand binding positions using probability distribution functions

We present two novel methods to predict native protein-ligand binding positions. Both methods identify the native binding position as the most probable position corresponding to a maximum of a probability distribution function (PDF) of possible binding positions in a protein active site. Possible binding positions are the origins of clusters composed, on the basis of root-mean square deviations (RMSD), from the multiple ligand positions determined by a docking algorithm. The difference between the methods lies in the ways the PDF is derived. To validate the suggested methods, we compare the averaged RMSD of the predicted ligand docked positions relative to the experimentally determined positions for a set of 135 PDB protein-ligand complexes. We demonstrate that the suggested methods improve docking accuracy by as much as 21-24% in comparison with a method that simply identifies the binding position as the energy top-scored ligand position.     

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.     

An all atom energy based computational protocol for predicting binding affinities of protein-ligand complexes

We report here a computationally fast protocol for predicting binding affinities of non-metallo protein-ligand complexes. The protocol builds in an all atom energy based empirical scoring function comprising electrostatics, van der Waals, hydrophobicity and loss of conformational entropy of protein side chains upon ligand binding. The method is designed to ensure transferability across diverse systems and has been validated on a heterogenous dataset of 161 complexes consisting of 55 unique protein targets. The scoring function trained on a dataset of 61 complexes yielded a correlation of r=0.92 for the predicted binding free energies against the experimental binding affinities. Model validation and parameter analysis studies ensure the predictive ability of the scoring function. When tested on the remaining 100 protein-ligand complexes a correlation of r=0.92 was recovered. The high correlation obtained underscores the potential applicability of the methodology in drug design endeavors. The scoring function has been web enabled at as binding affinity prediction of protein-ligand (BAPPL) server.     

Potential of mean force for protein-protein interaction studies.

Calculating protein-protein interaction energies is crucial for understanding protein-protein associations. On the basis of the methodology of mean-field potential, we have developed an empirical approach to estimate binding free energy for protein-protein interactions. This knowledge-based approach has been used to derive distance-dependent free energies of protein complexes from a nonredundant training set in the Protein Data Bank (PDB), with a careful treatment of homology. We calculate atom pair potentials for 16 pair interactions, which can reflect the importance of hydrophobic interactions and specific hydrogen-bonding interactions. The derived potentials for hydrogen-bonding interactions show a valley of favorable interactions at a distance of approximately 3 A, corresponding to that of an established hydrogen bond. For the test set of 28 protein complexes, the calculated energies have a correlation coefficient of 0.75 compared with experimental binding free energies. The performance of the method in ranking the binding energies of different protein-protein complexes shows that the energy estimation can be applied to value binding free energies for protein-protein associations.     

What determines the van der Waals coefficient beta in the LIE (linear interaction energy) method to estimate binding free energies using molecular dynamics simulations?

Recently a semiempirical method has been proposed by Aqvist et al. to calculate absolute and relative binding free energies. In this method, the absolute binding free energy of a ligand is estimated as deltaGbind = alpha + beta, where Vel(bound) and Vvdw(bound) are the electrostatic and van der Waals interaction energies between the ligand and the solvated protein from an molecular dynamics (MD) trajectory with ligand bound to protein and Vel(free) and Vel(free) and Vvdw(free) are the electrostatic and van der Waals interaction energies between the ligand and the water from an MD trajectory with the ligand in water. A set of values, alpha = 0.5 and beta = 0.16, was found to give results in good agreement with experimental data. Later, however, different optimal values of beta were found in studies of compounds binding to P450cam and avidin. The present work investigates how the optimal value of beta depends on the nature of binding sites for different protein-ligand interactions. By examining seven ligands interacting with five proteins, we have discovered a linear correlation between the value of beta and the weighted non-polar desolvation ratio (WNDR), with a correlation coefficient of 0.96. We have also examined the ability of this correlation to predict optimal values of beta for different ligands binding to a single protein. We studied twelve neutral compounds bound to avidin. In this case, the WNDR approach gave a better estimate of the absolute binding free energies than results obtained using the fixed value of beta found for biotin-avidin. In terms of reproducing the relative binding free energy to biotin, the fixed-beta value gave better results for compounds similar to biotin, but for compounds less similar to biotin, the WNDR approach led to better relative binding free energies.