Protein-ligand Binding Free Energy Research Articles

Protein–Ligand Binding Free Energies from Exhaustive Docking

We explore the use of exhaustive docking as an alternative to Monte Carlo and molecular dynamics sampling for the direct integration of the partition function for protein-ligand binding. We enumerate feasible poses for the ligand and calculate the Boltzmann factor contribution of each pose to the partition function. From the partition function, the free energy, enthalpy, and entropy can be derived. All our calculations are done with a continuum solvation model that includes solving the Poisson equation. In contrast to Monte Carlo and molecular dynamics simulations, exhaustive docking avoids (within the limitations of a discrete sampling) the question of "Have we run long enough?" due to its deterministic complete enumeration of states. We tested the method on the T4 lysozyme L99A mutant, which has a nonpolar cavity that can accommodate a number of small molecules. We tested two electrostatic models. Model 1 used a solute dielectric of 2.25 for the complex apoprotein and free ligand and 78.5 for the solvent. Model 2 used a solute dielectric of 2.25 for the complex and apoprotein but 1.0 for the free ligand. For our test set of eight molecules, we obtain a reasonable correlation with a Pearson r(2) = 0.66 using model 1. The trend in binding affinity ranking is generally preserved with a Kendall τ = 0.64 and Spearman ρ = 0.83. With model 2, the correlation is improved with a Pearson r(2) = 0.83, Kendall τ = 0.93, and Spearman ρ = 0.98. This suggests that the energy function and sampling method adequately captured most of the thermodynamics of binding of the nonpolar ligands to T4 lysozyme L99A.

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.

A very short history of structure-based design: how did we get here and where do we need to go?

