Predicting hydration free energies with a hybrid QM/MM approach: an evaluation of implicit and explicit solvation models in SAMPL4

The correct representation of solute-water interactions is essential for the accurate simulation of most biological phenomena. Several highly accurate quantum methods are available to deal with solvation by using both implicit and explicit solvents. So far, however, most evaluations of those methods were based on a single conformation, which neglects solute entropy. Here, we present the first test of a novel approach to determine hydration free energies that uses molecular mechanics (MM) to sample phase space and quantum mechanics (QM) to evaluate the potential energies. Free energies are determined by using re-weighting with the Non-Boltzmann Bennett (NBB) method. In this context, the method is referred to as QM-NBB. Based on snapshots from MM sampling and accounting for their correct Boltzmann weight, it is possible to obtain hydration free energies that incorporate the effect of solute entropy. We evaluate the performance of several QM implicit solvent models, as well as explicit solvent QM/MM for the blind subset of the SAMPL4 hydration free energy challenge. While classical free energy simulations with molecular dynamics give root mean square deviations (RMSD) of 2.8 and 2.3 kcal/mol, the hybrid approach yields an improved RMSD of 1.6 kcal/mol. By selecting an appropriate functional and basis set, the RMSD can be reduced to 1 kcal/mol for calculations based on a single conformation. Results for a selected set of challenging molecules imply that this RMSD can be further reduced by using NBB to reweight MM trajectories with the SMD implicit solvent model.     

Solutes probe hydration in specific association of cyclodextrin and adamantane

Using microcalorimetry, we follow changes in the association free energy of beta-cyclodextrin (CD) with the hydrophobic part of adamantane carboxylate (AD) due to added salt or polar (net-neutral) solutes that are excluded from the molecular interacting surfaces. Changes in binding constants with solution osmotic pressure (water activity) translate into changes in the preferential hydration upon complex formation. We find that these changes correspond to a release of 15-25 solute-excluding waters upon CD/AD association. Reflecting the preferential interaction of solute with reactants versus products, we find that changes in hydration depend on the type of solute used. All solutes used here result in a large change in the enthalpy of the CD-AD binding reaction. In one class of solutes, the corresponding entropy change is much smaller, while in the other class, the entropy change almost fully compensates the solute-specific enthalpy. For many of the solutes, the number of waters released correlates well with their effect on air-water surface tensions. We corroborate these results using vapor pressure osmometry to probe individually the hydration of reactants and products of association, and we discuss the possible interactions and forces between cosolute and hydrophobic surfaces responsible for different kinds of solute exclusion.     

Pathway dependence of the efficiency of calculating free energy and entropy of solute–solute association in water

In this study we investigate two alternative pathways to compute the free energy and the entropy of small molecule association (ΔF_assoc and ΔS_assoc) in water. The first route (direct pathway) uses thermodynamic integration as function of the distance R between the solutes. The mean force and the mean covariance of the force with the energy in solution are calculated from molecular dynamics simulation followed by integration of these quantities with respect to the reaction coordinate R. The alternative approach examined (solvation pathway) would first remove the solutes from the solution using thermodynamic integration as function of a solvation coupling parameter λ, change the solute–solute distance in vacuo and then solvate back the solute pair at the new separation distance. The system studied was a pair of CH₄ molecules in water. We investigate the influence of the CH₄–water interaction strength on the obtained ΔF_assoc and ΔS_assoc values by changing van der Waals and Coulomb interaction and evaluated the accuracy and efficiency for the two pathways. We find that the direct route seems more suitable for the calculation of free energies of hydrophobic solutes while the solvation pathway performs better when calculating entropy changes for solutes that have a stronger interaction with the solvent.     

Charge asymmetries in hydration of polar solutes

