Revisiting Free Energy Calculations: A Theoretical Connection to MM/PBSA and Direct Calculation of the Association Free Energy

The prediction of absolute ligand-receptor binding affinities is essential in a wide range of biophysical queries, from the study of protein-protein interactions to structure-based drug design. End-point free energy methods, such as the Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) model, have received much attention and widespread application in recent literature. These methods benefit from computational efficiency as only the initial and final states of the system are evaluated, yet there remains a need for strengthening their theoretical foundation. Here a clear connection between statistical thermodynamics and end-point free energy models is presented. The importance of the association free energy, arising from one molecule's loss of translational and rotational freedom from the standard state concentration, is addressed. A novel method for calculating this quantity directly from a molecular dynamics simulation is described. The challenges of accounting for changes in the protein conformation and its fluctuations from separate simulations are discussed. A simple first-order approximation of the configuration integral is presented to lay the groundwork for future efforts. This model has been applied to FKBP12, a small immunophilin that has been widely studied in the drug industry for its potential immunosuppressive and neuroregenerative effects.     

The “cratic correction” and related fallacies

The “cratic correction” for removing the contribution of translational motion from standard free energies and entropies of chemical reactions is examined. Although removal of translational effects is useful in interpreting free energy and entropy changes in molecular terms, this particular correction does not do so. No basis for it is found in thermodynamics or statistical mechanics. The appropriate correction for center-of-mass translational motion is derived from classical statistical mechanics. However, the importance of recognizing the proper molecular meaning of the remaining free energy is stressed. In any case, this latter correction, however legitimate, is shown to be futile, since it has long been recognized from classical statistical mechanics that all terms involving momenta (thermal wavelength terms) —be they classified as translational, rotational, vibrational, or whatever—necessarily cancel out informing the standard free energy change of any chemical reaction from the constituent standard chemical potentials. © 1995 John Wiley & Sons, Inc.     

Absolute binding free energy calculations using molecular dynamics simulations with restraining potentials

The absolute (standard) binding free energy of eight FK506-related ligands to FKBP12 is calculated using free energy perturbation molecular dynamics (FEP/MD) simulations with explicit solvent. A number of features are implemented to improve the accuracy and enhance the convergence of the calculations. First, the absolute binding free energy is decomposed into sequential steps during which the ligand-surrounding interactions as well as various biasing potentials restraining the translation, orientation, and conformation of the ligand are turned "on" and "off." Second, sampling of the ligand conformation is enforced by a restraining potential based on the root mean-square deviation relative to the bound state conformation. The effect of all the restraining potentials is rigorously unbiased, and it is shown explicitly that the final results are independent of all artificial restraints. Third, the repulsive and dispersive free energy contribution arising from the Lennard-Jones interactions of the ligand with its surrounding (protein and solvent) is calculated using the Weeks-Chandler-Andersen separation. This separation also improves convergence of the FEP/MD calculations. Fourth, to decrease the computational cost, only a small number of atoms in the vicinity of the binding site are simulated explicitly, while all the influence of the remaining atoms is incorporated implicitly using the generalized solvent boundary potential (GSBP) method. With GSBP, the size of the simulated FKBP12/ligand systems is significantly reduced, from approximately 25,000 to 2500. The computations are very efficient and the statistical error is small ( approximately 1 kcal/mol). The calculated binding free energies are generally in good agreement with available experimental data and previous calculations (within approximately 2 kcal/mol). The present results indicate that a strategy based on FEP/MD simulations of a reduced GSBP atomic model sampled with conformational, translational, and orientational restraining potentials can be computationally inexpensive and accurate.     

Loss of translational entropy in molecular associations

