Molecular Mechanical and Molecular Dynamic Simulations of

J. Phys. Chem. 1995,99, 3832-3846

3832

Molecular Mechanical and Molecular Dynamical Simulations of Glycoproteins and Oligosaccharides. 1. GLYCAM-93 Parameter Development Robert J. Woods,* Raymond A. Dwek, and Christopher J. Edge Glycobiology Institute, Department of Biochemistry, The University of Oxford, South Parks Road, Oxford, OX1 3QU, U.K.

Bert Fraser-Reid The Department of Chemistry, Paul M. Gross Chemical Laboratory, Duke University, North Carolina 27706 Received: August 15, 1994@

A new parameter set (GLYCAM-93) that is consistent with the all-atom AMBER force field has been developed for molecular dynamics simulations of glycoproteins and oligosaccharides. Torsional energy terms were derived from ab initio molecular orbital calculations performed at the restricted Hartree-Fock (RHF) level with the split valence 6-31G* basis set on a series of moderately large representative “fragments”, each a derivative of tetrahydropyran. High-level ab initio calculations have been performed that indicate that a gauche preference in 1,2dmethoxyethane, of particular relevance to 1-6 linked sugars, does not exist in the gas phase. Net atomic charges that reproduce molecular electrostatic potentials (ESPs) at the HF/6-31G* level were employed and have been determined for each relevant monosaccharide. These ESP-charges are consistent with the TIP3P water model. Two approaches that minimize the problem of the conformational dependence of partial atomic charges are discussed. The new parameters correctly predicted the subtle effects on the rotational properties of the glycosidic linkages in models for methyl a-D-gluco- and a-Dmannopyranoside.

Introduction The significant alterations in the biological activities of peptides and proteins that sometimes accompany the covalent attachment of a carbohydrate (glycosylation) to one or more of the amino acid residues are becoming increasingly The carbohydrate component of a glycoprotein may be present in a largely structural capacity (as in mucins, where it causes a 2-fold extension in the length of the p r ~ t e i n ) ,or~ it may alter the functioning of the protein (as in human chorionic gonadotropin, hCG, tissue plasminogen activator, t-PA, and ribonuclease B),“-6 or it may be recognized by receptors on other proteins, such as by the carbohydrate recognition domain (CRD) proteins.’ It is therefore essential that the spatial and dynamic properties of these molecules be accurately determined. Experimental techniques, such as nuclear magnetic resonance (NMR)8-1 and fluorescence energy-transfer~pectroscopy~~ have been valuable aids in structural studies of oligosaccharides and glycoproteins. While NMR spectroscopy has been widely applied to oligosaccharides, the conformations of the glycosidic linkages in these systems are particularly difficult to determine. This difficulty arises largely from two factors, namely, the characteristic paucity of nuclear Overhauser effects (nOes)13and the uncertain accuracy of Karplus-type relationships derived from heteronuclear J-coupling constants.l4 An alternative approach to the problem of glycoprotein and oligosaccharide conformational analysis is provided in principle by molecular dynamics (MD) simulations. That there is considerable interest in applying molecular mechanical and dynamical methods to oligosaccharides is illustrated by the numerous reports of force fields and parameter sets for *To whom correspondence should be addressed. Present address: Complex Carbohydrate Research Center, University of Georgia, 220 Riverbend Rd., Athens, GA 30602. Abstract published in Advance ACS Abstracts, January 15, 1995. @

carbohydrates. Among the force fields employed to date have been MM2,15 MM3,16-18 CHARMM,19 AMBER,2°,21GROMOS,22 HSEA,23and TRIPOS.24925However, none of these methods has as yet gained widespread acceptance for application to MD simulations of oligosaccharides. Moreover, at least one recent report suggests that no current method is suitable for fully solvated MD simulations performed without experimental restraints.26 One of the limitations common to each of the current carbohydrate parametrizations has been the lack of accurate torsional energy profiles. As we wished not only to determine the structures of oligosaccharides but also to examine their dynamic properties and eventually to apply the technique of free energy perturbation to the study of oligosaccharideprotein interactions, we undertook the development of a suitable set of carbohydrate parameters for use with the AMBER force field for proteins. We have recently applied these parameters in preliminary molecular dynamics simulations of oligosaccharides and report here full details of the parameter development .27-29 Regardless of the particular mathematical formulation, the force field must be accurately parametrized; the parameters must have been derived in a consistent manner for each species of interest. The ultimate aim in molecular modeling is both to explain experimental observations and to act in a predictive capacity. Thus, it is essential that the force field and parameter set should predict correctly the known experimental results. For a thorough discussion of the problems associated with forcefield development, see ref 30. In the following sections we derive and evaluate a consistent parameter set for application in the molecular modeling of oligosaccharidesand glycoproteins that is appropriate for use with the T P 3 P model for water and that is thoroughly integrated with the all-atom AMBER force field for proteins and nucleic acids. We selected the AMBER force field as the framework for our parameters for two reasons. First, AMBER has been

0022-3654/95/2099-3832$09.00/0 0 1995 American Chemical Society

MD Parameters for Glycoproteins and Oligosaccharides extensively employed in protein conformational analysis and in free energy perturbation calculations, and its strengths and weaknesses are well d ~ c u m e n t e d . ~ lSecond, - ~ ~ we believed that the presence of asymmetric phases in the cosinusoidal expansion of the torsional potentials in AMBER would facilitate an accurate description of the inherently asymmetric torsional energy profiles of the glycosidic linkages.

J. Phys. Chem., Vol. 99, No. 11, I995 3833 CHART 1 40-Me 3

1

OMe 1

2

Method We have performed ab initio molecular orbital calculations, using the G a ~ s s i a n 9 0and ~ ~ 9238 series of software packages, to estimate the energies associated with the rotation about the glycosidic (Cl-01) and C5-C6 bonds in a series of model glycosides derived from tetrahydropyran. The geometries of these molecules were fully optimized at the restricted HartreeFock (RHF) level with the 6-31G* basis set,39unless otherwise noted. In the determination of rotational energy profiles the relevant torsion angles were varied from 0-330' in 30' increments while all other structural variables were optimized. Ab initio calculations were performed on either a Cray-YMP or a Convex C210 computer. Semiempirical molecular orbital calculations were performed using the AM1@ Hamiltonian as present in a version of MOPAC 6.041modified by us to execute on an FPS P64/30 computer. All molecular mechanical calculations were performed with AMBER 4.P2on either a VAX 4000-200 or an FPS P64/30 computer. In general, the strategy for parameters derivation followed that reported by Hopfinger and Pearl~tein.4~Following this method the geometric parameters of interest were distorted from their equilibrium values to generate distortion-energy curves. Using a nonlinear least-squares curve-fitting algorithm and incorporating the relevant analytical forms of the energy functions, force constants were derived for the HF/6-31G* distortion curves. Torsional parameters were derived after all other parameters had been estimated by fitting the AMBER analytical form of the torsion energy function to the difference between the ab initio energies and the AMBER energies.

Discussion Glycoside Model Selection. To avoid an unnecessarily large or redundant parameter set, we sought to ensure parameter generality and transferability between sugars. This was achieved through the selection of model compounds that exhibited only the essential structural characteristics of pyranosides. To minimize the possibility of inaccurately modeling the rotational barriers associated with the glycosidic linkages in oligosaccharides, we chose a set of representative molecules all of which were analogs of tetrahydropyran (see Chart 1). While an earlier parametrization of AMBER for carbohydrates employed dimethoxymethane as a model for the glycosidic linkages,44we selected the larger structures, which reflect more accurately the torsional properties of pyranosides. In the discussion that follows, we have chosen to adopt the numbering scheme for carbohydrates rather than that for derivatives of tetrahydropyran. Basis Set Selection It was essential to select a basis set that could be expected to reproduce accurately the geometries and energies of carbohydrates. However, the comparatively large size of our model compounds imposed a constraint on the complexity of the basis set that we could reasonably expect to use in routine calculations. Previously we have reported the results of an examination of the abilities of several molecular orbital methods to predict the energies and geometries of conformations of dihydroxymethane, the simplest model for the 0-C-0 atomic sequence unique to

C=9-"Ac

6 OMe

4

3

4 OMe OMe

q

OMe OMe

6

5

7

8

carbohydrate^.^^ In that work we found that calculations performed at the RHF level with a split valence basis set with one set of polarization functions on the non-hydrogen atoms (6-3 1G*) gave results for relative conformational energies that were within 0.20 kcal mol-' of the corresponding values generated at the second order M~ller-Plesset46-~*(MP2) level with a significantly more sophisticated basis set that included polarization and diffuse functions on all atoms (6-311++G**).49 Moreover, the results of that study indicated that the relative changes in C - 0 bond length that are a recognized feature of carbohydrate conformation^^^-^^ were reproduced fully with the 6-31G* basis set. Consequently, we employed the 6-31G* basis set in all of the ab initio calculations. unless otherwise stated. Bond Length and Angle Parameters The well known variations in the C1-01 and C1-05 bond lengths with configuration at C1 in pyranose^^^-^* could be fully described by a stretch-torsion cross term in the force field. Cross terms may be necessary when a high level of accuracy in reproducing vibrational frequencies is desired53and have been incorporated into type II force such as MM3l6-I8 and CFF91.55*56However, a stretch-torsion term does not exist in any of the well-established macromolecular force fields, such as AMBER, CHARMM, or GROMOS. In the case of carbohydrates, an alternative approach for incorporating these variations is to define atom types that are specific for each anomer. Unique values for the force constants and equilibrium bond angles may then be specified for the a- and p-anomers. This method has the advantage of being mathematically simple and so is computationally rapid. However, it lacks the ability to describe the variations in geometry that may occur during

3834 J. Phys. Chem., Vol. 99, No. 11, 1995 chemical processes, such as glycoside hydrolysis, that involve large changes in ring c ~ n f o r m a t i o n .We ~ ~ have introduced only three new atom types, namely, a unique anomeric carbon atom for each anomer (labeled AC and EC) and a glycosidic oxygen atom (OG). The labels AC and EC denote the presence of an axial or equatorial substituent heteroatom at C1. We did not believe it was necessary to create new atom types for all of the atoms in a sugar and sought to introduce as few new parameters as possible. With the three new atom types we were able to characterize uniquely the bonds ( R ) and angles (8)involving C1, 0 1 , and 0 5 in both anomers, in terms of their equilibrium values ( R , and 8,) and stretching and bending force constants (Ksand Kb). The van der Waals parameters for the new atom types AC, EC, and OG were assumed to be the same as those employed in the current AMBER parameter set (PARM91) for the corresponding tetrahedral carbon atom (CT) and for an ethertype oxygen atom (OS). Not only do the glycosidic bond lengths and angles vary as a function of anomeric configuration but also as a function of the size of the aglycon.58 Although our rotational energies were to be based on the representative pyranoid molecules containing a methyl group as aglycon, our goal was to apply the parameters in the modeling of oligosaccharides. Consequently, values for R , and 8, at the anomeric center were selected either from a search of the Cambridge Structural Database (CSD)59 for 0-linked oligosaccharides by Schleifer and co-workers5*or, for N-linkages, from an analogous search of the CSD performed by us. Force constants for many of the bonds and angles in which we were interested have been published by Ha and co-workers for application in a CHARMM-type force field for carbohydrates.60 However, it is well-known that parameters are not generally transferable between force fields; particularly critical are the nonbonded and torsion parametem61 While such an amalgamation of parameters has been employed in other force fields for oligosaccharide^,^^^^^^^ we chose to derive parameters that were specific for use with the PARM91 parameter. Set as our interests are primarily directed toward macromolecular dynamics, we have not optimized the stretching force constants for the glycosidic bonds. Rather, it was assumed that the values for bonds involving the new atom types (AC, EC, or OG) would be approximately equal to those for the corresponding bonds involving atom types CT and OS in PARM91. Similarly, with the exception of the characteristic 0-C-0 and 0-C-N angles, the bending force constants in PARM91 were employed. Flexibility of the glycosidic angles could have a significant effect on the overall dynamics of an oligosaccharide and an accurate estimation of the bending force constants for these angles was required. Values for the bending force constants for the 0-C-0 and N-C-0 angles were estimated from distortion analyses of the relevant angles in dimethoxymethane (DMM) and (N-formy1amino)methoxymethane (FAMM). h e vious computational studies aimed at quantifying the anomeric effect in N-C-0 containing species have examined the potential energy surfaces of amino- or (N-methylamino)meth~xymethane.~f'~ As the N-linkage in glycoproteins always occurs through an amido nitrogen atom, we believed the N-formylated analogue would be a more appropriate model for this linkage. Fitting to the HF/6-31G* energies gave rise to preliminary values for K b for the 0-C-0 and N-C-0 angles of 123.0 and 118.8 kcal rad-2, respectively (see Figure 1). However, force constants derived from ab initio molecular orbital calculations have been reported to be 10-15% higher than those derived experimentally66and our values for Kb were scaled by a factor of 0.9, as recommended in ref 66. The scaled value of Kb for the 0-c-0 angle (110.7 kcal rd-*) is slightly

