Predicting aqueous free energies of solvation as functions of temperature

This work introduces a model, solvation model 6 with temperature dependence (SM6T), to predict the temperature dependence of aqueous free energies of solvation for compounds containing H, C, and O in the range 273-373 K. In particular, we extend solvation model 6 (SM6), which was previously developed (Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Chem. Theory Comput. 2005, 1, 1133) for predicting aqueous free energies of solvation at 298 K, to predict the variation of the free energy of solvation relative to 298 K. Also, we describe the database of experimental aqueous free energies of solvation for compounds containing H, C, and O that was used to parametrize and test the new model. SM6T partitions the temperature dependence of the free energy of solvation into two components: the temperature dependence of the bulk electrostatic contribution to the free energy of solvation, which is computed using the generalized Born equation, and the temperature dependence of first-solvation-shell effects which is modeled using a parametrized solvent-exposed surface-area-dependent term. We found that SM6T predicts the temperature dependence of aqueous free energies of solvation with a mean unsigned error of 0.08 kcal/mol over our entire database, whereas using the experimental value at 298 K produces a mean unsigned error of 0.53 kcal/mol.     

Extension of a temperature-dependent aqueous solvation model to compounds containing nitrogen, fluorine, chlorine, bromine, and sulfur

Most methods for predicting free energies of solvation have been developed or validated exclusively for room temperature. Recently, we developed a model called SM6T for predicting aqueous solvation free energies as a function of temperature for solutes composed of C, H, or O, and here we present solvation model 8 with temperature dependence (SM8T) for predicting the temperature dependence of aqueous free energies of solvation for compounds containing H, C, N, O, F, S, Cl, and Br in the range 273-373 K. We also describe the database of experimental aqueous free energies of solvation used to parametrize the model. SM8T partitions the temperature dependence of the free energy of solvation into two components: the temperature dependence of the bulk electrostatic contribution to the free energy of solvation, which is computed using the generalized Born equation, and the temperature dependence of first-solvation-shell effects, which is modeled by terms proportional to the solvent-exposed surface areas of atoms in functional groups determined entirely by geometry. SM8T predicts the temperature dependence of aqueous free energies of solvation with a mean unsigned error of 0.08 kcal/mol over a database of 4403 measurements on 348 compounds at various temperatures. We also discuss the accuracy of SM8T for predicting the temperature dependence of aqueous free energies of solvation for ions and present free energies of solvation as a function of temperature for two sample ions.     

Aqueous solvation free energies of ions and ion-water clusters based on an accurate value for the absolute aqueous solvation free energy of the proton

Thermochemical cycles that involve pKa, gas-phase acidities, aqueous solvation free energies of neutral species, and gas-phase clustering free energies have been used with the cluster pair approximation to determine the absolute aqueous solvation free energy of the proton. The best value obtained in this work is in good agreement with the value reported by Tissandier et al. (Tissandier, M. D.; Cowen, K. A.; Feng, W. Y.; Gundlach, E.; Cohen, M. J.; Earhart, A. D.; Coe, J. V. J. Phys. Chem. A 1998, 102, 7787), who applied the cluster pair approximation to a less diverse and smaller data set of ions. We agree with previous workers who advocated the value of -265.9 kcal/mol for the absolute aqueous solvation free energy of the proton. Considering the uncertainties associated with the experimental gas-phase free energies of ions that are required to use the cluster pair approximation as well as analyses of various subsets of data, we estimate an uncertainty for the absolute aqueous solvation free energy of the proton of no less than 2 kcal/mol. Using a value of -265.9 kcal/mol for the absolute aqueous solvation free energy of the proton, we expand and update our previous compilation of absolute aqueous solvation free energies; this new data set contains conventional and absolute aqueous solvation free energies for 121 unclustered ions (not including the proton) and 147 conventional and absolute aqueous solvation free energies for 51 clustered ions containing from 1 to 6 water molecules. When tested against the same set of ions that was recently used to develop the SM6 continuum solvation model, SM6 retains its previously determined high accuracy; indeed, in most cases the mean unsigned error improves when it is tested against the more accurate reference data.     

The solvation, partitioning, hydrogen bonding, and dimerization of nucleotide bases: a multifaceted challenge for quantum chemistry

