PEARLS: program for energetic analysis of receptor-ligand system

Analysis of the energetics of small molecule ligand-protein, ligand-nucleic acid, and protein-nucleic acid interactions facilitates the quantitative understanding of molecular interactions that regulate the function and conformation of proteins. It has also been extensively used for ranking potential new ligands in virtual drug screening. We developed a Web-based software, PEARLS (Program for Energetic Analysis of Ligand-Receptor Systems), for computing interaction energies of ligand-protein, ligand-nucleic acid, protein-nucleic acid, and ligand-protein-nucleic acid complexes from their 3D structures. AMBER molecular force field, Morse potential, and empirical energy functions are used to compute the van der Waals, electrostatic, hydrogen bond, metal-ligand bonding, and water-mediated hydrogen bond energies between the binding molecules. The change in the solvation free energy of molecular binding is estimated by using an empirical solvation free energy model. Contribution from ligand conformational entropy change is also estimated by a simple model. The computed free energy for a number of PDB ligand-receptor complexes were studied and compared to experimental binding affinity. A substantial degree of correlation between the computed free energy and experimental binding affinity was found, which suggests that PEARLS may be useful in facilitating energetic analysis of ligand-protein, ligand-nucleic acid, and protein-nucleic acid interactions. PEARLS can be accessed at http://ang.cz3.nus.edu.sg/cgi-bin/prog/rune.pl.     

Protein-ligand binding free energy estimation using molecular mechanics and continuum electrostatics. Application to HIV-1 protease inhibitors

A method is proposed for the estimation of absolute binding free energy of interaction between proteins and ligands. Conformational sampling of the protein-ligand complex is performed by molecular dynamics (MD) in vacuo and the solvent effect is calculated a posteriori by solving the Poisson or the Poisson-Boltzmann equation for selected frames of the trajectory. The binding free energy is written as a linear combination of the buried surface upon complexation, SASbur, the electrostatic interaction energy between the ligand and the protein, Eelec, and the difference of the solvation free energies of the complex and the isolated ligand and protein, deltaGsolv. The method uses the buried surface upon complexation to account for the non-polar contribution to the binding free energy because it is less sensitive to the details of the structure than the van der Waals interaction energy. The parameters of the method are developed for a training set of 16 HIV-1 protease-inhibitor complexes of known 3D structure. A correlation coefficient of 0.91 was obtained with an unsigned mean error of 0.8 kcal/mol. When applied to a set of 25 HIV-1 protease-inhibitor complexes of unknown 3D structures, the method provides a satisfactory correlation between the calculated binding free energy and the experimental pIC5o without reparametrization.     

Converging free energy estimates: MM-PB(GB)SA studies on the protein-protein complex Ras-Raf

Estimating protein-protein interaction energies is a very challenging task for current simulation protocols. Here, absolute binding free energies are reported for the complex H-Ras/C-Raf1 using the MM-PB(GB)SA approach, testing the internal consistency and model dependence of the results. Averaging gas-phase energies (MM), solvation free energies as determined by Generalized Born models (GB/SA), and entropic contributions calculated by normal mode analysis for snapshots obtained from 10 ns explicit-solvent molecular dynamics in general results in an overestimation of the binding affinity when a solvent-accessible surface area-dependent model is used to estimate the nonpolar solvation contribution. Applying the sum of a cavity solvation free energy and explicitly modeled solute-solvent van der Waals interaction energies instead provides less negative estimates for the nonpolar solvation contribution. When the polar contribution to the solvation free energy is determined by solving the Poisson-Boltzmann equation (PB) instead, the calculated binding affinity strongly depends on the atomic radii set chosen. For three GB models investigated, different absolute deviations from PB energies were found for the unbound proteins and the complex. As an alternative to normal-mode calculations, quasiharmonic analyses have been performed to estimate entropic contributions due to changes of solute flexibility upon binding. However, such entropy estimates do not converge after 10 ns of simulation time, indicating that sampling issues may limit the applicability of this approach. Finally, binding free energies estimated from snapshots of the unbound proteins extracted from the complex trajectory result in an underestimate of binding affinity. This points to the need to exercise caution in applying the computationally cheaper "one-trajectory-alternative" to systems where there may be significant changes in flexibility and structure due to binding. The best estimate for the binding free energy of Ras-Raf obtained in this study of -8.3 kcal mol(-1) is in good agreement with the experimental result of -9.6 kcal mol(-1), however, further probing the transferability of the applied protocol that led to this result is necessary.     

