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

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

Interaction of phospholipids (Lysophosphatidylethanolamines) with water and sodium cation

We have performed detailed ab initio SCF calculations on the intermolecular interaction energies for one Na⁺ ion and one water molecule with two molecular fragments, one exemplifying a phospholipid (PL) head (PLHD) and the other, a phospholipid tail (PLTL). A 6-12-1 atom-atom pair potential for the interaction of a Na⁺ ion and water with a lysophosphatidyl-ethanolamine (LPEA) was derived from these results by a fitting procedure. This fitted potential was used to obtain isoenergy maps that provide energy profiles of the Na⁺ ion and the water around the phospholipids. The interaction of the Na⁺ ion with PL, as well as the interaction of water with the PL, can be visualized from these maps, which, as expected, show regions of hydrophilicity and hydrophobicity for the water and indicate a very strong binding site for the Na⁺ ion on the phosphate. It appears to be a stationary site that would limit the Na⁺ ion mobility. This binding site is located near the double-bonded oxygen atom of the phosphate group; its binding energy for Na⁺ is 67 kcal/mol. On the other hand the NH⁺ group of PLHD ahows strong electrostatic repulsion of Na⁺ while interacting with water with a binding energy of 13 kcal/mol. This potential energy well region for water is separated from another of similar depth near the phosphate by a barrier and both regions are expected to act as binding sites for water.     

Calculation of binding affinities of HIV-1 RT and β-secretase inhibitors using the linear interaction energy method with explicit and continuum solvation approaches

The linear interaction energy (LIE) approach has been applied to estimate the binding free energies of representative sets of HIV-1 RT and β-Secretase inhibitors, using both molecular dynamics (MD) and tethered energy minimization sampling protocols with the OPLS-AA potential, using a range of solvation methodologies. Generalized Born (GB), ‘shell’ and periodic boundary condition (PBC) solvation were used, the latter with reaction field (RF) electrostatics. Poisson-Boltzmann (PB) and GB continuum electrostatics schemes were applied to the simulation trajectories for each solvation type to estimate the electrostatic ligand-water interaction energy in both the free and bound states. Reasonable agreement of the LIE predictions was obtained with respect to experimental binding free energy estimates for both systems: for instance, ‘PB’ fits on MD trajectories carried out with PBC solvation and RF electrostatics led to models with standard errors of 1.11 and 1.03 kcal mol⁻¹ and coefficients of determination, r ² of 0.76 and 0.75 for the HIV-1 RT and β-Secretase sets. However, it was also found that results from MD sampling using PBC solvation provided only slightly better fits than from simulations using shell or Born solvation or tethered energy minimization sampling. Figure Evolution of the running averages for compound H11 (binding to HIV-1RT) of the bound state ligand-water and ligand-protein interaction energies. The ligand-water electrostatic terms are twice the corresponding GB and PB electrostatic solvation free energies. The ligand-receptor van der Waals and Coulombic interaction energies are also shown, in addition to the ligand-water van der Waals interaction term. The terms were calculated (without application of a cut-off) from a trajectory sampled under PBC solvation with reaction field electrostatics Electronic supplementary material The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

The nature of intermolecular interactions between aromatic amino acid residues

The nature of intermolecular interactions between aromatic amino acid residues has been investigated by a combination of molecular dynamics and ab initio methods. The potential energy surface of various interacting pairs, including tryptophan, phenilalanine, and tyrosine, was scanned for determining all the relevant local minima by a combined molecular dynamics and conjugate gradient methodology with the AMBER force field. For each of these minima, single-point correlated ab initio calculations of the binding energy were performed. The agreement between empirical force field and ab initio binding energies of the minimum energy structures is excellent. Aromatic-aromatic interactions can be rationalized on the basis of electrostatic and van der Waals interactions, whereas charge transfer or polarization phenomena are small for all intermolecular complexes and, particularly, for stacked structures. Proteins 2002;48:117-125.     

Ab Initio Quantum Mechanics Analysis of Imidazole C-H…O Water Hydrogen Bonding and a Molecular Mechanics Forcefield Correction

Abstract

