Molecular Dynamics Simulations of Carbohydrates ... - ACS Publications

Chapter 5

Molecular Dynamics Simulations of Carbohydrates and Their Solvation 1

2

Downloaded by MONASH UNIV on April 11, 2016 | http://pubs.acs.org Publication Date: July 6, 1990 | doi: 10.1021/bk-1990-0430.ch005

L. J. Madsen , S. N. Ha, V. H. Tran , and J. W. Brady Department of Food Science, Cornell University, Ithaca, NY 14853-7201

Molecular dynamics (MD) simulations are a class of molecular mechanics calculation which directly model the motions of molecular systems, often providing considerable information which cannot be obtained by any other technique, theoretical or experimental. MD simulations have only recently been applied to problems of carbohydrate conformation and motions, but it is likely that this technique will be widely used for modeling carbohydrates in the future. This paper introduces the basic techniques of MD simulations and illustrates the types of information which can be gained from such simulations by discussing the results of several simulations of sugars. The importance of solvation in carbohydrate systems will also be discussed, and procedures for including solvation in molecular dynamics simulations will be introduced and again illustrated from carbohydrate studies. Molecular mechanics calculations are an attempt to understand the physical properties of molecular systems based upon an assumed knowledge of the way i n which the energy of such systems varies as a function o f the coordinates of the component atoms. While t h i s term i s most c l o s e l y associated with the conformational energy analyses of small organic molecules pioneered by A l l i n g e r (1), i n t h e i r more general applications molecular mechanics calculations include energy minimization studies, normal mode c a l c u l a t i o n s , molecular dynamics (MD) and Monte Carlo simulations, reaction path analysis, and a number o f r e l a t e d techniques (2). Molecular mechanics 1

Current address: Department of Physical Sciences, Mesa College, Grand Junction, CO 81502

2

Current address: Laboratoire de Physicochimie des Macromolécules, Institut National de la Recherche Agronomique, B.P. 527, Nantes, Cédex F-44026, France 0O97-6156/90/O43O-O069$06.5O/0 © 1990 American Chemical Society

French and Brady; Computer Modeling of Carbohydrate Molecules ACS Symposium Series; American Chemical Society: Washington, DC, 1990.


70

COMPUTER MODELING OF CARBOHYDRATE MOLECULES

calculations have long been used to analyze biopolymer conformations, perhaps beginning with the work by Ramachandran on the allowed conformations of polypeptides (3) . The h i s t o r y of molecular mechanics conformational energy studies of carbohydrates i s almost as old, with the f i r s t such studies of carbohydrates (4) coming very shortly a f t e r the i n i t i a l studies of peptides. Conformational energy calculations have become an i n t e g r a l part of polysaccharide studies (5), p a r t i c u l a r l y i n the i n t e r p r e t a t i o n of f i b e r d i f f r a c t i o n data (6,7). As the t h e o r e t i c a l study of protein and nucleic acid conformations matured, the inherently dynamical behavior of these molecules was recognized, and beginning with the simulations of Karplus and coworkers, these motions were studied t h e o r e t i c a l l y by d i r e c t l y modeling them i n molecular dynamics simulations (8). The systematic investigation of protein motions has revealed that i n many cases the physical properties and b i o l o g i c a l function of these polymers cannot be understood without taking i n t e r n a l motions into consideration (2). Although MD simulations of these classes of molecules have now become quite commonplace, p a r a l l e l studies of carbohydrates were not attempted u n t i l quite recently. Only i n the l a s t few years have MD studies of carbohydrates been reported (913). While the importance of f l e x i b i l i t y i n carbohydrate molecules has long been recognized (14-16), these simulations have served to d i r e c t l y i l l u s t r a t e the e f f e c t s of such dynamical f l e x i b i l i t y , and i t i s c e r t a i n that MD simulations of carbohydrates w i l l become routine i n the near future. I t i s also now generally understood that the structure and b i o l o g i c a l function of many biomolecules are affected by aqueous solvation. For t h i s reason, t h e o r e t i c a l models of biopolymers must include these solvent e f f e c t s as well as i n t e r n a l f l e x i b i l i t y . The unique s t r u c t u r a l behavior of water, p a r t i c u l a r l y immediately adjacent to solutes, makes i t d i f f i c u l t to apply continuum theories to aqueous solutions. However, with the development of high speed computers, i t i s now becoming f e a s i b l e to d i r e c t l y model the behavior of aqueous systems through MD and Monte Carlo simulations which s p e c i f i c a l l y include solvent water molecules. A wide v a r i e t y of such calculations have now been reported, including studies of pure water (12) and solutions of nonpolar atoms (18-20), ions (21,22), alcohols (23,24), urea (25,26), peptides (27,28), and even proteins (29). MD simulations of carbohydrate solvation could be p a r t i c u l a r l y useful, since there are a number of unanswered questions concerning the i n t e r a c t i o n of these molecules with water (30,31). Complex solution behavior might well be expected f o r carbohydrates, which contain a number of adjacent polar hydroxyl groups held i n r e l a t i v e l y f i x e d positions around the pyranoid rings, as well as nonpolar CH and CH groups and the e t h e r - l i k e r i n g oxygen atoms. This expectation of complex s o l u t i o n behavior has been f u l f i l l e d i n the solution simulations of carbohydrate molecules which have been reported (13,32, and Ha, S.; Gao, J . ; Tidor, B.; Brady, J.W.; Karplus, M. submitted to J_g. Am Chem Soc.), and demonstrates the need f o r further studies which include the e f f e c t s of aqueous solvation. In t h i s paper, the general technique of MD simulations w i l l be reviewed and examples of applications to 2


5. MADSEN ET AL

71

Molecular Dynamics Simulations

carbohydrate molecules w i l l be discussed, as well as the i n c l u s i o n of aqueous solvation into such studies.


