VOLUME 102, NUMBER 26, JUNE 25, 1998
© Copyright 1998 by the American Chemical Society
LETTERS A Molecular/Continuum Thermodynamic Model of Hydration Henry S. Ashbaugh† Department of Chemical Engineering, Center for Molecular and Engineering Thermodynamics, UniVersity of Delaware, Newark, Delaware 19716
Michael E. Paulaitis* Department of Chemical Engineering, Johns Hopkins UniVersity, Baltimore, Maryland 21218 ReceiVed: March 11, 1998; In Final Form: May 11, 1998
We propose an integrated molecular/continuum thermodynamic model for calculating electrostatic contributions to the free energy of hydration that is based on a correlation function expansion of the free energy. The model uses explicit water simulation results for simple monatomic and diatomic solutes to calculate lower order correlations in this expansion and systematically corrects for higher order correlations using the continuum description of hydration. The model is applied to calculating free energies of charging the tetramethylammonium (TMA) ion and idealized TMA-like solutes in water. We find that the effects of water structure on ion hydration are embodied in the lower order correlations and, thus, are accounted for by the explicit water simulation results. In addition, continuum corrections for higher order correlations are found to be independent of the Born radii chosen for oppositely charged ions. Thus, the model reduces the reliance of the continuum description on parametrization of the Born radii.
Continuum models of hydration are applied routinely and with little difficulty to describe the thermodynamic properties of complex macromolecules in aqueous solution.1-3 The continuum approximation, while accurate for long-range interactions, neglects the molecular details of solvent structure that govern short-range interactions and, as such, neglects the role of solvent structure in determining hydration properties.4,5 Moreover, assuming the solvent can be treated as a macroscopic dielectric on all length scales comes at the price of introducing phenomenological parameters (e.g., Born radii6) whose values have limited physical significance.7 The molecular details of solvent structure, as well as the connection between solvent structure and the macroscopic * To whom correspondence should be addressed. † Present address: Physical Chemistry 1, Center for Chemistry and Chemical Engineering, Lund University, P.O. Box 124, S-221 00 Lund, Sweden.
properties of hydration, can be probed by molecular simulations.8,9 However, the computational effort needed to obtain satisfactory statistical precision places inordinate demands on molecular simulations owing to the large number of degrees of freedom involved for even moderately complex solutes. Given the inherent complexities of evaluating multidimensional molecular potential functions in simulations and the questionable treatment of short-range interactions by the continuum model, alternative descriptions of hydration are needed that can overcome the deficiencies of both approaches. In this Letter, we present a thermodynamic model that selfconsistently combines both explicit water and continuum dielectric descriptions of hydration. One advantage of the model is that the free energy of hydration of complex molecules can be calculated using only a limited number of well-characterized correlation functions obtained from explicit water simulations
S1089-5647(98)01450-3 CCC: $15.00 © 1998 American Chemical Society Published on Web 06/05/1998
5030 J. Phys. Chem. B, Vol. 102, No. 26, 1998
Letters
of the hydration of simple monatomic and diatomic solutes. To test the utility of this hybrid model, we use it to calculate the free energy of charging the tetramethylammonium (TMA) ion and idealized TMA-like solutes in water. These results are compared to free energies calculated directly from molecular simulations. We express the reversible work, or equivalently the free energy, of hydrating a solute composed of n sites in a fixed configuration (r1, ..., rn) as an expansion in terms of one-, two-, three-site, and higher order correlations:10 n
ω(n)(1, ..., n) )
ω(1)(i) + ∑ δω(2)(i, j) + ∑ i)1 i,j
δω(3)(i, j, k) + ... + δω(n)(1, ..., n) ∑ i,j,k
(1)
In this expression ω(1)(i) is the reversible work of solvating a single solute site i, and δω(2)(i, j) is the potential of mean force (PMF) between solute sites i and j in solution:
δω(2)(i, j) ) ω(2)(i, j) - ω(1)(i) - ω(1)(j)
(2)
Three-site and higher order correlations correct for nonadditivity of the PMF. For example, the correction for triplet correlations is obtained by setting n ) 3 and solving for δω(3)(i, j, k) in eq 1:
δω(3)(i, j, k) ) ω(3)(i, j, k) - ω(2)(i, j) - ω(2)(i, k) ω(2)(j, k) + ω(1)(i) + ω(1)(j) + ω(1)(k) (3) Truncating eq 1 at the level of two-site and three-site correlations corresponds to the Kirkwood11 and Fisher-Kopeliovich12 superposition approximations, respectively. Owing to their high dimensionality, three-site and higher order correlations are difficult, if not impossible, to obtain from simulations. To circumvent this problem, we propose an approximation to eq 1 that involves evaluating the lower order correlation functions from molecular simulations and the higher order correlation functions from the continuum description of solvation: n
ω(1)(i) + ... + ∑ i)1
ω(n)(1, ..., n) ≈ ω(n|m)(1, ..., n) ) [
n
δω(m)]sim + [ω(n)(1, ..., n) - ∑ ω(1)(i) - ... ∑ m-tuples i)1
∑
δω(m)]cont (4)
m-tuples
The first term in brackets on the right-hand side of eq 4 is the m-body superposition approximation evaluated from simulation. The second term is a continuum correction to the superposition approximation for higher order correlations. We refer to ω(n|m)(1, ..., n) as the m-body solvation free energy. A useful feature of eq 4 is that the balance between simulation and the continuum description can be adjusted unambiguously by changing m. For example, ω(n|0)(1, ..., n) is the solvation free energy evaluated from the continuum description alone, while ω(n|n)(1, ..., n) is the solvation free energy obtained from molecular simulation. Canonical ensemble Monte Carlo (MC) simulations of 216 water molecules and 1 solute at 25 °C and a density of 0.997 g/cm3 were performed. Water was modeled using the simple point charge potential.13 TMA Lennard-Jones (LJ) interactions
TABLE 1: Free Energies of Charging Different Me Solutes in Water (q Is the Unit of Charge on Each Me Group and Is Equal to +0.25e)a continuum description solute (charges)
simulation
single Me groups (+q) -1.92 (0.03) (-q) -6.88 (0.03) Me pairs (+q, +q) -9.77 (0.05) (-q, -q) -20.55 (0.09) (+q, -q) -1.98 (0.03) c-TMA (+q, +q, +q, +q) -41.02 (0.09) a-TMA (-q, -q, -q, -q) -66.35 (0.13) p-TMA (+q, +q, -q, -q) -3.79 (0.03)
R+ ) 2.2 Å and R+ ) R- ) 2.0 Å R- ) 1.8 Å -5.14 -5.14
-4.66 -5.71
-16.77 -16.77 -2.95
-15.50 -18.25 -3.01
-56.99
-53.19
-56.99
-61.35
-5.38
-5.48
a
Free energies are reported in units of kcal/mol. The numbers in parentheses denote 1 standard deviation.
were modeled using the potential of Jorgensen and Gao.14 Water-TMA cross LJ interaction parameters were obtained using geometric mean combining rules. Long-range electrostatic interactions were evaluated using the Ewald summation technique with conducting boundary conditions.15,16 Appropriate Ewald self-interaction corrections were made to the charging free energies.4,17,18 The free energy was calculated using coupling parameter integration.19,20 The TMA molecule was modeled as a tetrahedral arrangement of four methyl (Me) groups bonded to a central nitrogen atom. The buried nitrogen atom was neglected in these calculations. The TMA ion carries a total charge of +1e, which was distributed equally among the Me groups. We also considered a TMA anion with charges of -0.25 e on each Me group and a polar TMA-like solute with charges of +0.25e on two Me groups and -0.25e on the other two Me groups. Hereafter, we denote the cation, anion, and polar solute as c-TMA, a-TMA, and p-TMA, respectively. Application of the one-body and two-body hydration models required additional simulations to calculate the free energy of charging a single Me (( 0.25e) and a pair of Me groups separated by 2.47 Å (the intramolecular Me-Me separation in TMA). All combinations of (0.25e charges on the two Me groups were considered for the Me-Me pairs (see Table 1). Complementary continuum electrostatic calculations of the charging free energy were also performed for each solute in water.21 In this treatment, the solvent contribution to the free energy of charging was obtained by numerical solution of Poisson’s equation1 for a molecular solute immersed in an aqueous dielectric medium. To test the sensitivity of the model predictions to the continuum model parametrization, we considered two cases: (i) equivalent Born radii for oppositely charged Me groups (R+ ) R- ) 2.0 Å) and (ii) different Born radii for oppositely charged Me groups (R+ ) 2.2 Å and R- ) 1.8 Å). The values of R+ and R- were selected as typical Born radii for simple ions in water. The boundary between the solute interior (dielectric constant of 1) and water (dielectric constant of 80) was defined by the molecular surface,25,27 which was calculated by rolling a 1.4 Å radius probe over the solute surface defined by the Born radii. The calculated free energies of charging the different Me solutes in water are given in Table 1. An estimate for the free energy of charging the TMA ion can be obtained from experiment as the difference between the measured free energy of transferring the solute from vacuum to water, -38.2 kcal/
Letters
J. Phys. Chem. B, Vol. 102, No. 26, 1998 5031
TABLE 2: Calculated Free Energies of Charging TMA-like Solutes in Water Compared to the Exact Simulation Results in Table 1a solute c-TMA a-TMA p-TMA
ac bd a b a b
ω(4|1)(1,2,3,4)
∆b
ω(4|2)(1,2,3,4)
∆
KSA
∆
-44.11 (0.11) -42.17 (0.11) -63.97 (0.10) -66.04 (0.10) -2.42 (0.07) -2.32 (0.07)
-3.08 (0.14) -1.15 (0.14) 2.38 (0.16) 0.31 (0.16) 1.36 (0.08) 1.47 (0.08)
-40.82 (0.35) -40.87 (0.35) -65.72 (0.57) -65.73 (0.57) -4.22 (0.21) -4.26 (0.21)
0.21 (0.37) 0.15 (0.37) 0.63 (0.58) 0.62 (0.58) -0.43 (0.21) -0.47 (0.21)
-43.34 (0.35)
-2.32 (0.37)
-68.25 (0.57)
-1.90 (0.58)
-3.07 (0.21)
0.72 (0.21)
a Free energies are reported in units of kcal/mole. The numbers in parentheses denote one standard deviation. b ∆ is the calculated charging free energy minus the value from simulation given in Table 1. c Continuum description results for R+ ) R- ) 2.0 Å. d Continuum description results for R+ ) 2.2 Å and R- ) 1.8 Å.
mol,28 and the measured free energy of transferring neopentane from vacuum to water, +2.7 kcal/mol.29 This value is -40.9 kcal/mol, which is in excellent agreement with the free energy of charging c-TMA obtained from simulation: -41.02 ( 0.09 kcal/mol. Limitations of the continuum description are clearly seen in Table 1 for the charging free energies of the oppositely charged solutes calculated using equivalent Born radii. In each case, the two free energies are equivalent. In contrast, both experiment30 and simulations4,18 show that ion hydration is asymmetric; i.e., anion hydration is favored over cation hydration for ions of similar size. This asymmetry is attributed to differences in water structure in the first hydration shell around positive and negative ions.4,18,31-33 In the continuum description these solvent specific structural effects are mimicked by assigning cations (anions) Born radii that are greater (less) than the van der Waals radius of the solute. Indeed, as shown in Table 1, the continuum description with R+ ) 2.2 Å and R- ) 1.8 Å captures this asymmetry, although the agreement with simulation is only qualitative. Finally, we note that the free energies of charging the two polar solutes in Table 1 are much less sensitive to the changes in Born radii than are the ionic solutes. Free energies of charging the TMA-like solutes in water calculated using the one- and two-body hydration models are reported in Table 2. Also included in this table are charging free energies calculated from simulation alone using the Kirkwood superposition approximation (KSA) to eq 1. Charging free energies obtained from the one-body hydration model with R+ ) R- are comparable to those obtained using KSA. Both predictions are in reasonably good agreement with the exact simulation results, although somewhat larger relative differences are evident for c-TMA and p-TMA compared to a-TMA. The results for ω(4|1)(1,2,3,4) can be improved significantly if different Born radii are used for positively and negatively charged Me groups in the continuum correction. For c-TMA, the difference between ω(4|1)(1,2,3,4) and the exact simulation results decreases by more than a factor of 2. For a-TMA, this difference decreases by nearly an order of magnitude. The results for p-TMA are, however, slightly worse. Overall, ω(4|1)(1,2,3,4) correctly predicts the more favorable hydration of a-TMA relative to c-TMA even when the continuum correction uses equivalent Born radii. We conclude, therefore, that the observed asymmetry in the hydration of positive and negative TMA ions is largely accounted for by differences in water structure around the individual charged Me groups alone. That is, the asymmetry in TMA ion hydration reflects in large part water organization in response to local interactions with the solute. The two-body hydration model significantly improves the agreement between the predicted charging free energies and the exact simulation results. For c-TMA and a-TMA, the differ-
ences are now comparable to the statistical uncertainties in the simulation results. For the polar solute, the error is within 2 standard deviations of the statistical uncertainty. It is apparent that adding the continuum correction for higher order correlations to KSA systematically improves the predicted charging free energies in all cases. Moreover, ω(4|2)(1,2,3,4) is insensitive to the Born radii chosen to calculate this correction. We conclude, therefore, that only the lower order continuum correlations are sensitive to the Born radii, which is consistent with the notion that solvent specific structural effects embodied in the Born radii are largely accounted for in the hydration of single solute sites. Although the simulations needed to calculate ω(4|3)(1,2,3,4) have not been performed, the continuum corrections to the Fisher-Kopeliovich superposition approximation to eq 1 were found to be negligible within the calculated statistical uncertainties of the two-body model. Two significant conclusions can be drawn from our results: First, the effects of water structure on ion hydration are accounted for by the leading (lower order) terms of eq 1, which are calculated from exact simulation results without the need for different Born radii for positively and negatively charged Me groups. Second, the continuum correction to the KSA makes a significant, systematic improvement to the free energies of charging the TMA ion and TMA-like solutes in aqueous solution. An important consequence of these conclusions is that reliance of continuum corrections on parametrization of the Born radii has been reduced. In addition, only the lower order correlation functions need to be evaluated from explicit water simulations. It should be noted that our model has been applied to solutes that are roughly spherically symmetric and have low surface charge densities. Additional work must be performed on solutes having a variety of shapes and conformational degrees of freedom to generalize the results presented here. It would also be interesting to calculate the PMF between two multisite solutes in aqueous solution. Our model could be easily extended to do so. In previous work, the KSA was successfully applied to predict DNA conformational equilibria as a function of salt concentration using the restricted primitive model of electrolytes.34 The model presented here would improve upon these calculations by including specific water structural effects, higher order correlations in the free energy expression, and no need for artificial parametrization of the restricted primitive model. Indeed, it has been shown that incorporating water structure through the PMF between ions in aqueous solution leads to the prediction of an attractive force between charged surfaces at molecular separations, in agreement with experimental observations.5 Finally, we note that applications of our model are not limited to the approach we have described herein, that is, combining explicit water simulations with the continuum description of
5032 J. Phys. Chem. B, Vol. 102, No. 26, 1998 hydration. Alternatively, our explicit water simulation results could be combined with higher order correlation functions obtained, for example, from integral equation theory. In general, the model presented here could be used to incorporate any approximate treatment of hydration with a more rigorous, molecular description to yield improved predictions of macromolecular hydration at a lower computational expense. Acknowledgment. We are indebted to Shekhar Garde, Gerhard Hummer, and Lawrence Pratt for numerous fruitful discussions. Financial support from the National Aeronautics and Space Administration (NAG3-1954), the National Science Foundation (Grants BES-9210401 and BES-9510420), and a National Science Foundation Fellowship for H.S.A. (GER9253850) is gratefully acknowledged. References and Notes (1) Jackson, J. D. Classical Electrodynamics, 2nd ed.; Wiley: New York, 1975. (2) Honig, B.; Nicholls, A. Science 1995, 268, 1144. (3) Nakamura, H. Q. ReV. Biophys. 1996, 29, 1. (4) Garde, S.; Hummer, G.; Paulaitis, M. E. J. Chem. Phys. 1998, 108, 1552. (5) Marcelja, S. Nature 1997, 385, 689. (6) Born, M. Z. Phys. 1920, 1, 45. (7) Pratt, L. R.; Hummer, G.; Garcı´a, A. E. Biophys. Chem. 1994, 51, 147. (8) Garde, S.; Hummer, G.; Garcı´a, A. E.; Pratt, L. R.; Paulaitis, M. E. Phys. ReV. 1995, E53, R4310. (9) Ashbaugh, H. S.; Paulaitis, M. E. J. Phys. Chem. 1996, 100, 1900. (10) Green, H. S. The Molecular Theory of Fluids; North-Holland Publishing Company: Amsterdam, 1952. (11) Kirkwood, J. G. J. Chem. Phys. 1935, 3, 300. (12) Fisher, I. Z.; Kopeliovich, B. L. Dokl. Akad. Nauk SSSR 1960, 133, 81. (13) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces: Proceedings of the Fourteenth Jerusalem Symposium on Quantum Chemistry and Biochemistry; Pullman,
Letters B., Ed.; D. Reidel Publishing Company: Dordrecht, Holland, 1981; pp 33142. (14) Jorgensen, W. L.; Gao, J. J. Phys. Chem. 1986, 90, 2174. (15) Ewald, P. P. Ann. Phys. 1921, 64, 253. (16) De Leeuw, S. W.; Perram, J. W.; Smith, E. R.; Proc. R. Soc. London 1980, A373, 27. (17) Hummer, G.; Pratt, L. R.; Garcı´a, A. E. J. Phys. Chem. 1995, 99, 14188. (18) Hummer, G.; Pratt, L. R.; Garcı´a, A. E. J. Phys. Chem. 1996, 100, 1206. (19) McQuarrie, D. A. Statistical Mechanics; Harper Collins: New York, 1976. (20) Simulations were performed for λ ) 0 (solute charges off), 1/2, and 1 (solute charges on) to evaluate ∂ωi/∂λ, the derivative of the free energy with respect to λ. The free energy of charging each solute was then calculated by numerically integrating ∂ωi/∂λ using Simpson’s rule. Each simulation was performed for a total of 8 × 105 MC cycles (1 MC cycle corresponds to one attempted move for each water molecule and four attempted moves for the solute). Statistical uncertainties in the free energies were estimated from block averages by dividing each simulation into five equal sized blocks that were assumed to be independent. (21) Poisson’s equation1 was solved numerically using the boundary element method (BEM).22-24 Surface discretization of the solute/solvent boundary was performed using Zauhar’s SMART (Smooth MoleculAR Triangulator) program.25 An angle parameter of 25° was used to triangulate the molecular surface. The SMART discretized surface was input into the BEM Poisson solver developed by Neal and Lenhoff.26 (22) Zauhar, R. J.; Morgan, R. S. J. Mol. Biol. 1985, 186, 815. (23) Rashin, A. A. J. Phys. Chem 1990, 94, 1725. (24) Yoon, B. J.; Lenhoff, A. M. J. Comput. Chem. 1990, 11, 1080. (25) Zauhar, R. J. J. Comput.-Aided Mol. Des. 1995, 9, 149. (26) Neal, B. L. Ph.D. Thesis, University of Delaware, 1997. (27) Richards, F. M. Annu. ReV. Biophys. Bioeng. 1977, 6, 151. (28) Marcus, Y. Biophys. Chem. 1994, 51, 111. (29) Ben-Naim, A.; Marcus, Y. J. Chem. Phys. 1984, 81, 2016. (30) Latimer, W. M.; Pitzer, K. S.; Slansky, C. M. J. Chem. Phys. 1939, 7, 108. (31) Hummer, G.; Pratt, L. R.; Garcı´a, A. E.; Berne, B. J.; Rick, S. W. J. Phys. Chem. B 1997, 101, 3017. (32) Ashbaugh, H. S.; Wood, R. H. J. Chem. Phys. 1997, 106, 8135. (33) Ashbaugh, H. S.; Kaler, E. W.; Paulaitis, M. E. Biophys. J., in press. (34) Soumpasis, D.-M. Proc. Natl. Acad. Sci. U.S.A. 1984, 81, 5116.
