Free energy difference calculations by thermodynamic integration: Difficulties in obtaining a precise value

Free energy difference calculations have been performed by the “slow growth” method of thermodynamic integration of the AMBER 3.0 molecular dynamics program for the mutation of a conformationally restricted threonine dipeptide, N-acetyl threonyl-N-methylamide, to the corresponding alanyl dipeptide. By varying the total simulation length, it has been determined that precise free energy values are obtained only for simulations of greater than 100 ps total simulation time length. By varying the starting configurations for simulations of the same length, it has been determined that averaging the free energies obtained from shorter simulations may not give precise answers. Possible reasons for this behavior are discussed.     

Free energy perturbation and molecular dynamics calculations of copper binding to azurin

Free energy perturbation/molecular dynamics simulations have been carried out on copper/azurin systems calculating the binding affinities of copper (II) ion to azurin either in the native or in the unfolded state. In order to test the validity of the strategy adopted for the calculations and to establish what force field is suitable for these kinds of calculations, three different force fields, AMBER, CVFF, and CFF, have been alternatively used for the calculations and the results have been compared with experimental data obtained by spectroscopic titrations of copper (II)/azurin solutions and denaturation experiments. Our findings have pointed out that only CFF gives satisfactory results, thus providing a reliable tool for copper binding simulations in copper protein.     

Free energy perturbation calculations on models of active sites: Applications to adenosine deaminase inhibitors

Free energy perturbation calculations were performed to determine the free energy of binding associated with the presence of perhaps an unusual hydroxyl group in the transition state analog of nebularine, an inhibitor of the enzyme adenosine deaminase. The presence of a single hydroxyl group in this inhibitor has been found to contribute ?9.8 kcal/mol to the free energy of binding, with a 10⁸-fold increase in the binding affinity by the enzyme. In this work, we calculate the difference in solvation free energy for the 1,6-dihydropurine complex versus that of the 6-hydroxyl-1,6-dihydropurine complex to determine if this marked increase in binding affinity is attributed to an unusually hydrophobic hydroxyl group. The calculated ΔG associated for the solvation free energy is ?11.8 kcal/mol. This large change in the solvation free energy suggests that this hydroxyl is instead unusually hydrophilic and that the difference in free energy of interaction for the two inhibitors to the enzyme must be at least ca. 20 kcal/mol. Although the crystal structure for adenosine deaminase is currently not known, we attempt to mimic the nature of the active site by constructing models which simulate the enzyme-inhibitor complex. We present a first attempt at determining the change in free energy of binding for a system in which structural data for the enzyme is incomplete. To do this, we construct what we believe is a minimal model of the binding between adenosine deaminase and an inhibitor. The active site is simulated as a single charged carboxyl group which can form a hydrogen bond with the hydroxyl group of the analog. Two different carboxyl anion models are used. In the first model, the association is modeled between an acetic acid anion and the modified inhibitor. The second model consists of a hydrophobic amino acid pocket with an interior Glu residue in the active site. From these models we calculate the change in free energy of association and the overall change in free energy of binding. We calculate the free energies of interaction both in the absence and presence of water. We conclude from this that the presence of a single suitably placed-CO^?₂ group probably cannot explain the binding effect of the-OH group and that additional interactions will be found in the adenosine deaminase active site.     

Free energy perturbation calculations on parallel computers: Demonstrations of scalable linear speedup