Molecular Dynamics Simulations Molecular dynamics calculations are a technique f o r modeling physical systems on the microscopic scale i n which atomic motions are d i r e c t l y simulated by numerically solving the c l a s s i c a l Newton's equations of motion f o r a l l of the atoms i n the system subject to the forces a r i s i n g from some given force f i e l d . MD simulations were o r i g i n a l l y developed to model simple physical systems such as monatomic rare gases (33), but they have now been applied to a wide v a r i e t y of problems, including water and aqueous solutions, and e s s e n t i a l l y a l l types of biopolymers (2). Several sophisticated MD programs are currently generally available, including the well-known molecular mechanics packages CHARMM (34), AMBER (35), and GROMOS (36), and many such programs are available commercially. A fundamental requirement of molecular mechanics studies of any type, including MD simulations, i s a complete description of the v a r i a t i o n of the t o t a l p o t e n t i a l energy of the system of i n t e r e s t as a function of the molecular coordinates. S t r i c t l y speaking, t h i s energy i s given by the complete solution of the Schrôdinger equation for the entire system, including solvent molecules, as t h e i r positions evolve with time. Since f o r macromolecules and condensed phases the accurate c a l c u l a t i o n of t h i s quantum mechanical energy i s not possible, i t i s common to employ analytic, semi-empirical energy expressions which have t h e o r e t i c a l l y reasonable functional forms and which have been parameterized to the results of experiments and quantum mechanical calculations of simple molecules. Extensive sets of such potentials have been developed to describe condensed phases of small molecules, including water (1Z»32,1S)» simple organic compounds (1,39), and various biopolymers (34-36). Such empirical energy surfaces have also been developed f o r sugars (40-42), since carbohydrates contain functional groups which do not occur i n other b i o l o g i c a l molecules and since p o t e n t i a l parameters are sensitive to chemical environment. Unfortunately, to date none of these carbohydrate parameterizations have proven to be completely s a t i s f a c t o r y . Work on developing improved force f i e l d parameters for various types of molecules i s continuing. Semiempirical p o t e n t i a l energy surfaces vary i n d e t a i l and even i n form from one set to the next (2), but most represent the intramolecular p o t e n t i a l as a sum of e l e c t r o s t a t i c and van der Waals interactions between non-bonded atoms and terms f o r hindered r o t a t i o n about molecular bonds. Bond stretching and angle bending forces are derived from harmonic restoring potentials, and s p e c i a l functions are sometimes included to account f o r hydrogen bonding. A t y p i c a l example of such a function of the i n t e r n a l coordinates q might be 2

V(q) - Z k ^ i b i - b i o ) + Z k ^ i - ^ o ) + Σ(Α^ΑΪ

2

2

+ Sk [l + cosOi^-^)] 0 i

- Bij/rfj + q ^ / r ^ )


(1)

72


where k and \α are the bond stretching and angle bending force constants, φ i s a t o r s i o n angle with a force constant k^ , p e r i o d i c i t y n, and phase factor 6 and q and qj are the atomic p a r t i a l charges. A and are van der Waals constants and r ^ i n t h i s equation represents the interatomic distance between atoms i and j . An energy function of t h i s type allows the c a l c u l a t i o n of the v a r i a t i o n of the system energy as a function of the positions of every atom i n the system, and because i t has an a n a l y t i c form, also permits the d i r e c t analytic evaluation of the derivatives of t h i s function with atomic positions, which gives the atomic forces. Of course, i n order f o r molecular mechanics calculations to be of use, i t i s necessary f o r the p o t e n t i a l energy functions to be s u f f i c i e n t l y r e a l i s t i c as to adequately mimic physical behavior. For t h i s reason considerable e f f o r t goes into the development of the parameters which appear i n the energy function (2,34)· This parameterization i s generally accomplished by the matching of calculated properties to experimental measurements, as a function of the parameter set f o r selected small model compounds. Representing the molecular p o t e n t i a l energy as an a n a l y t i c function of the nuclear coordinates i n t h i s fashion i m p l i c i t l y invokes the Born-Oppenheimer approximation i n separating the very f a s t e l e c t r o n i c motions from the much slower ones of the n u c l e i . This separation allows the electronic energy at each nuclear configuration to be represented as the p o t e n t i a l energy f o r the motions of those n u c l e i f o r that configuration (2). The energy function parameterization i s thus an empirical approximation to t h i s quantum mechanical energy, and i n t y p i c a l c l a s s i c a l molecular mechanics calculations, i s the only place where the quantum behavior of the molecular system i s included. The system i s subsequently considered to be governed s o l e l y by c l a s s i c a l mechanics, whether the simulation i s an energy minimization, normal mode analysis, or MD or Monte Carlo c a l c u l a t i o n . For most s t r u c t u r a l and dynamical properties of i n t e r e s t t h i s i s a v a l i d approximation, but biopolymers, including sugar, do contain hydrogen atoms, whose mass i s s u f f i c i e n t l y small that quantum e f f e c t s f o r these atoms can be s i g n i f i c a n t . Quantum e f f e c t s may be p a r t i c u l a r l y important i n chemical reactions, electron and proton transport, and other events which involve large changes i n e l e c t r o n i c d i s t r i b u t i o n s . From elementary c l a s s i c a l mechanics (43), the force acting on an atom i n a molecular mechanics system i s the negative gradient, or derivative, of the p o t e n t i a l energy function with respect to atomic position, θ±

bi

±

t

it

±


