Absolute pKa calculations with continuum dielectric methods - The

Jul 1, 1991 - Approach for Predicting the Standard Free Energy Solvation of H and Acidity Constant in Nonaqueous Organic Solvents. Hossein ... Role of...
0 downloads 7 Views 1MB Size
Download PDF
J. Phys. Chem. 1991,95, 5610-5620

5610

In equation IIA, the permittivity of the solvent has been calculatedIZ at t = -65 OC from the equation12 e = 1.61 + -0.2 X 10-3(-65) = 1.74. The arbitrary parameters r = 5.5 and r = 6.0 A as the electron-solvent distances have been used. The model envisages a tetrahedral trap formed by four hydrocarbon molecules surrounding the localized electron at the center of the cm, m = 0.64 X trap. For r = 5.5 X esu cm. With this value we can now calculate the polarization energy Up due to the electron at the center of the trap as em Up = -4= -2.33 X lo-" erg/molecule = -0.146 eV

This corresponds to increasin r by a factor of (lS)'/' in order to match, after conversion to I the values of E,. In other words, the same calculated values obtained with r = 5.5 X 10% cm and r = 6.0 X lo" cm, have been obtained (for an octahedral model) with r = 6.09 X cm and r = 6.64 X 10-8 cm, respectively. In the Introduction of the paper it has been pointed out that E,, is proportional to the molar polarization of the solvent. We advance the hypothesis that lOd E,,. This implies that the barrier of energy of activation for mobility of excess electrons, subjected to a gradient of potential of 1 V/cm, is comparable to the work necessary to extract the electron from the potential well created by the polarization or the solvent. Table V reports a breakdown of the above calculation for various straight-chain hydrocarbons. Figure 4 reports again E,, vs the molar polarization of the solvent. The calculated molar polarization energy lOpl is also reported for comparison, showing a surprisingly good correlation to exist, given the simplicity of the model assumed. The correlation between E,, and 10plis, however, lost if one tries the above calculation for non-linear or branched hydrocarbons. A discussion in the literatureI3 exists involving anisotropy in the polarizability of branched hydrocarbons and possible directional effects in electron localization.

%

f

69

which on a molar scale becomes (2.33 X 10-20)(6.02X loz3) = 3.35 kcal/md lop1 = 4.184 Repeating the calculation with r = 6.0 X 10" cm, 10pl= 2.36 kcal/mol, whereas the reported activation energy for electrical mobilityz E,, = 2.6 kcal/mol for excess electrons in propane. Notice that altering the arbitrary model of a tetrahedral trap to, for example an octahedral trap, changes Upto Up = -6em/e#. (12) Data from Handbook of Chemistry and Physics, 44th 4.; The Chemical Rubber Co.: Cleveland, OH, 1963.

(13) Dodelet, J.-P.; Freeman, G. R. Can. J . Chem. 1972, 50, 2667.

Absolute pK, Cakulattons with Continuum Dleiectric Methods Carmay Lim, Don Bashford, and Martin Karplus* Department of Chemistry, Haruard University, 12 Oxford Street, Cambridge, Massachusetts 02138 (Received: July 16, 1990) Solvation free energies and pK, values of models for ionizable side chains of amino acids are calculated by using continuum dielectric methods, integral equation techniques are also investigated. The dependence of the solvation free energies on the parameters (charges and van der Waals interactions used to describe the model compounds) is explored by comparing different sets that are being used in protein and liquid simulations. The solvation free energies calculated with both continuum and integral equation methods and various parameter sets agree qualitatively with experiment but are not accurate enough to yield absolute pK, values. To obtain the experimental solvation free energies and pK, values of the model compounds with the continuum dielectric method, an adjusted parameter set is introduced; only very small changes from the standard parameter values are required. The set of calibrated parameters is tested on some bifunctional compounds and yields pK, changes in reasonable agreement with experiment. However, the pK, changes are very sensitive to the solution conformation. This may result in large pK, errors if conformational changes (e.g., between the ionized and neutral species) are not taken into account.

1. Introduction Knowledge of the pK, values of ionizable groups is important for an understanding of many a r m of chemistry, both in the gas phase and solution. They are of particular interest for elucidating reaction mechanisms, especially those involving proton transfers, and for interpreting the binding of substrates or inhibitors to enzymes. However, experimental determinations of individual pK, values are difficult in complex systems. One case in point concerns the direct measurement of titrating group of catalytically important residues or of substrates in enzyme-substrate complexes. Kinetic assignments of pK,'s are often complicated by uncertainties in interpreting the pH dependence of the measured parameters.' It is useful, therefore, to have reliable and accurate means of calculating relative and/or absolute pK, values and to have an understanding of the factors involved. Such an understanding is essential for interpreting the measured effective pK, values in proteins and other complex polyions. Both microscopic and macroscopic theoretical methods are now available for the estimation of solvation free energie~.~JThis ( I ) Knowlea, J.

SCHEME I BH(g)

A%

~AGw)

BH(s)

We)

AGa

+

JAGW*)

AW)

B(s)

H*(g)

+ H*(s)

makes it possible, in principle, to determine theoretical relative or absolute pK, values from the thermodynamic cycle (Scheme I). In Scheme I, AG, and AGa are the gas-phase and solution free energy of ionization and the AG,'s are solvation free energies. Given these results, the pK, of BH(s) is

1

pK, = 2.303RTAGa

=- 1 {AG8 + AG,(B) - AG,(BH) 2.303R T

+ AG,(H+))

(2) Jorgensen, W. L. Acc. Chem. Res. 1989, 22, 184. (3) Gilwn, M. K.; Honig, B. H. Proteins 1988, I , 7 .

R. CRC Crfr. Reu. Blochem. 1975, 3, 165. 0022-3654/91/2095-5610$02.50/0

-1

Q

1991 American Chemical Society

(1.1)

Absolute pK, Calculations

The Journal of Physical Chemistry, Vol. 95, No.14, 1991 5611

and the difference in pK,’s of BH and AH is [AGAB) - AGs(A)I - [AGABH) - AGAAH)IJ (1.2) The determination of relative pK, leads to the important simplification that it does not require a knowledge of the proton solvation free energy, a quantity concerning which there exists considerable experimental uncertainty (see below) and for which there are no reliable calculations. In determining the relative or absolute pK, values, the AG, values can be obtained from experiments or ab initio methods, and the AG, values can be calculated by using Monte Carlo simulations, simulation methods combined with a reaction field approach, integral equation techniques, and continuum dielectric methods. A number of workers have employed one or another of these approaches in recent years with some success. Jorgensen et al.‘.’ used Monte Carlo simulations to obtain apriori predictions for the relative pK,’s of organic compounds in water. Ab initio quantum mechanics (Moller-Plesset theory to third order and the 6-31+G* basis set) was used to compute the gas basicities. The computed gas basicities differ from the experimental values by at most 3.3 kcal/mol or 0.9%. Monte Carlo simulations were performed to calculate the difference in solvation free energies using statistical mechanical perturbation theoryS5 The relatioe solvation free energies for the neutral molecules are within 1.3 kcal/mol of their experimental values. However, this can easily lead to a 50% error in the solvation free energies since these are typically a few kcal/mol for neutral molecules; e.g., the calculated AG,(MeSH) - AG,(MeCN), -1.3 f 0.2 kcal/mol, is half the experimental value, -2.6 k ~ a l / m o l . ~ The solvation free energies for ionic species cannot be compared to experiment since there are no known measurements of solvation free energies for the ions; the values usually referred to are estimated from eq 1.1. The calculated pK, of acetonitrile, methylamine, and ethane are consistent with experimental estimates, but the experimental range of pK, values is broad, especially in the case of ethane (i-e., the experimental pK, estimates range from 42 to 60). Jorgensen et ah5 attributed the main error in the calculations to the lack of explicit treatment of polarization effects in the anion-water interactions. Guissani et aL6 used molecular dynamics and test particle methods to calculate the excess chemical potentials of H20, H,O+, and OH- in water and used the results to estimate the absolute pK, of water at three different temperatures. Their calculated pH values were of the correct order of magnitude, the largest error occurring at the highest temperature of study, 593 K,where the calculated pK, (3.8) differs by 2.1 pK, units from the experimental value (5.9). They found that the polarization of water molecules by the ion and the back polarization of strong hydrogen bonds around the ion by water molecules were very important in the absolute pK, calculations. Rullmann and van Duijnen7 have implemented a polarizable water model in Monte Carlo simulations to calculate solvation energies and relative protonation energies of several amines in solution. The vacuum proton affinities were calculated at the Hart-Fock level with Rm-Siegbahn basis set with and without polarization functions. The calculated difference in proton affinities between ethylamine and propylamine agrees within 2 kcal/mol with the experimental value, although the individual values differ by 13-14 kcal/mol; this is due presumably to the neglect of correlation effects. Solvation energies were calculated with the Monte Carlo algorithm by using spherical boundary conditions and a continuum reaction field. The solute and two or three layers of solvent were embedded in a dielectric continuum; the oxygen centers of the water molecules were kept a minimum distance of 1.58 A from the dielectric boundary to prevent di(4) Jorgenren, W. L.; Briggs, J. M.;Gao, J. J. Am. Chem.Soc. 1987,109, 6857. (5) Jorgenren, W. L.; Briggs, J. M. J. Am. Chem. Soc. 1987,111,4190. (6) Guiwani, Y.; Guillot, B.; Brat=, S.J . Chem. Phys. 1988,88, 5850. (7) Rullmann, J. A. C.; van Duijnen, P. Th. Mol. Phys. 1988, 63, 451.

