DFTMD studies of glucose and epimers: anomeric ratios, rotamer populations, and hydration energies

Carbohydrate Research 345 (2010) 503–511

Contents lists available at ScienceDirect

Carbohydrate Research journal homepage: www.elsevier.com/locate/carres

DFTMD studies of glucose and epimers: anomeric ratios, rotamer populations, and hydration energies q U. Schnupf, J. L. Willett, F. Momany * Plant Polymer Research, USDA, ARS, National Center for Agricultural Utilization Research, 1815 N. University St., Peoria, IL 61604, USA

a r t i c l e

i n f o

Article history: Received 9 September 2009 Received in revised form 24 November 2009 Accepted 7 December 2009 Available online 4 January 2010 Keywords: Anomeric ratios Glucose Epimers COSMO DFTMD

a b s t r a c t Results are presented from density functional molecular dynamics (DFTMD) simulations, based on constant energy dynamics, of glucose and its cyclic form of 6-carbon epimers. Both in vacuo and an implicit solvent method (COSMO) were examined, including simulations of low-energy conformations of each molecule. Analysis of the DFTMD results includes the following: energies averaged over the simulation time, calculated anomeric ratios, hydroxyl and hydroxymethyl rotamer populations, and hydration energies. Hydrogen-bonding networks persistence times were examined, and the effects of solvation on rotamer populations were described. Anomeric ratios calculated from energy optimization of an ensemble of low-energy conformers are compared to those obtained from ensemble averages from molecular dynamics, with dynamics simulations giving populations in best agreement with experimental anomeric ratios. Ensemble results in vacuo were not in agreement with experimental anomeric ratios or hydroxymethyl populations, producing in some cases reversal of the a:b ratios. The difficulty in obtaining correct a:b ratios increases with the number of axial groups; the mono-axial epimers being best represented, epimers with two axial groups being more difficult, and the epimers with three axial hydroxyl groups being most difficult to analyze, the result of a large number of very strong hydrogen-bonding networks that form the ensemble of low-energy conformations in the multi-axial structures. Published by Elsevier Ltd.

1. Introduction 1.1. Background Glucose and its epimers are important biomolecules involved in a variety of processes such as supporting matrices, conversion to alcohol for energy, molecular recognition processes such as blood group incompatibility, glycoconjugate antibiotics, cell attachment and bonding, viral infections, antitumor agents, and control of many biological functions. In the bio-environment these molecules exist in a water environment, so the influence of solvent on flexibility and conformation must be of interest from a molecular point of view. In this work, structural and thermodynamic properties of the above-noted cyclic 6-carbon aldopyranoses are investigated by both optimization methods and density functional molecular dynamics (DFTMD). Figure 1 shows the conventional numbering system, and Figure 2 shows the different axial/equatorial positions of the hydroxyl groups for the eight epimers. The 4C1 chair confor-

q Names are necessary to report factually on available data; however, the USDA neither guarantees nor warrants the standard of the product, and the use of the name by USDA implies no approval of the product to the exclusion of others that may also be suitable. * Corresponding author. Tel./fax: +1 309 681 6362. E-mail address: [email protected] (F. Momany).

0008-6215/$ - see front matter Published by Elsevier Ltd. doi:10.1016/j.carres.2009.12.001

mation of the ring shown in Figure 1 is the lowest energy ring conformation by these DFT calculations applied throughout the epimer series (see Table 1). All epimers have an anomeric site with two different anomeric forms, denoted a:b, which differ only by the orientation (axial or equatorial, respectively, with respect to the ring Fig. 1) of the C1–O1 bond. Hydroxymethyl groups can adopt three orientations associated with the C5–C6 internal rotation, characterized by the dihedral angles O5–C5–C6–O6 and C4– C5–C6–O6, that can be gauche–gauche (gg), gauche–trans (gt), or trans–gauche (tg). Hydroxyl groups are designated as clockwise, ‘c’, or counterclockwise, ‘r’, around the ring counting from the C1 atom. The O6–H group is observed to rotate fairly freely, and during dynamics it behaves much like a rotating top. Experimental evidence1–14 suggests that the glucose abundance ratio of the axial a anomer with respect to the equatorial b anomer in aqueous solution at room temperature is 38:62%,3 corresponding to a free-energy difference of 0.4 kcal/mol favoring the b anomer. This number could be misleading by suggesting that one conformer of each anomer with this energy difference would be sufficient to give the anomeric ratio. However, there are in all the epimers many low-energy conformations, each with sufficiently low energy that it may contribute to the ensemble energy, and all low-energy conformers must be included in the analysis to obtain a correct anomeric ratio. Using cutting edge DFTMD calculations it is possible to answer important questions such as: What

504

U. Schnupf et al. / Carbohydrate Research 345 (2010) 503–511

O6 C6

C6

‘c’ O4

O3

C4

C5

C3 C2

C5 O5

anomeric carbon C1

C1 O1

O2

α -form

β -form (equatorial)

(axial)

‘r’ α -D-glucose,‘c’-gg (gauch+) conformation

β -D-glucose,‘r’-gg (gauch+) conformation

Figure 1. Shown are a ball and stick figures of D-glucose with conventional notions used. Rotation around the C6–O6 bond produces configurations labeled as gauch+, gauch , and trans with C5–C6–O6–H 60°,  60°, and 180°, respectively. Rotation around the C5–C5 bond produces configurations labeled as gg, gt, and tg with O5–C5–C6–O6  60°, +60°, and 180°, respectively.

Figure 2. Shown are a ball and stick figures of D-glucose epimers.

is the source of the anomeric ratio in water, is it from the solvent, an entropic contribution, or both? Of further interest is the flexibility of these molecules, and DFTMD gives explicit information on the flexibility and stability of specific conformations in water. Comparing results from vibrational analysis of DFT optimized in vacuo and implicitly solvated structures with the results from the DFTMD, one can find answers to questions concerning the source of the anomeric ratios, rotamer populations, conformational flexibility, and solvent-directed conformational effects.