i(j

Fi - - ViV

(2)

These forces can be substituted into Newton's equations of motion,

Fi - m ^

(3)

which can i n p r i n c i p l e be integrated numerically on a large d i g i t a l computer to provide a complete description of the motion of every atom i n the molecular system as a function of time, q(qo,v ,t), f o r 0


5.

MADSENETAL.

73


any given set of i n i t i a l positions q and v e l o c i t i e s v at the i n i t i a l time t . Because the forces are conservative (that i s , derivable from a p o t e n t i a l as i n equation (2)), the t o t a l energy of the system i s constant throughout the simulation, as are the l i n e a r and angular momenta. Properties calculated from such a simulation are thus microcanonical ensemble properties (constant energy, number of p a r t i c l e s , and volume) (44). The t o t a l energy, or Hamiltonian, i s made up of p o t e n t i a l and k i n e t i c energy terms, 0

0

0

E

T

- Εκ + V

(4)


which are not constant and which fluctuate as p o t e n t i a l energy i s converted into motion ( k i n e t i c energy) and then back again. The instantaneous temperature of the MD system i s c a l c u l a t e d from the v e l o c i t i e s v* of the atoms as

(5)

where Ν i s the t o t a l number of atoms i n the system, and t h i s quantity w i l l of course fluctuate as the k i n e t i c energy fluctuates. The thermodynamic temperature of the system i s the mean value of t h i s instantaneous temperature averaged over a s u f f i c i e n t l y long dynamics sequence such that the temperature converges to a stable value c h a r a c t e r i s t i c of a l l possible motions i n the system under the given conditions. In a t y p i c a l MD simulation, a s t a r t i n g structure q i s selected, usually from crystallographic data, or i n the case of s u f f i c i e n t l y simple small carbohydrates, stereotypical chair or boat forms, and i n i t i a l v e l o c i t i e s v f o r every atom are selected from a Boltzmann d i s t r i b u t i o n at the desired system temperature (2) . The integration of the equations of motion i s c a r r i e d out numerically using one of several integration procedures. The most commonly used method i s the V e r l e t algorithm (2,45,46), which calculates the p o s i t i o n at a time t'-t+At through a truncated Taylor series expansion i n the time step At about the p o s i t i o n at t, x(t) . More elaborate integration schemes, such as the Gear predictor-corrector algorithm are also sometimes used (2), but these require considerably more computer time. The integration time step size At i n the V e r l e t algorithm must be small (usually i n the femtosecond range) r e l a t i v e to the f a s t e s t motions of the system f o r the procedure to be r e l i a b l e , which means that the simulation of a lengthy period of time w i l l require quite a large number of time steps. Each time step of the integration requires the c o s t l y evaluation of the forces at that step, thus making the c a l c u l a t i o n expensive i n terms of computer time and l i m i t i n g the periods which may be p r a c t i c a l l y simulated. Since the s e l e c t i o n of s t a r t i n g conditions i n MD simulations can be somewhat a r b i t r a r y , and not necessarily r e a l i s t i c , i t i s necessary to " e q u i l i b r a t e " the system by integrating the equations of motion f o r some period of time during which the behavior may not 0

i 0



74


be t r u l y representative of the physical system. During t h i s period the integration i s c a r e f u l l y monitored to ensure that the energy i s being conserved and that the temperature i s stable and has the desired value. I f as the r e s u l t of relaxation processes the system temperature d r i f t s away from the desired value by more than some small tolerance, i t i s customary to either scale a l l of the atomic v e l o c i t i e s i n the system by an appropriate factor to bring the temperature back to the s p e c i f i e d value, or to again assign new values f o r a l l of the atomic v e l o c i t i e s by random s e l e c t i o n from a Boltzmann d i s t r i b u t i o n . Because i t may not be p h y s i c a l l y representative, the e q u i l i b r a t i o n portion of the trajectory, which must l a s t f o r a number of picoseconds, i s not analyzed f o r the c a l c u l a t i o n of system properties. A f t e r the system has s t a b i l i z e d s u f f i c i e n t l y to be considered "equilibrated" (that i s , a l l a r t i f i c i a l stresses produced by the s e l e c t i o n of i n i t i a l conditions have relaxed away), the integration i s continued f o r a much longer period without further intervention for the purpose of a c t u a l l y simulating equilibrium dynamical behavior. In p r i n c i p l e , when a simulation has been integrated f o r a s u f f i c i e n t l y long period of time, the mean properties computed as time averages over the entire MD simulation w i l l converge to the t r u e thermodynamic, canonical ensemble properties of the system. Unfortunately, i t i s d i f f i c u l t to know when t h i s "thermodynamic l i m i t " has been reached, and i n systems of even moderate complexity the thermodynamic l i m i t may not be attainable with r e a l i s t i c computer times. Reported MD simulations of carbohydrates have t y p i c a l l y ranged from 10-20 ps to 500 ps i n length, although supercomputers now permit simulations of modest-sized systems i n the nanosecond range. M

M

Molecular Dynamics Simulations of Carbohydrates The f i r s t reported molecular dynamics simulations of carbohydrates began to appear i n 1986, with the p u b l i c a t i o n of studies of the vacuum motions of a-D-glucopyranose (9), discussed below, and the dynamics of a hexa-NAG substrate bound to lysozyme (10), which are described i n greater d e t a i l i n the chapter by Post, et a l . i n t h i s volume. Since that time, simulations of the dynamics of many more carbohydrate molecules have been undertaken. A number of these studies are described i n subsequent chapters of t h i s volume. The introduction of t h i s well developed technique to problems of carbohydrate structure and function could contribute s u b s t a n t i a l l y to the understanding of t h i s class of molecules, as has been the case f o r proteins and related biopolymers. One of the most s i g n i f i c a n t p o t e n t i a l contributions of MD studies to the f i e l d of carbohydrate chemistry i s i n exploring the degree of f l e x i b i l i t y of these molecules and i n i l l u m i n a t i n g any possible b i o l o g i c a l or s t r u c t u r a l roles f o r such f l e x i b i l i t y . Although polysaccharides have been u s e f u l l y modeled i n the past as r i g i d monomer units (5), sugars are not completely r i g i d , and when MD simulations are applied to carbohydrates a wide v a r i e t y of i n t e r n a l motions, fluctuations and conformational t r a n s i t i o n s are found to occur. For example, i n MD simulations of glucopyranose (9), when MD t r a j e c t o r i e s were i n i t i a t e d i n the C conformation, A

X



5.

MADSENETAL.

75


the molecule remained i n t h i s stable chair conformation, but the r i n g a c t u a l l y executed numerous small-scale o s c i l l a t i o n s about t h i s average conformation. Figure 1 displays the h i s t o r y through a p a r t i c u l a r MD simulation of one of the r i n g torsion angles, C3-C4C5-C6. As can be seen, t h i s angle was c l e a r l y not constant i n time during the simulation, although i t s mean value was stable as the r i n g o s c i l l a t e d about the average value. Fluctuations i n r i n g torsion angles such as seen i n Figure 1 r e s u l t i n fluctuations i n the o v e r a l l shape, or pucker, of sugar rings. I t i s possible to represent the puckering away from p l a n a r i t y of a r i n g of atoms by a set of parameters developed by Cremer and Pople, which describe the degree of puckering and the conformational form of the r e s u l t i n g puckered structure (47). Normal thermal motions of the type represented i n Figure 1 r e s u l t i n a continual f l u c t u a t i o n i n t h i s puckering. Larger fluctuations i n r i n g pucker can lead to actual t r a n s i t i o n s i n r i n g conformation, as from chair to boat forms. Figure 2 i l l u s t r a t e s the evolution of one of the r i n g t o r s i o n angles, C1-C2-C3-C4, f o r a glucose trajectory which began with the molecule i n the C conformation (a discussion of t h i s conformational terminology can be found i n reference 48), but which underwent a t r a n s i t i o n to a twist boat form half-way through the simulation. Figure 3 displays the h i s t o r y of the Cremer-Pople pucker parameter Θ, which s p e c i f i e s conformational form, f o r the same trajectory. Values of θ around 0° correspond to the C conformation, and values around 90° correspond to the various twist-boat forms. As can be seen from the figures, the t r a n s i t i o n which occurred i n t h i s trajectory involved several changes i n the r i n g torsion angle, but because these o s c i l l a t i o n s were accompanied by correlated changes i n adjacent torsion angles, they produced only one change i n o v e r a l l r i n g conformation, which occurred d i r e c t l y and without returning to the higher energy 0 form. X

4

h

1

χ

4

MD simulations have also demonstrated that the pendant alcohol groups on carbohydrate rings are not r i g i d , but undergo frequent rotations and o r i e n t a t i o n a l t r a n s i t i o n s . For example, i n the same series of simulations of α-D-glucopyranose discussed above (9), the primary alcohol group was not found to be locked into the crystallographic GT conformation, but instead frequently jumped between the three major low energy conformers available to i t , and established an equilibrium d i s t r i b u t i o n between these three forms based upon t h e i r r e l a t i v e Boltzmann-weighted p r o b a b i l i t i e s . In addition to rotations about the C5-C6 bond, the hydroxyl groups themselves also frequently rotate, changing t h e i r orientation. Figures 4 and 5 displays the h i s t o r y of the t o r s i o n a l angle C4-C5C6-06, specifying the conformation of the primary alcohol group, and for the angle 05-C1-01-H, specifying the orientation of the anomeric hydroxyl group f o r a t y p i c a l α-D-glucopyranose t r a j e c t o r y i n vacuum calculated from an MD simulation using the Rasmussen PEF422 energy function (41). Such motions cannot be adequately studied by conformational energy studies alone, and require dynamics simulations f o r t h e i r characterization. For D-glucose, NMR studies have demonstrated that as a r e s u l t of s t e r i c crowding, the C form i s not present at room temperature X

4


COMPUTER MODELING OF CARBOHYDRATE MOLECULES 180.0

120.0 h 60.0


-60.0 Ρ

-120.0 I-

-180.0 Time (ps)

Figure 1. History of the dihedral angle C3-C4-C5-05 calculated from a t y p i c a l molecular dynamics simulation of a α-Dglucopyranose i n the C conformation i n vacuum. (Reproduced from Ref. 9. Copyright 1986 American Chemical Society.) ii

1

180.0

120.0 h

60.0

-60.0

-120.0

-180.0

Figure 2. History of the r i n g torsion angle C1-C2-C3-C4 calculated from a molecular dynamics simulation of the motions of an α-D-glucopyranose molecule i n vacuum which began i n the *0 conformation and which subsequently underwent a t r a n s i t i o n to a twist-boat conformation. (Reproduced from Ref. 9. Copyright 1986 American Chemical Society.)


4

MADSEN ET A L


77

150.0 h


120.0

Time (ps)

Figure 3. History of the Cremer-Pople pucker parameter θ calculated f o r the same t r a j e c t o r y i l l u s t r a t e d i n Figure 2. (Reproduced from Ref. 9. Copyright 1986 American Chemical Society.)

y>i ^ y · ν Γ V * Η

180.0

1

1

'' 1

1

— F T — Π Ι ψ^ί Ή | | Ί 1

120.0

60.0 h

0.0

-60.0

-120.0

-180.0 TIME (PICOSECONDS)

Figure 4. History of the torsion angle C4-C5-C6-06 c a l c u l a t e d from a t y p i c a l MD simulation of α-D-glucopyranose i n vacuum using the Rasmussen p o t e n t i a l energy function PEF422 (41). (Reproduced from Ref. 9. Copyright 1986 American Chemical Society.)


ν

'" Ι | !


78

to any great extent (49). At higher temperatures, however, t h i s form must occasionally a r i s e spontaneously, as the formation of levoglucosan (1,6-anhydro-0-D-glucopyranose) demonstrates (50). This molecule r e s u l t s from the 1,6 elimination of water from β-Ώglucopyranose, which can only happen when the inversion of the r i n g from the usual C conformation to the higher energy C conformer brings these two hydroxyl groups into close proximity. Such a spontaneous r i n g inversion has been observed i n molecular dynamics simulations of 0-D-glucopyranose (11) using the PEF 422 energy function (41). Figure 6 i l l u s t r a t e s the h i s t o r y of the Cremer-Pople pucker parameter θ from a simulation of 0-D-glucopyranose at a somewhat elevated temperature (320K) which was i n i t i a t e d i n the lowest energy C conformation. A f t e r approximately 21 ps of simulation, the molecule underwent a spontaneous t r a n s i t i o n to higher energy twist-boat forms, and a f t e r approximately 60 ps, b r i e f l y converted to the s t i l l less favorable *C conformation, where the elimination reaction i s possible, before returning to a twist-boat and f i n a l l y to the o r i g i n a l C conformation. Transitions such as t h i s i n monomer r i n g geometries, which can a r i s e spontaneously i n MD simulations but which cannot be e a s i l y incorporated into s t a t i c conformational energy c a l c u l a t i o n s , could have important consequences i f they occurred i n polysaccharides. Such chair-to-boat t r a n s i t i o n s have been suggested to play a r o l e i n determining polysaccharide structure (end-to-end lengths) i n aqueous s o l u t i o n (51). 4

X

X

A

4


X

4

4

X

Solvation Since the cost i n computer time of MD simulations increases r a p i d l y with the number of p a r t i c l e s (atoms) i n the system, many simulations of proteins and n u c l e i c acids have studied these biopolymers i n vacuum, due to the f a r larger cost of including solvent molecules. Many important b i o l o g i c a l processes, however, are the d i r e c t or i n d i r e c t r e s u l t of aqueous solvation, and simulating p h y s i c a l properties of t h i s type requires the i n c l u s i o n of solvent molecules (2). Although some generalized s o l v a t i o n e f f e c t s such as d i e l e c t r i c screening can be represented by continuum models, many of the most i n t e r e s t i n g consequences of aqueous s o l v a t i o n require an e x p l i c i t representation of the solvent s t r u c t u r i n g immediately adjacent to the solute. For t h i s reason, i t i s c l e a r l y desirable to include water molecules i n MD simulations i n such a way as to be p h y s i c a l l y reasonable. In p r i n c i p l e , i t i s a simple matter to include solvent water molecules d i r e c t l y i n MD simulations, since appropriate intermolecular p o t e n t i a l energy functions f o r water are a v a i l a b l e ( l Z » ! Z » M ) î one would j u s t surround the solute molecules with a s u f f i c i e n t number of water molecules to approximate a bulk s o l u t i o n . Unfortunately, a " s u f f i c i e n t number" of water molecules might be enormous, since many of the e f f e c t s of aqueous s o l v a t i o n are long range or are due to entropie contributions a r i s i n g from "structuring" of the solvent, which may be cooperative i n nature. In any "droplet" representation of a solution, the water molecules farthest from the solvent would constitute a boundary between a bulk l i q u i d phase and a vacuum. Since such a boundary would


5.

MADSEN ET A L

79


180.0

120.0 h

60.0


-60.0 h

-120.0 h

-180 Time (ps)

Figure 5. History of the t o r s i o n angle 05-C1-01-H c a l c u l a t e d from a t y p i c a l MD simulation of α-D-glucopyranose i n vacuum using the Rasmussen p o t e n t i a l energy function PEF422 (41). (Reproduced from Ref. 9. Copyright 1986 American Chemical Society.)

150.0 h 120.0

90.0 h60.0 h 30.0

60.0

45.0

30.0

75.0

T i m e (ps)

Figure 6. History of the Cremer-Pople pucker parameter θ calculated from molecular dynamics simulation of 0-Dglucopyranose i n the *C conformation i n vacuum using the Rasmussen p o t e n t i a l energy function PEF422. This molecule underwent four t r a n s i t i o n s i n conformation. (Reproduced with permission from Ref. 11. Copyright 1987 E l s e v i e r Science Publishers.) X



80


s u b s t a n t i a l l y a f f e c t the properties of those solvent molecules i n the boundary region, a large number of molecules would need to be included i n the c a l c u l a t i o n to ensure that t h i s perturbed boundary region was f a r from the solute-solvent region of i n t e r e s t . In order to avoid edge e f f e c t s , an approximation known as "periodic boundary conditions" (2,52) i s often used, i n which the solute i s placed i n the center of a cubic box of f i n i t e dimensions, surrounded by water molecules, and t h i s entire box i s then surrounded i n every d i r e c t i o n with exact images of i t s e l f . Atoms near to the surface of t h i s primary box interact with image, or "ghost" water molecules i n the image boxes rather than with a vacuum, eliminating d i r e c t edge e f f e c t s . The equations of motion of the image p a r t i c l e s are not integrated d i r e c t l y , but can be generated at each force evaluation by applying the appropriate symmetry operations on the coordinates of the atoms i n the c e n t r a l , primary box. Usually "minimum-image" periodic boundary conditions are employed i n such calculations, which means that only the closest i j i n t e r a c t i o n f o r two atoms i and j and a l l t h e i r possible images i s included i n the i n t e r a c t i o n forces; on the assumption that the forces decrease strongly with distance. This c r i t e r i o n means neglecting interactions between atoms greater than one h a l f of the box length apart, which implies that no molecule, including the solute, interacts with i t s own image. The solute molecule must be covered with several layers of water molecules, such that the outer layer experiences l i t t l e perturbation from b u l k - l i k e behavior due to the solute. Under these minimum-image periodic boundary conditions, the r e s u l t i n g system thus represents a p h y s i c a l l y unattainable state, an i n f i n i t e l y d i l u t e solution at f i n i t e concentration. Such a system i s a microcanonical ensemble (constant volume, constant energy, constant N). Once the boundary conditions have been implemented, the c a l c u l a t i o n of solution molecular dynamics proceeds i n e s s e n t i a l l y the same manner as do vacuum calculations, while the t o t a l energy and volume i n a microcanonical ensemble c a l c u l a t i o n remain constant, the temperature and pressure need not remain fixed. A variant of the periodic boundary condition c a l c u l a t i o n method keeps the system pressure constant by adjusting the box length of the primary box at each step by the amount necessary to keep the pressure calculated from the system second v i r i a l at a f i x e d value (46) . Such a procedure may be necessary i n simulations of processes which involve large volume changes or fluctuations. Techniques are also available, by coupling the system to a Brownian heat bath, f o r performing simulations d i r e c t l y i n the canonical, or constant T,N, and V, ensemble (2,46). To date, only a few solution calculations f o r carbohydrates have been attempted (one such study of mannitol and s o r b i t o l i s described i n the chapter by Grigera i n t h i s volume), but the r e s u l t s of these early studies bear out the expectation that solvation e f f e c t s i n carbohydrate systems can be both s i g n i f i c a n t and d i f f i c u l t to predict. In the case of pyranoid rings, molecular solvation i s further complicated by the close j u x t a p o s i t i o n of these groups i n e s s e n t i a l l y f i x e d r e l a t i v e orientations (assuming no conformational changes i n the r i n g ) . Under such circumstances, molecular stereochemistry could play important physical roles, as i s



5.

MADSENETAL.


81

i l l u s t r a t e d by the simple hexose monosaccharides. A l l of these molecules are s t r u c t u r a l isomers d i f f e r i n g only i n t h e i r stereochemistries at the various r i n g carbon atoms, yet they exhibit d i f f e r e n t physical properties i n aqueous solution, such as the c h a r a c t e r i s t i c equilibrium anomeric r a t i o s (31,50), i n d i c a t i n g that these properties are determined at l e a s t i n part by aqueous solvation. As an example of the possible importance of the s p e c i a l nature of aqueous solvation i n determining the properties of carbohydrate solutions, i t i s worthwhile to consider the r e s u l t s of recent MD simulations of α-D-glucopyranose i n aqueous solution (32)· In t h i s study, a minimum-image periodic boundary conditions MD simulation was conducted f o r a single α-D-glucopyranose molecule surrounded by 207 SPC water molecules (38) at a temperature of 300K. This number of water molecules i n the primary box i s s u f f i c i e n t f o r two complete solvation s h e l l s around the solute, and a portion of a t h i r d s h e l l . The simulation was run f o r 10.2 ps of e q u i l i b r a t i o n and an a d d i t i o n a l 32 ps of data c o l l e c t i o n , using an integration step size of 1 f s . Although t h i s run was apparently not long enough to allow complete thermodynamic convergence of a l l of the physical properties examined, i t was quite lengthy i n terms of computer time, requiring more than 40 days of DEC VAX 11/750 CPU time, demonstrating one of the d i f f i c u l t i e s of solution calculations. The presence of aqueous solvent was found to have l i t t l e e f f e c t upon the mean structure of the pyranoid r i n g i n these MD simulations, with only s l i g h t deviations i n the time-averaged structure away from that observed i n vacuum simulations or i n the crystallographic d i f f r a c t i o n experiments (32). However, the presence of the solute had substantial e f f e c t s upon the average "structuring" of the solvent. Figure 7 displays a p a i r d i s t r i b u t i o n function g(r), defined as (27.) 1

dN(r)

g(r)

(7) 4πρΓ

2

dr

which i s the normalized p r o b a b i l i t y of finding the oxygen atom of a solvent water molecule a given distance r from some p a r t i c u l a r atom i n the sugar solute molecule, where ρ i s the bulk number density. The g(r) i l l u s t r a t e d i n Figure 7 shows the p r o b a b i l i t y of f i n d i n g a water oxygen atom as a function of distance from the glucose C6 carbon atom as calculated from the simulations. The form of t h i s curve i s t y p i c a l of hydrophobic or nonpolar hydration, with a broad peak centered around 3.4 A, the approximate van der Waals contact distance, representing the f i r s t solvation s h e l l water molecules. This type of nonpolar d i s t r i b u t i o n function contrasts with that found i n the case of hydrogen bonding, as i l l u s t r a t e d i n Figure 8, which displays the p a i r d i s t r i b u t i o n function f o r water oxygen atoms around the 06 hydroxyl oxygen atom. The f i r s t peak i n t h i s curve occurs at a much closer distance, 2.7 Â, which i s t y p i c a l of hydrogen bonding, since the nearest neighbors of t h i s hydroxyl group w i l l be those water molecules which are hydrogen bonded to i t . This much higher and quite narrow f i r s t peak, with a deep f i r s t minimum, indicates that these hydrogen bonded water molecules are very




g(r)

2.0

4.0

6.0

r (A) Figure 7. Water oxygen-exocyclic methylene carbon p a i r d i s t r i b u t i o n function, calculated from a molecular dynamics simulation of α-D-glucopyranose i n aqueous solution, giving the normalized p r o b a b i l i t y of finding a water oxygen atom a given distance r from the C6 carbon atom. (Reproduced from Ref. 32. Copyright 1989 American Chemical Society.)



MADSEN ET AL.


2.0

g(r)

r (A) Figure 8. Water oxygen-hydroxyl oxygen 06 p a i r d i s t r i b u t i o n function, calculated from a molecular dynamics simulation of α D-glucopyranose i n aqueous solution. (Reproduced from Ref. 32 Copyright 1989 American Chemical Society.)



84


l o c a l i z e d s p a t i a l l y by the requirements of the hydrogen bond, which balances strong e l e c t r o s t a t i c attractions with van der Waals repulsions. Figure 9 i l l u s t r a t e s the p a i r d i s t r i b u t i o n function f o r the water molecules around the r i n g oxygen; the breadth of the f i r s t peak i n t h i s function, centered around 3.25 Â, and i t s low maximum height, indicate that t h i s atom i s not s i g n i f i c a n t l y hydrogen bonded to i t s nearest water molecule neighbors, due to i t s lower charge and lack of a hydrogen atom substituent. In addition to imposing s p a t i a l r e s t r i c t i o n s on the d i s t r i b u t i o n s of adjacent water molecules, the various groups i n the solute sugar molecule also impose o r i e n t a t i o n a l structuring upon these solvent molecules. Figure 10 displays the d i s t r i b u t i o n s of orientations f o r water molecules around the methylene carbon C6. This function i s the i n t e g r a l l y normalized p r o b a b i l i t y f o r each water molecule being oriented such that i t makes an angle θ between i t s OH bond vectors and the vector from the water oxygen to the carbon atom. This function i s calculated f o r those molecules within 4.9 Â of the carbon atom (nearest neighbors), as t h i s distance marks the f i r s t minimum i n the p a i r d i s t r i b u t i o n function f o r that atom. The curve i n Figure 10 i s t y p i c a l f o r hydrophobic hydration (22) . The peak at c o s ( 0 ) - l corresponds to one hydroxyl group pointed d i r e c t l y away from the nonpolar group, and the value of nearly zero at -1.0 indicates v i r t u a l l y no p r o b a b i l i t y of a hydroxyl group pointing d i r e c t l y at the methylene group. The broad peak around 0.33, the tetrahedral angle, i s a consequence of the tetrahedral structure of the SPC water molecule; i f one of the hydroxyl groups i s pointing d i r e c t l y away from the CH group, then the other must be making a tetrahedral angle with the CO vector. This type of o r i e n t a t i o n a l structuring, somewhat s i m i l a r to a f l u c t u a t i n g clatherate, i s adopted as the best way to solvate such a small nonpolar species, as each water molecule can s t i l l make hydrogen bonds to other water molecules and thus avoid the high energetic cost of l o s i n g a hydrogen bond (22)· Figure 11 displays the same type of o r i e n t a t i o n a l d i s t r i b u t i o n function f o r water molecules adjacent to the 03 hydroxyl oxygen atom. As can be seen, i n the case of normal hydrogen bonding, the inverse behavior i s seen, with one of the water hydroxyl groups pointing d i r e c t l y at the solute hydroxyl oxygen atom. The poorly hydrogen-bonding r i n g oxygen atom imposes l i t t l e o r i e n t a t i o n a l structuring, as can been seen from Figure 12. Surprisingly, however, t h i s o r i e n t a t i o n a l d i s t r i b u t i o n function f o r those water molecules around the 02 hydroxyl oxygen atom, i l l u s t r a t e d i n Figure 13, was s u b s t a n t i a l l y perturbed i n these simulations. This perturbation was apparently due to the interference i n o r i e n t a t i o n a l structuring caused by other adjacent groups, perhaps i n p a r t i c u l a r the requirements of the two CH groups 1 and 2 which constitute an extended region of hydrophobicity, which i s closer to the 02 hydroxyl group than to the a x i a l 01 group i n the alpha anomer (see Figure 14). Preliminary simulations using other water models indicate that t h i s o r i e n t a t i o n a l perturbation i s indeed r e l a t e d to the anomeric preference, and i s apparently involved i n the solventinduced anomeric d i s t r i b u t i o n s observed i n aqueous s o l u t i o n (H,50). The solvation requirements of the s p e c i f i c s p a t i a l arrangements r e s u l t i n g from p a r t i c u l a r stereochemistries thus give r i s e to the 2



MADSENETAL.


g(r)

2.0

4.0

6.0

r (A) Figure 9. Water oxygen-ring oxygen p a i r d i s t r i b u t i o n function, calculated from a molecular dynamics simulation of a-Dglucopyranose i n aqueous solution. (Reproduced from Ref. 32. Copyright 1989 American Chemical Society.)


86


1.00

h

0.75

h

0.50

h

0.25

h


P(Cos(0))

Cos(0) Figure 10. D i s t r i b u t i o n of orientations f o r water molecules adjacent to the exocyclic methylene carbon atom C6 as calculated from a molecular dynamics simulation of α-D-glucopyranose i n aqueous solution. The function plots the frequency of occurrence of an angle θ between the water OH bond vectors and the vector from the carbon atom to the water oxygen atom. A value of cos(0) of 1.0 corresponds to an OH bond vector pointing d i r e c t l y away from the carbon atom. (Reproduced from Ref. 32. Copyright 1989 American Chemical Society.)

1.00

H

0.75

h

P(Cos(0)) 0.50

0.25

h

Cos(0) Figure 11. D i s t r i b u t i o n of orientations f o r water molecules adjacent to the hydrogen-bonding 03 hydroxyl oxygen atom, as calculated from an MD simulation of α-D-glucopyranose, as i n Figure 10. (Reproduced from Ref. 32. Copyright 1989 American Chemical Society.)


MADSEN ET AL.


1.00

h

0.75

h

0.50

h

P(Cos(0))


0.25

-1.0

-0.5

0.0

0.5

1.0

Cos(0) Figure 12. D i s t r i b u t i o n of orientations f o r water molecules adjacent to the r i n g oxygen atom of α-D-glucopyranose i n aqueous solution, calculated as i n Figure 10. (Reproduced from Ref. 32. Copyright 1989 American Chemical Society.)

1.00

0.75 P(Cos(0)) 0.50

0.25

h

-1.0

-0.5

0.0

0.5

1.0

Cos(0) Figure 13. D i s t r i b u t i o n of orientations f o r water molecules adjacent to the 02 hydroxyl oxygen atom o f α-D-glucopyranose i n aqueous solution, calculated as i n Figure 10. (Reproduced from Ref. 32. Copyright 1989 American Chemical Society.)



88


Figure 14. A stereoview of a t y p i c a l "snapshot" from a molecular dynamics simulation of α-D-glucopyranose i n aqueous solution, showing selected water molecules close to the 01 and 02 hydroxyl groups. The hydrophobic hydration requirements of the a l i p h a t i c hydrogen atoms on CI and C2 impose additional s t r u c t u r i n g r e s t r i c t i o n s of the orientations of the water molecules hydrogen bonding to the 02 hydroxyl group. (Reproduced from Ref. 32. Copyright 1989 American Chemical Society.)

differences i n solution properties of the various sugars. For t h i s reason, i t may often be necessary to include solvent i n t h e o r e t i c a l treatments of carbohydrate molecules. Conclusions Although they have only recently been applied to carbohydrate problems, molecular dynamics simulations should contribute s u b s t a n t i a l l y to our understanding of a v a r i e t y of s t r u c t u r a l and conformational phenomena. MD calculations can help to quantify rates f o r a v a r i e t y of t r a n s i t i o n processes, which other s t a t i s t i c a l mechanical simulation techniques (Monte Carlo) cannot explore. MD simulations may also help resolve " v i r t u a l " or time-averaged structures i n disaccharides which o s c i l l a t e between one or more lowenergy conformations (described i n the chapter by Carver, et a l . i n t h i s volume). More ambitious MD simulations may help answer o l d questions about the conformation of polysaccharides such as amylose and any possible role which random changes i n r i n g conformation may play. Because of the importance of aqueous solvation i n many b i o l o g i c a l systems, including the carbohydrates, t h e o r e t i c a l studies which include solvent i n an attempt to understand the exact r o l e of solvent i n these systems w i l l be necessary. Again, MD simulations o f f e r an excellent method f o r including solvation d i r e c t l y and examining s o l u t i o n behavior at a l e v e l of d e t a i l which i s not possible i n experiment. Presumably, MD simulations of carbohydrates w i l l soon become as commonplace as are conformational energy studies of these molecules.


5.

MADSENETAL.


89

Acknowledgments The authors thank A. French, M. Karplus, S. Perez, and J . Carver f o r h e l p f u l discussions. This work was supported i n part by NIH grant GM 34970 and USDA Hatch project 143-433.

Literature Cited 1) 2)