vergence of the potential from the reaction field. This reaction field plus exclusion model in conjunction with a self-consistent model for electronic polarization of water gave solvation energies of ethylamine, propylamine, and ethylenediaminewhich arc within 2.2-5.3 kcal/mol (or 12-40%) of the experimental values. The magnitudes of the solution protonation energies of the amines relative to ethylene diamine exceed the experimental values by a factor of 2-14 and the protonation energies are not in the correct order. Russell and Warshe18 developed a protein dipoles Langevin dipoles (PDLD) model which includes charge-charge and charge-dipole interactions in the protein; the solvent molecules are modeled as Langevin dipoles and the average solvent polarization rather than the average of the actual polarization is calculated. The PDLD model was employed to calculate the intrinsic pK, values of the acidic groups of bovine pancreatic trypsin inhibitor (BPTI). The PDLD results are treated and calibrated asfree energies rather than enthalpies; the two Langevin dipole parameters, and C’, are calibrated to reproduce the solvation free energies of the ions in water. The calculated ionization energies of the four acidic groups in BPTI, when all the other ionizable groups are neutral or when each acid‘s nearestneighbor acid group is ionized, differ from the corresponding experimental estimates8 by up to 7 kcal/mol. As a test of the sensitivity of the results to the parametrized charges, they repeated the calculations with a charge set scaled by 20%. The energies obtained with the two sets of parameters are qualitatively similar. Alternatives to such computationally intensive simulation methods for calculating solvation free energies are provided by the extended RISM (XRISM) integral equation”z and macroscopic methods based on a continuum description of the solFrom an analysis of an asymptotic solution to the extended RISM integral equation in the limit of a large spherical solute, Roux et a1.16 showed that the Born model, which treats the solvent as a structurelesscontinuum, was a good approximation to microscopic models. Gilson and Honig3 have obtained finite difference solutions to the Poisson-Boltzmann equation to calculate solvation free energies of methanol, methylammonium, ammonium, and acetate. Their calculated values (using AG,(H+) = -262.5 kcal/mol) differ from the experimental values by 8-19%. The main advantage of using a continuum description of the solvent is that it speeds up computation considerably. Rashin and Namboodiri17 used a reaction field approach to calculate the hydration enthalpies of water, methanol, derivatives of ammonium, acetate, and nitrate. Using AH8(H+)= -262.2 kcal/mol to obtain absolute hydration enthalpies for the ions, the calculated hydration enthalpies, which include both electrostatic and nonpolar contributions, of all the molecules except acetate agree within 6% of the experimental values. The calculated hydration enthalpies of water and methyl ammonium are within 2% of the Monte Carlo simulation results.18 However, as noted by Roux et a1.,I6 the Born model is not appropriate for the calculation of solvation enthalpies since the latter are sensitive to the temperature dependence of the Born radius, which was not included in the work of Rashin and Na~nboodiri.’~ From the calculations described above, it is clear that, although the qualitative results obtained from the various methods are satisfactory, the quantitative errors are such that accurate ab initio pK, values are still very difficult to determine. To indicate what (8) Ruseell, S. T.; Warshel, A. J. Mol. Elol. 1985, 185, 389. (9) Chandler, D. In The Liquid State of Matter: Fluids, Simple and Complex; North-Holland New York, 1982; Vol. 8. (10) Hirata, F.; R a k y , P. J. Chem. Phys. Len. 1981, 83, 329. (1 1) Chandler, D.; McCoy,J. D.; Singer, S. J. J. Chcm. Phys. 1986,85, 5971. (12) Hirata, F.; Rmky, P. J.; Pcttitt, B. M.J. Chem. Phys. 1983,78,4133. (13) Warwicker, J.; Watson, H. C. J . Mol. Elol. 1982, 157, 671. (14) Sternbcrg, M.J. E.; Hays, F. R. F.; Russcll, A. J.; Thomas, P. 0.; Fersht, A. R. Nature 1987,330,86. (15) Bashford, D.; Karplus, M.J . Mol. Bfof. 1988, 203, 507. (16) Roux, B.; Yu, H. A.; Karplus, M. J . Phys. Chrm., in press. (17) Rashin, A. A.; Nambocdiri, K. J. Phys. Chem. 1987, 91, 6003. (IS) Jorgensen, W. L.; Gao, J. J . Phys. Chem. 1986, 90,2174.

5612 The Journal of Physical Chemistry, Vol. 95, No. 14, 199’1 is involved we refer to Scheme I and eq 1.1. The gas-phase basicities are in the range 2owoO kcal/mol, and the solvation free energies are in the range do to -150 kcal/mol for ions, and 0 to -10 kcal/mol for neutral molecules; the proton solvation free energy of -254 to -261 kcal/mol (depending on the method, see below) has already been mentioned. Thus, it is clear that the calculation of absolute pK, values in solution, which are in the range 4-16 pK, units (as compared to gas phase values of 150-273 pK, units), requirts extremely accurate gas basicities and solvation free energies. An error of 1-2 kcal/mol results in pK, shifts of 1.5-3 pK, units. For gas-phase basicities, accurate values can be achieved, in principle, by use of high-level ab initio methods for small model and from experimental data for more complex systems, though experimental gas basicities are accurate to only 2-5 kcal/mol. For the solvation free energies, the results obtained from all the methods discussed depend in an essential manner on the choice of empirical parameters (primarily charges and van der Waals radii) that are used for the system of interest. The inherent uncertainty in the parameters introduces a significant uncertainty in the calculated pK, values. To be able to assess the significance of pK, calculations, it is of considerable interest, therefore, to evaluate the parameter dependence of the solvation estimates. In this paper, we investigate various aspects of the calculation of absolute solution pK, with the thermodynamic cycle in Scheme I. The systems chosen for study include methanol, methylammonium, imidazolium, and acetic acid since they are good model systems for the ionizable side chains of serine, lysine, histidine, and aspartic acid, respectively. Experimental gas basicities were used. Solvation free energies of the protonated (BH) and unpmtonated (B) molecule for zero ionic strength were calculated by using integral equation techniques,e2o and continuum dielectric method^^*^*-'^ with ab initio geometries. The solvation free energy of the proton was taken from experiment. We assume that AG,(H+) is equal to -259.5 kcal/mol or 4.5 V which is an average of five separate measurements of the standard hydrogen potential (-4.43, -4.43, -4.44, -4.48, and -4.73 V)? the latter constitutes a 7 kcal/mol spread in AG,(H+) values ranging from -254 to -261 kcal/mol. The precision to which AG,(H+) can be known is limited by the fact that the standard hydrogen potential cannot be obtained by measurement alone, and some independent quantity is needed to determine an absolute half-cell potentiaLZ1 With AG,(H+), the gas basicities and solvation free energies, we can obtain absolute and relative pK, from eq 1.1 and eq 1.2, respectively. Both integral equation techniques and continuum dielectric methods were used to explore the parameter dependence in calculating solvation free energies. The input charges and radii are derived from the CHARMM22and 0PLSz3parameter sets that are used in liquid simulations. Since we are ultimately interested in determining relative pK, changes in more complex systems (e&, in enzymesubstrate interactions), we need to choose an optimized set of atomic charges and radii which yields the absolute solution pK, of ionizable groups of interest in aqueous solution. The continuum dielectric method was chosen to calibrate the parameters for two reasons. First, it is relatively fast, yet produces results of quality similar to the more compute intensive Monte Carlo free energy calculations. It is also possible to obtain absolute rather than relative values in a simple fashion. Integral equation techniques were found to yield solvation free energies of poorer quality than the continuum dielectric method (see below). Second, the dielectric method can be applied easily to nonhomogeneous many particle systems, such as protein^?^^.^^ while integral equations (19) Hehrc, W. J.; Radom,.L.; Schleycr, P. v. R.; Pople, J. A. Ab. Initio Moleculur Orbftal Theory; Wilcy: New York, 1986. (20) Yu, H. A.; Karplua, M. J . Chem. Phys. 1988,89,2366. (21) Reiu, H.; Hcller, A. J . Phys. Chem. 1985,89,4207. (22) Brooks, B. R.; Bruccoleri, R. D.; Olafson, B. 0.;States, D. J.; Swaminathan, S.; Karplur, M. J. Compur. Chem. 1983, I , 187. (23) Tirado-Rives, J.; Jorge~cn,W. L. J . Am. Chem. Soc. 1988, 110, 1657. (24) Barhford, D.; Karplw, M.Elochemlsrry 1990,29, 10219. (25) Gilron, M.K.; Honig, 8. H. Elopolymers 1986, 25, 2097.

Lim et al. are limited to molecules not significantly larger than a 25-atom peptide in aqueous solution. Even though the standard CHARMM and OPLS charges and radii are chosen to reproduce the interaction energies of model compounds with water molecules, this is not sufficient to obtain quantitative agreement of the calculated solvation free energies with experiment. Moreover, the van der Waals radius does not necessarily determine the effective Born radius in the continuum mode1.16*z6We exploit the sensitivity of the AG, values to the atomic charges and radii and tune them to reproduce the experimental pK, values of model compounds. Since the pK, values change by less than 0.2 pK, units in going from zero to 0.2 M ionic strength,2’ the set of adjusted charges and radii is appropriate for titration calculations in that range. This makes the results appropriate for use with proteins since the pK, values are generally measured at an ionic strength of 0.1 or 0.2 M. We have tested the set of adjusted parameters by determining the relative pK, values for several small molecules with bifunctional groups; viz., aminoethanol, 4hydroxymethyMdazole. and 3-hydroxypmpanoic acid. These results gave a consistent and accurate model that should be applicable to pK, shifts in proteins. Section 2 gives a brief discussion of the methods; for more details, the reader is referred to the original reference listed below. Section 3 presents the results and an analysis of their significance. In section 4 a summary of the conclusion is given and future work is outlined. The integral equation results are described in an Appendix. 2. Method 2.1. Calculatio~.The ab initio calculations were performed . ~geometries were with the programs GAUSSIAN82 and ~ 6 All fully optimized at the HartreeFock level by using the split valence 6-31GS basis set. The difference between the total energies of the model compounds with and without the ionizing proton yields the proton affinity. The ab initio geometries are employed in calculating the solvation free energies. Also the results obtained with Mulliken charges from the ab initio calculations are compared with those from the CHARMM and OPLS charges. Although experimental AG,, are used in calculating solvation free energies of the model systems, ab initio proton affinites are determined for comparison. 2.2. Continuum Dielectric Metbod. The continuum dielectric method is a generalization of the Born model for solvation free energies. The solute molecules are treated as irregularly shaped objects whose interior has an uniform dielectric constant el, with point charges q, at positions corresponding to the atomic nuclei, ri. The electrostatic contribution to AG, is the difference between the free energies of charging these objects in vacuo and in solution. The charges are taken to be the same in vacuo and in solution. Any tendency for solvation to alter the charge distribution by electronic polarization effects, for example, is presumed to be included in ci. The electrostatic field, 4, is determined by the Poisson equation with a dielectric constant that changes from el to the solvent dielectric constant, eo, at the solute/solvent interface; Le.,

