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.     

Computation of the contribution from the cavity effect to protein-ligand binding free energy

We present results of the investigation of the cavity creation/annihilation effect in view of formation of the protein-ligand (PL) complexes. The protein and ligand were considered as rigid structures. The change of the cavity creation/annihilation free energy DeltaG(cav) was calculated for three PL complexes using the thermodynamic integration procedure with the original algorithm for growing the interaction potential between the cavity and the water molecules. The thermodynamic cycle consists of two stages, annihilation of the cavity of the ligand for the unbound state and its creation at the active site of the protein (bound state). It was revealed that for all complexes under investigation, the values of DeltaG(cav) are negative and favorable for binding. The main contribution to DeltaG(cav) appears due to the annihilation of the cavity of the ligand. All computations were made using the parallel version of CAVE code, elaborated in our preceding work.     

Theory of binless multi-state free energy estimation with applications to protein-ligand binding

The weighted histogram analysis method (WHAM) is routinely used for computing free energies and expectations from multiple ensembles. Existing derivations of WHAM require observations to be discretized into a finite number of bins. Yet, WHAM formulas seem to hold even if the bin sizes are made arbitrarily small. The purpose of this article is to demonstrate both the validity and value of the multi-state Bennet acceptance ratio (MBAR) method seen as a binless extension of WHAM. We discuss two statistical arguments to derive the MBAR equations, in parallel to the self-consistency and maximum likelihood derivations already known for WHAM. We show that the binless method, like WHAM, can be used not only to estimate free energies and equilibrium expectations, but also to estimate equilibrium distributions. We also provide a number of useful results from the statistical literature, including the determination of MBAR estimators by minimization of a convex function. This leads to an approach to the computation of MBAR free energies by optimization algorithms, which can be more effective than existing algorithms. The advantages of MBAR are illustrated numerically for the calculation of absolute protein-ligand binding free energies by alchemical transformations with and without soft-core potentials. We show that binless statistical analysis can accurately treat sparsely distributed interaction energy samples as obtained from unmodified interaction potentials that cannot be properly analyzed using standard binning methods. This suggests that binless multi-state analysis of binding free energy simulations with unmodified potentials offers a straightforward alternative to the use of soft-core potentials for these alchemical transformations.     

Biophysical limits of protein-ligand binding

In classic work, Kuntz et al. (Proc. Nat. Acad. Sci. USA1999, 96, 9997-10002) introduced the concept of ligand efficiency. Though that study focused primarily on drug-like molecules, it also showed that metal binding led to the greatest ligand efficiencies. Here, the physical limits of binding are examined across the wide variety of small molecules in the Binding MOAD database. The complexes with the greatest ligand efficiencies share the trait of being small, charged ligands bound in highly charged, well buried binding sites. The limit of ligand efficiency is -1.75 kcal/mol·atom for the protein-ligand complexes within Binding MOAD, and 95% of the set have efficiencies below a "soft limit" of -0.83 kcal/mol·atom. On the basis of buried molecular surface area, the hard limit of ligand efficiency is -117 cal/mol·?(2), which is in surprising agreement with the limit of macromolecule-protein binding. Close examination of the most efficient systems reveals their incredibly high efficiency is dictated by tight contacts between the charged groups of the ligand and the pocket. In fact, a misfit of 0.24 ? in the average contacts inherently decreases the maximum possible efficiency by at least 0.1 kcal/mol·atom.     

Pairwise decomposition of residue interaction energies using semiempirical quantum mechanical methods in studies of protein-ligand interaction

Pairwise decomposition of the interaction energy between molecules is shown to be a powerful tool that can increase our understanding of macromolecular recognition processes. Herein we calculate the pairwise decomposition of the interaction energy between the protein human carbonic anhydrase II (HCAII) and the fluorine-substituted ligand N-(4-sulfamylbenzoyl)benzylamine (SBB) using semiempirical quantum mechanics based methods. We dissect the interaction between the ligand and the protein by dividing the ligand and the protein into subsystems to understand the structure-activity relationships as a result of fluorine substitution. In particular, the off-diagonal elements of the Fock matrix that is composed of the interaction between the ionic core and the valence electrons and the exchange energy between the subsystems or atoms of interest is examined in detail. Our analysis reveals that the fluorine-substituted benzylamine group of SBB does not directly affect the binding energy. Rather, we find that the strength of the interaction between Thr199 of HCAII and the sulfamylbenzoyl group of SBB affects the binding affinity between the protein and the ligand. These observations underline the importance of the sulfonamide group in binding affinity as shown by previous experiments (Maren, T. H.; Wiley: C. E. J. Med. Chem. 1968, 11, 228-232). Moreover, our calculations qualitatively agree with the structural aspects of these protein-ligand complexes as determined by X-ray crystallography.     

