The SPASIBA force field as an essential tool for studying the structure and dynamics of saccharides

Biochimie 85 (2003) 65–73 www.elsevier.com/locate/biochi

The SPASIBA force field as an essential tool for studying the structure and dynamics of saccharides G. Vergoten *, I. Mazur, P. Lagant, J.C. Michalski, J.P. Zanetta UMR CNRS 8576 « Glycobiologie Structurale et Fonctionnelle », Université des Sciences et Technologies de Lille, 59655 Villeneuve d’Ascq, France Received 21 November 2002; accepted 3 February 2003

Abstract The SPASIBA force field has been applied to the determination of the structure and dynamical properties of various disaccharides. It has been shown that the experimental properties (structure, dipole moment, conformational relative energies) are satisfactorily predicted. The anomeric and exo-anomeric effects are confidently reproduced without specific terms for the a and b anomers and the type of glycosidic linkages. © 2003 Éditions scientifiques et médicales Elsevier SAS and Société française de biochimie et biologie moléculaire. All rights reserved. Keywords: SPASIBA force field; Saccharides; Molecular vibrations; Molecular dynamics

1. Introduction It is now clear that the electronic computer is one of the single most important inventions that have led to the most dramatic changes in both experimental and theoretical parts of scientific disciplines. “Computer experiments” now represent an additional aspect to be considered besides the classical experimental and theoretical divisions of scientific disciplines. Glycobiology and, in particular, the physical chemistry of glycans and glycoconjugates are no exceptions. To understand many of the functions and properties of glycans and glycoconjugates, it is necessary to possess a molecular level knowledge of their structure and dynamics. Such knowledge cannot be obtained using the usual experimental techniques due to their high level of flexibility and conformational adaptability. Computer simulations (molecular dynamics, Monte Carlo simulations) are some of the most promising tools for achieving this purpose. Within the field of numerical simulations, the use of computers is quite remarkable. Computer simulations have come to play a dominant role in modern science. They fill the gap between theory and experiment. In many cases, they represent the only means of

* Corresponding author. Tel.: +33-320337150; fax: +33-320436555. E-mail address: [email protected] (G. Vergoten).

making quantitative progress in the study of complex systems or of a simple system with complex behaviour. Computer simulation methods can be divided into deterministic (molecular dynamics) and probabilistic (Monte Carlo simulations). They both have intrinsic advantages and disadvantages [1,2]. They, however, have a common limitation: the availability of the potential energy function and its specific parameters. For several years, we have been developing a molecular mechanics force field which can be used with confidence particularly in molecular dynamics simulations, where it has to be reliable over the whole potential energy surface rather than in the vicinity of the minima. Such confidence can only be attained if a good fit between the observed and calculated vibrational frequencies and intensities is obtained. This force field has been called SPASIBA (Spectroscopic Potential Algorithm for Simulating Biomolecular conformational Adaptability) [3]. It has been parameterized for a large number of model compounds of biological importance, recently including phospholipids and saccharides. It has also been introduced in the CHARMM software [4]. 2. The SPASIBA potential energy function The molecular potential energy function V is not known. The only certainty we have on it is that it must be a function of the relative atomic displacement coordinates ri, the number of which is 3 N, where N is the number of atoms in the

© 2003 Éditions scientifiques et médicales Elsevier SAS and Société française de biochimie et biologie moléculaire. All rights reserved. DOI: 1 0 . 1 0 1 6 / S 0 3 0 0 - 9 0 8 4 ( 0 3 ) 0 0 0 5 2 - X

66

G. Vergoten et al. / Biochimie 85 (2003) 65–73

molecule. If these displacements are small enough, V can be expanded in a Taylor series as V共 r 兲 = V共 r0 兲 +

兺共 3N i

⭸V ⭸ri

兲 dr

i

+1 2

兺 冉⭸r⭸ ⭸rV 冊 ⭸r ⭸r + 3N i,j

i

j

j

2

2

2

and where s = r0共 1 − cos U0 兲/q0, t = r0共 sin U0 兲/q0

2

i

q = r1 + r2 − 2r1 r2 cos U

···