We study the solvation of polar molecules in water. The center of water's dipole moment is offset from its steric center. In common water models, the Lennard-Jones center is closer to the negatively charged oxygen than to the positively charged hydrogens. This asymmetry of water's charge sites leads to different hydration free energies of positive versus negative ions of the same size. Here, we explore these hydration effects for some hypothetical neutral solutes, and two real solutes, with molecular dynamics simulations using several different water models. We find that, like ions, polar solutes are solvated differently in water depending on the sign of the partial charges. Solutes having a large negative charge balancing diffuse positive charges are preferentially solvated relative to those having a large positive charge balancing diffuse negative charges. Asymmetries in hydration free energies can be as large as 10 kcal/mol for neutral benzene-sized solutes. These asymmetries are mainly enthalpic, arising primarily from the first solvation shell water structure. Such effects are not readily captured by implicit solvent models, which respond symmetrically with respect to charge.     

Simulations of solvation free energies and solubilities in supercritical solvents

Computer simulations are used to study solvation free energies and solubilities in supercritical solvents. Solvation free energies are calculated using the particle insertion method. The equilibrium solvent configurations required for these calculations are based on molecular dynamics simulations employing model solvent potentials previously tuned to reproduce liquid-vapor coexistence properties of the fluids Xe, C(2)H(6), CO(2), and CHF(3). Solutes are represented by all-atom potentials based on ab initio calculations and the OPLS-AA parameter set. Without any tuning of the intermolecular potentials, such calculations are found to reproduce the solvation free energies of a variety of typical solid solutes with an average accuracy of +/-2 kJmol. Further calculations on simple model solutes are also used to explore general aspects of solvation free energies in supercritical solvents. Comparisons of solutes in Lennard-Jones and hard-sphere representations of Xe show that solvation free energies and thus solubilities are not significantly influenced by solvent density fluctuations near the critical point. The solvation enthalpy and entropy do couple to these fluctuations and diverge similarly to solute partial molar volumes. Solvation free energies are also found to be little affected by the local density augmentation characteristic of the compressible regime. In contrast to solute-solvent interaction energies, which often provide a direct measure of local solvent densities, solvation free energies are remarkably insensitive to the presence of local density augmentation.     

Implicit solvent models: Combining an analytical formulation of continuum electrostatics with simple models of the hydrophobic effect

The recent development of approximate analytical formulations of continuum electrostatics opens the possibility of efficient and accurate implicit solvent models for biomolecular simulations. One such formulation (ACE, Schaefer & Karplus, J. Phys. Chem., 1996, 100:1578) is used to compute the electrostatic contribution to solvation and conformational free energies of a set of small solutes and three proteins. Results are compared to finite-difference solutions of the Poisson equation (FDPB) and explicit solvent simulations and experimental data where available. Small molecule solvation free energies agree with FDPB within 1–1.5 kcal/mol, which is comparable to differences in FDPB due to different surface treatments or different force field parameterizations. Side chain conformation free energies of aspartate and asparagine are in qualitative agreement with explicit solvent simulations, while 74 conformations of a surface loop in the protein Ras are accurately ranked compared to FDPB. Preliminary results for solvation free energies of small alkane and polar solutes suggest that a recent Gaussian model could be used in combination with analytical continuum electrostatics to treat nonpolar interactions. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 322–335, 1999     

Monte Carlo simulations of the solution structure of simple alcohols in water-acetonitrile mixtures

Monte Carlo simulations have been performed to explore the solution structure of ethyl, isopropyl, isobutyl, and tertiary butyl alcohols in pure water, pure acetonitrile, and different mixtures of the two solvents. The explicit solvent studies in NpT ensembles at T = 298 K illustrate that the solute "discriminates" the solvent's components and that the composition of the first solvation shell differs from that of the bulk solution. Since the polarizable continuum dielectric method (PCM) does not presently model the solvation of molecules with both polar and apolar sites in mixed protic solvents, we suggest a direction for further program development wherein a continuum dielectric method would accept more than one solvent and the solute sites would be solvated by user-defined solvent components. The prevailing solvation model will be determined upon the lowest free energy calculated for a particular solvation pattern of the solute having a specific conformational/tautomeric state. Characterization of equilibrium hydrogen-bond formation becomes a complicated problem that depends on the chemical properties of the solute and its conformation, as well as upon the varying nature of the first solvation shell. For example, while the number of hydrogen bonds to secondary and tertiary alcohol solutes are nearly constant in pure water and in water-acetonitrile mixtures with at least 50% water content, the number of hydrogen bonds to primary alcohols gradually decreases for most of their conformations when acetonitrile content is increased. Nonetheless, the calculations indicate that O-H...O(water) hydrogen bonds are still possible in a small fraction of the arrangements for the solution models with water content of 30% or less. The isopentene solute does not form any observable hydrogen bonds, despite having an electron-rich, double-bond site.     