Vt(7) V4(7) = -4*p(r3

(2.1)

For the present application, it is appropriate to set the ionic strength of the solution equal to zero; for nonzero ionic strengths the Poisson equation is replaced by the Poisson-Boltzmann equation.29 In the interior, the solution to the Poisson equation can be decomposed into Coulombic terms and a term due to the dielectric interface, (26) Jayaram, B.;Fine, R.; Sharp,K.; Honig, B.J. Phys. Chem. 1989,93, 4320. (27) Smith, R. M.; Martell, A. E. Crftfcol Srabfliry Conrronrs: Second Suppkment; Plenum Prw: New York, 1989; Vol. 6.

(28),Binkley, J. S.; Frisch, M. J.; Defrccr. D. J.; Raghavachari, K.; Whiteside, R. A.; Schlcgel, B. H.; nuder, E. M.;Pople, J. A. Gousslon 82; Carnegie Mellon: Pittsburgh, PA 15213. (29) Mquarrie. D. A. Srorfsrlcol Mechonfcs; Harper and Row: New York, 1976.

Absolute pK, Calculations

The Journal of Physical Chemistry, Vol. 95, No. 14, 1991 5613

Although an integral formulation will not be used in the actual calculations, it is convenient to think of Q* as a potential arising from the induced surface charge, u, at the dielectric boundary,

Use of eq 2.3 provides an alternative approach to the solution of the praent problems by finite element methodsM The work of charging the system is (2.4)

where Q is now considered to be a function of the set of charges q = (q,, g2, ...) and can be separated into a self term for qj, j = i, and a cross term for ql, j # i, of the form

The Coulomb terms in eq 2.2 contributecharge-charge interaction terms and infinite point-charge self-energy terms to W. These terms are the same in vacuo and in solution and do not contribute to the solvation free energy. Dropping the Coulomb terms and using the fact that 4* is linear in the charges, one finds

where the subscripts s and g indicate potentials calculated in solvent and in the gas phase, respectively. As an example, the model is used to derive the Born formula for spherical monatomic ions. The ion is treated as a sphere of radius R with a charge Q at the center. In a solvent with zero ionic strength and dielectric constant the solution of the Poisson equation in the interior is

The second term corresponds to the 4* of eq 2.2. For the ion in vacuo,

Applying q 2.6 gives the Born formula, (2.9)

In this spherically symmetric case, the effect of the interior dielectric constant cancels out. This cancellation does not occur in general. For multiatom molecules, the low dielectric region is defined as the region inaccessible to contact by a 1 . 4 4 sphere rolling over a surface defined by the effective Born radii, RI,of the atoms i. This is consistent with Richards’” definition of the contact and reentrant surfaces. For monatomic ions, the dielectric boundary is independent of the solvent probe radius. van der Waals radii were initially used for RI,but they were ultimately adjusted to reproduce experimental solvation energies. For irregularly shaped molecules, Poisson’s equation is solved by the finite difference method. On a cubic lattice with spacing s, the finite difference form of the Poisson equation is 1

-s2Up’P # - Qp)e(PlP?

= -4*P#

(2.10)

where (p) indicates a lattice point, ( i , j, k),and the summation (30) Zauhar, R. J.; Morgan,R. S. J . Comput. Chem. 1988, 9, 171. (31) Richards, F. M. Annu. Rev. Biophys. Eiocng. 1977, 6, 151.

over (p3 is a summation over the six neighboring points, (i,j, k - l), ( i , j , k + l), (hi- 1, k), etc. The notation e@b’)indicates the average value of e in the region between lattice points p and p! This average is calculated from a lattice of c values that is offset from the 4 and p lattices by (s/2, s/2, s/2) in a manner described by Wachspress.)z Equation 2.10 is solved by a standard overrelaxation algorithm.” The values of Q at the outer boundary of the grid are the Dirichlet boundary condition of the pr~blem.’~These must be chosen to approximate the actual, infinite boundary conditions. Here, these 4 values are estimated by a Coulomb potential of the molecular charge distribution in a uniform medium of dielectric constant 6. The grid of 4 values obtained as the solution dots not provide a separation of the potential into Coulombic and Q* terms. It contains a finite difference analogue of the l/r singularity; Le., a charge q placed on a grid point p in a uniform medium of dielectric constant c produces a potential in the 4 grid at point p which is initially approximated by Q/(0.315~)(Born term with the Born radius set equal to a third of the grid spacing). This becomes infinite in the continuum limit as s 0. However, the Coulombic terms from the solvent and vacuum potentials cancel in eq 2.6, so that the potentials obtained from the finite difference 4 grids can be used directly in place of Q* in eq 2.6. The procedure used for calculating AG, values is as follows: (a) set up a 61 X 61 X 61 grid with 1.0-A spacing centered on the molecule; (b) assign partial charges to the atom centers and distribute the charge of each center between the eight nearest grid points by the procedure of Wachspress; (c) determine the corresponding e grid by a procedure that assigns t,,to points accessible to a solvent probe, and ei to inaccessible poiints; (d) determine the initial values of 4 at the outer boundary of the grid by a Coulomb potential as described above; to speed convergence, set the initial values of Q on the rest of the grid in the same way; (e) solve eq 2.1 for the values of 4 on the grid by the overrelaxation algorithm; (f) to improve the accuracy of the results, a calculation is made on a fine grid by setting up a 61 X 61 X 61 grid with 0.25-A spacing centered on the molecule; (g) repeat (b) and (c) on the fine grid; (h) set the initial and outer boundary values of the new grid by interpolation from the old coarse grid; (i) solve for the 4 gri4 (i) repeat (a) through (i) for the case of the molecule in vacuo by replacing t,, by 1.0; (k) use eq 2.6 and the appropriate 4 grids to calculate AG,. The procedure used here for calculating relative solution pK, for the test compounds is analogous to protein titration calculat i o n ~ ?A~ difference is that, in the protein calculations,u the total charge for the titrating group is localized on a single atom whereas distributed charges are used in the present work. Another difference lies in subtracting different models. For the protein titration calculation^?^ the model compounds are the Nformyl-N-methylamidederivatives of the amino acids whose coordinates are derived from the protein crystallographiccoordinates. In this work, the models for the test compounds, HO(CH&NH3+, HOCH21mH+,and HO(CH2)2COOH,are MeNH3+, ImH+, and MeCOOH, respectively; the coordinates for both the protonated and unprotonated state are obtained from ab initio calculations. In contrast to the relative pK, calculations, the absolute solution pK, of a molecule is determined with respect to the molecule in vacuo (see Scheme I) rather than a model compound in solution. Another difference is that, in determining absolute pK, values, AG, of the protonated and unprotonated molecule is calculated by using the ab initio geometry of the protonated and unprotonated molecule, respectively; in contrast, the ab initio geometry of the protonated or unprotonated molecule is employed in pK, shifts. Thus, conformational changes upon protonation/deprotonation in solution are partially accounted for in absolute pK, calculations (which assume corresponding geometries in vacuo and in solution)

-

(32) Wachsprcsr, E. I. Iterative Solution of Elliptic Systems; Rsntiw Hall: Englewood Cliffs, NJ, 1966. (33) P ~ M W. , H.; Flannery, B. P.; Teukohky, S. A.; Vetterling, W. T. Numerical Recipes. The Art of Scfentiflc Computins Cambridge University Press: Cambridge, U.K., 1986. (34) Jackson, J. D. CIaSsical E~ectmdynamics;Wiley: New York, 1975.

Lim et al.

5614 The Journal of Physical Chemistry, Vol. 95, No. 14, 1991 CH3

OG

TABLE I: CHARMM and OPE3 N o a W P8mmetm

CHI

CHARMM'

HG

?&iucQl

atom or group CH3

H a m 2

OG

HG

Famic r i d

CH3 NZ HZ1,2,3

M

c Add

-

OPLSb

e,

kcal/mol u, A Methanol and Methoxide 3.775 3.8576 0.1811 2.8509 0.1591 3.070 0.400 1.4254 0.0498 U, A

Methylammonium and Methylamine 3.8576 0.1811 3.775 2.8509 0.2384 3.250 1.4254 0.0498 0.400

f,

kcal/mol 0.207 0.170 4.598-2 0.207 0.170 4.598-2

CG NE2/ND1 CD2ICEl HE2/ND1

Imidazolium and Imidazole 3.7418 0.1200 3.750 2.8509 0.2384 3.250 3.7418 0.1200 3.750 1.4254 0.0498 0.400

0.145 0.170

CH3 CG OD 1 OD2 HD2

Acetic Acid and Acetate 3.8576 0.1811 3.910 3.7418 0.1200 3.750 2.8509 0.1591 2.960 2.8509 0.1591 3.000 1.4254 0.0498 0.400

0.160 0.105 0.210 0.170 4.598-2

H CG OD 1 OD2 HD2

Formic Acid and Formate 2.6157 0.045 2.500 3.7418 0.1200 3.750 2.8509 0.1591 2.960 2.8509 0.1591 3.000 1.4254 0.0498 0.400

0.050 0.105 0.210 0.170 4.598-2

0.145

4.598-2

m

CEl

HB2

Imidudium

4-hydroxy methyl imlduole

&

Hz1

HTI

Methyl Amine

=&

IU1 Bch.nolAmine

Figure 1. Fully optimized HF/6-31GS geometries: carbon, oxygen, nitrogen, and hydrogen atom are black, dark grey, light grey, and white circles, respectively.

but not in determining relative pK, values. 2.3. Solute Model and Input Parameters. For most of the calculations, except where noted otherwise, a polar hydrogen model is used; Le., the XRISM and dielectric calculations require the coordinates, charges, and nonbonded parameters of the solute as input. HF/6-31G0 optimized geometries were used in conjunction with various combinations of sets of charges and van der Waals (vdW) parameters. The HF/6-3 1G* geometries, which provide a good approximation to the experimental geom~tries,'~J~ and are similar to the standard geometries of the CHARMM22and O P L P force fields, are depicted in Figure 1. All parameter sets are denoted by two letters; the first refers to the charge set and the second to the set of nonbondcd parameters. For example, QC denotes quantum mechanical charges and CHARMM vdW parameters, CC denotes CHARMM charges and CHARMM vdW parameters, and 00 denotes OPLS charges and OPLS vdW parameters. Nonbonded parameters used in the XRISM calcu(35) Davidson, E. R.;Feller, D. Chcm. Rev. 1986, 86, 681.

CHARMM vdW parameters for methanol, methylammonium, (imidazolium, imidazole), (acetic acid, acetate) correspond to side chain of SER, LYS,HIS, and ASP residues, respectively. CHARMM vdW parameters for methoxide, methylamine, formic acid, and formate are assumed to be the same as those for methanol, methylammonium, acetic acid, and acetate, respectively; for H bonded to C (type HA atom) in formic acid and formate, a vdW radius of 2.6157 A and well depth of 0.045 kcal/mol are assumed. bOPLS vdW parameters for methanol, methylamine, and acetic acid were supplied by J. Gao (private communication). OPLS vdW parameters for methylammonium, and acetate are taken from: Jorgensen, W. L.,Gao, J. J. Phys. Chem. 1986, 90,2174, and parameters for imidazolium, imidazole are from: Jorgensen, W. L. and Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657. OPLS vdW parameters for methoxide, formic acid, and formate are assumed to be the same as those for methanol, acetic acid, and acetate; for H bonded to C (type HA atom) in formic acid and formate, a vdW radius of 2.5 A and well depth of 0.05 kcal/mol are assumed. OPLS hydrogens on heteroatoms require no vdW parameters but, for RISM calculations, they are set equal to the values for the hydrogen of water in the modified TIP3P model.

lations are given in Table I, and radii used in dielectriccalculations are shown in Table 11. The various sets of charges for model compounds are summarized in Table 111, and Mulliken charges for bifunctional compounds are given in Table IV. Dielectric measurements on dried and hydrated protein powders as well as protein solutions over a frequency range, 0-20 GHz, yield dielectric constants between 2 and 4.*% Gilson and H o n i e have computed a dielectric constant between 2.5 and 4 for a hypothetical isotropic protein and Nakamura et a1.39 have calculated dielectric constants ranging from l to 20 for different portions of BPTI. Generally, values ranging from 1 to 5 are used to estimate the dielectric constant of a protein, depending on the method of calculation. Given the uncertainty in the dielectric constant of a protein, we have explored calculating AG, with ci values ranging from 1 to 4. In calculating AG,,the protonated and unprotonated geometry of the model compounds are taken from ab initio calculations. Thus dipolar reorientation is not involved and the dielectric constant of the model compound may (36) Tahahima, S.;Schwan, H.P. J. Phys. Chem. 1965,69,4176. (37) Pennock, B. D.; Schwan, H. P.J. Phys. Chcm. 1969, 73, 2600. (38) Harvey, S.C.; Hoclutra, P. J. Phys. Chcm. 1972, 76, 2987. (39) Nakamura, H.;Sahmoto, T.; Wada, A. Protein Eng. 1988,2, 177.

The Journal of Physical Chemistry, Vol. 95, No. 14, 1991 5615

Absolute pK, Calculations TABLE II: Bora Radii (A)Used in Dkketric Cdculatioac atom or group type CHARMM. OPLSb RashinL adiusted" Methanol CH3 OG

HG

CH3E OH HO

2.165 1.600 0.800

2.119 1.723 0.m

2.780 1.500 1.160

2.165 1.600 0.800

Methoxide CH3 OG

CH3E

oc

2.165 1.600

2.119 1.723

2.780 1.500

CH3E NH3 HC

2.165 1.600 0.800

CH3 NZ HZ1.2

CH3E NH2 HC

2.165 1.600 0.800

2.119 1.824 0.m

2.780 1.500 1.160

2.165 1.600 1.200

2.780 1.500 1.160

2.165 1.600 0.800

Methylamine 2.119 1.824 0.m

Imidazolium CG NE2/NDl CD2)CEl HE2/HDl