Woods et al.

+HFl6-31G' ---b HF16-310' +GLYCAM-Q3 +GLYCAM-93

'"1 0.8 0.7

Relative Energy (kcal/mol)

P

A

__ 104

106

108

170

114

112

116

118

120

0-C-X (X=N,O) Bond Angle (deg)

iv\

Figure 1. Bond distortion energy for dimethoxymethane (0)and (N-

forrny1amino)methoxymethane(A).

'i Fi 6

(kcaVmoi) Relative Energy

--t HFl6-310'

2

4 GLYCAM-93

1

0 4 . I . 0

30

I

60

,

,

,

,

I

.

I

,

,

,

!

'

I

.=!. 1

.

/

90 120 150 180 210 240 270 300 330 360

05-C1-01-CTorsion Angle (deg) Figure 2. N-Acetyl group rotation in 3. The total energy for the minimum energy conformation at the HF/6-31G* level was -476.842 3% au

.

larger than that in either Jeffrey and Taylor's parametrization (80.6 kcal rad-2)67 of the MM2 force field15 or the parametrization of the C W M force field by Ha and co-workers (92.6 kcal and considerably larger than the values employed by Scarsdale and co-workers in their parametrization of AMBER (40.6 kcal rad-2)63or in GROMOS (68 kcal rad-2).68 The value for 8, for the 0-C-0 angle in DMM was calculated to be 113.8'. This value is larger than that observed in a-glycosides, by approximately 2°.58 Consequently, for application in pyranoid molecules the 8, values for the 0-C-0 angles in axial (1 12.0') and equatorial (107.6') glycosides were chosen based on the average of the values reported for The e,, value (111.1") for the 0-C-N angle in (N-formylamino)methoxymethane was also found to be larger than that observed in N-glycosides, by approximately 3°.69,70Therefore, the value of 0, in N-linked glycosides was selected from the average of the values reported for related pyranoside~.6~-'~ The differences between the 8, values in the small model compounds and those observed in pyranosides presumably reflects the influence of cyclization on these angles.

Torsional Parameters To incorporate torsion terms into the force field that were specific for the pglycosidic linkage present in 0- and Nglycosides (05-Cl-X-C,~,,,,, X = O, N) we chose to introduce three torsion parameters, namely, for the a-0-, /I-0-, and /I-N-linkages (see Figures 2-4). There has been only one report of an a-N-linked glycoside, and we will report suitable parameters for that linkage at a later date." A close agreement between the ab initio energies and the GLYCAM-93 energies was required to ensure that the parameters would reproduce both the conformation and dynamics of the carbohydrate-carbohy-

MD Parameters for Glycoproteins and Oligosaccharides

3 O.,,, CH3

12 1

Relative Energy (kcaVmol)

:h 5 4 3

--f

2 1 "

-

I

0

\

-e-

HFl6-310' GLYCAM-93

\

30 60 90 120 150 180 210 240 270 300 330 360

05-C1-01-C Torsion Angle (deg)

Figure 3. Methoxyl group rotation in 1. The total energy for the minimum energy conformation at the HF/6-31G* level was -383.907 631 au.

Relative Energy (kcaVmol)

13 12 11 10 9 8 7 6 5 4 3

+HFl6-310'

+GLYCAM-93

0

30 60 90 120 150 180 210 240 270 3 0 330 360

OB-C1-N-C Torsion Angle (deg) Figure 4. Methoxyl group rotation in 2. The total energy for the minimum energy conformation at the HF/6-31G* level was -383.909 932 au.

drate and carbohydrate-protein linkages. However, to maintain computational efficiency, the torsional expansions were restricted to a maximum of three terms. The energy profiles for rotation about the 0-C-0-C torsion angle in 1, a model for B-D-glycosides in the 4clconformation, calculated at the HF/6-31G* and AM1 levels shows two minima corresponding to Q, angles of approximately 60" and -60' with the latter preferred energetically over the former, by approximately 2.9 kcal mol-' in the case of the HF/6-31G* results (see Figure 2). Each of these orientations is expected to contain favorable exo-anomeric interactions. The exo-anomeric effect reflects to the preference for a gauche orientation of the aglycon 0-R bond with respect to the endocyclic C - 0 bond.72 The energetic difference between the two orientations presumably arose from the presence of repulsive steric interactions in the 60" orientation between the aglycon and the hydrogen atoms attached to C2. It is notable that no minimum was predicted in the region of 180" in contrast to the results obtained from studies employing the PCILO semiempirical molecular orbital method.73 Furthermore, the PCILO calculations predicted the maximum barrier to rotation to be approximately 4 kcal mol-', whereas our ab initio data predicted a value of approximately 8 kcal mol-'. A close similarity between the values calculated with full geometry optimization at the ab initio level (HF/63 1G*//6-31G*) and with those obtained with the same ab initio basis set at the AM1-optimized geometry (HFI6-31G*//AMl)

J. Phys. Chem., Vol. 99, No. 11, 1995 3835

was observed. This similarity is particularly intriguing as the AM1 energies for the AM1 geometries (AMWAh41) appeared to predict a rather different rotational profile than did the ab initio data. The AM1 rotational energy profile for 1 has been previously reported and ~ r i t i c i z e d .Three ~ ~ torsion coefficients (VI = 0.85, V2 = 0.74, V3 = 0.97, kcal mol-') and asymmetric phase corrections ( y l = 144.06", y2 = 353.55", y3 = 6.43') were employed to generate the GLYCAM-93 energies. For application to L sugars or sugars that exist in the lC4 conformation the phases for this torsion term must be reflected through the phase origin, while the coefficients remain unchanged ( y l = 215.94", y2 = 6.45', y3 = 353.57"). The energy profiles for rotation about the a-glycosidic torsion angle in 2 is presented in Figure 3. As in 1, a surprisingly good agreement was seen between the HF/6-31G*//6-31G* and HF/6-31G*//AM1 results. The ab initio results indicated the presence of only one minimum, at an 0-C-0-C torsion angle of approximately 60', corresponding to the preferred exoanomeric orientation. A comparison of the ab initio data with the results from PCILO calculations, in which minima at approximately 60" and 180" were reported, again indicated significant differences between the ab initio and the PCILO semiempirical results. However, in this case the maximum barrier to rotation predicted at the PCILO level (approximately 9 kcal mol-') was in reasonable agreement with that predicted at the HF//6-31G* level. In contrast, the AM1 Hamiltonian predicted the presence of two minima, at approximately 60" and -60°.74 Neither the results from the PCILO nor those from the AM1 Hamiltonians agreed with the ab initio results presented here, nor did they agree with each other. This observation is in contrast to conclusions based on studies of the model compound DMM, for which AM1 accurately reproduced the trends in the conformational energies and g e o m e t r i e ~ .As ~ ~in 1, a torsional expansion into three components was sufficient to reproduce the ab initio energy profiles (VI = 1.39, V2 = 0.70, V3 = 0.91, kcal mol-l, yl = 265.77', y2 = 312.04", y3 = 347.72'). For application to L sugars or sugars that exist in the 'C4 conformation the phases for this torsion term are y1 = 94.23", y2 = 47.96", y3 = 12.28'. An estimate of the energy of the anomeric effect may be deduced from the difference between the minimum energy conformations of 1 and 2. It should be noted that no attempt has been made to include entropic contributions to this en erg^.^^^^^ The anomeric effect refers to the tendency of an electronegative substituent at C1 of a pyranoid ring to assume the axial rather than equatorial orientation, in contrast to predictions based solely on steric g r o ~ n d s . ~We ~ , note ~ ~ ,that ~~ in our geometry optimizations of 1 and 2 the 0-C-0-C torsion angle was constrained, thus our value for the anomeric effect (1.44 kcalmol-') does not represent the exact HF/6-31G*/ /6-31G* value. It is, however, correct in sign,79 and in good agreement with both a previous ab initio result (1.33 kcal mol-') calculated at the HF/6-3 1G*//3-21G level,g0and experimental values determined in relatively nonpolar solutions (0.75 kcal mol-' 81 in CFC13/CDCl3, and 0.89 kcal mol-' 82 in CC4 and Cd-I6). Experimental and theoretical aspects of the anomeric effect have been reviewed recently.83 The computed SCF-dipole moments for the exo-anomeric conformations of 1 and 2 (1.8 and 0.4 D, respectively) are in reasonable agreement with the experimental values (2.1 and 0.8 D, respectively).82 In the case of both 1 and 2 the dipole moments of the 180" rotamers (2.8 and 2.2 D, respectively) are larger than those of the minimum energy conformers, in agreement with experiment (2.9 and 2.1 D, respectively), and indicated that these conformations would be preferentially stabilized in polar solvents. This result is supported by the

3836 J. Phys. Chem., Vol. 99, No. 11, 1995

i

Woods et al.

R-b? HH

I COW

RCH3,H

12 11 10 9

+R=Me, GLYCAM-93 --f R=Me, HFI6-31W -t- R=Me, MP2/531G'

8

0

4

A 0

-+-R=Me, AM1

Relative 7 Energy 6 (kcal/moi) 5

2t \

0-R

30 80 90 120 150 160 210 240 270 300 330 380

C1-CP-N-C Torsion Angle (deg)

Figure 5. N-Acetyl group rotation in 7. The total energy for the minimum energy confoxmation at the HF/631G* level was -476.834 990 au.

3 2 1 0 0

30 80 90 120 150 180 210 240 270 300 330 360

0-C-C-0 Torsion Angle (deg)

Figure 6. Rotation of the 0-C-C-0 torsion angle in DME and DHE (ZM~(H)-O-C-C = 180'). The total energy for the minimum energy conformation at the HF/6-31G*, MP2/6-31G*, and MP2/6-311++G** levels for DME were -306.980 196, -307.852 136, and -308.074 533 au, respectively. The heat of formation for the minimum energy conformation of DME calculated by AM1 was -99.449 kcal mol-'. For DHE the total energy for the minimum energy conformation at the HF/6-31G* level was -228.922 424 au.

observation that both the -60' and 180' conformers of 1 are present in solutions.84 On the basis of our results, we would also anticipate the presence of the 180' conformer of 2 in aqueous solution. Thus, the plateau regions in the gas-phase rotational profiles for 1and 2 (approximately 150") may become minima when the sugars are dissolved in polar solvents. Presented in Figure 4 is the rotational energy profile for the agreement with s o l i d - ~ t a t and e ~ ~recent ~ ~ ~NMR data in aqueous 0-C-N-C angle in 3 calculated at the HF/6-31G* level and solution,86 the carbonyl group at C2 was predicted nearly to the corresponding profile calculated with the GLYCAM-93 eclipse the proton at C2. However, unlike the case of the parameters (VI = 2.39, V2 = 1.40, V3 = -0.41, kcal mol-', y1 N-linkage, in which a high barrier to rotation in the 0" conformer = 200.31", y2 = 359.33", y3 = 175.02"). As this linkage is was found, the rotational barriers for the C2-N2 bond were critical to the conformation of N-glycosides, the ab initio calculated to be 4.34 and 4.83 kcal mol-', at 0 and 240°, calculations were performed without constraining the amide respectively. nitrogen atom to a planar geometry. This was found to be The final torsion angle that was examined was the O-Cextremely computationally intensive due to the occurrence of C-0 angle, which is relevant to 1 4 linked sugars. Parameters frequent inversions at the nitrogen atom during energy minifor this sequence exist in PARh491; however, they were derived mization. In agreement with 'H NMR studies,85at the minimum under the explicit assumption that a gauche rotamer (to-c-c-0 energy conformation the carbonyl group of the N-acetamido = f60.0') would be more stable than an anti, in the gas phase.20 moiety approximately eclipses the anomeric proton. The value A theoretical explanation for a gauche preference in atomic of the 0-C-N-C torsion angle in the minimum energy sequences of the type X-C-C-Y (X, Y = electronegative conformation (-90") is in good agreement with data for the atoms) has been pre~ented.'~*~'While it is known that in solid-state conformations of 4-N-(-#?-~-glucopyranosyl)-~-as- solution the hydroxymethyl group prefers a gauche orientation paragine (-83.6'),'O 4-N-(2-acetamido-2-deoxy-#?-~-glucopy- with respect to the ring oxygen atom,88-% an anti orientation ranosy1)-L-asparagine (- 100.2),7Oand 1-N-acetyl-2-acetamidohas been predicted to be preferred in the gas phase.91*92As we ~-D-g~ucopyranosy~amine (- 117.1°).69 A second stable wished to treat solvent explicitly, we sought to determine conformation was predicted to exist, in which the 0-C-N-C appropriate torsion terms for the 0-C-C-0 sequence by torsion torsion angle was approximately 60'. This conformer performing high-level ab initio calculations on models for the was less stable than the former, by approximately 3.5 kcal 1-6 linkage, namely, 1,2-dimethoxyethane (DME) and 1,2mol-'. The comparatively large barrier to rotation (12.5 kcal dihydroxyethane (DHE). Experimental evidence for the presmol-') that occurred when the carbonyl group eclipsed the ring ence of gauche conformers of 1,2-dihydroxyethane (DHE) in oxygen atom presumably arose from unfavorable Coulombic the gas phaseg3and 1,Zdimethoxyethane (DME) both neat and interactions between the carbonyl and ring oxygen atoms. In in chloroform solution has been r e p ~ r t e d ? ~The - ~ ~ gauche the presence of an N-acetyl group at C2, as found in N-linked conformation of the former compound may be stabilized by the glycosides, the barrier to rotation at 120" would be likely to presence of an intramolecular hydrogen b0nd.9~ Definitive gasincrease significantly. On the basis of these data we would phase data for DME have not been reported. To allow a direct anticipate that the 0-C-N-C torsion angle in N-linked comparison to be made between DHE and DME, the hydroxyl glycosides would have only very limited conformational flexand methoxyl groups were constrained to anti orientations, thus ibility. precluding the formation of internal hydrogen bonds in the Related to the 05-C1-N1-C glycosidic linkage is the C1former compound. The results of calculations at the W/6-3 lG* C2-N2-C linkage present in 2-deoxy-2-(N-acetylgluco)-and level indicated that the anti rotamers were preferred over the 2-deoxy-2-(N-acetylgalacto)pyranosides.For the sake of comgauche, by approximately 1.9 kcal mol-' for both DHE and putational efficiency, the energy profile for rotation of the DME (see Figure 6). Previous calculations have indicated that N-acetyl group in 7 was determined with the amide nitrogen both 2-chloro- and 2-fluoroethanol prefer gauche orientations atom constrained to planarity. The HF/6-31G* and only when an intramolecular hydrogen bond is present.% In GLYCAM-93 results are presented in Figure 5 . Analogous rotamers without a hydrogen bond the anti rotamer is again to the results for the 0-C-N-C glycosidic linkage, in the favored over the gauche, by approximately 1.29 and 1.13 kcal minimum energy conformation (n-1-C2-N2-C = 90.0°>, and in mol-I, for 2-chloro- and 2-fluoroethanol, r e s p e ~ t i v e l y .To ~~