A coarse-grain parallel implementation of the free energy perturbation (FEP) module of the AMBER molecular dynamics program is described and then demonstrated using five different molecular systems. The difference in the free energy of (aqueous) solvation is calculated for two monovalent cations ΔΔG_aq(Li⁺ Δ Cs⁺), and for the zero-sum ethane-to-ethane′ perturbation ΔΔG_aq(CH₃? methyl? X → X? methyl? CH₃), where X is a ghost methyl. The difference in binding free energy for a docked HIV-1 protease inhibitor into its ethylene mimetic is examined by mutating its fifth peptide bond, ΔG(CO? NH → CH?CH). A potassium ion (K⁺) is driven outward from the center of mass of ionophore salinomycin (SAL^?) in a potential of mean force calculation ΔG_MeOH(SAL^? · K⁺) carried out in methanol solvent. Parallel speedup obtained is linearly proportional to the number of parallel processors applied. Finally, the difference in free energy of solvation of phenol versus benzene, ΔΔG_oct(phenol → benzene), is determined in water-saturated octanol and then expressed in terms of relative partition coefficients, Δ log(P_o/w). Because no interprocessor communication is required, this approach is scalable and applicable in general for any parallel architecture or network of machines. FEP calculations run on the nCUBE/2 using 50 or 100 parallel processors were completed in clock times equivalent to or twice as fast as a Cray Y-MP. The difficulty of ensuring adequate system equilibrium when agradual configurational reorientation follows the mutation of the Hamiltonian is discussed and analyzed. The results of a successful protocol for overcoming this equilibration problem are presented. The types of molecular perturbations for which this method is expected to perform most efficiently are described. © 1994 by John Wiley & Sons, Inc.     

BLEEP—potential of mean force describing protein–ligand interactions: I. Generating potential

We have developed BLEEP (biomolecular ligand energy evaluation protocol), an atomic level potential of mean force (PMF) describing protein–ligand interactions. The pair potentials for BLEEP have been derived from high-resolution X-ray structures of protein–ligand complexes in the Brookhaven Protein Data Bank (PDB), with a careful treatment of homology. The use of a broad variety of protein–ligand structures in the derivation phase gives BLEEP more general applicability than previous potentials, which have been based on limited classes of complexes, and thus represents a significant step forward. We calculate the distance distributions in protein–ligand interactions for all 820 possible pairs that can be chosen from our set of 40 different atom types, including polar hydrogen. We then use a reverse Boltzmann methodology to convert these into energy-like pair potential functions. Two versions of BLEEP are calculated, one including and one excluding interactions between protein and water. The pair potentials are found to have the expected forms; polar and hydrogen bonding interactions show minima at short range, around 3.0 Å, whereas a typical hydrophobic interaction is repulsive at this distance, with values above 4.0 Å being preferred. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1165–1176, 1999     

Understanding the nitrate coordination to Eu3+ ions in solution by potential of mean force calculations

Coordination of nitrate anions with lanthanoid cations (Ln(3+)) in water, methanol and octanol-1 has been studied by means of molecular dynamics simulations with explicit polarization. Potential of mean force (PMF) profiles have been calculated for a mono-complex of lanthanoid nitrate (Ln(NO(3))(2+)) in these solvents using umbrella-sampling molecular dynamics. In pure water, no difference in the nitrato coordination to lanthanoids (Nd(3+), Eu(3+) and Dy(3+)) is observed, i.e. the nitrate anion prefers the monodentate coordination, which promotes the salt dissociation. Then, the influence of the nature of the solvating molecules on the nitrato coordination to Eu(3+) has been investigated. PMF profiles point out that both monodenate and bidentate coordinations are stable in neat methanol, while in neat octanol, only the bidentate one is. MD simulations of Eu(NO(3))(3) in water-octanol mixtures with different concentrations of water have been then performed and confirm the importance of the water molecules' presence on the nitrate ion's coordination mode.     

Theoretical calculations of homoconjugation equilibrium constants in systems modeling acid-base interactions in side chains of biomolecules using the potential of mean force

The potentials of mean force (PMFs) were determined for systems forming cationic and anionic homocomplexes composed of acetic acid, phenol, isopropylamine, n-butylamine, imidazole, and 4(5)-methylimidazole, and their conjugated bases or acids, respectively, in three solvents with different polarity and hydrogen-bonding propensity: acetonitrile (AN), dimethyl sulfoxide (DMSO), and water (H(2)O). For each pair and each solvent a series of umbrella-sampling molecular dynamics simulations with the AMBER force field, explicit solvent, and counterions added to maintain a zero net charge of a system were carried out and the PMF was calculated by using the Weighted Histogram Analysis Method (WHAM). Subsequently, homoconjugation-equilibrium constants were calculated by numerical integration of the respective PMF profiles. In all cases but imidazole stable homocomplexes were found to form in solution, which was manifested as the presence of contact minima corresponding to hydrogen-bonded species in the PMF curves. The calculated homoconjugation constants were found to be greater for complexes with the OHO bridge (acetic acid and phenol) than with the NHN bridge and they were found to decrease with increasing polarity and hydrogen-bonding propensity of the solvent (i.e., in the series AN > DMSO > H(2)O), both facts being in agreement with the available experimental data. It was also found that interactions with counterions are manifested as the broadening of the contact minimum or appearance of additional minima in the PMF profiles of the acetic acid-acetate, phenol/phenolate system in acetonitrile, and the 4(5)-methylimidazole/4(5)-methylimidzole cation conjugated base system in dimethyl sulfoxide.     