C NHl CRlE H

2.100 1.600 2.100 0.800

2.105 1A24 2.105 0.m

2.460 1.500 2.460 1.160

2.100 1.600 2.100 1.200

C NH 1 CRlE H

2.100 1.600 1.600 2.100 0.800

CH3 CG OD 1 OD2 HD2

CH3E C 0 OH 1 H

2.165 2.100 1.600 1.600 0.800

NR

2.105 1.824 1.824 2.105 0.m

2.460 1.500 1.500 2.460 1.160

2.100 1.600 1.600 2.100 0.800

0.m

2.780 2.460 1.500 1.500 1.160

2.165 2.100 1.600 1.600 0.800

2.194 2.105 1.661

2.780 2.460 1 .500

2.165 2.100 1.500

1.710 2.460 1.500 1.500 1.160

1.468 2.100 1.600 1.600 0.800

Acetic Acid 2.194 2.105 1.661 1.684

Acetate CH3 CG ODI/OD2

CH3E C OC

2.165 2.100 1.600

Formic Acid H CG OD 1 OD2 HD2

HA

C 0 OH 1 H

1.468 2.100 1.600 1.600 0.800

1.400 2.105 1.661 1.684 0.m

Formate H CG ODl/OD2

HA C OC

1.468 2.100 1.600

1.400 2.105 1.661

0.29

OG

-0.73

HG

0.44

CH3 OG

-0.1 1 -0.89

CH3 NZ HZ1 HZ2 HZ3

0.43 -0.85 0.47 0.47 0.47

CH3 NZ HZ1 HZ2

0.17 -0.83 0.33 0.33

CG NE2 ND 1 CD2 CEl HE2 HDl

0.74 -0.69 -0.69 0.35 0.35 0.47 0.47

CG NE2 NDl CD2 CE 1 HE2

0.47 -0.71

0.27

0.22

-0.70

-0.58

0.43

0.36

-0.02 -0.98

0.10 -1.10

0.3 1 -0.30 0.33 0.33 0.33

0.46 -0.30 0.28 0.28 0.28

0.20 -0.90 0.35 0.35

0.10 -0.70 0.30 0.30

0.50 -0.54 -0.54 0.33 0.33 0.46 0.46

0.36 -0.30 -0.30 0.20 0.20 0.42 0.42

0.10

0.41 -0.57 -0.49 0.13

0.10 0.30

0.42

0.41 -0.59 -0.49 0.13 0.10 0.44

0.25 0.1 5d -0.65 -0.406 0.40 0.25d

Mcthoxide -0.15 -0.85

Mcthvlammonium 0.25 -0.30 0.35 0.35 0.35

Methylamine 0.10 -0.70 0.30 0.30

Imidazolium

Imidazole CG NE2 NDl CD2/CEl HE2

Methanol CH3

2.165 1.500

Methylammonium CH3 NZ HZ1,2,3

TABLE IIk HF/6-31CL, CHARMM, OPLS, ud Adjmted C h g a RfOUD 6-31GS' CHARMM' OPLS' D2'

1.710 2.460 1.500

1.468 2.100 1.500

.From CHARMM u in Table 1. bFrom OPLS u in Table I with hydrogen radius set to 0. 'From Rashin, A. A.; Namboodiri, K. J . Phys. Chem. 1987, 91,6003. "HARMM radii adjusted to fit experimental solvation free energies of corresponding amino acid functional group. be ap roximated by the value of 2 for electronic polarization This value is used for the model compound in calibrating alone! the charges and radii to reproduce the experimental AG, values. 2.4. Nonpolar Contributiom to AG,. The continuum dielectric method yields only the electrostatic contribution to the solvation free energies. We have made an estimate of the nonpolar contribution by replacing a polar molecule with a nonpolar molecule of comparable size, as illustrated in Table V. For example, the nonpolar contribution to AG, of CH3NH3+is estimated from the experimental AG, of CH3CH3(where N is replaced by C). The nonpolar contributions are found to be less than 2 kcal/mol; the difference between the nonpolar contribution to AG,(BH) and AG,(B) is expected to be even smaller ( I 1 kcal/mol). Thus we have neglected comctions to the electrostatic AG, due to nonpolar contributions. In contrast, Rashin and Namboodiri" includes a nonpolar hydration enrhalpy term; this is estimated as a linear

0.45 -0.30 -0.30 0.225 0.225 0.35 0.35

