Analysis of interfacial water structure and dynamics in α-maltose solution by molecular dynamics simulation

29 March 1996

CHEMICAL PHYSICS LETTERS

~ ,~.,:l i,II ~,

ELSEVIER

Chemical PhysicsLetters 251 (1996) 268-274

Analysis of interfacial water structure and dynamics in a-maltose solution by molecular dynamics simulation Curl Xin W a n g a, Wei Zu Chen a, V. Tran h, R. Douillard c a Center for Fundamental Physics, University of Science and Technology of China, Hefei 230026, People's Republic o f China b Laboratoire de Physico-Chimie des Macromol~cules, lnstitut National de la Recherche Agronomique (INRA), BP 1627, 44316 Nantes Cedex 03, France c Laboratoire de Biochimie et de Technologie des Prot~ines, lnstitut National de la Recherche Agronomique (INRA), BP 1627, 44316 Nantes Cedex 03, France

Received 7 July 1995;in f'malform 12 January 1996

Abstract The structure and dynamic behavior of the interfacial water around an a-maltose molecule are studied by molecular dynamics (MD) simulations. It is found that 8 hydrogen bonds (H-bonds) are formed between water and hydroxyl groups of each ring in maltose. The results of the radial distribution function of water around maltose are consistent with the hydrogen bonding analysis. The calculated diffusion coefficient of water around the different atom types in maltose shows that no general trend can be found, which may suggest that the solvent mobility is not influenced significantly by the atom types in the carbohydrate molecule.

1. Introduction Because of the significant role of carbohydrates in biological systems, an understanding of the aqueous solvation behavior of carbohydrate polymers is of considerable importance [l]. Nowadays computer modeling techniques, such as MD, are widely used to explore the detailed structural and dynamic behavior of carbohydrates in aqueous solution [2-7]. Hardy and Sarko [6] used the Amber force field [8] and TIP3P water model [9] to study the conformational behavior of cellobiose in water. Brady and Schmidt [7] used the CHARMM force field [10] and TIP3P water model to examine the role of hydrogen bonding in a B-maltose solution. In the present paper our particular interest lies in

the determination of the dynamic behavior of the hydration water and its role on interfacial properties from a structural and dynamic point of view. Although our modeling molecule, o~-maltose, is somewhat similar to the one studied by Brady and Schmidt [7], the points of interest and the details studied in this work are quite different. Only quite recently it has been found that the diffusion coefficient of hydration water molecules near a polar atom of a protein [11] and polypeptide [12] is the smallest one. To explore the influence of the carbohydrate molecules on the structure and mobility of the solvent, we have selected the small disaccharide, etmaltose, as our research object and used the GROMOS force field [13] and S P C / E water model [14] to study the solvent structure and dynamic behavior

0009-2614/96/$12.00 © 1996 Elsevier Science B.V. All rights reserved Pll S0009-2614(96)00110-8

C.X. Wang et al./ Chemical Physics Letters 251 (1996) 268-274

of the hydration shell around an a-maltose molecule by MD simulations.

2. Method and procedures The initial coordinates of maltose were provided by the C conformation of the crystal structure [7,15]. Fig. 1 shows the schematic representation of the 1-4 linkage a-maltose molecule, in which all non-polar hydrogen atoms were included in the carbon atoms (the united atom model). The polar hydrogen atoms in the hydroxyl groups of maltose were treated explicitly. The simulation was performed using the GROMOS force field [13] with force field parameters for sugar atoms reported in Re/', [16]. The maltose was centered in a rectangular box (30.4784 × 32.5743 × 36.6010 ,~3). This box was filled using an equilibrated cubic box of a bulk sample of S P C / E water at 300 K as a building block. A distance of 13 ,~ between the maltose surface and the planes of the box edges and a minimum distance of 2.3 A between the non-hydrogen atoms and the closest water oxygen atoms were maintained. The resulting box with a volume of 36337.9 ,~3 contained 1202 water molecules and 3637 atoms in total to give an initial full hydration system with a density of 1.004 g / c m s. The concentration of the resulting solution was 0.046 M. In the beginning, two types of energy minimizations were successively performed. A first 100 steps minimization with the steepest descent method was done keeping the maltose harmonically constrained to its original conformation, followed by another one (80 steps of steepest descent) without these conH03'