Implementation of chemometric methodology in ACE: predictive investigation of protein-ligand binding

An ACE predictive investigation of protein-ligand binding using a highly effective chemometric response surface design technique is presented. Here, K(d) was estimated using one noninteracting standard which relates to changes in the electrophoretic mobility of carbonic anhydrase B (CAB, EC 4.2.1.1) on complexation with the ligand 4-carboxybenzenesulfonamide (CBSA) present in the electrophoresis buffer. Experimental factors including injection time, capillary length, and applied voltage were selected and tested at three levels in a Box-Behnken design. Statistical analysis results were used to create a mathematical model for response surface prediction via contour and surface plots at a given target response (K(d) = 1.19x10(-6) M). As expected, there were a number of predicted solutions that reached our target response based on the significance of each factor at appropriate levels. The adequacy of the model was validated by experimental runs with the predicted model solution (capillary length = 47 cm, voltage = 11 kV, injection time = 0.01 min) presented in detail as an example.     

The influence of variations of ligand protonation and tautomerism on protein-ligand recognition and binding energy landscape

We investigate the influence of variations of ligand protonation and tautomeric states on the protein-ligand binding energy landscape by applying the concept of structural consensus. In docking simulations, allowing full flexibility of the ligand, we explore whether the native binding mode could be successfully recovered using a non-native ligand protonation state. Here, we consider three proteins, dihydrofolate reductase, transketolase, and alpha-trichosanthin, complexed with ligands having multiple tautomeric forms. We find that for the majority of protonation and tautomeric states the native binding mode can be recovered without a great loss of accuracy.     

Sampling and convergence in free energy calculations of protein-ligand interactions: The binding of triphenoxypyridine derivatives to factor Xa and trypsin

The binding of a set of 10 triphenoxypyridine derivatives to two serine proteases, factor Xa and trypsin, has been used to analyze factors related to sampling and convergence in free energy calculations based on molecular dynamics simulation techniques. The inhibitors investigated were initially proposed as part of the Critical Assessment of Techniques for Free Energy Evaluation (CATFEE) project for which no experimental results nor any assessment of the predictions submitted by various groups have ever been published. The inhibitors studied represent a severe challenge for explicit free energy calculations. The mutations from one compound to another involve up to 19 atoms, the creation and annihilation of net charge and several alternate binding modes. Nevertheless, we demonstrate that it is possible to obtain highly converged results (+/- 5-10 kJ/mol) even for such complex multi-atom mutations by simulating on a nanosecond time scale. This is achieved by using soft-core potentials to facilitate the creation and deletion of atoms and by a careful choice of mutation pathway. The results show that given modest computational resources, explicit free energy calculations can be successfully applied to realistic problems in drug design.     

Solvated interaction energy (SIE) for scoring protein-ligand binding affinities. 1. Exploring the parameter space

We present a binding free energy function that consists of force field terms supplemented by solvation terms. We used this function to calibrate the solvation model along with the binding interaction terms in a self-consistent manner. The motivation for this approach was that the solute dielectric-constant dependence of calculated hydration gas-to-water transfer free energies is markedly different from that of binding free energies (J. Comput. Chem. 2003, 24, 954). Hence, we sought to calibrate directly the solvation terms in the context of a binding calculation. The five parameters of the model were systematically scanned to best reproduce the absolute binding free energies for a set of 99 protein-ligand complexes. We obtained a mean unsigned error of 1.29 kcal/mol for the predicted absolute binding affinity in a parameter space that was fairly shallow near the optimum. The lowest errors were obtained with solute dielectric values of Din = 20 or higher and scaling of the intermolecular van der Waals interaction energy by factors ranging from 0.03 to 0.15. The high apparent Din and strong van der Waals scaling may reflect the anticorrelation of the change in solvated potential energy and configurational entropy, that is, enthalpy-entropy compensation in ligand binding (Biophys. J. 2004, 87, 3035-3049). Five variations of preparing the protein-ligand data set were explored in order to examine the effect of energy refinement and the presence of bound water on the calculated results. We find that retaining water in the final protein structure used for calculating the binding free energy is not necessary to obtain good results; that is the continuum solvation model is sufficient. Virtual screening enrichment studies on estrogen receptor and thymidine kinase showed a good ability of the binding free energy function to recover true hits in a collection of decoys.     

Solvated interaction energy (SIE) for scoring protein-ligand binding affinities. 2. Benchmark in the CSAR-2010 scoring exercise