Imidazole -0.51

0.21 0.14 0.40

0.30 -0.40 -0.40

0.10

Acetic Acid 0.05 0.72

CH3 CG OD1 OD2 HD2

-0s -0.69 0.47

CH3 CG OD 1 OD2

-0.21 0.72 -0.76 -0.76

H CG OD 1 OD2 HD2

0.18 0.52 -0.51 -0.66 0.47

H CG OD 1 OD2

-0.07 0.53 -0.73 -0.73

0.00 0.75 -0.52 -0.61 0.38

0.08 0.55 -0.58 0.45

0.00 0.75 -0.52 4-61 0.38

-0.10

-0.04

0.70 -0.80 -0.80

0.70 -033 -0.83

0.08 0.55 -0.50 -0.58 0.45

0.00 0.75 -0.52 -0.61 0.38

-0.50

Acetate -0.16 0.36 -0.60 -0.60

Formic Acid 0.00 0.75 -0.52 -0.6 1 0.38

Formate -0.15 0.35 -0.60 -0.60

-0.04 0.50.

-0.75 -0.75

0.70 -0.83 -0.83

OHF/6-31GL Mulliken atomic charges. 'See footnotes under Table AG,. "Charges from OZ'H group of sugar ring. 'Charge for CH group. 1. CChargesadjusted to fit experimental

function of the cavity surface area of a hypothetical molecule which is obtained from the polar molecule by a formal switching off of all point charges, but retaining the hard-core radii of the hydrogen-bonding groups. Their calculated nonpolar enthalpy contributions range from -1 to -6 kcal/mol. We also have not introduced corrections for dielectric saturation. Monte Carlo simulations of Jayaram et al.% indicate that dielectric saturation is not important for monovalent cations.

5616 The Journal of Physical Chemistry, Vol. 95, No. 14, 1991 TABLE IV: HF/C31G*MulUken chufes group BH B group BH = HO(CHz)ZNH,+ CH2 0.41 0.32 NZ -0.77 -0.75 OG HZI 0.49 0.46 HG HZ2 CH2 0.31 0.15 HZ3

BH

B

-0.86 0.46 0.47 0.49

-0.87 0.34 0.35

CH2 OG HG CH2

0.32 -0.77 0.46 -0.03

BH

HO(CH2)zCOOH 0.25 CG 0.77 -0.81 OD 1 -0.54 0.49 OD2 -0.72 -0.17 HD2 0.51

0.76 -0.74 -0.78

CH2 OG HG CG NE2

0.44 -0.75 0.48 -0.74 -0.69

BH = HOCHzImH+ 0.33 NDI -0.74 CD2 0.46 CE 1 0.50 HE2 -0.72 HDI

-0.57 0.20 0.14 0.40

-0.70 0.30 0.23 0.47 0.48

TABLE V Nonpolar Contribution to AC, polar molecule CHIOH

(methanol)

CHINHI

(methylammonium) CHaCOOH

