Modeling Amino Acid Side Chains. 3. Influence of Intra - American

Aug 15, 1993 - of model amino acid side chains, revealing that charges depend on the ..... and the C@ carbon, in the side chain of Me-L-Ser-Me and in ...
0 downloads 0 Views 1MB Size
Download PDF
J. Phys. Chem. 1993,97, 9797-9807

9797

Modeling Amino Acid Side Chains. 3. Influence of Intra- and Intermolecular Environment on Point Charges Christophe Chipot, JBnos G. AngyBn,' and Bernard Maigret' Laboratoire de Chimie ThCorique, Unit8 de Recherche Ass0cik.e au CNRS No. 510, UniversitC de Nancy I, BP. 239, 54506 Vandoeuvre-12s-Nancy Cedex, France Harold A. Scheraga' Baker Laboratory of Chemistry, Cornell University, Ithaca, New York 14853-1301 Received: March 23, 1993'

Modifications of the molecular charge distributions of simplified models of amino acids due to intra- and intermolecular environment have been analyzed. Point charges of the extended conformers of Me-L-Ala-Me, Me-L-Asp-Me, Me-Gly-Me, Me-L-Ser-Me, and the extended (C5) as well as the C7-equatorial conformations of N-acetyl-N'-methyl-L-aspartylamide(NANMAsp) have been computed from a b initio optimized geometries, using the split-valence 6-31G** basis set. Point charge models were determined from the resulting selfconsistent-field (SCF) wave functions by means of least-squares fits to the exact a b initio electrostatic potential and to the approximate electrostatic potential created by distributed multipole moments (DMM). The results show that the charge distribution of the side chain is affected by that of the backbone. In turn, the conformation of the backbone, and hence its charges, are altered, depending on the polarity of the side chain. In addition, the distortion of the total molecular charge distribution caused by conformational changes strongly suggests that conventional models are unable to provide point charges that are conformationally invariant. Moreover, the computed quasi-conformation-independent charges obtained from D M M potentials support the view that it is impossible to obtain conformationally invariant atomic point charge models of good quality. The effects of the solvent on point charge models have also been considered a t two levels of approximation. First, a cavity model of solvation has been used to represent the bulk aqueous solution. 6-31G** a b initio polarized wave functions have been employed to study the influence of the continuum on the charge distribution of a series of model amino acid side chains, revealing that charges depend on the surroundings. Second, in the framework of a supermolecule model, and in order to understand the role of the first solvation shell, two hydrated complexes of the model aspartate side chain have been optimized in both vacuum and aqueous solution, using the 6-31G** basis set. Point charges were computed using the corresponding accurate wave functions. The significant charge transfer observed from these results indicates that it is necessary to employ a combined discretecontinuum model of solvation for the treatment of explicit hydration.

1. Introduction

The use of point charges as an essential ingredientof empirical interatomic potential energy functions14 has given rise to a substantial number of alternative computational procedures.5-29 Some of the most representative approaches to compute point charge models have been reviewed in the first two parts of this ~eries.~8JOIt should be observed that, unfortunately, many of these approaches are inappropriate to reproduce the molecular electrostatic properties, viz. multipole moments, electrostatic potentials, and fields. In the preceding paper^,^^^^^ we have investigated some aspects of the method which consists of fitting atom-centered charge models to the ab initio SCF potentials or fields,using a constrained numericalleast-squaresprocedure.14 Charges obtained from this technique do not reflect any consistent transferability from one molecule to another structurally similar 0ne.28,30,31In order to remedy this situation, we have explored an alternative method based on the derivation of the potentials created by distributed multipolemoments (DMM).32-35 Net atomic charges computed with this recent alg0rithm~~92~ are perfectly transferable but have the major drawback that they are less successful in mimicking the exact electrostatic potential and its successive derivatives. In an intermediate approach, a limited number of additional off-atom sites are introduced to relax the fitting procedure and thereby provide a better description of the high-order molecular *Abstract published in Advance ACS Abstracts, August 15, 1993.

0022-3654/93/2097-9797%04.00/0

multipolemoments. Such a scheme yields perfectly transferable charges and much more satisfactory electrostatic properties.29.30 As an application, we have studied the molecular charge distribution of the naturally occurring amino acid side chains at the ab initio SCF level of a p p r o x i m a t i ~ n .An ~ ~ extra hydrogen atom was added to each side chain R of the amino acid moiety -NH-CHR-CO-, hence constitutinga completemodelmolecule, which was fully optimized28 using the polarization split-valence type 6-31G** basis et.^^,^^ Electrostatic potential- and DMM potential-derivedpoint charges were computed from the corresponding wave functions. In the present contribution to this series, we attempt to answer the following questions regarding the role of the intra- and intermolecular environment on the charge distribution in the R moiety, which will now be considered in its realistic amino acid environment: (i) Is the model molecule R-H sufficient to characterize the correspondingamino acid; i.e. what is the direct influence of the backbone on the charges borne by the constituent atoms of the side chains, and vice versa? Moreover, are these charges significantlymodified by conformational changes of the backbone? These problems are discussed in section 2 by means of an ab initio study of theextended conformers of the terminally-blocked amino acid models H3C-NH-CHR-CQCH3 (R = H, CH,, C H t OH, and CH2-COO-) as well as the extended (CS) and C7-equatorial conformationsof the N-acetyl-N'-methyhaspartylamide (NANMAsp). For each molecule in each conformation, 0 1993 American Chemical Society

9798

Chipot et al.

The Journal of Physical Chemistry, Vol. 97,No. 38, 1993

Me-(L)-Als.Me

Me.(L).Asp-Me