On the origins of core-electron chemical shifts of small biomolecules in aqueous solution: insights from photoemission and ab initio calculations of glycine(aq)

The local electronic structure of glycine in neutral, basic, and acidic aqueous solution is studied experimentally by X-ray photoelectron spectroscopy and theoretically by molecular dynamics simulations accompanied by first-principle electronic structure and spectrum calculations. Measured and computed nitrogen and carbon 1s binding energies are assigned to different local atomic environments, which are shown to be sensitive to the protonation/deprotonation of the amino and carboxyl functional groups at different pH values. We report the first accurate computation of core-level chemical shifts of an aqueous solute in various protonation states and explicitly show how the distributions of photoelectron binding energies (core-level peak widths) are related to the details of the hydrogen bond configurations, i.e. the geometries of the water solvation shell and the associated electronic screening. The comparison between the experiments and calculations further enables the separation of protonation-induced (covalent) and solvent-induced (electrostatic) screening contributions to the chemical shifts in the aqueous phase. The present core-level line shape analysis facilitates an accurate interpretation of photoelectron spectra from larger biomolecular solutes than glycine.     

Calculations of the free energy of interaction of the c-Fos-c-Jun coiled coil: effects of the solvation model and the inclusion of polarization effects

The leucine zipper region of activator protein-1 (AP-1) comprises the c-Jun and c-Fos proteins and constitutes a well-known coiled coil protein-protein interaction motif. We have used molecular dynamics (MD) simulations in conjunction with the molecular mechanics/Poisson-Boltzmann generalized-Born surface area [MM/PB(GB)SA] methods to predict the free energy of interaction of these proteins. In particular, the influence of the choice of solvation model, protein force field, and water potential on the stability and dynamic properties of the c-Fos-c-Jun complex were investigated. Use of the AMBER polarizable force field ff02 in combination with the polarizable POL3 water potential was found to result in increased stability of the c-Fos-c-Jun complex. MM/PB(GB)SA calculations revealed that MD simulations using the POL3 water potential give the lowest predicted free energies of interaction compared to other nonpolarizable water potentials. In addition, the calculated absolute free energy of binding was predicted to be closest to the experimental value using the MM/GBSA method with independent MD simulation trajectories using the POL3 water potential and the polarizable ff02 force field, while all other binding affinities were overestimated.     

The Confine-and-Release Method: Obtaining Correct Binding Free Energies in the Presence of Protein Conformational Change

Free energy calculations are increasingly being used to estimate absolute and relative binding free energies of ligands to proteins. However, computed free energies often appear to depend on the initial protein conformation, indicating incomplete sampling. This is especially true when proteins can change conformation on ligand binding, as free energies associated with these conformational changes are either ignored or assumed to be included by virtue of the sampling performed in the calculation. Here, we show that, in a model protein system (a designed binding site in T4 Lysozyme), conformational changes can make a difference of several kcal/mol in computed binding free energies, and that they are neglected in computed binding free energies if the system remains kinetically trapped in a particular metastable state on simulation timescales. We introduce a general "confine-and-release" framework for free energy calculations that accounts for these free energies of conformational change. We illustrate its use in this model system by demonstrating that an umbrella sampling protocol can obtain converged binding free energies that are independent of the starting protein structure and include these conformational change free energies.     

Thermodynamic integration to predict host-guest binding affinities