While it is well established that classical hydrogen bonds play an important role in enzyme structure, function and dynamics, the role of weaker, but ‘activated’ C-H donor hydrogen bonds is poorly understood. The most important such case involves histidine which often plays a direct role in enzyme catalysis and possesses the most acidic C-H donor group of the standard amino acids. In the present study, we obtained optimized geometries and hydrogen bond interaction energies for C-H…O hydrogen bonded complexes between methane, ethylene, benzene, acetylene, and imidazole with water at the MP2-FC/6-31++G(2d,2p) and MP2-FC/aug-cc-pVDZ//MP2-FC/6-31++G(2d,2p) levels of theory. A strong linear relationship is obtained between the stability of the various hydrogen bonded complexes and both separation distances for H…0 and C—O. In general, these calculations indicate that C-H…0 interactions can be classified as hydrogen bonding interactions, albeit significantly weaker than the classical hydrogen bonds, but significantly stronger than just van der Waals interactions. For instance, while the electronic energy of stabilization at the MP2-FC/aug-cc-pVDZ//MP2-FC/6-31++G(2d,2p) level of theory of a water C-H…O water hydrogen bond is 4.36 kcal/mol more stable than the methane C-H…O water interaction, the water-water hydrogen bond is only 2.06 kcal/mol more stable than the imidazole C_e?H…O water hydrogen bond. Neglecting this latter hydrogen bonding interaction is obviously unacceptable. We next compare the potential energy surfaces for the imidazole C_e?H…O water and imidazole N_d?H…O hydrogen bonded complexes computed at the MP2/6-31++G(2d,2p) level of theory with the potential energy surface computed using the AMBER molecular mechanics program and forcefields. While the Weiner et al and Cornell et al AMBER forcefields reasonably account for the imidazole N-H…O water interaction, these forcefields do not adequately account for the imidazole C_e?H…O water hydrogen bond. A forcefield modification is offered that results in excellent agreement between the ab initio and molecular mechanics geometry and energy for this C-H…O hydrogen bonded complex.     

A CFF 91-based Continuum Solvation Model: Solvation Free Energies of Small Organic Molecules and Conformations of the Alanine Dipeptide in Solution

Abstract

For molecular mechanics simulations of solvated molecules, it is important to use a consistent approach for calculating both the force field energy and the solvation free energy. A continuum solvation model based upon the atomic charges provided with the CFF91 force field is derived. The electrostatic component of the solvation free energy is described by the Poisson-Bolzmann equation while the nonpolar comonent of the solvation energy is assumed to be proportional to the solvent accessible surface area of the solute. Solute atomic radii used to describe the interface between the solute and solvent are fitted to reproduce the energies of small organic molecules. Data for 140 compounds are presented and compared to experiment and to the results from the well-characterized quantum mechanical solvation model AM1-SM2. In particular, accurate results are obtained for amino acid neutral analogues (mean unsigned error of 0.3 kcal/mol). The conformational energetics of the solvated alanine dipeptide is discussed.     

Analysis of serine proteinase-inhibitor interaction by alanine shaving

We analyzed the energetic importance of residues surrounding the hot spot (the P(1) position) of bovine pancreatic trypsin inhibitor (BPTI) in interaction with two proteinases, trypsin and chymotrypsin, by a procedure called molecular shaving. One to eight residues of the structural epitope, composed of two extended and exposed loops, were mutated to alanine(s). Although truncation of the side chains of residues surrounding the P(1) position to methyl groups caused a decrease in Delta G(den) values up to 6.4 kcal mole(-1), it did not influence the overall conformation of the inhibitor. We found that the replacement of up to six residues with alanines was fully additive at the level of protein stability. To analyze the influence of the structural epitope on the association energy, we determined association constants for BPTI variants and both enzymes and applied the additivity analysis. Shaving of two binding loops led to a progressive drop in the association energy, more pronounced for trypsin (decrease up to 9.6 kcal mole(-1)) than chymotrypsin (decrease up to 3.5 kcal mole(-1)). In the case of extensively mutated variants interacting with chymotrypsin, the association energies agreed very well with the values calculated from single mutational effects. However, when P(1)-neighboring residues were shaved to alanine(s), their contribution to the association energy was not fully removed because of the presence of methyl groups and main chain-main chain intermolecular hydrogen bonds. Moreover, the hot spot had a different contribution to the complex stability in the fully shaved BPTI variant compared with the wild type, which was caused by perturbations of the P(1)-S(1) electrostatic interaction.     