There have been many experimental NMR and circular dichroism (CD) studies1–13 and a host of computational studies of glucose and its epimers at various levels of theory.14–18,11,19–96 Several individual studies will be noted here for background information. Isolated glucose and epimers have been examined computationally starting in the late 1960s and 1970s,13–18,11,19,20 and computational studies continue to be published today. In most cases, but not all, empirical force-field calculations have been carried out. Quantum mechanical (QM) studies at various levels of theory are also found throughout the literature.20,24,29,33–35,40,42,43,45,47,48,52,53,56,58–60,63,66,70,72–80,85,86, 91–94 Early work generally used energy optimization of several conformers, expanding to a few empirical dynamics simulations of glucose in a box with explicit water molecules, and in later years to quantum optimization of a few conformations. The ‘first principle’ molecular dynamics study of glucose by the Car–Parrinello method53,58 was the first to try treating both solute and solvent by quantum methods. In those works the BLYP gradient corrected functional and Troullier–Martins norm-conserving pseudopotentials were utilized. The electronic wave function was expanded in a plane-wave basis set with cutoff criteria. Constant energy dynamics were performed at room temperature for 6 ps. The hydroxymethyl group in all studies remained at the starting orientation for 3 ps of the simulation, after which it moved to the tg conformation and remained there throughout the remaining simulation. As will be seen, this result is not in agreement with the DFTMD studies reported here, nor for that matter, is it agreement with experimental NMR studies.3 An earlier QM optimization approach by Polavarapu and Ewing29 examined only the ‘r’ glucose conformations, and then only a small subset of low-energy in vacuo structures. This study29 did not use DFT or solvation methods, and carried out calculations using the 4-31G and 6-31G* basis sets. Many different empirical molecular dynamics simulations have been carried out, with Brady’s work23 in 1986 being the first in which solvation was included in the dynamics. At that time the force field used was that of Rasmussen and his co-workers.22 Although this work23 included both the 4C1 and 1C4 ring structures, today we understand that the 1C4 relative energy is, for the most part, too high to be of interest at room temperature (see Table 1). The favored hydroxymethyl form found in that work was tg, and the hydroxyl orientations favored the ‘r’ form with some rotations to the ‘c’ form taking place during the 20 ps simulation.23 Some studies found that once the ring transitioned into a boat or skew form, it never returned to the 4C1 form, a result of empirical potentials that were unable to drive the boat structure back to the chair

505

U. Schnupf et al. / Carbohydrate Research 345 (2010) 503–511 Table 1 Listed are the relative energies of low-energy glucose epimers with different ring conformations calculated at the B3LYP/6-311++G** level of theory Relative energiesa (kcal/mol)

Axial OH

a-D-Glucose b-D-Glucose a-D-Mannose b-D-Mannose a-D-Allose b-D-Allose a-D-Galactose b-D-Galactose a-D-Altrose b-D-Altrose a-D-Gulose b-D-Gulose a-D-Talose b-D-Talose a-D-Idose b-D-Idose

(none) (none) (2) (2) (3) (3) (4) (4) (2,3) (2,3) (3,4) (3,4) (2,4) (2,4) (2,3,4) (2,3,4)

4

Vacuum C1

1

Vacuum C4

Vacuum Skew/boat

1.17 2.13 1.97 2.49 1.34 1.60 0.58 2.77 1.93 3.72 2.52 2.46 0.43 0.00 2.43 1.42

6.52 7.65 5.40 3.92 5.55 3.95 6.91 6.93 5.71 5.91 5.45 4.19 6.14 2.46 5.74 2.88

6.49 9.89 6.13 7.25 4.63 7.34 7.40 4.57 6.19 8.81 4.84 6.12 7.90 5.00 5.39 4.15

(3S5) (1S5) (B14) (1S3) (3S5) (1S5) (3,OB) (B3O) (3S5) (B3O) (3,OB) (1S5) (B14) (14B) (3,OB) (1S5)

Solution C1

4

0.82 1.00 1.61 1.98 1.16 1.37 1.22 1.42 2.87 2.88 2.11 2.52 0.35 0.00 1.97 1.29

Solution values are calculated using COSMO in conjunction with B3LYP/6-311++G**. a Relative energies are calculated with respect to the lowest vacuum energy of 431354.6 kcal/mol and lowest COSMO energy of