An alchemical free energy method with explicit solvent molecular dynamics simulations was applied as part of the blind prediction contest SAMPL3 to calculate binding free energies for seven guests to an acyclic cucurbit-[n]uril host. The predictions included determination of protonation states for both host and guests, docking pose generation, and binding free energy calculations using thermodynamic integration. We found a root mean square error (RMSE) of 3.6 kcal mol(-1) from the reference experimental results, with an R(2) correlation of 0.51. The agreement with experiment for the largest contributor to this error, guest 6, is improved by 1.7 kcal mol(-1) when a periodicity-induced free energy correction is applied. The corrections for the other ligands were significantly smaller, and altogether the RMSE was reduced by 0.4 kcal mol(-1). We link properties of the host-guest systems during simulation to errors in the computed free energies. Overall, we show that charged host-guest systems studied here, initialized in unconfirmed docking poses, present a challenge to accurate alchemical simulation methods.     

Recognition of RNA by amide modified backbone nucleic acids: molecular dynamics simulations of DNA-RNA hybrids in aqueous solution

Thermodynamic and structural properties of a chemically modified DNA-RNA hybrid in which a phosphodiester linkage is replaced by a neutral amide-3 linkage (3'-CH(2)-CONH-5') were investigated using UV melting experiments, molecular dynamics simulations in explicit water, and continuum solvent models. van't Hoff analysis of the experimental UV melting curves suggests that the significant increase of the thermodynamic stability of a 15-mer DNA-RNA with seven alternated amide-3 modifications (+11 degrees C) is mainly due to an increased binding enthalpy. To further evaluate the origin in the observed affinities differences, the electrostatic contribution to the binding free energy was calculated by solving the Poisson-Boltzmann equation numerically. The nonelectrostatic contribution was estimated as the product of a hydrophobic surface tension coefficient and the surface area that is buried upon double strand formation. Structures were taken from 10 ns molecular dynamics simulations computed in a consistent fashion using explicit solvent, counterions, and the particle-mesh Ewald procedure. The present preliminary thermodynamic study suggests that the favorable binding free energy of the amide-3 DNA single strand to the complementary RNA is equally driven by electrostatic and nonpolar contributions to the binding compared to their natural analogues. In addition, molecular dynamics simulations in explicit water were performed on an amide-3 DNA single strand and the corresponding natural DNA. Results from the conformations cluster analysis of the simulated amide-3 DNA single strand ensembles suggest that the 25% of the population sampled within 10 ns has a pre-organized conformation where the sugar C3' endo pucker is favored at the 3'-flanking nucleotides. These structural and thermodynamic features contribute to the understanding of the observed increased affinities of the amide-3 DNA-RNA hybrids at the microscopic level.     

Calculations of the electrostatic free energy contributions to the binding free energy of sulfonamides to carbonic anhydrase

The interactions between biologically important enzymes and drugs are of great interest. In order to address some aspects of these interactions we have initiated a program to investigate enzymedrug interactions. Specifically, the interactions between one of the isozymes of carbonic anhydrase and a family of drugs known as sulfonamides have been studied using computational methods. In particular the electrostatic free energy of binding of carbonic anhydrase II with acetazolamide, methazolamide,p-chlorobenzenesulfonamide,p-aminobenzenesulfonamide and three new compounds (MK1, MK2, and MK3) has been computed using finite-difference Poisson-Boltzmann (FDPB) [1] method and the semimacroscopic version [2, 3] of the protein dipole Langevin dipole (PDLD) method [4]. Both methods, FDPB and PDLD, give similar results for the electrostatic free energy of binding even though different charges and different treatments were used for the protein. The calculated electrostatic binding free energies are in reasonable agreement with the experimental data. The potential and the limitation of electrostatic models for studies of binding energies are discussed.     

Electrostatic solvation free energy of amino acid side chain analogs: implications for the validity of electrostatic linear response in water

