Reference interaction site model study on the anomeric equilibrium of D-glucose in aqueous solution

The anomeric equilibrium of D-glucose in aqueous solution was studied by the extended reference interaction site model (XRISM) theory. In this study, all of the rotational degrees of freedom were considered upon the exocyclic hydroxyl and hydroxymethyl groups, namely 729 stereoisomers for each anomer. The free energy


Introduction
Saccharides play a variety of important roles in a biological system, such as metabolism and biological recognition [1].The structural information on saccharides is the basis for understanding their functions.The most important experimental techniques to analyze the conformation of saccharides are X-ray crystal structure analysis and NMR [1].In particular, the latter has a great advantage in investigating the conformation of saccharides in solution.Theoretical and/or computational studies regarding conformations of saccharides have also been reported [2], where much attention is still focused on monosaccharides such as D-glucose [3][4][5][6][7], D-mannose [8] and D-xylose [9,10] as well as oligosaccharides [11,12].
The computational schemes used in studying the saccharides have been mainly quantum mechanical (QM) calculations [4,6,12] and molecular simulations [3,5,[7][8][9][10][11] that have included solvent effects to model solutions, which are of biological interest.In describing the solvent environment in the above methods, continuum solvent models [4,8] or explicit solvents [3,[5][6][7][9][10][11][12] have been employed.In a continuum solvent model, microscopic details of solvation structure are entirely neglected, and additional macroscopic parameters must be introduced to describe a solvent environment.In an explicit solvent model, on the other hand, numerical results are contaminated by statistical errors originating from a finite sampling number since solvent configurations must be statistically averaged.Since the energetic differences among the lowest and other local-minimum energy conformations of a saccharide molecule are generally very small [1], a computational analysis of the thermodynamic stability (or conformational stability) of saccharides requires both an accurate method and a consideration of many candidate conformations at the same time.In practice, this is a very demanding task, especially in solution.This fact makes it difficult to study saccharides based on computational chemical approaches at an atomistic level.As a consequence, even monosaccharides, which are the simplest in structure among saccharides, are still active objects to attract many researchers in the field of computational chemistry.
One of the promising methodologies of studying the conformations of saccharides in solution is the reference interaction site model (RISM) theory [13], which is an integral equation theory based on statistical mechanics.The RISM theory gives not only a solvation structure around solute T.Miyata molecules but also a solvation free energy, which is needed when discussing a thermodynamic stability.It should also be noted that the RISM theory is free from any statistical errors.Therefore, the RISM theory is expected to be advantageous in studying the thermodynamic stability of saccharides in solution.Up till now, no application of the RISM theory to the "true" saccharide molecules has been reported.The only example that has been reported so far is its application to 4,6-dimethyl-2-methoxytetrahydropyran, which was studied as a simple model material of monosaccharides by Maw et al. [14].
We studied the anomeric equilibrium of D-glucose in aqueous solution using the RISM theory, as a first step to investigate thermodynamic stabilities of saccharides.The anomeric equilibrium of D-glucose is one of the major subjects that is intensively investigated [3][4][5]15].The physical quantity that determines the anomeric equilibrium in solution is the free energy difference between solvated α-D-glucose and β-D-glucose.The simplest expectation based on the geometry is that the hydrogen-bonding capacity of the anomeric hydroxyl group of β-D-glucose with surrounding water is relatively high since it is in an equatorial form, which directs toward a relatively open space from the pyranose ring.On the other hand, the anomeric hydroxyl group of α-D-glucose takes an axial form, and consequently accesses of water molecules to the group may be relatively hindered by the rest of the pyranose ring compared to the β anomer [3,6,16].This paper discusses the free energy difference between α-D-glucose and β-D-glucose in aqueous solution calculated by the RISM theory, in comparison with the experimental data.Comparisons of our results with the other computational methods such as quantum mechanics and molecular simulations are also reported.

Free energy calculation
Figure 1 shows typical molecular shapes of (a) α-D-glucose and (b) β-D-glucose, respectively.In molecular simulations, free energy perturbation (FEP), thermodynamic integration (TI), locally enhanced sampling (LES), and so on have been used in order to calculate the free energy difference between α-D-glucose and β-D-glucose in solution [3,5].Alternatively, we evaluated the free energy difference ∆G α→β via the partition function of aqueous solution of each anomer at an infinite dilution limit.The formulae we used are as follows: where k B , T , Q, E C,i , ∆µ S,i represent the Boltzmann constant, temperature, partition function, conformational energy (i.e.intramolecular interaction energy of solute molecule), and solvation free energy, respectively.When evaluating the free energy difference in gas phase between both anomers, a similar strategy to the equations ( 1) and ( 2) was used in the past [5].In our method, the solvation free energy ∆µ S,i is calculated using the RISM theory (See the next section for the details of the RISM theory).The index, i, seen in the equation ( 2) identifies the conformation sampled discretely in the conformational space.Note that the summation in the equation ( 2) is taken over all the conformations sampled for the corresponding anomer.Moreover, E C,i + ∆µ S,i is a potential of mean force: it should be noted that the exponential function form that appeared in the equation ( 2) is obtained merely by an "incomplete" integration of the original configuration integral, where the integration should be performed only for the degrees of freedom on the solvent molecules (The derivation process can be found in the reference [13]).
In what follows, we describe the conformations considered in this study, which are related to the summation in the equation (2).The intramolecular degrees of freedom of D-glucopyranose include the forms of pyranose ring (i.e.puckering), vibrations on each bond length and on each angle, the rotational freedom around the C-C bond of the exocyclic hydroxymethyl group (hereafter, this bond is referred to as C5-C6 bond), and the rotational freedom of five exocyclic hydroxyl groups (hereafter, the bonds corresponding to this freedom are referred to as C-O bond).Concerning the forms of pyranose ring, the 4 C 1 chair structure is known to be the most stable among the forms of D-glucopyranose ring: for instance, it was reported that the free energy of the 4 C 1 chair structure of β-D-glucopyranose is lower than that of 1 C 4 chair form by as much as 8 kcal/mol, which was evaluated by a quantum mechanical calculation [17].In this work, we considered only the 4 C 1 chair form, and the other ring forms were not treated since they are expected to be insignificant in the analysis of thermodynamic stability of D-glucose.All the vibrational degrees of freedom were frozen in this study, related with the bond lengths and angles.We considered the remaining rotational degrees of freedom around the C5-C6 bond and five C-O bonds, which become six degrees of freedom in total.Since each rotational degree of freedom produces three stable rotamers, i.e. gauche-trans conformers, the final number of conformations we treated for each anomer is 3 6 = 729: that is, in the equation (2), the summation was taken over 729 conformers for each anomer.

RISM theory
In the evaluation of solvation free energy, an extended reference interaction site model (XRISM) integral equation was employed in this study.The details of the XRISM theory have been described elsewhere [13].Here, we show only basic equations.The XRISM integral equation for solute-solvent correlation functions reads where h, c, w are the total correlation function, the direct correlation function and the intramolecular correlation function, respectively.χ represents the pure solvent site density pair correlation function, and is evaluated from the dielectrically consistent RISM (or DRISM) theory (DRISM formalism can be found in [13] and references therein)."u" and "v" stand for solute and solvent, respectively.The asterisk, * , represents the convolution integral.The equation ( 3) can be solved with coupling the closure equations.We used the HNC closure and the KH closure, which are given by the equations ( 4) and ( 5), respectively [13]: namely, and Here, u uv αγ (r) is the interaction potential energy, and β is defined as β = 1/k B T .Hereafter, the coupling between the equations (3) and ( 4) is referred to as the XRISM/HNC theory, and that between (3) and ( 5) as the XRISM/KH theory, respectively.
The solvation free energy for the XRISM/HNC theory and that for the XRISM/KH theory are expressed by the equations ( 6) and (7), respectively [13]: and where ρ γ represents the solvent site density, and Θ(x) is the Heaviside step function.The solvation free energy is calculated for each conformer from equations ( 6) or (7), and substituted into the equation (2).

Computational details
Both the conformational energy and the solvation free energy were calculated for 729 conformers in each anomer, molecular structures of which can be searched by rotating C5-C6 bond and five C-O bonds one by one with the grid size of 120 • .Then, they were substituted into the equation ( 2).For the calculation of conformational energy of D-glucose, OPLS-AA parameters [18] were used.As to the solvent water, SPC/E model [19] was employed.The temperature was set as 298.15K.The number of grid points and the grid separation we used are 512 and 0.05 Å, respectively.The convergence of the XRISM calculation was accelerated by means of the MDIIS [13].In this study, both the XRISM/HNC theory and the XRISM/KH theory were numerically solved 1458 times, which corresponds to two anomers, and all calculations were converged without using the techniques such as "charge up" or "temperature lowering."

Results and discussion
Table 1 shows the free energy differences ∆G α→β between α-D-glucose and β-D-glucose in aqueous solution calculated on the basis of the XRISM/HNC and XRISM/KH theories.Here, ∆G α→β indicates the free energy change when an α-D-glucose is transformed into a β-D-glucose, as can be seen from the equation (1).The free energy difference in gas phase is also shown in table 1, which was evaluated by the same procedure as equations ( 1) and ( 2) except for the assumption of ∆µ S,i = 0 in the equation (2).Experimental data and the other computational results reported to date are also listed in table 1.The difference between the XRISM/HNC theory and the XRISM/KH theory was insignificant with regard to the free energy difference.The XRISM theory combined with the equations ( 1) and ( 2) predicts that the β anomer is more stable than α in aqueous solution.This result agrees with the experimental data qualitatively [15].However, the present method overestimated the magnitude of the free energy difference in aqueous solution (i.e. the relative stability of the β anomer) by about 1.3 kcal/mol, compared with the experimental data.In gas phase, the β anomer was predicted to be more stable than α, which is contrary to the quantum mechanical result by Barrows et al. [4] and the one of molecular dynamics simulation by Simmerling et al. [5].As an integral tendency, our method overestimated the relative stability of the β-D-glucose.We shall point out, nevertheless, that the difference between ∆G α→β in gas phase and that in solution is about 1.0 kcal/mol in our study, which is close to the results reported by Barrows et al. (i.e.0.6 kcal/mol) [4] and Simmerling et al. (i.e.0.98 kcal/mol) [5].This fact implies that the solvation effect itself can be treated with a satisfactory accuracy in our calculation.Therefore, the above overestimation observed in our study may be attributed to an incomplete sampling in the conformational space, e.g. a lack of vibrational flexibilities of D-glucose.There remain the possibilities that a structural optimization related to the bond lengths and bond angles as well as dihedral angles may improve the numerical accuracy in describing the free energy surfaces.
Here, we show each component of potential of mean force (i.e.energy decomposition) such as van der Waals (i.e.vdw or Lennard-Jones) term, electrostatic term and so on, in order to explore physical origins of the fact that the β anomer is more favored in aqueous solution than α.Table 2 lists the conformational energy E C,α→β , intramolecular vdw energy E vdw,α→β , intramolecular electrostatic energy E elec,α→β , torsion energy E tors,α→β , and solvation free energy ∆µ S,α→β , where the numerical values shown in the table were evaluated as averages over all conformations using a Boltzmann weight factor, and presented in the form of an energy gain when an α-D-glucose was transformed into a β-D-glucose.In table 2, we observe again an insignificant difference between the XRISM/HNC theory and the XRISM/KH theory.It is found that the principal factor to stabilize the β anomer in aqueous solution is the solvation free energy, which stabilizes β by as much as 2.5 kcal/mol better than α in solution.It is also interesting to note that in terms of the intramolecular electrostatic energy, a slightly higher stability of β-D-glucose than α-D-glucose in gas phase (ca -0.4 kcal/mol) is drastically reversed when the surrounding environment is changed into aqueous solution (ca 1.3 kcal/mol): that is, the intramolecular electrostatic contribution favors the α anomer better than β in aqueous solution.On the one hand, both anomers are stabilized by forming intramolecular hydrogen bonds in gas phase.On the other hand, in aqueous solution, a competition can occur between the neighboring water molecules and intramolecular hydroxyl groups when a D-glucose molecule forms hydrogen bonds.The behavior of the intramolecular electrostatic contribution observed in table 2 may be interpreted as the evidence that the β anomer interacts with water molecules through hydrogen bonds more preferably than α in solution, which is regarded as a

Conclusion
The anomeric equilibrium of D-glucose was studied using the XRISM/HNC and XRISM/KH theories.All 729 possible torsional stereoisomers of the 4 C 1 chair form were taken into consideration in this study for each anomer.The free energy difference between the α and β anomers was calculated from partition functions, and a relative stability of each anomer was discussed on the basis of the free energy difference.Both the XRISM/HNC and XRISM/KH theories predicted that β-D-glucose was more stable in aqueous solution than α-D-glucose, which agreed with an experimental result qualitatively.It was found that the principal factor to stabilize the β anomer in solution compared with the α was the solvation free energy, which was revealed by the analysis of each component of the potential of mean force.Correspondingly, the intramolecular electrostatic energy of the β anomer was remarkably higher than that of α in solution.These findings imply that in aqueous solution, the β anomer favors an interaction with water molecules through hydrogen bonds, compared to α.Although our method succeeded in reproducing experimental data in a qualitative sense, it overestimated the relative stability of the β anomer significantly.This is still left as an open problem to be overcome, probably by introducing a better sampling algorithm in a conformational space.

Figure 1 .
Figure 1.Molecular shapes of (a) α-D-glucose and (b) β-D-glucose, respectively.Black, gray, and white spheres represent oxygen, carbon, and hydrogen atoms, respectively.The forms of pyranose ring depicted in the figures are 4 C1 chair structures (see the text).

Table 1 .
Free energy difference [kcal/mol] between α-D-glucose and β-D-glucose.The values are in the form of a free energy change when an α-D-glucose is transformed into a β-D-glucose.

Table 2 .
[3]ferences in averaged energy components [kcal/mol] between α-D-glucose and β-Dglucose.The values are in the form of a corresponding energy change when an α-D-glucose is transformed into a β-D-glucose.E C,α→β E vdw,α→β E elec,α→β E tors,α→β ∆µ S,α→β This interpretation is consistent with the above-mentioned fact that the solvation free energy stabilizes the β anomer better than α.Also, our results and interpretation qualitatively agree with the data of each energy component reported by Ha et al.[3].