straints. Maxwellian velocities corresponding to 300 K were then assigned to each atom in the system to start the MD simulation. The temperature of the maltose and solvent were separately coupled to an external bath of 300 K by using a temperature relaxation time of 0.01 ps for the first 500 steps and 0.1 ps for the remaining simulation [17]. The MD time step was set to 0.002 ps. A non-bonded cutoff radius of 8 ,~ was used for the van der Waals interactions and 11 ,~, for the long-range Coulombic interactions. The SHAKE constraint algorithm was used to fix all covalent bond lengths and the water angle with a relative tolerance of 0.0001 [18]. The neighbor pair list was updated every 10 steps. Periodic boundary conditions were imposed on the box to remove edge effects. The MD simulations consisted of 30 ps equilibration followed by 100 ps for data collection. Configurations of the trajectories were saved every 20 fs (10 MD time steps) for further analysis. Most of the simulations and the associated analysis were performed on a SUN Sparc station 10 and a Vax station 3200.

3. Results and discussions First, we have focused on the analysis of the hydrogen bond (H-bond) between maltose and water. The criteria used to define an H-bond are purely geometric [19,20]. According to the method developed by Berendsen et al. [21], the average number of H-bonds per hydroxyl group to water, (Nhb), is determined by ( NhD) -= NhJNstep,

'

06

Fig. 1. Schematic plot of ~t-maltose.

~HOI'

( l)

where Nhb is the total number of H-bonds that are found in the averaging period (50 ps) and Nstep is the number of trajectory frames which is equal to 2500 in the present work. The evaluation of the H-bond life time, r, is given by

~'= ( ( Nhb) /gw )tMD, H06

269

(2)

where N w is the number of different water molecules that form the H-bonds and tMD is the ~ simulation time for the averaging. Table 1 reports the results for the H-bonds between maltose and the water molecules, which are obtained by averaging over the 30-80 ps simulation

C.X. Wang et al. / Chemical Physics Letters 251 (1996) 268-274

270

period. In Table 1, for each hydroxyl group, the sum of (Nhb) is between 1.77 and 2.29, and the sum of the corresponding H-bond life time is more than 6 ps. It is clear that the O6-HO6 group has the strongest H-bond with water ( ( N h b ) = 2.29). However, all non-hydroxyl oxygen atoms O4', 0 5 and 05' have the lower average number of H-bonds and the shorter H-bond life time with water (Table 1). All hydroxyl groups participate in forming H-bonds with water, in which eight H-bonds can be generated in each ring. This result is consistent with that obtained by Hardy and Sarko for the simulation of cellobiose in water [6].

Next, we have analyzed the radial distribution functions which are calculated around the different exposed solute atoms presented in this work, such as polar oxygen atoms in hydroxyl groups and non-polar carbon atoms. The radial distribution function of the solvent, g(r), is defined as