We present M06-2X density functional calculations of the chloroform/water partition coefficients of cytosine, thymine, uracil, adenine, and guanine and calculations of the free energies of association of selected unsubstituted and alkylated nucleotide base pairs in chloroform and water. Both hydrogen bonding and π-π stacking interactions are considered. Solvation effects are treated using the continuum solvent models SM8, SM8AD, and SMD, including geometry optimization in solution. Comparison of theoretical results with available experimental data indicates that all three of these solvation models predict the chloroform-water partition coefficients for the studied nucleobases qualitatively well, with mean unsigned errors in the range of 0.4-1.3 log units. All three models correctly predict the preference for hydrogen bonding over stacking for nucleobase pairs solvated in chloroform, and SM8, SM8AD, and SMD show similar accuracy in predicting the corresponding free energies of association. The agreement between theory and experiment for the association free energies of the dimers in water is more difficult to assess, as the relevant experimental data are indirect. Theory predicts that the stacking interaction of nucleobases in water is more favorable than hydrogen bonding for only two out of three tested hetero-dimers.     

Performance of SM8 on a test to predict small-molecule solvation free energies

The SM8 quantum mechanical aqueous continuum solvation model is applied to a 17-molecule test set proposed by Nicholls et al. (J. Med. Chem. 2008, 51, 769) to predict free energies of solvation. With the M06-2X density functional, the 6-31G(d) basis set, and CM4M charge model, the root-mean-square error (RMSE) of SM8 is 1.08 kcal mol(-1) for aqueous geometries and 1.14 kcal mol(-1) for gas-phase geometries. These errors compare favorably with optimal explicit and continuum models reported by Nicholls et al., having RMSEs of 1.33 and 1.87 kcal mol(-1), respectively. Other models examined by these workers had RMSEs of 1.5-2.6 kcal mol(-1). We also explore the use of other density functionals and charge models with SM8 and the RMSE increases to 1.21 kcal mol(-1) for mPW1/CM4 with gas-phase geometries, to 1.50 kcal mol(-1) for M06-2X/CM4 with gas-phase geometries, and to 1.27-1.64 kcal mol(-1) with three different models at B3LYP gas-phase geometries.     

Continuum solvation models in the linear interaction energy method

The linear interaction energy (LIE) method in combination with two different continuum solvent models has been applied to calculate protein-ligand binding free energies for a set of inhibitors against the malarial aspartic protease plasmepsin II. Ligand-water interaction energies are calculated from both Poisson-Boltzmann (PB) and Generalized Born (GB) continuum models using snapshots from explicit solvent simulations of the ligand and protein-ligand complex. These are compared to explicit solvent calculations, and we find close agreement between the explicit water and PB solvation models. The GB model overestimates the change in solvation energy, and this is caused by consistent underestimation of the effective Born radii in the protein-ligand complex. The explicit solvent LIE calculations and LIE-PB, with our standard parametrization, reproduce absolute experimental binding free energies with an average unsigned error of 0.5 and 0.7 kcal/mol, respectively. The LIE-GB method, however, requires a constant offset to approach the same level of accuracy.     

Comparative assessment of computational methods for the determination of solvation free energies in alcohol‐based molecules

The determination of differences in solvation free energies between related drug molecules remains an important challenge in computational drug optimization, when fast and accurate calculation of differences in binding free energy are required. In this study, we have evaluated the performance of five commonly used polarized continuum model (PCM) methodologies in the determination of solvation free energies for 53 typical alcohol and alkane small molecules. In addition, the performance of these PCM methods, of a thermodynamic integration (TI) protocol and of the Poisson–Boltzmann (PB) and generalized Born (GB) methods, were tested in the determination of solvation free energies changes for 28 common alkane‐alcohol transformations, by the substitution of an hydrogen atom for a hydroxyl substituent. The results show that the solvation model D (SMD) performs better among the PCM‐based approaches in estimating solvation free energies for alcohol molecules, and solvation free energy changes for alkane‐alcohol transformations, with an average error below 1 kcal/mol for both quantities. However, for the determination of solvation free energy changes on alkane‐alcohol transformation, PB and TI yielded better results. TI was particularly accurate in the treatment of hydroxyl groups additions to aromatic rings (0.53 kcal/mol), a common transformation when optimizing drug‐binding in computer‐aided drug design. © 2013 Wiley Periodicals, Inc.     

Calculation of the solvation free energy of neutral and ionic molecules in diverse solvents

The solvation free energy density (SFED) model was modified to extend its applicability and predictability. The parametrization process was performed with a large, diverse set of solvation free energies that included highly polar and ionic molecules. The mean absolute error for 1200 solvation free energies of the 379 neutral molecules in 9 organic solvents and water was 0.40 kcal/mol, and for 90 hydration free energies of ions was 1.7 kcal/mol. Overall, the calculated solvation free energies of a wide range of solute functional groups in diverse solvents were consistent with experimental data.     

Extension of the platform of applicability of the SM5.42R universal solvation model