In this expansion, the first term corresponds to the equilibrium energy (i.e. the observed energy corrected by vibrational energy), r0 is the solution of the 3 N equation series

and r0, U0 and q0 are the equilibrium values of r1, r2, U and q. Since there is a redundancy relation among the four internal coordinates, the potential energy needs linear terms, as given in the following equation:

dV = 0 dri (equilibrium conformation). This indicates that a nonredundant set of displacement coordinates has to be used. The third term consists of the potential of infinitesimal vibrations around r0. The spectroscopic force field SPASIBA has been described in detail. It merges the AMBER [5,6] and UreyBradley-Shimanouchi [7] functional forms.

VAMBER consists of terms associated with torsional coordinates, Van der Waals and electrostatic interactions and hydrogen bonds. n and c are the multiplicity and phase of the dihedral angle φ, Vn is the energy barrier to internal rotation, r′ij is the non-bonded distance between atoms i and j, e is the dielectric constant, qi is the residual charge of atom i expressed in electronic unit, A and B are the Lennard-Jones parameters for each type of atom pairs, C and D are the parameters for the 10–12 potential for hydrogen bonds. VUREY − BRADLEY − SHIMANOUCHI has been introduced in order to impart “spectroscopic quality” to the force field. It is well known that cross-terms have to be taken into account in the potential energy function in order to reproduce the vibrational spectra. It is easily demonstrated that this purpose is achieved if the 1–3 non-bonded interaction in a valence angle is added and the corresponding redundancy problem is solved. Let us consider an H2O type molecule, where r1, r2, U and q are the two bond stretchings, the angle bending and the 1–3 non-bonded internal coordinates. Variations of these coordinates are related as seen in the following equation (redundancy relation):

which is derived from

The same dimensions are used for K′, H′, F′, K, H and F. The redundancy relation (Eq. (1)) is multiplied by Lagrange’s multiplier k and is added to Eq. (2):

Now the linear terms should be zero and we have

and

where F′q0 is the so-called gem internal tension. It is seen from Eq. (5) that if we introduce the 1–3 interaction force constant F and solve the redundancy problem, bond stretching-bond stretching and bond stretching-bond angle bending cross-terms clearly occur. The Urey-Bradley-Shimanouchi potential energy of the H2O type molecule is expressed by four force constants K, H, F and F′. The values of F and F′ are not independent. When the repulsion energy is expressed by V = a/q9, we have F′ = –F/10. When V = a/q12, we have F’ = –F/13. When the interaction energy is given by the Lennard-Jones 6–12 potential, we have F’ = –F[1–(q/qe)6]/13[1–(7/13)(q/qe)6], where

G. Vergoten et al. / Biochimie 85 (2003) 65–73

qe is the distance q at the position of minimum energy. From Eq. (5), it is also demonstrated that the linear terms in the potential energy function due to the redundancy relation among internal coordinates are eliminated when the gem internal tension is introduced to the second order. This procedure has to be performed each time the redundancy conditions are found among internal coordinates. If the basis redundancy is

67

This system of four equations has a solution if the determinant is equal to 0.

f共 Ui,j,Uik,... 兲 = 0 a function W of the first and second redundancies is defined: W = W ′ + W″ where W′ =

兺⭸U⭸f 共 DU ij

In order to eliminate the linear terms in the Urey-BradleyShimanouchi force field, the following term is added to the potential energy function:

ij 兲

and

″ V共 tetra 兲 = 1 R Kappa W tetra 2

where The corresponding potential energy term is ″ V = 1共 f.c. 兲W 2

where f.c. is the internal tension force tension (gem in the former example). Redundancy relations also exist in tetrahedral groups. For the CH4 type molecule, the redundancy relation has the form:

By differentiating A共 Uij 兲 the coefficients of the function W’ are obtained 2

aij = ⭸A bij = ⭸ A2 ⭸Uij ⭸U ij 2 2 k kl cij = ⭸ A d ij = ⭸ A ⭸Uij ⭸Uij ⭸Uij ⭸Ukl