A comparison of methods to compute the potential of mean force.

Most processes occurring in a system are determined by the relative free energy between two or more states because the free energy is a measure of the probability of finding the system in a given state. When the two states of interest are connected by a pathway, usually called reaction coordinate, along which the free-energy profile is determined, this profile or potential of mean force (PMF) will also yield the relative free energy of the two states. Twelve different methods to compute a PMF are reviewed and compared, with regard to their precision, for a system consisting of a pair of methane molecules in aqueous solution. We analyze all combinations of the type of sampling (unbiased, umbrella-biased or constraint-biased), how to compute free energies (from density of states or force averaging) and the type of coordinate system (internal or Cartesian) used for the PMF degree of freedom. The method of choice is constraint-bias simulation combined with force averaging for either an internal or a Cartesian PMF degree of freedom.     

Calculating thermodynamic properties from perturbation theory: I. An analytic representation of square-well potential hard-sphere perturbation theory

The Barker–Henderson macroscopic compressibility approximation of the second-order perturbation term is improved by assuming that the numbers of molecules in every two neighbour shells are correlated, based upon the original assumptions. The results are better than those for the original macroscopic compressibility and local compressibility approximation, especially at high densities. A simple analytic representation of square-well potential hard-sphere perturbation theory is derived based upon this improvement. The method is tested by calculating thermodynamic properties with the four-term truncated form, and the results are in good agreement with those of Monte Carlo and Molecular Dynamics simulation.     

Free energies of hydration from thermodynamic integration: Comparison of molecular mechanics force fields and evaluation of calculation accuracy

Four commonly used molecular mechanics force fields, CHARMM22, OPLS, CVFF, and GROMOS87, are compared for their ability to reproduce experimental free energies of hydration (ΔG_hydr) from molecular dynamics (MD) simulations for a set of small nonpolar and polar organic molecules: propane, cyclopropane, dimethylether, and acetone. ΔG_hydr values were calculated by multiconfiguration thermodynamic integration for each of the different force fields with three different sets of partial atomic charges: full charges from an electrostatic potential fit (ESP), and ESP charges scaled by 0.8 and 0.6. All force fields, except for GROMOS87, give reasonable results for ΔG_hydr · if partial atomic charges of appropriate magnitude are assigned. For GROMOS87, the agreement with experiment for hydrocarbons (propane and ethane) was improved considerably by modifying the repulsive part of the carbon-water oxygen Lennard-Jones potential. The small molecules studied are related to the chemical moieties constituting camphor (C₁₀H_l6O). By invoking force-field transferability, we calculated the ΔG_hydr for camphor. With the modified GROMOS force field, a ΔG_hydr within 4 kJ/mol of the experimental value of −14.8 kJ/mol was obtained. Camphor is one of the largest molecules for which an absolute hydration free energy has been calculated by molecular simulation. The accuracy and reliability of the thermodynamic integration calculations were analyzed in detail and we found that, for ΔG_hydr calculations for the set of small molecules in aqueous solution, molecular dynamics simulations of 0.8–1.0 ns in length give an upper statistical error bound of 1.5 kJ/mol, whereas shorter simulations of 0.25 nm in length given an upper statistical error bound of 3.5 kJ/mol. © 1997 by John Wiley & Sons, Inc.     

Applicability of multi-reference many-body perturbation theory to the determination of potential energy surfaces: A model study