We present eight new parameterizations of the SM5.42R solvation model: in particular we present parameterizations for HF/MIDI!, HF/6-31G^*, HF/6-31+G^*, HF/cc-pVDZ, AM1, PM3, BPW91/MIDI!, and B3LYP/MIDI!. Two of the new cases are parameterized using the reaction-field operator presented previously, and six of the new cases are parameterized with a simplified reaction-field operator; results obtained by the two methods are compared for selected examples. For a training set of 2135 data for 275 neutral solutes containing H, C, N, O, F, S, P, Cl, Br, and I in 91 solvents (water and 90 nonaqueous solvents), seven of the eight new parameterizations give mean unsigned errors in the range 0.43–0.46 kcal/mol, and the eighth – for a basis set containing diffuse functions – gives a mean unsigned error of 0.53 kcal/mol. The mean unsigned error for 49 ionic solutes (containing the same elements) in water is 3.5–3.9 kcal/mol for the Hartree–Fock, Becke–Perdew–Wang-1991 and Becke three-parameter Lee–Yang–Parr cases and 4.1 and 4.0 kcal/mol for parameterized model 3 and Austin model 1, respectively. The methods are tested for sensitivity of solvation free energies to geometry and for predicting partition coefficients of carbonates, which were not included in the training set. Received: 24 November 1998 / Accepted: 31 December 1998 / Published online: 7 June 1999     

Universal solvation model based on conductor‐like screening model

Atomic surface tensions are parameterized for use with solvation models in which the electrostatic part of the calculation is based on the conductor‐like screening model (COSMO) and the semiempirical molecular orbital methods AM1, PM3, and MNDO/d. The convergence of the calculated polarization free energies with respect to the numerical parameters of the electrostatic calculations is first examined. The accuracy and precision of the calculated values are improved significantly by adjusting two parameters that control the segmentation of the solvent‐accessible surface that is used for the calculations. The accuracy of COSMO calculations is further improved by adopting an optimized set of empirical electrostatic atomic radii. Finally, the electrostatic calculation is combined with SM5‐type atomic surface tension functionals that are used to compute the nonelectrostatic portions of the solvation free energy. All parameterizations are carried out using rigid (R) gas‐phase geometries; this combination (SM5‐type surface tensions, COSMO electrostatics, and rigid geometries) is called SM5CR. Six air–water and 76 water–solvent partition coefficients are added to the training set of air–solvent data points previously used to parameterize the SM5 suite of solvation models, thereby bringing the total number of data points in the training set to 2266. The model yields free energies of solvation and transfer with mean unsigned errors of 0.63, 0.59, and 0.61 kcal/mol for AM1, PM3, and MNDO/d, respectively, over all 2217 data points for neutral solutes in the training set and mean unsigned errors of 3.0, 2.7, and 3.1 kcal/mol, respectively, for 49 data points for the ions. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 340–366, 2000     

Calculating pKa values for substituted phenols and hydration energies for other compounds with the first‐order fuzzy‐border continuum solvation model

We have computed pK_a values for 11 substituted phenol compounds using the continuum Fuzzy‐Border (FB) solvation model. Hydration energies for 40 other compounds, including alkanes, alkenes, alkynes, ketones, amines, alcohols, ethers, aromatics, amides, heterocycles, thiols, sulfides, and acids have been calculated. The overall average unsigned error in the calculated acidity constant values was equal to 0.41 pH units and the average error in the solvation energies was 0.076 kcal/mol. We have also reproduced pK_a values of propanoic and butanoic acids within about 0.1 pH units from the experimental values by fitting the solvation parameters for carboxylate ion carbon and oxygen atoms. The FB model combines two distinguishing features. First, it limits the amount of noise which is common in numerical treatment of continuum solvation models by using fixed‐position grid points. Second, it uses either second‐ or first‐order approximation for the solvent polarization, depending on a particular implementation. These approximations are similar to those used for solute and explicit solvent fast polarization treatment which we developed previously. This article describes results of using the first‐order technique. This approximation places the presented methodology between the Generalized Born and Poisson‐Boltzmann continuum solvation models with respect to their accuracy of reproducing the many‐body effects in modeling a continuum solvent. © 2012 Wiley Periodicals, Inc.     

The effect of aqueous solvation upon alpha-helix formation for polyalanines