Molecular associations in solution are opposed by the loss of entropy (DeltaS) that results from the restriction of motion of each component in the complex. Theoretical estimates of DeltaS are essential for rationalizing binding affinities, as well as for calculating entropic contribution to enzyme catalysis. Recently a statistical-mechanical framework has been proposed for estimating efficiently the translational entropy loss (DeltaS(trsl)), while taking explicitly into account the complex intermolecular interactions between the solute and the solvent. This framework relates the translational entropy of a solute in solution to its "free volume," defined as the volume accessible to the center of mass of the solute in the presence of the solvent and calculated by using an extension of the cell model (CM) for condensed phases. The translational entropy of pure water, estimated with the CM algorithm, shows good agreement with the experimental information. The free volume of various solutes in water, calculated within the CM by using molecular dynamics simulations with explicit solvent, displays a strong correlation with the solutes' polar and total surface areas. This correlation is used to propose a parameterization that can be used to calculate routinely the translational entropy of a solute in water. We also applied the CM formalism to calculate the free volume and translational entropy loss (DeltaS(trsl)) on binding of benzene to a cavity in a mutant T4-lysozyme. Our results agree with previously published estimates of the binding of benzene to this mutant T4-lysozyme. These and other considerations suggest that the cell model is a simple yet efficient theoretical framework to evaluate the translational entropy loss on molecular association in solution.     

Thermodynamics of the isothermal denaturation of lysozyme by guanidinium chloride.

A thermodynamic analysis of the isothermal denaturation of lysozyme by guanidinium chloride has been performed. The analysis is based on the equation which relates the equilibrium constant for denaturation to the preferential binding of denaturant. The equation has been derived previously by thermodynamic methods, whereas in this article a derivation based on statistical mechanics is given. By application of the equation the free energy of denaturation is first calculated and from it, by subtracting the calorimetrically-determined enthalpy of denaturation, the entropy of denaturation is determined.     

Estimation of the hydrophobic effect in an antigen-antibody protein-protein interface

Antigen-antibody complexes provide useful models for analyzing the thermodynamics of protein-protein association reactions. We have employed site-directed mutagenesis, X-ray crystallography, and isothermal titration calorimetry to investigate the role of hydrophobic interactions in stabilizing the complex between the Fv fragment of the anti-hen egg white lysozyme (HEL) antibody D1.3 and HEL. Crystal structures of six FvD1.3-HEL mutant complexes in which an interface tryptophan residue (V(L)W92) has been replaced by residues with smaller side chains (alanine, serine, valine, aspartate, histidine, and phenylalanine) were determined to resolutions between 1.75 and 2.00 A. In the wild-type complex, V(L)W92 occupies a large hydrophobic pocket on the surface of HEL and constitutes an energetic "hot spot" for antigen binding. The losses in apolar buried surface area in the mutant complexes, relative to wild-type, range from 25 (V(L)F92) to 115 A(2) (V(L)A92), with no significant shifts in the positions of protein atoms at the mutation site for any of the complexes except V(L)A92, where there is a peptide flip. The affinities of the mutant Fv fragments for HEL are 10-100-fold lower than that of the original antibody. Formation of all six mutant complexes is marked by a decrease in binding enthalpy that exceeds the decrease in binding free energy, such that the loss in enthalpy is partly offset by a compensating gain in entropy. No correlation was observed between decreases in apolar, polar, or aggregate (sum of the apolar and polar) buried surface area in the V(L)92 mutant series and changes in the enthalpy of formation. Conversely, there exist linear correlations between losses of apolar buried surface and decreases in binding free energy (R(2) = 0.937) as well as increases in the solvent portion of the entropy of binding (R(2) = 0.909). The correlation between binding free energy and apolar buried surface area corresponds to 21 cal mol(-1) A(-2) (1 cal = 4.185 J) for the effective hydrophobicity at the V(L)92 mutation site. Furthermore, the slope of the line defined by the correlation between changes in binding free energy and solvent entropy approaches unity, demonstrating that the exclusion of solvent from the binding interface is the predominant energetic factor in the formation of this protein complex. Our estimate of the hydrophobic contribution to binding at site V(L)92 in the D1.3-HEL interface is consistent with values for the hydrophobic effect derived from classical hydrocarbon solubility models. We also show how residue V(L)W92 can contribute significantly less to stabilization when buried in a more polar pocket, illustrating the dependence of the hydrophobic effect on local environment at different sites in a protein-protein interface.     

Accounting for Ligand Conformational Restriction in Calculations of Protein-Ligand Binding Affinities