Beddell et al. [1] published the first structure-based design paper in 1976 on hemoglobin ligands using Kendrew wireframe models.The first protein–ligand docking paper was in 1982 [2]. The PDB contained about 200 protein X-ray crystal structures in 1982, but very few were drug discovery targets. Many medicinal chemists believed that relevant macromolecule crystal structures were not likely to ever become available early enough during the lifetime of a drug discovery project to matter. Many also thought that designing an optimal ligand would be trivially obvious if they were lucky enough to have the high-resolution structure of their target. The initial flurry of enthusiasm in the 80s around docking and de novo design yielded to the realization in the 90s that scoring is the hard problem, not design. It’s relatively easy to design a molecule to optimize any given scoring function; the hard part is recognizing which apparently complementary molecules actually will bind to their target. This is the crux of the still-unsolved scoring problem. Predicting free energies of binding in aqueous solution has proven to be far more difficult than most of us realized. Despite 30 years of effort with many different approaches, many different investigators, and massive amounts of computer time, we still cannot reliably predict relative binding free energies with sufficient accuracy to drive organic synthesis during hit-to-lead optimization. CASP (Critical Assessment of Protein Structure, http:// predictioncenter.org) showed the protein folders in 1994 that true, blind prediction is much harder than retrospective ‘‘prediction’’ [3]. Their field has advanced steadily since then. Our own field just began to learn this recently, thanks to the CCDC (Cambridge Crystallographic Data Centre) and SAMPL (http://sampl.eyesopen.com) ‘‘contests’’, plus other critical analyses (e.g., [4]). Blind predictions for small molecule crystal structures began with the CCDC in 1999 [5]. Anthony Nicholls, Vijay Pande, and others started SAMPL0 in 2007 to predict small molecule solvation free energies. SAMPL3 was held this year and included blind prediction of host–guest complexes and trypsin-fragment binding. The initial CCDC contest results were disappointing, but helped stimulate dramatic progress. The most recent CCDC contests proved that high-resolution, accurate, blind prediction of small molecule crystal structures is possible for some small organic molecules, albeit with heroic amounts of computing [6, 7]. This is really exciting work: it proves that current theory and its implementation are sufficient to solve this problem in several non-trivial cases. SAMPL0–3 showed that true blind predictions of solvation free energies, host–guest complexes, protein-fragment structures and relative binding affinities are still remarkably difficult and typically have surprisingly large errors [8–10]. Literature in our field is still cluttered with work claiming to predict what’s already known. It is not even clear where the major errors are (partial charges, dielectric model, fixed versus polarizable force fields, torsion terms, entropy changes due to torsional degrees of freedom, water model, etc.). We still don’t understand why protein–ligand binding free energy doesn’t continue to increase with increasing size of the binding interface [11]. Brute force computing has not helped. Major improvements in our methods are required. The PDB now has over 75,000 structures, of which about 50,000 are protein–ligand complexes (http://www.pdb.org). High-resolution X-ray crystal structures are now routinely J. Blaney (&) Computational Chemistry and Cheminformatics, Small Molecule Drug Discovery, Genentech, 1 DNA Way, South San Francisco, CA 94080, USA e-mail: blaney.jeff@gene.com

Evaluation of protein‐ligand binding free energy focused on its entropic components

The binding free energy for FK506-binding protein-ligand systems is evaluated as a sum of two entropic components, the water-entropy gain, and the configurational-entropy loss for the protein and ligand molecules upon the binding. The two entropic components are calculated using morphometric thermodynamics combined with a statistical-mechanical theory for molecular liquids and the normal mode analysis, respectively. We find that there is an excellent correlation between the calculated and experimental values of the binding free energy. This result is compared with those of several other binding-free energy calculation methods, including MM-PB/SA. The binding can well be elucidated by competition of the two entropic components. Upon the protein-ligand binding, the total volume available to the translational displacement of the coexisting water molecules increases, leading to an increase in the number of accessible configurations of the water. The water-entropy gain, by which the binding is driven, originates primarily from this effect. This study sheds new light on the theoretical prediction of the protein-ligand binding free energy.

Optimized Potential of Mean Force Calculations for Standard Binding Free Energies

The prediction of protein-ligand binding free energies is an important goal of computational biochemistry, yet accuracy, reproducibility, and cost remain a problem. Nevertheless, these are essential requirements for computational methods to become standard binding prediction tools in discovery pipelines. Here, we present the results of an extensive search for an optimal method based on an ensemble of umbrella sampling all-atom molecular simulations tested on the phosphorylated tetrapeptide, pYEEI, binding to the SH2 domain, resulting in an accurate and converged binding free energy of -9.0 ± 0.5 kcal/mol (compared to an experimental value of -8.0 ± 0.1 kcal/mol). We find that a minimum of 300 ns of sampling is required for every prediction, a target easily achievable using new generation accelerated MD codes. Convergence is obtained by using an ensemble of simulations per window, each starting from different initial conformations, and by optimizing window-width, orthogonal restraints, reaction coordinate harmonic potentials, and window-sample time. The use of uncorrelated initial conformations in neighboring windows is important for correctly sampling conformational transitions from the unbound to bound states that affect significantly the precision of the calculations. This methodology thus provides a general recipe for reproducible and practical computations of binding free energies for a class of semirigid protein-ligand systems, within the limit of the accuracy of the force field used.

A water-swap reaction coordinate for the calculation of absolute protein–ligand binding free energies

The accurate prediction of absolute protein-ligand binding free energies is one of the grand challenge problems of computational science. Binding free energy measures the strength of binding between a ligand and a protein, and an algorithm that would allow its accurate prediction would be a powerful tool for rational drug design. Here we present the development of a new method that allows for the absolute binding free energy of a protein-ligand complex to be calculated from first principles, using a single simulation. Our method involves the use of a novel reaction coordinate that swaps a ligand bound to a protein with an equivalent volume of bulk water. This water-swap reaction coordinate is built using an identity constraint, which identifies a cluster of water molecules from bulk water that occupies the same volume as the ligand in the protein active site. A dual topology algorithm is then used to swap the ligand from the active site with the identified water cluster from bulk water. The free energy is then calculated using replica exchange thermodynamic integration. This returns the free energy change of simultaneously transferring the ligand to bulk water, as an equivalent volume of bulk water is transferred back to the protein active site. This, directly, is the absolute binding free energy. It should be noted that while this reaction coordinate models the binding process directly, an accurate force field and sufficient sampling are still required to allow for the binding free energy to be predicted correctly. In this paper we present the details and development of this method, and demonstrate how the potential of mean force along the water-swap coordinate can be improved by calibrating the soft-core Coulomb and Lennard-Jones parameters used for the dual topology calculation. The optimal parameters were applied to calculations of protein-ligand binding free energies of a neuraminidase inhibitor (oseltamivir), with these results compared to experiment. These results demonstrate that the water-swap coordinate provides a viable and potentially powerful new route for the prediction of protein-ligand binding free energies.

Prediction of ligand‐binding sites of proteins by molecular docking calculation for a random ligand library

A new approach to predicting the ligand-binding sites of proteins was developed, using protein-ligand docking computation. In this method, many compounds in a random library are docked onto the whole protein surface. We assumed that the true ligand-binding site would exhibit stronger affinity to the compounds in the random library than the other sites, even if the random library did not include the ligand corresponding to the true binding site. We also assumed that the affinity of the true ligand-binding site would be correlated to the docking scores of the compounds in the random library, if the ligand-binding site was correctly predicted. We call this method the molecular-docking binding-site finding (MolSite) method. The MolSite method was applied to 89 known protein-ligand complex structures extracted from the Protein Data Bank, and it predicted the correct binding sites with about 80-99% accuracy, when only the single top-ranked site was adopted. In addition, the average docking score was weakly correlated to the experimental protein-ligand binding free energy, with a correlation coefficient of 0.44.

A fast empirical GAFF compatible partial atomic charge assignment scheme for modeling interactions of small molecules with biomolecular targets

We report here a new and fast approach [Transferable Partial Atomic Charge Model (TPACM4)-upto four bonds] for deriving the partial atomic charges of small molecules for use in protein/DNA-ligand docking and scoring. We have created a look-up table of 5302 atom types to cover the chemical space of C, H, O, N, S, P, F, Cl, and Br atoms in small molecules together with their quantum mechanical RESP fit charges. The atom types defined span diverse plausible chemical environments of each atom in a molecule. The partial charge on any atom in a given molecule is then assigned by a reference to the look-up table. We tested the sensitivity of the TPACM4 partial charges in estimates of hydrogen bond dimers energies, solvation free energies and protein-ligand binding free energies. An average error ±1.11 kcal/mol and a correlation coefficient of 0.90 is obtained in the calculated protein-ligand binding free energies vis-à-vis an RMS error of ±1.02 kcal/mol and a correlation coefficient of 0.92 obtained with RESP fit charges in comparison to experiment. Similar accuracies are realized in predictions of hydrogen bond energies and solvation free energies of small molecules. For a molecule containing 50-55 atoms, the method takes on the order of milliseconds on a single processor machine to assign partial atomic charges. The TPACM4 programme has been web-enabled and made freely accessible at http://www.scfbio-iitd.res.in/software/drugdesign/charge.jsp.

Mixed Quantum Mechanics/Molecular Mechanics Scoring Function To Predict Protein−Ligand Binding Affinity

Computational methods for predicting protein-ligand binding free energy continue to be popular as a potential cost-cutting method in the drug discovery process. However, accurate predictions are often difficult to make as estimates must be made for certain electronic and entropic terms in conventional force field based scoring functions. Mixed quantum mechanics/molecular mechanics (QM/MM) methods allow electronic effects for a small region of the protein to be calculated, treating the remaining atoms as a fixed charge background for the active site. Such a semi-empirical QM/MM scoring function has been implemented in AMBER using DivCon and tested on a set of 23 metalloprotein-ligand complexes, where QM/MM methods provide a particular advantage in the modeling of the metal ion. The binding affinity of this set of proteins can be calculated with an R(2) of 0.64 and a standard deviation of 1.88 kcal/mol without fitting and 0.71 and a standard deviation of 1.69 kcal/mol with fitted weighting of the individual scoring terms. In this study we explore using various methods to calculate terms in the binding free energy equation, including entropy estimates and minimization standards. From these studies we found that using the rotational bond estimate to ligand entropy results in a reasonable R(2) of 0.63 without fitting. We also found that using the ESCF energy of the proteins without minimization resulted in an R(2) of 0.57, when using the rotatable bond entropy estimate.

Test MM-PB/SA on True Conformational Ensembles of Protein−Ligand Complexes

The molecular mechanics Poisson-Boltzmann surface area (MM-PB/SA) method has been popular for computing protein-ligand binding free energies in recent years. All previous evaluations of the MM-PB/SA method are based upon computer-generated conformational ensembles, which may be affected by the defective computational methods used for preparing these conformational ensembles. In an attempt to reach more convincing conclusions, we have evaluated the MM-PB/SA method on a set of 24 diverse protein-ligand complexes, each of which has a set of conformations derived from NMR spectroscopy. Our results indicate that both MM-PB/SA and molecular mechanics generalized Born surface area (MM-GB/SA) are able to produce a modest correlation between their results and the experimentally measured binding free energies on our test set. In particular, both MM-PB/SA and MM-GB/SA produced better results by using a representative structure (R = 0.72-0.79) rather than averaging over the conformational ensemble of each given complex (R = 0.61-0.74). A head-to-head comparison with four selected scoring functions (X-Score, PLP, ChemScore, and DrugScore) on the same test set reveals that MM-PB/SA and MM-GB/SA results are marginally better than those produced by scoring funcitons, supporting the value of the MM-PB/SA method. Nevertheless, scoring functions are still more cost-effective options, especially for high-throughput tasks.

Prediction of Absolute Solvation Free Energies using Molecular Dynamics Free Energy Perturbation and the OPLS Force Field.

The accurate prediction of protein-ligand binding free energies is a primary objective in computer-aided drug design. The solvation free energy of a small molecule provides a surrogate to the desolvation of the ligand in the thermodynamic process of protein-ligand binding. Here, we use explicit solvent molecular dynamics free energy perturbation to predict the absolute solvation free energies of a set of 239 small molecules, spanning diverse chemical functional groups commonly found in drugs and drug-like molecules. We also compare the performance of absolute solvation free energies obtained using the OPLS_2005 force field with two other commonly used small molecule force fields-general AMBER force field (GAFF) with AM1-BCC charges and CHARMm-MSI with CHelpG charges. Using the OPLS_2005 force field, we obtain high correlation with experimental solvation free energies (R(2) = 0.94) and low average unsigned errors for a majority of the functional groups compared to AM1-BCC/GAFF or CHelpG/CHARMm-MSI. However, OPLS_2005 has errors of over 1.3 kcal/mol for certain classes of polar compounds. We show that predictions on these compound classes can be improved by using a semiempirical charge assignment method with an implicit bond charge correction.

Post Processing of Protein-Compound Docking for Fragment-Based Drug Discovery (FBDD): In-Silico Structure-Based Drug Screening and Ligand-Binding Pose Prediction

For fragment-based drug development, both hit (active) compound prediction and docking-pose (protein-ligand complex structure) prediction of the hit compound are important, since chemical modification (fragment linking, fragment evolution) subsequent to the hit discovery must be performed based on the protein-ligand complex structure. However, the naïve protein-compound docking calculation shows poor accuracy in terms of docking-pose prediction. Thus, post-processing of the protein-compound docking is necessary. Recently, several methods for the post-processing of protein-compound docking have been proposed. In FBDD, the compounds are smaller than those for conventional drug screening. This makes it difficult to perform the protein-compound docking calculation. A method to avoid this problem has been reported. Protein-ligand binding free energy estimation is useful to reduce the procedures involved in the chemical modification of the hit fragment. Several prediction methods have been proposed for high-accuracy estimation of protein-ligand binding free energy. This paper summarizes the various computational methods proposed for docking-pose prediction and their usefulness in FBDD.

Importance of Ligand Reorganization Free Energy in Protein−Ligand Binding-Affinity Prediction

Accurate prediction of the binding affinities of small-molecule ligands to their biological targets is fundamental for structure-based drug design but remains a very challenging task. In this paper, we have performed computational studies to predict the binding models of 31 small-molecule Smac (the second mitochondria-derived activator of caspase) mimetics to their target, the XIAP (X-linked inhibitor of apoptosis) protein, and their binding affinities. Our results showed that computational docking was able to reliably predict the binding models, as confirmed by experimentally determined crystal structures of some Smac mimetics complexed with XIAP. However, all the computational methods we have tested, including an empirical scoring function, two knowledge-based scoring functions, and MM-GBSA (molecular mechanics and generalized Born surface area), yield poor to modest prediction for binding affinities. The linear correlation coefficient (r(2)) value between the predicted affinities and the experimentally determined affinities was found to be between 0.21 and 0.36. Inclusion of ensemble protein-ligand conformations obtained from molecular dynamic simulations did not significantly improve the prediction. However, major improvement was achieved when the free-energy change for ligands between their free- and bound-states, or "ligand-reorganization free energy", was included in the MM-GBSA calculation, and the r(2) value increased from 0.36 to 0.66. The prediction was validated using 10 additional Smac mimetics designed and evaluated by an independent group. This study demonstrates that ligand reorganization free energy plays an important role in the overall binding free energy between Smac mimetics and XIAP. This term should be evaluated for other ligand-protein systems and included in the development of new scoring functions. To our best knowledge, this is the first computational study to demonstrate the importance of ligand reorganization free energy for the prediction of protein-ligand binding free energy.

Calculation of absolute protein-ligand binding free energy using distributed replica sampling

Distributed replica sampling [T. Rodinger et al., J. Chem. Theory Comput. 2, 725 (2006)] is a simple and general scheme for Boltzmann sampling of conformational space by computer simulation in which multiple replicas of the system undergo a random walk in reaction coordinate or temperature space. Individual replicas are linked through a generalized Hamiltonian containing an extra potential energy term or bias which depends on the distribution of all replicas, thus enforcing the desired sampling distribution along the coordinate or parameter of interest regardless of free energy barriers. In contrast to replica exchange methods, efficient implementation of the algorithm does not require synchronicity of the individual simulations. The algorithm is inherently suited for large-scale simulations using shared or heterogeneous computing platforms such as a distributed network. In this work, we build on our original algorithm by introducing Boltzmann-weighted jumping, which allows moves of a larger magnitude and thus enhances sampling efficiency along the reaction coordinate. The approach is demonstrated using a realistic and biologically relevant application; we calculate the standard binding free energy of benzene to the L99A mutant of T4 lysozyme. Distributed replica sampling is used in conjunction with thermodynamic integration to compute the potential of mean force for extracting the ligand from protein and solvent along a nonphysical spatial coordinate. Dynamic treatment of the reaction coordinate leads to faster statistical convergence of the potential of mean force than a conventional static coordinate, which suffers from slow transitions on a rugged potential energy surface.

Calculation of protein–ligand binding free energy by using a polarizable potential

The binding of charged ligands benzamidine and diazamidine to trypsin was investigated by using a polarizable potential energy function and explicit-water molecular dynamics simulations. The binding free energies were computed from the difference between the free energies of decoupling the ligand from water and protein environments. Both the absolute and the relative free energies from the perturbation simulations agree with experimental measurements to within 0.5 kcal.mol(-1). Comparison of free-energy components sampled from different thermodynamic paths indicates that electrostatics is the main driving force behind benzamidine recognition of trypsin. The contribution of electronic polarization to binding appears to be crucial. By computing the free-energy contribution caused by the polarization between the ligand and its surroundings, we found that polarization has the opposite effect in dissimilar environments. Although polarization favors ligand solvation in water, it weakens the protein-ligand attraction by screening the electrostatic interaction between trypsin and benzamidine. We also examined the relative binding free energies of a benzamidine analog diazamidine to trypsin. The changes in free energy on benzamidine-diazamidine substitution were tens of kilocalories in both water and trypsin environments; however, the change in the total binding free energy is <2 kcal.mol(-1) because of cancellation, consistent with the experimental results. Overall, our results suggest that the use of a polarizable force field, given adequate sampling, is capable of achieving chemical accuracy in molecular simulations of protein-ligand recognition.

Computation of entropy contribution to protein-ligand binding free energy

The entropy contribution DeltaS to protein-ligand binding free energy is studied for nine protein-lipid complexes. The entropy effect from the loss of the translational/rotational degrees of freedom (DeltaStr) is calculated using the ideal gas approach. The change in the vibrational entropy (DeltaSvib) is calculated using the effective quantum oscillator approach with frequencies derived from the coordinate covariance matrix, so the inharmonic effects are taken into account. The change in the entropy of solvation (DeltaSsolv) is considered using the binomial cell model (developed by the authors) for the hydrophobic effect. The entropy contribution from loss of conformations that are available for the free ligand (DeltaSconf) is also estimated. It is revealed that the negative in view of binding term DeltaStr is only partly compensated by increasing of DeltaSvib, so T(DeltaStr+DeltaSvib+DeltaSconf)<0 for all complexes under investigation, but taking into account DeltaSsolv leads to significantly increased DeltaS. For all complexes except biotin-streptavidin, the results are found to be in reasonable agreement with experimental data.

Continuum Solvation Models in the Linear Interaction Energy Method

The linear interaction energy (LIE) method in combination with two different continuum solvent models has been applied to calculate protein-ligand binding free energies for a set of inhibitors against the malarial aspartic protease plasmepsin II. Ligand-water interaction energies are calculated from both Poisson-Boltzmann (PB) and Generalized Born (GB) continuum models using snapshots from explicit solvent simulations of the ligand and protein-ligand complex. These are compared to explicit solvent calculations, and we find close agreement between the explicit water and PB solvation models. The GB model overestimates the change in solvation energy, and this is caused by consistent underestimation of the effective Born radii in the protein-ligand complex. The explicit solvent LIE calculations and LIE-PB, with our standard parametrization, reproduce absolute experimental binding free energies with an average unsigned error of 0.5 and 0.7 kcal/mol, respectively. The LIE-GB method, however, requires a constant offset to approach the same level of accuracy.

Efficient Generalized Born Models for Monte Carlo Simulations.

The Generalized Born Surface Area theory (GBSA) has become a popular method to model the solvation of biomolecules. While efficient in the context of molecular dynamics simulations, GBSA calculations do not integrate well with Monte Carlo simulations because of the nonlocal nature of the Generalized Born energy. We present a method by which Monte Carlo Generalized Born simulations can be made seven to eight times faster on a protein-ligand binding free energy calculation with little or no loss of accuracy. The method can be employed in any type of Monte Carlo or Hybrid Monte Carlo-molecular dynamics simulation and should prove useful in numerous applications.

Multiple target screening method for robust and accurate in silico ligand screening

We developed a new in silico multiple target screening (MTS) method, based on a multi-receptor versus multi-ligand docking affinity matrixes, and examined its robustness against changes in the scoring system. According to this method, compounds in a database are docked to multiple proteins. The compounds among these proteins that are likely bind to the target protein are selected as the members of the candidate-hit compound group. Then, the compounds in the group are sorted into descending order using the docking score: the first ( n-th) compound is expected to be the most ( n-th) probable hit compound. This method was applied to the analysis of a set of 142 receptors and 142 compounds using a receptor-ligand docking program, Sievgene [Y. Fukunishi, Y. Mikami, H. Nakamura, Similarities among receptor pockets and among compounds: analysis and application to in silico ligand screening, J. Mol. Graphics Modelling, 24 (2005) 34–45], and the results demonstrated that this method achieves a high hit ratio compared to uniform sampling. We prepared two new scores: the Δ G score, designed to reproduce the protein-ligand binding free energy, and the hit-optimized score, designed to maximize the hit ratio of in silico screening. Using the Sievgene docking score, Δ G score and hit-optimized score, the MTS method is more robust than the multiple active-site correction scoring method [G.P.A. Vigers, J.P. Rizzi, Multiple active site corrections for docking and virtual screening, J. Med. Chem., 47 (2004) 80–89].

Calculation of absolute protein–ligand binding free energy from computer simulations

A general methodology for calculating the equilibrium binding constant of a flexible ligand to a protein receptor is formulated on the basis of potentials of mean force. The overall process is decomposed into several stages that can be computed separately: the free ligand in the bulk is first restrained into the conformation it adopts in the bound state, position, and orientation by applying biasing potentials, then it is translated into the binding site, where it is released completely. The conformational restraining potential is based on the root-mean-square deviation of the peptide coordinates relative to its average conformation in the bound complex. Free energy contributions from each stage are calculated by means of free energy perturbation potential of mean force techniques by using appropriate order parameters. The present approach avoids the need to decouple the ligand from its surrounding (bulk solvent and receptor protein) as is traditionally performed in the double-decoupling scheme. It is believed that the present formulation will be particularly useful when the solvation free energy of the ligand is very large. As an application, the equilibrium binding constant of the phosphotyrosine peptide pYEEI to the Src homology 2 domain of human Lck has been calculated. The results are in good agreement with experimental values.