Predicting the rate enhancement of protein complex formation from the electrostatic energy of interaction

The rate of association of proteins is dictated by diffusion, but can be enhanced by favorable electrostatic forces. Here the relationship between the electrostatic energy of interaction, and the kinetics of protein-complex formation was analyzed for the protein pairs of: hirudin-thrombin, acetylcholinesterase-fasciculin and barnase-barstar, and for a panel of point mutants of these proteins. Electrostatic energies of interaction were calculated as the difference between the electrostatic energy of the complex and the sum of the energies of the two individual proteins, using the computer simulation package DelPhi. Calculated electrostatic energies of interaction were compared to experimentally determined rates of association. One kcal/mol of Coulombic interaction energy increased the rate of association by a factor of 2.8, independent of the protein-complex or mutant analyzed. Electrostatic energies of interaction were also determined from the salt dependence of the association rate constant, using the same basic equation as for the theoretical calculation. A Br?nsted analysis of the electrostatic energies of interactions plotted versus experimentally determined ln(rate)s of association shows a linear relation between the two, with a beta value close to 1. This is interpreted as the energy of the transition state varies according to the electrostatic interaction energy, fitting a two state model for the association reaction. Calculating electrostatic rate enhancement from the electrostatic interaction energy can be used as a powerful tool to design protein complexes with altered rates of association and affinities.     

Prediction of the binding energy for small molecules,peptides and proteins

A fast and reliable evaluation of the binding energy from a single conformation of a molecular complex is an important practical task. Knowledge‐based scoring schemes may not be sufficiently general and transferable, while molecular dynamics or Monte Carlo calculations with explicit solvent are too computationally expensive for many applications. Recently, several empirical schemes using finite difference Poisson–Boltzmann electrostatics to predict energies for particular types of complexes were proposed. Here, an improved empirical binding energy function has been derived and validated on three different types of complexes: protein–small ligand, protein–peptide and protein–protein. The function uses the boundary element algorithm to evaluate the electrostatic solvation energy. We show that a single set of parameters can predict the relative binding energies of the heterogeneous validation set of complexes with 2.5 kcal/mol accuracy. We also demonstrate that global optimization of the ligand and of the flexible side‐chains of the receptor improves the accuracy of the evaluation. Copyright © 1999 John Wiley & Sons, Ltd.     

Modeling interactions between a β‐O‐4 type lignin model compound and 1‐allyl‐3‐methylimidazolium chloride ionic liquid

Aiming at understanding the molecular mechanism of the lignin dissolution in imidazolium‐based ionic liquids (ILs), this work presents a combined quantum chemistry (QC) calculation and molecular dynamics (MD) simulation study on the interaction of the lignin model compound, veratrylglycerol‐β‐guaiacyl ether (VG) with 1‐allyl‐3‐methylimidazolium chloride ([Amim]Cl). The monomer of VG is shown to feature a strong intramolecular hydrogen bond, and its dimer is indicated to present important π‐π stacking and intermolecular hydrogen bonding interactions. The interactions of both the cation and anion of [Amim]Cl with VG are shown to be stronger than that between the two monomers, indicating that [Amim]Cl is capable of dissolving lignin. While Cl^– anion forms a hydrogen‐bonded complex with VG, the imidazolium cation interacts with VG via both the π‐π stacking and intermolecular hydrogen bonding. The calculated interaction energies between VG and the IL or its components (the cation, anion, and ion pair) indicate the anion plays a more important role than the cation for the dissolution of lignin in the IL. Theoretical results provide help for understanding the molecular mechanism of lignin dissolution in imidazolium‐based IL. The theoretical calculations on the interaction between the lignin model compound and [Amim]Cl ionic liquid indicate that the anion of [Amim]Cl plays a more important role for lignin dissolution although the cation also makes a substantial contribution.     