Electrostatic free energies of solvation for 15 neutral amino acid side chain analogs are computed. We compare three methods of varying computational complexity and accuracy for three force fields: free energy simulations, Poisson-Boltzmann (PB), and linear response approximation (LRA) using AMBER, CHARMM, and OPLS-AA force fields. We find that deviations from simulation start at low charges for solutes. The approximate PB and LRA produce an overestimation of electrostatic solvation free energies for most of molecules studied here. These deviations are remarkably systematic. The variations among force fields are almost as large as the variations found among methods. Our study confirms that success of the approximate methods for electrostatic solvation free energies comes from their ability to evaluate free energy differences accurately.     

Effects of system net charge and electrostatic truncation on all‐atom constant pH molecular dynamics

Constant pH molecular dynamics offers a means to rigorously study the effects of solution pH on dynamical processes. Here, we address two critical questions arising from the most recent developments of the all‐atom continuous constant pH molecular dynamics (CpHMD) method: (1) What is the effect of spatial electrostatic truncation on the sampling of protonation states? (2) Is the enforcement of electrical neutrality necessary for constant pH simulations? We first examined how the generalized reaction field and force‐shifting schemes modify the electrostatic forces on the titration coordinates. Free energy simulations of model compounds were then carried out to delineate the errors in the deprotonation free energy and salt‐bridge stability due to electrostatic truncation and system net charge. Finally, CpHMD titration of a mini‐protein HP36 was used to understand the manifestation of the two types of errors in the calculated pK_a values. The major finding is that enforcing charge neutrality under all pH conditions and at all time via cotitrating ions significantly improves the accuracy of protonation‐state sampling. We suggest that such finding is also relevant for simulations with particle mesh Ewald, considering the known artifacts due to charge‐compensating background plasma. © 2014 Wiley Periodicals, Inc.     

FAMBE-pH: a fast and accurate method to compute the total solvation free energies of proteins

A fast and accurate method to compute the total solvation free energies of proteins as a function of pH is presented. The method makes use of a combination of approaches, some of which have already appeared in the literature; (i) the Poisson equation is solved with an optimized fast adaptive multigrid boundary element (FAMBE) method; (ii) the electrostatic free energies of the ionizable sites are calculated for their neutral and charged states by using a detailed model of atomic charges; (iii) a set of optimal atomic radii is used to define a precise dielectric surface interface; (iv) a multilevel adaptive tessellation of this dielectric surface interface is achieved by using multisized boundary elements; and (v) 1:1 salt effects are included. The equilibrium proton binding/release is calculated with the Tanford-Schellman integral if the proteins contain more than approximately 20-25 ionizable groups; for a smaller number of ionizable groups, the ionization partition function is calculated directly. The FAMBE method is tested as a function of pH (FAMBE-pH) with three proteins, namely, bovine pancreatic trypsin inhibitor (BPTI), hen egg white lysozyme (HEWL), and bovine pancreatic ribonuclease A (RNaseA). The results are (a) the FAMBE-pH method reproduces the observed pK a's of the ionizable groups of these proteins within an average absolute value of 0.4 p K units and a maximum error of 1.2 p K units and (b) comparison of the calculated total pH-dependent solvation free energy for BPTI, between the exact calculation of the ionization partition function and the Tanford-Schellman integral method, shows agreement within 1.2 kcal/mol. These results indicate that calculation of total solvation free energies with the FAMBE-pH method can provide an accurate prediction of protein conformational stability at a given fixed pH and, if coupled with molecular mechanics or molecular dynamics methods, can also be used for more realistic studies of protein folding, unfolding, and dynamics, as a function of pH.     

Accurate prediction of protonation state as a prerequisite for reliable MM-PB(GB)SA binding free energy calculations of HIV-1 protease inhibitors

Binding free energies were calculated for the inhibitors lopinavir, ritonavir, saquinavir, indinavir, amprenavir, and nelfinavir bound to HIV-1 protease. An MMPB/SA-type analysis was applied to conformational samples from 3 ns explicit solvent molecular dynamics simulations of the enzyme-inhibitor complexes. Binding affinities and the sampled conformations of the inhibitor and enzyme were compared between different HIV-1 protease protonation states to find the most likely protonation state of the enzyme in the complex with each of the inhibitors. The resulting set of protonation states leads to good agreement between calculated and experimental binding affinities. Results from the MMPB/SA analysis are compared with an explicit/implicit hybrid scheme and with MMGB/SA methods. It is found that the inclusion of explicit water molecules may offer a slight advantage in reproducing absolute binding free energies while the use of the Generalized Born approximation significantly affects the accuracy of the calculated binding affinities.     