conformation. In a molecular dynamics study carried out in 1990,28 Kroon-Batenburg and Kroon found that the tg conformation of the glucose hydroxymethyl group was relatively unstable in water, appearing only as a minor conformational species. Their force field was GROMOS with CH and CH2 groups treated as united atoms. In 1993 further GROMOS32 work was carried out on glucose in a water bath, with reasonable results found for the, a:b ratio. Also in 1993 a QM optimization (4-31G and 6-31G*) with implicit solvation study was carried out33 on six glucose conformations, three a, and three b conformations. For these structures, only the ‘r’ form was examined, but an implicit solvation term from the SM2 model was included.33 A continued study using the AM1 and PM3 levels of theory that included 81 possible rotamers concluded that a-ggr was the lowest energy conformation. The authors33 concluded that there was no solvation effect on the a-D-glucose free energy of mutarotation, the final calculated value being  0.5 kcal/mol, which is compared to the experimental value of +0.3 kcal/mol. This result is not in agreement with the solvent effect presented here. These authors did suggest that solvation played a role in determining the hydroxymethyl rotamer populations, the solvation energy favoring the gt-r conformation. Since the 1994 empirical MM3 work of Dowd et al.,39 advances in computing power and software permit consideration of the monosaccharide conformational properties, including rotamers and ring shapes, at a much higher level of theory and without resort to any empirical potentials. In 1995 a semiempirical AM1 study43 using water as solvent showed that under those conditions the anti conformation (‘r’) of the exocyclic hydroxyl groups was favored. They suggested that water reduces the oscillation amplitude and lowered the barriers to rotation of the hydroxymethyl group. The year 1996 produced more carbohydrate studies.47–50 In 1997,52 ab initio studies (HF/ 6-311++G**//HF/6-31G**) of D-glucose and D-galactose were used to parameterize the AMBER force field. In the same year, molecular dynamics studies of D-glucose and D-galactose were carried out using the GROMOS force field.54 Molecular dynamics of both anomers of D-glucose61 were carried out in 1999 using GROMOS and the SPC/E water model to examine the glass transition temperature (calculated value 241 K vs experimental values of 232 K) of this solution. Hydration numbers and radial distribution functions for the oxygen atoms were examined. In the same year, both b-D-glucose and b-D-mannose62 were studied by molecular dynamics using empirical potentials.

Favored conformation vacuum/solution

4

C1 C1 C1 4 C1 4 C1 4 C1 4 C1 4 C1 4 C1 4 C1 4 C1 4 C1 4 C1 4 C1 4 C1 4 C1 4 4

431371.9 kcal/mol, respectively.

In 2002 a DFT optimization study73 using B3LYP/6-31G** was carried out on 264 conformers of methyl D-glucopyranoside and methyl D-galactopyranoside derivatives. Single point energies were obtained at the B3LYP/6-311++G** level of theory. Solvent was included using a self-consistent reaction field method. They concluded that the population of gg, gt, and tg rotamers were sensitive to solvent effects, and that there were large differences between the three hydroxymethyl populations of the in vacuo versus solution phase. Hydroxymethyl conformational properties of a-D-mannose2 by NMR spectroscopy and CD were described in 2003. In 2004, this laboratory published76,77 a series of DFT studies of a- and b-D-glucose and monohydrates of glucose with explicit water molecules. In 200579 structural properties of the glucose pentahydrates using high-level DFT methods were published. During the 2005–200780,86,89 period, mannose, galactose, and allose structures were examined in a series of DFT papers. A metadynamic study of the free energy landscape of glucose was carried out in 2009.96 As in a previous study20 the plane-wave method did not find hydroxymethyl populations in agreement with experiment.88 Finally, in 2009,93 glucose and the monohydrates were revisited while examining the OH stretch region, adding more glucose conformations to the earlier list. All the above were studied at the B3LYP/6-311++G** level of theory, and in some cases, the implicit solvation method, COSMO,97 was included. In the work presented here, all epimers of glucose have been studied at the same level of theory using density functional and basis sets described earlier. 2. Methodology 2.1. Computational methods Brief procedures for finding low-energy conformations have been published previously,80,86,89 and some expansion of the methods used is provided here. As before for glucose,76 mannose,80 galactose86, and allose,89 the initial chair, boat, and skew conformations were constructed using constraining potentials applied after optimizing all hydrogen-bond possibilities. All three low-energy hydroxymethyl group rotamers were included in the search, with O6–H hydroxyl directions tested at every conformer. Hydroxyl groups were chosen to be clockwise, ‘c’, or counterclockwise, ‘r’, or oriented to best allow for hydrogen-bonding patterns,

506

U. Schnupf et al. / Carbohydrate Research 345 (2010) 503–511

which sometimes required splitting the ‘c’ and ‘r’ orientations. As noted previously, it is not necessary or practical to consider all 729 possible combinations of staggered orientations of the hydroxyl and hydroxymethyl groups for each puckered ring conformation, those without optimized hydrogen bonding being of high relative energy. All combinations of the above were included in the search. Brief pre-optimization using semi-empirical (PM3)98 potentials was carried out, and the results for final unrestricted optimization utilized B3LYP/6-311++G**. Convergence criteria for optimizations were set at 1  10 6 Hartree for the energy and gradients of less than 3  10 4 au were demanded for both optimization with COSMO97 and those without. Constant energy dynamics, using the Verlet algorithm,99 B3LYP/6-31+G* DFTMD100 was carried out on the low-energy conformations (a 2 kcal/mol relative energy cutoff was used for the selection of low-energy conformations) for 1.5–5.0 ps, at 300 K, depending on energy convergence, with COSMO.97 The conditions for COSMO97 included the solvent as water, dielectric permittivity of 78.39, a scaling factor of 0.981, and a solvent sphere radius of 1.385 Å. Free energy is obtained from in vacuo and solvated optimized structures from calculation of the vibrational frequencies and subsequent entropies, and from the DFTMD results at 300 K by energy averaging (omitting the first 100 fs of the dynamics run for equilibration). The objective of the DFTMD dynamics is to better simulate the contribution of the flexibility of the hydroxyl groups to obtain free energies. This approach in conjunction with short simulation times requires that a DFTMD run needs to be performed for each low-energy conformer. The reason for this is that rotamer transitions require simulation times in the order of ns to obtain a reliable population distribution. Interpretation of the DFTMD results can be complicated by conformational transitions either in the hydroxymethyl group rotamers or by ring transition into boats, skews, or transition-state structures. When conformational transitions to other states occur, it is required to re-run the DFTMD dynamics with a different starting point for the trajectory. Rotation about the C–O bonds (hydroxyls) have much lower barriers to rotation than the C–C bond, and take up multiple orientations when free to do so without strong hydrogen bonding constraints, and this leads to simulations in which the conformation may undergo a transition to an equivalent and previously studied structure. This also suggests that because of the sampling problems, even using DFTMD results, simulations must be carried out on all low-energy conformations (relative energy <2 kcal/mol), and even then these runs are examined for structural transitions in order to derive usable Gibbs free energies for specific states, and this includes the entropy term intrinsic in the averaged energy. Anomeric ratios and rotomer populations are calculated using relative free-energy values. Nomenclatures are those published previously (see 80,86,89).

3. Results 3.1. Conformations The number of final optimized conformations for glucose and all the epimers are listed in Table 2. The number of optimized conformations is many times larger than the number of DFTMD conformations chosen because of the choice of low relative energy cutoffs (i.e., <2.0 kcal/mol) for the dynamics simulations. This choice occurs because those conformations with higher energies do not contribute significantly to the observable populations. For example, the lowest energy 1C4 conformations (Table 1) do not play significant roles in the calculation of a:b ratios for all epimers, the one exception being b-D-talose where the relative energy of the lowest energy 1C4 form is 2.5 kcal/mol. Even in this case this energy will be only a small part of the population calculation. Also in Table 1 are presented the lowest relative energies of the skew or boat conformations, with almost none of sufficiently low relative energy to make a contribution to the anomeric populations, although, as will be described, some twisted conformations appear during the dynamics simulations. In cases where many moderate relative energy conformations are available, these may contribute to the ensemble and so are included. Examination of the relative energies from in vacuo optimization, demonstrates that at this level of theory, b-D-talose is the lowest energy structure with a-D-talose next in relative energy. Upon optimization, a-D-glucose is found to be over 1 kcal/mol higher in energy than b-D-talose, with b-D-glucose over 2 kcal/mol higher in relative energy (see Table 1). The change in relative energies of this epimer series when COSMO97 is included in the optimization using the same basis sets is also presented in Table 1. Although the talose anomers remain of lowest energy, both glucose anomers move closer to talose with smaller relative energy, and the relative energy difference between a and b forms of both have become smaller upon application of COSMO during optimization. The lowest energy glucose anomer remains the a-form, but the difference between anomers is smaller when solvent is applied. From this optimization study it is clear that even with optimization that includes solvent, a conformation with the a form of glucose remains of lower energy than the best b anomer conformer. However, it is necessary to include all low-energy conformations in a population analysis before a true population can be deduced. DFTMD studies were carried out at the B3LYP/6-31+G* level of theory and include contributions from COSMO simulating solvent interactions. After each femptosecond dynamics step, the COSMO contributions were recalculated. Only those low relative energy conformations were included in the ensemble for dynamics simulations. The numbers included for glucose and each epimer are shown in Table 2. Clearly, some epimers, such as talose, only have

Table 2 Listed is the number of unique epimer structures investigated in the current work. Empirical structures are optimized with PM3 and then further optimized at the B3LYP/6311++G** level of theory

Glucose Mannose Allose Galactose Altrose Gulose Talose Idose

Empirical structures (PM3)

Optimized DFT vacuum structures

Optimized DFT COSMO structures

IR spectra and thermo data (vacuum)

IR spectra and thermo data (COSMO)

DFTMD

540 112 138 152 116 121 127 259

489 112 138 152 116 121 127 259

489 112 138 150 116 121 127 207

489 112 120 88 116 121 127 259

36 22 14 17 30 19 10 25

14 12 15 22 14 12 8 10

Solvent effects were simulated at the B3LYP/6-311++G**/COSMO level of theory. DFTMD calculations were carried out based on selected low-energy epimer structures at the B3LYP/6-31+G* level of theory.

507

U. Schnupf et al. / Carbohydrate Research 345 (2010) 503–511 Table 3 Calculated and experimental a:b anomeric ratios for glucose and its epimers Axial OH

b-gg:b-gt:b-tg

45:35:20 54:34:12 26:64:10 55:45:0 31:59:10

32:30:18 11:19:7 16:18:4 21:17:0

13:5:2 43:15:5 9:46:6 34:28:0

(2)

VACUUM COSMO DFTMDb

37:63 60:40 66:34 65:35 67:33

40:16:44 44:41:15 29:44:27 55:45:0

9:4:23 33:19:8 14:32:20 36:29:0

31:12:21 11:22:7 15:12:7 19:16:0

(3)

VACUUM COSMO DFTMDb

79:21 38:62 36:64 15:85 16:84

32:55:13 38:54:8 41:52:7

25:44:9 16:18:4 20:16:0

7:11:4 22:36:4 21:36:7

(4)

VACUUM COSMO DFTMDb

77:23 35:65 44:56 32:68 34:66

54:15:31 26:63:11 21:57:22 20:60:20

51:14:11 11:20:4 7:28:10 6:19:6

3:1:20 15:43:7 14:29:12 14:41:14

f

(2,3)

VACUUM COSMO DFTMDb

75:25 33:67 44:56 39:61 40:60

46:14:40 28:44:28 33:35:32

38:8:29 5:16:13 15:16:13

8:6:11 23:28:15 18:19:19

g

(3,4)

VACUUM COSMO DFTMDb

49:51 30:70 36:64 15:85 13:87

92:6:2 64:27:9 80:16:4

47:2:0 16:12:2 26:7:3

45:4:2 48:15:7 55:9:1

h

(2,4)

VACUUM COSMO DFTMD Experimental Zhu et al.

22:78 43:57 64:36 59:41 59:41

88:9:2 50:50:0 88:9:3

21:1:0 23:20:0 59:5:0

67:8:2 27:30:0 29:4:3

(2,3,4)

VACUUM COSMO DFTMD Experimental Zhu et al.

9:91 15:85 42:58 52:48 47:53

92:6:2 59:36:5 92:7:1

8:1:0 7:6:2 39:2:1

84:5:2 52:30:3 53:5:0

c

D-Galactose

D-Altrose

D-Gulose

i

a-gg:a-gt:a-tg

80:20 37:63 38:62 38:62

d D-Allose

D-Idose

gg:gt:tg

VACUUM COSMO DFTMDb

D-Mannose

D-Talose

a:b

(None)

D-Glucose

a

e

a For glucose the estimated contribution from conformations in solution other then 4C1 is <0.1%. The calculated average solvation energy for a-4C1 is 18.8 kcal/mol and for b-4C1 is 20.5 kcal/mol, respectively. b For experimental values see Refs. 3 (Zhu et al.), 10,12–14 (Angyal), and 88 (Suzuki et al.) c For mannose the estimated contribution from conformations in solution other the 4C1 is <0.1% (mostly a-1C4). The calculated average solvation energy for a-4C1 is 19.8 kcal/mol and for b-4C1 is 19.6 kcal/mol, respectively. d For allose the estimated contribution from conformations in solution other the 4C1 is <2.00% (mostly 1C4). The calculated average solvation energy for a-4C1 is 19.4 kcal/ mol and for b-4C1 is 20.5 kcal/mol, respectively. e For galactose the estimated contribution from conformations in solution other the 4C1 is <0.01%. The calculated average solvation energy for a-4C1 is 18.3 kcal/mol and for b-4C1 is 20.9 kcal/mol, respectively. f For altrose the estimated contribution from conformations in solution other the 4C1 is <3.00% (mostly 1C4). The calculated average solvation energy for a-4C1 is 18.9 kcal/ mol and for b-4C1 is 20.5 kcal/mol, respectively. g For gulose the estimated contribution from conformations in solution other the 4C1 is <1.50% (mostly 1C4 and OS2). The calculated average solvation energy for a-4C1 is 18.9 kcal/mol and for b-4C1 is 20.5 kcal/mol, respectively. h For talose the estimated contribution from conformations in solution other the 4C1 is <0.02%. The calculated average solvation energy for a-4C1 is 19.9 kcal/mol and for b-4C1 is 19.6 kcal/mol, respectively. i For idose the estimated contribution from conformations in solution other the 4C1 is <8.00% (there is a significant contribution by skew/boat forms totaling 6% from 3,OB, 2 SO, 1S5, and OS2, with an additional 2% from 1C4). The calculated average solvation energy for a-4C1 is 20.3 kcal/mol and for b-4C1 is 20.3 kcal/mol, respectively.

a few conformations of low energy, while others such as gulose, allose, mannose, and glucose have many low-energy conformations to include in the dynamics simulations. In general, the relative energy cutoff was 2.0 kcal/mol since relative energies larger than that contribute little to the total population. In some cases such as galactose and allose, simulations took on a different question. That is, it was necessary to take specific time segments out of the run in order to obtain free energies for specific confor-

mations from dynamic averages, since the structures assumed, for short periods of simulation time, boat, skew, or transition-state conformations for the ring. Altrose and gulose (see Table 1) have the highest relative energies, being higher than the 3-axial epimer, idose. In Table 1 the effect of adding solvent to the calculations is clearly observed. For example, the glucose a:b free energy difference is 1.0 kcal/mol, a being favored in vacuo, while it changes to 0.2 kcal/mol, a

508

U. Schnupf et al. / Carbohydrate Research 345 (2010) 503–511

being favored when COSMO is included in the energy optimization procedure. This would suggest that solvent plays a role in favoring the b form, but that is not the whole story. Table 3 shows the calculated a:b ratio, the overall gg:gt:tg population, and the hydroxymethyl rotamer populations for the individual a form and b form. One issue arises with the experimental data found in the literature. It has been found that the individual hydroxymethyl populations can change significantly with different studies, and some data (see 88) from work carried out in the 1970– 1980s has recently been found to be significantly in error. Further, it is usual in experimental studies to use blocked hydroxyl groups to maintain solely a or b forms. This chemical modification could lead to shifting of the populations relative to the pure sugars, but is not considered in this work. The most recent experimental a:b ratios of all the epimers are found in the NMR study of Serianni and co-workers3 and references therein, and these values are in excellent agreement with some older experimental studies. In comparing the results from the DFT studies described here to the experimental results from solution NMR studies, it becomes necessary to examine several epimers in detail. For example, allose and gulose have experimental a:b ratios 15:85, respectively.3 Both a epimers have structural features that make computational analysis difficult, a result of a very tight hydrogen-bonding network between the 3-hyroxyl axial group, the 1-hydroxyl axial group, and the 2-hydroxyl group that points down under the ring, as shown in Figure 2. This configuration persists throughout 5 ps of DFTMD, and contributes substantially to favoring the a form energetically; thus our prediction of nearly 34:66 population in these two epimers. Further, it has become clear that the b anomer, but not the a forms of several epimers, transition easily into higher energy ring conformations that are not 4C1 chairs. These transitions complicate the dynamics runs by short transitory structural changes leading to difficulty in establishing sufficient simulation times to get reliable b anomer averages for one specific conformational state. When these transitions occur, the approach used is to rerun the starting conformation. Since the first few femptoseconds of the simulation are so important in establishing the stable conformational state, one can get a different set of conformations by just restarting the run. Examination of these cases suggests that possibly very long simulation times might be required in the two cases described above, or possibly the addition of explicit hydration (in conjunction with implicit hydration) would help achieve a:b ratios in better agreement with NMR results. 3.2. Individual simulations 3.2.1. D-Glucose From Table 3 it is apparent that the different methods to obtain a:b ratios and hydroxymethyl rotamer ratios give, with the exception of the in vacuo case, excellent comparisons to experiment. In particular, the newer work of Suzuki et al.88 shows better agreement with the calculated hydroxymethyl rotamers than the earlier work. As will be seen in every case, applying the COSMO solvation method to the calculations results in changing a:b preferences in the correct direction. The relative difference in hydration energy (see Table 4) favors the b anomer by 1.7 kcal/mol. The relative populations of hydroxymethyl rotamers are given in Table 3 for the vacuum, COSMO, dynamics, and two experimental values. The change when going from vacuum to COSMO is fairly significant, with the values from the dynamics averaging very close to the most recent experimental and theoretical values of Suzuki et al.88 It is of interest that the experimental values have changed significantly from the earlier work, most difference being in the observation of a larger tg population at the expense of the gg con-

Table 4 Calculated hydration energies for glucose and its epimers, in kcal/mol Epimer

DFTa

Axial OH

a

a b c

CHARMMb

D

b

c

a

D

b

D-Glucose

(None)

18.8

20.5

1.7

23.3

24.2

0.9

D-Mannose

(2)

19.8

19.6

+0.2

22.4

22.8

0.4

D-Allose

(3)

19.4

20.5

1.1

23.0

23.9

0.9

D-Galactose

(4)

18.3

20.9

2.6

22.2

22.6

0.4

D-Altrose

(2,3)

18.9

20.5

1.6

23.3

22.5

+0.8

D-Gulose

(3,4)

18.9

20.5

1.6

22.7

22.5

+0.2

D-Talose

(2,4)

19.9

19.6

+0.3

22.2

22.0

+0.2

D-Idose

(2,3,4)

20.3

20.3

0.0

21.6

22.6

1.0

Calculated at the B3LYP/6-311++G**/COSMO level of theory. See Ref. 81. D is the energy value for b minus that for a.

formation. The individual populations for a and b anomers differ with the b form from the dynamics having significantly higher values for the gt than found in the a form. 3.2.2. D-Mannose Again, excepting the vacuum study which has the calculated a:b ratio reversed, application of solvation, that is, COSMO, to both the static and dynamics studies results in improvement when compared to experimental data. The a anomer is favored by 0.2 kcal/ mol of hydration energy, and that is in line with the observed a:b ratio which favors the a anomer. Dynamic averaging gives excellent agreement with the observed anomeric ratio. The hydroxymethyl group population for mannose is fairly consistent for COSMO and dynamics averages, with gg and gt being highly populated, but with some tg appearing in contrast to the older experimental results. The comparison between a and b forms are in reasonable agreement with the experimental analysis. 3.2.3. D-Allose This epimer was of particular difficulty because of a triad of strong hydrogen bonds in the a anomer case. The result of this hydrogen-bonding network was to lower the a anomer energy such that neither the solvent (hydration favored b anomer by 1.1 kcal/mol) contribution nor the favorable entropy of the b form could overcome the relative energy gap, and thus the populations obtained were biased toward the a form. In the allose case the optimized COSMO and dynamics give nearly identical hydroxymethyl rotamer populations with the dynamics average. There appears to be only slight differences in population between a and b forms with gg and gt being nearly equally populated. No experimental populations are available. 3.2.4. D-Galactose One might expect that this epimer would be relatively easy to obtain an agreement with the experimental results, unfortunately there are serious difficulties in choosing the set of conformations one wishes to use to carry out the dynamics. The result shows that COSMO drives the system in the right direction when compared to the vacuum for the optimized conformations, and dynamics moves the ratio in the wrong way making the a anomer population somewhat too large. It is of interest that the hydration energy does favor the b anomer by 2.6 kcal/mol. which is in the correct direction when compared to the experimental values. The hydroxymethyl rotamer populations obtained from the dynamics simulations are in very good agreement with the experimental results. The axial hydroxyl group in the 4-position appears to be responsible for the 22% tg rotamer population.

U. Schnupf et al. / Carbohydrate Research 345 (2010) 503–511

3.2.5. D-Altrose As noted previously, the vacuum optimization studies do not reproduce the experimental a:b ratio, strongly favoring the a anomer. Addition of COSMO alters the ratio to favor the b anomer, with favorable agreement with experiment. The hydration energy from COSMO favors the b anomer by 1.6 kcal/mol, and this again is in line with the experimental results for the anomeric a:b ratio. The dynamics average values give results in good agreement with experiment. The hydroxymethyl rotamer population dynamics results suggest that, in both a and b anomers, all three rotamers are significantly populated and of nearly the same amount. This distribution of rotamers is unique in the epimers studied here. 3.2.6. D-Gulose This epimer is doubly axial in the 3- and 4-positions and has proven to be a difficult molecule to model. The experimental a:b ratio is predominately b favored, and the hydration energy favors the b anomer by 1.6 kcal/mol. Even so, it is insufficient to overcome the very low relative energies of a few a anomer conformations. In this case, application of COSMO to the optimized conformations does not fully describe the system when compared to those from the vacuum, although this application of solvent is in the correct direction. The dynamics averages are not in better agreement with the experiment, even though they do predict a strong b anomeric preference. The primary problem with this molecule, as well as others, is that the conformation about the C5–C6 hydroxymethyl bond is not stable, moving in one simulation from gt to tg to gg in the course of 5 ps. From this simulation one has only short time segments in which one hydroxymethyl conformation is stable, and this time is insufficient in order to get a reasonable average free energy for one rotamer. The solution to this is either to run the simulation for a very long time, which is very expensive, or to run many more short simulations until one is found that maintains the starting rotamer conformation throughout the run, providing a good average energy. Long runs pose a different problem in that for the multiply axial epimers, there is a high probability of ring distortions to the transition-state structure, boats, or skew forms appearing to make analysis more complicated. Neither solution is practical since these are expensive and time-consuming runs. Dynamics strongly favors the gg rotamer population, particularly in the b anomer, with only small populations of gt and tg appearing. Optimization studies with COSMO gave somewhat more favor to the gt form but less to the tg rotamer population. 3.2.7. D-Talose For this doubly axial epimer, solvation with COSMO to the optimized conformations gave a hydration energy that slightly favored the a anomer (0.3 kcal/mol). In this case the 1C4 ring conformation makes up <0.02% and will not interfere with the statistical analysis. As before, the vacuum calculations give almost the reverse of the observed anomeric ratio, while including COSMO in the static optimized structures reverses this affect and brings the anomeric ratio into better agreement with experimental values. Dynamic averaging appears to slightly overstate the a anomeric preference. There is a dramatic enhancement of the gg hydroxymethyl population when it is calculated from the dynamics runs over that calculated from the optimized structures with COSMO, whereby it gains over 30% in gg at the expense of the gt rotamer. 3.2.8. D-Idose This triply axial epimer (axial in the 2-, 3-, and 4-positions) models by dynamic averaging compared very well with the experimental results, producing a 42:58% ratio. However, the vacuum-

509

and COSMO-optimized conformations do not reflect the experimental values, as the b anomer is highly preferred. Dynamic averaging appears to correct this large difference by opening the molecule to other ring and hydroxyl conformations. It is of interest that in <6.00% of the dynamics simulation the structure took up skew/boat forms including 3,OB, 2SO, 1S5, and OS2 forms, with 2% being in the 1C4 conformation. The average solvation energy favors neither a or b anomers (energy difference of 0.0 kcal/mol). Hydroxymethyl rotamer population is, in this case, highly favored toward the gg conformation. 3.3. Hydroxyl rotamer populations During dynamics the hydroxyl groups undergo significant rotational motion jumping about from one position to another. In order to examine this motion, analysis was carried out for each position on the ring, and several of the more important rotations are described here. First it is clear from Table 1 that the anomeric configuration strongly influences the hydroxymethyl rotamer population. For example, glucose, and for that matter all the epimers exhibit significant changes in population that entirely reverse the populations from the gg state preferred to the gt state preferred, depending upon the presence or absence of solvent. In some cases the gg form may be of largest population in vacuo for the a anomer, while in solvent, the b anomer has the largest gg population. These changes occur time and again and show that solvent does play a significant role in determining the hydroxymethyl populations. Importantly, application of DFTMD again shifts the populations, sometimes reducing one form and increasing another, as shown for glucose where the optimized with COSMO overall results for both anomers gives a gg:gt:tg ratio of 54:34:12. But when DFTMD+COSMO is run the populations changed to 26:64:10, in close agreement with recent experimental88 results of 31:59:10. 3.4. Hydration energies There are no reliable experimental hydration energies for glucose and its epimers. However, some hydration values have been estimated using empirical potentials and explicit water molecules within a molecular dynamics simulation.81 Table 4 gives calculated results from CHARMM simulations81 and compares the values obtained here with results from our energy-averaged DFTMD methods. The agreement is generally good, with some differences on preferences depending upon the anomeric configuration. The empirical values are more negative in free energy of hydration than the DFT values by 2 kcal/mol, the differences being fairly consistent. The anomer with the largest negative hydration energy is most usually found to be the most populated in water solutions, and this is in agreement with our DFTMD hydration values. The most highly b anomer populated epimers are allose and gulose, and the DFT hydration predictions are consistent with this, while the empirical results81 show allose to also be b favored, but gulose to be a favored and idose to be highly b favored by hydration energy. 4. Conclusions In previous studies of glucose and several epimers, a large number of low-energy conformations have been described, and exhaustive tables of conformations have been presented. Here these ensembles of aldohexopyranose conformations are condensed, and results of many optimizations and DFTMD simulations are reduced to simplified tables of populations and hydration energies. As learned from this work, one can obtain very close agreement with experimental results of the monomers using these high-level

510

U. Schnupf et al. / Carbohydrate Research 345 (2010) 503–511

DFT methods with no variable parameters. The advantages of these calculations over empirical studies is now apparent, as particular empirical studies cannot obtain close relative energies over all the epimers in order to predict anomeric ratios or hydroxyl rotamer populations of this series. The advantage of being able to carry out DFTMD with implicit solvent is also apparent, since for the rotamer and anomeric populations to be correct, the solvent is absolutely necessary, and the relative energies with solvent must be very close to the correct values. The absolute values obtained for the solvation energies are probably only approximate; however, the trends in hydration energy when mapped do correlate with the experimental data, the hydration energy negative values pointing to the preferred solution populations in the a:b ratios. To answer the question of the source of the anomeric ratios, one must only look at the data. The free-energy results in vacuo were not in agreement with the experimental values, reversing the a:b ratio in some cases. In fact one must even go beyond optimization with COSMO to obtain results that reflect reality, that is, DFTMD + COSMO. The conclusions are that both solvent and entropy play roles in these energetic problems. In some cases where the complex hydrogen bonding is very strong, the calculated results are in greatest error. This could be a result of the short dynamics runs used or the possibility of having missed a low-energy conformation, but this latter is unlikely. Unfortunately, longer DFTMD runs pose a serious problem by producing complex structural transitions to boat or skew forms, particularly in the multiple axial cases and those that are highly b preferred in solution. It is a serious challenge to computational chemists to find ways to treat these systems, but the results can be productive when correctly carried out.

60. 61. 62. 63. 64. 65.

References

66. 67.

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. 30. 31. 32. 33. 34. 35.

Grindley, T. B. Chem. Chem. Biol. 2001, 1, 3–51. Nobrega, C.; Vazquez, J. T. Tetrahedron: Asymmetry 2003, 14, 2793–2801. Zhu, Y.; Zajicek, J.; Serianni, A. S. J. Org. Chem. 2001, 66, 6244–6251. Astley, T.; Birch, G. G.; Drew, M. G. B.; Rodger, P. M.; Wilden, R. H. J. Comput. Chem. 1998, 19, 363–367. Prdron, J. I.; Moreles, E. Q.; Vazquez, J. T. J. Org. Chem. 1998, 63, 8247–8258. Rockwell, G. D.; Grindley, T. B. J. Am. Chem. Soc. 1998, 120, 10953–10963. Morales, E. Q.; Padron, J. I.; Trujillo, M.; Vazquez, J. T. J. Org. Chem. 1995, 60, 2537–2548. Bock, K.; Duus, J. J. Carbohydr. Chem. 1994, 13, 513–543. De Vries, N. K.; Buck, H. M. Carbohydr. Res. 1987, 165, 1–16. Angyal, S. J. Adv. Carbohydr. Chem. Biochem. 1984, 42, 15–69. Marchessault, R. H.; Perez, S. Biopolymers 1979, 18, 2369–2374. Angyal, S. J.; Pickles, V. A. Aust. J. Chem. 1972, 25, 1695–1710. Angyal, S. J. Aust. J. Chem. 1968, 21, 2737–2746. Angyal, S. J. Angew. Chem., Int. Ed. Engl. 1969, 8, 157–166. Rao, V. S. R.; Vijayalakshmi, K. S.; Sundararajan, P. R. Carbohydr. Res. 1971, 17, 341–352. Kildeby, E.; Melberg, S.; Rasmussen, K. Acta Chem. Scand. 1977, A31, 1–13. Pérez, S.; St-Piere, J.; Marchessault, R. H. Can. J. Chem. 1978, 56, 2866–2871. Melberg, S.; Rasmussen, K. Acta Chem. Scand. 1978, A32, 187–188. Joshi, N. V.; Rao, V. S. R. Biopolymers 1979, 18, 2993–3004. Melberg, S.; Rasmussen, K. Carbohydr. Res. 1979, 76, 23–37. Tvaroska, I.; Kozar, T. J. Am. Chem. Soc. 1980, 102, 6929–6936. Rasmussen, K. Acta. Chem. Scand. 1982, A36, 323–327. Brady, J. W. J. Am. Chem. Soc. 1986, 108, 8153–8160. Tvaroska, I.; Kozar, T. Theor. Chim. Acta 1986, 70, 99–114. Brady, J. W. Carbohydr. Res. 1987, 165, 306–312. Ha, S. N.; Giammona, A.; Field, M.; Brady, J. W. Carbohydr. Res. 1988, 180, 207– 221. Brady, J. W. J. Am. Chem. Soc. 1989, 111, 5155–5165. Kroon-Batenburg, L. M.; Kroon, J. Biopolymers 1990, 29, 1243–1248. Polavarapu, P. L.; Ewig, C. S. J. Comput. Chem. 1992, 13, 1255–1261. Sekkal, M.; Legrand, P.; Vergoten, G.; Dauchez, M. Spectrochim. Acta 1992, 48A, 959–973. Daucher, M.; Derreumaux, P.; Vergoten, G. J. Comput. Chem. 1992, 14, 263– 277. Kouwijzer, M. L. C. E.; van Eijck, B. P.; Kroes, S. J.; Kroon, J. J. Comput. Chem. 1993, 14, 1281–1289. Cramer, C. J.; Truhlar, D. G. J. Am. Chem. Soc. 1993, 115, 5745–5753. Van Eijck, B. P.; Hooft, R. W. W.; Kroon, J. J. Phys. Chem. 1993, 97, 12093– 12099. Forster, M. J.; Mulloy, B. Biopolymers 1993, 33, 575–588.

36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56. 57. 58. 59.

68. 69. 70. 71. 72. 73. 74. 75. 76. 77. 78. 79. 80. 81. 82. 83. 84. 85. 86. 87. 88. 89. 90. 91. 92. 93.

Hajduk, P. J.; Horita, D. A.; Lerner, L. E. J. Am. Chem. Soc. 1993, 115, 9196–9201. Grootenhuis, D. J.; Haasnoot, A. G. Mol. Simul. 1993, 10, 75–95. Adams, B.; Lerner, L. E. Magn. Reson. Chem. 1994, 32, 225–230. Dowd, M. K.; French, A. D.; Reilly, P. J. Carbohydr. Res. 1994, 264, 1–19. Salzner, U.; von Rague Schleyer, P. J. Org. Chem. 1994, 59, 2138–2155. Galema, S. A.; Howard, E.; Engberts, J. B. F. N.; Grigera, J. R. Carbohydr. Res. 1994, 265, 215–225. Tvaroska, I.; Gajdos, J. Carbohydr. Res. 1995, 271, 151–162. Zuccarello, F.; Buemi, G. Carbohydr. Res. 1995, 273, 129–145. Asensio, J. L.; Martin-Pastor, M.; Jiménez-Barbero, J. Int. J. Biol. Macromol. 1995, 17, 137–148. Barrows, S. E.; Dulles, F. J.; Cramer, C. J.; French, A. D.; Truhlar, D. G. Carbohydr. Res. 1995, 276, 219–251. Kouwuzer, M. L. C. E.; van Eijck, B. P.; Kooijman, H.; Kroon, J. Acta Crystallogr., Sect. B 1995, 51, 209–220. Jebber, K. A.; Zhang, K.; Cassady, C. J.; Chung-Phillips, A. J. Am. Chem. Soc. 1996, 118, 10515–10524. Brown, J. W.; Wladkowski, B. D. J. Am. Chem. Soc. 1996, 118, 1190–1193. Csonka, G. I.; Elias, K.; Csizmadia, I. G. Chem. Phys. Lett. 1996, 257, 49–60. Cheetham, N. W. H.; Lam, K. Carbohydr. Res. 1996, 282, 13–23. Jonsdottir, S. O.; Klein, R. A. Fluid Phase Equilib. 1997, 132, 117–137. Senderowitz, H.; Still, W. C. J. Org. Chem. 1997, 62, 1427–1438. Molteni, C.; Parrinello, M. Chem. Phys. Lett. 1997, 275, 409–413. Klewinghaus, P.; van Eijck, B. P.; Kouwijzer, M. L. C. E.; Kroon, J. J. Mol. Struct. (Theochem) 1997, 395–396, 289–295. Caffarena, E. R.; Grigera, J. R. Carbohydr. Res. 1997, 300, 51–57. Csonka, G. I.; Kolossvary, I.; Csaszar, P.; Elias, K.; Csizmadia, I. G. J. Mol. Struct. (Theochem) 1997, 395–396, 29–40. Wladkowski, B. D.; Chenoweth, S. A.; Jones, K. E.; Brown, J. W. J. Phys. Chem. 1998, A102, 5086–5092. Molteni, C.; Parrinello, M. J. Am. Chem. Soc. 1998, 120, 2168–2171. Barrows, S. E.; Storer, J. W.; Cramer, C. J.; French, A. D.; Truhlar, D. G. J. Comput. Chem. 1998, 19, 1111–1129. Ma, B.; Schafer, H. F., III; Allinger, N. I. J. Am. Chem. Soc. 1998, 120, 3411–3422. Caffarena, E. R.; Grigera, J. R. Carbohydr. Res. 1999, 315, 63–69. Roberts, C. J.; Debenedetti, P. G. J. Phys. Chem. 1999, B103, 7308–7318. Gregurick, S. K.; Kafafi, S. A. J. Carbohydr. Chem. 1999, 18, 867–890. Andersson, C.; Engelsen, S. B. J. Mol. Graph. Model. 1999, 17, 101–105. Spieser, S. A. H.; van Kuik, A.; Kroon-Batenburg, L. M. J.; Kroon, J. Carbohydr. Res. 1999, 322, 264–273. Hoffman, M.; Rychlewski, J. Comput. Methods Sci. Technol. 2000, 6, 61–64. Stenger, J.; Cowman, M.; Eggers, F.; Eyring, E. M.; Kaatze, U.; Petrucci, S. J. Phys. Chem. 2000, B104, 4782–4790. Palma, R.; Himmel, M. E.; Brady, J. W. J. Phys. Chem. 2000, B104, 7228– 7234. Fuchs, K.; Kaatz, U. J. Phys. Chem. 2001, B105, 2036–2042. Lewis, B. E.; Schramm, V. L. J. Am. Chem. Soc. 2001, 123, 1327–1336. Kirschner, K. N.; Woods, R. J. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 10541–10545. Csonka, G. I. J. Mol. Struct. (Theochem) 2002, 584, 1–4. Tvaroska, I.; Taravel, F. R.; Utille, J. P.; Carver, J. P. Carbohydr. Res. 2002, 337, 353–367. Lii, J.-H.; Chen, K.-H.; Allinger, N. L. J. Comput. Chem. 2003, 24, 1504–1513. Hemmingsen, L.; Madsen, D. E.; Esbensen, A. L.; Olsen, L.; Engelsen, S. B. Carbohydr. Res. 2004, 339, 937–948. Appell, M.; Strati, G.; Willett, J. L.; Momany, F. A. Carbohydr. Res. 2004, 339, 537–551. Momany, F. A.; Appell, M.; Strati, G.; Willett, J. L. Carbohydr. Res. 2004, 339, 553–567. Corchado, J. C.; Sanchez, M. L.; Aquilar, M. A. J. Am. Chem. Soc. 2004, 126, 7311– 7319. Momany, F. A.; Appell, M.; Willett, J. L.; Bosma, W. B. Carbohydr. Res. 2005, 340, 1638–1655. Appell, M.; Willett, J. L.; Momany, F. A. Carbohydr. Res. 2005, 340, 459– 468. Dashnau, J. L.; Sharp, K. A.; Vanderkooi, J. M. J. Phys. Chem. 2005, B109, 24152– 24159. Behrends, R.; Kaatze, U. ChemPhysChem 2005, 6, 1133–1145. Mason, P. E.; Neilson, G. W.; Enderby, J. E.; Saboungi, M.-L.; Brady, J. W. J. Phys. Chem. 2005, B109, 13104–13111. Gonzalez-Outeirino, J.; Kirschner, K. N.; Thobhani, S.; Woods, R. J. Can. J. Chem. 2006, 84, 569–579. Kurihara, Y.; Ueda, K. Carbohydr. Res. 2006, 341, 2565–2574. Momany, F. A.; Appell, M.; Willett, J. L.; Schnupf, U.; Bosma, W. B. Carbohydr. Res. 2006, 341, 525–537. Mason, P. E.; Neilson, G. W.; Enderby, J. E.; Saboungi, M-L.; Cuello, G.; Brady, J. W. J. Chem. Phys. 2006, 125, 224505-1–224505-9. Suzuki, T.; Kawashima, H.; Sota, T. J. Phys. Chem. 2006, B110, 2405–2418. Schnupf, U.; Willett, J. L.; Bosma, W. B.; Momany, F. A. Carbohydr. Res. 2007, 342, 196–216. Biarnes, X.; Ardevol, A.; Planas, A.; Rovira, C.; Laio, A.; Parrinello, M. J. Am. Chem. Soc. 2007, 129, 10686–10693. Barnett, C. B.; Naidoo, K. J. J. Phys. Chem. 2008, B112, 15450–15459. Roslund, M. U.; Tahtinen, P.; Niemitz, M.; Sjoholm, R. Carbohydr. Res. 2008, 343, 101–114. Bosma, W. B.; Schnupf, U.; Willett, J. L.; Momany, F. A. J. Mol. Struct. (Theochem) 2009, 905, 59–69.

U. Schnupf et al. / Carbohydrate Research 345 (2010) 503–511 94. Csonka, G. I.; French, A. D.; Johnson, G. P.; Stortz, C. A. J. Chem. Theor. Comput. 2009, 5, 679–692. 95. Hansen, H. S.; Hunenberger, P. H. J. Comput. Chem. 2009. doi:10.1002/ jcc.21253. 96. Spiwok, V., Tvaroska, I., Carbohydr. Res. doi:10.1016/j.carres.2009.05.019. 97. Klamt, A.; Schuurmann, G. J. J. Chem. Soc., Perkin Trans. 2 1993, 799–805.

511

98. HyperChem 7.5, Hypercube, Inc., 115 NW 4th Street, Gainesville, FL 32601, USA. 99. Verlet, L. Phys. Rev. 1967, 159, 98–103; . Phys. Rev. 1968, 165, 201– 214. 100. PQS Ab Initio Program Package, Parallel Quantum Solutions, 2013 Green Acres,Suite E, Fayetteville, AR 72703, USA.