The recently written CI -based multi-reference many-body perturbation theory (MR-MBPT ) program package is exploited to study a simple ab initio minimum basis set model involving four hydrogen atoms in a rectangular configuration. This model was examined earlier by several authors using both coupled cluster (CC ) and finite-order MBPT approaches. Here we present the MR-MBPT results up to the 50th order and examine the effect of various shifting techniques on the convergence behavior of this approach. It is shown that in contrast with CC methods, both single and MR finite-order MBPT potential energy calculations are plagued with convergency and intruder state problems, which can be particularly severe when the latter approach is employed for non-degenerate situations.     

Toward the description of electrostatic interactions between globular proteins: potential of mean force in the primitive model

Monte Carlo simulations are used to calculate the exact potential of mean force between charged globular proteins in aqueous solution. The aim of the present paper is to study the influence of the ions of the added salt on the effective interaction between these nanoparticles. The charges of the model proteins, either identical or opposite, are either central or distributed on a discrete pattern. Contrarily to Poisson-Boltzmann predictions, attractive, and repulsive direct forces between proteins are not screened similarly. Moreover, it has been shown that the relative orientations of the charge patterns strongly influence salt-mediated interactions. More precisely, for short distances between the proteins, ions enhance the difference of the effective forces between (i) like-charged and oppositely charged proteins, (ii) attractive and repulsive relative orientations of the proteins, which may affect the selectivity of protein/protein recognition. Finally, such results observed with the simplest models are applied to a more elaborate one to demonstrate their generality.     

Exceptional thermodynamic stability of DNA duplexes modified by nonpolar base analogues is due to increased stacking interactions and favorable solvation: Correlated ab initio calculations and molecular dynamics simulations

The geometries of DNA hexamer (5'-GGAACC-3') and DNA 13-mer (5'-GCGTACACATGCG-3') have been determined by molecular dynamics (MD) simulations using an empirical force field. The central canonical base pair was replaced by a pair of nonpolar base analogues, 2,2'-bipyridyl and 3-methylisocarbostyril. The stabilization energy of the model system (model A) consisting of a central base pair (base-analogue pair) and two neighboring base pairs was determined by the RI-MP2 method using an extended aug-cc-pVDZ basis set. The geometry of the model was averaged from structures determined by MD simulations. The role of the solvent was covered by the COSMO continuum solvent model and calculations were performed for a larger model system (model B) which also contained a sugar-phosphate backbone. The total stabilization energies of the unperturbed system and the system perturbed by a base-analogue pair (model A) were comparable to the stability of both duplexes experimentally determined. This is due to large stacking interaction energy of the base-analogue self-pair which compensates for the missing hydrogen-bonding energy of the replaced adenine...thymine base pair. The selectivity of the base-analogue pair was reproduced (model B) when their desolvation energy was included with the interaction energy of both strands determined by the approximate SCC-DFTB-D method.     

Free energy profiles for penetration of methane and water molecules into spherical sodium dodecyl sulfate micelles obtained using the thermodynamic integration method combined with molecular dynamics calculations

The free energy profiles, ΔG(r), for penetration of methane and water molecules into sodium dodecyl sulfate (SDS) micelles have been calculated as a function of distance r from the SDS micelle to the methane and water molecules, using the thermodynamic integration method combined with molecular dynamics calculations. The calculations showed that methane is about 6-12 kJ mol(-1) more stable in the SDS micelle than in the water phase, and no ΔG(r) barrier is observed in the vicinity of the sulfate ions of the SDS micelle, implying that methane is easily drawn into the SDS micelle. Based on analysis of the contributions from hydrophobic groups, sulfate ions, sodium ions, and solvent water to ΔG(r), it is clear that methane in the SDS micelle is about 25 kJ mol(-1) more stable than it is in the water phase because of the contribution from the solvent water itself. This can be understood by the hydrophobic effect. In contrast, methane is destabilized by 5-15 kJ mol(-1) by the contribution from the hydrophobic groups of the SDS micelle because of the repulsive interactions between the methane and the crowded hydrophobic groups of the SDS. The large stabilizing effect of the solvent water is higher than the repulsion by the hydrophobic groups, driving methane to become solubilized into the SDS micelle. A good correlation was found between the distribution of cavities and the distribution of methane molecules in the micelle. The methane may move about in the SDS micelle by diffusing between cavities. In contrast, with respect to the water, ΔG(r) has a large positive value of 24-35 kJ mol(-1), so water is not stabilized in the micelle. Analysis showed that the contributions change in complex ways as a function of r and cancel each other out. Reference calculations of the mean forces on a penetrating water molecule into a dodecane droplet clearly showed the same free energy behavior. The common feature is that water is less stable in the hydrophobic core than in the water phase because of the energetic disadvantage of breaking hydrogen bonds formed in the water phase. The difference between the behaviors of the SDS micelles and the dodecane droplets is found just at the interface; this is caused by the strong surface dipole moment formed by sulfate ions and sodium ions in the SDS micelles.     