MD Parameters for Glycoproteins and Oligosaccharides Me

--t HFI6-31G' -8AM1

7

Relative Energy (kcaVmol)

-+- GLYCAM-93

5

3 2

1 0

0

30

80 90 120 150 180 210 240 270 300 330 360

06-C6-C5-05 Torsion Angle (deg) Figure 7. Methoxymethyl group rotation in 4. The total energy for the minimum energy conformation at the HF/6-31* level was -535.823 267 au. The heat of formation for the minimum energy conformation calculated by Ah41 was -154.514 kcal mol-'.

examine the extent to which electron correlation altered the rotational properties, ab initio calculations were performed on DME at the MP2 level with both the 6-31G* and the 6-31 l++G** basis sets. Wiberg and co-workers have reported that this level of correlation correction was necessary to correctly predict a gauche preference in 1,2-difluor~ethane.~~~~~ Although calculations at the MP2/6-31G* level lowered the energy difference between the gauche and anti rotamers in DME, the anti conformer was still preferred, by approximately 1.0 kcal mol-'. Calculations performed with the more flexible 6-31l++G** basis set reduced the gauche/anti energy difference by only 0.16 kcal mol-', to approximately 0.83 kcal mol-' in favor of the anti rotamer. The complete failure of the AM1 semiempirical molecular orbital calculations to predict the existence of a stable gauche rotamer of DME is evident from an inspection of Figure 6. The poor performance of AM1 may be attributed to its lack of polarization functions. On the basis of these limited studies we concluded that in the gas-phase DME would adopt an anti orientation. Although we had parametrized the preceding torsion angles based on data generated at the HF/6-31G* level, we chose to parameterize the 0-C-C-0 torsion angle on the basis of the MP2/6-31G* data. The molecular symmetry associated with DME resulted in torsional parameters that contained no asymmetric phase corrections (VI = 1.34, VZ = -1.15, V3 = 0.77, kcal mol-', y l = y3 = O.O", y2 = 180.0'). A comparison of the ab initio and GLYCAM-93 energies for rotation about the C5-C6 bond in 4 (see Figure 7) indicated that the symmetric parameters accurately reproduced this asymmetric rotational profile. These results and the close similarity between the rotational profiles of DHE and DME led us to conclude that the new torsional term should be applied for all 0-C-C-0 torsion angles when explicit solvation is to be employed.

Partial Atomic Charges One of the most complex issues in developing a carbohydrate force field pertains to the choice of algorithm used to calculate and assign partial atomic charges. The possibility for hydroxyl and hydroxymethyl group rotation leads to ambiguity in the choice of conformation best suited for charge deviation. Moreover, a typical monosaccharide contains four or five chiral atoms. Epimerization at any of these centers can be expected to lead to a variation in the atomic charges in the vicinity of that atom. As the Coulombic interaction between partial atomic

J. Phys. Chem., Vol. 99, No. 11, I995 3837 charges on proximal sugar residues may influence the conformations adopted by an oligosaccharide, a detailed examination of the effects of monosaccharide configuration on the values of the charges was undertaken. Previously we have reported net atomic charges for ketopyranoses derived from fitting a monopole approximation to the molecular electrostatic potential generated with a HF wave function (ESP charges).99 One of the conclusions from that research was that in asymmetric molecules, such as sugars, it may be necessary to sample the quantum mechanical elecrostatic potential at a large number of sites to ensure charge convergence. Similar conclusions were also reached in a study related to charge variation as a function of bond rotation in formamide.'" However, even with a high sampling density, ESP charges may depend on molecular conformation.'0'J02 In some cases the charge fluctuations between neighboring atoms may be minimized by fitting the charge to a subset of the atoms. For example, fluctuations in the charges on the carbon and hydrogen atoms in methyl and methylene groups may appear quite large, but when these groups are treated as a single charge site, either through summation of the individual charges99J01,102 or through charge fitting only to the central carbon atom,lo3the united atom (UA) ESP charges appear much less sensitive to either conformation or method of derivation. Furthermore, in a previous study of oxocarbenium ions, of relevance to glycoside hydrolysis, we found that charges derive by fitting directly to the non-hydrogen atoms, rather than by summation, led to interaction potentials between those ions and water that were virtually indistinguishable from the same interactions derived from all-atom charge models.lo3 The calculation of the ab initio UA charges was performed with a modified version of subroutine L602 in Gaussian 92. To examine the effects of basis set and ESP-fitting algorithm, the ab initio ESP-charges were computed for methyl a-Dmannopyranoside (S), a sugar which occurs commonly in N-linked oligosaccharides. The STO-3G and 6-31G* ESP charges were calculated with a variety of ESP sampling algorithms, including the radially distributed array of points from the CHELP method,lM a high-density grid (CHELPG),lm and a high-density array of random points (PDQP).99 The ESP charges for 8 are presented in Tables 1 and 2. Between each of the all-atom models, in which monopoles were placed at all of the nuclear positions, certain trends appeared constant. As expected, hydroxyl oxygen atoms carried partial negative charges while the hydroxyl protons were partially positive. The carbon atoms were positive and the aliphatic protons carried the smallest charges. However, the magnitudes and relative orderings of the partial charges appeared to depend strongly on the ESP sampling method. For example, the CHELP ESP charges indicated that the methyl carbon atom is the most positive atom (0.674 and 0.772 au at the STO-3G and 6-31G* levels, respectively), whereas both PDQP (0.136 and 0.200 au) and CHELPG (0.056 and 0.196 au) indicated much smaller charges on this atom. Moreover, in the case of the aliphatic protons, CHELP predicted that all but H1 (in the case of the STO-3G data only) were partially negative, while both PDQP and CHELPG predicted that only H4 was slightly negative. The UA charge model had the least effect on the charges of the hydroxyl atoms, although there was somewhat more variation among these charges as calculated by CHELP than found by either PDQP or CHELPG. As anticipated, the methy4 group now carried a similar charge, independent of the calculational method. More significant to a carbohydrate force field was the observation that the UA model resulted in ESP charges

3838 J. Phys. Chem., Vol. 99,No. 11, 1995

Woods et al.

