J. Phys. Chem. 1993, 97, 958-966
958
The Role of Hydrogen Bonding in Carbohydrates: Molecular Dynamics Simulations of Maltose in Aqueous Solution J. W.Brady’ and R. K. Schmidt Department of Food Science, Stocking Hall, Cornell University, Ithaca, New York 14853- 7201 Received: October 5, 1992
Molecular dynamics simulations of maltose in two different conformations in vacuum and aqueous (TIP3P) solution have been used to examine the types of hydrogen bonds made by this carbohydrate molecule. Maltose was found to be extensively hydrogen bonded to solvent molecules in aqueous solution, and these hydrogen bonds were found to have potential conformational consequences. The exchange of an intramolecular 02-03’ hydrogen bond found in the crystal structure for hydrogen bonds to solvent was observed to produce significant changes in the solvation of the sugar molecule. Solvation was found to increase root mean square conformational fluctuations, in contrast to the results found in a previous simulation of a mannose dimer in aqueous solution. A substantial fraction of the first hydration layer water molecules around the polyfunctional maltose solute was found to make simultaneous hydrogen bonds to more than one sugar hydroxyl group. The nature and extent of the conformational effects of solvation in model systems were found to be sensitive to the molecular mechanics force fields used.
potential energy functions (such as the so-called HSEA model).8 Perhaps the most seriouslimitation of much previous carbohydrate work was the inadequate treatment or even complete neglect of hydrogen bonding, along with the absence of solvent water molecules. Modern simulation methods and currently available computers no longer require such approximations,however, and molecular dynamics simulations of carbohydrates in aqueous solution are becoming increasingly As a result, it is now possible to determine directly the nature of hydrogen bonding in model carbohydrate solutions. We report here such MD simulations of the aqueous solvation of maltose and an analysis of the hydrogen bonding which occurs in this system. Maltose (Figure 1) is the a(l+4)-linked dimer of D-glucose and, as the basic repeat unit of amylose, is an important prototypical molecule. The crystal structure of &maltose has been determined to high resolution by neutron diffraction.16 In addition, there is some limited NOE data available for maltose in aqueous solution,l7 and its conformation in solution has been analyzed by optical rotation measurements.18 Maltose has long been the subject of molecular mechanics studies using simplified potential energy functionsand/or rigid monomer ~ n i t s ; ~more J~,~ recently, relaxed conformational energy maps for maltose have been prepared using several different force fields.21*22 The crystal conformation has been found to be stabilizedby an intramolecular hydrogen bond between the 0 2 hydroxyl group of the nonreducing residue and the 03’ hydroxyl group of the reducing residue.16 The relaxed conformational energy studies have generally found that maltose has several low energy conformations in vacuum, some of which possessthe intramolecular 02-03’ hydrogen bond, but with at least one important low-energy structure which, for geometric reasons, cannot form this hydrogen bond. In the relaxed vacuum energy maps for maltose,22 the hydrogen bonded conformers have been labeled ‘B”and ‘C” and the non-hydrogen bonded conformer is referred to as the conformation (Figures 2 and 5, and Figure 3 of ref 21). The fate of this intramolecular hydrogen bond in aqueous solution, where alternate hydrogen bond partners are provided, could have an important influence on thesolutionconfonnation of thecarbohydrate. Both theoptical rotation studies and the reported NOE data imply that in solution there is a shift away from the crystal structure (the C conformationl6) to the A conformation which does not contain the intramolecular hydrogen bond.I7J8
I. Introduction Hydrogen bonds are often quite important in determining the conformationsof many types of biopolymers. These interactions, usually 3-6 kcal/mol in magnitude, are much weaker than typical chemical bonds; however they are at least an order of magnitude larger than the strongest van der Waals interactions, and substantially larger than kTat room temperature. This strength makes them quite resistant to thermal disruption. Furthermore, these interactions are both short range and angle-dependent.’ This combination of strength, directionality, and short range of influence allows hydrogen bonds to play a significant role in structuring polymers that contain hydrogen-bondinggroups. The various secondary structures of proteins, such as 8-sheets, a-helices,and reverse turns are well-known examplesof structures dominated by hydrogen bonding. Hydrogen bonding might also be expected to be important in the conformational behavior of carbohydrates,2since they contain many hydroxyl groups which can serve both as hydrogen bond donors and acceptors. For a variety of reasons, the conformations of carbohydrates are less well understood than are those of peptides. Efforts to obtain diffracting single crystals of polysaccharides for crystallographic analysis have failed, and results from fiber diffraction studies have often proven ambiguous, requiring additional information for their interpretation.3 Although crystallographic experiments have produced high resolution structures for a variety of disaccharides, including the basic repeat units of many larger polymers, it is far from certain that these molecules adopt the same conformations in aqueous solution. While useful, NMR NOE studies of carbohydrates generally do not produce enough distance constraints to allow a complete structure determination. Recent limited NMR data for several important disaccharides in aqueous solution have also sparked controversies concerning whether the molecules adopt a single conformation in solution or whether conformational averaging between multiple conformers occum4 Computer modelings of carbohydrate conformational behavior could potentially be quite useful in helping to resolve some of these experimental uncertainties.6 Unfortunately, the results of such molecular mechanics calculations have also proven to be ambiguous. Often carbohydrate calculationshave been limited by severe approximations, including the neglect of internal flexibility7and the use of primitive To whom correspondence should be addressed. (B
1993 American Chemical Society
The Journal of Physical Chemistry, Vol. 97, No. 4, 1993 959
The Role of Hydrogen Bonding in Carbohydrates
I
I
I
I
I
I
5.0
10.0
I 15.0
120.0
\\
W
60.0
d
0.0 -60.0 -120.0
b
-
Time (ps)
03& 02'
Figure 1. &Maltose in the C conformation.
120.0
I
120 -
l
160
l
-
60.0
II/
0.0 -60.0
- 120.0
80
c t
I 40
I 5.0
1
10.0
I 15.0
Time (ps)
Figure 3. History of C#J and $ for the first 20 ps of data collection in the A-well vacuum trajectory.
0
TABLE I: Partial Atomic Charges for @-Maltose
-40
-a0
- 120 - 160 -120
-80
-40
0
40
80
120
$ Figure 2. Two vacuum trajectories starting from the A well and the C well superimposed on the adiabatic energy map. Contours are in 2 kcal/ mol intervals from 2.0 to 10.0 kcal/mol above the energy minimum.
II. Procedures In the present studies, several molecular dynamics simulations were performed for a single &maltose molecule surrounded by 488 water molecules using the molecular mechanics program CHARMM.23 The intramolecular energy of the solute was represented by a preliminary CHARMM-type force field proposed for carbohydrate^.^^ This all-atom potential function contains no term to explicitly represent the exo-anomericeffect,25although a revision of this force field to include such a contribution has recently been proposed.15 The energy function was thus the same as that used in the previously reported adiabatic energy mapping for this molecule.21 Although some empirical force fieldscontain a specific term to represent hydrogen bonding, this type of interaction can be modeled adequately by appropriatelyadjusting the atomic partial charges and van der Waals parameters for hydrogen bond donor and acceptor atoms.Z6 This approach is consistentwith the way in which hydrogen bonding is represented in most empirical water force fields, including the TIP3P model used here.Z7 The charges used for each of the solute atoms are listed in Table I. An additional simulation in solution was also computed using the alternate charge set included in Table I. Each simulation was initiated by superimposing previously identified minimum-energy vacuum conformations for maltose2' upon the center of an equilibrated cubic box containing 512 TIP3P27water molecules. Those water molecules with oxygen atoms that overlapped with any of the maltose heavy atoms were eliminated, leaving 488 solvent molecules, and a weight con-
atoms in atoms in nonreducing charge reducing charge residue standard alternate residue standard alternate 0.25 0.35 C1' 0.3 c1 0.3 0.1 H1' 0.1 0.1 H1 0.1 -0.55 -0.65 01' -0.4 01 -0.4 0.4 0.4 OH1' 0.05 0.15 C2' 0.05 c2 0.15 0.1 0.1 H2' H2 0.1 0.1 -0.55 -0.65 02' -0.55 02 -0.65 0.4 0.4 OH2' 0.4 OH2 0.4 0.15 0.05 0.05 C3' c3 0.15 0.1 0.1 H3' 0.1 0.1 H3 -0.55 -0.65 -0.55 03' 03 -0.65 0.4 OH3' 0.4 0.4 0.4 OH3 0.1 0.1 C4' 0.15 0.05 c4 0.1 0.1 H4' 0.1 H4 0.1 04 -0.65 -0.55 0.4 0.4 OH4 0.1 0.1 C5' 0.1 0.1 c5 0.1 0.1 H5 H5' 0.1 0.1 -0.4 -0.4 05' -0.4 05 -0.4 0.05 -0.05 -0.05 C6' C6 0.05 0.1 0.1 H6' 0.1 H6 0.1 0.1 0.1 H7' 0.1 H7 0.1 -0.55 -0.65 -0.55 06' -0.65 06 0.4 0.4 OH6' 0.4 0.4 OH6
centration of 3.75%. The box length was then adjusted to 24.64 A to give a density of 1.013 g/cm3, the approximate experimental density of an aqueous maltose solution of this weight percentage. The concentration of the resulting solution was 0.12 M. The box was subjected to minimum-image periodic boundary conditions to eliminate edge effects. Long-range interactions were adjusted to smoothly approach zero by using a switching function applied on a group-by-groupbasis between 10.0 and 12.0A.23 The groups corresponded to entire water molecules and electrically neutral collections of adjacent atoms in the solute. Initial velocities for all atoms in these systems were selected at random from a Boltzmann distribution at 300 K. The equations of motion were then integrated as a microcanonical ensemble using a Verlet integrator28 with a step size of 1 fs. The water molecules were kept rigid using the SHAKE constraint which was also used to fix the bond lengths for solute bonds involving
Brady and Schmidt
960 The Journal of Physical Chemistry, Vol. 97, No. 4, I993 hydrogen atoms. Each simulation was equilibrated for 20 ps, during which velocities were reassigned or rescaled if the average temperature in a 2-ps interval deviated from 300 K by more than 3 K. Following equilibration, the trajectories were continued without further interference, and complete phase space points were saved every 10 fs for subsequent analysis. A total of two vacuum and three solution simulations were conducted. One vacuum trajectory was calculated with the starting structure for the maltose in the A conformation (using the nomenclature of Tran et a1.22), and the second vacuum trajectory was started in the C conformation. Both trajectories consisted of 30 ps of equilibration and 80 ps of data collection. Two primary solution simulations were conducted. The first used as its starting structure maltose in the A conformation; this simulation was equilibrated for 20 ps and was then continued for an additional 80 ps of dynamics. The second primary solution trajectory was started from the C conformation and comprised 20 ps of equilibration and 50 ps of data collection. An additional solution simulation was started in the C conformation and employed thealternate set of maltose atomic partial charges listed in Table 1. This simulation comprised 20 ps of equilibration and 40 ps of data collection. The main difference in this alternate set of charges was a reduction in the charge of the hydroxyl oxygen. These short simulation times were not adequate to allow a conformational equilibrium to develop between the wells but were sufficientfor convergence of "V"-state averaging30of solvent properties in each well. The length of the simulations were limited by computer speed; the A-well solution simulation, for example, required 447 h of CPU time on an IBM RISC/6000 computer.
III. Results and Discussion Figure 2 displays trajectories for the two different vacuum dynamics simulations projected onto the vacuum adiabatic potential energy map for this energy function.21 The energy is plotted as a function of the glycosidic linkage torsion angles 4, defined as Hl-Cl-Ol-C4', and $, defined as C1-01-C4'-H4'. These new vacuum trajectories were qualitatively like those calculated previously for this molecule.21 In the absence of alternate hydrogen bond partners provided in solution by water molecules, the maltose hydroxyl groups are forced to make internal hydrogen bonds to satisfy their hydrogen-bonding requirements. In general, each hydroxyl group orients so as to make a hydrogen bond to an adjacent hydroxyl on the ring, with each of these hydrogen bonds pointing in the same direction around the sugar rings (see Figure 1). The crystallographically observed 02-03' hydrogen bond must be broken in any transition from the B or C well to the A well, requiring reorientation of the other internal hydrogen bonds, thereby contributing to a barrier to such transitions in vacuum. In both of the illustrated simulations, the conformational space within the low-energy contours is wellsampled by the vacuum trajectory, and the adiabatic energy surface appears to be an excellent description of the likely motions of the molecule in vacuum. Figure 3 displays the individual histories of 4 and $ for motion in the A well as calculated from the vacuum trajectory. For both angles,there is a relatively simple low-frequency oscillationabout the mean values, with similar periods. Figure 4 illustrates the correlation functions for fluctuations in these two angles, which are sufficiently simple in form that one can easily estimate quantitatively the primary frequencies revealed in the power spectra of these correlation functions,displayed in the insets. The peak in the power spectrum for 4 is at 4.7 ps-I, and the peak for $ is at 6.3 ps-I. The weak damping evident in the slow decay of these correlation functions presumably results primarily from long-range nonbonded interactions between the two residues. However, for both angles the significant feature of the motion is long-term correlation due to the weakness of interactions with the environment, which in vacuum consists solely of the remainder
100
-025
I
I
I
025
050
075
,
I
I
I
100
125
1.50
171
-
-
Time (P)
WO IWO Is00 2WO 2500 Y
-OZ5
tI
(P."
)
v I
I
I
I
I
I
I
I
025
050
075
I W
125
150
175
I
Time (ps)
Figure 4. Correlation functions and spectral densities (insets) of fluctuations in (a, top) $I and (b,bottom) in the A-wellvacuumtrajectory.
+
of the molecule. Both histories and correlation functions also exhibit a high-frequency, small-amplitude oscillation which is a result of defining the torsion angles in terms of the instantaneous positions of four specific atoms. C-0-C and H-C-O bond bending motions, which will also cause these atoms to move, can thus register by this definition as small changes in 4 and $, with corresponding peaks in the power spectra around 25Ops-I. Similar correlation functions and power spectra were obtained for the C-well simulation. The presence of solvent molecules can affect the behavior of the maltose solute in several ways. The "bath" of surrounding molecules should exert a frictional damping on both the overall motions of the molecule and on its internal conformational fluctuations. The solvent might also induce a conformational shift in the molecule, either locally or by causing a transition to another conformation altogether, and the magnitude of the fluctuations of the solute motions might also be affected. Figure 5 displays the trajectories for the two primary solution simulations, plotted on the vacuum adiabatic energy map, and Table I1 lists the mean values of the linkage torsion angles for these two trajectories and their root mean square (rms) fluctuations. Both trajectories were stable in their respective wells for the entire period of the simulations, with no fluctuations bringing the molecules close to apparent transition states for conformational changes. It is not possible toevaluate the relative stabilities (free energies) of the two wells from such short simulations, but it would appear from the lackof transitionsthat the barriers between them are moderately large in solution with this energy function (see below). As can be seen from the table, solvation produced only small conformational shifts within the wells. In contrast to previously reported results for another disaccharide, Mana( 1+2)Mana, in solution14 the rms fluctuations in 4 and $ were not found to be smaller than in vacuum for either simulation but were actually larger. In both conformations, the moleculeexplores a larger region of conformational space as a result of interaction
The Role of Hydrogen Bonding in Carbohydrates 180 120
I
-
l
l
-
The Journal of Physical Chemistry, Vol. 97, No. 4, 1993 961 I
-
-
I
I
I
I
I
I
I
I
120.0 60.0
qj
0.0
-80.0 -120.0
t20.0
40.0
60.0
60.0
Time (pa)
120.0
,
nn -.-
-60.0
-120
-80
-40
40
0
80
120 -120.0
$ Figure 5. Two solution trajectories superimposed on the adiabatic energy map.
simulation
mean
vacuum A wellz1 vacuum C we1121 present vacuum A well present vacuum C well solution A well solution C well
-58.9 32.1 -57.5 32.7 -62.3 28.5
8.5 11.8 10.0 11.5 10.5 16.0
5.0
I .o
mean -43.8 15.5 -44.2 15.7 -49.4 19.8
8.3 8.1 8.3 8.7 11.0 10.7
1
F 20 0
60 0
400
60 0
80.0
Figure 7. History of 0 and $Jfor the A-well solution simulation over the entire trajectory, including equilibration. rms fluctuation
t
40.0
Ti” (p.)
TABLE Ik Ensemble-Averaged Mean Values and Fluctuations for 4 and 9 rms fluctuation
20.0
800
Tunc (pa)
Figure 6. History of the distance between 0 2 and 03’ for the C-well solution simulation over the entire trajectory, including equilibration.
with the solvent. In addition, in vacuum, the C conformation is constrained by the intramolecular 02-03’ hydrogen bond, but in aqueous solution water molecules can compete as alternate hydrogen bond partners for these hydroxyl groups. Figure 6 displays the history of this 02-03’ hydrogen bond, which generally persisted throughout the simulation. However, brief disruptions of this hydrogen bond did occur, the longest lasting about 7 ps, and these disruptions presumably contributed to the increase in fluctuations in 4 and $ for the C-well trajectory. For the most part, however, the internal hydrogen bond appears to have been favored over the mixed solvent-solute hydrogen bonds which temporarily replaced it. These results are almost certainly a function of the choice of empirical force field, since parameters which allow the water molecules to effectively compete with the internal hydrogen bond will allow larger fluctuationsin the torsion angles or even produce conformational shifts (to the A well, for example). The mannose disaccharide modeled by Edge et al.I4 does not possess any hydrogen bonds between the two rings due to its
stereochemicalstructure. Thus, hydrogen bond competition from the solvent would not necessarily be expected to affect the conformational equilibriumfor that molecule. However, solventsolute hydrogen bonding would probably not have been significant in that simulation regardless of stereochemistry,due to the choice ofatomic partialcharges. In that study, theatomic partialcharges were taken directly, without any scaling,from an AM 1calculation, which are known to produce charges inappropriately small for solution simulations.31 These hydroxyl oxygen charges, which ranged from -0.3e to -0).34e,are less than half the values used in the present study, and substantially smaller than the water oxygen charge. These charges are so small that the hydroxyl groups would not be expected to hydrogen bond to the solvent (see below). Although these authors report no detailed information concerning the water structure in their study, it is possible that the carbohydrate solute may have behaved essentially as an hydrophobic species. Since the overall conformation in that simulation agrees satisfactorily with the NOE data, it would appear that solvent does not significantly perturb the conformational equilibrium in this mannose disaccharide, perhaps due to the absence of intramolecular hydrogen bonds. However, the de-emphasis of solvent-solute hydrogen bonds may have sufficiently weakened the coupling between the carbohydrate and solvent resulting in the small rms fluctuationsin the linkage torsion angles. While in the present solution simulations the magnitudes of conformational fluctuations are larger than in vacuum, their temporal correlation is stronglydamped by friction with theviscous solvent. Figures 7 and 8 display the histories of 4 and 9 calculated from the simulation in the A well and the correlation functions for fluctuations in these angles, along with their power spectra. The low-frequency oscillations are strongly damped, and the correlation functions exhibit a slow overdamped decay, with a substantial zero-frequency component as the only important peak in the power spectra. There was little difference in this relaxation behavior between the solution simulations in the two conformational wells. In these simulations, the majority of the hydrogen bonds made by the maltose molecules were with the solvent water. Figure 9 displays the pair distribution function for water molecule oxygen atoms around the 0 3 hydroxyl oxygen of the nonreducing group of the C-well simulation. As can be seen, this function, which
962
Brady and Schmidt
The Journal of Physical Chemistry, Vol. 97, No. 4, 1993 I
OZ5
-OZ5
,
I
1
,
I
I
t
1
t
i 025
OW
I
I
075
100
125
I50
I
I
I
00
-0.5
0.5
1.0
Cos(@)
Figure 10. Distributionoforientations for those water moleculesadjacent to 0 3 calculated from the data collection period of the C-well solution
175
Time (ps)
I
-1 0
I
simulation.
TABLE IIk Avera e Number of SoluteSolvent Hydrogen for S e v e d ddtose ~ o ~ u t i o~imuiatiow n
I 00
0 75
0 50 0 25
Oo0
R atom
A well geointemetric grated 1.237 0.365 4.388 2.709 4.607 2.539 3.817 2.580 1.831 0.704 4.202 2.836 4.5 19 2.800 4.323 2.549 3.592 2.282 2.482 0.722 3.605 2.574
C well geometric
integrated
' 050
025
IW
075
125
150
175
Time (ps)
Figure 8. Correlation functions and spectral densities (insets) of fluctuations in (a, top) and 6 and (b, bottom) t j in the A-well solution simulation.
2.0
1.5
eW 1.o
0.5
4.0
2.0
r
6.0
(4
Figure 9. Water oxygen-hydroxyl oxygen 0 3 pair distribution function calculated from thedata collection period of theC-well solution simulation.
was typical of these distributions for hydroxyl groups in all of the simulations, displays characteristic hydrogen-bonding behavior, with a sharp, narrow first peak around 2.8 A and a deep first minimum around 3.6 A. It is well-known however, that theTIP3P water model used in the present simulations exhibits less 'structuring" than other models,27such as the ST2 or SPC m o d e l ~ . ~As~aJresult, ~ the pair distribution functions calculated from these simulations also display somewhat less structuring than observed in a previous simulation of D-glucose in SPC solution,'O where the first maxima were at 2.7 A and the first minima generally were around 3.4 A. The integral of the curve in Figure 9 out of 3.5 Agives 4.2 nearest-neighbors, substantially
01 02 03 04 05 06 01' 02' 03' 05' 06'
0.285 2.164 2.610 2.396 0.933 2.999 2.787 2.976 2.209 0.453 2.649
3.746 4.21 1 4.150 2.130 4.791 4.325 4.585 3.707 3.91 1
C well (alternate charges) integee: metric grated 0.364 2.333 2.472 1.945 0.736 2.243 2.726 2.355 2.095 1.044 2.723
3.533 4.205 3.480 1.636 3.981 4.338 4.356 3.686 2.365 4.125
The first column of values for each simulation was calculated using the geometric criteria mentioned in text, and the second column of values was calculated by integrating the pair distribution function out to the average first minimum of 3.5 A. a
more than the three observed in our previous glucose simulations in SPC water.1° This increase results from first neighbors crowding in around the hydroxyl group at less than optimal angles of hydrogen bonds. Figure 10 displays the orientational distribution function for the nearest-neighbor water molecules around this same hydroxyl group; as can be seen, there is substantially less tetrahedral ordering of these neighboring solvent molecules than in the previous simulation of the glucosbsPC solution (Figure 13 of ref 10). Table 111lists the number of near neighbors for each OH group in the maltose solutes as calculated by integration of the pair distribution functions from both the Aand C-well simulations. The pair distribution function in Figure 9 also exhibits a prominent second peak around 5.7 A which was also typical of these functions for the maltose hydroxyl group. This peak is at a much larger distance than the weak second peak in the pure TIP3P water and the first minimum is much deeper than t h a t for pure TIP3P water. This second peak apparently results from the first shell water molecules around the adjacent sugar hydroxyl groups, and the deeper minimum is caueed by partial exclusion of water molecules from the regions between the hydroxyl groups. This additional structuring thus apparently results from the mutual interference of the hydration behavior of adjacent functional group in a way which would not arise in monofunctional solutes such as methanol. Since not all of the near neighbors contributing to the first hydration shell around a hydroxyl group using this water potential function are in classic hydrogen-bonding geometries, it is useful to calculate the number of such hydrogen-bonding near neighbors using a geometric definition. Table 111 also lists the average number of hydrogen bonded water molecules for each hydroxyl
The Journal of Physical Chemistry, Vol. 97, No. 4,1993 963
The Role of Hydrogen Bonding in Carbohydrates TABLE Iv: Average Number of Solute-solrent Hydrogen Bonds Calculated from the Data Collection Period of the Solution Simulation Starting from the C Well Conformation (a) Number of SolutbSolvent Hydrogen Bonds as a Function of Distance Cutoff with an Angular Cutoff of 120' atom
0 1 02 03 04 05 06 01' 02' 03' 05' 06'
rm = 3.4 A 0.25 2.05 2.48 2.29 0.86 2.87 2.66 2.85 2.12 0.42 2.54
rm = 3.5 A 0.29 2.16 2.6 1 2.40 0.93 3.OO 2.79 2.98 2.2 1 0.45 2.65
rm = 3.6 A
rm = 3.7 A
0.3 1 2.27 2.73 2.50 1 .oo 3.12 2.9 1 3.10 2.30 0.49 2.74
0.33 2.36 2.87 2.62 1.07 3.25 3.07 3.23 2.39 0.53 2.84
(b) Number of SolutbSolvent Hydrogen Bonds as a Function of Angular Cutoff with a Distance Cutoff of 3.5 A atom
0 1 02 03 04 05 06 01' 02' 03' 05' 06'
8=
loo0
0.53 3.56 3.85 3.73 1.45 4.60 4.18 4.50 3.45 0.92 3.68
8 = 120'
0.29 2.16 2.61 2.40 0.93 3.00 2.79 2.98 2.21 0.45 2.65
8 = 140° 0.13 1.09 1.46 1.29 0.40 1.56 1.57 1.63 1.1 1 0.18 1.64
group in the maltose molecules using such a geometric definition which required an oxygen-oxygen distance of less than 3.5 A and a donor oxygen-hydrogen-acceptor oxygen angle of greater than 120°. Using thesegeometriccriteria, those hydroxyl groups fully exposed to the solventand not engaged in intramolecular hydrogen bonds were indeed found to make approximately three hydrogen bonds to the solvent as expected. These geometric criteria for a hydrogen bond are somewhat arbitrary, despite being based on the concept of a hydrogen bond as a relatively linear first-neighbor interaction at a distance less than the minimum in the pair distribution function. In order to test the sensitivity of the calculated number of hydrogen bonds to the values of these cutoffs, Table IV presents the number of water-solute hydrogen bonds calculated for each solute group in the A-well simulation as a function of cutoff distance and angle. As can be seen, the number of near neighbors is not strongly sensitive to small changes in the cutoff distance criterion but is very sensitiveto the angular cutoff. This data and theorientational distribution functions such as Figure 10 confirm that the larger number of near neighbors in the pair distribution functions in these TIP3P simulations compared to our previous SPC calculations results from nonlinear hydrogen-bonding water molecules crowding around the solute hydroxyl group in the generally less structured TIP3P solvent. It would be desirable to be able to examine the number of water-maltose hydrogen bonds in terms of an energetic definition, as is often done in the case of water-water hydrogen bonds. Unfortunately, such a definition is difficult to define in any meaningful way for such an extended solute. In the case of two water molecules, one can take the total nonbonded interaction energy between the two molecules as the hydrogen bond energy, but for a soluteas large as maltose, such a definitionwould include atoms more than 10 A away, which clearly have little to do with the local nature of hydrogen bonding. However, using a nonneutral subset of atoms from the solutecan substantially affect the computed hydrogen bond energies, making comparisons between different hydroxyl groups in the solute problematic, as the computed hydrogen bond energy will oscillate significantly
(C)
(B)
(A)
(0)
figure 11. Schematic depiction of four possible double hydrogen bond arrangementsbetween solute hydroxylsand a water molecule: (A) two solute donor-water acceptorhydrogen bonds; (B) one solutedonor-water acceptor and one solute acceptor-water donor hydrogen bond; (C) two solute acceptor-water donor hydrogen bonds (involving two different water hydrogens); (D) two solute acceptor-water donor hydrogen bonds (involving only one water hydrogen).
TABLE V: Percent of All SoluteSolvent Hydrogen Bonds Classified as One of the Double Hydrogen Bond Arrangements Depicted in Figure 11, Calculated from the Data Collection Periods of the Solution Simulations
type
percent of all solute-solvent hydrogen bonds C well A well C well (alternate charges)
~~
A
D
3.32 13.58 5.62 15.29
total
30.46
B C
2.66 8.69 4.44 14.67 32.87
3.83 7.99 4.33 16.72 37.81
with the number of included atoms.34 Particularly in a large solutewith adeepcleft, such as the region between the twoglume rings, one water molecule could also be interacting with more than one hydroxyl group, obscuring the meaning of such an energetic definition. For these reasons, no unambiguous energetic criterion for water-solute hydrogen bonds could be developed. van Eijcket al.13 have studied the hydrogen bonding of hydroxyl groups from similar MD simulations of glucose and fructose in SPC solution and reported finding fewer hydrogen bonds for each hydroxyl than the three found here and in our previous simulation.IO This result is surprisingsince each hydroxyl group should be able to form two hydrogen bonds as an acceptor and one as a donor. The reason for this differenceapparently arises primarily from thedifferent definitions of hydrogen bonds employed by the two groups. In our previous paper, the number of hydrogen bonds was calculated by integrating the pair distribution function out to the first minimum. (This nongeometric definition apparently worked well for the SPC solution but, as seen above, is less satisfactory for the less structured TIP3P solutions.) In their analysis, van Eijck et al. used a geometric criterion similar to the one employed here except that it excluded the possibility that water molecules might be simultaneously hydrogen bonded to more than one hydroxyl group. In their classification scheme, each water was considered to be hydrogen bonded to the hydroxyl group to which it was closest, even if it was also simultaneously close enough to another group to make a second hydrogen bond. Such shared hydrogen bonds were found to arise frequently in the present maltose simulations, however. Hydrogen bonds in the present study were analyzed using geometric criteria which considered all such interactions that a solvent molecule might make. Using this definition, a variety of doubly hydrogen bonded arrangements, illustrated schematically in Figure 11, were observed, and their percentage of the total number of maltosesolvent water hydrogen bonds are given in Table V. As can be seen from the table, in the present simulationsthesevarious shared water molecules account for approximately one-third of all of the hydrogen bonds between the solute and water molecules, which accounts for the discrepancy between the number of hydrogen bonds counted by the two methods. The most common type of
964
Brady and Schmidt
The Journal of Physical Chemistry, Vol. 97, No. 4, 1993
120.0
x
0.0
tI-
1-I
I
-
1
20.0
40.0
60.0
I
80.0
1.0
1
Time (pa)
Figure 12. History of dihedral angle C4-CS-C6-06 for the C-well
solution simulation over the entire trajectory, including equilibration. shared hydrogen bond in these simulations was type D, in which a single proton from a water molecule bridges two hydroxyl oxygen atoms. Not all of the hydroxyl groups in the present simulations formed three hydrogen bonds to solvent, as can be seen from Table 111. In the simulation of the C-well geometry, both the 0 2 and 03' hydroxyl groups formed only about two hydrogen bonds to solvent molecules. This smaller number can readily be understood, however,since these twogroups were engaged in an intramolecular hydrogen bond between the rings during the majority of the simulation, as noted above. In the A-well geometry, the02 group did indeed form closer to three hydrogen bonds to solvent,although the number of solvent bonds for the 03' group was still small. The number of hydrogen bonds to water molecules for hydroxyl group 04 was also smaller than expected in the simulation of the C-well geometry. This smaller number may also be the result of an intramolecular hydrogen bond. Due to the rotational freedom of the C4-C5-C6-06 torsional angle, there are three possible orientations of the primary alcohol groups of both rings, usually designated GT, GG, and TG,3Sspecifying the position of the06 atomrelative to theringoxygenandC4atoms,respectively. In the crystal structure, the primary alcohol of the nonreducing residue is in the GT conformation, and that of the reducing residue is in the GG conformation. The TG conformation places the 0 6 and 0 4 hydroxyl groups sufficiently close to form an intramolecular hydrogen bond, and for this reasons this conformation is favored in vacuum calculations for the nonreducing ring. However, NMR experiments indicate that this conformation is not the dominant one in aqueous solution.36 It might again be expected that this conformation will be sensitive to relative hydrogen bond strengths, and previous simulations of glucose in solution have found that this conformational equilibrium is sensitive to the force field e m p l ~ y e d . ~For ~ ~the ~ ,maltose ~~ molecules considered here, the orientation of this primary alcohol group for the reducing residue will not generally be significantly affected by any tendency to form an intramolecular hydrogen bond, since its 04 group is involved in the glycosidic linkage and thus has a much lower partial charge. All of the simulations described here employed a starting structure with the primary alcohol of the nonreducing residue in the TG conformation, the lowest energy vacuum conformer. In the C-well simulations, this group rotated to the GT conformer and remained for the majority of the simulation, as shown in Figure 12, but returned to the TG conformation for the final 15 ps. The intramolecular hydrogen bond during this final interval thus reduced the mean number of hydrogen bonds to water for the 04 group. In the reducing residue, this exocyclic group remained in the GG conformation throughout. In the simulation of the A-well geometry, this primary alcohol group on the nonreducing residue remained in the TG conformation for 50 ps. It then made a transition to the GT conformation for 30 ps, followed by a transition to the GG conformation for the remainder of the simulation. The primary alcohol group on the reducing residue stayed in the GG conformation throughout the simulation.
2.0
4.0
6.0
r (A)
Figure 13. Water oxygen-linkage oxygen 01 pair distribution function calculated from the data collection period of the A-well solution simulation.
Since so few transitions were observed, it would seem likely that much longer simulations would be needed to calculate meaningful relative probabilities for these different conformers. A third type of intramolecular hydrogen bond is possible in the maltose solutes in these simulations. These are the strained hydrogen bonds between adjacent hydroxyl groups around the periphery of each ring, which tend to point in the same direction to avoid steric clashes (see Figure 1). These near neighbor hydrogen bonds have been observed in vacuum calculations,2I where there are no alternate hydrogen bond partners available, and have also been observed in NMR studies of carbohydrates in dimethyl sulfoxide (DMSO) solution.38Little evidence exists, however, to indicate that they are present in aqueous solution. In all of the simulations reported here, adjacent hydroxyl groups on both rings made hydrogen bonds to one another for extended periods, although numerous transitions occurred which resulted in the exchange of donor and acceptor roles between the groups, and there were many cases where these interactionsweredisrupted for extended periods by rotation without replacement. The presence of these intramolecular hydrogen bonds may have been responsible for the number of hydrogen bonds to water molecules being slightly less than 3 for many of the hydroxyl groups. If this type of hydrogen bond is not physically realistic, then the nonbonded potential function used here would appear to weight this type of intramolecular interaction too strongly. In severalreported carbohydrate crystal structures, ring oxygen atoms or glycosidic linkage oxygen atoms are engaged in intramolecular hydrogen bonds.lJ9 These hydrogen bonds may be the result of packing considerations which preclude the possibility of stronger hydrogen bond partners for the hydroxyl groups involved due to structural restrictions. In the present MD simulations, the ether-like ring and glycosidic oxygen atoms had a substantially lower partial charge (Table I) than the hydroxyl oxygen atoms, significantly reducing their ability to compete for hydrogen bond donors. As a result of this smaller charge, as well as steric restrictions, these groups were found to make few hydrogen bonds in these solution simulations. Figure 13 displays the pair distribution function for water oxygen atoms around the glycosidic linkage oxygen atom of the A-well solution simulation, which is typical of these functions for all of the linkage and ring oxygen atoms. The small peak indicates that there is some weak tendency to serve as a hydrogen bond acceptor. The first peak in this distribution is at 2.9 A, and the first minimum is at 3.5 A. By integrating the pair distribution function out to 3.5 A, 1.24 nearest-neighbors are obtained; this atom made an average of 0.37hydrogen bonds using the geometric definition (see Table
The Journal of Physical Chemistry, Vol. 97, No. 4, 1993 965
The Role of Hydrogen Bonding in Carbohydrates
I
120.0
120
80.0
I
I
I
1
I
I
I
t
I
1
80
40 -120.0
*
1
o 20.0
-80 ‘20’o 60.0
-120
q
-180 -80
-40
0
40
80
120
$ Figure 14. Solution trajectory (using alternate set of charges) starting from the C-well conformation superimposedon the adiabatic energy map. Contours are in 2 kcal/mol intervals from 2.0 to 10.0 kcal/mol above the energy minimum.
111). In general, the ring oxygen atoms made somewhat more hydrogen bonds than did the linkage oxygen atoms, perhaps due to greater solvent accessibility. Counting all types of hydrogen bonding groups, the total average number of water molecules hydrogen bonded to the maltose solute was 18.5 for the A-well simulation and 19.0 for the C-well simulation. Mean square displacements for the solute center of mass were computed, and diffusion coefficients for maltose were estimated from the limiting slopes of these curves. The diffusion coefficient was found to be 9.01 X 10” cm2 s-I for the C-well simulation. The experimentally measured value for this quantity is 4.825 X 10” cmz s-1 for 0.124 M solution.40 TIP3P water, which is less structured than real water, has been found to have a diffusion coefficient almost twice the experimentally observed value,41 so it is perhaps not surprising that in this less viscous solvent the calculated diffusion coefficient for maltose would also be higher than experiment by approximately the same factor. In order to test the sensitivity of the trajectory results to the selection of potential energy function, an additional simulation was performed using the alternate set of charges listed in Table I. This simulation began with the maltose solute in the C-well conformation. In contrast to the first such simulation, several conformation transitions were observed (Figure 14), including an extended excursion in the A well during which the internal 02-03’ hydrogen bond was broken and exchanged for hydrogen bonds to the solvent. Figure 15 illustrates the evolution of 4 and and Figure 16 displays the history of the 02-03’ distance in this simulation. The molecule spent the majority of the simulation in the B conformation, in which the 02-03’ hydrogen bond is geometrically possible; however, this bond exchanged for solvent hydrogen bonds for much of the time that the molecule was in this conformation (Figure 15) as well. Apparently, lowering the hydroxyl oxygen charge, without lowering the magnitude of the hydrogen atomic charge, decreases the strength of the intramolecular interaction relative to a mixed hydrogen bond to a TIP3P water molecule (which has a much larger charge of -0.834 for the oxygen) by a sufficiently large amount to favor disruption of this interaction. Lowering the hydroxyl oxygen charge also apparently resulted in lowering the barrier between the wells. Over the course of this second simulation, all three wells were sampled as the result of severalconformational transitions, while no transitions were observed in theother two solution simulations. Averaged over the entire trajectory, the maltose molecule made
+,
80.0
Time (pa)
-40
-120
60.0
40.0
t
1
0.0
-60.0 -120.0
20.0
40.0
80.0
80.0
Time (ps)
Figure 15. History of 4 and II.for the alternatecharge solution simulation starting from the C well over the entire trajectory, including equilibration.
,
t
20.0
40 0
60.0
80 0
Time (pa)
Figure 16. History of the distance between 0 2 and 03’ for the alternate charge solution simulation starting from the C well over the entire trajectory, including equilibration.
17.6 hydrogen bonds to water molecules in this simulation, which is fewer than was observed with the standard charge set.
IV. Conclusions These simulations reveal that a wide variety of hydrogen bonds are formed in carbohydrate solutions, both internally and with the solvent, and that these hydrogen bonds can potentially have significant influences upon the conformation of the solute as well as the structure of the solvent. Some of the effects of solvation appear to be the direct consequence of specific hydrogen bonds, rather than dielectric screening or mean field effects from distant water molecules. Such hydrogen-bonding effects would thus not be present in simulations employing mean field treatments of solvation4*or in vacuum calculations which modeled the absent water molecules with a large dielectric constant. However, there may be many carbohydrate systems in which solute-water hydrogen bonds play little role in the conformational structure, as is apparently true for the Maria( 1+2)Mana di~accharide.’~ Solvation was found to have a significant damping effect upon disaccharide conformationalfluctuations,as well as the magnitude of these fluctuations. In contrast to previous simulations, fluctuations in solution were observed to be larger in the present simulations than in vacuum, possibly due in part to water competition for a stabilizing intramolecular hydrogen bond. As expected, solute hydroxyl groups were found to tend to form
966
The Journal of Physical Chemistry, Vol. 97, No. 4, 1993
three hydrogen bonds to solvent, in the absence of competing internal hydrogen bonds. Care must be used in the selectionof solute-solvent force fields for carbohydrate simulations, and in particular the solute atomic partial charges, due to the role that these charges play in hydrogen bonding. The choice of water model used in molecular mechanics calculationsalso appears tobe important in determining the extent of hydrogen bond structuring around the solutes,due to previously known differences in the structural and dynamical behavior of the different water models. The balance of the relative strengths of water-water, water-solute, and solute-solute (intramolecular) hydrogen bonds can be quite important in determining physical properties, and presents a particular challenge in modeling simulations. Due to the similarlity of hydroxyl groups in sugars and in water, modeling this balance accurately may demand greater precision in molecular mechanics potential energy functions that many are capableof producing. It is quite possible that a simple static empirical function may not be adequate for reproducing this type of behavior, and in the future it may be necessary to develop dynamically polarizable potential functions such as those being developed for proteins and pure water.43In the absence of definitive experimental data concerning the solution conformations of carbohydrates, caution would seem to be warranted in the use of presently available potential functions which did not consider such hydrogen bond exchange data in their initial parameterization.
Acknowledgment. We thank K. Tasaki, A. D. French, A. D. MacKerell, and M. Karplus for helpful discussions. This work was supported in part by grants from NIH and the National Dairy Promotion and Research Board and Hatch Project 143410, USDA. References and Notes (1) Jeffrey, G. A.; Saenger, W. Hydrogen Bonding in Biological Structures; Springer-Verlag: Berlin, 1991. (2) Jeffrey, G. A. In Computer Modeling of Carbohydrate Molecules; French, A. D., Brady, J. W., Eds.; ACS Symposium Series 430; American Chemical Society: Washington, D.C., 1990. (3) French, A. D., Gardner, K. H., Eds. Fiber Diffraction Methods; ACS SymposiumSeries 141; AmericanChemicalSwiety: Washington, D.C., 1980. (4) Carver, J. P. Curr. Opin. Struct. Biol. 1991, I , 716. ( 5 ) Brooks, C. L.; Karplus, M.; Pettitt, B. M. Proteins: A Theoretical Perspective of Dynamics, Structure, and Thermodynamics; Advances in Chemical Physics, Vol. LXXI; Wiley-Interscience: New York, 1988. (6) French, A. D., Brady, J. W., Ed. ComputerModelingof Carbohydrate Molecules; ACS Symposium Series 4 3 0 American Chemical Society: Washington, D.C., 1990. (7) Brandt, D. A. Annu. Rev. Biophys. Bioeng. 1972,l. 369. Jordan, R. C.; Brant, D. A.; Cesaro, A. Biopolymers 1978, 17, 2617. (8) Lemieux, R. U.; Bock, K.; Delbaere, L. T. J.; Koto, S.; Rao, V. S. Can. J . Chem. 1980, 58, 631. (9) Koelher, J. E. H.; Saenger, W.; van Gunsteren, W. F. J . Mol. Biol. 1988, 203, 24 1.
Brady and Schmidt (10) Brady, J. W. J. Am. Chem. Soc. 1989,111, 5155. (1 1) Ha, S.; Gao, J.; Tidor, B.; Brady, J. W.; Karplus, M. J . Am. Chem. SOC.1991, 113, 1553. (12) Kroon-Batenburg, L. M. J.; Kroon, J. Biopolymers 1990, 29, 1243. (13) Van Eijck, B. P.; Kroon-Batenburg, L. M. J.; Kroon, J. J. Mol.Strucr. 1990, 237, 315. (14) Edge, C. J.; Singh, U. C.; Bauo, R.; Taylor, G. L.; Dwek, R. A.; Rademacher, T. W. Biochemistry 1990, 29, 1971. (15) Homans, S. W. Biochemistry 1990, 29, 9110. (16) Gress, M. E.; Jeffrey, G. A. Acta Crystallogr. 1977, 833, 2490. (17) Lipkind, G.M.; Verovsky, V. E.; Kochetkov, N. K. Carbohydr. Res. 1984,133,l. Shashkov, A. S.; Lipkind,G. M.; Kochetkov, N. K. Carbohydr. Res. 1986, 147. 175. (18) Stevens, E. S.;Sathyanarayana. B.K. J. Am. ChemSoc. 1989, I l l , 4149. (19) Blackwell, J.;Sarko, A,; Marchessault, R. H. J. Mol. Biol. 1969,42, 379. (20) Rees, D. A.; Smith, P. J. C. J. Chem. Soc. Perkin Tram. 2 1975,836. (21) Ha, S.N.; Madsen, L. J.; Brady, J. W. Biopolymers 1988,27,1927. (22) Tran, V.; Buleon, A.; Imberty, A.; Perez, S.Biopolymers 1989, 28, 679. (23) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.;Karplus, M. J. Comput. Chem. 1983, 4, 187. (24) Ha, S. N.; Giammona, A.; Field, M.; Brady, J. W. Curbohydr. Res. 1988, 180, 207. (25) Perez,S.;Marchessault,R. H. Carbohydr.Res. 1978,65,114. Jeffrey, G.A.;Taylor, R. J. Compur. Chem. 1980, I, 99. Tvaroska, I.; Bleha, T. Ado. Carbohydr. Chem. Biochem. 1989,47, 45. (26) Reiher, W. E. Ph.D. Thesis, Harvard University, Cambridge, MA, 1985. (27) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (28) Verlet, L. Phys. Reu. 1967, 159, 98. (29) van Gunsteren, W. F.; Berendsen, H. J. C. Mol. Phys. 1977, 34, 1311. (30) Hirata, F.; Rossky, P. J. J. Chem. Phys. 1981, 74, 6867. (31) Gao, J. J. Phys. Chem. 1992, 96, 6432. (32) Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1974, 60, 1545. (33) Berendsen, H. J. C.; Postma, J. P. M.;vanGunsteren, W. F.;Hermans, J. In Intramolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981; p 331. (34) Lau, W. F.; Pettitt, B. M. J. Med. Chem. 1989, 32, 2542. (35) Marchessault, R. H.; Perez, S.Biopolymers 1979, 18, 2369. (36) Nishida, Y.; Ohrai, H.; Meguro, H. Tetrahedron Lrrr. 1984, 25, 1575. Nishida, Y.; Hori, H.; Meguro, H.J . Carbohydr. Chem. 1988,7,239. (37) In a previous simulation of glucose in SPC water,Io the exocyclic group was found to adopt the TG conformation, while in a subsequent series of studies of the same solute sugar in TIP3P water," the exocyclic group was predominantly in the G T conformation, with a small fraction of time spent in the GG conformation, and with no TG conformations surviving the equilibration periods of the simulations. This difference was tentatively attributed in ref 11 to the difference in water models used, but subsequent simulations in our laboratory of this same system would seem to indicate that the differences between the long-range smoothing algorithms employed were actually the cause of the differences. (38) Christofides, J. C.; Davies, D. B. J. Chem. SOC.,Chem. Commun. 1985, 1533; J. Chem. Soc. Perkin Tram. 2 1987, 97. (39) Chu, S. S.C.; Jeffrey, G.A. Acta Crystallogr., Sect. B 1968,24,830. (40) Ueadaira, H.; Uedaira, H. Bull. Chem. Soc. Jpn. 1969, 42, 2140. (41) Tasaki, K.; McDonald, S.;Brady, J. W. J . Compur. Chem., in press. (42) Tvaroska, I.; Perez, S.Carbohydr. Res. 1986, 149, 389. (43) Sprik, M.; Klein, M. L.J. Chem. Phys. 1988,89,7556. Cieplak. P.; Kollman, P.; Lybrand, T. J. Chem. Phys. 1990,92,6755. Caldwell, J.; Dang, L.X.; Kollman, P. A. J . Am. Chem. Soc. 1990,112,9144. Straatsma, T. P.; McCammon, J. A. Mol. Simul. 1990, 5, 181.