3) 4) 5) 6) 7) 8) 9) 10) 11) 12) 13) 14) 15) 16) 17) 18) 19) 20) 21) 22) 23) 24) 25) 26) 27) 28) 29)

Burkert, U.; Allinger, N.L. Molecular Mechanics. ACS Monograph 177, American Chemical Society, Washington, 1982. Brooks, C . L . ; Karplus, M.; Pettitt, B.M. Proteins: A Theoretical Perspective of Dynamics, Structure, and Thermodynamics. Advances in Chemical Physics. Wiley Interscience: New York, 1988, Vol. LXXI. Ramachandran, G.N.; Ramakrishnan, C.; Sasisekharan, V. J. Mol. Biol. 1963, 7, 95. Rao, V.S.R.; Sundararajan, P.R.; Ramakrishnan, C.; Ramachandran, G.N. in Conformation in Biopolymers. Vol. 2, G.N. Ramachandran, ed., Academic Press, London, 1967. Brant, D.A. Ann. Rev. Biophys. Bioeng. 1972, 1, 369. Wu, H.-C.H.; Sarko, A. Carbohydr. Res. 1978, 61, 7. Imberty, Α.; Perez, S. Biopolymers 1988, 27, 1205. McCammon, J . A . ; Gelin, B.R.; Karplus, M. Nature 1977, 267, 585. Brady, J.W. J. Am. Chem. Soc. 1986, 108, 8153. Post, C.B.; Brooks, B.R.; Karplus, M.; Dobson, C.M.; Artymiuk, P . J . ; Cheetham, J . C . ; Phillips, D.C. J. Mol. Biol. 1986, 190, 455. Brady, J.W. Carbohyd. Res. 1987, 165, 306. Prabhakaran, M.; Harvey, S . C . Biopolymers 1987, 26, 1087. Grigera, J.R. J . Chem. Soc., Faraday Trans. I, 1988, 84, 2603. Goebel, C.V.; Dimpfl, W.L.; Brant, D.A. Macromolecules 1970, 3, 644. Joshi, N.V.; Rao, V.S.R. Biopolymers 1979, 18, 2993. French, A.D.; Murphy, V.G. Carbohydr. Res. 1973, 27, 391; Polymer 1977, 18, 489. Stillinger, F . H . ; Rahman, A. J. Chem. Phys. 1974, 60, 1545. Dashevsky, V . G . ; Sarkisov, G.N. Mol. Phys. 1974, 27, 1272. Owicki, J . C . ; Scheraga, H.A. J. Am. Chem. Soc. 1977, 99, 7413. Swaminathan, S.; Harrison, S.W,; Beveridge, D.L. J. Am. Chem. Soc. 1978, 100, 5705. Palinkas, G.; Riede, W.O.; Heinzinger, Κ. Z. Naturforsch. 1977, 32a, 1137. Dang, L . X . ; Pettitt, B.M. J. Chem. Phys. 1987, 86, 6560. Okazaki, S.; Nakanishi, K.; Touhara, H. J. Chem. Phys. 1984, 81, 890. Tanaka, H . ; Nakanishi, K.; Touhara, H. J. Chem. Phys. 1984, 81, 4065. Kuharski, R.A.; Rossky, P.J. J . Am. Chem. Soc. 1984, 101, 5794. Tanaka, H . ; Nakanishi, K.; Touhara, H. J. Chem. Phys. 1985, 82, 5184. Rossky, P . J . ; Karplus, M. J. Am Chem. Soc. 1979, 101, 1913. Ravishanker, C.; Mezei, M.; Beveridge, D.L. J. Comput. Chem. 1986, 7, 345. Ahlstrom, P.; Teleman, O.; Jönsson, B. J. Am. Chem. Soc. 1988, 110, 4198.