TABLE 1: Ab Initio STO-3G ESP-Partial Atomic Charges for 8 PDQP CHELPG CHELP atom all united all united all united c- 1 0.184 0.494 0.117 0.473 0.450 0.555 c-2 0.160 0.145 0.181 0.142 0.313 0.150 c-3 0.225 0.240 0.145 0.224 0.360 0.278 c-4 0.426 0.359 0.517 0.359 0.354 0.357 c-5 0.083 0.178 0.013 0.143 0.140 0.194 C-6 0.145 0.214 0.160 0.227 0.345 0.220 C-Me 0.136 0.216 0.056 0.203 0.674 0.235 0-1 -0.379 -0.473 -0.324 -0.442 -0.597 -0.536 0-2 -0.545 -0.567 -0.547 -0.560 -0.612 -0.592 0-3 -0.552 -0.548 -0.546 -0.546 -0.582 -0.541 0-4 -0.599 -0.601 -0.594 -0.585 -0.610 -0.617 0-5 -0.344 -0.447 -0.306 -0.425 -0.475 -0.472 -0.5 12 -0.537 -0.506 -0.533 -0.544 -0.523 0-6 H- 1 0.100 0.115 0.021 H-2 0.025 0.026 -0.066 H-3 0.002 0.014 -0.042 H-4 -0.002 -0.023 -0.016 0.040 -0.001 H-5 0.037 H-6A 0.030 0.033 -0.027 H-6B 0.017 0.014 -0.038 H1-Me 0.014 0.033 -0.122 H2-Me 0.020 0.038 -0.123 H3-Me 0.023 0.041 -0.103 HO- 1 0.337 0.346 0.342 0.346 0.343 0.354 HO-2 0.313 0.312 0.321 0.316 0.3 11 0.295 HO-3 0.336 0.342 0.326 0.334 0.338 0.337 HO-4 0.319 0.327 0.314 0.324 0.308 0.305 dipole" 1.556 1.552 1.555 1.5 14 1.536 1.645 rmS* 0.982 1.057 0.988 1.122 1.889 2.072 lTlllSC 13.81 14.87 12.61 14.33 22.11 24.25 NPTSd 32645 32645 14889 14889 1617 1617 a SCF dipole moment = 1S70 D. Root-mean-square deviation between quantum mechanically calculated and point charge derived electrostatic potentials in kcal mol-'. Relative rms in percent. Total number of points used in the charge-fitting routine. for the ring carbon atoms that could be ranked in the same decreasing order, namely, C1 > C4 > C3 > C5 > C2, regardless of the computational method or the basis set.

Dipole Moments A quantitative measure of the ability of the ESP charge models to reproduce the molecular electrostatic potential may be obtained by comparing the dipole moments calculated from the all atom and UA charge models with that calculated quantum mechanically (see Tables 1 and 2). The best agreement between all atom and UA dipole moments is found for the PDQP ESP charges. The CHELPG ESP charges also performed reasonably well, underestimating the UA dipole by only approximately 0.05 D. The difference in dipole moments was largest for the CHELP ESP charges, being approximately 0.10 D, in this case overestimating the dipole moment in the UA model. Similar conclusions as to the goodness of fit to the electrostatic potentials may be reached by examining the root-mean-square (rms) deviation between the quantum mechanically calculated and point-charge-derived electrostatic potentials and the relative rms (rrms) values for the all atom and UA models. The rms values for all three methods indicated that a UA model gave rise to a marginally worse reproduction of the electrostatic potentials than did the all atom model. While both the PDQP and CHELPG ESP charges gave rise to similar values for the Rh4S and RRMS, the corresponding values based on CHELP ESP charges were larger, by a factor of approximately 2. As the UA model contained approximately one-third fewer monopoles than the all-atom case, the degree to which the dipole moments and rms values were reproduced, as compared to the all atom model, was noteworthy.

Conformational Charge Dependence Presented in Figure 8 are the net atomic charges for selected atoms during rotation of the hydroxymethyl group in 4 calculated using the CHELP algorithm and a UA model. It is readily apparent that there is no simple relationship between the 06-C6-C5-05 torsion angle and the charges. Furthermore, the magnitudes of the charges vary by approximately 0.1 au for each atom with the exception of the charge on the 06methyl group which shows less dependence on conformation. Thus, it was not possible to decompose the structure unambiguously into subsets of atoms whose charges were affected uniquely by the rotation. Moreover, since rotation of the hydroxymethyl group altered the charges on the atoms in the pyran ring, incorporation of average charge values would require that the rotational variations be determined for each individual monosaccharide. While this is feasible, the variations in charge due to the orientation of each hydroxyl group would also have to be determined, and this is clearly impractical. While other studies have employed average charge values for monosaccharides, the extreme variations in the charges of the ring atoms clearly illustrates the approximate nature of such approaches."'@ Consequently, we elected to determine net atomic charges for the solid-state conformation of each monosaccharide or, in the case of model structures, for the HF/6-31G* optimized geometries. We believe the CHELPG algorithm offers a sufficient degree of charge convergence for general application to carbohydrates.

Anomericity and Partial Atomic Charges Related to the question of conformational charge dependence is the effect of anomeric configuration. For free sugars in aqueous solution, anomerization is an equilibrium process

MD Parameters for Glycoproteins and Oligosaccharides

J. Phys. Chem., Vol. 99, No. 11, 1995 3839

TABLE 2: Ab Initio 6-31G* Partial Atomic Charges for 8 PDQP CHELPG CHELP atom all united all united all united c-1 0.375 0.694 0.333 0.686 0.717 0.738 c-2 0.141 0.127 0.175 0.119 0.379 0.167 c-3 0.203 0.308 0.093 0.291 0.374 0.329 c-4 0.658 0.465 0.763 0.494 0.448 0.437 c-5 0.066 0.254 0.042 0.224 0.186 0.291 C-6 0.234 0.265 0.241 0.282 0.5 11 0.263 C-Me 0.200 0.295 0.196 0.289 0.772 0.303 0- 1 -0.527 -0.614 -0.510 -0.593 -0.772 -0.675 0-2 -0.703 -0.715 -0.722 -0.715 -0.774 -0.735 0-3 -0.740 -0.745 -0.750 -0.768 -0.729 -0.690 0-4 -0.806 -0.784 -0.824 -0.788 -0.788 -0.792 0-5 -0.490 -0.622 -0.482 -0.621 -0.648 -0.633 0-6 -0.693 -0.708 -0.708 -0.722 -0.713 -0.664 H- 1 0.112 0.126 -0.007 H-2 0.036 0.033 -0.094 H-3 0.027 0.054 -0.033 H-4 -0.028 -0.048 -0.024 H-5 0.055 0.046 0.006 H-6A 0.019 0.021 -0.064 H-6B 0.016 0.018 -0.068 H1-Me 0.002 0.003 -0.139 H2-Me 0.041 0.043 -0.112 H3-Me 0.027 0.027 -0.118 HO- 1 0.448 0.448 0.464 0.455 0.436 0.303 HO-2 0.439 0.435 0.460 0.455 0.407 0.438 HO-3 0.455 0.461 0.460 0.465 0.450 0.387 0.444 0.448 0.399 0.446 HO-4 0.434 0.439 2.242 2.226 dipole" 2.261 2.206 2.179 2.261 llllSb 1.220 1.341 1.132 1.338 2.361 2.584 lm-lSC 12.91 14.19 10.66 12.60 21.50 23.53 32645 32645 NFTSd 14889 14889 1617 1617 " SCF dipole moment = 2.236 D. Root-mean-square deviation between quantum mechanically calculated and point charge derived electrostatic potentials in kcal mol-'. Relative rms in percent. Total number of points used in the charge-fitting routine. TABLE 3: Effect of Anomericity on the HF/6-31G* ESP-Partial Atomic Charges for 1 and 2 1 2 ""Me

0.1

Net

7 0.8 '"'Me

0.2

atom

all atom" united atoma all atomc united atomd 0.597 0.326 0.489 0.456 0.006 c-2 -0.204 0.086 -0.257 c-3 0.126 -0.045 0.216 0.023 c-4 -0.064 0.043 -0.176 -0.022 c-5 0.217 0.185 0.342 0.281 C-Me 0.146 0.251 0.135 0.251 0- 1 -0.503 -0.458 -0.444 -0.482 0-5 -0.480 -0.387 -0.546 -0.513 H- 1 -0.018 0.032 H-2a 0.076 0.067 H-2e 0.034 0.054 H-3a -0.030 -0.028 H-3e -0.014 -0.028 H-4a 0.027 0.054 H-4e -0.001 0.011 H-Sa -0.020 -0.003 H-5e 0.031 0.001 H1-Me 0.025 0.011 H2-Me -0.003 0.026 H3-Me 0.057 0.043 dipole' 1.8169 1.9545 0.3682 0.3905 " R m s = 1.100 kcal mol-', NPTS = 12368. b R m s = 1.418 kcal mol-'. 'Rms = 1.191 kcal mol-l, NFTS = 12 131. dRms = 1.363 kcal mol-'. SCF dipole moments for 1 and 2, 1.8314 and 0.3937 D, respectively.

c- 1

0

Atomic -o,l Charge (ad.) -0.2

+Q(Me6) -e-

Q(06) +Q(C6)

0

30

t-Q(C5)

+Q(05)

---t Q(C4)

60 90 120 150 180 210 240 270 300 330 360

06-C6-C5-05 Torsion Angle (deg) Figure 8. United atom HF/6-31G* ESP-charges (Q) in 4 as a function of methoxymethyl group rotation. Numbering as in hexopyranoses.

leading to the interconversion of the a- and p-forms of the sugar. Any differences in atomic charge that are associated with anomeric configuration may be important for accurate aqueous simulations of sugars. We wished to determine the extent to which a- and B-anomers displayed similar charge distributions and chose to examine the net atomic charges calculated for 1 and 2 at their HF/6-31G* optimized geometries. While the same general trends in charge values were present in both anomers, the data in Table 3 indicated that epimerization at C1 led to a variation in all of the atomic charges and in the molecular dipole moments. The charge variations were largest for the ring atoms, both in the all atom and in the UA models. A direct comparison

of the differences between the charges on 1 and 2 was provided by the UA charges. That the variation between the all atom

3840 J. Phys. Chem., Vol. 99,No. 11, 1995

Woods et al.

TABLE 4: Effect of Hydroxyl Configuration on the HF/6-31G* ESP Partial Atomic Charges for 5 and 6 5

6

4

%

Energy (kcalhol)

Q,,,

a,,

atom c-1 c-2 c-3 c-4 c-5 C-Me 0- 1 0-2 0-5 HO-2 H-1 H-2a H-2e H-3a H-3e H-4a H-4e H-5a H-5e H1-Me H2-Me H3-Me dipole'

Relative

"Me

"'Me

all atomu 0.511 0.240 0.032 -0.070 0.256 0.156 -0.509 -0.732 -0.509 0.452 0.051

united atomb 0.564 0.200 -0.007 0.002 0.271 0.253 -0.525 -0.695 -0.506 0.443

-0.018 -0.001 -0.005 0.044

-0.002 0.013 0.017 0.004

0.021 0.049 1.8286

1.8232

all atomc 0.167 0.325 -0.021 -0.128 0.231 0.082 -0.361 -0.701 -0.467 0.408 0.100 0.055 0.016 0.029 0.034 0.038 0.022 0.033 0.039 0.052 0.046 2.2613

united atomd 0.420 0.406

-0.007 -0.038 0.294 0.273 -0.499 -0.718 -0.535 0.404

2.2775

R m s = 1.212 kcal mol-', NPTS = 12 430. R m s = 1.414 kcal mol-'. 'Rms = 1.166 kcal mol-', NPTS = 12882. d m= 1.374 kcal mol-1. SCF dipole moments for 5 and 6, 1.8646 and 2.2424 D,