Constant pH molecular dynamics in generalized Born implicit solvent

A new method is proposed for constant pH molecular dynamics (MD), employing generalized Born (GB) electrostatics. Protonation states are modeled with different charge sets, and titrating residues sample a Boltzmann distribution of protonation states as the simulation progresses, using Monte Carlo sampling based on GB-derived energies. The method is applied to four different crystal structures of hen egg-white lysozyme (HEWL). pK(a) predictions derived from the simulations have root-mean-square (RMS) error of 0.82 relative to experimental values. Similarity of results between the four crystal structures shows the method to be independent of starting crystal structure; this is in contrast to most electrostatics-only models. A strong correlation between conformation and protonation state is noted and quantitatively analyzed, emphasizing the importance of sampling protonation states in conjunction with dynamics.     

Energetic basis for drug resistance of HIV-1 protease mutants against amprenavir

Amprenavir (APV) is a high affinity (0.15 nM) HIV-1 protease (PR) inhibitor. However, the affinities of the drug resistant protease variants V32I, I50V, I54V, I54M, I84V and L90M to amprenavir are decreased 3 to 30-fold compared to the wild-type. In this work, the popular molecular mechanics Poisson-Boltzmann surface area method has been used to investigate the effectiveness of amprenavir against the wild-type and these mutated protease variants. Our results reveal that the protonation state of Asp25/Asp25′ strongly affects the dynamics, the overall affinity and the interactions of the inhibitor with individual residues. We emphasize that, in contrast to what is often assumed, the protonation state may not be inferred from the affinities but requires pK_a calculations. At neutral pH, Asp25 and Asp25′ are ionized or protonated, respectively, as suggested from pK_a calculations. This protonation state was thus mainly considered in our study. Mutation induced changes in binding affinities are in agreement with the experimental findings. The decomposition of the binding free energy reveals the mechanisms underlying binding and drug resistance. Drug resistance arises from an increase in the energetic contribution from the van der Waals interactions between APV and PR (V32I, I50V, and I84V mutant) or a rise in the energetic contribution from the electrostatic interactions between the inhibitor and its target (I54M and I54V mutant). For the V32I mutant, also an increased free energy for the polar solvation contributes to the drug resistance. For the L90M mutant, a rise in the van der Waals energy for APV-PR interactions is compensated by a decrease in the polar solvation free energy such that the net binding affinity remains unchanged. Detailed understanding of the molecular forces governing binding and drug resistance might assist in the design of new inhibitors against HIV-1 PR variants that are resistant against current drugs.     

Calculations of pH-dependent binding of proteins to biological membranes

Binding of proteins to membranes is often accompanied by titration of ionizable residues and is, therefore, dependent on pH. We present a theoretical treatment and computational approach for predicting absolute, pH-dependent membrane binding free energies. The standard free energy of binding, DeltaG, is defined as -RTln(P(b)/P(f)), where P(b) and P(f) are the amounts of bound and free protein. The apparent pK(a) of binding is the pH value at which P(b) and P(f) are equal. Proteins bind to the membrane in the pH range where DeltaG is negative. The components of the binding free energy are (a) the free energy cost of ionization state changes (DeltaG(ion)), (b) the effective energy of transfer from solvent to the membrane surface, (c) the translational/rotational entropy cost of binding, and (d) an ideal entropy term that depends on the relative volume of the bound and free state and therefore depends on lipid concentration. Calculation of the first term requires determination of pK(a) values in solvent and on the membrane surface. All energies required by the method are obtained from molecular dynamics trajectories on an implicit membrane (IMM1-GC). The method is tested on pentalysine and the helical peptide VEEKS, derived from the membrane-binding domain of phosphocholine cytidylyltransferase. The agreement between the measured and the calculated free energies of binding of pentalysine is good. The extent of membrane binding of VEEKS is, however, underestimated compared to experiment. Calculations of the interaction energy between two VEEKS helices on the membrane suggest that the discrepancy is mainly due to the neglect of protein-protein interactions on the membrane surface.     