(AN(r)> g(r)

(3)

47rprAr '

where r is the distance between the solute atom and the solvent, A N ( r ) is the number of solvent molecules in the region between r - A r / 2 and r +

4.0

4.0 --.Ccg-Ow ..... O 2 - H w

3.0

--C2-Ow ..... C2-Hw

3.0

~ 1.0

2.0 1.0

0.0 0

2

4

6

~t

0.0

10

0

2

r (A) 4.0

,i 6 r (A)

lO

4.O

--03..Ow . . . . .

3.0

~

a

O3-Hw

--C"3-1~v ..... C 3 - H w

3.O

~

2.0

2.0

-,\ 1.0

1.0

O.O 0

2

4

6

8

10

0.0

0

2

4

r (A)

6

8

10

r (A)

4.O

4.0 --O4-Ow ..... O 4 - H w

3.0

I

--C4-Ow ..... C 4 - H w

3.0

2.0

1.0

,L

/\.

0.O O

2

4

6

s

lO

1.0~ o . o o~

. 2 ,/

4

r (A) 4.0

'~" 2.0

1.0

1.0

0

2

4

6 r (~,)

F i g . 2 . Radial distribution

10

--c6,ow ..... c 6 . . H w

2.0

O.0

8

4.0 --O6-O~ ..... O ~ H w

3.0

~

6

r (A)

function

g(r)

8

io

0.0

.:/'

2

4

6 r

8

1o

(A)

o f w a t e r a r o u n d the h y d r o x y l g r o u p s o f o n e r i n g i n m a l t o s e . T h e solid c u r v e s are plotted for the

distribution o f w a t e r o x y g e n around solute a t o m s and the d a s h e d c u r v e s for the distribution o f w a t e r h y d r o g e n .

C.X. Wang et a l . / Chemical Physics Letters 251 (1996) 268-274 Table 1 Average number of hydrogen bonds to SPC/E water molecules, (Nhb), and the hydrogen bond life time, v, for maltose in water solution over the simulation period 30-80 ps Atom name

Nhb

(N,b)

N.

'r (ps)

HOI' O1' HO2' 02' HO3' 03' 04' 05' HO6' O6'

2483 2529 2453 2566 2429 2198 689 1363 2450 3056

0.99 1.01 0.98 1.03 0.97 0.88 0.28 0.55 0.98 1.22

9 20 10 26 11 17 9 12 10 15

5.5 2.5 4.9 2.0 4.4 2.6 1.6 2.3 4.9 4.1

HO2 02 HO3 03 HO4 04 05 HO6 06

2458 2904 2377 2040 2426 2454 1262 2441 3268

0.98 1.16 0.95 0.82 0.97 0.98 0.50 0.98 1.31

9 16 10 22 10 24 12 14 23

5.4 3.6 4.8 1.9 4.9 2.0 2.1 3.5 2.8

Ar/2

from the solute atom, Ar is a small increment and p is the density of bulk water. Since the radial distribution functions for both water oxygen (O w) and water hydrogen (H w) can 4.0

271

provide a picture of the hydrogen bonding patterns, we have analyzed their distributions around different surface atoms of maltose. The shapes of the g(r) curves for both rings of maltose are very similar so that here we have only reported the distributions for one ring (Figs. 2 and 3). As seen in Fig. 2, the first peaks in the g(r) curves are sharper around the polar oxygen atoms in comparison with the non-polar carbon atoms. However, the intensities of the first peaks of the g(r) curves around the 04' and 0 5 atoms (Fig. 3) are smaller than that around the hydroxyl oxygen atoms (Fig. 2). This indicates that the average number of water molecules near a hydroxyl group exceeds the number near the non-hydroxyl atoms (O4', 0 5 and O5'). The first peak positions in g(r) o f O w and H w are listed in the 2rid and 3rd column of Tables 2 and 3, respectively. The intensities of the peaks are reported in the 4th and 5th column of Tables 2 and 3, respectively. We have found that the first sharp peaks in og(r) of O w around polar atoms locate at 2.7-3.0 A (Fig. 2, Table 2). It is clear that the sharper the first peak in g(r) of O , appears, the more H-bonds to water exist. This result is consistent with the hydrogen bonding analysis mentioned above. For instance, atom 0 6 has the sharpest first peak (Fig. 2). This correlates with the highest number of H-bonds to water (Table 1). The distributions of O w 4.0. --Cl-Ow ..... C I - H w

--O4'-Ow ..... O 4 ' - H w

3.O

3.0 ¸

2.0

I.O

1.0

. . ? . . ~ ~ 0.0

0

' 0.0

2

4r ( A ) 6

s

~ (A)6

s

to

4.0

4.0

--O5-Ow .....O5-Hw

3.oi

--C5-Ow .....C S - H w

3.0

2.0

,~ 2.0, 1.0. 0.0

2

lo

I

1.0-

2

4

6

r (/~)

8

10

0.0 4

.. . " ' ~ 2

4r

(A)6

s

to

Fig. 3. Same as Fig. 2 except for different solute atoms of C1, C5, 04' and 05 in maltose.

C.X. Wang et a l . / Chemical Physics Letters 251 (1996) 268-274

272

Table 2 Solvent characteristics around polar atoms of maltose

Table 3 Solvent characteristics around non-polar atoms of maltose a

Atom

Atom

First peak position in g(r) (/~) a

Intensity of the peak a

Ow

Hw

Ow

Hw

O1' 02'

2.7 2.7

2.44 2.76

03' 04' 05' 06' 02

2.7 3.0 3.0 2.7 2.7

03 04 05 06

2.7 2.7 3.0 2.7

3.3 1.8 3.2 3.3 3.2 1.8 3.2 3.2 3.2 1.8 3.1

1.44 1.24 1.42 1.34 1.45 1.35 1.27 1.56 1.46 1.25 1.53

2.57 1.18 1.53 2.88 2.72 2.48 2.58 1.41 3.02

Diffusion constant (,~2/ps) b

Intensity of the peak

g(r) (~.)

0.21 (0.03) 0.20 (0.04) 0.18 0.17 0.18 0.17 0.19

First peak position in

CI' C2' C3' C4' C5'

(0.04) (0.03) (0.03) (0.04) (0.02)

C6' CI C2 C3 C4 C5 C6

0.20 (0.02) 0.18 (0.02) 0.17 (0.02) 0.20 (0.03)

Diffusion constant

(~2/ps)

Ow

Hw

Ow

Hw

3.7 3.7 3.5 4.0 3.9 3.5 3.9 3.6 3.6 3.5 3.8 3.7

4.3 4.1 4.2 5.0 4.7 4.0 4.3 4.3 4.3 4.6 4.2

1.96 1.94 1.63 1.21 1.68 2.23 1.72 2.09 2.08 2.83 1.60 2.50

1.22 1.21 1.18 1.06 1.13 1.15 1.14 1.29 1.28 1.27 1.25

0.19 (0.02) 0.18 (0.04) 0.16 (0.03) 0.15 (0.02) 0.15 (0.03) 0.16 (0.03) 0.16 (0.02) 0.17 (0.02) 0.17 (0.02) 0.18 (0.02) 0.17 (0.03) 0.18 (0.03)

a See footnotes of Table 2.

a The mark " - " means an absence of a peak. b The values in parentheses are the statistical errors of the diffusion coefficient constant.

around most non-polar atoms have broad peaks near 3.7 ,~ (Figs. 2 and 3), which demonstrates the presence of a hydrophobic shell. Such a solvent structure

was also observed in our recent simulations of a polypeptide [12]. It should be pointed out that the distributions of H w around hyd,roxyl oxygen atoms have two peaks at 1.8 and 3.2 A (Fig. 2), in which the first peak corresponds to water hydrogens do-

0.4"

0.4

~ 5

~

ps

0.3.

~

10p6

0.3

0.2 '

0.1

o., 1

o-0 o

5

10

15

Distartee

20

0.0

25

* 0

5

10

15

20

25

20

25

(.,/IL)

0.4"

0.4

0.3"

0.3

~ 2 0

~

0.2 -

~0.2

0.1

0.0

0.1t

0

5

10

Distmae-.e

IS

20

(.~)

25

O.O.

5

10

Distance

15

(A)

Fig. 4. Calculated diffusion coefficients of water around maltose as a function of their starting distance to the closest solute atom.

C.X. Wang et al./ Chemical Physics Letters 251 (1996) 268-274

nated to the solute oxygen atoms for forming Hbonds. However, a broad peak of Hw around nonpolar atoms appears around 4.2 ,~ (Fig. 2, Table 3). But in some cases, for instance in Fig. 3, no prominent H w peak can be observed. Finally, the self-diffusion coefficient (D) of water has been calculated from the mean square displacement by using the Einstein relation lim At---* ~

([r(t+At)-r(t)]2)=6DAt,

(4)

where r(t) is the location of a water molecule at an initial time t. The value of D can be obtained from the asymptotic slope of the mean square displacement versus At. In order to obtain a more complete dynamic behavior of the solvent, the mean square displacements have been calculated for each water molecule and then averaged as a function of the distance from the closest solute atom for 1 ,~ shells. The results are shown in Fig. 4, in which different values of A t (5, 10, 15, 20 ps) in Eq. (4) have been selected to calculate D. We have averaged over all choices of the initial time t. The error bars are obtained from a standard statistical analysis. From Fig. 4 it is found that the overall line shapes of the four different cases are quite similar except the slight change occurring in the shells around 5 A. The diffusion coefficient is low at small distance from the maltose surface (0.15 .~2/ps for the first shell). This result indicates that the water molecules within a sphere of 5 A from the maltose surface have rather low mobility. On the other hand, the water molecules separated by more than 8 ,~ from the maltose surface have an approximate diffusion coefficient of 0.23 ,~2/ps (Fig. 4), which agrees with the experimentally determined value for pure water at 298 K [22] and is close to the value of 0.25 ,~2/ps found in the simulation of bulk S P C / E water at 300 K [14]. In order to compare our calculated value of D with the experimental data for maltose in solution, we have evaluated the average value of D for all the water molecules around maltose for 20 different simulation intervals (10 ps for each time period). The average value of D is 0.21 ,~2/ps, which is higher than the experimentally observed value (0.05 A2/ps) for 0.124 M solution [23]. It is found that the concentration of our model system (0.046 M) is much lower

273

than that of the experimental system, which might be the reason that the average value of D is higher than experiment. It is noted that in the previous simulation of [3-maltose in TIP3P water with a concentration similar to experiment, Brady and Schmidt [7] obtained an average diffusion coefficient of 0.09 A2/ps which is less than our result but higher than the experimentally observed value. This means that the calculated value of D around a carbohydrate molecule depends on the system concentration and the solvent model. To examine if the trend in the solvent mobility obtained from recent work [11,12] is reproduced in carbohydrates in aqueous solution, we have calculated the diffusion coefficients for water molecules within a 6 ~, sphere around the different kinds of solute atoms. For this case, the value of D was averaged over 20 time intervals in which 5 ps was taken for each time period. The values of D and their errors around polar and non-polar atoms are reported in the last column of Tables 2 and 3, respectively. It can be found that for each hydroxyl group the values of Dw~ar are somewhat larger than those of Dnonpolar. But there is no general trend for all cases except t h a t Dpola r and Dnonpolar a r e all smaller than the calculated value of Dbulk for bulk S P C / E water. Thus, no consistent conclusion with our recent work [12] can be drawn for the influence of different kinds of solute atoms on the solvent mobility. Because of the more flexibility and the similarity of hydroxyl groups in disaccharide, the dynamic behavior of disaccharide could be much different from that of protein and polypeptide molecules. Our result obtained from the present study suggests that the effect of neighboring enyironments and the particular structural property of a carbohydrate molecule could be a key point in the restriction of the mobility of the solvent by comparing the nature of specific atoms interacting with the solvent molecules.

4. Conclusions

In summary, the results described above demonstrate the following points: (1) The hydroxyl groups of maltose all participate in intermolecular hydrogen

274

C.X. Wang et al. / Chemical Physics Letters 251 (1996) 268-274

bonding with water. They form approximately eight H-bonds in total for each ring, in which the hydroxyl 0 6 group has the strongest H-bond with water. Other non-hydroxyl oxygen atoms have the lower average number of H-bonds with water as well as shorter H-bond life time. (2) The sharp peaks in the g(r) curves of water oxygen (O w) locate at 2.7-3.0 ,~, from polar hydroxyl oxygen, and the broad peaks of g(r) curves for Ow are found at 3.5-4.0 ,~ from non-polar atoms. The results of radial distribution function of water around maltose are consistent with the hydrogen bonding analysis described in the present and recent work [12]. The sharper the first peak in g(r) of O,~, the more H-bonds to water exist. (3) The calculated self diffusion coefficient of water molecules within a sphere of 5 ,~ from the surface of maltose has a relatively low value, which corresponds to lower mobility of water in this region. Examination of the diffusion coefficient of water around different types of atoms shows that no general trend can be found. This result indicates that the solvent mobility is not influenced significantly by the different atom types of maltose. This study will be continued with other simulations of oligosaccharide in solution to investigate the hydration behavior around the larger carbohydrate molecule.

Acknowledgement We thank M. Lahaye, H. Bizot, M.M. Delage and B. Leroux for helpful discussions. CXW and WZC thank INRA for the research grant and kind hospitality during their stay at the INRA Center of Nantes. This work was supported by the Chinese National Natural Science Foundation and the China-France Advance Research Project.

References [1] G.A. Jeffrey, in: Computer modeling of carbohydrate molecules, eds. A.D. French and J.W. Brady, ACS Symposium Series 430 (Am. Chem. Soc., Washington, DC, 1990), [2] J.E.H. Koehler, W. Saenger and W.F. van Gunsteren, J. Mol. Biol. 203 (1988) 241. [3] J.W. Brady, J. Am. Chem. Soc. 111 (1989) 5155. [4] C.J. Edge, U.C. Singh, R. Bazzo, G.L. Taylor, R.A. Dwek and T.W. Rademacher, Biochemistry 29 (1990) 1971. [5] S. Ha, J. Gao, B. Tidor, J.W. Brady and M. Karplus, J. Am. Chem. Soc. 113 (1991) 1553. [6] B.J. Hardy and A. Sarko, J. Comput. Chem. 14 (1993) 848. [7] J.W. Brady and R.K. Schmidt, J. Phys. Chem. 97 (1993) 958. [8] S.J. Weiner, P.A. Kollman, D.T. Nguyen and D.A. Case, J. Comput. Chem. 7 (1986) 230. [9] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey and M.L. Klein, J. Chem. Phys. 79 (1983) 926. [10] B.R. Brooks, R.E. Bruccoleri, B.D. Olafson, D.J. Sta~es, S. Swaminathan and M. Karplus, J. Comput. Chem. 4 (1983) 187. [11] R.M. Brurme, E. Liepinsh, G. Otting, K. Wi)thrich and W.F. van Gunsteren, J. Mol. Biol. 231 (1993) 1040. [12] C.X. Wang, R. Douillard and V. Tran, Chem. Phys. 189 (1994) 511. [13] W.F. van Gunsteren and H.J.C. Berendsen, Groningen Molecular Simulation (GROMOS) Library Manual (Biomos, Groningen, 1987). [14] H.J.C. Berendsen, J.R. Grigera and T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. [15] M.E. Gress, G.A. Jeffrey, Acta Cryst. B 33 (1977) 2490. [16] J.E.H. Koehler, W. Saenger and W.F. van Gunsteren, Eur. Biophys. J. 15 (1987) 197. [17] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola and J.R. Haak, J. Chem. Phys. 81 (1984) 3684. [18] J.P. Ryckaert, G. Ciccotti and H.J.C. Berendsen, J. Comput. Phys. 23 (1977) 327. [19] J. Tirado-Rives and W.L. Jorgensen, J. Am. Chem. Soc. 112 (1990) 2773. [20] C.X. Wang, Y.Y. Shi, F. Zhou and L. Wang, Proteins Struct. Funct. Gene 15 (1993) 5. [21] H.J.C. Berendsen, W.F. van Gunsteren, H.R.J. Zwinderman and R.G. Geurtsen, Ann. NY Acad. Sci. 482 (1986) 269. [22] R. Mills, J. Phys. Chem. 77 (1973) 685. [23] H. Ueadaira, Bull. Chem. Soc. Japan 42 (1969) 2140.