respectively. charges in 1 and 2 was not an artifact of ESP sampling was indicated by the observation that the methyl groups in each anomer had identical UA charges. Notably, the UA charges predicted that epimerization at C1 had the least effect on the charges of the anomeric oxygen atom and aglyconic methyl group. The largest UA charge variations were predicted for ring atoms C1, C5, and 05.

Hydroxyl Group Codlguration Related to anomericity is a variation in hydroxyl group configuration at nonanomeric carbon atoms. In naturally occurring oligosaccharides both axial (manno) and equatorial (gluco) hydroxyl groups are present at C2. As these groups may influence the properties of the glycosidic linkages, we undertook an examination of the effects of configuration at C2 on charges and rotational properties. We selected two simple models, with a hydroxyl group at C2 in either an axial (5) or equatorial (6) configuration, that were closely related to 2. Presented in Table 4 are the partial atomic charges and dipole moments for the HF/6-3 lG* optimized geometries of 5 and 6. The calculations indicated that the configuration of the hydroxyl group at C2 had a relatively small effect on the molecular dipole moments. Analogous to the case of 1 and 2, the methyl groups and the acyclic oxygen atoms, including 0 2 , were largely unaffected by epimerization at C2. However, the charges at both C1 and C2 differed significantly between the two epimers. Moreover, the charges at C1 in 5 and 6 were found to differ from the corresponding charge in 2. To examine further the effect of configuration at C2 on the properties of the glycosidic linkage, we computed the HF/631G* energies for rotation of the 0-C-0-C torsion angles in 5 and 6. The energy profile for rotation of the axial methoxyl

:h E 4

+HFIE-310' +GLYCAM-93 ) "gluco'

2

0

30 60 SU 120 150 150 210 240 270 300 330 360

0541-01-C Torsion Angle (deg)