DU12 + DU13 + DU14 + DU23 + DU24 + DU34 = 0 where indices 1–4 are associated with the four equivalent atoms. For the CH3X type molecule, we have a1共 DU12 + DU13 + DU14 兲 + b1共 DU23 + DU24 + DU34 兲 = 0 where the index 1 refers to atom X and indices 2, 3 and 4 refer to the three equivalent atoms. For the CH2X2 type molecule, we have a2 DU12 + b2共 DU13 + DU14 + DU23 + DU24 兲 + c2 DU34 = 0 where 1 and 2 refer to the H type atom and 3 and 4 refer to atoms X. To solve the general tetrahedral redundancy problem, four unit vectors s s s es1, e 2, e 3, e 4 along the four bonds are considered. In the 3D space, four vectors are linearly dependent: aes1 + bes2 + ces3 + des4 = s 0

It has been claimed [8] that if the Cartesian coordinates are used, the problems of redundant coordinates are eliminated. This is true only if the potential energy function is expanded in a Taylor series in terms of these coordinates and if Cartesian force constants are determined. Of course, chemists prefer internal coordinates. In the internal coordinate space: V=

兺

a=1

Fa Ra + 1 2

n

兺F

a,b = 1

ab

Ra Rb + ···

In the 3N Cartesian displacement coordinate space V=

3N

3N

兺 f x + 12 兺 i=1

i

i

i,j = 1

fij xi xj + ···

There is a relation between the two sets of force constants: fi =

n

兺b

a=1

a i

Fa fij =

n

兺b

a=1

a ij

Fa +

where

(6)

To determine the coefficients a, b, cand d, four vectors are multiplied to Eq. (6)

n

bi = a

and

⭸Ra ⭸xi

n

兺b

a,b = 1

a i

b

b j Fab

68

G. Vergoten et al. / Biochimie 85 (2003) 65–73 2

b ij = a

⭸ Ra ⭸xi ⭸xj

The linear internal coordinate force constants contribute to the quadratic Cartesian force constants through the coeffia cients b ij VUREY − BRADLEY − SHIMANOUCHI consists of four terms:

The trans, gauche and angular interactions have been described elsewhere [3]. The SPASIBA force field has been parameterized in our group for all the non-aromatic amino acids and also for a large number of model compounds such as alkanes, alkenes, aliphatic acids, aliphatic ethers, alcohols, esters, mono- and di-saccharides. This force field is able to reproduce at the same time as the molecular structures, conformational energy differences, wavenumbers of vibrational normal modes in addition to their potential energy distribution (PEDs).

3. The SPASIBA force field as introduced into the CHARMM software In the CHARMM program, the SPASIBA force field is activated using a particular parameter file including the “SPAS” keyword. The CHARMM parameters are used by default when the corresponding SPASIBA value cannot be found. The program was implemented in such a way that different values for the bond stretching, valence angle bending, torsional or improper torsional parameters could be used for the same kind of internal coordinate according to the nature of the residues. Introduction of a slow subroutine (ebaspas.src) into the CHARMM program enables the calculation of the specific SPASIBA bonding and in plane bending terms, including the potential energy and its first and second derivatives with respect to internal coordinates. The specific parameters are entered via a parameter file modified in order to include the K, H and F terms and the specific kappa, lCH2, trans and gauche force constants which are calculated for every potential energy evaluation. The kappa, lCH2, trans and gauche parameters are primitively evaluated in the “setup_spa.src” subroutine.