The conformation adopted by a ligand on binding to a receptor may differ from its lowest-energy conformation in solution. In addition, the bound ligand is more conformationally restricted, which is associated with a configurational entropy loss. The free energy change due to these effects is often neglected or treated crudely in current models for predicting binding affinity. We present a method for estimating this contribution, based on perturbation theory using the quasi-harmonic model of Karplus and Kushick as a reference system. The consistency of the method is checked for small model systems. Subsequently we use the method, along with an estimate for the enthalpic contribution due to ligand-receptor interactions, to calculate relative binding affinities. The AMBER force field and generalized Born implicit solvent model is used. Binding affinities were estimated for a test set of 233 protein-ligand complexes for which crystal structures and measured binding affinities are available. In most cases, the ligand conformation in the bound state was significantly different from the most favorable conformation in solution. In general, the correlation between measured and calculated ligand binding affinities including the free energy change due to ligand conformational change is comparable to or slightly better than that obtained by using an empirically-trained docking score. Both entropic and enthalpic contributions to this free energy change are significant.     

Conformational dynamics and thermodynamics of protein-ligand binding studied by NMR relaxation

Protein conformational dynamics can be critical for ligand binding in two ways that relate to kinetics and thermodynamics respectively. First, conformational transitions between different substates can control access to the binding site (kinetics). Secondly, differences between free and ligand-bound states in their conformational fluctuations contribute to the entropy of ligand binding (thermodynamics). In the present paper, I focus on the second topic, summarizing our recent results on the role of conformational entropy in ligand binding to Gal3C (the carbohydrate-recognition domain of galectin-3). NMR relaxation experiments provide a unique probe of conformational entropy by characterizing bond-vector fluctuations at atomic resolution. By monitoring differences between the free and ligand-bound states in their backbone and side chain order parameters, we have estimated the contributions from conformational entropy to the free energy of binding. Overall, the conformational entropy of Gal3C increases upon ligand binding, thereby contributing favourably to the binding affinity. Comparisons with the results from isothermal titration calorimetry indicate that the conformational entropy is comparable in magnitude to the enthalpy of binding. Furthermore, there are significant differences in the dynamic response to binding of different ligands, despite the fact that the protein structure is virtually identical in the different protein-ligand complexes. Thus both affinity and specificity of ligand binding to Gal3C appear to depend in part on subtle differences in the conformational fluctuations that reflect the complex interplay between structure, dynamics and ligand interactions.     

Targeting the cell wall of Mycobacterium tuberculosis: a molecular modeling investigation of the interaction of imipenem and meropenem with L,D-transpeptidase 2

The single crystal X-ray structure of the extracellular portion of the L,D-transpeptidase (ex-Ldt_Mt2 – residues 120–408) enzyme was recently reported. It was observed that imipenem and meropenem inhibit activity of this enzyme, responsible for generating L,D-transpeptide linkages in the peptidoglycan layer of Mycobacterium tuberculosis. Imipenem is more active and isothermal titration calorimetry experiments revealed that meropenem is subjected to an entropy penalty upon binding to the enzyme. Herein, we report a molecular modeling approach to obtain a molecular view of the inhibitor/enzyme interactions. The average binding free energies for nine commercially available inhibitors were calculated using MM/GBSA and Solvation Interaction Energy (SIE) approaches and the calculated energies corresponded well with the available experimentally observed results. The method reproduces the same order of binding energies as experimentally observed for imipenem and meropenem. We have also demonstrated that SIE is a reasonably accurate and cost-effective free energy method, which can be used to predict carbapenem affinities for this enzyme. A theoretical explanation was offered for the experimental entropy penalty observed for meropenem, creating optimism that this computational model can serve as a potential computational model for other researchers in the field.     

Nucleic acid base-pairing and N-methylacetamide self-association in chloroform: affinity and conformation.

A recently developed computational method, 'mining minima', is used to examine the hydrogen-bonding interactions of nucleic acid base-pairs and of the N-methylacetamide homodimer in chloroform. The mining minima algorithm aggressively samples molecular conformations, identifies the most important local minima, and computes their contributions to the overall free energy of the system. Here, the CHARMM 98 parameter set is used for the potential energy and the generalized Born/surface area solvent model is used to account for the influence of the solvent. Good agreement with experiment is obtained for the non-covalent binding affinities of a series of complexes. The computational approach used here is applicable to a range of molecular systems.     

Direct measurement of protein binding energetics by isothermal titration calorimetry.