Alpha7 nicotinic acetylcholine receptor agonists: prediction of their binding affinity through a molecular mechanics Poisson-Boltzmann surface area approach

A group of agonists for the alpha7 neuronal nicotinic acetylcholine receptors (nAChRs) was investigated, and their free energies of binding DeltaG(bind) were calculated by applying the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) approach. This method, based on molecular dynamics simulations of fully solvated protein-ligand complexes, allowed us to estimate the contribution of both polar and nonpolar terms as well as the entropy to the overall free energy of binding. The calculated results were in a good agreement with the experimentally determined DeltaG(bind) values, thereby pointing to the MM-PBSA protocol as a valuable computational tool for the rational design of specific agents targeting the neuronal alpha7 nAChR subtypes.     

Calculation of absolute binding free energies for charged ligands and effects of long-range electrostatic interactions

A recently proposed molecular dynamics method for estimating binding free energies is applied to the complexation of two charged benzamidine inhibitors with trypsin. The difficulties with calculations of absolute binding energies for charged molecules, associated with long-range electrostatic contributions, are discussed and it is shown how to deal with these effectively. In particular, energetic effects caused by the trunction of dipole-dipole interactions in the medium surrounding the charged ligand are examined and found to be significant. Calculations of the absolute binding energy for benzamidine using the free energy perturbation approach are also reported. These simulations illustrate the typical problems associated with annihilation transformations of molecules bound inside proteins. © 1996 by John Wiley & Sons, Inc.     

Calculation of ligand binding free energies from molecular dynamics simulations

A recently developed method for predicting binding affinities in ligand–receptor complexes, based on interaction energy averaging and conformational sampling by molecular dynamics simulation, is presented. Polar and nonpolar contributions to the binding free energy are approximated by a linear scaling of the corresponding terms in the average intermolecular interaction energy for the bound and free states of the ligand. While the method originally assumed the validity of electrostatic linear response, we show that incorporation of systematic deviations from linear response derived from free energy perturbation calculations enhances the accuracy of the approach. The method is applied to complexes of wild-type and mutant human dihydrofolate reductases with 2,4-diaminopteridine and 2,4-diaminoquinazoline inhibitors. It is shown that a binding energy accuracy of about 1 kcal/mol is attainable even for multiply ionized compounds, such as methotrexate, for which electrostatic interactions energies are very large. © 1998 John Wiley & Sons, Inc. Int J Quant Chem 69: 77–88, 1998     

Prediction of pKa shifts in proteins using a combination of molecular mechanical and continuum solvent calculations

The prediction of pKa shifts of ionizable groups in proteins is of great relevance for a number of important biological phenomena. We present an implementation of the MM-GBSA approach, which combines molecular mechanical (MM) and generalized Born (GB) continuum solvent energy terms, to the calculation of pKa values of a panel of nine proteins, including 69 individual comparisons with experiment. While applied so far mainly to the calculation of biomolecular binding free energies, we show that this method can also be used for the estimation of protein pKa shifts, with an accuracy around 1 pKa unit, even for strongly shifted residues. Our analysis reveals that the nonelectrostatic terms that are part of the MM-GBSA free energy expression are important contributors to improved prediction accuracy. This suggests that most of the previous approaches that focus only on electrostatic interactions could be improved by adding other nonpolar energy terms to their free energy expression. Interestingly, our method yields best accuracy at protein dielectric constants of epsilonint = 2-4, which is in contrast to previous approaches that peak at higher epsilonint > or = 8. An important component of our procedure is an intermediate minimization step of each protonation state involving different rotamers and tautomers as a way to explicitly model protein relaxation upon (de)protonation.