(acetic acid) (CH)I(NHh

nonwlar molecule CHZCHZ (ethene) CH3CH3

AG.:.. kcalhol , 1.27

TABLE M: tal dcllcvhted Cas-wuC B.dcitks and Proton AtlWks (kd/apd) at 298 K4 BH MeOH McCOOH HCOOH ImH+

B

PA expt

calcb

MeOMeCOO-

408.5 367.5 HCOO363.8 240.0 Im MeNH3+ MeNHz 228.2 HO(CHJ2HO(CH2)Z- 226.0 NH]" NH2 HO(CH2)ZHO(CHz)z- 350.1 COOH COOHOCH21mH+ HOCH21m 241.8

379.2d 3483 347.6 222.1* 214.11 221.31

TAS 6.6 7.0 7.0 7.8h 8.4 7.9

AG, calCe expt 401.9 360.5 356.8 232.2 219.8 218.1

372.6' 341.5' 340.6 214.31 205.71 213.4

7.8h 342.3

-

7.8h 234.0

#AG8 is negative for the free energy change and PA is negative for the free enthalpy change for B H BH. bProton affinities are calculated at the HF/6-31G* level. 'AG (calc) = PA(ca1c) - T U . dBartmtss, E.; Scott, J. A.; McIver, R. 'f J . Am. Chem. Soc. 1979, 101,6046. Bartmess, J. E.; McIver, R. T. Gas Phase Ion Chemistry; Bowers, M. T., Ed.; Academic Press: 1979, p 87. fKoppel, I. A.; Molder, U. H.; Palm, V. Org. Reactiuity 1985,22, 77. #Lias, S. G.; Liebman, J. F.; Levin, D. J . Phys. Chem. ReJ Data 1984,13. "roton entropy taken from Stull, D. R.; Prophet, H. JANAF Thermochem. Tables 1971.

+

1.80

(ethane) CH3CCH2CH2H

1.20

(CH)Z(CH2)3

0.56

(2-methylpropene)

(imidazolium) (cyclopentene) "Ben-Naim, A.; Marcus, Y.J. Chem. Phys. 1984,81,2016. Cabani, S.; Gianni, P.; Mollica, V.; Lepori, L. J . Solution Chem. 1981,10, 563.

2.5. Analysis and Application of Experimental Data. The calculated and experimental AG, correspond to the process of transferring the molecule from a fixed position in the gas phase into the same fixed position in the liquid at 298 K and 1 atm." Since the experimental AG, of imidazole has not been measured, we have used the experimental AGEof 4methylimidazole measured by Wolfenden et al." to approximate that of imidazole. The AG, value reported by them corresponds to pH 7 where 50% of 4methylimidazole (pK, = 7) is ionized. However, the actual measurement was made for a higher pH where only the uncharged molecule is present; we use the latter value to obtain AGa. Experimental solvation free energies for the charged species are not available and were derived from the experimental AGEvalues of their neutral counterpart, the experimental values for AG8, AG,(H+) = -259.5 kcal/mol, and the experimental pK,; Le., eq 1.1 was used to obtain AG,(BH). For charged species, only the difference in AG, between two different ions (e.g., AG,(ImH+) AG,(MeO-)) can be determined directly from experiment by means of eq 1.2 due to the uncertainty in AG,(H+) (see Introduction). The experimental pK, values of the isolated model compounds are close (within a pK, unit) to the pK, values of the comsponding amino acid side chains," suggesting that the present analysis and parametrization should be useful for proteins. The calculated solution pK,'s are obtained from eq 1.1 with the experimental AGg and AG,(H+) values and the calculated solvation free energies. Care must be used in applying experimental data to yield solvation free energies as defined above. Experimentally, the partition coefficient, (ptlp:),, of solute u between the gas phase and dilute solution at equilibrium is measured. If the gas and solution are assumed to be ideal, AGs can be evaluated from the relation

-

AG, = Apou,rol= -RT In

Lim et al.

(a),

(2.11)

(40) Ben-Naim, A. J . Phys. Chem. 1978,82, 782. (41) Wolfenden, R.;Andemson, L.; Cullis, P.M.Bfochcmisrry1981, 20, 849.

where Apourolis the difference in standard chemical potential of the solute between the two phases using the molarity scale and p i and p i are the number densities in solution and in the gas phase, Experimental data corresponding to a different standard state must be transformed to yield the proper molarity scale standard state. For example, the data of Jones and Amett4z are in the form of Henry's law constants, K, (mmHg/molal) = Pu/mu,where Puis the measured partial pressure (mmHg) of the solute u above an aqueous solution of molality mu (mol/kg of solvent). The solution standard state is the hypothetical unit molal solution and the gaseous standard state is the real gas at 1 atm. Their reported solvation free energy for methylamine (denoted by the subscript JA), (AG,)jA = R T In (K,/760) is -2.68 kcal/mol. The correct AGa can be determined by substituting p: = NAmuand pf = Pu/kOT= K,mu/kBT into eq 2.1 1 to obtain AGrolv= -RT In

NAmu

-

~mmu/k0~ -RT In (RT/760)

+ (AGrolv)jA(2.12)

Note that K, and R T in parentheses in eq 2.12 should have the same units; thus for RTin J/mol, K, (mmHg/molal) is O.1333Km (J/mol). The AG, calculated from eq 2.12 for methylamine is -4.6 kcal/mol, instead of (AGJjA. 3. Results 3.1. Cas-Flmse Proton AfRnitk and BLpidtiea Table VI gives the experimental and calculated proton affinities and gas-phase basicities in kcal/mol. The calculated proton affinities for methoxide, formate, imidazole, and methylamine were taken from Brooks et al.43 The calculated gas-phase basicities are obtained from the proton affinities by including an estimate of the entropic contribution; it is equal to the difference between the experimental gas-phase basicity and proton affinity. When the latter is not available, the entropy change is approximated by the entropy of a proton, 26.012 eu," (TAS = 7.8 kcal/mol at 298 K) which is calculated from statistical mechanical formulas for the translation and electronic entropy; only the ground state is included for the electronic entropy since the first excited state lies at 82 258 cm-' above the ground state. The remaining entropic contribution,S(B) - S(BH+), is usually relatively small, ranging from 0 to 5 eu, and the methods (statistical mechanics or isoelectronic approximation) used to estimate this difference are only good to f 2 eu.45*46

Jonts, F. M.;Amett, E. M.Prog. Phys. Org. Chcm. 1914, I I , 2b3. (43) Brooks,111, C.; Brunger, A,; Fnncl, M.;Haydock, K.;A h , L. C.; Karplua, M.Ann. N.Y.Acad. Sei. 1966,471,295. (44)Stull, D.R.;Prophet, H. JANAF Thermochcm. Tables 1971. (42)

The Journal of Physical Chemistry, Vol. 95, No. 14, 1991 5617

Absolute pK, Calculations

------------------

BH + B. Me.0.H Me.0 Me.0 Me.0.H Mc.0.H Me.0 Me.0 Mc.0.H Mc.0.H Mc.0 Me.0 Me.0.H Me.0.H -. Me.0 Me.0 Me.0.H C.H.3.0 Mc.0.H C.H.3.0 Me.0.H MeOH MeOMcOMcOH Me.N.H.2 Me.N.H.3 Me.N.H.2 Me.N.H.3 Me.N.H.2 Me.N.H.3 Me.N.H.2 Me.N.H.3 Me.N.H.2 Me.N.H.3 Me.N.H.3 Me.N.H.2 Me.N.H.2 Me.N.H.3 MeNHS' McNHZ ImH.7 Im.6 Im.6 ImH.7 Im.6 ImH.7 ImH.7 Im.6 Irr.Y.7 Im.6 Im.6 ImH.7 ImH.7 -. Im.6 ImH+ Im ImH+ Im H.C.O.0 H.C.O.0.H H.C.O.0 H.C.O.0.H H.C.O.0.H H.C.O.0 H.C.O.0 H.C.O.0.H HCOOHCOOH Mc.C.O.0 Me.C.0.O.H Me.C.O.0.H -. Me.C.O.0 Me.C.O.0 Me.C.0.O.H MeCOOMeCOOH

parameters*

* oc cc cc cc cc

CC' CC'

00 00

Qc QC

cc cc cc

00 00

+

QC QC

cc cc cc

00 00 QC QC

cc cc cc cc cc

methode D1 D4 D1 D2 D3 D4 D1 D4 D1 D4 E E D1 D4 D1 D2 D4 D1 D4 E D1 D4 D1 D2 D4 D1 D4 E E D1 D4 D1 D4 E " D1 D2 D4

AG,(BH)~ -1 1.8 -4.7 -9.8 -6.4 -4.8 -3.8 -3.8 -1.5 -9.8 -3.9 -5.13 -5.13 -103.0 -83.6 -102.0 -92.3 -83.9 -88.6 -77.2 -72.9 -76.7 -65.2 -73.7 -68.4 -64.1 -71.3 -62.4 -64.6 -65.0 -15.1 -6.4 -12.4 -5.5

AGSB)'

-9.9 -6.7 -4.2 -6.7"'

-75.2 -71.1 -67.3 -82.2

E

-92.6 -79.2 -90.3 -83.9 -80.5 -78.5 -90.3 -78.5 -88.4 -74.9 -97.7 -96.3 -8.8 -3.5 -7.3 -4.8 -2.9 -1.2 -2.7 -4@

-10.6 -4.4 -7.7 -5.3 -3.3 -10.0 -4.4 -10.2' -10.2' -81.9 -73.8 -77.7 -71.8

PK,/ 23.7 28.2 23.9 26.1 27.4 28.2 19.5 26.5 25.3 30.8 15.0' 16.0' 29.3 19.3 30.0 24.7 20.0 20.2 15.2 lo.@& 15.3 11.4 15.3 13.1 11.4 11.8 9.4 6.7" 7.e 10.5 10.1 11.6 10.8 3.7Y 12.3 12.9 13.8 4.81-

.Different interacting sites are separated by periods, H.n means that n H sites are included, and 1m.n or 1mH.n denotes an n-site polar hydogen model. *QC denotes quantum mechanical charges and CHARMM vdW parameters, CC denotes CHARMM charges and vdW parameters, CC' denotes CHARMM charges and vdW parameters for the hydroxyl group of a sugar ring, 00 denotes OPLS charges and vdW parameters. CDn denotes dielectric method using el = n, E denotes experimental values. "Free energy of solvation for BH(g) BH(s) at 298.15 K and 1 atm. 'Free energy of solvation for B(g) B(s) at 298.15 K and 1 atm. fSolution pK, for BH(s) B(s) + H+(s) at 298.15 K and 1 atm calculated from experimental AG, and calculated AG,. #Ben-Naim, A.; Marcus, Y. J . Chem. Phys. 1984,81, 2016. "abani, S.; Gianni, P.; Mollica, V.; Lepori, L. J. Solution Chem. 1981, 10, 563. 'Bartmess, J. E.; Scott, J. A.; McIver R. T. J. Am. Chem. Soc. 1979, 101, 6046. jEbel, H. F.Die Acidifdr der Clf-Saiiren; Georg Thieme Verlag: Stuttgart 1969. kSmith,R. M.; Martell, A. E. Critical Stability Constants; Plenum Press: New York, 1989; Vol. 2, Suppl. 2. 'Wolfenden, R.; Anderson, L.; Cullis, P. M.; Southgate, C. C. B. Biochemistry 1981, 20, 849. mPearson,R. G. J . Am. Chem. Soc.

-

-

-

1986, 108,6109.

Experimental gas basicities are not available for the compounds with two functional groups except for amino ethanol. The order of the calculated proton affinities and gas basicities for monofunctional groups is in agreement with experiment; i.e., MeOH > MeCOOH > HCOOH > ImH+ > MeNH3+. The relative proton affinities differ from the experimental values by 2-15 kcal/mol but the absolute values differ by as much as 14-29 kcal/mol.. In all cases the calculated value is larger than the experimental result. This difference could be reduced to within experimental error b using extended basis sets and including electron correlation;4r49 e.g., a linear regression analysis of the proton affinitiescalculated at the MP2/6-31GS level for 11 neutral compounds gives PA,lp.= 0.48 i 0.995PAak. For our present purposes, we are mainly interested in relative changes, as described above. 3.2 Sdntion Fne bergla and sohrtion pK,. Table VI1 gives the solvation free energies AG, and pK, at 298 K as a function (45) Bartmcu. J. E.;Scott, Jr., J. A,; McIver, R. T. J . Am. Chcm. Soc. 1979,101,6046. (46) Bartmcu, Jr., J. E.;Mclver, R. T. In Gas Phase Ion Chemistry; Bowers, M. T., Ed.; Academic Prcu: New York, 1979. (47) Del Bene, J. E. J . Compur. Chem. 1985.6, 296. (48) DeFrca, D. J.; Mcban, A. D. J . Comput. Chem. 1986, 7, 321. (49) Koppel, 1. A.; Molder, U. H.; Palm, V. A. Org. Reactfvlty 1985, 22, 3.

of the parameters and the methods employed. The abbreviations used in columns 2 and 3 of Table VI1 are explained in the footnote and in section 2. Overall, AG8 for the charged species with the various models are in qualitative agreement with experiment. As expected, increasing the dielectric constant decreases the solvation. For the uncharged species, qualitatively correct values are also obtained, although the percentage error is larger. The dielectric method generally yields the correct ordering of pK, such that pK,(MeOH) > pK,(MeNH3+) > pK,(ImH+) > pK,(MeCOOH) > pK,(HCOOH); exceptions are the parameter sets QC and CC using cI = 1 where the order of the pK, of MeOH and MeNH,+ is reversed, and for the all CHARMM parameter set with eI = 4 where the pK, of ImH+ is less than that of MeCOOH. The dielectric method overestimatesthe pK, values by roughly a factor of 2 and, in most but not all cases, overestimates AG, of cations by HI%, but underestimates AG, of anions by 523%. The AG, of methylamine and acetic acid obtained by using cf = 2 and CHARMM charges and radii are in agreement with experiment. AG,(MeOH) is overestimated by 1.3 kcal/mol whereas AG,(Im) is underestimated by almost a factor of 2 but is close to its experimental value for OPLS charges and radii with q = 1. As noted in the Introduction, a relatively small error in the calculated AG, can lead to a large error in the computed pK,. For example, even though the dielectric method using an cf of 2 reproduces the experimental AG, of acetic acid and undmtimates AG, of acetate

Lim et al.

5618 The Journal of Physical Chemistry, Vol. 95, No. 14. 1991

by only 13.596, the calculated pK, of 12.9 for acetic acid is 8.1 pK, units larger than the experimental value of 4.8. For small values of cf (Dl),there is an overestimate of AG,. This d d be due to dielectric saturation effects which would lower the dielectric constant in regions near the ion and thus decrease the calculated solvation free energies. However, the results of a Monte Carlo based free energy simulation of cation solvationx show that the effect of dielectric saturation is opposed by the tendency of waters to crowd closer to the ion as the charge is increased (electrostriction) which tends to increase the solvation free energy. As a result, dielectric saturation was found to be insignificant for charges less than or equal to unity.% A possible rationalization of the observed overestimation of AG, for cations is based on the effective Born radius, which separates the low and high dielectric media in the simple Born model. Roux et a1.I6 derived an asymptotic solution to the XRISM equations in the limit that a monatomic solute ion is large compared to solvent, and interpreted the Born radius in terms of the first peak in the solute-solvent radial distribution function. Thus the low dielectric region of a spherical monatomic cation is underestimated since the Born radius is equated to the atom's van der Waals radius instead of the equilibrium cation-water oxygen distance. Using a similar argument, one would expect solvation free energies of monoanions to be overestimated, but to a laser extent due to the smaller hydrogen radius of water. Instead, the calculated anion solvation free energies are underestimated in most cases. It should be noted that for arbitrarily shaped molecules, the "effective" Born radius may not be simply related to the first peak in the solutesolvent radial distribution function. Using essentially the same macroscopic dielectric method and an cf equal to 2,Gilson and Honig' have calculated AGs of -4.3, -77, and -70 kcal/mol for methanol, methylammonium, and acetate, respectively. We have calculated AG, for the same set of parameters (CHARMM charges in Table I11 and CHARMM radii in Table I1 except r(H) = 1.2A and r(0)= 1.5 A) and found AG, for methanol and acetate in exact agreement with the values obtained by Gilson and Honig; AGs for methylammonium, -79.0 kcal/mol, is close to their calculated value, -77.0 kcal/mol. Jorgensen and Briggs5have performed Monte Carlo simulations and computed the difference in solvation free energies of methylamine and methanol to be 1.7 f 0.1 kcal/mol compared to the experimental difference of 0.5 kcal/mol. For the same set of charges and van der Waals radii, AG,(MeNH,) AG,(MeOH) calculated from continuum dielectric methods is 1.2kcal/mol for cf = 4,and 2.6 kcal/mol for cf = 1; XRISM with nonzero van der Waals parameters on hydrogen gives 3.1 kcal/mol. As noted by Gilson and Honig,3 the dielectric method yields results similar to those obtained by more accurate free energy calculations; both methods give values which are in close but not exact agreement with the experimental values. However, the dielectric method is relatively fast; the calculation of AG, for each small molecule takes only 10 to 15 minutes on a CONVEX C1, as compared to 8 h for a free energy simulation? In summary, the results in Table VI1 and those cited in the Introduction show that both microscopic methods and macroscopic methods 4ve solvation free energies for ionic species that are close to "experimental" values. Nevertheless, when these values are employed to calculate solution pK,, large errors in pK, result with various standard choices of parameters as illustrated in Table VII. To obtain a formulation that can be used for calculating relative pK, in other systems, e.g., in enzyme-substrate complexes, we choose a consistent set of parameters that yields the absolute solution pK, of the model compounds. We use continuum dielectric methods and adjust radii and charges to reproduce the experimental AG, and relative AG,. We assume q q u a l to 2 to approximate the electronic polarizability of the solute in solution where reorientational polarizability effects are presumably negligible. In trying to fit the experimental AG,, radii of polar atoms and hydrogens are first adjusted to obtain as close an agreement with the experimental values without adjusting the charges. If necessary, charges are subsequently adjusted to reach agreement with the experimental values. The adjusted radii and the adjusted

-

TABLE MU SdvatiOa Free Empks (kcal/mol) ud Cu ud PK, 8t 298 Kg

---

AGs-

BH+B

(BH)

MeOH MeOMeNH3+ MeNH2

-5.15 -73.3 -64.4 -6.7 -7.0

ImH+

Im

AG,-

(B) PK,W pK.(db APK,' -96.6 -4.8 -10.1 -82.5 -84.5

15.9 10.8 6.7 4.5 2.6

273.1 150.8 157.1 250.3 249.6

257.2 140.0 150.4 245.8 247.0

MeCOOH MeCOOHCOOH HCOO"Solvation free energies calculated with adjusted charges and radii (see Tables I1 and 111) using the dielectric method with t, = 2. Experimental AG values from Table VI divided by 1.3644 to convert to pK, units. CDif&rencebetween experimental gas pK, and calculated solution pK,.

8

charges are shown in the last column of Table I1 and of Table 111, respectively. The effective Born radii of all atoms of neutral molecules are unchanged from their CHARMM values as well as the radius of all carbon types. The radius of 0 is decreased from 1.6to 1.5 A for anions and the radius of H bonded to cations is increased from 0.8 to 1.2 A. These adjustments are consistent with the oxygen and hydrogen cavity radii of Rashin and Honigm who defined the Born (cavity) radius as the radius of a sphere which contains a negligible electron density contribution from the surrounding solvent (seeTable 11). The adjusted charges are all within f0.15 of either their CHARMM or OPLS values. The largest change in the charges is required for methylammonium. Although increasing the hydrogen radius of MeNH3+ further may yield better agreement with experiment without changing the charges on MeNH,', it would worsen the agreement of AG,(ImH+) with experiment. Consequently, the charges on the hydrogens of MeNH3+were reduced and the remaining charge was added to the methyl group since the other alternative of decreasing the nitrogen chaige led to an increase rather than a decrease in AG,. Table VI11 gives the optimized AG, and solution pK,, the gas-phase experimental pK,, and the difference between the gas-phase and solution pK, values. The gas-phase pK, values are lowered by 247-257 pK units (335-351 kcal/mol) for ionization of acids, and 140-150 pK units (191-200 kcal/mol) for ionization of bases. 3.3. Results for Bifunctional Systems. Table IX contains the test results for the three bifunctional compounds. In each, the pK, shift, ApK,, relative to the model compound without the CH20H group was examined. Ideally, we would have tested the adjusted parameters for a hydroxyl group as well but we could not find data for the ionization of a hydroxyl group in a bifunctional compound. The geometry of either the protonated or the unprotonated state was used in the calculation of relative pK, values (see Table IX and section 2.2);the various geometries are shown in Figure 1. The self-contributions (Born) and the cross contributions (Coulomb) to the calculated pK, shift, as defined by eqs 2.2,2.6,and 1.2,are also given in Table IX. The calculated pK, difference between MeNH3+ and HO(CH2)2NHp+is -0.85, in close agreement with experiment. The observed decrease in pK, is due mainly to the Born term which disfavors protonated aminoethanol since the volume of low dielectric medium is increased relative to methylammonium. In principle, the Coulomb term should be more negative since the conformational change upon deprotonation of HO(CH2)2NH3+,which was not allowed for in these calculations, would contribute a favorable Coulomb term due to the internal hydrogen bond (see Figure 1). The sensitivity of the relative pK, to the conformation of the molecule becomes even more apparent when the change in pK, is determined with two different geometries. When the protonated geometry of 4hydroxymethylimidamle is used, the calculated pK,, -0.64,is in excellent agreement with experiment. However, with the geometry of 4-hydroxymethylimidazole, the calculated pK, is twice as large due to the Born and Coulomb terms, both of which have doubled. With the geometry of the unprotonated state, the Coulomb contribution results in a greater decrease in the relative (50) Raahin, A. A,; Honig, B. J . Am. Chcm. Soc. 1985, 89, 5588.

The Journal of Physical Chemistry, Vol. 95, No. 14, 1991 5619

Absolute pK, Calculations

TABLE IX: Rchtive sdutioa pK, at 298 K Uaiog Dkkctric Methods@ BH-B HO(CH2)2NH,+ HO(CH2) NHz" HOCH21mH+ H O C H 2 1 ~ c * ~ HOCH21mH+ HOCHzImd*f HO(CH2)2COOH HO(CHJ2C00-dJ HO(CH2)zCOOH HO(CH2)2CW-'J

--

+

+

+

Born -0.83 -0.22 -0.48 0.44 -0.94

calculated ApK, Coulomb -0.02 -0.42 -0.95 -2.37 1.60

total -0.85 -0.64 -1.44 -1.93 0.66

exptlb ApK, -1.10 -0.64 -0.64 -0.40 -0.40

'Using adjusted charges and radii and e, = 2. bSmith, R. M.; Martell, A. E. Critical Stability Constants; Plenum Press: New York, 1989; Vol. HF/6-31G* geometry of unprotonated state. eRelative to pK, of methyl2, Suppl. 2. Wsing HF/6-31GS geometry or protonated state. ammonium. /Relative to pK, of imidazolium. 'Relative to pK, of acetic acid. TABLE X Solvation AG, (kal/mol) as a Function of Total Cbarge -AG, species Q XRISM D1 or D4 expt 100-80 -1 132-92 93-72 anions 78-58 1 53-32 103-62 cations 15-2 11-3 neutrals 0 5 to -14

pK, since the interaction between the hydroxyl hydrogen (q(H) = 0.36) and the deprotonated nitrogen (q(N) = 4 . 4 9 ) favors 4-hydroxymethylimidazole, but the interaction between the hydroxyl hydrogen (q(H) = 0.36) and the protonated nitrogen (q(NH) = 0.12) favors 4-hydroxymethylimidazolium.In the case of 3-hydroxypropanoic acid, the experimental difference is 4 . 4 0 compared to the calculated pK, change of -1.93 and 0.66 for the unprotonated and protonated geometry of 3-hydroxypropanoic acid, respectively. The negative and positive signs in the Coulombic pK, change indicate that the unprotonated and protonated states are favored, respectively, due to the presence of an intemal hydrogen bond in these states. The discrepancy in calculated pK, shifts may be due to the fact that, in solution, intervening water molecules will attenuate the internal hydrogen bond in the cyclic gas-phase structure. This is supported by gas-phase experiments5' which show that the exothermicity of ring closure and the temperature needed for thermal opening of the hydrogen-bonded ring of H2N(CH2)2NH3tis greater than when the ion is solvated by 1-4 water molecules. Thus, the calculated pK, difference is very sensitive to the conformation of the molecule in solution, which may not be the same in vacuo. In a continuum dielectric calculation of apparent pK, values for geometries of lysozyme obtained from the triclinic and tetragonal crystal forms, the results differ by up to 3.6 pK, units; in several cases, this could be attributed to conformational effects.2' 4. Concluding Mscussion This paper is concerned both with the qualitative aspects of solvation and with quantitative models needed for accurate pK, evaluations. Table X summarizes the calculated and experimental solvation free energies in water as a function of the total charge. The order of magnitude of the solvation free energies is as expected from the simple Born equation, eq 2.9, which predicts that AG, for a neutral atom is zero and AG,for a spherical cation or anion is -163.95/R, where R is the Born radius. Continuum dielectric methods yield solvation free energies which are often within 10-20% of the experiment estimates. However, these relatively small errors can result in calculated pK, which are roughly twice the experimental value. Because of the difficulty of absolute pK, calculations, there is a need to parametrize the atomic charges and radii to reproduce the experimental solvation free energies. These can then be applied to proteins. A set of optimized charges and radii was obtained for four commonly occurring functional groups, viz., a hydroxyl, an ammonium, an imidazolium, and a carboxyl group. These were tested on three bifunctional compounds. The adjusted parameters, which differ only slightly from standard values, yielded pK, changes for the test system which were in reasonable agreement with experiment. However, the pK, changes were found to be (51) Meot-Ner, M.;Hamlet, P.;Hunter, E. P.;Field, F.H.J . Am. Chem.

Soc. 1980, 102.6393.

TABLE XI: Solvation AC, (kul/mol) and Solution pK, at 298 K Calculated Using XRISM' BH+B parameters AG.(BH) AG,(B) pK. -5.1 -131.9 -10.1 Me.0.H Me.0 QC 6.3 -127.3 -15.0 cc Me.0 Me.0.H 1.8 -2.1 -1 12.7 00 Me.0 Me.0.H -3.7 -50.3 -1.6 Me.N.H.3 Me.N.H.2 QC 0.2 -5 1.7 2.4 cc Me.N.H.3 Me.N.H.2 0.1 -52.9 1.o Me.N.H.3 Me.N.H.2 00 -5.1 1.2 -37.1 ImH.7 Im.6 QC 0.9 14.0 -32.3 cc ImH.7 Im.6 6.9 12.3 -42.3 00 ImH.7 Im.6 -1.2 -101.0 -1 3.8 H.C.O.0 H.C.O.0.H QC -0.2 -91.7 -7.6 H.C.O.0 cc H.C.0.O.H -0.8 -105.4 -17.2 00 HC.O.0.H. HC.O.0

-------

@Seefootnotes under Table VII.

sensitive to the conformation in solution. We plan to employ the set of optimized parameters in titration calculations of the RNase A-CpA complex. We hope that the results will be more reliable than those obtained with charges and radii which have not been optimized for such dielectric calculations. Further, the protein calculations together with the results for ion-pair systems, such as MeNH3+-HCOO- representing a lysine interacting with an aspartic acid, should aid in elucidating the possible roles of catalytically important residues in ribonuclease A.

Acknowledgment. We thank J. Gao for many useful references,

K.Haydock for providing the coordinates for some of the ab initio calculations, M. Sommer for helping with the dielectric program, and H.-A. Yu for helping with the XRISM program and for stimulating discussions. We are grateful to Cray Research for the use of their computing facilities. We thank E. Wimmer and J. Mertz for their assistance at Cray Research. This work was supported in part by a grant from the National Science Foundation. Appendix A. RISM Integral Equation Method. We briefly outline the extended reference interaction site model (XRISM) integral equation formalism here and rely on the original references to provide the details. For a solute (u) dissolved in a solvent (v) at infinite dilution (p, 0, denoted by superscript (0)), the w and uv parts of the XRISM integral in Fourier k space (denoted by a caret) are

-

=

a$)+[fi,e$)

+ Q$)$$)][fi,,+ ,&$)I