Of all the techniques that are currently available to measure binding, isothermal titration calorimetry is the only one capable of measuring not only the magnitude of the binding affinity but also the magnitude of the two thermodynamic terms that define the binding affinity: the enthalpy (AH) and entropy (AS) changes. Recent advances in instrumentation have facilitated the development of experimental designs that permit the direct measurement of arbitrarily high binding affinities, the coupling of binding to protonation/deprotonation processes and the analysis of binding thermodynamics in terms of structural parameters. Because isothermal titration calorimetry has the capability to measure different energetic contributions to the binding affinity, it provides a unique bridge between computational and experimental analysis. As such, it is increasingly becoming an essential tool in molecular design.     

Calculation of absolute protein-ligand binding affinity using path and endpoint approaches

A comparative analysis is provided of rigorous and approximate methods for calculating absolute binding affinities of two protein-ligand complexes: the FKBP protein bound with small molecules 4-hydroxy-2-butanone and FK506. Our rigorous approach is an umbrella sampling technique where a potential of mean force is determined by pulling the ligand out of the protein active site over several simulation windows. The results of this approach agree well with experimentally observed binding affinities. Also assessed is a commonly used approximate endpoint approach, which separately estimates enthalpy, solvation free energy, and entropy. We show that this endpoint approach has numerous variations, all of which are prone to critical shortcomings. For example, conventional harmonic and quasiharmonic entropy estimation procedures produce disparate results for the relatively simple protein-ligand systems studied in this work.     

Predicting the conformational states of cyclic tetrapeptides

Biologically active cyclic tetrapeptides, usually found among fungi metabolites, exhibit phytotoxic or cytostatic activities that are likely to be governed by specific conformations adopted in solution. For conformational studies and drug design, there is a strong interest in using fast and reliable methods to determine correctly the conformational population of cyclotetrapeptides. We show here that standard molecular mechanics computational approach gives satisfactory results. The method was validated step by step by experimental data either obtained after synthesis and NMR analysis, or found in the literature. The cyclo(Gly)(4), cyclo(Ala)(4), cyclo(Sar)(4), and cyclo(SarGly)(2) peptides were used to evaluate the prediction of the peptide backbone conformation, and the detailed conformational analysis of tentoxin, a natural phytotoxic cyclotetrapeptide in which N-alkylated peptide bonds alternate with regular secondary ones, was used to validate the computation of conformers proportions. From the knowledge of an initial cyclic primary structure and of the D or L configuration of the amino acids, we show that it is possible to determine the exact orientation of carbonyl groups and to predict the nature of conformers present in solution. The proportion of each conformer can be inferred from a statistical thermodynamics approach by using the potential energy values of each conformer, computed by molecular mechanics methods with the TRIPOS force field, which allowed us to account for the solvent. The solvent contribution was processed by two different methods according to the nature of the interactions: whether through the dielectric constant introduced in the electrostatic potential, when interaction with solute molecules are weak or negligible, or through the computation of free energy of solvation using the algorithm SILVERWARE for solvents explicitly interacting with the solute. When applied to tentoxin, this conformational analysis yielded results in very good agreement with the experimental data reported by Pinet et al. (Biopolymers, 1995, Vol. 36, pp. 135-152), on both the nature of existing conformers and their relative proportions, whatever the nature of the considered solvent.     

Calculation of cyclodextrin binding affinities: energy, entropy, and implications for drug design

The second generation Mining Minima method yields binding affinities accurate to within 0.8 kcal/mol for the associations of alpha-, beta-, and gamma-cyclodextrin with benzene, resorcinol, flurbiprofen, naproxen, and nabumetone. These calculations require hours to a day on a commodity computer. The calculations also indicate that the changes in configurational entropy upon binding oppose association by as much as 24 kcal/mol and result primarily from a narrowing of energy wells in the bound versus the free state, rather than from a drop in the number of distinct low-energy conformations on binding. Also, the configurational entropy is found to vary substantially among the bound conformations of a given cyclodextrin-guest complex. This result suggests that the configurational entropy must be accounted for to reliably rank docked conformations in both host-guest and ligand-protein complexes. In close analogy with the common experimental observation of entropy-enthalpy compensation, the computed entropy changes show a near-linear relationship with the changes in mean potential plus solvation energy.     

Fatty acid interactions with native and mutant fatty acid binding proteins