N- Acetyl-N’-methyl-(Lbaspartylamide lIH/O2

H4

a

Me.Gly.Me

\

,HIo

(extendedconformation)

a, Me.(L).Ser.Me

Figure 1. Terminally blocked models of four naturally occurring amino

acids. a full geometry optimization was carried out, and the net point charges were determined from accurate SCF wave functions. (ii) What are the effects of the surroundings on the molecular charge distribution? This point seems to be most important, due to the fact that classical molecular dynamics (MD) or Monte Carlo (MC) simulations of polypeptides are frequently carried out in a solvent, typically aqueous solution. If, indeed, use is madeof effective pair potentials (Le. noexplicit polarization terms in the global energy expression), the charge distribution models should reflect the solvated state of the fragments, in a potential of mean force sense. Ab initio calculations on a supermolecule consisting of an aspartate side chain and one or two water molecules have been used to elucidate the role of the first solvation shell. The influence of the bulk solvent on the electronic and molecular structure of the solute has been examined, using the self-consistent-reaction field (SCRF) model of solvation.384 The details of the results obtained for five representative side chains and for the hydrogen-bonded complexes, CH3C00-(H20), ( n = 1, 2), a t the 6-31G** level of approximation are presented in section 3.

2. Influence of Intramolecular Environment on Point Charges In this section, we consider the influence of the backbone as well as its conformational variations on the molecular charge distribution. The geometry of each of the four molecules (see Figure l), involving a side chain and a simplified model of the backbone, in its extended conformation, has been fully optimized. In order to obtain a more reliable description of the backbone, the representative example of the N-acetyl-N’-methyl-L-aspartylamide (NANMAsp) has been treated in two different conformations (Le. extended and C7-equatorial). The geometries in both conformations have also been optimized. In each case, electrostatic potential- and DMM potential-derived charges were computed using 6-3 1G** wave functions. The constrained numerical fit of net point charges to the a b initio S C F potential was carried out with the GRID employing a finite set of regularly spaced points (grid step Ar = 0.5 A) surrounding the molecule. Points lying inside the van der Waals envelope or further than the arbitrarily chosen exclusion radius (Le. rmx = 4.0 A) from any atom were discarded. This manner of selectingthe grid points does not constitute a limitation, and other algorithms (e.g. GEPOLW8) may easily be employed without modifyingthe fiial results.30 A description of this method is available in the first part of this study.28330 The electrostatic potential was evaluated on each point of the grid with the ab initio package Gaussian 92.49

N-Acetyl-N’-methyl-(Lbaspartylamide

(c’- equatorial conformation) Figure 2. Extended (Cs) and C7-equatorialconformationsof theN-acetylN’-methyl-L-aspartylamide(NANMAsp).

The analytical least-squares fit of net atomic charges to the potential generated by distributed multipole series was carried out with the MULFIT program.50 The procedure introduced by F e r e n ~ z yeliminates ~~ the explicit evaluation of the exact electrostatic potential or its derivatives on a grid of regularly or irregularly spaced points and replaces it by a direct integration over given regions of the space surrounding the molecule. An analysis of this approach is presented in previous paper~.~5929 The one-electron density matrix over Gaussian type orbitals (GTO) was provided by Gaussian 92.49 The distributed multipole analysis (DMA)32v34was carried out with a program written in the OSIPE e n v i r ~ n m e n tusing , ~ ~ the weight functions introduced by VignQ Maeder and C l a ~ e r i eand , ~ ~the multipole expansion was limited to the hexadecapole (I = 4); the boundaries of integration were fixed at P I = rvdw + 2.0 A and p2 = rvdW 5.0 A.25.29 Our criteria for the quality of the fit are the root-mean-square deviation (rmsd) between the quantum mechanically calculated and the point charge-derived electrostatic potentials together with the corresponding mean absolute deviation (A) respectively defined by

+

- P(r,,))’

Ill2

where NPnt denotes the number of points in the grid. VSCF(r,) is the reference ab initio S C F electrostatic potential evaluated at a point p, located at r,,, whereas VQ(rr) is the potential created by the set of net atomic charges. The values of the rmsd are given in atomic units. However, in order to correlate our results with those available in the literature, which are generally expressed in kcal mol-’, these values should be multiplied by a factor of 627.5095 (hartrees to kcal mol-1). The mean absolute deviations A are given in percent. Innuence of the Backbone on the MolecularChargeDistribution. In order to assess the mutual influence of the charge distribution

Modeling Amino Acid Side Chains of the backbone of an amino acid residue and its side chain, we have studied a set of three molecules H3C-NH-CHR-CO-CH3, involving three different substituents correspondingto increasing polarity, uiz. R = CH3, CH2-OH, and CH2-COO-. Supplemented by an extra hydrogen atom, each functional R group constitutes a model amino acid side chain, which was studied separately in the first two parts of this series of papers.28~30 In the present work, we also have considered the case, where R = H, correspondingto an -unperturbed” backbone that may, hence, be regarded as a reference molecule. The geometry of the four molecules, referred to as Me-L-AlaMe, Me-L-Asp-Me, Me-Gly-Me, and Me-L-Ser-Me and considered in their extended conformation (see Figure l), have been fully optimized using the extended split-valence type 6-31G** basis set. Their point charges were derived from the ab initio SCFelectrostatic potential as well as from the potential generated by distributed multipole series. From these point charge models, the corresponding molecular dipole moments were computed. Our results are reported in Tables I and 11. It may be observed that the DMM potential-derivednet atomic charges of the backbone are transferable (Ferenczy’s method25), and for the polar parts of this set of molecules,they are reasonably close to those obtained from the SCF electrostaticpotential,which, nevertheless,are influenced to a great extent by the environment. This agreement with the two kinds of fitting procedure (e.g.the charges borne by the CO fragment in Me-Gly-Me are qCDMM = 0.612 ecu [electronic charge units] and qoDMM = -0.603 ecu, whereas qcv = 0.633ecu and qov = -0.575 ecu) is essentially due to the fact that each molecule involves several polar sites along its backbone, which counterbalance the effects of the nonpolar parts (typically, terminal methyl groups). Consequently, the magnitude of the errors in the fitting procedure is generally limited (albeit still large in most cases). The key question is, now, whether the previously used model molecule R-H provides an acceptable representation of the corresponding true side chain R group placed in its natural peptide or protein environment? Comparison with the results reported in refs 28 and 30 shows that, when R = CH3 [i.e.Me-~-Ala-Me], the charge borne by the C@carbon is qcv = -0.638 ecu, which is fairly close (within a 10% range of error) to that in the model side chain, i.e. CH4 (uiz. qcv = -0.561 ecu). It is also worth noting that, in the case of DMM potential-derived charges, a consistent C6+-Hb polarity is observed, unless the carbon is directly connected to a polar fragment. Moreover, since the side chain is not isolated,**JOthe charges of the three hydrogens do not reflect the tetrahedral symmetry appearing in CH, and, as a result, are substantially different from one another (uiz. 0.130 5 qHv 5 0.209 ecu to be compared to qHv = ecu in CH4). These remarks are also valid for Me-L-Asp-Me, where the charges borne by the oxygensin the-COO-fragment arelogically equivalent, if the intramolecular environment is not taken into account. This is practically the case (i.e. qcv = 0.847 ecu and qov = -0.797 and -0.804 ecu, to be compared with those in the model R-H side chain: qcv = 0.935 ecu and qov = -0.854 ecu). The electrostatic potential-derived charge of the C@carbon (qcv = -0.473 ecu) is slightly larger than that in the model side chain (qcv = -0.303 ecu). Incontrast, thechargesobtainedfromDMM potentials are almost identical (e.g.qCDMM = -0.122 and -0.124 ecu in Me-L-Asp-Me and in the isolated side chain). Concerning Me-L-Ser-Me, the charges of the OH fragment are qov = -0.532 ecu and qHv = 0.378 ecu, which agrees with those in the corresponding model side chain (uiz. qov = -0.685 ecu and qHv = 0.426 ecu). However, in the case of the c@ carbon, there is a significant deviation (qcv = -0.193 and 0.278 ecu, in the model amino acid and in the simplified side chain, respectively). This discrepancy is noticeably reduced when additional sites are added in the fitting procedure (viz. qcv = -0.623 and -0.337 ecu, in Me-L-Ser-Me and in the isolated side chain, respectively). As

The Journal of Physical Chemistry, Vol. 97, No. 38, 1993 9799 for both Me-L-Ala-Me and Me-L-AsgMe, the net atomiccharges derived from the potential created by distributed multipole series are very similar to those in the side chain model (uiz. q$MM = -0.7 15 and -0.723 ecu, qHDMM = 0.414 and 0.404 ecu and qCDMM = 0.458 and 0.628 ecu, for the oxygen, the hydroxyl hydrogen, and the C@carbon, in the side chain of Me-L-Ser-Me and in the model side chain, respectively). In the present context, the best results are reached with point charge models involving a limited number of off-atom sites, from both a transferability and an accuracy point of view. With these models, we observe that the charges borne by the backbone of the four terminally blocked amino acid models are comparable (see Table I) and that the charges of the different side chains are in better agreement with those in the simplified R-H models**.30 than when use is made of a classical atom-centered charge model. This may be explained by a better description of the high-order contributions to the multipole expansion of the potential, hence leading toa smaller number of contradictory fits. However,unlike the analytical fit to the DMM potentials, this procedure does not overcomethe short-range influenceof the side chain which greatly modifies the charge borne by the Ca carbon of the extended backbone (uiz. qcv = -0.775, 4 2 5 1 , 4 4 9 3 , and 0.020ecu, to be compared with qCDMM = 0.364,0.335,0.338,and 0.279 ecu, in Me-L-Ala-Me, Me-L-Asp-Me, Me-Gly-Me, and Me-L-Ser-Me, respectively). In addition to these observations, one should consider the reciprocal influence of the side chain on the charge distribution of the backbone. As an attempt to elucidate this aspect of the problem, we used the followingapproach, in which the functional group R of each of the three complete amino acid models MeL-Ala-Me, Me-L-Asp-Me, and Me-L-Ser-Me was removed and replaced by a hydrogen atom. In all cases, the C e H bond length was corrected to dc-H = 1.09 A, which roughly corresponds to the optimized interatomic distance in the true Me-Gly-Me. The ab initio SCF electrostatic potential of the three resulting pseudoMe-Gly-Me nonreoptimized structures (Le. reflecting the geometry of the backbones in the presence of the alanine, aspartate, and serine side chains) has been evaluated on a grid of regularly spaced points, using the 6-31G** basis set and the net atomic charges have been derived, applying the conventionalconstrained least-squares meth0d.1~28 The charges obtained in each case reveal a few discrepancies with the reference molecule Me-GlyMe (e.g.the charge borne by the Ca carbon is qcv = 0.067,O.127, and 0.142 ecu for the transformed Me-L-Ala-Me, Me-L-AspMe, and Me-L-Ser-Me, respectively,which is compared with the original value of qcv = 0.175 ecu in the reference Me-Gly-Me). The fitting error (uiz. 0.562 X Irmsd I0.572 X lo-’ au) as well as the recalculated dipole moment (uiz. 2.948 5 cc I3.11 8 D) is reasonably close to the reference values (see Tables I and 11). In brief, since the modified molecules are not geometryoptimized after the transformation consisting of replacing the side chain by a simple hydrogen atom, the small deviations encountered in these calculations may be explained by slight conformationalchanges induced by the side chain R in the original models. The magnitude of these modifications is evidently dependent on the nature of the side chain, and large effects are generally expected with very polar groups. This demonstrates clearly the difficulties that will be encountered in a future use of point charge models obtained from small fragments of macromolecules. Taking into account the importance of inductive effects, manifested by the mutual influence of the side chain and backbone charge distributions,inevitablyimplies that polarization should be treated explicitly in order to provide a completely satisfactory p r o c e d ~ r e . ~ ~ From these examples, we conclude that (i) the charges borne by the amino acid side chain are certainly influenced by the molecular charge distribution of the associated backbone, (ii) the representation of the true side chain by a simplified model

9800 The Journal of Physical Chemistry, Vol. 97, No. 38, 1993

Chipot et al.

TABLE I: Cartesian Coordinates* and Potential-Derived (4")and DLptdbuted Multipole-Derived ((IDMM) Net Point Charges' of Several Model Peptides, Using the Split-Valence 6-31C** Basis Set; Influence of the Inclusion of Additional Sites 4v atomC

Hs C1

c;

H6 H7

N

atomic atomic centers and bond only centersd

X

Y

1.311 0.224 1.040 2.579 3.309 2.948 2.539 0.319 -1.095 -2.376 -2.270 -2.641 -3.168

-0.45 1 0.322 -1.410 -0.307 -0.958 0.709 -0,548 0.437 -0.425 0.300 0.883 0.985 -0.425

-0.339 -0.748 0.21 1 0.365 -0.403 0.375 0.340 -0.096 0.108 -0.126 0.080 0.246 1.406 0.039 1.299 0.021 0.616 0.021 0.372 -0.532 1.281 0.122 -0.428 0.166 0.487 0.159

-1.234 -0.715 0.415 -0.546 0.138 0.134 0.098 0.156 1.935 -0.502 0.167 0.176 0.182

-1.944 -0.656 -2.623 -2.382 -3.355 -1.692 -2.464 -0.475 -0.667 0.633 1.094 1.340 0.435 -1.682 0.5 18 0.308

-0.643 -0.141 0.082 -1.845 -2.141 -2.659 -1.747 -0.382 1.383 2.106 1.734 1.905 3.169 1.995 -0.751 -1.811

-0.389 0.071 -0.291 0.278 -0.107 0.08 1 1.366 1.120 0.020 0.276 1.183 -0.521 0.340 -0.175 -0.709 -0.831

-0.808 0.351 0.389 -0,143 0.093 0.093 0.043 0.041 0.621 -0,441 0.120 0.153 0.093 -0.618 -0.473 0.165

-1.168 -0.25 1 0.412 -0,543 0.106 0.136 0.092 0.139 2.655 -0.446 0.203 0.190 0.150 -0.579

1.334 0.250 2.642 3.394 2.758 2.849 1.169 0.213 0.336 -1.082 -2.320

0.293 -0.569 -0.199 0.524 -0.397 -1.123 1.213 -1.437 -0.966 0.157 -0.704

-0.299 -0.742 0.077 0.175 0.064 -0.223 -0.231 0.132 1.133 0.081 -0.467 0.126 0.053 0.41 3 -0.577 0.043 1.099 0.015 0.028 0.633 -0.047 -0.397

-1.315 -0.493 -0.533 0.144 0.090 0.131 0.420 0.179 0.128 2.348 -0.569

1.323 0.228 1.062 2.597 3.335 2.945 2.566 0.346 -1.077 -2.375 -2.307 -2.612 -3.162 -1.060

-0,501 -0.046 -1.352 -0.637 -1.008 0.327 -1.316 -0.347 -0.692 -0.159 0.034 0.778 -0.873 -1.606

-0.346 0.482 -0.798 0.329 -0.373 0.685 1.183 1.530 0.010 0.566 1.632 0.072 0.370 -0.756

1.577 0.23 1 1.628 0.218 -0.639 -0.217 -0.268 0.531 -0,439 -1.531 -1.889 -2.037

-0.436 -0.044 -1.268 0.217 -1.288 -2.261 1.172 1.905 0.893 1.902 2.857 1.454

-0.240 0.141 -0.781 1.190 -0.072 -0.668 -0.661 -0.669 -1.699 -0.127 -0.799 0.924

Z

-0.550

0.127

-0.525 -1.200 -0.100 0.020 0.433 0.320 -0.538 0.009 0.145 0.056 0.138 0.045 0.106 0.030 0.098 0.097 2.532 0.739 -0.543 -0,538 0.132 0.174 0.180 0.210 0.151 0.200 -0.554 -0.468 N-Acetyl,."-methyl, -0.801 -0.858 0.620 0.088 0.305 0.426 -0.114 0.156 0.519 3.441 -0.609 -0.518 -0,173 -0.850 0.024 0.193 -0.02 1 0.164 2.348 0.899 -0.806 -0.603 -0.653 -0.864

atomC X Y Z Me-L-Ala-Me; Npte = 4041 0 -1.117 -1.554 -0.364 -0,930 1.725 -0.408 0.186 0.364 CS 2.360 0.057 -0.560 Hg 0.381 1.148 2.206 -0.284 0.599 Hio 1.659 -1.472 -0.016 -0.088 HI1 0.768 -0.065 -0,064 -0.096 XNCl 1.945 -0,379 0.000 -0.120 XNCl -0.435 -0.051 0.116 -0.094 XClCl -1.735 -0.062 0.197 0.614 XClC4 -1.106 -0.989 -0,171 -0.029 XCl0 1.023 -0.098 0.205 -0,015 XCIC, 0.009 rmsd ( 10-3 au) A8 (%) 0.021 Me-L-Asp-Me; Npt = 4045 0.577 -0.334 -1.710 -0.931 Hio 1.899 -0.648 -0.001 0.335 c 6 2.865 -0.5 15 -0.750 0 2 0.362 1.235 1.866 -0.723 0.646 00 -1.300 -0.392 -0.159 -0.135 XNCl -2.162 -1.244 -0.055 -0.098 XClCI 0.621 0.045 -0,661 -0.138 XCA 1.744 0.148 -0,017 -0.060 XClC4 -1.174 1.689 -0,078 0.643 XClOl -0.069 -0.446 -0,319 -0.069 XClCJ 1.208 -0,699 -0.355 0.035 XCJC, 2.382 -0.58 1 -0.375 0.026 XC602 1.883 -0.685 0.617 -0.020 xc601 -0.672 rmsd ( au) -0.122 A (%I -0.003 Me-L-Gly-Me; Npt = 3830 -2.253 -1.548 0.634 -0.846 H7 3.194 -0,111 0.183 0.338 Ha -2.415 -1.104 -1.053 0.583 Hg 0 -0.088 -1.141 1.348 0.073 0.792 -0,138 -0.111 -0.120 XNCl 1.988 0.047 -0.117 -0,089 XClC2 -0.416 -0,206 0.053 0.346 XClCl -1.701 -0.274 -0.009 -0.054 XC,C, -1.111 0.752 0.050 -0.065 xc,o 0.6 12 rmsd (1 au) -0.047 A (%I Me-L-Ser-Me;Npt = 4084 0.151 1.481 0.419 -0.909 cs -0.63 1 1.866 1.061 0.279 Hg 1.897 0.772 1.092 HII 0.389 1.900 -0.880 -0.131 0.525 0 2 1.464 -1.454 0.484 -0.068 Hi2 0.776 -0.273 0.068 -0,077 XNCl 1.960 -0.569 -0.009 -0.099 XClC2 -0.424 -0.369 0.246 -0.079 XClCl -1.726 -0.425 0.288 0.625 xctc4 -1.069 -1.149 -0.373 -0,027 xc301 0.190 0.7 18 0.451 -0.022 XClCJ 0.596 0.700 -0.613 0.024 XCJO2 0.020 rmsd ( l t 3 au) -0.6 13 A (W) -L-aspartylamide(Extended Confcirmation);hNPt = 4791 2.693 0.125 0.238 -0,709 CS 2.718 1.060 1.000 0 4 0.240 -1.866 -1.217 0.438 0.321 Nz -2.154 -0.304 0.773 -0.001 Hs -2.825 -2.264 0.209 0.859 c6 -3.139 -2.014 0.735 -0.777 H6 -3.055 -2.386 -0.846 -0.178 H7 -2.458 -3.215 0.579 0.048 H8 3.992 -0.485 -0.262 -0.012 c7 4.513 0.258 -0.855 1.021 Hg 4.614 -0.718 0.594 -0,857 Hio 3.851 -1.379 -0.859 -0.904 HI1 ~DMM

9v atomic atomic centers and bond only centersd

~DMM

-0.547 -0.638 0.130 0.209 0.171

-0,543 -0.654 0.118 0.150 0.143 0.399 0.686 -0.358 -0.595 -0,345 0.654 0.109 10.36

-0.612 0.093 -0.059 -0.007 -0.03 1

0.099 2.463 -0.649 -0.621 0.069 0.679 -0,821 -0,944 -0.554 0.640 -0.498 -0.7 15 -0.819 0.105 0.09

-0.03 1 1.020 -0.889 -0,900

0.179 0.189 0.185 -0.433 0.629 0.652 -0.668 -0.558 -0.704 0.134 6.78

-0.009 0.028 0.014 -0.603

0.458 -0.067 -0,057 -0.7 15 0.414

0.419 22.32

-0.623 0.165 0.126 -0.844 0.449 0.190 0.630 -0.901 -0,747 -0.652 0.324 0.571 0.199 8.11

0.9 12 -0.658 -0.442 0.395 -0.592 0.203 0.194 0.197 -0.785 0.201 0.196 0.200

2.407 -0.467 -0.744 0.465 -0.507 0.129 0.128 0.131 -0.595 0.172 0.167 0.153

0.848 -0.727 -0,769 0.386 0.531 -0.090 -0.086 -0.087 -0.039 0.006 0.003 -0.026

0.327 20.72 0.074 0.847 -0.797 -0.804

0.335 0.28

0.07 1 0.138 0.108 -0.575

0.508 20.36 -0,193 0.140 0.170 -0.532 0.378

0.769 63.33

0.88 1 0.75

0.926 41.32

0.801 31.58

Modeling Amino Acid Side Chains

The Journal of Physical Chemistry, Vol. 97, No. 38, 1993 9801

TABLE I (Continued) IV

UV

atomic atomic centers andbond only centersd

atomic atomic centers and bond only centersd

~

X

Y

Z

0.904 -0.204 -0.428 -0.019 -0.900 -1.710 -1.784

-0.240 -0.666 -1.775 0.564 1.537 2.379 1.678

-0.049 0.035 -0.370 -0.260 -0,394 -0.463 0.398

0.233 -0.157 0.357 -0,320 -0.485 0.338 -1.676 -1.985 -1.901 -2.586 -2.052 -3.758 1.384 2.308 1.664 2.219 2.313 2.432 1.735

1.066 -0.120 -1.416 -2.402 1.323 -0.080 -0.118 -1.073 0.635 0.207 0.800 -0.104 1.706 1.335 -1.394 -0.630 -2.586 -3.339 -3.026

0.06 1 0.802 0.158 0.031 -0.608 1.769 1.012 1.414 1.765 -0.210 -1.171 -0.066 0.187 0.893 -0.196 0.123 -0.678 0.099 -1.479

qDMM atomC X Y Z N-Acetyl-N’-methyl-L-aspartylamide (Extended Conformation);* Nmt= 479 1

0.547 -1.542 -1.009 0.473 -0.295 -0.731 -0.721

XN~Cs XC@4 xCi”I XNlG

xc5c7

2.135 2.706 -1.253 -2.346 3.342

-0.155 -0.001 0.593 0.619 -1.253 0.183 -1.740 0.324 -0.180 -0.012

rmsd ( 104 au) A ($1 N-Acetyl-N’-methyl-L-aspartylamide (C; Conformation);‘ Npt = 4753 -0.507 -0.776 -0,720 Ha 3.294 -2.322 -1.057 c 7 1.519 2.977 -0.627 0.297 -0,033 0.214 0.614 3.271 0.875 Hg 0.647 3.186 -1.232 -0.576 -0,459 -0.723 Hio 1.692 3.803 0.055 0.285 0.449 0.371 Hi1 2.392 2.890 -1.263 XNiCi 0.038 0.473 0.432 -0.038 0.138 -0.015 -0.247 -0.563 0.100 -0,768 0.480 XClCl -0.144 0.032 0.181 0.017 XClOl 0.018 -1.909 0.094 0.036 0.139 -0.012 0.907 XCICI -0,916 -0,119 0.843 2.505 1.039 XCIC, -2.131 0.044 0.401 XClO, -2.319 0.504 -0.690 -0.812 -0.642 -0.914 -0,858 XC401 -3.172 0.052 -0.138 -0.769 -0.609 0.800 2.410 0.850 XN~Cs 0.808 1.386 0.124 XCP4 1.846 1.520 0.540 -0.698 -0.535 -0.792 -0.457 -0.782 xCiNi 1.011 -1.405 -0.019 -0.803 0.362 0.444 1.989 -1.990 -0.437 XNiC( 0.371 -0.542 -0.494 XCFl 1.451 2.342 -0.220 0.531 0.169 0.115 -0.093 mud ( l k 3 au) 0.206 0.155 -0.061 A (%)

0.283 0.21 0.175 -0.821 0.226 0.216 0.205

0.231 0.18

~~

-0.736 -0.851 -1.027 0.501 -0.377 0.090 0.07 0.112 -0.578 0.174 0.168 0.164 0.514 -1.103 -0.968 0.120 4457 -0.802 -0,744 -0.723 -0.799 -1.076 0.521 -0.434 0.086 0.07

qDMM

1.140 1.01 -0.103

0.O00 -0.013 -0,006 -0.012

1.538 1.45

0 In angstroms. b In electronic charge units (ecu). Indices given in Figures 1 and 2. Bonds involving heavy atoms only. e Number of points requested for the fitting procedure. /Root-mean-square deviation between the ab initio computed and the point charge derived electrostaticpotentials and fields. g Mean absolute deviations. Side chain conformation: x1 = -165.5’, xz = 0.8’. ‘Side chain conformation: XI = 38.8’, x~ = 158.0D.

*

TABLE II: Dipole Moment’ of Modified Amino Acids from Potential- and Distributed Multipole-Derived Net Atomic Charges, Using the Split-Valence 6-316** Basis Set; Comparison with SCF Values and Influence of the Inclusion of Additional Sites PV

atomic atomic centers and bond only moleculeb centersc PDMM PSCF 2.916 2.909 2.881 2.901 Me+- Ala-Me 7.439 7.536 7.436 7.442 Me+ AspMed Me-Gly-Me 2.994 2.981 2.928 2.979 3.065 3.057 2.974 3.054 Me-L-Ser-Me 9.158 N-acetyl-N’-methyl-L-aspartyl9.759 9.851 9.758 amide (extended conformation) N-acetyl-N’-methyl+ 10.385 10.388 10.301 10.388 aspartylamide (C& conformation) 0

In debye. b See Figures 1 and 2. Bonds involving heavy atoms only.

d Charged species;dipole moment calculated in the standardorientation’3

framework of Gaussian 92.49 molecule (vir. R-H)constitutes an economical and satisfactory approximation, and (iii) the polarity of the side chain contributes to a slight structural alteration of the backbone and, as a result, to a distortion of its molecular charge distribution. Such a modification appeared to be relatively minor in the latter few examples. T h a t this is not always the case, however, may be seen in the details in the following section. Influence of Conformational Modifications on the Molecular Charge Distribution. In this paragraph, we consider the charge distribution of two different conformers of the N-acetyl-N’-methylL-aspartylamide (NANMAsp). This compound constitutes a case example for two major reasons: (i) it is a charged species, and (ii) it involves peptide bonds (one on each side of the central amino acid) that are responsible for induction effects within the molecule. T h e geometry of both conformations was optimized, using the extended split-valence type 6-31G** basis set. T h e

extended, or C5, conformer, with 4 = -153.6O and $ = 168.7O, is characterized by an intramolecular O1-.HI hydrogen bond ( d o , H , N 2.10 A), and the C7-equatorial conformer, with 4 = -82.6’ and $ = 59.3’, shows a strained O p H 5 hydrogen bond

(doa,- 2.11A)andslightlyrotatedpeptidegroups(ol =-175,4O and w2 = -179.1’). T h e backbone structures are, somehow, very close t o those computed for t h e N-acetyl-”-methyl-L-alaninamide, using a double-rand triple-fplus polarization basis set (DZP).S5 Our a b initio computations a t the 6-31G** level revealed that the extended, or C5,conformation of N A N M A s p is approximately 4.05 kcal/mol lower in energy than the C7equatorial conformation, in good agreement with empirical calculations, using the ECEPP/356 force field. I t is worth noting that conformational changes are a source of occasionally significant variations of the geometry. In particular, from the extended to the C7-equatorial conformers of NANMAsp, r(NCaC’) increases by 7.0°, r(C’NCa) increases by 0.9O, and r(C’CaCe) decreases by 2.6O. Simultaneously, the peptide bond length dC5N, decreases by 0.014 A, whereas d c a N a increases by 0.024 A. A similar trend has been observed recently for other oligopeptides,57 using the 4-2 1G basis set.58 Consequently, it has been suggestedS7 that the geometrical changes due to variations of the conformation should be taken into account in force field parametrization. From the above-minimized conformations, a b initio SCF electrostatic potential derived- and DMM potential-derived point charges were computed using the correspoding wave functions. T h e geometries in Cartesian coordinates, together with the charges, are listed in Table I. T h e molecular dipole moments computed from the fitted point charge models are presented in Table 11. T h e first comment suggested by our results concerns the problem of conformational invariance. A s expected, potentialderived net atomic charges show a limited transferability from the extended to the C7-equatorial conformer of the N A N M A s p

Chipot et al.

9802 The Journal of Physical Chemistry, Vol. 97, No. 38, 199‘3

(viz. the charge borne by the Ca carbon is qcv = 0.620 and 0.297 ecu, for the extended and C7-equatorial conformations, respectively). On the contrary, DMM potential-derived net atomic charges are consistently transferable, with no major deviation (uiz. the net charge of the Ca carbon is qCDMM = 0.240 and 0.214 ecu, for the extended and C7conformations, respectively),hence indicating that, in this particular case, conformational changes affect the charge distribution moderately. The origin of this observation probably lies in the background of the method, in which partial charges are fitted to individual potentials evaluated at different locations of the molecule, rather than net charges to global potentials.25329 The distortion of the charge distribution caused by conformational modifications is, hence, likely to have a lesser effect on the net atomic charges determined from distributed multipole series, as far as (i) the multipoles of each atomic site can be described satisfactorily by a simple set of net atomic charges, located on the nearest neighboring atoms, and (ii) no significant variations occur in the local fields experienced by any particular atom. On the other hand, in the case of a global fit (e.g. to the SCF electrostatic potential or field), a more extended use of the furthermost sites and the couplings of the partial fitting proceduresz9may occur to a much larger extent, so that the resulting net atomic charges are much more dependent on the conformation. It may also be noted that, owing to the relatively large rmsd obtained with the charges computed from the distributed multipole series, contradictory partial fits have apparently occurred. It, therefore, suggests that the description of the exact potential by the present point charge models is obviously quite poor. The inclusion of a limited number of off-atom sites located at the center of the bonds involving heavy atoms exclusively, which admittedly constitutes a reasonable compromise between charge transferabilityand reliable reproduction of electrostaticproperties, indeed contributes to a better transferability of the charges from one conformer to the other. However, as may be observed in Table I, this is not completely achieved, therefore indicating that such a model is not fully adequate to yield conformationally invariant charges. A second remark concerning the transferability within the side chain is suggested by a variation of the charges from Me-LAsp-Me to the extended conformer of NANMAsp (see Table I). The alteration of the charge distribution is a probable consequence of small conformational changes within the backbone of NANMAsp, together with the induction effects resulting from the dipole moments generated by its peptide groups. From this, we conclude that (i) the molecular charge distribution is undoubtedly affected by conformational changes, so that the point charge models derived from electrostatic quadtities are altered accordingly, (ii) these effects may be damped to a certain extent by including a small number of additional sites in the course of the least-squares fitting procedure to the ab initio SCF potential, and (iii) fully conformationally invariant charges cannot be obtained by simple point charge models59-52 which suggeststhat explicit polarizabilities should be included, as it has been demonstrated in ref 54. 3. Influence of Intermolecular Environment on Point Charges

Computational Details. Although quantum chemistry computations have reached a relatively high degree of accuracy, they deal mainly with isolated molecules, correspondingto low-pressure gaseous states. However, numerous physical, chemical, and biological properties are studied experimentally in condensed phases, especially in liquids. It is therefore of great importance to take the solvent into account when studying the conformational behavior of biopolymers. The physical reality of the liquid may be approached by means of a large scale statistical mechanical study of model molecules of fixed charge distribution in order to reproduce the thermo-

Alanine

Aspartate

Arginine

Cystine

Histidine D

Figure 3. Side chains of five naturally occurring amino acids, with an extra hydrogen atom to create a complete model molecule.

dynamical properties of the system. Such an approach has the major drawback that it does not correctly describe the influence of the surroundings on the molecular electron distribution. In contrast, the self-consistent-reactionfield (SCRF) methods put the emphasis on the quantum chemical representation of the electronic structure of a solute molecule embedded in a rudimentary model of the solvent. In the present study, the bulk solvent effects on some naturally occurring amino acid side chains will be taken into account by means of the cavity model proposed by Rivail et al.39.63-66 The solute molecule is placed in a cavity surrounded by a contin~um6~ characterized by specificdielectric properties. This approach has been widely applied at different levels of s o p h i s t i c a t i ~ n The . ~ ~molecular ~ ~ ~ ~charge ~ ~ ~ distribution ~~ within the cavity polarizes the dielectric medium, which, in turn, generates an electrostatic potential in the cavity, hence perturbing the electronicand molecular structureof the solute. The potential created in return by the continuum is referred to as a reaction field, defined by

where Qlm denotes the mth component of the multipole moment of rank 1 and f y ‘ i s the reaction field factor, a tensor depending on the dielectric characteristics of the medium as well as the geometry of the cavity. This factor can be evaluated analytically in the case of spherical or ellipsoidal cavities.65J3 Introducing the perturbation due to the solute-solvent interactions into the classical HartreeFock equation, one obtains

(3) is an element of the unperturbed Fock matrix for the Gaussian basis functions (x,) and 1 ~ ” ) . Finding the equilibrium geometry of the solute implies a minimization of the total energy

with respect to the nuclear coordinates. 19)stands for thegroundstate determinant. Details of the complete geometry optimization procedure in the ellipsoidal cavity approximation are reported elsewhere.74 This formalism has been implemented in the Gaussian 88 program75 as a set of independent links.76 The following computations were performed with a multipole expansion up to the order 1 = 6.

Modeling Amino Acid Side Chains

The Journal of Physical Chemistry, Vol. 97, No. 38, 1993 9803

TABLE III: Muence of the Solvent8on the Potential-Derived (qv) and Distributed Multipole-Derived (-) Net Point Chargesb of Several Amino Acid Side Cham, Using the Split-Valence 6316** Basis Set; Comideration of a Limited Number of Additional Sites 4v

4v

atomic atomic and atomic atomic and centers only centers only bond centersc 4DMM bond ccntersc in in in in in in in in in in atomd vacuo solution vacuo solution vacuo solution atomd vacuo solution vacuo solution Alanine: Nmte= 2263 0.172 0.156 H4 0.140 0.145 -0,561 -0.580 RMSW (lo-’ au) 0.049 0.050 -0.043 -0.039 0.140 0.145 r--

0.140 0.140

0.145 0.145

-0,425 0.138 0.116 0.116 0.203 0.007 0.007 -0.005 0.080 0.080 -0.488 0.322 0.909 -1.085

-0.436 0.113 0.112 0.112 0.249 -0.006 -0.006 0.009 0.075 0.075 -0.484 0.326 0.866 -1.076

-0,043 -0.039 -0,043 -0,039 -0.585 0.156 0.133 0.133 -0.651 0.126 0.126 -0.763 0.183 0.183 -0.722 0.433 3.234 -1.021

-0.630 0.135 0.134 0.134 -0.728 0.133 0.133 -0.778 0.184 0.184 -0.738 0.448 3.432 -1.006

0.126 -0.003 -0.033 -0,033 0.127 -0.037 -0,037 0.362 -0.041 -0.041 -0.779 0.365 1.322 -1.117

A#(%)

Arginine; Npnt= 4434 N3 Hg Hio Hii Hiz

0.146 -0.033 -0.043 -0,043 0.124 4.037 -0,037 0.383 -0,044 -0,044 -0.798 0.374 1.321 -1.129

37.34

43.77

QDMM

in vacuo

in solution

-0.043 -0.039 0.558 0.561 133.16 130.78

-0.902 -0.886 -0.880 0.507 0.527 0.525 0.515 0.509 0.528 0.466 0.478 0.505 0.438 0.443 0.501 0.460 0.460 0.623 -1.038 -0.755 -0398 0.219 0.231 0.108 0.16 0.15 0.08

-0.836 0.547 0.526 0.518 0.505 0.536 0.483 0.663 -1.128 -0,836 -1.018 0.100 0.08

-1.114 0.477 0.484 0.485 0.487

-1.124 0.500 0.484 0.503 0.496

0.764 0.55

0.774 0.57

-0.854 -0,903 -0.718 -0.083 XCIOl -0.492 -0.492 XC202 rmsd ( au) 0.325 0.269 0.130 A(%) 0.22 0.18 0.09 Cystine (+ Lone Pairs); NPt = 2824 0.192 Xqs, 0.072 -0.002 0.786 XSlS2 0.072 0.786 xs2c2 0.192 xi -0.806 -0.810 -0.760 Xi) -0.779 -0.794 -0.738 0.008 -0,025 xz -0.806 -0.810 -0.160 -0.001 XZl -0.779 -0.794 -0.738 rmsd (lt3 au) 0.170 0.195 0.154 0.008 A(%) 7.75 -0,025 8.32 7.68 -0.001 Histidine D; NPt = 3295 0.179 H5 0.121 0.134 0.224 -0.012 H6 0.123 0.129 0.225 0.530 -0.012 XClCZ -0.019 -0.025 xc2c4 0.302 0.031 XC~N~ 0.560 -0.530 XN~C3 1.028 0.513 XC~Ni 1.042 -0.708 xN2C4 0.019 rmsd (lk3 au) 0.391 0.401 0.104 7.97 3.63 A(%) 12.52 0.396

-0.767 -0.147 -0.509 -0.509 0.1 13 0.08

-0,918 -0.957

XCIC2 XC2C3 XC& XNIC, XC4N2

XW3 rmsd (10-3 au) A(%)

Aspartate; NPt = 2447 -0.303 -0.342 -0.628 -0.597 -0.124 -0.174 0.020 0.052 0.099 0.118 -0.035 -0.001 0.036 0.075 0.112 0.140 4.022 0.019 0.020 0.052 0.099 0.118 -0.035 -0.001 0.935 0.969 1.822 1.918 1.053 1.072 -0.854 -0.903 -0.7 18 -0.767 -0,918 -0.957 -0.576 -0.565 -0.632 -0.639 0.205 1.668 1.644 1.532 1.484 0.811 1.668 1.644 1.532 1.484 0.811 -0.576 4.565 -0.632 -0.639 0.205 0.185 0.195 0.194 0.207 -0.008 0.148 0.160 0.158 0.173 -0.044 0.160 0.170 0.174 0.189 -0.018 0.185 0.195 0.194 0.207 -0.008 0.148 0.160 0.158 0.173 -0,044 0.160 0.170 0.174 0.189 -0.018 -0.693 -0.721 -0.681 -0.709 0.205 0.227 0.145 0.166 0.205 0.227 0.145 0.166 0.188 0.191 0.175 0.176 0.229 0.230 -0.301 -0.490 -0.478 -0.464 -0,931 -1.007 0.247 0.272 -0.733 -0.769 -0,573 -0,663 -1.430 -1.519 0.046 0.032 -0.715 -0.851 0.380 0.406 0.434 0.468

0.217 -0,035 -0.035 -0.029 -0.002 -0.523 0.495 -0.639 0.062 0.362

0 2

XClC2

0.090 -0.014 0.090 -0.751 -0,747 -0.751 -0.747 0.173 7.60

0.247 0.244 0.571 0.145 0.417 0.629 1.034 1.081 0.108 5.87

0.681 0.46

0.531 0.33

-0.480 -0.466 -0.480 -0.466 0.875 54.89

-0.48 1 -0.479 -0.48 1 -0.479 0.901 49.07

0.024 0.101

0.038 0.1 12

0.901 44.66

0.859 31.14 = 78.3). In electronic charge units (ecu). Bonds involving heavy atoms only. Indices given in Figure 3.

0 Continuum model of solvation (e Number of points requested for the fitting procedure. f Root-mean-squaredeviation between the ab initio computed and the point charge derived electrostatic potentials. 8 Mean absolute deviations.

e

Results and Discussion. The goal of this section is to examine several aspects of the solvation process and its effects on point charge models. As a first approach, we have considered the direct influence of the bulksolvent (uiz. water, in this work), represented by a dielectric continuum (i.e. e = 78.3) on the molecular charge distribution of a variety of model naturally occurring amino acid side chains (see Figure 3). The geometry of these molecules was kept fixed, corresponding to the minimized structures in vacuum, a t the 6-31G** level of approximation.28 In other words, we have studied the effects of a perturbation due to the s o l u t e solvent interactions on the ground-state wave function. This procedure is not limited to the HartreeFock level of approxi-

mation and may also be carried out if intramolecular electron correlation is taken into account.77J8 Ab initio SCF molecular electrostatic potential derived- and DMM potential-derived point charges were computed using the perturbed 6-31G** wave function (see Table 111). Dipole moments recalculated from thecorresponding point chargemodels are presented in Tavle IV. From these results, we may note that (i) there is a general trend for the absolute value of the charges to be slightly increased when the solute is surrounded by the dielectric medium. This is obviously related to polarization effects induced in the solute by the reaction field, which distort the original charge distribution.

9804

The Journal of Physical Chemistry, Vol. 97, No. 38, 1993

Chipot et al.

TABLE I V Effects of the Solvent' on the Dipole Moment* of Several Amino Acid Side Chains, Using the Split-Valence 6-31G** Basis Set; Influence of the Inclusion of a Limited Number of Additional Sites U"

atomic centers only atomic and bond centersc PDMM PSCF side chain in vacuo in solution in vacuo in solution in vacuo in solution in vacuo in solution 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 alanine 6.666 5.838 6.692 5.783 5.785 6.663 6.665 arginined 5.782 4.700 3.794 4.697 3.901 4.726 3.793 4.697 3.799 aspartated 2.391 2.383 3.1 17 2.391 3.083 3.123 3.124 2.387 cystine (+ lone pairs) 5.101 4.039 5.104 4.041 5.073 4.040 4.037 5.106 histidine D Continuum model of solvation (e = 78.3). In debye. Bonds involving heavy atoms only. Charged species; dipole moment calculated in the standard orientations3 framework of Gaussian 92.49 We also note that (ii) the resulting molecular dipole moment is increased accordingly, as expected from Onsager's original model,@and confirmed by other quantum chemical reaction field calculations, e.g. ref 7 1. However, modeling the solvent by a polarizable continuum has a major drawback in that the detailed structure of the model solvent is ignored, thereby preventing an accurate description of anisotropic intermolecular interactions, e.g. hydrogen bonding. As an attempt to understand the role of the first shell of solvation, an intermediate approach, consisting of a combined discrete-continuum model of was applied to the negatively charged aspartate side chain. The mono- and dihydrated complexes of the aspartate side chain, which, in the framework of our approximation,28 reduce to CH$20O-(H20), (n = 1,2), have been studied, in both vacuum and bulk aqueous solution (Le. e = 78.3). Their respective geometries were fully optimized in both cases,74,76using the extended split-valence type 6-31G** basis set. The starting point of these ab initio S C F geometry optimization was taken from the 1 i t e r a t ~ r e . Q ~ ~ In order to assess quantitatively the influence of the explicit water molecules surrounding the solute on the molecular charge distribution, additional computations on isolated water and on the model aspartate side chain were carried out a t the same level of approximation. It is worth noting that the effects of the bulk solvent, in the course of the geometry optimization, increases the SCF molecular dipole moments noticeably. This is due mainly to a slight decrease of the bond angles, together with the polarization of the electron cloud of the solute by the reaction field. Point charges derived from the electrostatic potential and from distributed multipole series are listed in Table V. The dipole moments calculated from these point charge models are presented in Table VI. As previously mentioned, the absolute value of the charges increases quite significantly when molecules are solvated by the dielectric medium. This is particularly evident for the charged species, which polarize the solvent significantly, so that the long-rangesolute-solvent interactions become the predominant contribution to the solvation term.42*69*84 It is interesting to observe that, in the case of water, DMM potential derived- and SCF electrostatic potential-derived net atomic charges are very similar (see Table V). This is obviously due to the fact that the low-order multipole moments are large enough to conceal the influence of site-site couplings.29 In other words, at each polar site, the representation of the potential is acceptable, and the chances to have contradictory partial fits are extremely limited. However, such fits generally yield high errors because atom-centered point charge models are inappropriate for the description of the electrostatic properties of molecules involving lone pair~.~4,27-30385.86Several approaches have been proposed to remedy this d e f i c i e n ~ y . ~The ~ - use ~ ~of~ additional ~~-~~ charges is economical, although not totally satisfactory.32 In the present work, we consider an extra charge located along the bisector of the water molecule,w+'1 below the oxygen (i.e. d ~ x 0.10 A), toward the hydrogens. This model offers a much better description of the potential, implying a substantial decrease

of the rmsd, and is in good agreement with the fact that the atomic dipole moments of the hydrogens are not exactly collinear with the 0-H bonds. From Table V, it may be observed that the inclusion of one or two water molecules with the model aspartate side chain modifies its charge distribution significantly. If only one water molecule is added, the conformation adopted by the resulting supermolecule is bifurcated83 (see Figure 4), H20 being coplanar with the COOgroup. The two oxygens of the carbonyl group are practically equivalent in both vacuum (Le. do,& = dolHs = 2.08 A and LC2O1H5N 104.7") and solution (Le. do,& = 2.03 A, dolHs= 1.92 A, and LC201HS 104.6"). As expected, a charge transfer occurs from the negatively charged ion toward the water molecule; i.e. the total charge of the model aspartate side chain decreases from -1 to -0.925 ecu, in vacuum, whereas in the bulk solvent, it decreases from -1 to -0.927 ecu, hence indicating that, beyond the first water molecule, the effects of the bulk solvent are relatively minor.81 In the case of the dihydrated complex, one of the water molecules remains in a quasi-bifurcated conformation, while the second one is located above the aspartate ion, roughly perpendicular to the carbonyl (i.e. in vacuum: dolb 1.84 A, do,^^ N 2.13 A, do,& N 2.08 A, LC201H6 N 124.3', LC202H4 N 106.0°, and ~C202HsN 105.2O; in aqueous solution: do1h N 1 . 8 4 A , d o l H s2.16A,do,g= ~ 1.92A,LC201Hs 121.2', LCzO2H4 N 102.9O, and LCzOzHs 104.7"). Potential-derived net atomic charges show that the contribution of the latter molecule to the global charge transfer is slightly more pronounced, i t . the total charge borne by (i) the bifurcated water molecule is -0.029 and -0.027 ecu, in vacuo and in solvent, respectively, (ii) the second water molecule is -0.063 and -0.041 ecu, in vacuo and in solvent, respectively. It should be pointed out that the low-entropy bifurcated (or quasi-bifurcated) carboxylate-water complexes seldom occur in Monte Carlo (MC) simulations.92~93This indicates that more water molecules should have been introduced explicitly in our combined discrete-continuum model of solvation in order to achieve better agreement with the hydrated structures obtained from statistical simulations. So far, we have learned that (i) a nonnegligible charge transfer occurs when the model aspartate side chain is complexed with a single water molecule, which preferentially adopts a bifurcated conformation, (ii) the presence of a second water molecule reinforces this effect to a muchlesser extent, and (iii) thedielectric medium representing the bulk solvent appears to damp the charge transfer, probably as a result of a shortening of the H-0 bond lengths. 4. Conclusions

In the two first parts of the present series of papers,28JO we have built point charge models of the naturally occurring amino acid. As a first approach, we derived net atomic charges from SCF molecular electrostatic potential and field, employing a constrained least-squares proced~re.1~ However, an efficient

The Journal of Physical Chemistry, Vol. 97, No. 38, 1993 9805

Modeling Amino Acid Side Chains

TABLE V: Influence of the Solvent*on the Potential-Derived (4")and Distributed Multipole-Derived (am) Net Point Charges' of Two Hydrated Complexes of the Aspartate Side Wins, Using the Split-Valence 6316** Basis Set; Inclusion of a Limited Number of Supplementary Sites 4v

atomic atomic and centers only virtual centersc qDMM in in in in in in atomd vacuo solution vacuo solution vacuo solution HzO, NPte= 1869 (vacuum), NPt = 1871 (solvent) 0 H X

-0.797 0.398

-0,872 0.436

0.786 0.571 -0.787 -0,861 0.549 0.574 0.394 0.430 -1.884 -1.718 0.195 0.188 0.784 0.722 4.62 3.50 61.32 21.28

rm~df(lO-~au)0.775 0.708 (96) 61.94 21.40 CH,COO-, NPt = 2447 (vacuum), Npt = 2451 (solvent) c1

HI H2

H3 C2 01

0 2 XClCl

-0.303 0.020 0.036 0.020 0.935 -0,854 -0,854

XC.0.

X&

rmsd(lW au) A (96)

0.325 0.22

-0,409 -0.628 0.073 0.099 0.094 0.112 0.073 0.099 0.993 1.822 -0,912 -0.718 -0,912 -0.718 -0.083 -0.492 -0.492 0.229 0.130 0.15 0.09

-0.587 0.119 0.141 0.119 2.077 -0.739 -0.739 -0.137 -0.627 -0.627 0.127 0.09

-0.124 -0.222 -0.035 0.014 -0.022 0.033 -0.035 0.014 1.053 1.094 -0.918 -0.966 -0.918 -0.966

0.681 0.46

0.770 0.52

CH&OO-,HzO, NPt = 3000 (vacuum), Nmt= 2975 (solvent) -0.462 -0.439 -0.626 -0.590 -0.086 -0.114 HI 0.081 0.093 0.102 0.125 -0.029 -0.005 H 2 0.079 0.093 0.100 0.123 -0.032 -0.006 H 3 0.092 0.117 0.112 0.147 -0.019 0.016 C2 0.874 0.884 1.252 1.505 1.057 1.057 01 -0.793 -0,835 -0.815 -0.821 -0.925 -0.947 0 2 -0.796 -0.840 -0.804 -0.837 -0.927 -0.954 H4 0.421 0.444 0.609 0.631 0.449 0.461 00 -0.919 -0.954 0.844 0.753 -0.936 -0.960 H5 0.422 0.438 0.610 0.601 0.449 0.452 XClCl 0.096 -0.063 XCA -0.172 -0.273 XCh -0.197 -0.267 XO, -2.109 -2.035 rmsd(lO-3au) 0.461 0.408 0.118 0.122 0.858 0.777 0.58 0.51 0.07 0.08 0.26 A 0.31

c1

CH3COO-,ZHzO Npt = 3560 (vacuum), NPt = 3533 (solvent) -0.388 -0.430 -0.576 -0.564 -0.100 -0,120 0.062 0.092 0.096 0.133 -0.024 -0,003 H2 0.090 0.081 0.107 0.127 -0.008 0.001 Ha 0.060 0.127 0.099 0.166 -0.013 0.027 C2 0.885 0.940 0.964 2.114 1.056 1.055 01 -0,838 -0.883 -0.937 -0.710 -0.957 -0.948 0 2 -0.780 -0,860 -0.808 -0.798 -0.894 -0.946 H4 0.433 0.499 0.627 0.640 0.452 0.494 0 3 -0,909 -1,019 1.006 0.335 -0,926 -0.996 0.447 0.468 0.593 HJ 0.413 0.479 0.625 0.474 0.516 0.609 0.642 H6 0.480 0.523 0 4 -0.901 -1.043 0.597 0.174 -0,886 -1.027 0.380 0.478 0.590 H7 0.392 0.494 0.510 XClCl -2.295 -0.258 XCIOl -1.750 -0,662 XCDl 0.118 -0.491 0.073 -1.598 XOl -0.069 -1.433 XO, rmsd(1P3au) 0.465 0.366 0.135 0.179 0.794 0.711 A (96) 0.34 0.25 0.10 0.13 0.61 0.55 CI HI

a Continuum mcdelof solvation(e = 78.3). In electronicchargeunits (ecu). H20 supplementary charge below the oxygen at do-x H 0.1 A, toward the hydrogens. CH3COW supplementary charges at the center of bonds involving heavy atoms only. d Full geometry optimization in vacuo and in solution, using the split-valence6-31G** basis set; indicts given in Figures 3 and 4. e Number of points requested for the fitting procedure. /Root-mean-squaredeviationbetween the ab initiocomputed and the point charge derived electrostatic potentials. Mean absolute deviations.

fitting of atom-centeredpoint charges to the potential is practically impossible, unless the molecule is highly polar and, as a result, the low-order contributions to the multipole expansion of the potential are predominant. This constitutesan intrinsicdeficiency of the method and forced us to reconsider the problem of point charge models from a different point of view. That distributed multipoleanalysis(DMA)32-35 is an interesting tool to investigatethe information contained in the wave function, offers a correct reproduction of the electrostatic potential outside the actual charge d i s t r i b ~ t i o n , ~and ~ - ~has ~ enabled the development of new methods to compute point charge models. A procedure, recently proposed by Ferenczy25 and which consists of fitting atomic charges to individual potentials generated by distributed multipole moments, has been investigated in the case of saturated hydrocarbon^^^ and model amino acid side chains.30 This method provides consistently transferable charges from one molecule to another structurally related one, together with an unfortunately poorer, but still reasonable, reproduction of electrostatic properties. In a parallel way, we have observed that the inclusion of a small number of adequately located off-atom sites leads to transferable ab initio SCF potential-derived point charge model~?~.M together with a better descriptionof electrostaticproperties. This may easily be explained by the fact that charges are not constrained any more on the atoms, so that the global distribution is more acceptable and the reproduction of the molecular multipole moments and, as a consequence, the electrostatic potential, is very satisfactory. Having obtained transferable point charge models, a key question remains: to what extent is the simplifiedmodel of amino acid sidechains acceptable? In other words, what is the magnitude of the perturbation induced by the presence of a backbone? In the first three examples presented in this work, Le. the extended conformations of Me-L-Ala-Me, Me-L-AspMe, and Me-MerMe, we have seen that the inclusion of a model backbone evidently modifies the charge distribution within the side chain but that our simplistic representation of the side chain is still acceptable, owing to the small magnitude of the alteration of the molecular charge distribution (especiallyif additional off-atom point charges are included). However, in our systematically simplified representation of the naturally occurring amino acid side chains,28 the charges on some locally symmetrical groups do not reflect the influence of the backbone environment adequately; e.g. the hydrogens are equivalent in the model alanine side chain (Le. C H 3 but not in the side chain of Me-L-Ala-Me (Le. CH3). In addition, the modification of the charges borne by the backbone in these three models, and in comparison with the reference "unperturbedwMe-Gly-Me, suggests that the side chain is at the origin of some conformational changes within the backbone. This is confirmed by a study of the molecular charge distribution in the extended and C7-equatorial conformers of NANMAsp, in which we have observed that point charges are altered by conformational modifications. This alteration may be partially counterbalanced by including additional sites in the fitting procedure. It is, however, obvious that the induction effects generated by the variation of the torsional angles cannot be taken into consideration by simple point charge models and that completely conformationally invariant charges requires the involvement of polarizabilities.54 Because most simulationsof oligo- and polypeptides are carried out in a solvent (typically water), the next logical question was, how is the molecular charge distribution of our model amino acid side chains influenced by the surroundings? As a first approach, we have considered the effects of a polarizable medium representing a bulk aqueous solution. In this model of solvation, the solute is placed in an ellipsoidal cavity surrounded by a continuum67 which is polarized by the solute. Its charge distribution is modified in return by a reaction field created by

The Journal of Physical Chemistry, Vol. 97, No. 38, 1993

9806

Chipot et al.

TABLE VI: Effects of the Solvent. on the Dipole Moment* of Two Hydrated Complexes of the Aspartate Side Chaii Using the Split-Valence6-316** Basis Set: Influence of the Inclusion of a Limited Number of Additional Sites ~

PLY

moleculed

HzO CH3COO- e CHaC00-,H2Oe CH3C00-,2H20C

atomic centers only in vacuo in solution 2.173 3.799 3.856 3.525

2.399 4.992 4.948 7.017

atomic and virtual centersC in vacuo in solution 2.147 3.794 3.832 3.523

2.376 4.987 4.929 7.012

PEP

PDMM

in vacuo

in solution

in vacuo

in solution

2.147 3.901 3.843 3.561

2.369 4.829 4.905 6.978

2.147 3.793 3.832 3.519

2.376 4.988 4.930 7.012

a Continuum model of solvation (e = 78.3). In debye. H20: supplementary charge below the oxygen at d G x N 0.1 A, toward the hydrogens. CH,COO-: supplementary charges at the center of bonds involving heavy atoms only. Full geometry optimization in vacuo and in solution, using the split-valence 6-31G** basis set. e Charged species; dipole moment calculated in the standard orientations3 framework of Gaussian 92.49

doctoral fellowship. We thank the Groupement Scientifique “Moddisation Molkulaire” IBM-CNRS for their generous support. The computations presented in this paper were carried out on the IBM-3090/600 andon the Serial IBM-RS/6000 cluster of the Cornel1 National Supercomputer Facility, Cornell University, Ithaca, NY, and on the IBM-3090/600 of the Centre Inter-Regional de Calcul Electronique, Orsay (France). This work was alsosupported by the U S . Nationalscience Foundation (Grants DMB-90-15815 and INT-91-15638) and the action incitative CNRS-NSF LET92 No. 21 MDRI.

References and Notes

(b) Figure4. Monohydrated (a) anddihydrated (b) complexesof theaspartate side chain, placed in the ellipsoidal cavity model of solvation.

the medium. This phenomenon has been observed in a series of five amino acid side chains and has led to the conclusion that point charges, as well as multipole moments, are affected by the presence of a solvent and that such effects should be considered when parametrizing thecharges for Monte Carlo (MC) or classical molecular dynamics (MD) simulations. A second approach, which consists of a mixed discretecontinuum model of has been applied for the study of two hydrated complexes of the negatively charged aspartate side chain. In the case of the monohydrate, the preferred position of the water molecule leads to a bifurcated structure. This is also observed with the dihydrate, in which the second water molecule is placed above the carboxylate group, perpendicular to it. In both supermolecules, a clear charge transfer occurs, slightly damped by the solvent effects. This strongly suggests that, at least, one water molecule should be taken into consideration explicitly for the molecular simulations of strongly hydrogenbonded complexesand that the utility of the dielectric continuum is limited to the bulk of solvent.

Acknowledgment. We thank Dr. G. G. Ferenczy for providing his program, Dr. F. Colonna for many stimulating discussions, and the continuous interest and encouragement of Prof. J.-L. Rivail. C.C. is indebted to Roussel Uclaf Co. (France) for his

(1) Momany, F. A,; McGuire, R. F.; Burgess, A. W.; Scheraga, H. A. J . Phys. Chem. 1975, 79, 2361. (2) Brooks, B. R.; Bmccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Compur. Chem. 1983,4, 187. (3) Hermans, J.; Berendsen, H. J. C.; van Gunsteren, W. F.; Postma, J. P. M. Biopolymers 1984, 23, 1513. (4) (a) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, Jr., S.; Weiner, P. J. Am. Chem. Soc. 1984, 106, 765. (b) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J . Comput. Chem. 1986, 7, 230. (5) Ltiwdin, P. 0. J. Chem. Phys. 1950, 18, 365. (6) Mulliken, R. S. J. Chem. Phys. 1955, 23, 1833. (7) Hinze, J.; Jaff6, H. H. J . Am. Chem. SOC.1%2,84, 540. (8) Bader, R. F. W.; Beddall,P. M.; Cade, P. E. J . Am. Chem.Soc. 1971, 93, 3095. (9) Gasteiger, J.; Marsili, M. Tetrahedron 1980, 36, 3219. (10) Cox, S. R.; Williams, D. E. J. Compur. Chem. 1981, 2, 304. (11) Guillen, M. D.; Gasteiger, J. Tetrahedron 1983, 39, 1331. (12) Nalewajski, R. F. J. Am. Chem. SOC.1984, 106, 944. (13) Mortier, W. J.; Genechten, K. V.;Gasteiger, J. J. Am. Chem. Soc. 198~,107,a29. (14) Chirlian, L. E.; Francl, M. M. J. Compur. Chem. 1987, 8, 894. (15) Williams, D. E.; Yan, J. M. Adu. At. Mol. Phys. 1987, 23, 87. (16) (a) Kim, S.; Jhon, M. S.; Scheraga, H. A. J. Phys. Chem. 1988,92, 7216. (b) Wee, S. S.; Kim, S.; Jhon, M. S.; Scheraga, H.A. J . Phys. Chem. 1990, 94, 1656. (17) Dinur, U.; Hagler, A. T. J . Chem. Phys. 1989, 91, 2949. (18) Baler, B. H.; Merz. Jr., K. M.; Kollman, P. A. J. Compur. Chem. 1990, 11, 431. (19) Woods, R. J.; Khalil, M.; Pell, W.; Moffat, S. H.; Smith, Jr., V. H. J . Compur. Chem. 1990,11, 297. (20) Breneman, C. M.; Wiberg, K. B. J . Comput. Chem. 1990,lI, 361. (21) Ferenczy, G. G.; Reynolds, C. A.; Richards, W. G. J . Compur. Chem. 1990, 11, 159. (22) Luque, F. J.; Illas, F.;Orozco, M. J. Compur. Chem. 1990,II, 416. (23) (a) No, K. T.; Grant, J. A.; Scheraga, H. A. J. Phys. Chem. 1990, 94,4732. (b) No, K. T.; Grant, J. A.; Jhon, M. S.; Scheraga, H. A. J. Phys. Chem. 1990, 94, 4740. (24) Colonna, F.; Angyh, J. G.; Tapia, 0. Chem. Phys. Lett. 1990,172, 55. (25) Ferenczy, G. G. J . Comput. Chem. 1991, 12, 913. (26) Aida, M.; Corongiu, G.; Clementi, E. Int. J. Quantum Chem. 1992, 42, 1353. (27) Colonna, F.; Evleth, E.; k g y h , J. G. J. Compur. Chem. 1992, 13, 1234. (28) Chipot, C.; Maigret, B.; Rivail, J.-L.;Scheraga, H. A.J. Phys. Chem. 1992, 96, 10276. (29) Chipot, C.; AngyBn, J. G.; Ferenczy, G. G.; Scheraga, H. A. J. Phys. Chem. 1993, 97, 6628. (30) Chipot, C.; hgydn, J..G.; Maigret, B.; Scheraga, H. A. J . Phys. Chem., preceding paper in this issue. (31) Merz, K. M. J . Compur. Chem. 1992, 13, 749. (32) Stone, A. J. Chem. Phys. Lett. 1981, 83, 233.

Modeling Amino Acid Side Chains (33) (a) Price, S.L.; Stone, A. J. Mol. Phys. 1982,47,1457. (b) Price, S. L.; Stone, A. J. Chem. Phys. Lett. 1983,98,419. (c) Price, S. L.; Stone, A. J. Mol. Phys. 1984,51,569.(d) Price, S.L. Chem. Phys. Lett. 1985,114, 359. (e) Price, S.L.; Stone, A. J.; Alderton, M. Mol. Phys. 1984,52, 987. (34) Stone, A. J.; Alderton, M. Mol. Phys. 1985,56, 1047. (35) Stone, A. J.; Price, S.L. J. Phys. Chem. 1988,92,3325. (36) (a) Hariharan, P. C.; Pople, J. A. Chem. Phys. Lett. 1972,16,217. (b) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.;Gordon, M. S.;DeFrees, D. J.; Pople, J. A. J . Chem. Phys. 1982, 77, 3654. (37) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio Molecular Orbital Theory; Wiley: New York, 1986. (38) Tapia, 0.; Goscinski, 0. Mol. Phys. 1975,29, 1653. (39) Rivail, J. L.; Rinaldi, D. Chem. Phys. 1976,18, 233. (40) Hylton McCreery, J.; Christoffersen, R. E.; Hall, G. C. J . Am. Chem. SOC.1976,98,7191. (41) (a) Miertus, S.;Scrocco, E.; Tomasi, J. Chem. Phys. 1981,55,117. (b) Miertus, S.;Tomasi, J. Chem. Phys. 1982.65, 239. (42) Tapia, 0.In Molecularlnreractior;Ratajczak, H., Orville-Thomas, W. J., Eds.; Wiley: New York, 1982;Vol. 3, p 47. (43) Karlstrbm, G. J . Phys. Chem. 1988, 92, 1318. (44) AngyBn, J. G. J. Math. Chem. 1992,10, 93. (45) Chipot, C. GRID (Fortran program performing charge fitting to molecular electrostatic potentials or fields;available from the authors), 1992. (46) Pascual-Ahuir, J. L.; Silla, E. J . Compur. Chem. 1990, 11, 1047. (47) Silla, E.;Villar, F.; Nilsson, 0.;Pascual-Ahuir, J. L.; Tapia, 0. J . Mol. Graphics 1990,8, 168. (48) Silla, E.; Tufion, I.; Pascual-Ahuir, J. L. J . Comput. Chem. 1991,12,

1077. (49) Frisch, M. J.; Trucks, G. W.; Head-Gordon, M.; Gill, P. M. W.;

Wong, M. W.; Foresman, J. B.; Johnson, B. G.; Schlegel, H. B.; Robb, M. B.; Replogle, E. S.; Gomperts, R.; Andres, J. L.; Raghavachari, K.; Binkley, J.S.;Gonzalez,C.;Martin,R.L.;Fox,D. J.;Defrees,D. J.;Baker,J.;Stewart, J. J. P.; Pople, J. A. GAUSSIAN 92;Gaussian Inc.: Pittsburgh, PA, 1992. (50). Ferenczy,G. G. MULFIT (Fortranprogramperformingchargefitting to multipole series; available from the author: Chemical Works of Gedeon Richter Ltd., H-1475Budapest, P.O. Box 27,Hungary), 1992. (51) (a) AngyAn, J. G.;Colonna, F.; Jansen, G. To be submitted to Chem. Phys. Lett. (b) Colonna, F. OSIPE (Open Structured Interfaceable Programming Environment): Dynamique des Interactions Molkulaires, Univexsite Pierre et Marie Curie. 4,Place Jussieu, Tour 22,75005Paris, France. (52) VignbMaeder, F.; Claverie, P. J . Chem. Phys. 1988,88, 4934. (53) Frisch, M.; Foresman, J.; Frisch, A. E. Gaussian 92User’s Guide and Gaussian 92Programmer’s Guide; Gaussian, Inc.: Pittsburgh, PA, 1992;pp

47,341. (54) Colonna, F.; Evleth, E. Submitted to Chem. Phys. Lett. (55) Bbhm, H. J.; Brcde, S. J. Am. Chem. SOC.1991, 113, 7129. (56) Nbmethy, G.; Gibson, K. D.; Palmer, K. A.; Yoon, C. N.; Paterlini, G.; Zagari, A.; Rumsey, S.;Scheraga, H. A. J. Phys. Chem. 1992,96,6472. (57) SchHfer, L.; Newton, S.Q.; Cao, M.; Peeters, A.; Van Alsenoy, C.; Wolinski, K.; Momany, F. A. J . Am. Chem. Soc. 1993,115,272. (58) Pulay, P.; Fogarasi, G.; Pang, F.; Boggs, J. E. J . Am. Chem. SOC.

1979,101,2550. (59) Jorgensen, W. L.; Gao, J. J . Am. Chem. SOC.1988,110,4212. (60) Reynolds, C. A,; Essex, J. W.; Richards, W. G. J . Am. Chem. SOC. 1992,114,9075. (61) Stouch, T. R.; Williams, D. E. J . Comput. Chem. 1992,13,622. (62) Urban, J. J.; Famini, G. R. J . Comput. Chem. 1993,14,353.

The Journal of Physical Chemistry, Vol. 97, No. 38, 1993 9807 (63) Rinaldi, D.; Rivail, J. L. Theor. Chim. Acta 1973,32,57. (64) Rivail, J. L.; Terryn, B. J . Chim. Phys. Phys.-Chim. Biol. 1982,79, 1.

(65) Rinaldi, D.;Ruiz-Lopez, M. F.; Rivail, J. L. J. Chem. Phys. 1983, 78. 834. (66) Rivail, J. L.; Terryn, B.; Rinaldi, D.; Ruiz-Lopez, M. F. J. Mol. Struct.: THEOCHEM 1985,120,387. (67) Kirkwood. J . G. J . Chem. Phvs. 1934,2,351. (68) Onsager, L. J. Am. Chem. S&. 1936,58, 1486. (69) (a) Rashin, A. A.; Honig, B. J . Phys. Chem. 1985, 89,5588. (b) Rashin, A. A.; Namboodiri, K. J . Phys. Chem. 1987,91, 6003. (70) (a) Sharp, K. A.; Honig, B. Annu. Rev. Biophys. Biophys. Chem. 1990,19,301.(b) Sharp, K.A,; Honig, B. J. Phys. Chem. 1990,947684. (71) Grant, J. A.; Williams, R. L.; Scheraga, H. A. Biopolymers 1990, 30,929. (72) Vorobjev, Y. N.; Grant, J. A,; Scheraga, H. A. J. Am. Chem. Soc. 1992,114,3189. (73) Rinaldi, D. Comput. Chem. 1982,6,155. (74) Rinaldi, D.; Rivail, J. L.; Rguini, N. J. Comput. Chem. 1992,13, 675. (75) Frisch, M. J.; Head-Gordon, M.; Schlegel, H. B.; Raghavachari, K.; Binkley, J. S.;Gonzalez, C.; Defrees, D. J.; Fox, D. J.; Whiteside, R. A.; Seeger, R.; Melius, C. F.; Baker, J.; Martin, R. L.; Kahn, L. R.; Steward, J. J. P.; Fluder, E. M.; Topiol, S.;Pople, J. A. GAUSSIAN 88; Gaussian Inc.: Pittsburgh, PA, 1988. (76) Rinaldi, D.; Pappalardo, R. R. QCPE No. 622;SCRFPAC: a selfconsistent reaction field package. (77) Olivares del Valle, F. J.; Tomasi, J. Chem. Phys. 1991, 150, 139. (78) Chipot, C.; Rinaldi, D.; Rivail, J. L. Chem. Phys. Lett. 1992,191,

287. (79) (a) Beveridge, D. L.; Schnuelle, G. W. J. Phys. Chem. 1974, 78, 2064. (b) Schnuelle, G. W.; Beveridge,D. L. J. Phys. Chem. 1975,79,2566. (80) Claverie, P.; Daudey, J. P.; Langlet, J.; Pullman, B.; Piazzola, D.; Huron, M. J. J. Phys. Chem. 1978,82,405. (81) Hodes, Z.I.; Nkmethy, G.; Scheraga, H. A. Biopolymers 1979,18, 1565. (82) Thanki, N.; Thornton, J. M.; Goodfellow, J. M. J. Mol. Biol. 1988, 202, 637. (83) Jinfeng, C.; Topom, R. D. J . Mol. Struct.: THEOCHEM 1989, 188,45. (84) Born, M. Z . Phys. 1920,I, 45. (85) Williams, D. E. J . Comput. Chem. 1988,9,745. (86) Williams, D. E. Biopolymers 1990,29, 1367. (87) Bernal, J. D.; Fowler, R. H. J . Chem. Phys. 1933,I, 515. (88) Stillinger, F. H.; Rahman, A. J. J. Chem. Phys. 1974, 60, 1545. (89) Berendsen,H. J.C.;Postma, J.P. M.;VanGunsteren, W.F.;Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrtcht. 1981;p 331.

(90) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79,926. (91) (a) Lie, G. C.; Clementi, E.; Yoshimine, M. J. Chem. Phys. 1976, Clementi, E.; Yoshimine, M. J . Chem. Phys. 64,2314. (b) Matsuoka, 0.; 1976,64,1351. (92) Jorgensen, W. L.; Gao, J. J. Phys. Chem. 1986, 90,2174. (93) Alagona, G.; Ghio, C.; Kollman, P. J. Am. Chem. Soc. 1986,108, 185.