(A.1)

Qg) = fi#@,[41rV + ,,Q$]

(A4 The site-site HNC closure in r space adopted in this study is hh:)(r) = exp(-@trc,(r) h$)(r) - c$) (a+: A&)+,(r)l - 1 (A.3) where {a)= {u,v). In eqs A.I-A.3,6 is the matrix of intermolecular total correlation functions; is the chain sum of and @ = A$ with &,(k) = 4 r @ q j 7 / k where @ = l/kBTand q,, is the charge on site a; 41r is the matrix of intramolecular correlation

+

+

B

(52) Rossky, P. J.; Dale,

W.D.T.J. Chem. Phys. 1980, 73, 2457.

J. Phys. Chem. 1991,95, 5620-5624

5620

functions; for rigid molecules, $&) = sin (kLq)/kLq, where Lw is the distance between sites a acd y ; S is the matrix of direct correlation functions and S* = 3 - k’;kvis the matrix of solvent density; and V is a short-range potential which usually, but not necessarily, takes the Lennard-Jones 6-12 form. The symbol A is a numerical factor introduced to impose consistency between the calculated solvent correlation functions and a given value of the solvent dielectric constant, e, It is given bysFss

and p ~ =p0.997 g/cm3, the A of eq A.4 is 0.959 to match the experimental water dielectric constant of = 78.4.60 B. Results a d Discussion. The XRISM results in Table XI show that the order of the solution pK, is pK,(MeNH3+) > pK,(ImH+) > pK,(MeOH) > pK,(HCOOH) for parameter set QC; pK,(ImH+) > pK,(MeNH3+) > pK,(HCOOH) > pK,(MeOH) for parameter set CC; and pK,(ImH+) > pK,(MeOH) > pK,(MeNHf) > pK,(HCOOH) for parameter set 00.Thus for all three parameter sets, XRISM fails to predict the order of the absolute and relative pK, values; this results from XRISM ovemtimating AG,of anions by 16-36% and underestimating AG, of cations by 27-51%. It also underestimates AG, of neutral compounds and often yields positive AG, values. These trends are also manifested in AG, of water, Cl- and Na+ obtained by Yu and Karplus,20though the errors are smaller (-10%). A reason for the discrepancy between the calculated and experimental AG, can be inferred from the work of Yu and Karplus.2o They have used XRISM to obtain the solute-solventradial and a distribution functions for a cation (g+), an anion (g), with water oxygen and Lennard-Jones nonpolar particle (BO), hydrogen sites, respectively, at 300 K the nonbonded parameters for cation, anion, and nonpolar particle are the same and correspond to those of a chloride ion. The nearest-neighbor peak in the pair distribution function of a cation with the oxygen site of water, g,-,+, is positioned 1.28 A further than that of an anion with the hydrogen site of water, gH-. Thus XRISM would predict a more favorable solvation energy for an anion relative to a cation, which lacks favorable hydrogen-bonding interactions. The overestimation of anion AG,and underestimation of cation AG, may be mainly due to the inadequacy of the HNC closure.llJ” Furthermore, a comparison of the enthalpic and entropic contribution to AG, of the model chloride neutral atom shows that both terms are unfavorable with the entropy term dominant. Since solvation free energies for neutral molecules are underestimated, and their experimental values are all negative, it appears that the magnitude of the entropy term is overestimated and/or there is a lack of stabilizing energetic interactions with the solvent which would make the enthalpy negative. The latter is supported by an integral equation calculation” of the solvation free energy of methane using the Percus-Yevick closure where a large c(CH4-O) of 0.412 was needed in order to obtain the correct experimental value.