The determination of the non-bonded and electrostatic potential energies have not been modified in order to maintain specificities of the CHARMM program. Potential energy minimization and normal mode analyses have been performed for small molecules which are characteristic of several chemical groups, i.e. peptides (Nmethylacetamide, N-acetyl-L-alanine), alkanes, ethers (methylethylether), thiols (ethanethiol), esters (methylacetate), alkenes, alcohols (methanol, ethanol) and saccharides. As possible, a direct comparison of the calculated wavenumbers using the original CHARMM and SPASIBA (as introduced into CHARMM) force fields, together with the PED, was performed. Methylacetate will be given as the first example. Esters form an important class of molecules for which mechanical studies can be extended to check the effect of their insertion into the ester groups of glycolipids. Chhiba et al. [9]determined the associated SPASIBA force field parameters for several molecules taken as models (methyl and ethyl formate, methylacetate) and compared the dipole moments and moments of inertia to the experimental values. In the present work, we have directly transferred the parameter values given by these authors for bonds, valence angles and torsions. The CHARMM program was modified to accept n-fold terms for improper torsions and we used those values given by Chhiba et al. (n = 2 twofold rotational barrier). The topology file (OET) was created and the associated SPASIBA force constants have been included in the specific parameter file. Methylacetate has a hindered rotation around the O=C-O bond and experiments have shown that the most stable structure is obtained for the planar trans structure (with an energy difference of 8.5 kcal mol–1 relative to the cis conformer) [10]. In their work, Chhiba et al. used atomic charges derived from DFT methods (Becke3LYP/6–31G*); however, in the present work, we retained the values given in the CHARMM22 (C26a2) library. Table 1 displays the vibrational frequencies and PED for the original SPASIBA force field, for SPASIBA introduced into CHARMM and, as a matter of comparison, for the original CHARMM22 force field. Optimized geometries obtained in this work are quite comparable to the values given by Chhiba et al., but the dipole moment was calculated here to be 1.45D (1.68D for original SPASIBA force field). Examination of the various atomic charges used for each specific force field can easily explain the differences obtained for the two series of calculations; however, the relative root mean square (RMS) deviations show a good transferability for the SPASIBA force field. The second example is ethanethiol (C2H5SH). For ethanethiol, three molecular conformations can be put in evidence with staggered conformations about the C-S bond (two gauche enantiomers g, g’ with C1 group symmetry and one trans (t) conformer with Cs symmetry). Smith et al. [11] have shown that at room temperature, the infrared spectra cannot distinguish between the three conformers. These authors performed infrared and normal mode analyses in order to evaluate the relative contribu-

G. Vergoten et al. / Biochimie 85 (2003) 65–73

69

Table 1 Calculated vibrational wavenumbers (cm–1) for methylacetate and PED for original SPASIBA, SPASIBA in CHARMM and original CHARMM force fields Experiment [9] 110 136 187 303 429 607 639 844 980 1036 1060 1159

SPASIBA [9] 117 137 184 312 432 604 625 849 993 1060 1099 1158

SPASIBA into CHARMM 100 s C–CH3 109 s C–CH3 184 s C–O 303 d COC + d CO 418 d CCO + d CO 610 c C=O 628 m C–C + d CO 846 m C–O + m C=O + d COC 983 q CH3ac + m C–O + m C=O 1058 q CH3ac + c C=O 1097 m C–O + m C–CH3 1149 q CH3e

1187 1248 1375 1430 1430 1440 1460 1460 1771 2964 2966 2994 3005 3031 3035 RMS (cm–1)

1169 1278 1357 1454 1457 1459 1478 1479 1771 2940 2945 3003 3006 3008 3027 15.4

1171 1262 1348 1439 1452 1454 1471 1476 1770 2944 2964 3004 3006 3021 3023 15.08

q CH3e + d COC + m C=O + m C–O q CH3ac + m C–O + m C–CH3 ds CH3ac + c C–O ds CH3e da CH3ac da CH3ac da CH3e da CH3e m C=O ms C–H(CH3ac) ms C–H(CH3e) ma C–H(CH3ac) ma C–H(CH3ac) ma C–H(CH3e) ma C–H(CH3e)

Original CHARMM 28 ? 96 s C–O + g C=O 221 q CH3e 188 d COC 361 d CC–O + d CC=O 566 c C=O 495 m C–C + d CO + d CCO + m C–O 598 m C–O + m –C + r CH3 871 q CH3ac + c C=O 1103 d CH3e 980 m O–CH3e + d CCO + q CH3 1114 q CH3ac + q CH3e +m C–O + d CCO + d OCO 1135 q CH3e 1285 q CH3ac + m C–C ? 1365 d CH3ac 1366 d CH3ac 1415 da CH3e 1417 ds CH3e 1545 da CH3e 1722 m C=O 2787 ms C–H(CH3ac) 2793 ms C–H(CH3e) 2899 ma C–H(CH3ac) 2901 ma C–H(CH3ac) 2902 ma C–H(CH3e) 2905 ma C–H(CH3e) 104.6