The incremental free energies of aqueous solution for acetyl(ala)NNH2 in its extended unfolded and alpha-helical conformations are compared using the SM5.2 solvation method of Cramer and Truhlar. A combination of density functional theory (DFT) at the B3LYP/D95(d,p) and AM1 has been employed using the ONIOM method. The incremental solvation energies of alpha-helical structures are very similar for both ONIOM and AM1 optimized structures as these structures do not significantly change upon solution. However, the conformations of the unfolded peptides change from extended beta-strand to polyproline II conformations upon aqueous solution. The incremental solvation free energy per residue of the polyproline II structure is about 2 kcal/mol/residue greater than that for the alpha-helix, representing an upper limit for the difference between the solvation energies. However, most of this difference disappears when the energy required to distort the optimized gas-phase extended beta-strand structure to the optimized polyproline II solution structure is included in the analysis, leaving an estimated difference in incremental solvation free energy of 0.3-0.5 kcal/mol favoring the unfolded structure. The solution structure sacrifices the stability derived from the intramolecular C5 H-bonds for more favorable interactions with the aqueous solvent.     

Accurate pK(a) calculations for carboxylic acids using complete basis set and Gaussian-n models combined with CPCM continuum solvation methods

Complete Basis Set and Gaussian-n methods were combined with CPCM continuum solvation methods to calculate pK(a) values for six carboxylic acids. An experimental value of -264.61 kcal/mol for the free energy of solvation of H(+), DeltaG(s)(H(+)), was combined with a value for G(gas)(H(+)) of -6.28 kcal/mol to calculate pK(a) values with Cycle 1. The Complete Basis Set gas-phase methods used to calculate gas-phase free energies are very accurate, with mean unsigned errors of 0.3 kcal/mol and standard deviations of 0.4 kcal/mol. The CPCM solvation calculations used to calculate condensed-phase free energies are slightly less accurate than the gas-phase models, and the best method has a mean unsigned error and standard deviation of 0.4 and 0.5 kcal/mol, respectively. The use of Cycle 1 and the Complete Basis Set models combined with the CPCM solvation methods yielded pK(a) values accurate to less than half a pK(a) unit.     

Density-functional theory and hybrid density-functional theory continuum solvation models for aqueous and organic solvents: universal SM5.43 and SM5.43R solvation models for any fraction of Hartree-Fock exchange

Hybrid density functional theory, which is a combined Hartree–Fock and density functional method, provides a simple but effective way to incorporate nonlocal exchange effects and static and dynamical correlation energy into an orbital-based theory with affordable computational cost for many important problems of gas-phase chemistry. The inclusion of a reaction field representing an implicit solvent in a self-consistent hybrid density functional calculation provides an effective and efficient way to extend this approach to problems of liquid-phase chemistry. In previous work, we have parameterized several models based on this approach, and in the present article, we present several new parameterizations based on implicit solvation models SM5.43 and SM5.43R. In particular, we extend the applicability of these solvation models to several combinations of the MPWX hybrid-density functional with various one-electron basis sets, where MPWX denotes a combination of Barone and Adamos modified version of Perdew and Wangs exchange functional, Perdew and Wangs correlation functional, and a percentage X of exact Hartree–Fock exchange. SM5.43R parameter optimizations are presented for the MPWX/MIDI!, MPWX/MIDI!6D, and MPWX/6-31+G(d,p) combinations with X=0 (i.e., pure density functional theory), 25, 42.8, and 60.6, and for MPWX/6-31G(d) and MPWX/6-31+G(d), with X=0, 42.8, and 60.6; this constitutes a total of 18 new parameter sets. [Note that parameter optimizations using MPW25/6-31G(d) and MPW25/6-31+G(d) were carried out in a previous SM5.43R parameterization.] For each of the five basis sets, we found no significant loss in the accuracy of the model when parameters averaged over the four values of X are used instead of the parameters optimized for a specific value of X. Therefore for each of the five basis sets used here, the SM5.43R and SM5.43 models are defined to have a single parameter set that can be used for any value of X between 0 and 60.6. The new models yield accurate free energies of solvation for a broad range of solutes in both water and organic solvents. On the average, the mean-unsigned errors, as compared with those from experiment, of the free energies of solvation of neutral solutes range from 0.50 to 0.55 kcal/mol and those for ions range from 4.5 to 4.9 kcal/mol. Since the SM5.43R model computes the free energy of solvation as a sum of bulk-electrostatic and non-bulk-electrostatic contributions, it may be used for detailed analysis of the physical effects underlying a calculation of the free energy of solvation. Several calculations illustrating the partitioning of these contributions for a variety of solutes in n-hexadecane, 1-octanol, and water are presented.Acknowledgement This work was supported by a Department of Defense Multidisciplinary University Research Initiative (MURI) grant managed by the Army Research Office and by the National Science Foundation.     

Simple models for nonpolar solvation: Parameterization and testing