The application of thermodynamic perturbation theory to the calculation of excess free energy of small systems: 1. The study of the dimensional dependence of specific free energy of small droplets

Fundamentals of the application of the method of thermodynamic perturbation theory to finding the excess free energy of small objects were analyzed. The size dependence of the specific surface free energy of small droplets of a simple Lennard-Jones fluid, water, and metal melts was studied in the first and second approximations of the perturbation theory. It was established that, for all studied systems, the size dependence of the surface tension was satisfactorily described by Tolmans formula. At small particle radii, Rusanovs formula was well fulfilled.Translated from Kolloidnyi Zhurnal, Vol. 66, No. 6, 2004, pp. 844–849.Original Russian Text Copyright © 2004 by Shcherbakov, Samsonov, Bazulev.Deceased.     

Temperature dependence of three-body hydrophobic interactions: potential of mean force, enthalpy, entropy, heat capacity, and nonadditivity

Temperature-dependent three-body hydrophobic interactions are investigated by extensive constant-pressure simulations of methane-like nonpolar solutes in TIP4P model water at six temperatures. A multiple-body hydrophobic interaction is considered to be (i) additive, (ii) cooperative, or (iii) anti-cooperative if its potential of mean force (PMF) is (i) equal to, (ii) smaller than, or (iii) larger than the corresponding pairwise sum of two-methane PMFs. We found that three-methane hydrophobic interactions at the desolvation barrier are anti-cooperative at low to intermediate T, and vary from essentially additive to slightly cooperative at high T. Interactions at the contact minimum are slightly anti-cooperative over a wider temperature range. Enthalpy, entropy, and heat capacity are estimated from the computed PMFs. Contrary to the common expectation that burial of solvent-accessible nonpolar surface area always leads to a decrease in heat capacity, the present results show that the change in heat capacity upon three-methane association is significantly positive at the desolvation barrier and slightly positive at the contact minimum. This suggests that the heat capacity signature of a hydrophobic polymer need not vary uniformly nor monotonically with conformational compactness. Ramifications for protein folding are discussed.     

Hydrophobic interactions between methane and a nanoscopic pocket: three dimensional distribution of potential of mean force revealed by computer simulations

We consider a model system of methane molecule and a hemispherical, hydrophobic pocket of an 8 A radius, remaining together in aqueous environment. A spatial map of potential of mean force acting on methane molecule due to presence of pocket is constructed, based on a series of explicit solvent molecular dynamics simulations. A relation between free energy variations associated with methane translocations and accompanying changes in solvent density distribution is analyzed. A funnel-like area where free energy is diminished with respect to bulk is present over the pocket entrance and extends up to 9 A toward the bulk solvent. In order to get into the pocket, however, methane has to cross a free energy barrier, which is more prominent around the circumferential part of pocket entrance, while achieving bulklike free energy values at the very center. As a methane molecule crosses this barrier, the pocket gets completely dehydrated, which leads to "hydrophobic collapse," manifested by a sharp decrease in free energy. We find that the observed free energy changes are closely related to interactions between the methane hydration shell and the surrounding solvent. Results presented here are a continuation of our previous studies of methane-pocket systems.