Figure 9. Methoxyl group rotation in "manno-" 5 (0)and "gluco-" 6 (A) analogues of 2. The total energy for the minimum energy conformations of 5 and 6 at the HF/6-31G* level were -458.759 232 and -458.763 931 au, respectively. group in 5 (a model for mannopyranosides) appeared to be essentially equivalent to that found for 2 with only a modest increase in the rotational barrier at 240' of approximately 1.5 kcal mol-' (see Figure 9). However, in the case of 6, the maximum barrier to rotation (at TO-C-0-c= 240") was larger than that in 2, by approximately 3.5 kcal mol-'. For each rotation, the ab initio and GLYCAM-93 data were in good agreement. The larger barrier in 6 than in 5 reflects the presence of steric repulsions between the equatorial-hydroxyl and the axial-methoxyl groups in 6. Furthermore, the plateau in the energies found for 2 and 5, corresponding to 0-C-0-C torsion angles of approximately 150', was absent in 6. This observation should manifest itself as a higher preference in polar solvents for the exo-anomeric orientation in glucopyranosides than in mannopyranosides. We drew two conclusions from these charge studies. First, it was apparent that each monosaccharide would require a separate set of partial charges, calculated using a high sampling density (CHELPG) to account correctly for the dependence of the charges on molecular configuration and conformation. Both the UA and all-atom charge models performed well in our analyses. However, to be able to adopt the standard AMBER practice of scaling all 1-4 electrostatic interactions by a factor of 0.5, the more consistent method is to fit charges to all of the atoms. The close agreement between the GLYCAM-93 and the ab initio rotational energy profiles for both 5 and 6 indicated that the parameters (see Table 5 ) derived from our choices of representative fragments and ESP charges were capable of describing the subtle effects on glycosidic linkage conformation between different monosaccharides.

Hydrogen Bonding To examine the suitability of the HF/6-31G* ESP-charges for modeling hydrogen bonding, we calculated the interaction energy between water and methanol. Such an interaction provided a simple model for the aqueous solvation of a sugar hydroxyl group. Interaction energies were determined for complexes in which water was either the proton donor (type I interaction) or the proton acceptor (type II). As the ESP charges computed for the water monomer (q(H) = 0.410) were essentially equivalent to those of the TIP3P water model (q(H) = 0.417),'05 and as we intended to employ the TIP3P water in molecular dynamics simulations with GLYCAM-93, we chose the TIP3P monomer geometry and charges for these calculations. In the case of methanol monomer, the experimental gas-phase equilibrium geometry was employed with the HF/6-31G* ESP charges.'" The interaction between water and methanol has been the subject of several molecular orbital calculations, both with fmed monomer g e ~ m e t r i e s ' ~ ~ -and ~ ' l more recently with complete

MD Parameters for Glycoproteins and Oligosaccharides

J. Phys. Chem., Vol. 99, No. 11, 1995 3841

TABLE 5: GLYCAM-93 Parameter# for Oligosaccharides and Glycoproteins atom types

KS

r,

KS

atom types

rM

Stretching AC-CT AC-HC AC-OS AC-OG AC-OH EC-CT AC-CT-CT AC-CT-HC AC-CT-OH AC -CT- OG AC-CT-OS AC-OG-CT AC-OH-HO AC-OS-CT EC-CT-CT EC-CT-HC EC-CT-N EC-CT-OH EC-CT-OG EC-CT-OS EC -OG- CT CT-EC-OS CT-CT-OG HC-AC-OG HC-AC-OH HC-AC-OS HC-EC-HC HC-EC-N HC-EC-OG

310.0 331.0 320.0 320.0 320.0 310.0 40.0 35.0 50.0 50.0 50.0 60.0 55.0 60.0 40.0 35.0 80.0 50.0 50.0 50.0 60.0 50.0 50.0 35.0 35.0 35.0 35.0 35.0 35.0

1.527 1.090 1.416 1.405 1.396 1.519 Bending 111.80 109.50 109.50 109.50 109.50 113.80 108.50 113.50 109.40 109.50 109.70 109.50 109.50 109.50 114.70 109.50 109.50 109.50 109.50 109.50 109.50 109.50 109.50

EC-HC EC-N EC-OS EC-OG EC-OH CT-OG

331.0 337.0 320.0 320.0 320.0 320.0

EC-OH-HO EC-OS -CT EC-N-C EC-N-H CT-AC-HC CT-AC-OG CT-AC-OH CT- AC -OS CT-EC-HC CT-EC-N CT-EC-OG CT-EC-OH HC -EC -OH HC -EC -OS HC-CT-OG N-EC-OS OG-AC-OS OH-AC-OS OG-EC-OS OH-EC-OS

55.0 60.0 50.0 38.0 35.0 50.0 50.0 50.0 35.0 80.0 50.0 50.0 35.0 35.0 35.0 106.9 110.7 110.7 110.7 110.7

1.090 1.460 1.425 1.389 1.387 1.435 108.50 111.90 121.90 118.40 109.50 107.80 109.50 110.50 109.50 109.70 108.70 109.50 109.50 109.50 109.50 107.90 112.00 111.60 107.60 107.20

Torsion atom types X-CT-AC-X X-CT-EC-X X-CT-OG-X X-AC-OS-X X-EC-OS-X CT-C-N-EC CT-EC-N-C AC-CT-N-C EC-CT-N-C EC-N-C-0 H-N-EC-HC H-N-EC-CT H-N-EC-OS H-N-CT-AC H-N-CT-EC HC-EC-N-C OS-EC-N-C 6-31G* ESP-Charges OS-EC-N-C

STO-3G ESP-Charges CT-AC-OG-CT CT-EC-OG-CT CT-OG-AC-OS CT-OG-EC-OS HC-AC-OG-CT HC-EC-OG-CT OH-CT-CT-OG CT-CT-OG-EC AC-CT-OG-AC

V"l2 1.30 1.30 1.15 1.15 1.15 10.0 0.00 0.62 0.43 -0.18 0.62 0.43 -0.18 10.00 0.00 0.00 0.00 0.00 0.00 0.00 2.39 1.40 -0.41 0.64 1.91 -0.87 0.00 0.00 1.39 0.70 0.91 0.85 0.74 0.97 0.00 0.00 1.34 0.77 0.20 0.383 0.20 0.383

Y ,* 0.00,3 0.00,3 0.00,3 0.00,3 0.00.3 180.00, 2 0.00, 1 20.76, 1 26.06,2 339.31, 3 20.76, 1 26.06, 2 339.31, 3 180.00,2 0.00,l 0.00, 1 0.00.1 0.00, 1 0.00, 1 0.00,l 200.31, 1 359.33.2 175.02, 3 341.97, 1 3.37, 2 173.13, 3 0.00,3 0.00,3 275.77, 1 312.04,2 347.72, 3 144.06, 1 353.55, 2 6.43, 3 0.00,3 0.00,3 0.00,l 0.00,3 180.00,2 0.00,3 180.00,2 0.00,3

atom types AC-CT-OG-EC EC-CT-OG-

AC

EC-CT-OG-EC OH-CT-CT-OS OH-CT-AC-OS OS-CT-CT-OS OH-CT-AC-OG OH-AC-CT-OH OH-EC-CT-OH OH-CT-EC-OG OS-CT-AC-OG OS-CT-EC-OG OH-CT-EC-OS AC-CT-OH-HO EC-CT-OH-HO CT-AC-OH-HO HC-AC-OH-HO HC-EC-OH-HO OS-AC-OH-HO OS-EC-OH-HO CT-CT-OG-AC

VJ2 0.20 0.383 0.20 0.383 0.20 0.383 1.34 -1.15 0.77 1.34 -1.15 0.77 1.34 -1.15 0.77 1.34 -1.15 0.77 1.34 -1.15 0.77 1.34 -1.15 0.77 1.34 -1.15 0.77 1.34 -1.15 0.77 1.34 -1.15 0.77 1.34 -1.15 0.77 0.50 0.50 0.50 0.50 0.50

0.50 0.50 0.20 0.383

Y,n 180.00,2 0.00,3 180.00,2 0.00,3 180.00,2 0.00,3 0.00,l 180.00,2 0.00,3 0.00,l 180.00,2 0.00,3 0.00,l 180.00,2 0.00,3 0.00,l 180.00,2 0.00,3 0.00,l 180.00,2 0.00,3 0.00,l 180.00,2 0.00,3 0.00,l 180.00,2 0.00,3 0.00,l 180.00,2 0.00,3 0.00,l 180.00.2 0.00.3 0.00,l

180.00,2 0.00,3 0.00,3 0.00,3 0.00,3 0.00,3 0.00,3 0.00,3 0.00,3 180.00,2 0.00,3

3842 J. Phys. Chem., Vol. 99, No. 11, I995

Woods et al.

TABLE 5 (Continued) Atomic Mass atom types

mass

AC EC

12.01 12.01 16.00

OG

van der Waals

atom types

R*

E

AC

1.80 1.80 1.65

0.06 0.06

EC OG

0.15

Hydrogen Bond atom pairs H, OG

C

D

7557

H W ,OG

7557

2385 2385 2385

HO, OG 7557 Units are as follows: bond lengths (A), angles (deg), K,(kcal mol-' A2), Kb(kcal mol-' rad2), V (kcal mol-'), y (deg), mass (au), C (kcal mol-' e (kcal mol-'). A'Z),D (kcal mol-' geometry optirnization.'12 No experimental data have been reported with regard to either the geometries or interaction energies for the methanol-water complexes. Gas-phase interaction energies for the water dimer (4.5-5.0 kcal mol-') and for the methanol dimer (3.2-7.3 kcal mol-') provide only a rough estimate for the water-methanol interaction energy.' l3 At the ab initio HF/6-3 1G* level the interaction energies appear to be comparatively insensitive to the relative orientation of the donor and acceptor molecules. In the type I complex several values for the HF/6-31G* interaction energy have been reported, namely, -5.55,'12 -5.50,'09 and -5.73110 kcal mol-' at 0.0 separations of 2.95, 2.92, and 2.97 A, respectively. Corresponding values for the type II complex were reported as -5.59112 and -5.55'1° kcal mol-' at 2.96 and 2.97 8, respectively. None of these values included corrections for basis set superposition errors (BSSE). (For a discussion of the importance of these corrections in hydrogen bond calculations, see ref 114.) In the type I1 complex correction for the BSSE was reported to lead to a reduction in the interaction energy of approximately 1.5 kcal We selected the geometries of the complexes from ref 110. The approach trajectories for complex formation were determined by varying the 0.0 separations in 0.2 A increments while maintaining all other geometrical parameters constant (see Figure 10). For the type I complex the GLYCAM-93 results indicated a minimum in the interaction energy corresponding to a hydrogen bond strength of 6.69 kcal mol-' at an 0 0 separation of approximately 2.76 A. In the type II complex the minimum was also predicted to occur at an 0 0 separation of approximately 2.76 8, and had an energy of 6.98 kcal mol-l. These values compare favorably with the TIP3P water dimerization energy (6.50 kcal mol-') and 0.0separation (2.74 A).io5 Furthermore, the relative strengths of the type I and 11interactions were correctly predicted. However, the GLYCAM-93 value for the 0 0 separation in these com lexes is shorter than the ab initio values, by approximately 0.2 and the computed interaction energies are intermediate between those determined at the HF/6-3 lG* and those at the MP2/6-31G* levels (MP2 -7.56 kcal mol-' type I, -7.80 kcal mol-' type II).l12 Further insight into the suitability of the HF/6-31G* ESP charges for modeling hydrogen bonding in solution was provided through computation of the radial distribution functions for oxygen-oxygen (goo)and oxygen-hydrogen (gH0) interactions from an MD simulation of methanol in water. Radial distribution functions provide estimates of the degree of structural ordering of a solvent relative to a randomly oriented

1

Type I: R=H, R'=CH3 Type 11: R=CH3, R=H 01

i

-l -2

Interaction Energy (kcal/rnol)

-3 -4-5

U T +Type

W 1: HF/6-310'ESP II:HF/6-310* ESP

-6 ~

-7-8 J 2

2.5

3

3.5

4

4.5

5

5.5

6

O...O Separation (Angstroms) Figure 10. Interaction energies for water and methanol calculated with HF/6-31G* ESP-charges for type I (0)and type II (A) complexes.

arrangement of solvent molecules. The comparatively short (35 ps) NPT simulation was performed at 300 K and approximately 1 atm. The methanol molecule was solvated by a total of 296 previously equilibrated water molecules. A nonbonded cutoff value of 10 8, was employed. Bond lengths were constrained using the SHAKE algorithm, which is compatible with the use of the TIP3P water model.li5 The first 10 ps of the simulation was allocated for equilibration. The computational formalism has been previously described in f~ll.~'9*~ In close agreement with the results for pure TIP3P water,105J16 0 the maximum in the methanol-water goo occurred at an 0. separation of 2.80 8,(see Figure 1 la). The characteristic second peak in the goo at 4.5 8, was separated from the first by a minimum at approximately 3.5 A. Integration of the f i s t peak to 3.5 A provided an estimate of the number of nearest neighbors (3.73). This value is smaller than that for pure water (5.0), which is consistent with the loss of a hydrogen bond donor atom in methanol relative to water. A estimate of 4.2 nearest neighbors has been computed recently for a hydroxyl group in the disaccharide maltose in TIP3P water.'17 The position and intensity of the second peak varied only slightly from that computed for pure TIP3P ~ a t e r , ' ~ indicative ~J'~ of a comparatively small solute influence.

MD Parameters for Glycoproteins and Oligosaccharides 41

O

i

1

2

3

l

4

5

8

r (Angstroms)

r (Angstroms)

r (Angstroms)

a

b

a=-%

than the value for the number of nearest water molecules (3.73), a result that is also found for TIP3P water and may be expected to arise from the dynamic nature of the solvent.105 That water and methanol should display similar hydrogen bond patterns in aqueous solvent, within the limits of their structural similarities, would be predicted on the basis of their computed gas-phase interaction energies. That this expectation has been verified through the computed radial distribution functions supports considerably the choice of partial atomic charges, as well as the suitability of l V 3 P water in dynamics simulations employing our parameters.

Oligosaccharide Generation

C

Figure 11. Radial distribution functions for methanol in water: (a) goo, (b) gH0 methanol acting as a proton acceptor, (c) gHO methanol acting as a proton donor.

o.0w

J. Phys. Chem., Vol. 99, No. 11, 1995 3843

o=om

Figure 12. Linkage formation methodology for typical 0-and N-linked glycoproteins showing net charges (Q) on the component fragments. Methanol may act as either a proton donor or acceptor and consequently there are two possible gHO curves. For TIP3P water, the gHO displays two peaks of approximately equal intensity at approximately 1.9 and 3.3 A very similar function was computed for the interaction between TIP3P water and the methanolic oxygen atom, in which the maxima in gHO occurred at 1.85 and 2.80 8, (see Figure llb). Integration of the f i s t peak to 2.5 A indicated that 1.7 hydrogen bonds were formed with this oxygen atom. When methane acted as a proton donor, the gH0 was calculated to contain a single sharp peak at 1.90 A, which when integrated to 2.5 A indicated a single hydrogen bond was formed (see Figure 1IC). The total of these integrations (2.7 hydrogen bonds/methanolmolecule) is smaller

Within the standard AMBER program each amino acid or nucleoside residue is stored in a structural database. The three dimensional geometry for each residue is defined by a z-matrix, and associated with each atom is its net atomic charge and atom type. We have created an analogous database for monosaccharides. However, unlike a protein backbone, in which the linkage sites between amino acids are unambiguous, the possibility of several linkage positions between monosaccharides in a carbohydrate polymer necessitated a modified strategy for the linking of the individual monosaccharides. One alternative that has been employed previously by us involves the creation of a single structural file for the complete oligosaccharide.''* That approach avoided the problem of defining linkage sites and allowed complete flexibility in the assignment of net atomic charges and atom types. However, the generation of a z-matrix for a complete oligosaccharide is extremely time consuming. Moreover, such an approach is not readily amenable to automation of the sequence generation. In this work we chose to apply a method related to that used by AMBER to link amino acid residues together. With this method we created individual structural files for each monosaccharide based on the solidstate conformations. Ab initio HF/6-3 lG* ESP-charges were derived as described above for each residue and assigned to all of the atoms. In the chemical formation of a glycosidic linkage (see Figure 12) a condensation reaction between two monosaccharides occurs in which one molecule of water is expelled. Of relevance to this work was the choice of model for each chemical component of this condensation, that is, a choice had to be made

TABLE 6: HF/6-31G* ESP Charges on Oligosaccharide Components Related to 8 charges (au) atom MA 2MA 3MA 4MA 6MA c1 c2 c3 c4 c5 C6 H1 H2 H3 H4 H5 H6A H6B 02 03 04 06 H02 H03 H04 H06

0.333 0.175 0.093 0.762 0.042 0.241 0.126 0.033 0.054 -0.048 0.046 0.022 0.0 19 -0.722 -0.750 -0.824 -0.482 -0.708 0.464 0.460 0.460 0.444

0.333 0.175 0.093 0.762 0.042 0.241 0.126 0.033 0.054 -0.048 0.046 0.022 0.019 -0.498 -0.750 -0.824 -0.482 -0.708

total

0.240

05

0.460 0.460

0.333 0.175 0.093 0.762 0.042 0.241 0.126 0.033 0.054 -0.048 0.046 0.022 0.019 -0.722 -0.530 -0.824 -0.482 -0.708 0.464

0.333 0.175 0.093 0.762 0.042 0.241 0.126 0.033 0.054 -0.048 0.046 0.022 0.019 -0.722 -0.750 -0.604 -0.482 -0.708 0.464 0.460

0.444

0.460 0.444

0.444

O.Oo0

O.Oo0

O.Oo0

0.333 0.175 0.093 0.762 0.042 0.241 0.126 0.033 0.054 -0.048 0.046 0.022 0.019 -0.722 -0.750 -0.824 -0.482 -0.504 0.464 0.460 0.460

O.Oo0

26MA

36MA

0.333 0.175 0.093 0.762 0.042 0.241 0.126 0.033 0.054 -0.048 0.046 0.022 0.019 -0.498 -0.750 -0.824 -0.482 -0.504

0.333 0.175 0.093 0.762 0.042 0.241 0.126 0.033 0.054 -0.048 0.046 0.022 0.019 -0.722 -0.530 -0.824 -0.482 -0.504 0.464

0.460 0.460

0.460

-0.240

-0.240

3844 J. Phys. Chem., Vol. 99,No. 11, 1995 as to which component would be assumed to have contributed a proton and which a hydroxyl group. Our choice of linkage structure was based on the observation that, when combined with the LINK module in AMBER, it led to oligosaccharides in which the glycosidic linkages adopted the exo-anomeric orientation.

F'rocedure for Assigning Charges In calculations involving proteins and oligosaccharides, compatibility with the current AMBER structural database necessitated employing HF/STO-3G ESP charges, despite the better performance of 6-31G* ESP charges in terms of reproducing experimental properties, such as dipole and quadrupole moments and hydrogen bond energies.llg Our torsional parameters were derived for HF/6-31G* ESP charges. However, with the exception of the 0-C-N-C torsion term for N-linked glycoproteins, we found that the rotational profiles were essentially unaltered when STO-3G ESP charges were employed. Thus, we have included in Table 5 a set of torsion parameters for 0-C-N-C linkages for use with STO-3G ESP-charges in N-linked glycoproteins. The next release of AMBER will employ HF/6-31G* ESP charges and thus will be completely compatible with our parameters.120 For each monosaccharide a structural file containing partial atomic charges was assembled according to the following protocol. In the case of terminal residues, the ESP charges for the solid-state conformation of the methyl glucoside were calculated. The total molecular charge excluding the charge on the 0-Me group was then determined. In the case of intermediate residues, the charge on the linking hydroxyl proton was subtracted from Qrem to give the charge on the intermediate residue The difference between Qmtand 0.0 au was determined and added to the charge on the linking oxygen atom (see Table 6). This approach facilitates oligosaccharide generation by ensuring that each nonterminal sugar residue has zero net charge and is analogous to that employed in AMBER for amino acid residues. The difference in charge between the deprotonated aglycon and Qremwas then added to the charge on the heteroatom of the aglycon. This method of charge partitioning accounts directly for the observation that the glycosidic oxygen atom carries less charge than a free hydroxyl group and in general is based on the assumption that the largest changes in partial charge during glycoside formation would be localized in the vicinity of the glycosidic linkage. Furthermore, it leads to the correct molecular charge while enabling the incorporation of charge variations arising from changes in hydroxyl group configurations.

(ere&

(emt).

Conclusions The AMBER force field has been extended to include parameters that reproduce the structural features and conformational preferences, as predicted by ab initio calculations, of a series of tetrahydropyran derivatives related to glycopyranosides. A method for the application of these parameters in the molecular modeling of oligosaccharides and glycoproteins has been presented. The effects of conformation and configuration on partial charges have been discussed. The use of partial atomic ESP charges derived at the ab initio HF/6-31G* level has been proposed and shown to be consistent with the TIP3P model for water. We have presented a protocol for deriving partial atomic charges for monosaccharides that facilitates the incorporation of any monosaccharide into our parameter set. High level ab initio calculations have been performed that indicate that the gauche effect in 1,2-dimethoxyethanedoes not exist in the gas phase and our parameters have been derived to fit that data. Our parametrization correctly predicts the subtle

Woods et al. effects on the rotational properties of the glycosidic linkage that result from hydroxyl group epimerization at C2.

Acknowledgment. We thank the North Carolina Supercomputing Center and the Oxford Centre for Molecular Sciences for use of their facilities. This work was supported in part by a grant from Monstanto Co. and by the National Sciences and Engineering Research Council of Canada, in the form of a Fellowship (to R.J.W.). Supplementary Material Available: An AMBER-style database containing geometries and HF/6-3lG* ESP-partial atomic charges for sugars commonly found in glycoproteins, including a-and /3-D-manno-, galacto-, and glucopyranose and j3-~-fuco-and 2-deoxy-2-N-acetyl-/3-~-glucopyranose (34 pages). Ordering information is given on any current masthead page. A modified version of Gaussian 92 subroutine L602 will be supplied upon request to RWOODS @MOND1.CCRC.UGA. EDU. References and Notes (1) Rademacher, T. W.; Parekh, R. B.; Dwek, R. A. Annu. Rev. Biochem. 1988,57, 785-838. (2) Varki, A. Glycobiology 1993,3,97-130. (3) Gerken, T.A.; Butenhof, K. J.; Shogren, R. Biochemistry 1989, 28,5536-5543. (4) Kalyan, N. K.; Bahl, 0. P. J. Biol. Chem. 1983,258,67-74. (5) Wittwer, A. J.; Howard, S. C.; Cam, L. S.;Harakas, N. K.; Feder, J.; Parekh, R. B.; Rudd, P. M.; Dwek, R. A.; Rademacher, T. W. Biochemistry 1989,28,7662-7669. (6) Rudd, P. M.; Joao, H. C.; Coghill, E.; Fiten, P.; Saunders, M. R.; Opdenakker, G.; Dwek, R. A. Biochemistry 1994,33,17-22. (7) Weis, W. I.; Drickamer, K.; Hendrickson, W. A. Nafure 1992,360,

127- 134. (8)Lemieux, R. U.;Bock, K. Arch. Biochem. Biophys. 1983,221,125134. (9) Homans, S.W. Prog. NMR Spectrosc. 1990,22,55-81. (10) Tvaroska, I.; Kozfir, T.; Hricovlni, M. In Computer Modeling of Carbohydrate Molecules; French, A. D., Brady, J. W., Eds.; Symposium Series 430;American Chemical Society: Washington, DC, 1990 pp 162176. (11) Carver, J. P.; Mandel, D.; Michnick, S. W.; Imberty, A.; Brady, J. W. In Computer Modeling of Carbohydrate Molecules; French, A. D., Brady, J. W., Eds.; Symposium Series 430;American Chemical Society: Washington, DC, 1990; pp 266-280. (12) Rice, K. G.; Wu, P.; Brand, L.; Lee, Y. C. Biochemistry 1991,30, 6646-6655. (13)Wooten, E. W.; Edge, C. J.; Bazzo, R.; Dwek, R. A.; Rademacher, T. W. Carbohydr. Res. 1990,203, 13-17. (14) Cano, F. H.; Foces-Foces, C.; Jimknez-Barbero, J.; Alemany, A.; BerbaM, M.; Martin-Lomas, M. J . Org. Chem.1987,52,3367-3372. (15)(a) Allinger, N. L. J . Am. Chem. SOC.1977,99,8127-8134. See, for example: (b) Woods, R. J.; Andrews, C. W.; Bowen, J. P. J. Am. Chem. SOC. 1992,114, 850-858. (c) Woods, R. J.; Andrews, C. W.; Bowen, J. P. J . Am. Chem. SOC.1992,114,859-864. (d) Imberty,A.; Tran, V.; Pkrez, S.J . Comput. Chem. 1989,11, 205-216. (e) Jeffrey, G. A.; Taylor, R. J . Comput. Chem. 1980,1,99-109. (f)Norskov-Lauritsen, L.; Allinger, N. L. J . Comput. Chem. 1984, 5, 326-335. (g) Tvaroska, I.; Pkrez, S. Carbohyd. Res. 1986,149,389-410.(h) French, A. D.; Tran, V. H.; Pkrez, S.In Computer Modeling of Carbohydrate Molecules; French, A. D., Brady, J. W., Eds.; Symposium Series 430 American Chemical Society: Washington, DC, 1990;pp 191-212. (i) Homans, S.W.; Dwek, R. A.; Boyd, J.; Mahmoudian, M.; Richards, W. G.; Rademacher, T. W. Biochemistry

1986,25,6342-6350. (16)Allinger, N. L.; Yuh,Y. H.; Lii, J.-H. J . Am. Chem. SOC. 1989, 111, 8551-8566. (17) Lii, J.-H.; Allinger, N. L. J. Am. Chem. Soc. 1989,111, 85668576. (18) (a) Lii, J.-H.; Allinger, N. L. J . Am. Chem. SOC.1989,111,85768582. See, for example: (b) French, A. D. J . Mol. Srrwcr. (THEOCHEM) 1993, 286, 183-201. (c) Dowd, M. K.; Reilly, P. J.; French, A. D. Biopolymers 1994,34,625-638. (19) (a) Brooks, B. R.; Bmccoleri, R. E.; Olafson,B. D.; States, D. J.; Swaminathan, S.;Karplus, M. J . Comput. Chem. 1983,4,187-217. See, for example: (b) Carver, J. P.; Mandel, D.; Michnick, S. W.; Imberty,A.; Brady, J. W. In Computer Modeling of Carbohydrate Molecules; French, A. D., Brady, J. W., Eds.; Symposium Series 430; American Chemical Society: Washington, DC, 1990;pp 266-280. (c) Ha, S.N.; Giammona, A.; Field, M.; Brady, J. W. Carbohydr. Res. 1988,180, 207-221. (d)

MD Parameters for Glycoproteins and Oligosaccharides Rasmussen, K. Acta Chem. Scand. 1982, A36, 323. (e) Yan, Z.-Y.; Bush, A. C. Biopolymers 1990, 29, 799-811. (20) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.;Rofeta, S., Jr.; Weiner, P. J. Am. Chem. SOC.1984,106,765784. (21) (a) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J. Cornput. Chem. 1986, 7, 230-252. See, for example: (b) Edge, C. J.; Singh,U. C.; Bazzo, R.; Taylor, G. L.; Dwek, R. A.; Rademacher, T. W. Biochemistry 1990,29,1971- 1974. (c) Homans, S. W. Biochemistry 1990, 29,9110-9118. (d) Scarsdale, J. N.; Ram. P.; Prestegard, J. H.; Yu, R. K. In Computer Modeling of Carbohydrate Molecules; French, A. D., Brady, J. W., Eds.; Symposium Series 430 American Chemical Society: Washington, DC, 1990; Chapter 15. (22) (a) Hermans, J.; Berendsen, H. J. C.; Van Gunsteren, W. F.; Postma, J. P. M. Biopolymers 1984,23, 1513-1518. See, for example: (b) Koehler, J. E. H.; Saenger, W.; van Gunsteren, W. F. Eur. Biophys. J . 1987, 15, 197-210. (c) Koehler, J. E. H.; Saenger, W.; van Gunsteren, W. F. Eur. Biophys. J. 1987, 15, 211-214. (d) Van Eijck, B. P.; Kroon, J. J . Mol. Struct. 1989, 195, 133-146. (e) Homans, S. W.; Pastore, A.; Dwek, R. A.; Rademacher, T. W. Biochemistry 1987, 26, 6649-6655. (23) (a) Thogersen, H.; Lemieux, R. U.; Bock, K.; Meyer, B. Can. J . Chem. 1982,60,44-57. See, for example: (b) Cumming, D. A.; Carver, J. P. Biochemistry 1987, 26, 6676-6683. (c) Carver, J. P.; Cumming, D. A. Pure Appl. Chem. 1987,59, 1465-1476. (d) Bock, K.; Josephson, S.; Bundle, D. R. J . Chem. Soc., Perkin Trans. 2 1982, 59-70. (e) Bock, K.; Lemieux, R. U. Carbohydr. Res. 1982, 100, 63-74. (f) Lemieux, R. U.; Koto, S. Tetrahedron 1974, 30, 1933-1944. (g) Lemieux, R. U.; Bock, K.; Delbaere, L. T. J.; Koto, S.; Rao, V. S. Can. J . Chem. 1980, 58, 631653. (h) Lemieux, R. U.; Wong, T. C.; Thorgersen, H. Can. J. Chem. 1982, 60, 81. (i) Lemieux, R. U.; Bock, K. Arch. Biochem. Biophys. 1983, 221, 125-134. (i)Thogersen, H.; Lemieux, R. U.; Bock, K.; Meyer, B. Can. J. Chem. 1982,60,44-57. (k) Alvarado, E.; Nukada, T.; Ogawa, T.; Ballou, C. E. Biochemistry 1991, 30, 881-886. (1) Peters, T.; Meyer, B.; StuikePrill, R.; Somarjai, R.; Brisson, J.-R. Carbohydr. Res. 1993, 238, 49-73. (24) Clark, M.; Cramer, R. D., III; van Opdenbosch, N. J. Comput. Chem. 1989, 10, 982-1012. (25) (a) Imberty, A.; Hardman, K. D.; Carver, J. P.; Pkrez, S. Glycobiology 1991, l, 631-642. See, for example: (b) Imberty, A.; Pkrez, S.; Hricovini, M.; Shah, R. N.; Carver, J. P. lnt. J. Biol. Macromol. 1993, 15, 17-23. (26) Rutherford, T. J.; Partridge, J.; Weller, C. T.; Homans, S . W. Biochemistry 1993,32, 12715- 12724. (27) Woods, R. J.; Edge, C. J.; Dwek, R. A. In Modelling the Hydrogen Bond; Smith, D., Ed.; Symposium Series 569; American Chemical Society: Washington, DC,1994; Chapter 16. (28) Woods, R. J.; Edge, C. J.; Wormald, M. R.; Dwek, R. A. In Complex Carbohydrates in Drug Research; Bock, K., Clausen, H., KrogsgaardLarsen, P., and Kofod, H., Eds.; Munksgaard: Copenhagen, Denmark, 1993; pp 15-26. (29) Woods, R. J.; Edge, C. J.; Dwek, R. A. Nature Struct. Biol. 1994, 1,499-501. (30) van Gunsteren, W. F.; Berendsen, H. J. C. Angew. Chem., Znt. Ed. Engl. 1990, 29, 992-1023. (31) Hall, D.; Pavitt, N. J. Comput. Chem. 1984, 5, 441-450. (32) McCammon, A. J.; Harvey, S. C. Dynamics of Proteins and Nucleic Acids; Cambridge University Press: Cambridge, UK, 1987; p 234. (33) Reynolds, C. A.; King, P. M.; Richards, G. W. Mol. Phys. 1992, 76, 251-275. (34) Komeiji, Y.; Uebayasi, M.; Someya, J.-I.; Yamato, I. Protein Eng. 1992, 5, 759-767. (35) Mimshima, N.; Spellmeyer, D.; Hirono, S.; Pearlman, D.; Kollman, P. J. Biol. Chem. 1991, 266, 11801-11809. (36) Ferguson, D. M.; Pearlman, D. A.; Swope, W. C.; Kollman, P. A. J. Comput. Chem. 1992, 13, 362-370. (37) Frisch, M. J.; Head-Gordon, M.; Trucks, G. W.; Foresman, J. B.; Schlegel, H. B.; Raghavachari, K.; Robb, M. A.; 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.; Stewart, J. J. P.; Topiol, S.; Pople, J. A. Gaussian 90; Gaussian, Inc.: Pittsburgh PA, 1990. (38) Frisch, M. J.; T N C ~ SG. , W.; Head-Gordon, M.; Gill, P. M. W.; Wong, M. W.; Foresman, J. B.; Johnson, B. G.; Schlegel, H. B.; Robb, M. A.; 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. (39) Hariharan, P. C.; Pople, J. A. Theor. Chim. Acta 1973, 28, 213222. (40) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J. Am. Chem. SOC. 1985, 107, 3902-3909. (41) Stewart, J. J. P. MOPAC 6.0; Frank J. Seiler Research Laboratory: United States Air Force Academy, CO 80840, 1990. (42) Pearlman, D. A.; Case, D. A.; Caldwell, J. C.; Seibel, G. L.; Singh, U. C.; Weiner, P.; Kollman, P. A. AMBER 4.0; University of California: San Francisco, 1991. (43) Hopfinger, A. J.; Pearlstein, R. A. J . Comput. Chem. 1984.5.486499.

J. Phys. Chem., Vol. 99, No. 11, 1995 3845 (44) Homans, S.W. Biochemistry 1990,29,9110-9118. (45) Woods, R. I.; Szarek, W. A.; Smith, V. H., Jr. J . Chem. Soc., Chem. Commun. 1991, 334-337. (46) Mgiller, C.; Plesset, M. S. Phys. Rev. 1934, 46, 618. (47) Pople, J. A.; Binkley, J. S.; Seeger, R.Znt. J. Quantum Chem. Symp. 1976, 10, 1. (48) Krishnan, R.; Pople, J. A. Znt. J . Quantum Chem. 1976, 14, 91. (49) Spitznagel, G. W.; Clark, T.; Chandrasehar, J.; Schleyer, P. v. R. J. Comput. Chem. 1983, 3, 363. (50) Deslongchamps, P. Stereoelectronic Effects in Organic Chemistry; Pergamon Press: New York, 1983. (51) Jeffrey, G. A.; Pople, J. A.; Radom, L. Carbohydr. Res. 1972, 25, 117-131. (52) Jeffrey, G.A.; Pople, J. A.; Binkley, J. S.; Vishveshwara, S. J. Am. Chem. SOC.1978, 100, 373-379. (53) Bowen, J. P.; Allinger, N. L. In Reviews in Computational Chemistry; Lipkowitz, K. B., Boyd, D. B., E&.; VCH Publishers: New York, 1991; Chapter 3. (54) Dinur, U.; Hagler, A. T. In Reviews in Computational Chemistry; Lipkowitz, K. B., Boyd, D. B., Eds.; VCH Publishers: New York, 1991; Chapter 4. (55) Maple, J. R.; Thacher, T. S.; Dinur, U.; Hagler, A. T. Chem. Design Automat. News 1990, 5, 5-10. (56) Maple, J.; Dinur, U.; Hagler, A. T. Proc. Natl. Acad. Sci. CJ.SA. 1988, 85, 5350-5354. (57) Andrews, C. W.; Fraser-Reid, B.; Bowen, J. P. J. Am. Chem. SOC. 1991, 113, 8293-8298. (58) Schleifer, L.; Senderowitz, H.; Aped, P.; Tartakovsky, E.; Fuchs, B. Carbohydr. Res. 1990, 206, 21-39. (59) Allen, F. H.; Davies, J. E.; Galloy, J. J.; Johnson, 0.; Kennard, 0.; Macrae, C. F.; Mitchell, E. M.; Mitchell, G. F.; Smith, J. M.; Watson, D. G.J. Chem. In& Compnt. Sci. 1991, 31, 187-204. (60) Ha, S. N.; Giammona, A.; Field, M.; Brady, J. W. Carbohydr. Res. 1988, 180, 207-221. (61) van Gunsteren, W. F.; Mark, A. E. Eur. J. Biochem. 1992, 204, 947-961. (62) Yan, Z.-Y.; Bush, A. C. Biopolymers 1990, 29, 799-811. (63) Scarsdale, J. N.; Ram, P.; Prestegard, J. H.; Yu, R. K. In Computer Modeling of Carbohydrate Molecules; French, A. D., Brady, J. W., Eds.; Symposium Series 430; American Chemical Society: Washington, DC, 1990; Chapter 15. (64)Krol, M. C.; Huige, C. J. M.; Altona, C. J . Comput. Chem. 1990, 11, 765-790. (65) Senderowitz, H.; Aped, P.; Fuchs, B. Helv. Chim. Acta 1990, 73, 2113-2128. (66) Reviews in Computational Chemistry; Lipkowitz, K. B., Boyd, D. B., Eds.; VCH Publishers, Inc.: New York, 1991. (67) Jeffrey, G. A.; Taylor, R. J. Comput. Chem. 1980, 1, 99-109. (68) Koehler, J. E. H.; Saenger, W.; van Gunsteren,W. F. Eur. Biophys. J . 1987, 15, 197-210. (69) Rush, C. A.; Blumberg, K.; Brown, J. N. Biopolymers 1982, 21, 1971-1977. (70) Delbaere, L. T. J. Biochem. J. 1974, 143, 197-205. (71) Shibata, S.; Nakanishi, H. Carbohydr. Res. 1980, 86, 316-320. (72) Lemieux, R. U.; Koto, S.; Voisin, D. In Anomeric Effect, Origin and Consequences; Szarek, W. A., Horton, D., Eds.; Symposium Series 87; American Chemical Society: Washington, DC, 1979; Chapter 2. (73) Tvaroska, I.; Kozar, T. Carbohydr. Res. 1981, 90, 173-185. (74) Tvaroska, I.; Carver, J. P. J. Chem. Res. 1991, 0123-0144. (75) Booth, H.; Grindley, T. B.; Khedhair, K. A. J. Chem. Soc., Chem. Commun. 1982, 1047-1048. (76) Booth, H.; Khedhair, K. A.; Readshaw, S. A. Tetrahedron 1987, 43,4699-4723. (77) Anomeric Effect. Origin and Consequences;Szarek, W. A., Horton, D., Eds.; Symposium Series 87; American Chemical Society: Washington, DC, 1979. (78) Kirby, A. J. The Anomeric Effect and Related Stereoelectronic Effects at Oxygen; Springer-Verlag: New York, 1983. (79) Praly, J.-P.; Lemieux, R. U. Can. J . Chem. 1987, 65, 213-223. (80) Wiberg, K. B.; Murcko, M. A. J. Am. Chem. SOC.1989,111,48214828. (81) Booth, H.; Khedhair, K. A. J . Chem. Soc., Chem. Commun.1985, 467-468. (82) de Hoog, A. J.; Buys, H. R.; Altona, C.; Havinga, E. Tetrahedron 1969, 25, 3365-3375. (83) Juaristi, E.; Cuevas, G. Tetrahedron 1992, 48, 5019-5087. (84) Booth, H.; Dixon, M.; Khedhair, K. A.; Readshaw, S. A. Tetmhedron 1990, 46, 1626-1652. (85) Wormald, M. R.; Wooten, W. E.; Bazzo, R.; Edge, C. J.; Feinstein, A.; Rademacher, T. W.; Dwek, R. A. Eur. J . Biochem. 1991, 198, 131139. (86) Davis, J. T.; Hirani, S.; Bartlett, C.; Reid, B. R. Biochemistry 1994, 269, 3331-3338. (87) Wolfe, S. Acc. Chem. Res. 1972, 5, 102-111. (88) Nishida, Y.; Ohrui, H.; Meguro, H. Tetrahedron Lett. 1984, 25, 1575-1578.

3846 J. Phys. Chem., Vol. 99,No. 11, 1995 (89) Nishida, Y.; Hori, H.; Ohrui, H.;Megum, H.Carbohydr. Res. 1987, 170, 106-111. (90) Perkins, S. J.; Johnson, L. N.; Phillips, D. C.; Dwek, R. A. Carbohydr. Res. 1977,59, 19-34, (91) Kroon-Batenburg, I. M. J.; Kroon, J. Biopolymers 1990.29, 12431248. (92) Ha,S.;Gao, J.; Tidor, B.; Brady, J. W.; Karplus, M. J . Am. Chem. SOC.1991,113,1553-1557. (93) Caminati, W.; Giorgio, C. J. Mol. Spectrosc. 1981,90,572-578. (94) Snyder, R.G.; Zerbi, G. Spectrochim. Acta 1%7,23A, 391-437. (95) Connor, T. M.; McLaughlin, K. A. J . Phys. Chem.1%5,69,18881893. (96) Wiberg, K.B.; Murcko, M. A. J. Mol. Struc. (THEOCHEM) 1988, 163,1-17. (97) Wiberg, K. B.;Murcko, M. A. J. Phys. Chem. 1987,91,36163620. (98) Wiberg, K.B.; Murcko, M. A.; Laidig, K. E.; MacDougall, P. J. J . Phys. Chem. 1990,94,6956-6959. (99) Woods, R. J.; Khalil, M.; Pell, W.; Moffat, S. H.; Smith, V.H., Jr. J . Cornput. Chem. 1990,11,297-310. (100)Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990,11, 361-373. (101)Stouch, T. R.;Williams, D. E. J. Comput. Chem. 1992,13,622632. (102) Williams, D. E. Biopolymers 1990,29, 1367-1386. (103) Woods, R. J.; Andrews,C. W.; Bowen, J. P. J. Am. Chem. SOC. 1992,114,850-858. (104)Chirlian, L. E.; Francl, M. M. J. Comput. Chem. 1987,8, 894905.

Woods et al. (105) Jorgensen, W. L.;Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Phys. Chem.1983, 79,926-935. (106)Kimura, K.;Kubo, M.J. Chem.Phys. 1959,30,151-158. (107)Dannenberg, J. J.; Evleth, E. M. lnr. J . Quantum Chem. 1992,

44,869-885. (108) Kim, S.:Jhon, M. S.;Scheraga, H.A. J. Phys. Chem. 1988,92, 7216-7223. (109)hn-Batenburg, L.M. J.; van Duijneveldt, F. B. J . Phys. Chem. 1986,90,5431-5435. (110)Tse, Y.-C.; Newton, M. D.; Allen, L. C. Chem. Phys. Le#. 1980, 75,350-356. (111) Del Bene, J. E. J. Chem.Phys. 1971,55,4633-4636. (112) Zheng, Y.4.; Me&, K. M., Jr. J . Comput. Chem.1992,13,11511169. (113) Pimentel, G. C.; McClellan, A. L. The Hydrogen B o d , W. H. Freeman and Company: New York, 1960,p 475. (114) Scheiner, S.In Reviews in Computational Chemistry; Lipkowitz, K. B., Boyd, D. B., Eds.;VCH Publishers: New York, 1991;Chapter 5. (115) Ryckaert, J.-P.; Ciccott, G.; Berendsen, H. J. J . Comput. Phys. 1977,23, 327-341. (116) Jorgensen, W. L. J. Am. Chem. SOC.1981,103,335-340. (117)Brady, J. W.; Schmidt, R. K. J. Phys. chem. 1993,97,958-966. (118) Edge, C. J.; Singh, U. C.; Bazzo, R.; Taylor, G. L.; Dwek, R. A.; Rademacher, T. W. Biochemistry 1990,29, 1971-1974. (119)Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1984,5, 129145. (120) Kollman, P.A. Personal communication. Jp9421803