Analysis of intermolecular interaction energy inputs in benzene-imidazole and imidazole-imidazole systems in parallel displaced and T-configuration

Intermolecular interactions in several dimer aromatic systems were analyzed to determine how various energy contributions (electrostatic, exchange, repulsion, and polarization) change depending on the value of monomers separation. Different contributions to the intermolecular energy interactions between imidazole-imidazole and benzene-imidazole dimers are studied using the aug-cc-pVDZ basis set in the framework of ab initio Hartree-Fock and second-order Møller-Plesset perturbation theory methods. Special attention is paid to the exchange and dispersion energy binding contributions.     

Model assembly study of the ligand binding by p-hydroxybenzoate hydroxylase: Correlation between the calculated binding energies and the experimental dissociation constants

The energies of binding of seven ligands by p-hydroxybenzoate hydroxylase (PHBH) were calculated theoretically. Direct enzyme–ligand interaction energies were calculated using the ab initio quantum mechanical model assembly of the active site at the 3-21G level. Solvation energies of the ligands needed in the evaluation of the binding energies were calculated with the semiempirical AM1–SM2 method and the long-range electrostatic interaction energies between the ligands and the protein matrix classically using the static charge distributions of the ligands and the protein. Energies for proton-transfer between the ligands OH or SH substituent at position 4 and the active-site tyrosine within the ab initio model assemblies were calculated and compared to the corresponding pK_as in aqueous solution. Excluding 3,4-dihydroxybenzoate, the natural product of PHBH, a linear relationship between the calculated binding energies and the experimental binding free energies was found with a correlation coefficient of 0.90. Contributions of the direct enzyme–ligand interaction energies, solvation energies and the long-range electrostatic interaction energies to the calculated binding energies were analyzed. The proton-transfer energies of the ligands with substituents ortho to the ionized OH were found to be perturbed less in the model calculations than the energies of their meta isomers as deduced from the corresponding pK_as. © 1995 Wiley-Liss, Inc.     

Interactions in pseudorotoxanes based on crown ether-secondary ammonium motifs. A theoretical study

A theoretical analysis of the nature of the interactions in dibenzo[24]crown-8 (DB24C8)-n-dibutylammonium (DBM)—pseudorotaxane complex at the MP2 and DFT levels shows that the main contribution to the binding energy is the electrostatic interaction with moderate (20–25%) correlation stabilization. The total binding energy in the DB24C8-DBM complex represents a sum of the binding energies of two NH–O and one CH–O hydrogen bonds and the latter constitutes about 25% of the total interaction energy, giving the total binding energy of −41.2 kcal mol⁻¹ at the BHandHLYP/6-311++G** level. Deprotonation of the DB24C8-DBM complex reduces the binding energy by some 50 kcal mol⁻¹, giving metastable complexes DB24C8-DBA-1 or DB24C8-DBA-2, which will dissociate to give free crown ether and n-dibutylamine because of the strong exchange repulsion that prevails in neutral complexes. Figure Formation of DB24C8-DBM pseudorotoxane complex     

Insights into protein-protein binding by binding free energy calculation and free energy decomposition for the Ras-Raf and Ras-RalGDS complexes