Predicting hydration free energies using all-atom molecular dynamics simulations and multiple starting conformations

Molecular dynamics simulations in explicit solvent were applied to predict the hydration free energies for 23 small organic molecules in blind SAMPL2 test. We found good agreement with experimental results, with an RMS error of 2.82 kcal/mol over the whole set and 1.86 kcal/mol over all the molecules except several hydroxyl-rich compounds where we find evidence for a systematic error in the force field. We tested two different solvent models, TIP3P and TIP4P-Ew, and obtained very similar hydration free energies for these two models; the RMS difference was 0.64 kcal/mol. We found that preferred conformation of the carboxylic acids in water differs from that in vacuum. Surprisingly, this conformational change is not adequately sampled on simulation timescales, so we apply an umbrella sampling technique to include free energies associated with the conformational change. Overall, the results of this test reveal that the force field parameters for some groups of molecules (such as hydroxyl-rich compounds) still need to be improved, but for most compounds, accuracy was consistent with that seen in our previous tests.     

Compensation Between the Entropic and Enthalpic Changes Arising from Changes in Solvent–Solvent Interactions in Mixed Solvents

A theorem presented by Professor Ben-Naim (J Phys Chem 82:874–885, 1978) states that the standard state enthalpy and entropy changes arising from changes in the solvent structure that are induced by solvation of a solute cancel exactly in the standard state Gibbs energy. In this paper this is explored by consideration of the thermodynamics of transfer of electrolytes in mixed solvents, using previously developed models of the solvation process. Two cases are considered. One is random solvation, where curvatures in plots of the transfer enthalpies and entropies, which arise from changes in solvent–solvent interactions, exactly compensate in the transfer Gibbs (free) energies, which are sensibly linear with solvent composition. The second type of system are those with strong preferential solvation where it is found that the transfer Gibbs energies can be accounted for quantitatively in terms of changes in the solute–solvent interactions, with no contribution from changes in solvent–solvent interactions. The results are entirely consistent with the Ben-Naim theorem.     

Simple physics-based analytical formulas for the potentials of mean force for the interaction of amino acid side chains in water. 1. Approximate expression for the free energy of hydrophobic association based on a Gaussian-overlap model

A physics-based model is proposed to derive approximate analytical expressions for the cavity component of the free energy of hydrophobic association of spherical and spheroidal solutes in water. The model is based on the difference between the number and context of the water molecules in the hydration sphere of a hydrophobic dimer and of two isolated hydrophobic solutes. It is assumed that the water molecules touching the convex part of the molecular surface of the dimer and those in the hydration spheres of the monomers contribute equally to the free energy of solvation, and those touching the saddle part of the molecular surface of the dimer result in a more pronounced increase in free energy because of their more restricted mobility (entropy loss) and fewer favorable electrostatic interactions with other water molecules. The density of water in the hydration sphere around a single solute particle is approximated by the derivative of a Gaussian centered on the solute molecule with respect to its standard deviation. On the basis of this approximation, the number of water molecules in different parts of the hydration sphere of the dimer is expressed in terms of the first and the second mixed derivatives of the two Gaussians centered on the first and second solute molecules, respectively, with respect to the standard deviations of these Gaussians, and plausible analytical expressions for the cavity component of the hydrophobic-association energy of spherical and spheroidal solutes are introduced. As opposed to earlier hydration-shell models, our expressions reproduce the desolvation maxima in the potentials of mean force of pairs of nonpolar solutes in water, and their advantage over the models based on molecular-surface area is that they have continuous gradients in the coordinates of solute centers.     

On the transferability of hydration-parametrized continuum electrostatics models to solvated binding calculations