Implicit solvent models are important for many biomolecular simulations. The polarity of aqueous solvent is essential and qualitatively captured by continuum electrostatics methods like Generalized Born (GB). However, GB does not account for the solvent‐induced interactions between exposed hydrophobic sidechains or solute‐solvent dispersion interactions. These “nonpolar” effects are often modeled through surface area (SA) energy terms, which lack realism, create mathematical singularities, and have a many‐body character. We have explored an alternate, Lazaridis–Karplus (LK) gaussian energy density for nonpolar effects and a dispersion (DI) energy term proposed earlier, associated with GB electrostatics. We parameterized several combinations of GB, SA, LK, and DI energy terms, to reproduce 62 small molecule solvation free energies, 387 protein stability changes due to point mutations, and the structures of 8 protein loops. With optimized parameters, the models all gave similar results, with GBLK and GBDILK giving no performance loss compared to GBSA, and mean errors of 1.7 kcal/mol for the stability changes and 2 Å deviations for the loop conformations. The optimized GBLK model gave poor results in MD of the Trpcage mini‐protein, but parameters optimized specifically for MD performed well for Trpcage and three other small proteins. Overall, the LK and DI nonpolar terms are valid alternatives to SA treatments for a range of applications. © 2017 Wiley Periodicals, Inc.     

Treating entropy and conformational changes in implicit solvent simulations of small molecules

Implicit solvent models are increasingly popular for estimating aqueous solvation (hydration) free energies in molecular simulations and other applications. In many cases, parameters for these models are derived to reproduce experimental values for small molecule hydration free energies. Often, these hydration free energies are computed for a single solute conformation, neglecting solute conformational changes upon solvation. Here, we incorporate these effects using alchemical free energy methods. We find significant errors when hydration free energies are estimated using only a single solute conformation, even for relatively small, simple, rigid solutes. For example, we find conformational entropy (TDeltaS) changes of up to 2.3 kcal/mol upon hydration. Interestingly, these changes in conformational entropy correlate poorly (R2 = 0.03) with the number of rotatable bonds. The present study illustrates that implicit solvent modeling can be improved by eliminating the approximation that solutes are rigid.     

Generalized conductor-like screening model (GCOSMO) for solvation: An assessment of its accuracy and applicability

We present an assessment on the accuracy of a dielectric continuum solvation model, the generalized conductor-like screening model (GCOSMO), for predicting hydration free energies, tautomeric equilibria, and reaction profiles in solution. © 1996 John Wiley & Sons, Inc.     

GBr(6): a parameterization-free, accurate, analytical generalized born method

The Poisson-Boltzmann (PB) equation is widely used for modeling electrostatic effects and solvation for macromolecules. The generalized Born (GB) model has been developed to mimic PB results at substantial lower computational cost. Here, we report an analytical GB method that reproduces PB results with high accuracy. The analytical approach builds on previous work of Gallicchio and Levy (J. Comput. Chem. 2004, 25, 479), and incorporates an improvement, proposed by Grycuk (J. Chem. Phys. 2003, 119, 4817), of the Coulomb-field approximation used in most GB methods. Tested against PB results, our GB method has an average unsigned relative error of only 0.6% for a representative set of 55 proteins and of 0.4% and 0.3%, respectively, for folded and unfolded conformations of cytochrome b562 sampled in molecular dynamics simulations. The dependencies of the electrostatic solvation free energy on solute and solvent dielectric constants and on salt concentration are fully accounted for in our method.     

Mathematical analysis of the boundary-integral based electrostatics estimation approximation for molecular solvation: exact results for spherical inclusions

We analyze the mathematically rigorous BIBEE (boundary-integral based electrostatics estimation) approximation of the mixed-dielectric continuum model of molecular electrostatics, using the analytically solvable case of a spherical solute containing an arbitrary charge distribution. Our analysis, which builds on Kirkwood's solution using spherical harmonics, clarifies important aspects of the approximation and its relationship to generalized Born models. First, our results suggest a new perspective for analyzing fast electrostatic models: the separation of variables between material properties (the dielectric constants) and geometry (the solute dielectric boundary and charge distribution). Second, we find that the eigenfunctions of the reaction-potential operator are exactly preserved in the BIBEE model for the sphere, which supports the use of this approximation for analyzing charge-charge interactions in molecular binding. Third, a comparison of BIBEE to the recent GBε theory suggests a modified BIBEE model capable of predicting electrostatic solvation free energies to within 4% of a full numerical Poisson calculation. This modified model leads to a projection-framework understanding of BIBEE and suggests opportunities for future improvements.