11339
J, Phys. Chem. 1995,99, 11339-11343
Use of Umbrella Sampling in the Calculation of the Potential of Mean Force for Maltose in Vacuum from Molecular Dynamics Simulations R. K. Schmidt, B. Teo, and J. W. Brady* Department of Food Science, Stocking Hall, Cornel1 University, Ithaca, New York 14853 Received: March 15, I995@
Molecular dynamics simulations have been used to calculate the potential of mean force for conformational transitions in Ramachandran space for the disaccharide /?-maltose in vacuum. The negative of the adiabatic energy map for maltose in vacuum was used as an umbrella potential in this calculation to increase sampling in high-energy regions. The calculated potential of mean force was not identical to the vacuum adiabatic surface. Entropic effects were found to lower the barrier between two of the wells on the vacuum surface, leading to facile transitions between these conformations, as has been observed in standard molecular dynamics simulations of this molecule, and resulting in an increased probability for that region of the Ramachandran map, in accord with experimental evidence.
I. Introduction The conformations of carbohydrate oligomers and polymers are, in general, poorly understood. This situation is unfortunate given the importance of carbohydrates in many biological processes. For example, glycoconjugateshave been identified as playing an important role in many molecular recognition events, with the oligosaccharide component of the molecule serving as the actual recognition site.' Given this role of carbohydrates in molecular recognition in particular, it would be useful to develop an understanding of the conformational determinants for such molecules. The analysis of oligosaccharide conformations has long been limited by a severe lack of high-quality experimental information, due to a variety of problems. Very few single-crystal diffraction studies of oligosaccharide structures are available because of the difficulty of obtaining crystals,* and until recently, reliable NMR data sufficient to specify structures has also been limited and difficult to ir~terpret.~ A number of disaccharide crystal structures have been determined, but it is not always certain that their conformation in aqueous solution is the same as that found in the c r y ~ t a l . ~Fiber diffraction data is available for some polysaccharides but the information obtained from such studies is also quite limited and must be supplemented with molecular mechanics calculation^.^ For these reasons, theoretical studies of oligosaccharides could potentially be useful. Generally, molecular mechanics analyses of oligosaccharide conformations have employed Ramachandran-like conformational energy maps which plot the molecular energy as a function of the glycosidic linkage torsion angles 41and ly (Figure 1) to identify low-energy conformation^.^.^ Initially such maps neglected monomer flexibility: but more recent studies have calculated relaxed or adiabatic energy maps in which the molecular energy is minimized for each (@,$J)value by allowing the other degrees of freedom of the molecule to adjust through energy minimi~ation.~,~ Figure 2 illustrates an adiabatic energy map for the prototypical disaccharide maltose, the basic repeat unit of starch, calculated using a previously reported CHARMMtype force field.8 Such adiabatic maps of the mechanical energy have been quite useful in the analysis of carbohydrate conformations and dynamics. However, molecular conformationsare actually determined by the free energy, which also includes such
* Abstract published in Advance ACS Abstracts, July 1, 1995. 0022-3654/95/2099-11339$09.00/0
Figure 1. @-Maltosein the energy-minimized A conformation. 120
PI
I'
1
'
' I" ' I '
"
I
'
" " "
I
" '
I
"
' I'
"
I , "~7-
I
"
'
1'4
80 40
y
o -40
-80
-120
-80 -40
0
40
80
120
Q Figure 2. Adiabatic conformational energy map for @-maltose in vacuum, calculated on a 20' x 20" grid, using the energy function of ref 16. Energy contours are indicated at 2,4,6, and 8 kcal/mol above the global minimum in the C well. The dotted line illustrates the reaction path used in the calculation of Figure 4. 0 indicates the conformation of the crystal structure. The shaded areas indicate grid regions in which various close contacts occur; within 5ach region the atoms of the indicated atomic pair are closer than 2.0 A in the lowest energy structure.
contributions as entropic effects and environmental influences. Thus it would be desirable to obtain conformational maps of the total free energy, or potentials of mean force (PMF), for important disaccharides in aqueous solution. 0 1995 American Chemical Society
Schmidt et al.
11340 J. Phys. Chem., Vol. 99, No. 29, 1995 In principle, the PMF, W($,W), can be calculated from the relative probabilities of occurrence for different configurations during molecular dynamics (MD) simulations or Monte Carlo (MC) since
P(4,W) = c exp[-PW(4,V)I
W(@,q)= -kT In P ( 4 , V )
(1)
+ C'
(2)
where
(3) In practice, this approach does not generally succeed since MD and MC calculations tend to avoid high-energy regions of conformation space, frustrating statistically meaningful estimates of the normalizing configuration integral. So-called "umbrella sampling" techniques have been developed to circumvent this d i f f i c ~ l t y . ' ~In . ' ~this approach, the normal potential energy of a system is supplemented with an umbrella potential VU which has the effect of increasing the probability of the system having high-energy conformations. The probability density for the system using this biasing function, P*($,Y), can be calculated by direct averaging from the MD simulation. The effects of this additional forcing or umbrella potential function can then be eliminated from the calculation of the relative conformational probability density for the unbiased potential function, P($,*), by employing a computational trick. Multiplying both the numerator and denominator of eq 3 by l / e - ~ v ~ ,
where ( )U indicates averaging over the ensemble using the umbrella potential VU. Substituting eq 5 into eq 2 gives W(4,V) = - k T h P * ( 4 , q ) - Vu(4,V)
+C
(6)
Umbrella sampling calculations have been used extensively in the study of one-dimensional system^.'^-'^ Often simple harmonic constraining potentials can be used to keep a system within a given range of some conformational transition coordinate E, and a series of such separate calculations are subsequently pieced together to span the range of 4 using appropriately chosen scaling factors in regions of overlap for the separately normalized cal~ulations.'~Far fewer studies of two-dimensional systems have been reported, in part due to difficulties in normalization. To avoid possible problems in normalization, one single VU which spans the entire range of E could alternatively be used. The negative of the PMF would be the ideal choice for this purpose since it would, on average, smooth out all variations in energy with conformation, so that each conformation was equally probable. However, since the PMF is unknown, and its determination is in fact the objective of the calculation, some other umbrella function which is a good estimate of the PMF will be needed. As a first approximation, the negative of the adiabatic energy map should be a reasonable choice, since this energy surface must contribute significantly to the total free energy. Deviations could result from errors in the adiabatic map (that is, deviations from the actual effective potential energy experienced at each ($,*) point) or from
entropic effects. Environmental influences, such as solvent interactions, might also be expected to affect the free energy surface. It would thus be interesting to use such a vacuum relaxed energy map to attempt to calculate PMFs in aqueous solution, in order to determine the effects of solvent. In practice, such a calculation is expected to be both difficult and expensive. Due to the difficulties expected in the calculation of the solution surface, it would be instructive to first calculate a PMF in vacuum both for comparison with the solvent surface and to explore the operational difficulties which might be encountered in the calculation of the solution case. The vacuum calculation can be extended to far longer simulation times than the solution case due to the much smaller number of atoms. More importantly, the vacuum case has far fewer important conformational states than does the molecule in aqueous solution. In vacuum, the only possible hydrogen bond partners for the numerous hydroxyl groups of the carbohydrate are either the adjacent groups of the same rings or, in certain cases, hydroxyl groups on the other ring of the dimer. Previous analysis of maltose found that in vacuum only a few hydroxyl conformational states of low energy exist, which generally involve all of the hydroxyl groups pointing in the same direction around the periphery of the rings in a network of hydrogen bonds between adjacent hydroxyl group^.^ In addition, the distribution of primary alcohol conformations is distorted since only a few allow hydrogen bonds either to the same ring or to hydroxyl groups on the opposite ring. In solution simulations these restrictions are relaxed because the solvent water molecules provide hydrogen bond partners and allow numerous local conformational transitions. In this paper we report the results of applying umbrella sampling to the calculation of the PMF for &maltose in vacuum from MD simulations, in preparation for the more difficult solution study. 11. Methods The relaxed potential energy surface for maltose using a preliminary CHARMM-type energy function for carbohydratesI6 has previously been reported.8 The central portion of this Ramachandran surface contains three low-energy conformations, labeled A, B, and C, separated by barriers of 1-2 kcal/mol (Figure 2). The global minimum-energy structure is in the C well, which also contains the crystal structure for ma1t0se.l~ The negative of this surface was used as the umbrella potential for a series of three estimates of the PMF in vacuum calculated using eq 6 from molecular dynamics simulations. Significant problems were encountered in preliminary simulations using this umbrella potential function. These problems arose from the fact that for regions away from the central portion of the map the energy rises steeply due to steric repulsions resulting from van der Waals overlaps of atoms on different rings. Figure 2 also schematically illustrates the principal steric clashes arising for various combinations of glycosidic angles. These respulsive forces rise so steeply that the umbrella potential cannot smooth out the variations effectively, and even augmented with the umbrella potential the trajectories do not sample these regions with adequate frequency. To overcome this problem the calculation was restricted to the region between -120.0" and +120.0', illustrated in Figure 2, by imposing an additional quadratic restraining force to prevent excursions into the biologically uninteresting high-energy regions outside of this range. The same steric clashes result in conformational transitions in the reducing ring for some regions of the map even within these bounds, from the 4Cl chair to various twist-boat forms, in order to relieve the strains. Our previously reported adiabatic map contains several regions in which the lowest energy
Potential of Mean Force for Maltose
J. Phys. Chem., Vol. 99, No. 29, 1995 11341
TABLE 1
120
trajectory percentages for primary alcohol conformations TG:GG:GT starting run conform
I I1
B/TG
111
BRG
ClGG
80
nonreducing
no. of transitions
reducing
no. of transitions
94.9:0.8:4.3 79.9:2.5:17.7 91.1:1.7:7.2
14 10 16
44.3:48.2:7.5 28.4:64.4:7.2 49.5:42.1:8.4
22 16 22
40
-40
structure corresponds to a twist-boat geometry for the reducing ring.8 These ring transitions involve energy barriers as a function of coordinates other than I$ and v, and as a result, once the ring "flips" to the twist-boat forms, it does not return to a chair form even when it traverses regions of glycosidic angle space where the steric clashes which occasioned the transition do not occur. For this reason, an additional constraining force was applied to maintain the rings in the 4C1 chair conformation. During such transitions the H2-H4, H3-H5, and 01-H5 distances in the nonreducing ring and the H2'H4', H3'-H5', and Hl'-H5' distances in the reducing ring exceed 3.6 A, as determined empirically from trajectories which exhibited such transitions. In the PMF calculation an additional constraining force as a function of these atomic separations was applied which was zero for interatomic distances less than 3.6 and which increased quadratically with a force constant of 1000 kcaVA2 for distances greater than 3.6 A.18 Both of these additional constraining forces were treated as additional terms in the umbrella potential. The adiabatic map illustrated in Figure 2 was generated using these constraints but differs only slightly from the previously reported map, and only in the boat regions. Three Langevin dynamics simulations of maltose in vacuum were calculated using the program CHARMM,I8 two starting from the B conformation and one starting from the C conformation. Each was integrated for 6 ns at a temperature of 300 K. The umbrella potential energy was calculated from a cubic spline of the adiabatic surface, which was calculated on a 20" x 20" grid. The probability density P*(I$,v)was calculated for 5" x 5" bins over the entire surface by counting the number of configurations in each bin over the duration of each trajectory. These three trajectories were then used to produce three separate estimates of the PMF, which were subsequently averaged together after adding an arbitrary constant so that all three surfaces matched at the B-well minimum energy geometry (-35", -15"). 111. Results and Discussion
The three simulations were started using minimized geometries corresponding to the B and C well structures. For the potential energy function employed in these calculations, the C well minimum has the GG primary alcohol conformation for the reducing ring (as does the A well), while the B well structure has the TG conformation.8 For all three wells the primary alcohol of the nonreducing ring is in the TG conformation. However, infrequent transitions in the reducing ring primary alcohol conformation did occur in the simulations, sampling all three possible conformers in all three simulations. Table 1 lists the proportions of each rotamer from each simulation and the number of transitions which occurred. The disproportionate representation of TG is in part due to the absence of solvent, since the TG conformation allows an intramolecular hydrogen bond for the nonreducing ring, which will be favored in vacuum due to the absence of competition from solvent. However, the large proportion of TG even for the reducing ring indicates that this explanation is not the sole factor involved. The potential energy function used here is known to favor the TG conforma-
-80 -120 -120 -80 -40
80
40
0
120
0 Figure 3. Potential of mean force for p-maltose as calculated from the present simulations. The dashed contour is at 1 kcal/mol above the lowest energy structure on the surface in the B well, and the solid line contours are at 2 kcaVmol intervals above this minimum energy. No contouring is presented for regions for which no sampling occurred during the trajectories, giving rise to the irregular outer boundary. 141.0
140.0 139.0
Ba
138.0
1
137.0 136.0 L
f
,
,
-80
:
.. l
,
,
-60
,
l
,
,
-40
,
l
8
-20
~
,
1
0
8
,
20
,
1
,
~
~
1
,
~
40
4 Figure 4. Potential of mean force calculated along the dotted conformational transition path illustrated in Figure 2. The dotted line represents the mechanical energy of the adiabatic map; the other three lines indicated by squares, circles, and triangles correspond to the three individual PMFs calculated from each of the three simulations.
tion over the GT and GG forms, in contrast to e~periment.'~ The results of the two simulations begun in the B well (I and In) were more similar to each other in their rotameric distributions than they were to the one begun from the C well geometry (II), indicating that even in these long calculations the rotameric distributions were not completely converged. The differences, however, are small. Hydroxyl rotations also occurred, but due to the much larger number of possible states, it is less likely that all were sampled appropriately. Presumably the resulting calculated PMF was not affected significantly by this reduced sampling, however, since most hydroxyl orientations which differed from the all-clockwise or all-counterclockwise patterns should have been of much higher energy. Figure 3 displays the average vacuum PMF calculated from these simulations. Contouring is only presented for those regions sampled by the trajectories, giving rise to the irregular outer boundary of the map. Each of the three individual estimates was very similar to this average surface. As can be seen, the PMF is qualitatively similar to the adiabatic surface in most regions, as would be expected. Interestingly, the only significant differences between these two surfaces are that for the PMF the barrier between the A and B wells has nearly disappeared and the C well is no longer the global minimum. The lowering of this barrier can be easily seen in Figure 4, in which the calculated PMF is plotted along the dotted reaction path illustrated in Figure 2 for all three individual calculations,
Schmidt et al.
11342 J. Phys. Chem., Vol. 99, No. 29, 1995 TABLE 2 local
energy
minimum
lP
u a
A B C
-61
-43 -38
-13 33
(kcal/moW 0.4 0.0 0.1
13
Approximate, estimated from nearby grid points. Relative to the energy of the B minimum.
80
40
y
o -40
-80 -120 -120 -80 -40
0
40
80
120
Q Figure 5. An example of a vacuum molecular dynamics trajectory projected onto the adiabatic map of the mechanical energy, illustrating transitions between the A and B wells, which occur frequently in ensembles of such trajectories, in contrast to transitions between the B and C wells.
along with the potential energy from the vacuum adiabatic map. The barrier between the B and C wells is essentially unaffected, while the A and B wells almost coalesce. The global minimum on this composite surface is in the B well, but the ordering of the three wells, which are all close in energy, was different for each of the three separate calculations. Table 2 lists the relative energies of these three local minima for the average PMF surface. In general, the vacuum adiabatic surface has been found to explain quite well the trajectories of maltose in vacuum calculated from molecular dynamics simulations. The relatively large barriers between minima result in few transitions between the wells in such simulations. The main exception to this statement is that more transitions between the A and B wells have been observed than might be expected,8 based on the barrier between these conformations on the adiabatic map (Figure 5, and also Figure 14 of ref 8). However, these transitions seem quite reasonable based on the negligible free energy barrier on the PMF surface seen in Figure 3. Presumably this attenuation results from entropic effects involving degrees of freedom of the molecule other than the glycosidic angle coordinates of the Ramachandran map. A greater entropy for the A-B region might be expected, since there are apparently more alternate states of low energy in this region than in the C region (see Figure 11 of ref 8).
IV. Conclusions As can be seen from the results presented here, the use of potentials of mean force to analyze conformational stability and dynamics is superior to simple energy surface analysis even for the vacuum case. The observed coalescence of the A and B wells for the maltose surface in vacuum, apparently due to entropic effects, explains the anomalous results of MD trajectories which were difficult to understand in terms of the Ramachandran map alone. The known limitations of the
potential energy surface employed here, as well as the artificial nature of the calculation, which excludes solvent effects, make it difficult to interpret the results directly in terms of the physical behavior of actual maltose. It is interesting, however, that the increased probability of the conformations to the lower left on the map, away from the crystal structure at (5", 13") in the vicinity of the C well, is in accord with recent chiroptical studies of maltose in solution,20which suggest a shift to this region of conformational space. Since it is quite possible that the presence of solvent would further affect the conformational behavior of carbohydrates,*' it would be desirable to prepare an analogous map for maltose in solution, either with the present potential function or with the revised and improved energy parameters currently being developed in this laboratory. However, the present studies indicate that the calculation of PMFs for disaccharides in solution will be difficult for several reasons. In addition to the dramatic increase in cost of the MD simulations due to the far larger number of atoms in the condensed phase studies, difficulties arise from the complexity of carbohydrate structure which would not necessarily arise in a molecule with a simpler structure, such as the model compound alanine dipeptide.22 These problems arise from the difficulty of representing the energy, which is a complex function of a number of intemal variables, as an explicit function of just two angles. The umbrella potential energy used depends explicitly only upon the glycosidic torsions 4 and tp, which are calculated from the positions of only 5 atoms, H1, C1, 01, C4', and H4'. The adiabatic potential is the lowest possible energy for a particular (4,W)point and does not necessarily represent the instantaneous energy arising from the clash of two hydroxyl groups, for example, which could be relieved at the same values of glycosidic torsion angles by the rotation of a hydroxyl group, the deformation of an angle, the compression of a bond length, or a combination of several deformations or ring torsions. One possible solution to this problem might be to add additional terms to the umbrella potential which are explicit functions of the interatomic separations involved in the clashes illustrated in Figure 2 and which serve to attenuate these interatomic repulsions somewhat. In addition to the problem of steric clashes, each hydroxyl group in the molecule has multiple conformational states, each of different energy, as do both primary alcohol groups, and the barriers between these states inhibit conformational sampling, further increasing the length of simulations required for the calculation of meaningful results. While the conformational energy may not depend significantly on all of these "hidden" degrees of freedom, for those that do significantly affect the conformational energy, many transitions will be required to achieve a Boltzmann distribution of occupancies. Until this is achieved, the calculated energies will be skewed by the kinetic trapping of the system in specific conformational states between the infrequent jumps, with the resulting calculated P * ( ~ , w ) reflecting these kinetic bottlenecks rather than the true conformational probabilities. Overcoming this difficulty will require longer simulations, or possibly different approaches altogether, such as has been used in the study of peptide PMFs.*~ Acknowledgment. The authors thank A. Kuki, K. Naidoo, and G. Liang for helpful discussions. This work was supported by NSF grant CHE-9307690. R.K.S. thanks the NSF for a graduate student fellowship. References and Notes (1) Lemieux, R. U. Chem. SOC.Rev. 1989, 18, 347-374. (2) Imberty, A,; Chanzy, H.;Ptrez, S.;Bulkon, A,; Tran, V. J . Mol. B i d . 1988, 201, 365-378.
Potential of Mean Force for Maltose (3) Engelsen, S. B.; PCrez, S.; Braccini, I.; du Penhoat, C. H. J . Comput. Chem., in press. (4) Pkrez, S.; Taravel, F.; Vergelati, C. Nouv. J. Chim. 1985, 9,561-
_S h_ 4...
(5) French, A. D., Gardner, K. H., Eds. Fiber Dijfraction Methods; ACS Symposium Series 141; American Chemical Society: Washington, DC, 1980. (6) Brant, D. A. Annu. Rev. Biophys. Bioeng. 1972, I , 369-408. (7) Brady, J. W. Adv. Biophys. Chem. 1990, I , 155-202. (8) Ha, S. N.; Madsen, L. J.; Brady, J. W. Biopolymers 1988,27, 19271952. (9) Tran, V.; BulCon, A,; Imberty, A.; PCrez, S. Biopolymers 1989, 28, 679-690. (10) Brooks, C. L.; Karplus, M.; Pettitt, B. M. Proteins: A Theoretical Perspective of Dynamics, Structure, and Thermodynamics; Adv. Chem. Phys., Vol. LXXI; Wiley-Interscience: New York, 1988. (1 1) McCammon, J. A.; Harvey, S. C. Dynamics of Proteins and Nucleic Acids; Cambridge University Press: Cambridge, U.K., 1987. (12) Valleau, J. P.; Tome, G. M. In Statistical Mechanics, Part A: Equilibrium Techniques; Beme, J., Ed.; Plenum: New York, 1977. (13) King, P. M. In Computer Simulation of Biomolecular Systems, Vol. 2 : Theoretical and Experimental Applications; van Gunsteren, W. F.,
J. Phys. Chem., Vol. 99, No. 29, 1995 11343 Weiner, P. K., Wilkinson, A. J., Eds.; ESCOM: Leiden, 1993. (14) Pangali, C.; Rao, M.; Beme, B. J. J . Chem. Phys. 1979, 71,2975298 1. (15) van Eijck, B. P.; Hooft, R. W. W.; Kroon, J. J . Phys. Chem. 1993, 97, 12093-12099. (16) Ha, S. N.; Giammona, A.; Field, M.; Brady, J. W. Carbohydr. Res. 1988, 180, 207-221. (17) Gress, M. E.; Jeffrey, G. A. Acta Crystallogr. 1977, 833, 24902495. (18) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J . Comput. Chem. 1983, 4, 187-217. (19) Nishida, Y.; Hori, H.; Ohrui, H.; Meguro, H. J . Carbohydr. Chem. 1988, 7, 239-250. (20) Stevens, E. S.; Sathyanarayana, B. K. J . Am. Chem. SOC.1989, 111, 4149-4154. (21) Brady, J. W.; Schmidt, R. K. J . Phys. Chem. 1993, 97,958-966. (22) Hemans, J. In Computer Simulations of Biomolecular Systems; van Gunsteren, W. F., Weiner, K., Eds.; ESCOM: Leiden, 1989. JF'9507354