Using molecular mechanics force field partial atomic charges, we show the nonuniqueness of the parametrization of continuum electrostatics models with respect to solute atomic radii and interior dielectric constant based on hydration (vacuum-to-water transfer) free energy data available for small molecules. Moreover, parameter sets that are optimal and equivalent for hydration free energy calculations lead to large variations of calculated absolute and relative electrostatic binding free energies. Hence, parametrization of solvation effects based on hydration data, although a necessary condition, is not sufficient to guarantee its transferability to the calculation of binding free energies in solution.     

Calculations of solute and solvent entropies from molecular dynamics simulations

The translational, rotational and conformational (vibrational) entropy contributions to ligand-receptor binding free energies are analyzed within the standard formulation of statistical thermodynamics. It is shown that the partitioning of the binding entropy into different components is to some extent arbitrary, but an appropriate method to calculate both translational and rotational entropy contributions to noncovalent association is by estimating the configurational volumes of the ligand in the bound and free states. Different approaches to calculating solute entropies using free energy perturbation calculations, configurational volumes based on root-mean-square fluctuations and covariance matrix based quasiharmonic analysis are illustrated for some simple molecular systems. Numerical examples for the different contributions demonstrate that theoretically derived results are well reproduced by the approximations. Calculation of solvent entropies, either using total potential energy averages or van't Hoff plots, are carried out for the case of ion solvation in water. Although convergence problems will persist for large and complex simulation systems, good agreement with experiment is obtained here for relative and absolute ion hydration entropies. We also outline how solvent and solute entropic contributions are taken into account in empirical binding free energy calculations using the linear interaction energy method. In particular it is shown that empirical scaling of the nonpolar intermolecular ligand interaction energy effectively takes into account size dependent contributions to the binding free energy.     

Atomic decomposition of the protein solvation free energy and its application to amyloid-beta protein in water

We report the development of an atomic decomposition method of the protein solvation free energy in water, which ascribes global change in the solvation free energy to local changes in protein conformation as well as in hydration structure. So far, empirical decomposition analyses based on simple continuum solvation models have prevailed in the study of protein-protein interactions, protein-ligand interactions, as well as in developing scoring functions for computer-aided drug design. However, the use of continuum solvation model suffers serious drawbacks since it yields the protein free energy landscape which is quite different from that of the explicit solvent model and since it does not properly account for the non-polar hydrophobic effects which play a crucial role in biological processes in water. Herein, we develop an exact and general decomposition method of the solvation free energy that overcomes these hindrances. We then apply this method to elucidate the molecular origin for the solvation free energy change upon the conformational transitions of 42-residue amyloid-beta protein (Aβ42) in water, whose aggregation has been implicated as a primary cause of Alzheimer's disease. We address why Aβ42 protein exhibits a great propensity to aggregate when transferred from organic phase to aqueous phase.     

Comparison of two simulation methods to compute solvation free energies and partition coefficients

The thermodynamic integration (TI) and expanded ensemble (EE) methods are used here to calculate the hydration free energy in water, the solvation free energy in 1‐octanol, and the octanol‐water partition coefficient for a six compounds of varying functionality using the optimized potentials for liquid simulations (OPLS) all‐atom (AA) force field parameters and atomic charges. Both methods use the molecular dynamics algorithm as a primary component of the simulation protocol, and both have found wide applications in fields such as the calculation of activity coefficients, phase behavior, and partition coefficients. Both methods result in solvation free energies and 1‐octanol/water partition coefficients with average absolute deviations (AAD) from experimental data to within 4 kJ/mol and 0.5 log units, respectively. Here, we find that in simulations the OPLS‐AA force field parameters (with fixed charges) can reproduce solvation free energies of solutes in 1‐octanol with AAD of about half that for the solute hydration free energies using a extended simple point charge (SPC/E) model of water. The computational efficiency of the two simulation methods are compared based on the time (in nanoseconds) required to obtain similar standard deviations in the solvation free energies and 1‐octanol/water partition coefficients. By this analysis, the EE method is found to be a factor of nine more efficient than the TI algorithm. For both methods, solvation free energy calculations in 1‐octanol consume roughly an order of magnitude more CPU hours than the hydration free energy calculations. © 2012 Wiley Periodicals, Inc.     