ac, acetate; e, ester; m, stretching; ms, ma, symmetric and antisymmetric stretching; d, in plane bending; ds, dd, symmetric and antisymmetric in plane deformations; q, rocking deformation; c, out of plane deformation; s, torsion.

tions between the gauche and trans conformers. They studied the ethanethiol molecule in different physical states (gas, liquid and crystal states). Using calorimetric data and infrared experiments, these authors found the C1 (g and g’) as being the most stable by 0.2–0.3 kcal mol–1 with respect to the trans conformer value.Table 2 shows the calculated vibrational frequencies and PED among the internal coordinates (gauche), both for original SPASIBA (RMS = 7 cm–1) and for the SPASIBA in CHARMM. In the present work, the electrostatic potential was calculated using atomic charges derived for DFT calculations using the B3LYP functionals and a * * 6–31G basis set. Molecular dynamical studies using pure SPASIBA force field and the Verlet algorithm were performed on ethanethiol under the gas phase at 300 K to study the relative proportions in the g, g′ and t conformers vs. time. The relative amount given in percent for each conformer is displayed in Table 3 (starting from a gauche conformer). After 60 ps, the gauche conformer proportion decreases from 1.0 to 0.23, while the g’ conformer appears to an extent of 0.75. From 90 to 150 ps, the g conformer becomes the predominant form until the appearance of the three conformers in equal proportions. The energy difference obtained from the SPASIBA force field between the g (g’) and t conformers was calculated to be

0.2 kcal mol–1; this value being of the same order as the experimental value obtained from calorimetric data. This result clearly shows that this small difference would lead to the appearance of the same relative amounts of each conformer and to the predominance of the gauche forms.

4. The SPASIBA force field of the various glycosidic linkages Table 4 gives the refined force constants obtained from normal mode analysis performed on di- and tri-saccharides with various glycosidic linkages. Atomic charges were obtained for aa-trehalose (Fig. 1) (with a dipole moment of 1.75D) from DFT/B3LYP calculations using a 6–31G* * basis set. They were also used in other disaccharides. The present minimization procedures used at each step were performed on the disaccharide molecules using increments of 15° for the constrained phi (H–C–O–C′) (UH) and psi (C–O–C′–H′) (WH), while retaining all the remaining degrees of freedom. Hydrogen bonds have been described via electrostatic contributions and all energetic terms have been determined as currently used in the CHARMM program. Along each minimization, 5000 steps of steepest descent

70

G. Vergoten et al. / Biochimie 85 (2003) 65–73

Table 2 Vibrational wavenumbers (cm–1) and PED for ethanethiol (CH3-CH2-SH) Experiment (cm–1) [11] 191 246 331

Calculated (cm–1)a 191 251 335

658

653

738

742

871

872

971

980

1051

1052

1092

1057

1253

1255

1274 1379 1431 1449 1457 2571 2875 2902 2928 2967 2980 m(Debye) RMS (cm–1)

1277 1371 1433 1454 1457 2575 2879 2909 2916 2970 2971 1.45 7.3

PEDa sC–S(0.99) sC–C(0.99) dCCS(0.54), mC–S(0.27), mCC(0.07) mC–S(0.82), dCCS(0.05), dCSH(0.03) qCH2(0.45), mC–C(0.19), qCH3(0.17), sC–S(0.09) dCSH(0.57), mCC(0.25), mC–S(0.14) mC–C(0.64), dCSH(0.14), dCCH(0.05), mC–S(0.05) qCH3(0.28), mC–C(0.18), qCH2(0.12), dCSH(0.08) qCH3(0.56), qCH2(0.28), sC–C(0.10) swC2H2(0.61), dCH3(0.12), mC–S(0.07) wagC2H2(0.6), dCH3(0.33) dsCH3(0.95), mC–C(0.03) scC2H2(0.82) daCH3(0.91) daCH3(0.94) mS–H(1.0) ms CH2(0.84) ms CH3(0.87) ma CH2(0.98) ma CH3(0.99) ma CH3(1.0)

Calculated (cm–1)b 202 245 319

796

PEDb sC–S(0.97) sC–C(0.98) dCCS(0.78), mC–S(0.11), mC–C(0.04) mC–S(0.66), dCSH(0.17), dCCS(0.07) qCH2(0.70), qCH3(0.12)