Absolute binding free energy calculations and free energy decompositions are presented for the protein-protein complexes H-Ras/C-Raf1 and H-Ras/RalGDS. Ras is a central switch in the regulation of cell proliferation and differentiation. In our study, we investigate the capability of the molecular mechanics (MM)-generalized Born surface area (GBSA) approach to estimate absolute binding free energies for the protein-protein complexes. Averaging gas-phase energies, solvation free energies, and entropic contributions over snapshots extracted from trajectories of the unbound proteins and the complexes, calculated binding free energies (Ras-Raf: -15.0(+/-6.3)kcal mol(-1); Ras-RalGDS: -19.5(+/-5.9)kcal mol(-1)) are in fair agreement with experimentally determined values (-9.6 kcal mol(-1); -8.4 kcal mol(-1)), if appropriate ionic strength is taken into account. Structural determinants of the binding affinity of Ras-Raf and Ras-RalGDS are identified by means of free energy decomposition. For the first time, computationally inexpensive generalized Born (GB) calculations are applied in this context to partition solvation free energies along with gas-phase energies between residues of both binding partners. For selected residues, in addition, entropic contributions are estimated by classical statistical mechanics. Comparison of the decomposition results with experimentally determined binding free energy differences for alanine mutants of interface residues yielded correlations with r(2)=0.55 and 0.46 for Ras-Raf and Ras-RalGDS, respectively. Extension of the decomposition reveals residues as far apart as 25A from the binding epitope that can contribute significantly to binding free energy. These "hotspots" are found to show large atomic fluctuations in the unbound proteins, indicating that they reside in structurally less stable regions. Furthermore, hotspot residues experience a significantly larger-than-average decrease in local fluctuations upon complex formation. Finally, by calculating a pair-wise decomposition of interactions, interaction pathways originating in the binding epitope of Raf are found that protrude through the protein structure towards the loop L1. This explains the finding of a conformational change in this region upon complex formation with Ras, and it may trigger a larger structural change in Raf, which is considered to be necessary for activation of the effector by Ras.     

Calculations of electrostatic energies in proteins. The energetics of ionized groups in bovine pancreatic trypsin inhibitor

The importance of including different energy contributions in calculations of electrostatic energies in proteins is examined by calculating the intrinsic pKa values of the acidic groups of bovine pancreatic trypsin inhibitor. It appears that such calculations provide a powerful and revealing test; the relevant solvation energies of the ionized acids are of the order of -70 kcal/mol (1 cal = 4.184 J), and microscopic calculations that do not attempt to simulate the complete protein dielectric effect (including the surrounding solvent) can underestimate the solvation energy by as much as 50 kcal/mol. Reproducing correctly, by the same set of parameters, the solvation energies of ionized acids in different sites of a protein cannot be accomplished by including only part of the key energy contributions. The problems associated with macroscopic calculations are also considered and illustrated by the specific case of bovine pancreatic trypsin inhibitor. A promising approach is shown to be provided by a refinement of the previously developed Protein Dipoles Langevin Dipoles model. This model seems to represent consistently the microscopic dielectric of the protein and the surrounding water molecules. The model overcomes the problems associated with the macroscopic models (by treating explicitly the solvent molecules) and avoids the convergence problems associated with all-atom solvent models (by treating the average solvent polarization rather than averaging the actual polarization energy). This paper describes in detail the actual implementation of the model and examines its performance in evaluating intrinsic pKa values. Preliminary microscopic considerations of charge-charge interactions are presented.     

Strong interaction between disulfide derivatives and aromatic groups in peptides and proteins

The intermolecular interaction energy of complexes of dimethyldisulfide with benzene and cyclohexane, respectively, was computed as function of the relative distance and orientation within each pair of molecules. The energy of the most stable orientation of the dimethyldisulfide-cyclohexane complex is ?2.57 kcal/mol, while that of the most stable orientation of the dimethyldisulfide-benzene complex is ?3.33 kcal/mol. The energy difference of ～0.8 kcal/mol is due to favorable specific nonbonded interactions between the sulfur atoms and the atoms of the aromatic ring. Proper parameterization of empirical interatomic energies, used in computations in this laboratory, accounts for these interactions without the need for a $\text{special}$ sulfur-aromatic potential energy function.     

Calculation of the magnitude and orientation of electrostatic interactions between small aromatic rings in peptides and proteins: implications for angiotensin II

Electrostatic interactions between aromatic rings are believed to contribute to the intramolecular structuring and biological function of peptides and proteins. The modes of interaction of benzene with 1) benzene (Phe:Phe), 2) imidazole (Phe:His), and 3) imidazole anion (Phe:His), were calculated using all-electron ab initio wavefunctions. In all cases parallel-plate stackings were found to be purely repulsive, whereas perpendicular-plate geometries were attractive with interaction energies of -0.82 (Phe:Phe), -1.19 (Phe:His) and -3.39 (Phe:His) kcal/mole. These data show that small ring interactions in peptides and proteins will prefer perpendicular-plate geometry. For the proposed His:Phe interaction in angiotensin II, the attraction will be three times greater when the imidazole ring carries a negative charge.