Evaluation of the contribution to hydration of nonelectrolytes from the hydrophobic effect

A new method was suggested for estimating the hydrophobic effect of contributions to the Gibbs energies and enthalpies of hydration of hydrocarbons, inorganic gases and rare gases. In accordance with this method the hydrophobic effect contribution to the Gibbs energy was evaluated from the difference between the hydration Gibbs energy of a solute and the non hydrophobic contribution. To estimate the latter value, the known dependence connecting the Gibbs energies of solvation of a solute in a number of aprotic solvents to the Hildebrand solubility parameter for these solvents was used. The non hydrophobic contribution to the Gibbs energy of hydration was calculated for various solutes from such dependences extended to water as solvent. The Hildebrand solubility parameter for water used in the calculation was corrected for the effect of association through hydrogen bonding. This correction was made by subtraction of the water self-association enthalpy from the enthalpy of vaporization of water. The evaluated Gibbs energies of the hydrophobic effect are positive for saturated hydrocarbons, inorganic gases and rare gases and linearly depend on the solute molecular refraction. The hydrophobic contribution to the hydration enthalpies of the solutes was calculated in the same manner as was made to calculate the hydrophobic contribution to Gibbs energies of hydration. Enthalpies of the hydrophobic effect for the solutes under study are negative.     

MST continuum study of the hydration free energies of monovalent ionic species

In this study, we revisit the protocol previously proposed within the framework of the Miertus-Scrocco-Tomasi (MST) continuum model to define the cavity between the solute and solvent for predicting hydration free energies of univalent ions. The protocol relies on the use of a reduced cavity (around 10-15% smaller than the cavity used for neutral compounds) around the atom(s) bearing the formal charge. The suitability of this approach is examined here for a series of 47 univalent ions for which accurate experimental hydration free energies are available. Attention is also paid to the effect of the charge renormalization protocol used to correct uncertainties arising from the electron density located outside the solute cavity. The method presented here provides, with a minimum number of fitted parameters, reasonable estimates within the experimental error of the hydration free energy of ions (average relative error of 4.7%) and is able to reproduce solvation in water of both small and large ions.     

Solvation enthalpy and the thermodynamics of hydration of trans-cyclohexyl-1,4-diamine and cis-cyclohexyl-1,2-diamine

The enthalpy of solution of trans-cyclohexyl-1,4-diamine and cis-cyclohexyl-1,2-diamine in water was determined by calorimetry. The enthalpy of hydration was determined from this quantity and from the enthalpy of sublimation/vaporization presented in another paper by the authors. Considering the solvation process resulting from cavity creation in the solvent and variation of solute conformation transfer steps, the enthalpy corresponding to solute–solvent interaction was estimated. The entropies of solvation and interaction were calculated from the values given for the enthalpies in the present paper and those available for the Gibbs free energies.     

The importance of excluded solvent volume effects in computing hydration free energies

Continuum dielectric methods such as the Born equation have been widely used to compute the electrostatic component of the solvation free energy, DeltaG(solv)(elec), because they do not need to include solvent molecules explicitly and are thus far less costly compared to molecular simulations. All of these methods can be derived from Gauss Law of Maxwell's equations, which yields an analytical solution for the solvation free energy, DeltaG(Born), when the solute is spherical. However, in Maxwell's equations, the solvent is assumed to be a structureless continuum, whereas in reality, the near-solute solvent molecules are highly structured unlike far-solute bulk solvent. Since we have recently reformulated Gauss Law of Maxwell's equations to incorporate the near-solute solvent structure by considering excluded solvent volume effects, we have used it in this work to derive an analytical solution for the hydration free energy of an ion. In contrast to continuum solvent models, which assume that the normalized induced solvent electric dipole density P(n) is constant, P(n) mimics that observed from simulations. The analytical formula for the ionic hydration free energy shows that the Born radius, which has been used as an adjustable parameter to fit experimental hydration free energies, is no longer ill defined but is related to the radius and polarizability of the water molecule, the hydration number, and the first peak position of the solute-solvent radial distribution function. The resulting DeltaG(solv)(elec) values are shown to be close to the respective experimental numbers.