90


30) Suggett, A. in Water: A Comprehensive Treatise. Vol. 4, F. Franks, ed., Plenum, New York, 1975, pp. 519-567. 31) Franks, F. Pure and Appl. Chem. 1987, 59, 1189. 32) Brady, J.W. J . Am Chem. Soc. 1989, 111, 5155. 33) Rahman, A. Phys. Rev. 1964, 136, A405. 34) Brooks, B.R.; Bruccoleri, R.E.; Olafson, B.D.; States, D.J.; Swaminathan, S.; Karplus, M. J . Comput. Chem. 1983, 4, 187. 35) Weiner, S.J.; Kollman, P.A.; Nguyen, D.T.; Case, D.A. J . Comp. Chem. 1986, 7, 230. 36) van Gunsteren, W.F.; Berendsen, H.J.C.; Hermans, J.; Hol, W.G.J.; Postma, J.P.M., Proc. Natl. Acad. USA 1983, 80, 4315. 37) Jorgensen, W.L.; Chandrasekhar, J.; Madura, J.D.; Impey, R.W.; Klein, M.L. J . Chem. Phys. 1983, 79, 926. 38) Berendsen, H.J.C.; Postma, J.P.M.; van Gunsteren, W.F.; Hermans, J . in Intermolecular Forces. B. Pullman, ed., Reidel, Dordrecht, 1981, 331. 39) Jorgensen, W.L., J. Am. Chem. Soc., 1981, 103, 335; Ibid. 1981, 103, 341. 40) Ha, S.N.; Giammona, Α.; Field, M.; Brady, J.W. Carbohydr. Res., 1988, 180, 207. 41) Rasmussen, K. Acta Chem. Scand., Ser. A 1982, 36, 323, and references therein. 42) Koeler, J.E.H.; Saenger, W.; van Gunsteren, W.F. Eur. Biophs. J. 1987, 15, 197. 43) Goldstein, H. Classical Mechanics. 2nd Edition, Addison-Wesley, Reading, MA, 1980. 44) McQuarrie, D.A. Statistical Mechanics, Harper and Row, New York, 1976. 45) Verlet, L. Phys. Rev. 1967, 159, 98. 46) McCammon, J.A.; Harvey, S.C. Dynamics of Proteins and Nucleic Acids, Cambridge University Press, Cambridge, 1987. 47) Cremer, D.; Pople, J.A. J . Am Chem. Soc. 1975, 92, 1354. 48) Stoddart, J.F. Stereochemistry of Carbohydrates, WileyInterscience, New York, 1971. 49) Perkins, S.J.; Johnson, L.N.; Phillips, D.C.; Dwek, R.A. Carbohydr. Res. 1977, 59, 19. 50) Shalleriberger, R.S. Advanced Sugar Chemistry, AVI Publishing, Westport, CT, 1982. 51) Goebel, K.D.; Harvie, C.E.; Brant, D.A. Appl. Polym. Symp. 1976, 28, 671. 52) Metropolis, N.; Rosenbluth, A.W.; Rosenbluth, M.N.; Teller, A.H.; Teller, E. J . Chem. Phys. 1953, 21, 1087. RECEIVED March 21, 1990


Molecular Dynamics Simulations of Carbohydrates ... - ACS Publications

Recommend Documents