Molecular dynamics simulations of hexopyranose ring distortion in different force fields

Introduction

The flexibility of pyranose ring is a degree of freedom that often is neglected when studying the three-dimensional structure of carbohydrates. This is because in many cases the ring-inversion free energies are high enough to assume that the ring shape is ‘conformationally locked’ in the energetically favorable state (which, for hexopyranoses of the d series, is the ⁴C₁ conformation, as a rule). However, there is a growing number of reports pointing out that the ring shape and the related conformational rearrangements may play an important role in various, biologically-relevant processes [1], [2], [3], [4]. Furthermore, there exists a significant collection of various hexopyranoses for which a large degree of ring flexibility has been proven in terms of the NMR-based studies, performed for molecules in solution [5], [6], [7], [8], [9]. In particular, such studies usually concern idose [5], altrose [6] and their derivatives [8], [9] as well as some functionalized pyranoses such as mannuronic [7] or iduronic [8] acids. Even if the population of non-standard ring geometries is not high, it may significantly increase under influence of external factors. The example are atomic force microscopy (AFM) experiments [10], [11] during which the carbohydrate chain is pulled by the mechanical force. Contrary to the case of monomers, there are no experiental studies that would quantitatively describe the ring flexibility for residues in carbohydrate chain and the available data are either of qualitative nature [12], [13], [14] or concern the static crystal structures (see discussion in Ref. [15]). Currently, there is no quantitative interpretation of the AFM ‘stretching curves’ in terms of population of particular ring conformers in the pulled chain. Thus, the results of such experiments (although giving an important assessment of the favorable ring conformation), provide only indirect evidences which can hardly be translated into the quantitative ring-inversion properties [10], [11].

Because of scarce number of experimental reports that provide the quantitative description of hexopyranose ring flexibility, the theoretical approaches based on the molecular dynamics (MD) technique seem to be an attractive alternative to investigate this subject. The results of recent theoretical studies suggest that the ring flexibility may be enhanced when inserting a given residue in a chain of specific linkage type or upon functionalization (β-O₁-methylation; see Fig. 1a for atom numbering) [17]. Furthermore, the ring shape may influence the favorable orientation of glycosidic linkages while the significance of this effect depends not only on the ring-inversion properties but also on the linkage type [19], [20].

Fig. 1:

(a) Molecular structure of the d-glucopyranose (Glc) molecule (α- or β-anomer) illustrating the atom labeling used in the article. (b) Illustrative free energy (ΔF) profile (FEP) along the Cremer-Pople θ parameter [16] corresponding to the unfunctionalized α-Glc monomer within the 56A6_CARBO force field [17], [18]; the dashed vertical lines at θ=60° and 120° separate the three free-energy basins (⁴C₁, B/S, ¹C₄).

The accuracy of the MD-based methods relies primarily on the quality of potentials that approximate various interactions present in the molecular system of interest. In spite of several cases of successful application of the ab initio potentials into MD calculations aimed at studying the carbohydrate conformation [21], [22], the high computational cost makes such attempts limited to molecular systems of small size. Larger number of studies focused on puckering of the hexopyranose ring were conducted by using classical, biomolecular force fields dedicated to study the conformational properties of carbohydrates [15], [17], [18], [19], [20], [23], [24], [25], [26]. So far, the accuracy of such force fields has been comparatively studied only in the context of small sets of molecular systems (e.g. β-d-glucopyranose) [26] or only for unfunctionalized hexopyranose monomers [24]. The collection of all hexopyranoses of the d series either in monomeric forms or connected by α(1→4) and β(1→4) linkages has never been considered.

The main aims of our investigation are:

1. Describing and comparing the ring inversion properties of unfunctionalized monomers of all hexopyranoses of the d series calculated within four different force fields (GROMOS 56a6_CARBO/56a6_{CARBO_R} [17], [18], GROMOS 53a6_GLYC [24], CHARMM [27], [28] and GLYCAM06 [29]). Contrary to already existing studies (see e.g. table 4 in Ref. [24]) these properties will be expressed here as the two separate quantities: (i) the free energy change associated with the full chair-chair ring inversion from ⁴C₁ to ¹C₄ (ΔF_i); (ii) the free energy change associated with the ‘half’ ring inversion from ⁴C₁ to boat/skew-boat conformers (B/S) (ΔF_h).

