- Email: [email protected]
The effect of hydrogen bonding on protonation energies and ion pair formation
Chemical Physics ELSEVIER
Chemical Physics 204 (1996) 403-410
The effect of hydrogen bonding on protonation energies and ion pair formation Markus P. Filscher a Depurtment ofTheoretical b Depurtment
ofPhysiology
Chemistry,
University
and Biophysics,
a, Ernest L. Mehler b,* ofLund,
Chemical
Mount Sinai School
Centre, P.O.B. 124, 2210 Lund, Sweden
of Medicine,
CUNY, New York, NY 10029,
USA
Received February 21, 1995
Abstract The effects of environmental hydrogen bonding by solvent molecules on the stability of charged groups and ion pairs have been analyzed by ab initio calculations on molecular models representing acidic and basic amino acid residues. By constructing closed reactions cycles, the energetics of the isolated and solvated groups can be directly compared and the origin of the changes in relative stability exhibited. The results of the calculations showed that the solvent molecules stabilize interacting pairs in both the neutral and charged states, but that the effects shift the equilibrium towards the latter. The analysis further indicated that the main contribution to this effect is due to the changes in the work required to form the charges attendant upon solvation. Moreover, although the relative stability of the neutral and charged states is primarily controlled by electrostatic effects, it can be modulated by additional short range interactions when the species are close to each other. In such cases the relative stability of the states predicted by the total interaction energy may be different to what
is obtain from just the electrostatic interactions.
1. Introduction The stability of ion pairs is of fundamental importance in protein structure and function. The majority of ionizable groups are found at or near the protein surface where they can form dynamic networks of hydrogen bonds involving solvent molecules and possibly proximal polar residues. However, some titratable groups are found in the protein interior, mostly in pairs, but sometimes as isolated charges [ 1,2], and they are often involved in stabilizing the active site region or participate directly in the catalytic mechanism [3,4]. To explain their effects in
??
Corresponding author.
0301-0104/96/$15.00
0 1996 Elsevier Science B.V. All rights resewed
SSDI 0301-0104(95)00423-S
terms of the attendant free energy changes at the molecular level, and to evaluate them quantitatively by computational methods are presently major challenges to theory. Molecular simulation techniques are one way to obtain free energies, but it is necessary to include a large part of the environment in the model for such calculations to approach an acceptable degree of reliability. These methods also have the disadvantage that it becomes very difficult to identify the individual interactions which lead to the observed macroscopic effects or to assess their significance [5,61. In spite of the obvious limitations, analysis of the interactions leading to energy changes seems to be most meaningful when carried out on models that are not dependent on experimentally determined parameters, and are flexible enough to
404
M.P. Fiilscher, E.L. Mehler/Chemical
account for the relevant contributions to the interaction energy. Due to the long range nature of Coulombic forces, the intermolecular interactions which maintain protein structure are dominated by electrostatics. At shorter intermolecular separations, and especially those where chemical processes take place, additional effects of quantum chemical origin become important and may even be dominant. These include charge transfer, exchange and correlation effects that sometimes lead to unexpected results not predictable by a consideration of electrostatic effects alone. Therefore, to provide a proper theoretical accounting of protein structure and function, it is essential that these interactions are described in a reasonably reliable way. Here the methods of ab initio molecular quantum mechanics have been used to achieve this goal. Current methods of molecular quantum chemistry have been used for electronic structure calculations to study chemical reactions, photochemistry, interactions, structure and spectroscopy. However, attempts to apply these paradigms to biological systems or condensed phase chemistry have been limited due to their excessive computational costs. In these types of systems it has therefore become common to divide them into regions, and only solve the electronic problem for the region of primary interest and incorporate the effects of the rest of the system by means of operators that modify the electronic Hamiltonian of the region of primary interest [7]. Even then the electronic structure has to be modeled by very simple systems that may be unable to properly represent the phenomena being studied. More complex models can be used by recognizing that most such systems can be constructed from chemically identifiable substructures. This feature of larger molecular aggregates can be incorporated into the wave function by writing it as an antisymmetrized product of group functions. Such a theory was proposed earlier and was based on relaxing the orthogonality constraints between the group functions and reformulating the variational problem for determining the orbitals in terms of suitably defined group energies 183. Applying this method of nonorthogonal group functions (NOGF) to various systems has demonstrated its reliability and its potential for reducing computational cost 19,101.
Physics 204 (19961403-4/O
The purpose of the present analysis is to investigate the effects of polar solvent molecules, which can hydrogen bond to the groups of interest, on the stability of isolated charged residues and the relative stability of acid-base pairs in the charged and neutral states. Since the aim of these studies is to apply the methodology and findings to large biological systems, we have used the NOGF approximation to carry out the calculations. The functional side chains of the amino acids aspartate, glutamate, lysine, arginine and histidine will be considered here, and they are modeled by formic acid, ammonia, guanidine (GD) and imidazole (IMI), respectively. The environment is modeled by two water molecules hydrogen bonded to each of these fragments. In spite of the simplicity of these models, they appear to be capable of exhibiting the most important features affecting the interactions of ion pairs.
2. Methods The general formulation of the NOGF approximation for closed shell systems has been described in detail elsewhere [8,9,11], and applications of the method to a variety of molecular aggregates was recently reviewed [lo]. Here we consider the influence of the part of the first solvation shell, comprising two water molecules, on the protonation energies of HCOOH, NH,, GD and IMI. When these systems interact with the water molecules the corresponding NOGF wave functions are of the form *x =~‘[@x%,0%,0]
9
(1)
where X represents any one of the ionizable species in the protonated or deprotonated state, and d’ is an antisymmetrizer which only exchanges electrons between groups. Each group function, QR, is an antisymmetrized product of NR doubly occupied orthonormal orbitals, but no orthogonality constraints are assumed between the orbitals belonging to different groups [IO]. The interactions of the isolated acid-base pairs are represented by two group wave functions TPJJa=-pP’[@“*@oB], and
when
the pair
(2) is surrounded
by four
water
M.P. Fiilscher, E.L. Mehler/Chemical
b_u )_H_P_.%:::
p--i”,-. H-C
H-C 'JO ,' H
\
H
0'\ charged
H neulral
pair
fi H
,I"-II ,O
H’
pair
H '0'
H
-N’
'H
B
H Heo’H
Fig. 1. Schematic representations of H-bonding patterns of complexes in the neutral and charged states. (A) Singly H-bonded NH,-HOOCH and IMI-HOOCH; Note that for IMI one water is H-bonded to the H of N,, the other to the H of C, and N, participates in H-bonding with formic acid. (B) Doubly H-bonded GD-HOOCH. R denotes proton donor-proton acceptor distance for calculation of the interaction potentials. For geometries see references given in the text.
molecules in the first solvation tion has the form T*a=.&[@
A @ B @ H,O
shell, the wave func-
@ H,O @ H,O @ H,O 1.
(3)
Because of the relaxation of the orthogonality constraints between groups, the basis expansions for the orbitals of a group need only span a part of the Hilbert space defining the complete system. The systems described by the wave functions of Eqs. (2) and (31, consist of 2 and 6 groups, respectively, that only interact non-covalently, as shown in Fig. 1. In such cases it seems most natural to expand the orbitals of each groups in terms of the basis set used to expand the orbitals of the corresponding isolated molecule. However, it has been shown that such orbitals cannot directly account for charge transfer effects, but the modifications of the orbital expansions required to account for these effects has been discussed elsewhere 19,111. For the hydrogen-bonding interactions considered here, it is only necessary
Physics 204 (1996) 403-410
405
to extend the orbital expansion of the proton acceptor by the basis functions of the proton that defines the H-bond. By allowing the basis space of the proton to be shared by both proton acceptor and proton donor, charge transfer effects are completely accounted for, and, in addition, this procedure has been shown to minimize basis set superposition errors [ 10,111. In the NOGF approximation the orbitals of different groups can be represented at different levels of approximation. The form of the orbital equations ensures that no variational distortion results from using such mixed basis sets and that reliable results are achieved with substantial reduction in computing costs [9]. Since the water molecules are part of the environment, they do not participate directly in the interactions of interest, but modify them. For such environmental groups it is reasonable to expand the orbitals in smaller basis sets than the bases used for the orbitals of the species involved in the interactions of interest. Therefore, for formic acid and the bases, double zeta (DZ) basis sets are constructed from the (7~3~) atomic bases of Roos and Siegbahn [12] and Dunning’s (4s) basis for hydrogen [ 131, but the orbitals of the water molecules are expanded in a minimum basis (MB) contracted from (5~2~) atomic bases due to Graf and Mehler [14] and Huzinaga’s (3s) basis for hydrogen [15]. To compute intermolecular interaction potentials of the isolated acid-base pairs, the monomer geometries were kept frozen throughout. The structure of the monomers, and relative orientations and the positions of the H-bonded water molecules were determined previously [9,10,16]. To compute the wave functions for the acid-base pairs interacting with the water molecules transferability of the intermolecular separations and orientations is assumed. The conformational arrangement of the acid-base pairs is given schematically in Fig. 1. The interaction of ammonia and imidazole with formate comprises a single H-bond (Fig. lA), whereas for guanidine a double H-bonded configuration was used (Fig. 1B) since this is the structure of lowest energy [17]. The protonation energies of the isolated and the solvated monomers were calculated, and the optimum acidbase separation was determined for the isolated pairs in the charged and neutral states; these distances were then also used for the single point calculations
406
M.P. Fiilscher, E.L. Mehier / Chemical Physics 204 (1996) 403-410
on the solvated pairs. All computations were performed with the program described in Ref. [9].
respectively. From the reaction cycle the protonation energy, PES, can be written in the alternative form PES=PE+SxH+-S,.
3. Results and discussion 3.1. Protonation
Eq. (6) shows that in a polar solvent the electrostatic interactions will always stabilize XH+ relative to X since the charge-dipole interaction will always be stronger than the interaction of a dipole with a dipole or higher order pole. A way to calculate the full effect of bulk solvent is by the methods of macroscopic continuum electrostatics [22,23] which, of course, also shows that the high dielectric environment favors the charged state. Eq. (6) shows that this result is inherent in the relation between PES and PE, and does not depend on the method used to calculate these quantities. The PE of the isolated species are given in Table 1. Comparison of the calculated PEs with the observed values shows that they are in error by 22-26 kcal/mol. The fact that these errors are quite close is fortunate as will become clear below. The remaining columns in Table 1 show how the solvation energies contribute to the stabilization of the charged species, which results in the protonation energy of formate increasing whereas it decreases for the bases. Perhaps most surprising is the magnitude of the effect given the fact that only two waters are bound to each species, although the magnitude may be somewhat overemphasized because the screening of the electrostatic interactions due to the remaining system is neglected. These waters comprise part of the nearest solvation shell, and the results simply reflect the fact that the first solvation shell will have by far the strongest modulating effect on the stability of any species. Moreover, in our calculations the full polarizing and charge transfer effects of the solvent molecules on the protonated and deprotonated species are included.
energies
The calculation of the protonation energy is of fundamental importance for the determination of acid-base equilibria [ 18-211. Moreover, this quantity also plays a key role in the way solvent (or environment) modulates ion pair formation, as will be seen below. The protonation energy for the reaction X+H+=XH+ is defined by PE= E,,+-
E,,
(da)
with E,,+ and E, being the total energies of the protonated and deprotonated systems XH+ and X, respectively. For the solvated system, S . X, the protonation energy for the reaction S.X+H+oS.XH+ can be similarly
defined as
PES = Es.xH+-
E,.,.
(4b)
The relationship between PE and PES can be obtained from the following closed reaction cycle: S
S.X+H+
&
s.XH+ S XH+ II S+XH+
Sx
II S+X+H+
pE
(5)
where S, and S,,+ are the interaction (solvation) energies between solvent and species X and XH+,
Table I Decomposition
(6)
of the protonation
energy due to salvation
Molecule
-PE
- exp(err) b
-sX
-SXH+
-PES
6PE ’
HCOONH, Guanidine Imidazole
368.7 227.7 263.2 248.6
342.2 ( - 26.5) 203.6(-24.1)
36.0 5.5 7.1 9.6
8.7 49.7 35.7 34.4
341.5 272.0 291.8 273.5
27.3 -44.3 - 28.6 - 24.9
a Energies in kcal/mol.
b Ref. [25].
a
226.4 (- 22.2) ’ 6PE = PES - PE.
M.P. Fiilscher, E.L. Mehler/Chemical
3.2. Environmental
modulation
407
Physics 204 (1996) 403-410
of ion pair fonnation
In proteins ion pairs are surrounded by amino acid residues and, depending on the extent of solvent exposure, water. All these groups will interact with the ion pair, modulating the relative stability of the charged and neutral states and the proton transfer potential, and depending on the proximity of these molecules to the ion pair, the non-electrostatic interactions will be of increasing importance. Ion pair formation is described by the proton transfer reaction AH-B
4
Fig. 2. Interaction potential of ammonia-formic curve: NH: -_ OOCH; light curve NH,-HOOCH.
* A--BH+,
acid.
Heavy
and the affinity to adopt the ion pair state is defined by AE,,=
EA--BH+-
EAH-B.
(7)
A decomposition of A&,, which is just the energy required to change the neutral pair (NP) into the charged pair (CP), can be obtained from the reaction cycle AH-B
I
<
A Ear,
> A-HB+
II
AEA-B+
AEAt3
AH+B
M
(8)
A- + BH+
APE = PE, - PE,and allows Eq. (7) to be written form AE,,=
AE,-,++ = AE,-,+(net)
APE-
in the alternative
AEAB
- AE,,.
are on the same energy scale and therefore are also the appropriate quantities to be used for calculating proton transfer potentials [ 101. For each isolated, interacting pair, calculations were carried out using the group wave function structure of Eq. (2) at several intermolecular separations (see Fig. 1). The intermolecular potentials for ammonia-formic acid (Fig. 2) and IMI-formic acid (Fig. 3) exhibit several similar features: In the gas phase the NP is lower in energy than the CP for all intermolecular distances, and the separation between the energy minima of the NP and CP is comparable in magnitude, e.g., 15.6 kcal/mol for the ammonia-formic acid and 17.0 kcal/mol for the IMI-formic acid pairs. In contrast, the GD-formic acid interaction potentials behave quite differently as shown in Fig. 4. At
(9)
APE is the difference in protonation energy of the base and the acid, and can be considered to represent the work required to form the charges BHf and Afrom the neutral systems, when they are at infinite separation. In the supermolecule approach this work is automatically removed when calculating A EA B+ from the difference in energy between the charged species in the bound and unbound states. However, from Eq. (8) it is seen that the quantity APE is a required part of the energy of ion-pair formation and therefore A EA-B+ (net) represents the proper value of the interaction energy of the charged pair. It should be further noted that A E,- B+ (net> and A EAB
-10
1
22
i 24
26
3
Fig. 3. Interaction potential of imidazole-formic curve: IMIH+ -_ OOCH; light curve IMI-HOOCH.
3.2
acid.
3.4
Heavy
M.P. Fiilscher,
E.L. Mehler/Chemical I
Physics 204 (1996)
Table 2 Electrostatic and total interaction energies of the guanidine-formic acid pair in the neutral and charged states a R(C-C)
Neutral pair
(.Q
ES
A Ebb
AES(net)
AE,-,+(net)
4.00 4.50
-31.6 - 16.2
- 26.6 - 20.4
- 20.6 6.2
- 29.0 -9.1
a Energies
3’4
3’6
318
R&C)
412
4:4
4:6
4’6
A
Fig. 4. Interaction potential of guanidine-formic curve: GDH+-OOCH; light curve GD-HOOCH.
acid.
Heavy
intermolecular distances (R in Fig. 1B) greater than about 4.1 A the NP is lower in energy than the CP, but at 4.1 A the two curves cross and for shorter distances the CP is more stable than the NP. The vertical separation of the minima is 4.3 kcal/mol. This effect is due to the resonance stabilization of guanidinium resulting from the delocalization of the C-N double bond that also accounts for the very high pK of guanidinium. For the GD-formic acid pair the additional stabilization resulting from the resonance effect is clearly not electrostatic in origin as mty be seen from Table 2. For a C-C separation of 4 A the CP is about 3.5 kcal/mol more stable than the NP, but the electrostatic component of the interaction energy would predict that the NP is more stable than the CP. At a
Table 3 Acid-base
ion pair formation
and effects of solvation
System
A EAB
NH,-HOOCH solvated NH:--0OCH solvated GD-HOOCH solvated GDH+OOCH solvated IMI-HOOCH solvatcd IMIH+- - OOCH solvated
- 14.7 - 25.2
* Energies in kcal/mol.
403-410
Ion pair
in kcal/mol.
distance of 4.5 A both electrostatic and total energy predict the NP to be more stable. However, the electrostatic component gives a difference of about 22 kcal/mol, while the difference in total interaction energy between the two species is only about 11 kcal/mol. These results show that electrostatic interactions alone, do not always predict even qualitatively correct results when interactions of an essentially quantum mechanical origin make important contribution to the interaction energy. In such cases it is essential that all components of the interaction be included in the calculation. If the interacting species are solvated a reaction cycle similar to Eq. (8) can be formulated [ 181. The affinity towards the solvated ion pair then reads AE& = AE&++
APES - AE,
= AEi-,+(net)
- AEA.
(10)
For each solvated pair (see Fig. l), the wave function given in Eq. (3) was used to calculate the intermolecular energies of the CP and NP at the values of R,,,
a
AEAma+
APE
A EA -B +(net)
- 140.2 - 93.6
141.0 69.5
0.8 - 24.2
15.5 0.9
- 133.6 -99.1
102.7 49.7
- 30.9 - 49.5
- 4.2 - 19.8
- 116.4 - 93.0
120.1 68.0
3.7 - 25.0
17.1 - 3.4
A E,rr
- 26.7 - 29.7
- 13.4 -21.6
M.P. Fiilscher, E.L. Mehler/Chemical
obtained from the potential curves of the corresponding isolated pairs. The resulting Al& and the decomposition defined by Eq. (10) are given in Table 3. The analysis for the isolated pairs is also included. Comparing the values of AE, with A E: and shows that the effect of AEA-B+ with hEi-,+ solvation is different for the two states. The interaction energy of the NP decreases upon solvation whereas for the CP it increases. This behavior is primarily due to the fact that the increase in the exchange energy upon solvation is much larger for the charged than the neutral pairs, although the electrostatic contribution also increases slightly in the CP whereas for the NP it decreases. In any case the expected increase in stability ofcharged relative to neutral pairs is clearly not due to the effect of polar solvent molecules on the electrostatic interactions between the pairs. The values for APE can be calculated from the protonation energies given in Table 1, and it is seen that the effect of solvation on these quantities is dramatic. In each case the vacuum value is approximately halved due to the interaction of the ion pair with only four water molecules. From the values of the shifts in PE due to solvation (6PE in Table l), which are of opposite sign for formate and the bases, and the fact that APE is the difference between the PE, the 6PE values add, resulting in the very large decrease in APE due to the interaction with water. This effect is clearly dominant as evidenced by the values of the net interaction energies of the CPs that are 20-30 kcal/mol lower in the solvated pair than in the isolated pair, in spite of the fact that A Ei- B+ is about 40 kcal/mol higher than AE,-,+. The net result is that the affinity to form the CP increases due to the solvation, and in the present calculations amounts to 4-5 kcal/mol-solvent molecules in the first solvation shell that are H-bonded to the members of the interacting pairs. The enhanced stability of the CP relative to the NP is sufficient to stabilize solvated IMIH+-formate relative to IMI-formic acid whereas for the ammonia-formic acid the CP is about 1 kcal/mol less stable than the NP as compared to 15 kcal/mol for the isolated system. The GDH+--0OCH pair, which in the region of optimal separation is already stable by about 4 kcal/mol, is further stabilized by the water molecules to nearly 20 kcal/mol.
Physics 204 (1996) 403-410
409
4. Conclusions The results of the calculations presented above show that proximal polar groups have a large stabilizing effect on isolated charges, and they increase the stability of interacting pairs in both the neutral and charged state. The effect on the latter is always larger than on the former and therefore tends to shift the equilibrium towards the CP. However, the origin of this effect is not due to interactions of the polar molecules with the interacting pair, which tend to destabilize the CP relative to the NP, but is due to the large changes in the energy of charge formation caused by the presence of the polar groups. Several caveats should be kept in mind when assessing these calculations: In computing the interaction potentials of the pairs, the position of the protons participating in the H-bonding were not reoptimized at each point of the potential curves. The consequence of such optimization would be to shift each point to a slightly lower energy and possibly yield small shifts of the minima. It seems unlikely, however, that any of the qualitative features of the interaction potentials or their consequences on the stabilization of ion pairs or isolated charged groups would substantially change. Other weaknesses of the present calculations are the small basis sets used and the neglect of the electron correlation contribution to the total energies. Fortunately, for many applications of the type discussed here, correlation effects can be neglected in the first approximation. On the other hand, extended basis sets including diffuse functions would improve, in particular, the protonation energies. However, as it is only the differences in the protonation energies that enters the analysis presented above, and in the present calculations the error in the protonation energies are similar for all the systems considered, the effect of improving the basis sets is expected to be small for the calculations discussed here. An intriguing result of the present calculations is the fundamental difference in behavior of the interaction potentials of GD-HCOOH and the other two acid-base pairs. The fact that the curve for the CP, GDH+--OOCH, crosses the potential curve of the NP suggests that even in a low dielectric environment this pair could be in the charged state, while both NH,-HOOCH and IMI-HCOOH would be in
410
M.P. Fiilscher, E.L. Mrhler/
Chemical Physics 204 (1996) 403-410
the neutral state. The resonance effect would be somewhat smaller in arginine, but its pK of 12.7 suggests that it is still operative. Thus the result found for GDH+--0OCH appears to be applicable to the interaction of Arg with acidic residues in proteins and suggests that this residue might assume functions not just attributable to its charge. This possibility is supported by the finding that the charge conserving mutation Arg to Lys can affect activity [24]. Also of interest is the fact that the C-C distance at which the two curves cross (see Fig. 4) is typical for the separation between Arg and either Asp or Glu in many proteins. Given the dynamical nature of protein structure and the structural changes that often occur during reaction, it does not seem unreasonable to speculate about the functional role of the curve crossing. If one considers a situation where an Arg-acid residue pair are approaching each other in a low dielectric environment there will be an abrupt change in the gradient of the potential at the distance where the charged state becomes more stable. Moreover, the change in the forces will tend to bring the pair still closer together. Such a mechanism could easily be tested by mutating Arg to Lys, which is not expected to exhibit such a curve crossing (see Fig. 2). Thus for enzymes where this property of Arg is important for function, this mutation would be expected to lead to lower activity in the resulting mutant.
Acknowledgement The research reported in this paper was carried out at the Biocenter of the University of Basle, Switzerland, and was supported by Swiss National Science Foundation Grants 3 l-8840.86 and 3 l-26261.89 and the ‘Stipendium des Kantons Basel-Land’ (no. 11343). The support of University of Base1 Computing Center and the Computing Center BaselCity for generous grants of computing time is gratefully acknowledged. Finally, the partial support of the Swedish Natural Science Research Council
(MPF), and the National Institutes (ELM) is gratefully acknowledged.
of Health,
USA
References []I A.A. Rashin, B. Honig, J. Mol. Biol. 173 (1984) 515. 121F.A. Quiocho, J.S. Sack and N.K. Vyas, Nature 329 (1987) 561. [31 J.N. Jansonius and M.G. Vincent, in: Biological macromolecules and assemblies, Vol. 3, eds. F.A. Jumak and A. McPherson (Wiley, New York, 1987) pp. 187-285. [41 J.F. Kirsch, G. Eichele, G.C. Ford, M.G. Vincent, J.N. Jansonius, H. Gehring and P. Christen, J. Mol. Biol. 174 ( 1984) 497. [51 A.E. Mark and W.F. van Gunsteren, J. Mol. Biol. 240 (1994) 167. [61 S. Yun-yu, A.E. Mark, W. Cun-xin, H. Fuhua, H.J.C., B. and W.F. van Gunsteren, Prot. Eng. 6 (1993) 289. 171 M.L.J. Drummond, Prog. Biophys. Mol. Biol. 47 (1986) I. [81 E.L. Mehler, J. Chem. Phys. 67 (1977) 2728. [91 M.P. Wlscher and E.L. Mehler, J. Comput. Chem. 12 (1991) 811. [lOI EL. Mehler, J. Math. Chem. 10 (1992) 57. 1111 EL. Mehler, J. Chem. Phys. 74 (1981) 6298. [121 B. Roos and P. Siegbahn, Theoret. Chim. Acta 17 (1970) 209. [I31 T.H. Dunning, J. Chem. Phys. 53 (1970) 2823. 1141 P. Graf and E.L. Mehler, Intern. J. Quantum Chem. QBS8 (1981) 49. [I51 S. Huzinaga, J. Chem. Phys. 42 (1965) 1293. [I61 M.P. Wlscher and E.L. Mehler, J. Mol. Strut. THEOCHEM 165 (1988) 319. [171 A.M. Sapse and C.S. Russel, Intern. J. Quantum Chem. (1984) 91. [IS] D. Bashford and M. Karplus, Biochemistry 29 (1990) 10219. [ 191 D. Bashford and M. Karplus, J. Phys. Chem. 95 (1991) 9556. 1201 C. Lim, D. Bashford and M. Karplus, J. Phys. Chem. 95 (1991) 5610. [21] J. Antosiewics, J.A. McCammon and M.K. Gilson, J. Mol. Biol. 238 (1994) 415. 1221 M. Born, Z. Physik 1 (1920) 45. [231 K.A. Sharp and B. Honig, Amm. Rev. Biophys. Biophys. Chem. 19 (1990) 301. [241 H. Friesen and P.D. Sadowski, J. Mol. Biol. 225 (1992) 313. [251 HCOOH: T.B. McMahon and P. Kabarle, J. Am. Chem. Sot. 99 (1977) 2222; NH,: S.T. Ceyer, P.W. Tiedemann and B.H. Mahan, J. Chem. Phys. 70 (1979) 14; Met-Imidazole: M.R. Ellenberg, R.A. Eades, M.W. Thomsen, WE. Fameth and D.A. Dixon, J. Am. Chem. Sot. 101 (1979) 7151.
Chemical Physics 204 (1996) 403-410
The effect of hydrogen bonding on protonation energies and ion pair formation Markus P. Filscher a Depurtment ofTheoretical b Depurtment
ofPhysiology
Chemistry,
University
and Biophysics,
a, Ernest L. Mehler b,* ofLund,
Chemical
Mount Sinai School
Centre, P.O.B. 124, 2210 Lund, Sweden
of Medicine,
CUNY, New York, NY 10029,
USA
Received February 21, 1995
Abstract The effects of environmental hydrogen bonding by solvent molecules on the stability of charged groups and ion pairs have been analyzed by ab initio calculations on molecular models representing acidic and basic amino acid residues. By constructing closed reactions cycles, the energetics of the isolated and solvated groups can be directly compared and the origin of the changes in relative stability exhibited. The results of the calculations showed that the solvent molecules stabilize interacting pairs in both the neutral and charged states, but that the effects shift the equilibrium towards the latter. The analysis further indicated that the main contribution to this effect is due to the changes in the work required to form the charges attendant upon solvation. Moreover, although the relative stability of the neutral and charged states is primarily controlled by electrostatic effects, it can be modulated by additional short range interactions when the species are close to each other. In such cases the relative stability of the states predicted by the total interaction energy may be different to what
is obtain from just the electrostatic interactions.
1. Introduction The stability of ion pairs is of fundamental importance in protein structure and function. The majority of ionizable groups are found at or near the protein surface where they can form dynamic networks of hydrogen bonds involving solvent molecules and possibly proximal polar residues. However, some titratable groups are found in the protein interior, mostly in pairs, but sometimes as isolated charges [ 1,2], and they are often involved in stabilizing the active site region or participate directly in the catalytic mechanism [3,4]. To explain their effects in
??
Corresponding author.
0301-0104/96/$15.00
0 1996 Elsevier Science B.V. All rights resewed
SSDI 0301-0104(95)00423-S
terms of the attendant free energy changes at the molecular level, and to evaluate them quantitatively by computational methods are presently major challenges to theory. Molecular simulation techniques are one way to obtain free energies, but it is necessary to include a large part of the environment in the model for such calculations to approach an acceptable degree of reliability. These methods also have the disadvantage that it becomes very difficult to identify the individual interactions which lead to the observed macroscopic effects or to assess their significance [5,61. In spite of the obvious limitations, analysis of the interactions leading to energy changes seems to be most meaningful when carried out on models that are not dependent on experimentally determined parameters, and are flexible enough to
404
M.P. Fiilscher, E.L. Mehler/Chemical
account for the relevant contributions to the interaction energy. Due to the long range nature of Coulombic forces, the intermolecular interactions which maintain protein structure are dominated by electrostatics. At shorter intermolecular separations, and especially those where chemical processes take place, additional effects of quantum chemical origin become important and may even be dominant. These include charge transfer, exchange and correlation effects that sometimes lead to unexpected results not predictable by a consideration of electrostatic effects alone. Therefore, to provide a proper theoretical accounting of protein structure and function, it is essential that these interactions are described in a reasonably reliable way. Here the methods of ab initio molecular quantum mechanics have been used to achieve this goal. Current methods of molecular quantum chemistry have been used for electronic structure calculations to study chemical reactions, photochemistry, interactions, structure and spectroscopy. However, attempts to apply these paradigms to biological systems or condensed phase chemistry have been limited due to their excessive computational costs. In these types of systems it has therefore become common to divide them into regions, and only solve the electronic problem for the region of primary interest and incorporate the effects of the rest of the system by means of operators that modify the electronic Hamiltonian of the region of primary interest [7]. Even then the electronic structure has to be modeled by very simple systems that may be unable to properly represent the phenomena being studied. More complex models can be used by recognizing that most such systems can be constructed from chemically identifiable substructures. This feature of larger molecular aggregates can be incorporated into the wave function by writing it as an antisymmetrized product of group functions. Such a theory was proposed earlier and was based on relaxing the orthogonality constraints between the group functions and reformulating the variational problem for determining the orbitals in terms of suitably defined group energies 183. Applying this method of nonorthogonal group functions (NOGF) to various systems has demonstrated its reliability and its potential for reducing computational cost 19,101.
Physics 204 (19961403-4/O
The purpose of the present analysis is to investigate the effects of polar solvent molecules, which can hydrogen bond to the groups of interest, on the stability of isolated charged residues and the relative stability of acid-base pairs in the charged and neutral states. Since the aim of these studies is to apply the methodology and findings to large biological systems, we have used the NOGF approximation to carry out the calculations. The functional side chains of the amino acids aspartate, glutamate, lysine, arginine and histidine will be considered here, and they are modeled by formic acid, ammonia, guanidine (GD) and imidazole (IMI), respectively. The environment is modeled by two water molecules hydrogen bonded to each of these fragments. In spite of the simplicity of these models, they appear to be capable of exhibiting the most important features affecting the interactions of ion pairs.
2. Methods The general formulation of the NOGF approximation for closed shell systems has been described in detail elsewhere [8,9,11], and applications of the method to a variety of molecular aggregates was recently reviewed [lo]. Here we consider the influence of the part of the first solvation shell, comprising two water molecules, on the protonation energies of HCOOH, NH,, GD and IMI. When these systems interact with the water molecules the corresponding NOGF wave functions are of the form *x =~‘[@x%,0%,0]
9
(1)
where X represents any one of the ionizable species in the protonated or deprotonated state, and d’ is an antisymmetrizer which only exchanges electrons between groups. Each group function, QR, is an antisymmetrized product of NR doubly occupied orthonormal orbitals, but no orthogonality constraints are assumed between the orbitals belonging to different groups [IO]. The interactions of the isolated acid-base pairs are represented by two group wave functions TPJJa=-pP’[@“*@oB], and
when
the pair
(2) is surrounded
by four
water
M.P. Fiilscher, E.L. Mehler/Chemical
b_u )_H_P_.%:::
p--i”,-. H-C
H-C 'JO ,' H
\
H
0'\ charged
H neulral
pair
fi H
,I"-II ,O
H’
pair
H '0'
H
-N’
'H
B
H Heo’H
Fig. 1. Schematic representations of H-bonding patterns of complexes in the neutral and charged states. (A) Singly H-bonded NH,-HOOCH and IMI-HOOCH; Note that for IMI one water is H-bonded to the H of N,, the other to the H of C, and N, participates in H-bonding with formic acid. (B) Doubly H-bonded GD-HOOCH. R denotes proton donor-proton acceptor distance for calculation of the interaction potentials. For geometries see references given in the text.
molecules in the first solvation tion has the form T*a=.&[@
A @ B @ H,O
shell, the wave func-
@ H,O @ H,O @ H,O 1.
(3)
Because of the relaxation of the orthogonality constraints between groups, the basis expansions for the orbitals of a group need only span a part of the Hilbert space defining the complete system. The systems described by the wave functions of Eqs. (2) and (31, consist of 2 and 6 groups, respectively, that only interact non-covalently, as shown in Fig. 1. In such cases it seems most natural to expand the orbitals of each groups in terms of the basis set used to expand the orbitals of the corresponding isolated molecule. However, it has been shown that such orbitals cannot directly account for charge transfer effects, but the modifications of the orbital expansions required to account for these effects has been discussed elsewhere 19,111. For the hydrogen-bonding interactions considered here, it is only necessary
Physics 204 (1996) 403-410
405
to extend the orbital expansion of the proton acceptor by the basis functions of the proton that defines the H-bond. By allowing the basis space of the proton to be shared by both proton acceptor and proton donor, charge transfer effects are completely accounted for, and, in addition, this procedure has been shown to minimize basis set superposition errors [ 10,111. In the NOGF approximation the orbitals of different groups can be represented at different levels of approximation. The form of the orbital equations ensures that no variational distortion results from using such mixed basis sets and that reliable results are achieved with substantial reduction in computing costs [9]. Since the water molecules are part of the environment, they do not participate directly in the interactions of interest, but modify them. For such environmental groups it is reasonable to expand the orbitals in smaller basis sets than the bases used for the orbitals of the species involved in the interactions of interest. Therefore, for formic acid and the bases, double zeta (DZ) basis sets are constructed from the (7~3~) atomic bases of Roos and Siegbahn [12] and Dunning’s (4s) basis for hydrogen [ 131, but the orbitals of the water molecules are expanded in a minimum basis (MB) contracted from (5~2~) atomic bases due to Graf and Mehler [14] and Huzinaga’s (3s) basis for hydrogen [15]. To compute intermolecular interaction potentials of the isolated acid-base pairs, the monomer geometries were kept frozen throughout. The structure of the monomers, and relative orientations and the positions of the H-bonded water molecules were determined previously [9,10,16]. To compute the wave functions for the acid-base pairs interacting with the water molecules transferability of the intermolecular separations and orientations is assumed. The conformational arrangement of the acid-base pairs is given schematically in Fig. 1. The interaction of ammonia and imidazole with formate comprises a single H-bond (Fig. lA), whereas for guanidine a double H-bonded configuration was used (Fig. 1B) since this is the structure of lowest energy [17]. The protonation energies of the isolated and the solvated monomers were calculated, and the optimum acidbase separation was determined for the isolated pairs in the charged and neutral states; these distances were then also used for the single point calculations
406
M.P. Fiilscher, E.L. Mehier / Chemical Physics 204 (1996) 403-410
on the solvated pairs. All computations were performed with the program described in Ref. [9].
respectively. From the reaction cycle the protonation energy, PES, can be written in the alternative form PES=PE+SxH+-S,.
3. Results and discussion 3.1. Protonation
Eq. (6) shows that in a polar solvent the electrostatic interactions will always stabilize XH+ relative to X since the charge-dipole interaction will always be stronger than the interaction of a dipole with a dipole or higher order pole. A way to calculate the full effect of bulk solvent is by the methods of macroscopic continuum electrostatics [22,23] which, of course, also shows that the high dielectric environment favors the charged state. Eq. (6) shows that this result is inherent in the relation between PES and PE, and does not depend on the method used to calculate these quantities. The PE of the isolated species are given in Table 1. Comparison of the calculated PEs with the observed values shows that they are in error by 22-26 kcal/mol. The fact that these errors are quite close is fortunate as will become clear below. The remaining columns in Table 1 show how the solvation energies contribute to the stabilization of the charged species, which results in the protonation energy of formate increasing whereas it decreases for the bases. Perhaps most surprising is the magnitude of the effect given the fact that only two waters are bound to each species, although the magnitude may be somewhat overemphasized because the screening of the electrostatic interactions due to the remaining system is neglected. These waters comprise part of the nearest solvation shell, and the results simply reflect the fact that the first solvation shell will have by far the strongest modulating effect on the stability of any species. Moreover, in our calculations the full polarizing and charge transfer effects of the solvent molecules on the protonated and deprotonated species are included.
energies
The calculation of the protonation energy is of fundamental importance for the determination of acid-base equilibria [ 18-211. Moreover, this quantity also plays a key role in the way solvent (or environment) modulates ion pair formation, as will be seen below. The protonation energy for the reaction X+H+=XH+ is defined by PE= E,,+-
E,,
(da)
with E,,+ and E, being the total energies of the protonated and deprotonated systems XH+ and X, respectively. For the solvated system, S . X, the protonation energy for the reaction S.X+H+oS.XH+ can be similarly
defined as
PES = Es.xH+-
E,.,.
(4b)
The relationship between PE and PES can be obtained from the following closed reaction cycle: S
S.X+H+
&
s.XH+ S XH+ II S+XH+
Sx
II S+X+H+
pE
(5)
where S, and S,,+ are the interaction (solvation) energies between solvent and species X and XH+,
Table I Decomposition
(6)
of the protonation
energy due to salvation
Molecule
-PE
- exp(err) b
-sX
-SXH+
-PES
6PE ’
HCOONH, Guanidine Imidazole
368.7 227.7 263.2 248.6
342.2 ( - 26.5) 203.6(-24.1)
36.0 5.5 7.1 9.6
8.7 49.7 35.7 34.4
341.5 272.0 291.8 273.5
27.3 -44.3 - 28.6 - 24.9
a Energies in kcal/mol.
b Ref. [25].
a
226.4 (- 22.2) ’ 6PE = PES - PE.
M.P. Fiilscher, E.L. Mehler/Chemical
3.2. Environmental
modulation
407
Physics 204 (1996) 403-410
of ion pair fonnation
In proteins ion pairs are surrounded by amino acid residues and, depending on the extent of solvent exposure, water. All these groups will interact with the ion pair, modulating the relative stability of the charged and neutral states and the proton transfer potential, and depending on the proximity of these molecules to the ion pair, the non-electrostatic interactions will be of increasing importance. Ion pair formation is described by the proton transfer reaction AH-B
4
Fig. 2. Interaction potential of ammonia-formic curve: NH: -_ OOCH; light curve NH,-HOOCH.
* A--BH+,
acid.
Heavy
and the affinity to adopt the ion pair state is defined by AE,,=
EA--BH+-
EAH-B.
(7)
A decomposition of A&,, which is just the energy required to change the neutral pair (NP) into the charged pair (CP), can be obtained from the reaction cycle AH-B
I
<
A Ear,
> A-HB+
II
AEA-B+
AEAt3
AH+B
M
(8)
A- + BH+
APE = PE, - PE,and allows Eq. (7) to be written form AE,,=
AE,-,++ = AE,-,+(net)
APE-
in the alternative
AEAB
- AE,,.
are on the same energy scale and therefore are also the appropriate quantities to be used for calculating proton transfer potentials [ 101. For each isolated, interacting pair, calculations were carried out using the group wave function structure of Eq. (2) at several intermolecular separations (see Fig. 1). The intermolecular potentials for ammonia-formic acid (Fig. 2) and IMI-formic acid (Fig. 3) exhibit several similar features: In the gas phase the NP is lower in energy than the CP for all intermolecular distances, and the separation between the energy minima of the NP and CP is comparable in magnitude, e.g., 15.6 kcal/mol for the ammonia-formic acid and 17.0 kcal/mol for the IMI-formic acid pairs. In contrast, the GD-formic acid interaction potentials behave quite differently as shown in Fig. 4. At
(9)
APE is the difference in protonation energy of the base and the acid, and can be considered to represent the work required to form the charges BHf and Afrom the neutral systems, when they are at infinite separation. In the supermolecule approach this work is automatically removed when calculating A EA B+ from the difference in energy between the charged species in the bound and unbound states. However, from Eq. (8) it is seen that the quantity APE is a required part of the energy of ion-pair formation and therefore A EA-B+ (net) represents the proper value of the interaction energy of the charged pair. It should be further noted that A E,- B+ (net> and A EAB
-10
1
22
i 24
26
3
Fig. 3. Interaction potential of imidazole-formic curve: IMIH+ -_ OOCH; light curve IMI-HOOCH.
3.2
acid.
3.4
Heavy
M.P. Fiilscher,
E.L. Mehler/Chemical I
Physics 204 (1996)
Table 2 Electrostatic and total interaction energies of the guanidine-formic acid pair in the neutral and charged states a R(C-C)
Neutral pair
(.Q
ES
A Ebb
AES(net)
AE,-,+(net)
4.00 4.50
-31.6 - 16.2
- 26.6 - 20.4
- 20.6 6.2
- 29.0 -9.1
a Energies
3’4
3’6
318
R&C)
412
4:4
4:6
4’6
A
Fig. 4. Interaction potential of guanidine-formic curve: GDH+-OOCH; light curve GD-HOOCH.
acid.
Heavy
intermolecular distances (R in Fig. 1B) greater than about 4.1 A the NP is lower in energy than the CP, but at 4.1 A the two curves cross and for shorter distances the CP is more stable than the NP. The vertical separation of the minima is 4.3 kcal/mol. This effect is due to the resonance stabilization of guanidinium resulting from the delocalization of the C-N double bond that also accounts for the very high pK of guanidinium. For the GD-formic acid pair the additional stabilization resulting from the resonance effect is clearly not electrostatic in origin as mty be seen from Table 2. For a C-C separation of 4 A the CP is about 3.5 kcal/mol more stable than the NP, but the electrostatic component of the interaction energy would predict that the NP is more stable than the CP. At a
Table 3 Acid-base
ion pair formation
and effects of solvation
System
A EAB
NH,-HOOCH solvated NH:--0OCH solvated GD-HOOCH solvated GDH+OOCH solvated IMI-HOOCH solvatcd IMIH+- - OOCH solvated
- 14.7 - 25.2
* Energies in kcal/mol.
403-410
Ion pair
in kcal/mol.
distance of 4.5 A both electrostatic and total energy predict the NP to be more stable. However, the electrostatic component gives a difference of about 22 kcal/mol, while the difference in total interaction energy between the two species is only about 11 kcal/mol. These results show that electrostatic interactions alone, do not always predict even qualitatively correct results when interactions of an essentially quantum mechanical origin make important contribution to the interaction energy. In such cases it is essential that all components of the interaction be included in the calculation. If the interacting species are solvated a reaction cycle similar to Eq. (8) can be formulated [ 181. The affinity towards the solvated ion pair then reads AE& = AE&++
APES - AE,
= AEi-,+(net)
- AEA.
(10)
For each solvated pair (see Fig. l), the wave function given in Eq. (3) was used to calculate the intermolecular energies of the CP and NP at the values of R,,,
a
AEAma+
APE
A EA -B +(net)
- 140.2 - 93.6
141.0 69.5
0.8 - 24.2
15.5 0.9
- 133.6 -99.1
102.7 49.7
- 30.9 - 49.5
- 4.2 - 19.8
- 116.4 - 93.0
120.1 68.0
3.7 - 25.0
17.1 - 3.4
A E,rr
- 26.7 - 29.7
- 13.4 -21.6
M.P. Fiilscher, E.L. Mehler/Chemical
obtained from the potential curves of the corresponding isolated pairs. The resulting Al& and the decomposition defined by Eq. (10) are given in Table 3. The analysis for the isolated pairs is also included. Comparing the values of AE, with A E: and shows that the effect of AEA-B+ with hEi-,+ solvation is different for the two states. The interaction energy of the NP decreases upon solvation whereas for the CP it increases. This behavior is primarily due to the fact that the increase in the exchange energy upon solvation is much larger for the charged than the neutral pairs, although the electrostatic contribution also increases slightly in the CP whereas for the NP it decreases. In any case the expected increase in stability ofcharged relative to neutral pairs is clearly not due to the effect of polar solvent molecules on the electrostatic interactions between the pairs. The values for APE can be calculated from the protonation energies given in Table 1, and it is seen that the effect of solvation on these quantities is dramatic. In each case the vacuum value is approximately halved due to the interaction of the ion pair with only four water molecules. From the values of the shifts in PE due to solvation (6PE in Table l), which are of opposite sign for formate and the bases, and the fact that APE is the difference between the PE, the 6PE values add, resulting in the very large decrease in APE due to the interaction with water. This effect is clearly dominant as evidenced by the values of the net interaction energies of the CPs that are 20-30 kcal/mol lower in the solvated pair than in the isolated pair, in spite of the fact that A Ei- B+ is about 40 kcal/mol higher than AE,-,+. The net result is that the affinity to form the CP increases due to the solvation, and in the present calculations amounts to 4-5 kcal/mol-solvent molecules in the first solvation shell that are H-bonded to the members of the interacting pairs. The enhanced stability of the CP relative to the NP is sufficient to stabilize solvated IMIH+-formate relative to IMI-formic acid whereas for the ammonia-formic acid the CP is about 1 kcal/mol less stable than the NP as compared to 15 kcal/mol for the isolated system. The GDH+--0OCH pair, which in the region of optimal separation is already stable by about 4 kcal/mol, is further stabilized by the water molecules to nearly 20 kcal/mol.
Physics 204 (1996) 403-410
409
4. Conclusions The results of the calculations presented above show that proximal polar groups have a large stabilizing effect on isolated charges, and they increase the stability of interacting pairs in both the neutral and charged state. The effect on the latter is always larger than on the former and therefore tends to shift the equilibrium towards the CP. However, the origin of this effect is not due to interactions of the polar molecules with the interacting pair, which tend to destabilize the CP relative to the NP, but is due to the large changes in the energy of charge formation caused by the presence of the polar groups. Several caveats should be kept in mind when assessing these calculations: In computing the interaction potentials of the pairs, the position of the protons participating in the H-bonding were not reoptimized at each point of the potential curves. The consequence of such optimization would be to shift each point to a slightly lower energy and possibly yield small shifts of the minima. It seems unlikely, however, that any of the qualitative features of the interaction potentials or their consequences on the stabilization of ion pairs or isolated charged groups would substantially change. Other weaknesses of the present calculations are the small basis sets used and the neglect of the electron correlation contribution to the total energies. Fortunately, for many applications of the type discussed here, correlation effects can be neglected in the first approximation. On the other hand, extended basis sets including diffuse functions would improve, in particular, the protonation energies. However, as it is only the differences in the protonation energies that enters the analysis presented above, and in the present calculations the error in the protonation energies are similar for all the systems considered, the effect of improving the basis sets is expected to be small for the calculations discussed here. An intriguing result of the present calculations is the fundamental difference in behavior of the interaction potentials of GD-HCOOH and the other two acid-base pairs. The fact that the curve for the CP, GDH+--OOCH, crosses the potential curve of the NP suggests that even in a low dielectric environment this pair could be in the charged state, while both NH,-HOOCH and IMI-HCOOH would be in
410
M.P. Fiilscher, E.L. Mrhler/
Chemical Physics 204 (1996) 403-410
the neutral state. The resonance effect would be somewhat smaller in arginine, but its pK of 12.7 suggests that it is still operative. Thus the result found for GDH+--0OCH appears to be applicable to the interaction of Arg with acidic residues in proteins and suggests that this residue might assume functions not just attributable to its charge. This possibility is supported by the finding that the charge conserving mutation Arg to Lys can affect activity [24]. Also of interest is the fact that the C-C distance at which the two curves cross (see Fig. 4) is typical for the separation between Arg and either Asp or Glu in many proteins. Given the dynamical nature of protein structure and the structural changes that often occur during reaction, it does not seem unreasonable to speculate about the functional role of the curve crossing. If one considers a situation where an Arg-acid residue pair are approaching each other in a low dielectric environment there will be an abrupt change in the gradient of the potential at the distance where the charged state becomes more stable. Moreover, the change in the forces will tend to bring the pair still closer together. Such a mechanism could easily be tested by mutating Arg to Lys, which is not expected to exhibit such a curve crossing (see Fig. 2). Thus for enzymes where this property of Arg is important for function, this mutation would be expected to lead to lower activity in the resulting mutant.
Acknowledgement The research reported in this paper was carried out at the Biocenter of the University of Basle, Switzerland, and was supported by Swiss National Science Foundation Grants 3 l-8840.86 and 3 l-26261.89 and the ‘Stipendium des Kantons Basel-Land’ (no. 11343). The support of University of Base1 Computing Center and the Computing Center BaselCity for generous grants of computing time is gratefully acknowledged. Finally, the partial support of the Swedish Natural Science Research Council
(MPF), and the National Institutes (ELM) is gratefully acknowledged.
of Health,
USA
References []I A.A. Rashin, B. Honig, J. Mol. Biol. 173 (1984) 515. 121F.A. Quiocho, J.S. Sack and N.K. Vyas, Nature 329 (1987) 561. [31 J.N. Jansonius and M.G. Vincent, in: Biological macromolecules and assemblies, Vol. 3, eds. F.A. Jumak and A. McPherson (Wiley, New York, 1987) pp. 187-285. [41 J.F. Kirsch, G. Eichele, G.C. Ford, M.G. Vincent, J.N. Jansonius, H. Gehring and P. Christen, J. Mol. Biol. 174 ( 1984) 497. [51 A.E. Mark and W.F. van Gunsteren, J. Mol. Biol. 240 (1994) 167. [61 S. Yun-yu, A.E. Mark, W. Cun-xin, H. Fuhua, H.J.C., B. and W.F. van Gunsteren, Prot. Eng. 6 (1993) 289. 171 M.L.J. Drummond, Prog. Biophys. Mol. Biol. 47 (1986) I. [81 E.L. Mehler, J. Chem. Phys. 67 (1977) 2728. [91 M.P. Wlscher and E.L. Mehler, J. Comput. Chem. 12 (1991) 811. [lOI EL. Mehler, J. Math. Chem. 10 (1992) 57. 1111 EL. Mehler, J. Chem. Phys. 74 (1981) 6298. [121 B. Roos and P. Siegbahn, Theoret. Chim. Acta 17 (1970) 209. [I31 T.H. Dunning, J. Chem. Phys. 53 (1970) 2823. 1141 P. Graf and E.L. Mehler, Intern. J. Quantum Chem. QBS8 (1981) 49. [I51 S. Huzinaga, J. Chem. Phys. 42 (1965) 1293. [I61 M.P. Wlscher and E.L. Mehler, J. Mol. Strut. THEOCHEM 165 (1988) 319. [171 A.M. Sapse and C.S. Russel, Intern. J. Quantum Chem. (1984) 91. [IS] D. Bashford and M. Karplus, Biochemistry 29 (1990) 10219. [ 191 D. Bashford and M. Karplus, J. Phys. Chem. 95 (1991) 9556. 1201 C. Lim, D. Bashford and M. Karplus, J. Phys. Chem. 95 (1991) 5610. [21] J. Antosiewics, J.A. McCammon and M.K. Gilson, J. Mol. Biol. 238 (1994) 415. 1221 M. Born, Z. Physik 1 (1920) 45. [231 K.A. Sharp and B. Honig, Amm. Rev. Biophys. Biophys. Chem. 19 (1990) 301. [241 H. Friesen and P.D. Sadowski, J. Mol. Biol. 225 (1992) 313. [251 HCOOH: T.B. McMahon and P. Kabarle, J. Am. Chem. Sot. 99 (1977) 2222; NH,: S.T. Ceyer, P.W. Tiedemann and B.H. Mahan, J. Chem. Phys. 70 (1979) 14; Met-Imidazole: M.R. Ellenberg, R.A. Eades, M.W. Thomsen, WE. Fameth and D.A. Dixon, J. Am. Chem. Sot. 101 (1979) 7151.