Solvated interaction energy (SIE) is an end-point physics-based scoring function for predicting binding affinities from force-field nonbonded interaction terms, continuum solvation, and configurational entropy linear compensation. We tested the SIE function in the Community Structure-Activity Resource (CSAR) scoring challenge consisting of high-resolution cocrystal structures for 343 protein-ligand complexes with high-quality binding affinity data and high diversity with respect to protein targets. Particular emphasis was placed on the sensitivity of SIE predictions to the assignment of protonation and tautomeric states in the complex and the treatment of metal ions near the protein-ligand interface. These were manually curated from an originally distributed CSAR-HiQ data set version, leading to the currently distributed CSAR-NRC-HiQ version. We found that this manual curation was a critical step for accurately testing the performance of the SIE function. The standard SIE parametrization, previously calibrated on an independent data set, predicted absolute binding affinities with a mean-unsigned-error (MUE) of 2.41 kcal/mol for the CSAR-HiQ version, which improved to 1.98 kcal/mol for the upgraded CSAR-NRC-HiQ version. Half-half retraining-testing of SIE parameters on two predefined subsets of CSAR-NRC-HiQ led to only marginal further improvements to an MUE of 1.83 kcal/mol. Hence, we do not recommend altering the current default parameters of SIE at this time. For a sample of SIE outliers, additional calculations by molecular dynamics-based SIE averaging with or without incorporation of ligand strain, by MM-PB(GB)/SA methods with or without entropic estimates, or even by the linear interaction energy (LIE) formalism with an explicit solvent model, did not further improve predictions.     

The origins of femtomolar protein-ligand binding: hydrogen-bond cooperativity and desolvation energetics in the biotin-(strept)avidin binding site

The unusually strong reversible binding of biotin by avidin and streptavidin has been investigated by density functional and MP2 ab initio quantum mechanical methods. The solvation of biotin by water has also been studied through QM/MM/MC calculations. The ureido moiety of biotin in the bound state hydrogen bonds to five residues, three to the carbonyl oxygen and one for each--NH group. These five hydrogen bonds act cooperatively, leading to stabilization that is larger than the sum of individual hydrogen-bonding energies. The charged aspartate is the key residue that provides the driving force for cooperativity in the hydrogen-bonding network for both avidin and streptavidin by greatly polarizing the urea of biotin. If the residue is removed, the network is disrupted, and the attenuation of the energetic contributions from the neighboring residues results in significant reduction of cooperative interactions. Aspartate is directly hydrogen-bonded with biotin in streptavidin and is one residue removed in avidin. The hydrogen-bonding groups in streptavidin are computed to give larger cooperative hydrogen-bonding effects than avidin. However, the net gain in electrostatic binding energy is predicted to favor the avidin-bicyclic urea complex due to the relatively large penalty for desolvation of the streptavidin binding site (specifically expulsion of bound water molecules). QM/MM/MC calculations involving biotin and the ureido moiety in aqueous solution, featuring PDDG/PM3, show that water interactions with the bicyclic urea are much weaker than (strept)avidin interactions due to relatively low polarization of the urea group in water.     

Thermodynamics of buried water clusters at a protein-ligand binding interface

The structure of the complex of cyclophilin A (CypA) with cyclosporin A (CsA, 1) shows a cluster of four water molecules buried at the binding interface, which is rearranged when CsA is replaced by (5-hydroxynorvaline)-2-cyclosporin (2). The thermodynamic contributions of each bound water molecule in the two complexes are explored with the inhomogeneous fluid solvation theory and molecular dynamics simulations. Water (WTR) 133 in complex 1 contributes little to the binding affinity, while WTR6 and 7 in complex 2 play an essential role in mediating protein-ligand binding with a hydrogen bond network. The calculations reveal that the rearrangement of the water molecules contributes favorably to the binding affinity, even though one of them is displaced going from ligand 1 to 2. Another favorable contribution comes from the larger protein-ligand interactions of ligand 2. However, these favorable contributions are not sufficient to overcome the unfavorable desolvation free energy change and the conformational entropy of the hydroxylpropyl group of ligand 2 in the complex, leading to a lower binding affinity of ligand 2. These physical insights may be useful in the development of improved scoring functions for binding affinity prediction.     

Measurement of the protein-ligand bond energy of carboxymyoglobin by pulsed photoacoustic calorimetry

Pulsed photoacoustic spectroscopy has been used to study two horse heart myoglobin derivatives, deoxymyoglobin and carboxymyoglobin, in the temperature range of 0–25° C. The myoglobin-CO bond energy has been measured to be 13.4±0.5 kcal/mole. The volume change associated with the photodissociation of carboxymyoglobin into deoxymyoglobin and CO has also been observed.     