766

dCSH(0.80), mC–S(0.15)

986

qCH3(0.37), mC–C(0.38), qCH2(0.07), mC–S(0.05) qCH3(0.71), qCH2(0.21)

636

1026

qCH3(0.46), mC–C(0.36), dCCS(0.06), dCSH(0.03) swCH2(0.94)

1040 1196

wagCH2(0.82), mC–C(0.09) dsCH3(0.90), mC–C(0.09) scCH2(0.85) daCH3(0.82) daCH3(0.95) mS–H(1.0) msCH2(0.99) maCH2(0.99) msCH3(0.99) maCH3(1.0) maCH3(1.0)

1374 1413 1447 1422 1428 2575 2851 2882 2902 2958 2959 1.45 19.1

m, stretching; d, in plane bending; s, torsion; sc, scissor; wag, wagging; tw, twisting; q, rocking. a From SPASIBA parameters. b From CHARMM parameters.

were required. A dielectric constant of 4 was used in each calculation; this value representing a correct estimation of the local environment acting on the disaccharide [12–14]. The cut-off for non-bonded atom pairs and electrostatic potentials were taken using a switch function (cutnb = 10 Å,

ctofnb = 8 Å, ctonnb = 7 Å) and the NBXMOD = 5 option in the CHARMM program. The 1–4 non-bonded and electrostatic interactions have been scaled to 0.5. The trehalose disaccharide is hereafter taken as an example. aa-Trehalose, ab-trehalose and bb-trehalose exhibit the 1→1 type of link-

Table 3 Ethanethiol: dihedral percent of gauche and trans conformers vs. time (ps) Angle 0–30 30–60 60–90 90–120 120–150 150–180 180–210 210–240 240–270 270–300 300–330 330–360

Time 25–50 ps 13.0 36.3 38.1 12.7

50–60 ps 25.7 24.0 24.5 8.

60–70 ps 7.2 3.2 3.8 8.0

70–80 ps 12.5 4.5 4.2 8.2

80–90 ps 7.2 3.5 4.2 8.0

90–100 ps 21.0 19.8 20.7 19.5

20.2 14.5 12.5 30.7

19.2 14.0 13.2 24.2

21.0 15.5 14.25 26.5

4.0 4.2 3.7 7.2

100–150 ps 9.35 39.9 43.5 7.3

150–200 ps 16.5 17.7 19.0 16.8 5.1 1.4 2.2 2.8 6.3 3.6 3.3 5.4

200–250 ps 8.1 7.3 7.8 7.9 9.3 8.9 9.6 9.0 8.55 7.8 7.5 8.2

250–500 ps 8.0 8.2 8.3 8.4 8.4 8.4 8.8 8.4 8.35 8.25 7.9 8.7