The interactions of long chain fatty acids (FA) with wild type (WT) fatty acid binding proteins (FABP) and engineered FABP mutants have been monitored to determine the equilibrium binding constants as well as the rate constants for binding and dissociation. These measurements have been done using the fluorescent probes, ADIFAB and ADIFAB2, that allow the determination of the free fatty acid (FFA) concentration in the reaction of FA with proteins and membranes. The results of these studies indicate that for WT proteins from adipocyte, heart, intestine, and liver, Kd values are in the nM range and affinities decrease with increasing aqueous solubility of the FA. Binding affinities for heart and liver are generally greater than those for adipocyte and intestine. Moreover, measurements of the rate constants indicate that binding equilibrium at 37øC is achieved within seconds for all FA and FABPs. These results, together with the level of serum (unbound) FFA, suggests a buffering action of FABPs that helps to maintain the intracellular concentration of FFA so that the flux of FFA between serum and cells occurs down a concentration gradient. Measurements of the temperature dependence of binding reveal that the free energy is predominately enthalpic and that the enthalpy of the reaction results from FA-FABP interactions within the binding cavity. The nature of these interactions were investigated by determining the thermodynamics of binding to engineered point mutants of the intestinal FABP. These measurements showed that binding affinities did not report accurately the changes in protein-FA interactions because changes in the binding entropy and enthalpy tend to compensate. For example, an alanine substitution for arginine 106 yields a 30 fold increase in binding affinity, because the loss in enthalpy due to the elimination of the favorable interaction between the FA carboxylate and Arg106, is more than compensated for by an increase in entropy. Thus understanding the effects of amino acid replacements on FA-FABP interactions requires measurements of enthalpy and entropy, in addition to affinity.     

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.     

Simulating Metabolism with Statistical Thermodynamics

New methods are needed for large scale modeling of metabolism that predict metabolite levels and characterize the thermodynamics of individual reactions and pathways. Current approaches use either kinetic simulations, which are difficult to extend to large networks of reactions because of the need for rate constants, or flux-based methods, which have a large number of feasible solutions because they are unconstrained by the law of mass action. This report presents an alternative modeling approach based on statistical thermodynamics. The principles of this approach are demonstrated using a simple set of coupled reactions, and then the system is characterized with respect to the changes in energy, entropy, free energy, and entropy production. Finally, the physical and biochemical insights that this approach can provide for metabolism are demonstrated by application to the tricarboxylic acid (TCA) cycle of Escherichia coli. The reaction and pathway thermodynamics are evaluated and predictions are made regarding changes in concentration of TCA cycle intermediates due to 10- and 100-fold changes in the ratio of NAD⁺:NADH concentrations. Finally, the assumptions and caveats regarding the use of statistical thermodynamics to model non-equilibrium reactions are discussed.     

Thermodynamics of T-cell receptor-peptide/MHC interactions: progress and opportunities

alphabeta T-cell receptors (TCRs) recognize peptide antigens presented by class I or class II major histocompatibility complex molecules (pMHC). Here we review the use of thermodynamic measurements in the study of TCR-pMHC interactions, with attention to the diversity in binding thermodynamics and how this is related to the variation in TCR-pMHC interfaces. We show that there is no enthalpic or entropic signature for TCR binding; rather, enthalpy and entropy changes vary in a compensatory manner that reflects a narrow free energy window for the interactions that have been characterized. Binding enthalpy and entropy changes do not correlate with structural features such as buried surface area or the number of hydrogen bonds within TCR-pMHC interfaces, possibly reflecting the myriad of contributors to binding thermodynamics, but likely also reflecting a reliance on van't Hoff over calorimetric measurements and the unaccounted influence of equilibria linked to binding. TCR-pMHC binding heat capacity changes likewise vary considerably. In some cases, the heat capacity changes are consistent with conformational differences between bound and free receptors, but there is little data indicating these conformational differences represent the need to organize disordered CDR loops. In this regard, we discuss how thermodynamics may provide additional insight into conformational changes occurring upon TCR binding. Finally, we highlight opportunities for the further use of thermodynamic measurements in the study of TCR-pMHC interactions, not only for understanding TCR binding in general, but also for understanding specifics of individual interactions and the engineering of TCRs with desired molecular recognition properties.