Information theory-based scoring function for the structure-based prediction of protein-ligand binding affinity

The development and validation of a new knowledge based scoring function (SIScoreJE) to predict binding energy between proteins and ligands is presented. SIScoreJE efficiently predicts the binding energy between a small molecule and its protein receptor. Protein-ligand atomic contact information was derived from a Non-Redundant Data set (NRD) of over 3000 X-ray crystal structures of protein-ligand complexes. This information was classified for individual "atom contact pairs" (ACP) which is used to calculate the atomic contact preferences. In addition to the two schemes generated in this study we have assessed a number of other common atom-type classification schemes. The preferences were calculated using an information theoretic relationship of joint entropy. Among 18 different atom-type classification schemes "ScoreJE Atom Type set2" (SATs2) was found to be the most suitable for our approach. To test the sensitivity of the method to the inclusion of solvent, Single-body Solvation Potentials (SSP) were also derived from the atomic contacts between the protein atom types and water molecules modeled using AQUARIUS2. Validation was carried out using an evaluation data set of 100 protein-ligand complexes with known binding energies to test the ability of the scoring functions to reproduce known binding affinities. In summary, it was found that a combined SSP/ScoreJE (SIScoreJE) performed significantly better than ScoreJE alone, and SIScoreJE and ScoreJE performed better than GOLD::GoldScore, GOLD::ChemScore, and XScore.     

A general mass spectrometry-based assay for the quantitation of protein-ligand binding interactions in solution

A new method that utilizes matrix-assisted laser desorption/ionization (MALDI) mass spectrometry and exploits the hydrogen/deuterium (H/D) exchange properties of proteins was developed for measuring the thermodynamic properties of protein-ligand complexes in solution. Dissociation constants (Kd values) determined by the method for five model protein-ligand complexes that included those with small molecules, nucleic acids, peptides, and other proteins were generally in good agreement with Kd values measured by conventional methods. Important experimental advantages of the described method over existing methods include: the ability to make measurements in a high-throughput and automated fashion, the ability to make measurements using only picomole quantitities of protein, and the ability to analyze either purified or unpurified protein-ligand complexes.     

Calculations of protein-ligand binding entropy of relative and overall molecular motions

In the context of virtual database screening, calculations of protein-ligand binding entropy of relative and overall molecular motions are challenging, owing to the inherent structural complexity of the ligand binding well in the energy landscape of protein-ligand interactions and computing time limitations. We describe a fast statistical thermodynamic method for estimation the binding entropy to address the challenges. The method is based on the integration of the configurational integral over clusters obtained from multiple docked positions. We apply the method in conjunction with 11 popular scoring functions (AutoDock, ChemScore, DrugScore, D-Score, F-Score, G-Score, LigScore, LUDI, PLP, PMF, X-Score) to evaluate the binding entropy of 100 protein-ligand complexes. The averaged values of binding entropy contribution vary from 6.2 to 9.1 kcal/mol, showing good agreement with literature. We calculate positional sizes and the angular volume of the native ligand wells. The averaged geometric mean of positional sizes in principal directions varies from 0.8 to 1.4 A. The calculated range of angular volumes is 3.3-11.8 rad(2). Then we demonstrate that the averaged six-dimensional volume of the native well is larger than the volume of the most populated non-native well in energy landscapes described by all of 11 scoring functions.     

Theory and simulation on the kinetics of protein-ligand binding coupled to conformational change

Conformational change during protein-ligand binding may significantly affect both the binding mechanism and the rate constant. Most earlier theories and simulations treated conformational change as stochastic gating with transition rates between reactive and nonreactive conformations uncoupled to ligand binding. Recently, we introduced a dual-transition-rates model in which the transition rates between reactive and nonreactive conformations depend on the protein-ligand distance [H.-X. Zhou, Biophys. J. 98, L15 (2010)]. Analytical results of that model showed that the apparent binding mechanism switches from conformational selection to induced fit, when the rates of conformational transitions increase from being much slower than the diffusional approach of the protein-ligand pair to being much faster. The conformational-selection limit (k(CS)) and the induced-fit limit (k(IF)) provide lower and upper bounds, respectively, for the binding rate constant. Here we introduce a general model in which the energy surface of the protein in conformational space is coupled to ligand binding, and present a method for calculating the binding rate constant from Brownian dynamics simulations. Analytical and simulation results show that, for an energy surface that switches from favoring the nonreactive conformation while the ligand is away to favoring the reactive conformation while the ligand is near, k(CS) and k(IF) become close and, thus, provide tight bounds to the binding rate constant. This finding has significant mechanistic implications and presents routes for obtaining good estimates of the rate constant at low cost.     

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.