The atomic charges were derived from DFT (B3LYP/6–31G** calculations for the dynamical studies. C(CH3) –0.326, H(CH3) 0.1153, 0.1165, 0.1327, C(CH2) –0.318, H(CH2) 0.1464, 0.1431, S –0.081, H(SH) 0.072.

G. Vergoten et al. / Biochimie 85 (2003) 65–73

71

Table 4 SPASIBA parameters for the glycosidic linkage (O5,H,C)–C–OG–C–(C,C,H) OG is the oxygen atom of the glycosidic linkage. Neighbouring atoms are given in parentheses. R0 and h0 are the equilibrium values for bonds and valence angles, respectively Stretching C–OG

K 253.0

r0 1.418

Bending H–C–OG C–OG–C C–C–OG O5–C–OG C–C–OG

H 17.0 32.0 21.0 30.0 21.1

h0 109.0 115.7 108.5 110.0 110.0

F 80.0 30.0 60.0 55.0 55.0

Torsion C–C–OG–C O5–C–OG–C H–C–OG–C

Vn/2 0.70 0.70 0.70

n 3 3 3

c 0.0 0.0 0.0

K in kcal mol–1 Å–2, r0 in Å, H in kcal mol–1 rad–2, F in kcal mol–1 Å–2, h0 in degrees. Barrier to internal rotation (Vn in kcal mol–1), order (n) and phase in the Fourier series (c).

age. The aa and bb conformers of trehalose have their phi and psi angles opposite in sign and the distances between the anomeric hydrogen are quite comparable. The aa form of trehalose was shown to have significant properties for the protection of organisms against dehydration (yeast cells, plants, fungal spores), as they are able to recover their metabolism upon hydration after severe dryness. The low energetic conformer of aa trehalose exhibits an “open state” in its anhydrous state allowing a reversible hydration. These features have been well studied using numerous experiments such as NMR and molecular dynamics [15], vibrational

analysis [16], thermodynamical and electrical studies [17,18] and neutron scattering experiments [19,20]. aa-Trehalose is an axial-axial conformer, and the global minimum in energy was found to lie in the –70° < UH < –20° and –60° < WH < –20° range [13], correctly overlapping the experimental data observed for this disaccharide and its analogues in the crystal state while predicting the global minimum at UH = –44°, WH = –44°. The present work predicts nearly the same region (–75° < UH < –30°, –75° < WH < –30°) and practically the same values for the global minimum (UH = –51.5°, WH = –51.6°) (Fig. 2).

Fig. 1. CHELP charges for aa-trehalose as determined from DFT (B3LYP/6–31G* *) calculations. The H1–C1–O1–C1and C1–O1–C1–H1 are the (UH, WH) dihedrals relative to the glycosidic linkage.

72

G. Vergoten et al. / Biochimie 85 (2003) 65–73

Fig. 2. aa-Trehalose (O-(a-D-glucopyranosyl)-(1→1)-a-D-glucopyranose).

The present phi-psi map, however, does not appear to be completely symmetrical. This may be due to problems in the minimization procedures but we did not attempt to presently improve the calculations. The two hydroxymethyl groups appear to have g’g conformations, while only g or g’ conformations were detected for the exo-cyclic hydroxyl groups for the global minimum.

5. Atomic charge effects When introducing the SPASIBA force field in the CHARMM program, the electrostatic and Van der Waals contributions were mostly retained as given by atomic charges in the CHARMM library. Electrostatic contributions (EEL) to the total strain energy (ET) may be important as found in the case of methylacetate (ET= –14.8 kcal mol–1, EEL = –15.7 kcal mol–1) and ethanethiol (ET = 18.7 kcal mol–1, EEL = 16.6 kcal mol–1). In other cases, the contributions remain small as for aa-trehalose (bb-trehalose, respectively) (ET = 11.4 (10.0) kcal mol–1, EEL = 0.7 (1.0) kcal mol–1). In the CHARMM library, no charges were available for atoms involved in the various glycosidic linkages. Therefore, electrostatic charges were computed using the Density Functional Theory at the B3LYP/6–31G** level.

6. Conclusion The SPASIBA force field is available for a large variety of chemical groups. Its specific functional form and parameters were introduced into the CHARMM modelling software. It has been shown that the refinement of the force field param-

eters using the spectroscopic data leads to a better approach of the various physical properties than those obtained from other available force fields. The SPASIBA force field has shown its ability to correctly reproduce the empirical potential energy for disaccharides as the experimental values for UH and WH dihedral angles can be retrieved using unified parameters and charges that are transferred directly in the prediction of the global minima for various linkages. For all types of linkages, some differences in the localization of the global minimum can be noted when comparing our results to other works [13]. This is particularly true for the WH region extends, which appear more limited in the present work than in other studies [13]. This may be due to the particular analytical form of the force field used in this work. On the other hand, no particular attention was paid to the starting conformations of the hydroxyl groups but their orientations would certainly play a central role in the global minimum localizations. This work will enable to refine the SPASIBA force field and to extend it to polysaccharides including explicit solvation. A better confidence may be attained in results given by molecular dynamics and Monte Carlo simulations performed on biomolecular systems containing peptides, saccharides, peptidoglycans and/or phospholipids.

References [1]

D.P. Tieleman, S.J. Marrink, H.J.C. Berendsen, A computer perspective of membranes: molecular dynamics studies of lipid bilayer systems, Biochimica et Biophysica Acta 1331 (1997) 235–240. [2] K. Kubica, Computer simulation studies on significance of lipid planar head orientation, Computers Chem. 26 (2002) 351–356. [3] P. Derreumaux, G. Vergoten, A new spectroscopic molecular mechanic force field. Parameters for proteins, J. Chem. Phys. 102 (1995). [4] Lagant, D. Nolde, R. Stote, G. Vergoten, M. Karplus, Increasing normal mode analysis accuracy: the SPASIBA spectroscopic force field introduced into the CHARMM program (submitted for publication). [5] S.J. Weiner, P.A. Kollman, D.T. Nguyen, A.D. Case, An all atom force field for simulations of proteins and nucleic acids, J. Comp. Chem. 7 (1986) 230–252. [6] S.J. Weiner, P.A. Kollman, D.A. Case, U.C. Singh, C. Ghio, S. Profeta, P. Weiner, A new force field for molecular mechanical simulations of nucleic acids and proteins, J. Am. Chem. Soc. 106 (1984) 765–784. [7] T. Shimanouchi, Force constants of small molecules, Pure Appl. Chem. 7 (1963) 131–145. [8] A. Warshel, The QCFF/PI+MCA program package efficiency and versability in molecular mechanics, Computers Chem. 1 (1977) 195–202. [9] M. Chhiba, F. Tristram, G. Vergoten, The SPASIBA force field of esters, J. Mol. Struct. 405 (1997) 113–122. [10] C.E. Blom, H.S.H. Gunthard, Rotational isomerism in methyl formate and methylacetate: a low-temperature matrix infrared study using thermal molecular beams, Chem. Phys. Lett. 84 (2) (1981) 267–271. [11] D. Smith, J.P. Devlin, D.W. Scott, Conformational analysis of ethanethiol and 2-propanethiol, J. Mol. Spect. 25 (1968) 174–184. [12] A.D. French, A.M. Kelterer, G.P. Johnson, M.K. Dowd, C.J. Cramer, Constructing and evaluating energy surfaces of crystalline disaccharides, J. Mol. Graph Modelling 18 (2000) 95–108.

G. Vergoten et al. / Biochimie 85 (2003) 65–73 [13] A.D. French, A.M. Kelterer, C.J. Cramer, G.P. Johnson, M.K. Dowd, A QM/MM analysis of the conformations of crystalline sucrose moieties, Carbohydrate Res. 326 (2000) 305–322. [14] M.K. Dowd, A.D. French, P.J. Reiley, Conformational analysis of the anomeric forms of sophorose, laminarabiose and cellobiose using MM3, Carbohydrate Res. 233 (1992) 15–34. [15] S.B. Engelsen, C. Monteiro, C. Hervé du Penhoat, S. Pérez, The dilute aqueous solvation of carbohydrates as inferred from molecular dynamic simulations and NMR spectroscopy, Biophys. Chem. 93 (2001) 103–127. [16] K.I. Akao, Y. Okubo, N. Sakawa, Y. Inoue, M. Sakurai, Infrared spectroscopic study on the properties of the anhydrous form II of trehalose. Implications for the functional mechanism of trehalose as a biostabilizer, Carbohydrate Res. 334 (2001) 233–241.

73

[17] F. Sussich, C. Skopec, J. Brady, A. Cesaro, Reversible dehydration of trehalose and anhydrobiosis: from solution state to an exotic crystal, Carbohydrate Res. 334 (2001) 165–176. [18] T. Matsuoka, T. Okada, K. Murai, S. Koda, H. Nomura, Dynamics and hydration of trehalose and maltose in concentrated solutions, J. Mol. Liquids 98–99 (2002) 317–329. [19] H.D. Middendorf, S. Magazù, C. Branca, A. Faraone, G. Maisano, P. Migliardo, V. Villari, Molecular dynamics of disaccharides by inelastic light scattering, Physica B. Condensed Matter 267–278 (2000) 526–527. [20] S. Magazù, R.E. Lechner, S. Longeville, G. Raisano, D. Mapolino, P. Migliardo, V. Wanderlingh, Diffusive dynamics in trehalose aqueous solutions by QENS, Physica B. Condensed Matter 276–278 (2000) 475–476.