2. Describing and comparing the ring inversion properties of residues in a chain composed of uniform monomers connected by glycosidic linkages of the α(1→4) and β(1→4) types calculated within three different force fields (GROMOS 56a6_CARBO/56a6_{CARBO_R}, GROMOS 53a6_GLYC and CHARMM). In this stage of the study the ring inversion properties are also investigated in the context of both the ΔF_i and ΔF_h values.

Methods

All MD simulations were carried out with the GROMACS 4.5.5 program [30], using the following four force fields: GROMOS 56a6_CARBO/56a6_{CARBO_R} [17], [18], GROMOS 53a6_GLYC [24], CHARMM [27], [28] and GLYCAM06 [29]. Note that one of these force fields exists in two versions, namely: GROMOS 56a6_CARBO and GROMOS 56a6_{CARBO_R}; the latter version contains parameters applicable for simulating the alkylated and glycosylated carbohydrates whereas parameters for unfunctionalized monomers are identical in both variants. For the sake of simplicity, in the present article we will be using only the 56a6_CARBO term. The GROMACS-compatible CHARMM parameters were downloaded from the official CHARMM website (http://mackerell.umaryland.edu/charmm_ff.shtml); the GROMOS 56a6_{CARBO_R} parameters were prepared during the procedure of the force field parameterization (see Ref. [17]); the GROMOS 53a6_GLYC parameters were downloaded from the website of Biomaterial Modeling Group of Federal University of Pernambuco (http://dqfnet.ufpe.br/biomat); the GLYCAM parameters for monomers were provided by Dr V.H. Rusu. In the case of trisaccharides, the pdb2gmx program was used to produce the molecular topologies from respective monomeric building blocks.

For each of these force field an appropriate water model was applied to describe the explicit solvent present in the system, i.e. SPC [31] (GROMOS force fields) or TIP3P [32] (GLYCAM06 and CHARMM). The systems considered consisted of one solute molecule solvated by water molecules within a cubic computational box simulated under periodic boundary conditions at constant volume and at a temperature of 298 K. Given the selected numbers of water molecules (~900 to ~1500 for the different systems), the box edges (3.1–3.6 nm for the different systems) were all preoptimized by a 1 ns constant-pressure MD equilibration at 1 bar and 298 K, ensuring an effective solvent density appropriate for these conditions in the subsequent isochoric production simulations.

The solutes considered were mono- or trisaccharides. The monosaccharides included the eight hexopyranoses of the d-series, namely allose (All), altrose (Alt), glucose (Glc), mannose (Man), gulose (Gul), idose (Ido), galactose (Gal) and talose (Tal), considering both the α- and β-anomers. The trisaccharides considered systematically were α(1→4)- and β(1→4)-linked homooligomers of the eight monosaccharides without any functionalization, the anomeric configuration at the reducing end being the same as that of the linkage.

The calculations involved the determination of one-dimensional free-energy profiles (FEPs) along the Cremer-Pople θ parameter [16] describing the conformation of specific rings. For the FEP analyses along θ, the intervals [0–60°], [60–120°] and [120–180°] were ascribed to regular chair (⁴C₁), boat/skew-boat (B/S) and inverted chair (¹C₄) conformations, respectively, as illustrated in Fig. 1. The free-energy changes (ΔF_i and ΔF_h) associated with the ring inversions (⁴C₁→¹C₄ and ⁴C₁→B/S, respectively) were defined as the difference between the values at the FEP minima in the two corresponding regions. These analyses were carried out systematically for the monomers and for the central residue in the trimers.

The calculation of FEPs relied on an enhanced-sampling scheme combining [33] parallel tempering [34] and well-tempered metadynamics [35] as implemented in the PLUMED 1.3 plug-in [36] of the GROMACS program [30]. The well-tempered metadynamics [35] relied on Gaussian local functions of widths 2.86°, an initial deposition rate of 0.01 kJ mol⁻¹ ps⁻¹ and a temperature parameter ΔT (see eq. (2) in Ref. [35]) of 7152 K (simulations within the GROMOS 53a6_GLYC force field) or 1788 K (remaining force fields). The parallel-tempering [34] relied on 16 metadynamics simulations carried out in parallel at different temperatures ranging from 298.0 to 363.2 K in steps of about 4.3 K, along with replica-exchange attempts performed at 2 ps intervals.

After equilibration, all simulations were carried out until the convergence of ΔF_i and ΔF_h values were detected which resulted in total simulation time varying from 10 ns (GROMOS 56a6_CARBO, GLYCAM06) to 110 ns (GROMOS 53a6_GLYC). The convergence of the ΔF_i and ΔF_h values was checked during the course of simulations by using the PLUMED built-in sum_hills tool and hand-written scripts.

In all simulations, the equations of motion were integrated using the leap-frog scheme [37] with a timestep of 2 fs. The full rigidity of the water molecules was enforced by application of the SETTLE procedure [38]. The temperature was maintained close to its reference value by application of the V-rescale thermostat [39] with a relaxation time of 0.1 ps. The translational center-of-mass motion was removed every timestep. The solute bond lengths were constrained by application of the LINCS procedure [40], [41] with a relative geometric tolerance of 10⁻⁴ to all covalent bonds (the GROMOS force fields) or to hydrogen-containing bonds (GLYCAM06 and CHARMM). In the case of GROMOS 56a6_CARBO, the non-bonded interactions were calculated using a twin-range scheme [42], with short- and long-range cutoff distances set to 0.8 and 1.4 nm, respectively, and an update frequency of five timesteps for the short-range pairlist and intermediate-range interactions. A reaction-field correction [43], [44] was applied to account for the mean effect of the electrostatic interactions beyond the long-range cutoff distance, using a relative dielectric permittivity of 61 [45]. In the case of GROMOS 53a6_GLYC, long-range electrostatic interactions were treated by the reaction-field method with a relative dielectric permittivity of 66. A 1.2 nm cutoff was used for the short range electrostatics and Lennard-Jones interactions. Within GLYCAM06-based simulations long-range electrostatic interactions were treated by using the particle mesh Ewald summation [46] (a 0.9 nm distance cutoff was used for direct space nonbonded calculation). In the case of the CHARMM force field, a force-switched smoothing function [47] was applied to Lennard-Jones interactions in the range of 1–1.2 nm. Long-range Coulomb interactions were handled using particle mesh Ewald with a real-space cutoff of 1.2 nm.

Results

The results of simulations are gathered in Tables 1 and 2 as the free energy values ΔF_i and ΔF_h. The correlation between selected sets of data is presented graphically in Figs. 2–5. All possible monomers of d-hexopyranoses as well as homotrimers composed of corresponding residues were considered, except of the case of GLYCAM06 where, due to the lack of available parameters, investigations were limited to non-Ido monomers.

Fig. 2:

Chair-chair-inversion (⁴C₁→¹C₄) free energies ΔF_i for the central residue of α(1→4)- and β(1→4)-linked homotrimers in water calculated using three different force fields, compared to corresponding values for monomers.

Fig. 3:

Regular chair-boat/skew-boat-inversion (⁴C₁→B/S) free energies ΔF_h for the central residue of α(1→4)- and β(1→4)-linked homotrimers in water calculated using three different force fields, compared to corresponding values for monomers.

Fig. 4:

Chair-chair-inversion (⁴C₁→¹C₄) free energies ΔF_i for the hexopyranose monomers in water calculated using three different force fields, compared to regular chair-boat/skew-boat-inversion (⁴C₁→B/S) free energies ΔF_h for the same compounds.

Fig. 5:

Chair-chair-inversion (⁴C₁→¹C₄) free energies ΔF_i for the central residue of α(1→4)- and β(1→4)-linked homotrimers in water calculated using three different force fields, compared to regular chair-boat/skew-boat-inversion (⁴C₁→B/S) free energies ΔF_h for the same compounds.

Let us mention that both force fields belonging to the GROMOS family (56a6_CARBO and 53a6_GLYC) have been parametrized with attention paid to possible deformations of pyranose rings. The respective papers (Refs. [17], [18], [24]) contain a detailed description of the hexopyranose ring-inversion properties which were studied in terms of all possible ring conformers [18] or just the chair-chair inversion [17], [24]. Initially, these analyses have been performed only at the level at which the ring flexibility-related features were adjusted, i.e. for unfunctionalized monomers; therefore, in the case of GROMOS 56a6_CARBO, the additional patch was needed to account for the more accurate description of the ring distortion in O₁-functionalized pyranoses [17]. Conversely, GROMOS 53a6_GLYC, GLYCAM06 and CHARMM have never been systematically tested with respect to the ring-puckering properties in the context of varying local environment of the distorted ring (i.e. unfunctionalized monomer vs. functionalized monomer vs. residue in a chain).

As the experimental data providing the desired quantities (i.e. the relative populations of at least the two most abundant ring conformers) are limited to idose [5] and altrose [6], for the rest of hexopyranoses the semi-empirical Angyal scheme [48] was used for comparative purposes (Table 1). Note that this scheme concerns only unfunctionalized hexopyranose monomers and the chair-chair inversion energy. Thus, for the second type of considered data, i.e. free energies of chair-boat/skew-boat transition, there is neither experimental data nor semi-empirical schemes that could be used to verify the results obtained by MD simulations within different force fields. The same is true for any type of ring inversions occurring for residues in a chain.

Both the GROMOS 56a6_CARBO and CHARMM force fields show an excellent agreement with both the Angyal scheme prediction and (where available) experimental data. The largest divergences are observed for compounds that are assumed to exhibit very rigid rings, independently of the method of calculations, e.g. β-Gal, β-Tal or α-Man. In the case of experimental data the largest divergence corresponds to the 56a6_CARBO force field and β-anomers of Ido and Alt and equals to ~6 kJ/mol. Note that for each compound the sign of the ring-inversion free energy remains uniform with that of the Angyal scheme and/or experimental data.

A similar level of qualitative agreement is offered by the GROMOS 53a6_GLYC force field. Also in this case the sign of ΔF_i values remains in agreement with the Angyal scheme prediction and the NMR data. However, now a much broader dispersion of the ring-inversion free energies is suggested. This results in both more extreme values of ΔF_i observed for most of compounds (|ΔF_i|>20 kJ/mol for 10 out of 16 compounds) and larger divergences between 53a6_GLYC and Angyal scheme which, in the extreme case of β-Glc, can reach ~50 kJ/mol. Moreover, in spite of the correct sign of ΔF_i for α-Ido, the absolute difference between the calculated and NMR-inferred values equals to ~25 kJ/mol. The related difference for β-Alt is also large (~12 kJ/mol).

The least degree of agreement between MD simulations and the available experimental data and/or the Angyal scheme was exhibited by the GLYCAM06 force field. GLYCAM06 is clearly at odd with the experimental data concerning α-Alt and suggests strongly dominating population of ¹C₄ conformers (ΔF_i<−25 kJ/mol). Furthermore, GLYCAM06 predicts the negative values of ΔF_i for α-All and α-Gul which is rather surprising in the context of both the semi-empirical scheme-inferred data [48] and the ab initio calculations [49]. The lack of parameters for Ido additionally hinders the comparison with the second set of existing experimental data.

The reasons for relatively worse GLYCAM06 performance in comparison with other force fields are most likely associated with the method of parametrization accepted during its development. Contrary to both GROMOS force fields and CHARMM, each molecule of unfunctionalized hexopyranose modeled by GLYCAM is associated with its own unique set of parameters (e.g. partial atomic charges) developed on the basis of ensemble-averaged conformation. As the dominating conformer is usually ⁴C₁, the derived charges are suited very well to describe the ensemble of states associated with the regular chair conformation but may fail to represent alternative ring shapes as well as their relative populations. GROMOS and CHARMM are designed to keep the set of parameters identical or nearly identical with respect to all hexopyranoses, thus, they suffer less from this effect.

Let us emphasize that the divergence between predictions of the given force field and the other types of data may influence the force field applicability only in some very special cases. This is due to the following reasons: (i) the behavior of the most frequently studied hexopyranose monomers (Glc, Man, Gal) is roughly independent of the applied force field and corresponds to strongly dominating population of the ⁴C₁ conformers; (ii) the free energy barriers separating particular ring conformers are high enough to prevent spontaneous ring ‘flips’ at the timescale of tens or even hundreds of ns of unbiased MD simulation.

Finally, let us refer to the ring-inversion free energies reported elsewhere [24]. The differences between our calculations and the data shown in table 4 of Ref. [24] are of quantitative nature and can be ascribed to different schemes of metadynamics-based calculations. Namely, calculations reported in Ref. [24] employed ‘plain’ metadynamics, without bias scaling, whereas in the present case, the well-tempered metadynamics scheme was used in combination with the replica-exchange method. Appreciating the existence of these minor differences, let us note that also in this case, the signs of ΔF_i values remain uniform.

The tested force fields, compared in the context of correlation between ring-inversion properties of monomers and those of residues in a chain, reveal a series of differences in their predictions. As discussed in details in Ref. [17], in the case of GROMOS 56A6_CARBO, the values of ΔF_i are tendentially reduced in magnitude upon insertion of the free residue into a chain (for 12 out of 18 trisachcarides). Furthermore, the correlation between the ΔF_i values for the monomers and for the corresponding residues in a chain is very poor along the d-hexopyranose series. In contrast, there exists a much higher correlation between the same quantities calculated within both CHARMM and GROMOS 53a6_GLYC force fields. While the monomer- and trimer-related ΔF_i values for β-anomers are very close to each other (except of decrease of ΔF_i by ~20 kJ/mol for β-Alt and β-Tal simulated within GROMOS 53a6_GLYC), the tendential decrease of the trimer-related ΔF_i values is observed for 4–6 out of 8 α-anomers (i.e. for Glc, Gal, Man and Tal in the case of CHARMM and Glc, Gal, Man, Tal, Gul and Alt in the case GROMOS 53a6_GLYC). Interestingly, Glc, Gal, Man and Tal are the same compounds for which also the GROMOS 56A6_CARBO force field predicts the largest decrease of ΔF_i upon inserting residue into a α(1→4)-linked chain. All of these compounds contain equatorial hydroxyl group at C₃ in the ⁴C₁ chair (see Fig. 1 for notation). Note that both the CHARMM and GROMOS 56A6_CARBO force fields agree in their predictions that the dominating ring conformation at the level of (non-Ido) monomer is the same as that of residue in a chain. On the contrary, in the case of GROMOS 53a6_GLYC, there are three compounds for which the monomer- and trimer-related ΔF_i values differ in their sings, i.e. α-Alt, α-Gul and α-Tal. The reason for that is a tendency to systematic decrease of ΔF_i exhibited by GROMOS 53a6_GLYC in the case of α(1→4)-linked residues. Such decrease in combination with the residues of relatively small, monomer-related, ΔF_i values results in the negative values of ring-inversion free energies for residues in a chain.

Analogous investigations focused on the possible correlation between the ΔF_h values calculated for the monomer and for the central residue in trimer give similar results, i.e. a poor correlation for GROMOS 56A6_CARBO and much better ones for CHARMM and GROMOS 53a6_GLYC. Interestingly, ΔF_h values characteristic of residues in a chain are again decreased in comparison to those of monomers for the same group of four α anomers containing equatorial hydroxyl group at C₃ in the ⁴C₁ chair conformer (the related decrease is additionally exhibited by α-Gul in the case of GROMOS 53a6_GLYC).

There is no quantitative experimental data that would allow to assign the ring-inversion free energy value to the residue being a part of a larger chain. The existing qualitative data determining the preferred conformation of such residues concern only selected systems and agree with the predictions of all tested force fields (e.g. Refs. [12], [13], [14] suggesting the ⁴C₁ chair conformation as dominating in residues of the amylose chain). Similarly, there is no semi-empirical schemes generalized for the case of glycosylated hexopyranoses and no ab initio methods have been used so far to systematically study this issue. For these reasons, currently there is no possibility to compare the quality of the force fields on the basis of the data collected in Table 2. However, the most striking difference between force fields is the presence or the absence of the correlation between ring-inversion properties of monomers and those of residues in a chain. This is especially important in the context of the lack of quantitative data concerning the ring conformation of residues in a chain and the possibility of predicting of such properties on the basis of relatively better recognized properties of respective monomers.

The mentioned difference can be explained when considering the monomer-inherent flexibility of the rings predicted by each of the tested force fields. The chair-chair ring-inversion free energies suggested by GROMOS 56A6_CARBO for monomers vary from ~−2 to ~28 kJ/mol (i.e. in the range of ~30 kJ/mol). The corresponding ranges for the CHARMM and GROMOS 53a6_GLYC force fields are larger and equal to ~40 kJ/mol and ~100 kJ/mol, respectively. Let us assume that the non-systematic alterations of the ring flexibility upon inserting residue into a chain, as suggested by GROMOS 56A6_CARBO, originate from the residue-residue interactions and that the strength of such interactions is of comparable magnitude for each of the force fields. Then, the ‘decorrelation effect’ will be more significant in the case of force field that predicts larger flexibility of rings, i.e. GROMOS 56A6_CARBO. The residue-residue interactions of the same strength are less important for rings of greater degree of rigidity, i.e. those described in terms of CHARMM and, especially, GROMOS 53a6_GLYC. Such hypothesis would link the possible decorrelation of ring-inversion free energies with the capability of the tested force fields to reflect very subtle differences between ring puckering at the level of monomers. As most of the currently available data allow only to assign a particular ring conformation as the dominating one and are not capable to characterize this domination by providing the ring-inversion free energy (or equivalent quantity), the possible alternative methods of choice would include MD-related techniques employing ab initio potentials of increased accuracy (see e.g. Ref. [21], [22]). However, due to the high significance of intramolecular hydrogen bonding [50] that strongly influences the conformation of hexopyranoses in the absence of solvent, such simulations would require the presence of explicit solvent, possibly modeled at the molecular-mechanics level of theory.

Interestingly, there exists a relatively good correlation of the chair-chair inversion free energies (ΔF_i) values with those of chair-boat/skew-boat ones (ΔF_h). The correlation is present for all tested force fields, however, it becomes very poor in the case of residue in a chain simulated within GROMOS 56A6_CARBO and is of (better) comparable quality for the remaining cases (see Figs. 4 and 5). In the case of monomers, with the exception of very rigid rings (ΔF_i>20 kJ/mol), the relation ΔF_i≤ΔF_h is always satisfied. For residues in a chain the same inequality holds except of α- and β-Ido simulated within the CHARMM force field. This means that the moderately flexible rings should be described in terms of the ⁴C₁ → ¹C₄ transition rather than the ⁴C₁ → B/S one. For more rigid rings the population of B/S conformers may be larger in comparison to ¹C₄, however both of them will be very scarce in relation to ⁴C₁. Finally, the most flexible rings (All and Ido, according to predictions of all force fields) may require a more complex description involving both types of transitions or, possibly, introducing additional coordinate(s) accounting for the ring shape within the B/S set of conformers (e.g. Cremer-Pople φ).

Conclusions

Four different classical force fields have been tested in the context of their capabilities to simulate the ring distortion in hexopyranose molecules (in monomers and in residues in a chain exploiting either α(1→4) or β(1→4) glycosidic linkages). Both the GROMOS 56a6_CARBO and CHARMM force fields offer very similar predictions for unfunctionalized monomers which, at the same time, remain in agreement with the available experimental data and with the semi-empirical Angyal scheme. GROMOS 56a3_GLYC also reflects this agreement at the qualitative level but differs quantitatively from GROMOS 56a6_CARBO and CHARMM by predicting a much broader dispersion of ring-inversion free energies varying significantly from one compound to another. GLYCAM06 exhibits the worst accuracy at the level of unfunctionalized monomers, being at odd with the NMR data for α-Alt and predicting the dominating population of inverted chairs for α-Gul and α-All. The most frequently studied hexopyranoses (Glc, Man, Gal) exhibit similar ring-inversion properties, independently of the force field used. All force fields predict alterations of the chair-chair inversion free energies upon inserting residue into a chain. However, the scale and direction of such alterations are dependent on the force field: GROMOS 56a3_GLYC and CHARMM suggest a strong correlation between monomer- and trimer-related properties and increased flexibility of some of the α(1→4)-linked residues, relative to the properties of corresponding monomers. According to GROMOS 56a6_CARBO, the ring-inversion free energies for residues in a chain are typically smaller in magnitude compared to those of the monomers (for both the α(1→4)- and β(1→4)- linkages) and correlate poorly with the latter. Except of very rigid rings, chair-boat/skew-boat inversion free energies are higher than those corresponding to the chair-chair inversion. For all tested force fields, there exists a correlation between chair-chair and chair-boat/skew-boat inversion free energy for both monomers and residues in a chain. However, it becomes rather poor for trimers simulated within GROMOS 56a6_CARBO.

Finally, as the availability of the reliable quantitative data related to the ring distortion of residues in a chain is very limited, the use of the force field-independent, ab initio potentials in combination with the presence of explicit solvent may be suggested as an alternative method to study this issue.