where (d,?) is the average squared dipole moment of the solvent molecule. By solving eqs A.l-A.3 self-consistently one can obtain hz)(r) and c$(r) from which the excess solvation free energy, AG,,of introducing a solute into the solvent at infinite dilution can be calculated by using the e x p r e s ~ i o n ~ ~ , ~ ’ P”“

4

AG, = -E E J d i {‘/2[h$?,!(r)]2 - @(r) - f/2hLv(r)c$(r)) Bu-lT-I

(A.5) where nu is the number of solute sites and n, is the number of solvent sites, equal to 3 in the present case. The solvent model used in the XRISM method is the rigid threasite TIP3P water with minor modifications of the hydrogen Lennard-Jones interaction parameters to avoid a Coulombic singularity in the interaction potential due to the absence of repulsive cores on the hydrogen sitesas9 The empirical potential is characterized by the experimental water geometry (ROH= 0.9572 LHOH = 104.52O) and charges of -0.834 on the oxygen and 0.417 on the hydrogen. The Lennard-Jones parameters are = 3.15 A, 0.152 kcal/mol) for the oxygen site and (0.40 0 046 (kcal/mol) for the hydrogen site. For the Lennard-Jones interaction between oxygen and hydrogen sites, the combination and uq = (au u,)/2, are used. At 25 OC rules, em =

tC!

+

(53) Cummings, P. J.; Stell, G.Mol. Phys. 1981, 44, 529. (54) Cummings, P. J.; Stell, G. Mol. Phys. 1982.46, 383. (55) Rouky, P. J.; Pettitt, B. M.;Stell, 0.Mol. Phys. 1983, 50, 1263. (56) Zichi, D. A.; Rouky, P. J. J. Chrm. Phys. 1986,84, 1712. (57) Singer, S. J.; Chandler, D. Mol. Phys. 1985, 55, 621. (58) Jotgemen, W. L.; Chandrasekbr, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J . Chcm. Phys. 1983, 79, 926. (59) Pettitt, B. M.;Rwky, P. J. J . Chcm. Phys. 1982, 77, 1451.

(60)Owen, B. B.; Miller, R. C.; Milner, C. E.; Cogan, H. L. J. Chcm. Phys. 1961, 65, 2065.

Transport Phenomena and Asymmetry Effects In Membranes wlth Asymmetric Fixed Charge Distributions J& A. Manzanares, Salvador MaG, and Julio PeUicer* Departamento de Termodinbmica, Facultad de I;fiFica,Universidad de Valencia, 46100 Burjarsot, Spain (Received: October 5, 1990) Transport phenomena through membranes with asymmetric fixed charge distributionsvarying linearly with position inside the membrane are theoretically studied. The limits of applicability of previous models based on Henderson’s assumption have been established, and the effects that the asymmetry exerts on the flux and the membrane potential have been found to be of minor importance for the physical model employed. The theory is basal on the Nemst-Planck equations, and the numerical solution dcriwd doss not show the doubtful results arising from the uw of the approaimated Hcndcr”s amumption: the nonzero, steady-statevalues of the flux and the membrane potential for a situation in which two identical external solutions are separated by an asymmetric membrane, and the Occurrence of steady-state reverse transport. Some comments on the temporary nature of the experimentally observed asymmetry effects are also included. the existence of a structural or chemical inhomogeneity in a membrane requires the revision and extension of the classical

Introduction The problem of ion transport through membranes with i n h e mogeneous charge diatributions has received some attention recently.’+ The description of the physical phenomena related to OO22-3654/91/2095-5620$02.50/0

( 1 ) Man, s.; Manzanarca, J. A.; Compsfi, v. An. Fh.,in preu. (6

1991 American Chemical Society