Subscriber access provided by Iowa State University | Library
Molecular Mechanics
A GROMOS Force Field for Furanose-Based Carbohydrates Karina Nester, Karolina Gaweda, and Wojciech Plazinski J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00838 • Publication Date (Web): 04 Jan 2019 Downloaded from http://pubs.acs.org on January 4, 2019
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
A GROMOS Force Field for Furanose-Based Carbohydrates Karina Nester, Karolina Gaweda, Wojciech Plazinski1 Jerzy Haber Institute of Catalysis and Surface Chemistry, Polish Academy of Sciences, Niezapominajek Str., 8, 30-239 Cracow, Poland
Abstract The article describes a GROMOS force field parameter set for molecular dynamics simulations of furanose carbohydrates. The proposed united-atom force field is designed and validated with respect to the conformational properties of furanose mono-, di-, oligo- and polymers in aqueous solvent. The set accounts for the possibility of arbitrary glycosidic linkage connectivity between units, O-alkylation as well as of different anomery. The compatibility with the already existing, pyranose-dedicated GROMOS 56A6CARBO/CARBO_R set allows to use the presently proposed extension for studying more diverse and biologicallyrelevant carbohydrates that exploit both pyranose and furanose units. The validation performed against the quantum-mechanical and experimental data concerning the structural and conformational features shows that the newly-developed set is capable to reproduce conformational equilibrium within the furanose ring, relative free energies of anomers, hydroxymethyl rotamers, and glycosidic linkage conformers. Additionally, the results concerning the conformation of the furanose ring with relation to the two-state model as well as other conformational features of furanose-containing saccharides are discussed.
1Corresponding
author. E-mail address:
[email protected] ACS Paragon Plus Environment
1
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 2 of 47
1. Introduction Carbohydrates exploiting furanose units play an important role in biology with a special emphasis put on their structural role as a building blocks of nucleic acids.1–3 However, they are also central for other biology-related fields, such as nutrition, medicine, infection, and immune defense.4–8 The most common furanose-containing saccharide, sucrose, is often treated as a model compound to study the conformational preferences of carbohydrates and for refining the experimental and computational methods toward that aim.9–11 Other, popular examples of biologically-relevant furanoses are: ribose and deoxyribose (components of RNA and DNA, respectively), raffinose, lactulose, melezitose, 1-kestose, 6-kestose, isomaltulose, planteose, nystose and arabinoxylan. There also exists a class of polysaccharides called fructans,12 that almost exclusively exploit the furanose-furanose glycosidic linkages. Atomistic molecular dynamics (MD) simulations are a powerful tool that can be used to provide insight into molecular details of the carbohydrate structure, conformation, as well as of intra- and intermolecular interactions that are not accessible for experiment. The underlying molecular model determines the accuracy of the MD simulations and the quality of the final results. Thus, carefully designed and properly validated force fields are of a great use for biomolecular community. Recent attempts in the development of the carbohydrate-dedicated force fields for explicit-water MD simulations were focused mostly on pyranoses and include the GROMOS family force fields for unfunctionalized aldohexopyranoses13,14 and pyranoses in general15,16 as well as the extension for uronic acids.17 The other types of force field families, such as CHARMM,18,19 OPLS20 and GLYCAM21 also provide parameter sets for simulation of pyranose-based carbohydrates. Considering furanoses, there exists a much smaller number of available fixed-charge force fields. Namely, the CHARMM sets for furanose monomers22 and furanose-containing di- and oligomers23 as well as the GLYCAM set for aldopentofuranoses24 were proposed. Both these sets have their shortcomings. More precisely, the parameterization of the furanose ring flexibility in the CHARMM force field was performed against the PSEUDOROT-originating25 secondary data concerning the conformational populations in the north and south regions of the pseudorotation itinerary1,26 (see Fig. 1). Such procedure assumes the full applicability of the so-called two-state model for pseudorotation;24,26 however, there is a growing number of evidences that this model is often oversimplified.27–30 The GLYCAM furanose-dedicated force field24 is focused on only single group of furanoses,
ACS Paragon Plus Environment
2
Page 3 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
i.e. aldopentoses, and the procedure of its generalization for other types of furanosides is not straightforward. It is also worth noting that very recently the CHARMM polarizable force field for saccharides (including furanoses) has been proposed.31 The main aim of this article is to present a GROMOS force field that would be free of the above-mentioned shortcomings and, additionally, compatible with the existing, pyranosededicated GROMOS set called 56A6CARBO15 as well as with its revised version, 56A6CARBO_R.16 Such combination of the pyranose- and furanose-related parameters permits to perform computational studies concerning a wide group of biologically-relevant and structurally-heterogeneous di-, oligo- and polysaccharides that contain both furanose and pyranose building blocks. The presently proposed force field is designed with respect to conformational properties of furanoses in aqueous solution (conformational equilibrium within the ring, relative free energies of anomers, hydroxymethyl rotamers, and glycosidic linkage conformers). The force field build-up rules (in their large part shared with the 56A6CARBO set) allow for constructing molecular models for diverse furanoses, structurally exceeding beyond the set of compounds used by us in the validation procedure.
Fig. 1. (A) The pseudorotation itinerary exploiting the Altona-Sundaralingam coordinates1 and describing the possible conformations of the five-membered rings. (B) Atom numbering shown on the example of -D-fructofuranose and the chemical structures of disaccharides used in the MD simulations; the definitions of the glycosidic dihedral angles are also shown.
ACS Paragon Plus Environment
3
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 4 of 47
The article is organized as follows: the “Methods” section contains: the nomenclature used in this work (section 2.1), the brief description of the force field functional form (section 2.2), the force field build-up rules and the final parameters, including the description of both modified and unchanged elements (section 2.3), the parameterization stages, including the description of model fitting procedures (section 2.4) and the detailed description of the computational methodology used throughout the study (section 2.5). The “Results and Discussion” section, divided into several subsections with respect to the considered aspect, reports the results of force field validation procedure (including the comparison with the available experimental data) as well as the discussion on the force field applicability and its possible limitations. This section also contains the description of several findings regarding the conformational features of furanose-containing saccharides. Finally, the “Conclusion” section summarizes the main findings and provides concluding remarks.
2. Methods 2.1 Nomenclature and Definitions The atom numbering is in accordance with the IUPAC nomenclature and is presented on the example of D-fructofuranose in Fig. 1. Furthermore, Fig. 1 contains the graphical illustration of the exocyclic dihedral angles describing the conformation of glycosidic linkages (i.e. , and ). Note that atom numbering as well as the corresponding dihedral angle definitions may vary from one compound to another [e.g. the anomeric carbon atom can be either C1 (for aldoses) or C2 (for ketoses)]. The ring geometry was described in terms of the AltonaSundaralingam notation1 which exploits the two coordinates: P (corresponding to the pseudorotation phase) and m (corresponding to the pseudorotation amplitude). The alternative Cremer-Pople coordinates32 were used in the case of oxolane, in order to subsequently apply the Wu-Cremer variation of the Karplus equation.33 The most significant dihedral angles subjected to deeper, quantitative analysis included: the dihedral angle, describing the rotation of the exocyclic hydroxymethyl group (defined by the O6−C6−C5−O5 and O5−C5−C4−O4 quadruplets of atoms for ketohexofuranoses and aldopentofuranoses, respectively), , and dihedral angles, describing the orientation of the glycosidic linkage (defined for (2→1)-linked di--D-fructofuranose by the O5-C2-O2C1’, C2-O2-C1’-C2’and O2-C1’-C2’-O5’quadruplets of atoms, respectively, whereas for (1→5)linked di--D-arabinofuranose by the O4-C1-O1-C5’, C1-O1-C5’-C4’and O1-C5’-C4’-O4’ quadruplets of atoms, respectively) and the dihedral angle, describing the rotation of lactol
ACS Paragon Plus Environment
4
Page 5 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
group (defined for unfunctionalized and O1-methylated aldopentofuranoses by the O4−C1−O1−H1 and O4−C1−O1−C quadruplets of atoms, respectively). 2.2 Force field functional form The functional form of the GROMOS force field can be divided into non-bonded and bonded interactions34. The bonded interactions are the sum of bond, bond angle, improper dihedral angle, and torsional dihedral angle terms. The nonbonded interactions are the sum of van der Waals and electrostatic (Coulomb) interactions between all pairs of atoms. All individual terms are described in detail elsewhere34. Here we limit ourselves to provide the basic functional form of the GROMOS force field which is applied in the presently proposed extension. The electrostatic interactions between atom pairs consist of three contributions. The first is a sum over all interacting (i,j) pairs, using a Coulomb potential which depends on both the partial charges q on the atoms and on the atom positions (r): 𝑞𝑖𝑞𝑗
𝑉C(𝐫) = ∑pairs 𝑖.𝑗4𝜋𝜀0𝜀1𝑟𝑖𝑗 + 𝑉RF(𝐫) + 𝑉RFc,
(1)
where rij is the distance, 0 is the dielectric permittivity of vacuum and 1 the relative permittivity of the medium in which the atoms are embedded. Additional, cutoff-dependent terms, VRF(r) and VRFc, represent the reaction-field contribution and a distance-independent reaction-field term, respectively. More details on these terms can be found in ref. 34. The nonbonded van der Waals interactions are calculated as a sum over all interacting nonbonded atom pairs (i,j) using a Lennard–Jones (LJ) 12/6 interaction function with parameters C12 and C6: 𝑉LJ(𝐫) = ∑pairs 𝑖.𝑗
(
C12,𝑖𝑗 𝑟𝑖𝑗
―
C6,𝑖𝑗 𝑟𝑖𝑗
),
(2)
The parameters C12,ij and C6,ij depend on the type of atoms involved and the character of the interaction. In particular, the special parameters are applied to atoms that are separated by three covalent bonds (so-called third neighbors and the corresponding 1-4 LJ interactions). The functional form of the potential-energy term, associated with the stretching of bond and applied to all unique pairs of covalently linked atoms is given by eq. (3): 1
2
𝑉(𝑏) = 4𝑘𝑏(𝑏2 ― 𝑏20) ,
(3)
where b is the bond-length distance, b0 its reference value, and kb the corresponding force constant. The functional form of the potential-energy term, associated with the bending of bond angle and applied to all unique triplets of covalently linked atoms, is given by eq. (4):
ACS Paragon Plus Environment
5
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 6 of 47
1
𝑉(𝜃) = 2𝑘𝜃(cos 𝜃 ― cos 𝜃0)2,
(4)
where θ is the bond-angle value, θ0 its reference value, and kθ the corresponding force constant. The functional form of the potential-energy term, associated with the deformation of improper-dihedral angle and applied to the subset of atom quadruplets in order to control outof-plane or out-of-tetrahedron distortions, is given by eq. (5): 1
𝑉(𝜉) = 2𝑘𝜉(𝜉 ― 𝜉0)2,
(5)
where ξ is the improper-dihedral angle value, ξ0 its reference value, and kξ the corresponding force constant. Finally, the functional form of the potential energy term, associated with the torsion around dihedral angle is given by eq. (6): 𝑉(𝜙) = 𝑘𝜙[1 + cos(𝑚𝜙 ― 𝜙0)],
(6)
where ϕ is the dihedral angle value, m the multiplicity of the term, 0 the associated phase shift, and kϕ the corresponding force constant. This term is applied to a subset of dihedral angles as specified in section 2.3. The presented extension relies mainly on the modifications within the terms accounting for the deformation of torsional dihedral angles (eq. (6)) and for the non-bonded interactions expressed by the LJ potentials (eq. (2)). The parameters present within the remaining terms (bond-stretching, bond-angle bending, improper dihedral deformation and Coulombic interactions) are unchanged with respect to the pyranose-dedicated set15; however, the additional tests have been performed in order to check the correctness of the previously accepted atomic partial charges. Section 2.3 contains the detailed description of both newlyintroduced and unmodified elements of the force field. 2.3 Force field build-up rules The force field build-up rules given below are restricted to compounds that contain furanose rings and: (i) involve only the C, O and H elements; (ii) involve only the single (saturated) CO, C-C, O-H, and C-H (treated implicitly according to the united-atoms approach) bonds; (iii) do not involve oxygen functions other than alcohol, ether, hemiacetal or acetal; (iv) may exhibit structural topologies consisting of furanose rings functionalized by arbitrary patterns of exocyclic groups, as well as any furanose-furanose or furanose-pyranose rings separated by at least one intermediate atom. The latter limitation relies on the compatibility of the presently proposed set with the GROMOS 56A6CARBO set15 and the use of the corresponding, pyranosededicated parameters. Furthermore, the above restrictions closely resemble those inherent to 56A6CARBO, in order to remain the compatibility of these two parameter sets.
ACS Paragon Plus Environment
6
Page 7 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Tab. 1. Parameters for bond stretching, bond-angle bending and improper-dihedral distortion in the presently proposed force fielda. Bond stretching O‒H C‒O C‒C Bond-angle bending C‒C‒C C‒C‒O, O‒C‒O C‒O‒C C‒O‒H Improper-dihedral distortion CH1 center a
kb [kJ∙mol-1∙nm-4] 15.70 ∙ 106 6.10 ∙ 106 5.43 ∙ 106 kθ [kJ∙mol-1] 285 320 380 450 kξ [kJ∙mol-1∙deg-2] 0.102
bo [nm] 0.1000 0.1435 0.1520 θ0 [deg] 109.5 109.5 109.5 109.5 ξ0 [deg] 35.2644
The parameters refer to the standard GROMOS functional forms for bond stretching, bond-angle bending and
improper-dihedral distortion (see eqs. (3), (4) and (5)). These terms are included for all unique atom sequences matching the indicated topological pattern (bond, bond angle), or only once for each chiral center (improper dihedral; atom sequence determined by the chirality of the center). The following notation applies: any explicit hydrogen atom (H); any oxygen atom (O); any carbon atom (C); trisubstituted united carbon atom (CH1).
ACS Paragon Plus Environment
7
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 8 of 47
Tab. 2. Atom types, atomic partial charges, and charge group definitions for the presently proposed set (in agreement with those present in 56A6CARBO) shown on the example of both unfunctionalized and O2-methylated ketohexofuranose.a Atom Charge group Unfunctionalized C1 1 O1 1 H1 1 C2 2 O2 2 H2 2 C3 3 O3 3 H3 3 C4 4 O4 4 H4 4 C5 2 O5 2 C6 5 O6 5 H6 5 Alkylated Oalk (b) Calk (c) Cchn (d) Ogly (d)
Charge
Atom type
0.232 -0.642 0.410 0.464 -0.642 0.410 0.232 -0.642 0.410 0.232 -0.642 0.410 0.232 -0.464 0.232 -0.642 0.410
CH2 OA HO CH0r OA HO CH1r OA HO CH1r OA HO CH1r Or CH2 OA HO
-0.464 0.232 0.000 -0.464
OE CH3 CH* OE
a
Atom numbering in agreement with Fig. 1.
b
Oalk and Calk are used for an alkylated sugar to replace On and Hn (n = 1, 2, 3, 4, or 6); these atoms are included
into charge group n, which remains neutral. c
Cchn is used for any united carbon atom of an alkyl chain that is not bound to an oxygen atom; CH* stands for
CH0, CH1, CH2, or CH3 depending on the chain connectivity at this atom; each atom of this type belongs to its own (single-atom) neutral charge group. d
Ogly is used for a non-reducing residue in an oligosaccharide to replace an anomeric hydroxyl group. For
glycosylation at the exocyclic group (e.g. levan or inulin) this atom is included together with the functionalized carbon atom C'n, the ring oxygen atom and atoms neighboring ring oxygen atoms into single, neutral charge group. For glycosylation at the furanose ring (e.g. sucrose), the group additionally contains the ring oxygen atom and atoms neighboring ring oxygen from the second residue.
ACS Paragon Plus Environment
8
Page 9 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Tab. 3. Parameters for torsional interactions in the presently proposed force field. The examples of topological patterns are given for both unfunctionalized and O2-methylated -D-fructofuranose
Purpose
3
0
Generic C-C torsion
4.9
3
0
O5-C2-O2-C7
3.77
3
0
Generic C-O torsion C-O torsion for exocyclic groups
C2-C1-O1-H1, C5-C6-O6-H6 O5-C2-O2-H2, C2-C3-O3-H3, C3-C4-O4-H4
C2-C1-O1-H1, C5-C6-O6-H6
2.4
3
0
Generic hydroxyl torsion
C2-C3-O3-H3, C3-C4-O4-H4
1.86
3
0
Hydroxyl torsion exocyclic groups
Or‒CH*r‒O‒H (c)
O5-C2-O2-H2
-
1.00
1
180
T4F2
Or‒CH*r‒O‒C (c,e)
-
O5-C2-O2-C7
2.00
1
180
T5
CH*r‒Or‒CH*r‒O (a)
C5-O5-C2-O2
-
2.00
1
180
T5R1
CH*r‒Or‒CH*r‒O (b)
-
C5-O5-C2-O2
4.50
1
180
T5R2
CH*r‒Or‒CH*r‒O (b)
-
-
6.50
1
180
T6 T7 T8 T9
Xr‒CH*r‒O‒Cx (a,e,g) Xr‒CH*r‒O‒Cx (a,e,h) Xr‒CH*r‒CH*‒Cx (a,e,i) Xr‒CH*r‒CH*‒Cx (a,e,j)
k [kJ/mol]
m
5.92
X-C-O-X, X ≠ H (a,d,e)
O5-C2-C1-O1, O5-C5-C6-O6 -
Me---Dfructose O5-C2-C1-O1, O5-C5-C6-O6 -
T2F
Xr-CH*r‒O-C (c,d,e)
-
T3
X‒C‒O‒H (a,d,f)
T3F
Xr‒CH*r‒O‒H (c,d)
T4F1
Type
Pattern
T1
X-C-C-X (a,d,e)
T2
-D-fructose
[deg]
for
Exo-anomeric torsion for hydroxyl group Exo-anomeric torsion for methoxy group Anomeric torsion for unfunctionalized residue Anomeric torsion for methylated residue Anomeric torsion for glycosylated residue
C3-C2-O2-C7 2.50 1 60 Exocyclic methoxy torsion O5-C2-O2-C7 2.50 1 -60 1.50 1 60 Exocyclic ethyl torsion 1.50 1 -60 C4-C5-C6-O6, C4-C5-C6-O6, T10 Xr‒CH*r‒CH*‒O (a,e,i) 1.00 1 60 C3-C2-C1-O1 C3-C2-C1-O1 Exocyclic oxymethoxyl torsion O5-C5-C6-O6, O5-C5-C6-O6, T11 Xr‒CH*r‒CH*‒O (a,e,j) 1.00 1 -60 O5-C2-C1-O1 O5-C2-C1-O1 O3-C3-C4-O4, O3-C3-C4-O4, O2-C2-C3-O3, O2-C2-C3-O3, Ox‒X‒X‒Ox; Oxygen-oxygen gauche T12 4.50 1 180 O6-C6-C5-O5, O6-C6-C5-O5, not T13 (a,k) torsion O1-C1-C2-O5, O1-C1-C2-O5, O1-C1-C2-O2 O1-C1-C2-O2 O5-C5-C4-O4, O5-C5-C4-O4, Oxygen-oxygen intracyclic T13 Ox‒Xr‒Xr‒Or (a) 1.00 1 180 O5-C2-C3-O3 O5-C2-C3-O3 torsion C2-C3-C4-C5, C2-C3-C4-C5, TFR1 CH*r-Xr-CH*r-CH*r (c) C2-O5-C5-C4, C2-O5-C5-C4, 4.00 5 0 Intrinsic conformational C5-O5-C2-C3 C5-O5-C2-C3 properties of unsubstituted ring O5-C5-C4-C3, O5-C5-C4-C3, TFR2 Or-CH*r-CH*r-CH*r (c) 2.90 5 0 O5-C2-C3-C4 O5-C2-C3-C4 a present in the 56A6 b CARBO set; introduced in the 56A6CARBO_R set; c introduced in the current set; d one per bond; e not applicable for interactions within the ring; f not applicable for exocyclic hydroxyl groups; g to be used if the improper dihedral angle O-Xr-Xr'-CH*r is negative in the reference structure, where Xr' is the other ring neighbor of CH*r; h to be used if the improper dihedral angle O-Xr-Xr'-CH*r is positive in the reference structure, where Xr' is the other ring neighbor of CH*r; i to be used if the improper dihedral angle CH*-Xr-Xr'-CH*r is negative in the reference structure, where Xr' is the other ring neighbor of CH*r; j to be used if the improper dihedral angle CH*-Xr-Xr'-CH*r is positive in the reference structure, where Xr' is the other ring neighbor of CH*r; k to be applied only if the atom sequence does not match the topological pattern for T13.
ACS Paragon Plus Environment
9
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 10 of 47
The following build-up rules were used in the presently proposed parameter set: A. The covalent bond stretching, bond-angle bending, and improper-dihedral distortion parameters are identical to those used in the previous GROMOS sets of parameters (i.e. 43A1, 45A3, 45A4, 53A6, 56A6CARBO, 56A6CARBO_R and 54A7).13,15,34 The functional forms of the corresponding potentials are given by eqs. (3), (4) and (5), respectively. The values of the associated parameters are collected in Table 1. They were obtained on the basis of the experimental spectroscopic (force constants) and X-ray diffraction (optimal distance and angle values) data for small molecules. B. The atomic partial charges (being the parameters for potential given by eq. (1)) are derived by means of a bond-increment approach.35,36 The details of this method in the context of carbohydrates are given in ref.15 Increments are used for covalent bonds involving an oxygen atom and are described by the parameters qC→O (charge transfer from a neighboring carbon atom) and qH→O (charge transfer from a neighboring hydrogen atom). We accepted the qC→O = −0.232 and qH→O = −0.410 e values, identical to those introduced in 56A6CARBO. Such choice facilitates the combination of the presently proposed force field with the 56A6CARBO set, aimed at description of the more complex saccharides, involving both the furanose and pyranose residues. Accepting these values is also justified on the basis of the RESP charge fitting following the electrostatic potential calculations performed at the HF/631G* level of theory. The corresponding calculations involved six different model compounds (shown in Fig. 2) and the largest deviation between the RESP-fitted atomic charges and their bond-increment-derived analogues does not exceed 0.09 e. The charge group definitions (inherent to the GROMOS force field) are identical to those used in the 56A6CARBO force field (see ref.15 and Fig. 4 therein). Namely, they are defined so as to include the first covalent neighbors of all oxygen atoms it encompasses, whereas other (neutral) atoms define their own charge groups. As a result, all charge groups are neutral. Tab. 2 includes the definition of the charge groups for a specific case of both unfunctionalized and O2-methylated ketohexofuranose, whereas Fig. 2 shows the charge groups for model compounds, used in the RESP calculations. The graphical illustration of the bond-increment method is given in Fig. 4 in ref.15 The resulting set of possible atomic charges in the most common furanosides is rather narrow and limited to only several possible values for each element.
ACS Paragon Plus Environment
10
Page 11 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Fig. 2. The partial atomic charges derived from the RESP calculations (red numbers) or from the bond-increment approach (black numbers) for six compounds containing furanose rings. The dotted lines limit the neutral charge groups. The details are given in the text.
C. The Lennard–Jones interaction parameters for atoms that are not covalently linked or are separated by more than three covalent bonds are unchanged with respect to those existing in the 53A6, 56A6CARBO and 56A6CARBO_R force fields. See Tab.7 in ref.34 and Tab. 4 in ref.15 for their values. The latter reference contains the saccharide-specific atom types (e.g. Or, CH0r, CH1r) that do not appear in the earlier versions of the GROMOS force fields as well as the slight modifications of some of the LJ existing parameters. Note that the Xr atom type (where X = CH0, CH1, CH2, O) corresponds to the atoms that create a furanose ring.
ACS Paragon Plus Environment
11
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 12 of 47
The assignment of these atom types into the atoms of a specific molecule is exemplary shown in Tab. 2 for both unfunctionalized and O2-methylated ketohexofuranose. Eq. (2) describes the associated potential of interactions. D. Following the GROMOS conventions,34 first and second covalent neighbors are excluded from any nonbonded (Coulomb and Lennard–Jones) interaction, whereas third covalent neighbors interact with normal electrostatic interactions but LJ interactions determined by a special set of parameters (see Tab. 9 in ref.34 and Tab. 4 in ref.15). The GROMOS 56A6CARBO/CARBO_R force field introduced a number of exceptions which do not obey the above rules. These exceptions concern specific topological patterns that can be found in pyranose molecules (see Tab. 5 in ref.15). The presently proposed force field also contains several exceptional LJ interactions between selected atom pairs. More specifically, the N1F term has been introduced to describe the third-neighbor lactol repulsion between lactol hydrogen (H) and ring oxygen (Or) atoms. This term is analogous to the N1 one, present in the 56A6CARBO/CARBO_R set and applies to the Or-CH*r-O-H topological pattern; however, presently the C12 LJ term has been reduced to 0.25 ∙ 10-6 kJ∙nm12/mol, whereas the C6 term is equal to 0. Thus, in the case of furanose-inherent interactions, N1 is replaced by N1F. The second topological pattern ascribed to the exceptional LJ interactions is CH*-Xr-Xr-CH*r (N3, present also in the same form in the 56A6CARBO/CARBO_R set) described by the following values of the LJ parameters: C6 = 5.689469 ∙ 10-3 kJ∙nm6/mol and C12 = 3.33986∙ 10-6 kJ∙nm12/mol. Regarding the remaining 56A6CARBO-related exceptional LJ interactions (i.e. N2, N4 and N5; see Tab. 5 in ref.15), they were removed from the present, furanose-dedicated set and replaced by the corresponding, standard LJ interactions for the respective atom pairs (which may include the parameters for the 1-4 pairs). Note that in the case of heterogeneous compounds containing both furanose and pyranose units, the exceptional LJ interactions at the furanose/pyranose interface should be assigned separately for each unit, i.e. depending to which unit the corresponding topological pattern belongs. E. The covalent torsional dihedral parameters, listed in Tab. 3, are reoptimized and may differ from those characteristic of the 53A6 and 56A6CARBO/CARBO_R sets, especially when considering the terms involving ring atoms. However, some terms (e.g. those describing the rotation of glycosidic linkages as well as of the hydroxymethyl group) are unchanged with respect to 56A6CARBO. The functional form of the potential used is given by eq. (6). The significant part of the parameters in the present set is unchanged in comparison to the 56A6CARBO/CARBO_R set (T1-T3, T5-T13). This includes the terms describing the rotation around the glycosidic linkage which makes the description of the furanose/pyranose interface
ACS Paragon Plus Environment
12
Page 13 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
straightforward. However, a series of other parameters are significantly changed and the flexibility of furanose ring is described in a way completely different from that applied in 56A6CARBO/CARBO_R for pyranoses rings. The combination of the TFR1 and TFR2 terms describes the intrinsic conformational properties of the furanose ring, independently of the presence or absence of any ring-substituents. This combination is inherent to all compounds exhibiting the topological pattern of furanose ring and remains unchanged throughout all possible furanosides. The T1 and T2 terms are kept; however, they are no longer used to describe torsions within the ring and their applicability is reduced to characterize the rotation of exocyclic groups. The applicability of the T3 term is also reduced to the case of hydroxyl torsions that do not involve any two ring atoms. The exocyclic torsions of such type (i.e. Xr‒CH*r‒O‒H) are described by the newly-introduced T3F term instead. The universal T4 term was removed and replaced by the two separate T4F1 and T4F2 terms, designed to reflect the exo-anomeric torsions for the hydroxyl and methoxyl groups, respectively. The T5, T5R1 and T5R2 terms are retrieved from 56A6CARBO/CARBO_R and follow the distinction into the three possible cases of anomeric torsions in unfunctionalized, O-methylated and O-glycosylated anomeric center, respectively. The remaining terms (T6-T13) are unchanged with respect to the 56A6CARBO/CARBO_R set. The presently proposed set of parameters, in analogy to all GROMOS sets, is designed for use with reaction-field electrostatics37,38 with an effective cutoff distance of 1.4 nm and an adjusted reaction-field dielectric permittivity equal to 61 and is compatible with the SPC water model.39 These conditions were accounted for during the validation procedure and other simulations reported in this work.
ACS Paragon Plus Environment
13
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 14 of 47
Fig. 3. Comparison of the energy profiles calculated at the classical (the presently proposed parameters) and QM (DFT/B3LYP/6-311++G**) levels of accuracy. The energy profiles correspond to the full rotation of the single exocyclic group with the ring shape conformationally locked in the 3T 2
shape. The QM- and force field-derived values are given in red and green, respectively.
ACS Paragon Plus Environment
14
Page 15 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Fig. 4. Comparison of the energy profiles calculated at the classical (the presently proposed parameters) and QM (DFT/B3LYP/6-311++G**) levels of accuracy. The energies correspond to the different ring geometries, expressed by the value of the pseudorotation phase (P). The QM- and force field-derived values are given in red and green, respectively.
2.4 Parameterization stages The parameterization was performed roughly incrementally, starting from the simplest compound (i.e. oxolane) to more complex compounds being the mono- or multisubstituted oxolane derivatives. Such procedure allowed for the refinement of successive parameter while freezing the previously adjusted ones. For simplicity, the intermediate attempts that were required to reach the final parameter set are not discussed, and only the procedure that led to final set is presented in more detail. The successive steps involved in this refinement were the following: A. Conformational properties of unsubstituted ring. After assigning atomic partial charges (derived from the bond-increment approach) and the LJ parameters (CH2r and Or atom types), the parameters for torsions within the oxolane ring were optimized on the basis of the QM/MM-MD-derived free energy profiles, describing the ring-distortion properties. This approach differs significantly from those usually accepted
ACS Paragon Plus Environment
15
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 16 of 47
during the GROMOS force field parameterization. Its justification, as well as other details related to this step, is given in section 3.1. Here, we only mention that the final set of parameters is capable of excellent reproduction of the experimental data (NMR-inferred 3JHH coupling constant values) related to the oxolane conformation in the solvent. At this stage the procedure of adjusting the parameters relied on minimizing the difference between the compared 2D (P vs. m) free energy profiles. Initially, the random search was performed in the space of the k and m values (see eq. (6)) in order to narrow the possible number of combinations. Then, the systematic search in the space of the two k parameters was performed (~120 tested combinations). This step required performing short (~150 ps) MD simulations for each of the tested parameter combinations to compare the resulting P vs. m free energy profiles with the reference one (i.e. derived from the QM/MM-MD calculations). B. Rotation of exocyclic groups. The initial torsional and LJ parameters were those adopted from the 56A6CARBO set. They were slightly refined (compare Tab. 6 in ref.15 and Tab. 3 in the present work) in order to improve the agreement with the QM energy scans (see Fig. 3). In particular, the N1 term for special LJ interactions was replaced by N1F whereas the T3 and T4 parameters were modified to become T3F, T4F1 and T4F1. Thus, instead of one universal term (T4), the distinction into the two separate cases of unfunctionalized and alkylated/glycosylated anomeric center was introduced to describe the exo-anomeric torsions in the corresponding compounds. Initially, the rotation of hydroxyl groups in 1-hydroxyoxolane, 1-methoxyoxolane, 2-hydroxyoxolane and 2-methoxyoxolane was considered. The two different orientations of the exocyclic group with respect to the conformationally frozen ring (3T2 shape) in the first two compounds were additionally taken into account (this is equivalent to the analysis of the same group orientation and the two limiting 3T2 and 2T3 ring geometries). After achieving a good agreement with the QM data, the resulting set was successfully tested in the context of rotation of the hydroxyl groups in the 1,2-dihydroxy- and 2,3-dihydroxyoxolane (see Fig. 3 for details). The procedure of adjusting the parameters relied on minimizing the difference between the energy profiles calculated either at the QM or classical levels. The systematic search in the space of the one or two k parameters was performed (~100-1000 tested combinations for each type of torsional angle). This only step required calculating the single-point energies for the same structures as those obtained during QM calculations.
ACS Paragon Plus Environment
16
Page 17 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
C. Relative free energies of ring conformers. The parameters concerning the LJ interactions of ring-linked hydroxyl oxygen atoms with the ring atoms were refined at this stage (i.e. the interactions involving Xr-OA and Xr-OE pairs). The torsional potential involving the CH*r-Or-CH*r-O topological term was unchanged with respect to the 56A6CARBO/CARBO_R set. The previously proposed16 refinements include the distinct parameters for both O-methylated, O-glycosylated and unfunctionalized compounds. The main difference with respect to the 56A6CARBO set it that the special LJ parameters describing the ring-substituent interactions (N2, N4 and N5 terms in ref.15) were replaced by the standard LJ parameters for the corresponding 1-4 or 1-5 neighbors. The validation at this stage relied mainly on the experimental data concerning ring conformational equilibria for various furanoses in water (quantitative analysis against the NMR-inferred 3JHH and 3JCH coupling constant values) and in crystal state (qualitative analysis against the compilation of crystal structures taken from the literature). Additionally, the results provided by the final set were compared with the QM energy scans performed along the pseudorotation phase coordinate for the three mono-substituted oxolane derivatives (see Fig. 4). At this stage the validation procedure concerned nearly whole group of considered compounds and required performing relatively long (~20 ns) MD simulations for all group members. Due to this fact, the final parameters were selected among a more limited number of combinations in comparison to the procedures described in points A and B. Instead of using the systematic search based on a grid of possible values of the k and/or LJ parameters, we preselected the small number of combinations describing the substituent-substituent and ring-substituent interactions basing on chemical intuition, the knowledge about the considered system and the existing pyranose-related parameters. The parameters adjusted at this stage were simultaneously validated with respect to the hydroxymethyl group rotation, anomeric equilibria and glycosidic linkage conformation (points D-F). D. Hydroxymethyl group rotation in furanoses. Seven different values of the k force constant of the T10/T11 term were tested in the context of agreement with the experimental data regarding the hydroxymethyl rotamer populations (more precisely: 0.5, 0.75, 1.0, 1.25, 1.5, 1.75 and 2.0 kJ/mol). Finally, the parameters for the corresponding exocyclic hydroxymethyl torsion were left unchanged in comparison to the 56A6CARBO set. This choice resulted in the best agreement with the hydroxymethyl rotamer populations of a series of aldopentofuranoses, either α- or β-anomers, in water.
ACS Paragon Plus Environment
17
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 18 of 47
E. Anomeric equilibria of furanoses. The parameters for the corresponding anomeric torsion were unchanged in comparison to the 56A6CARBO/CARBO_R set. This includes the T5, T5R1 and T5R2 terms which are expected to significantly affect both the anomeric equilibria and ring conformational properties.16 The final combination of parameters provided the satisfactory agreement with the experimental data concerning the anomeric populations of aldohexofuranoses, ketohexofuranoses and ketopentofuranoses in water. F. Glycosidic linkage conformation. The parameters adjusted/refined at the previous stages are sufficient to provide the description of the glycosidic linkage conformational properties. The additional validation was performed against the available 3JHH and 3JCH coupling constants measured for sucrose and the dimer of -D-arabinofuranose in the context of the pyranose-furanose α(1→2)- and furanose-furanose (1→5)-linkages, respectively. Due to the lack of the analogous NMR data for other types of furanose-related linkages, we qualitatively compared the force field predictions about the favorable glycosidic linkage geometry with the available crystal structures containing various di- and oligomers containing furanose units. 2.5 Computational details The compounds under consideration included: (i) the - and -anomers of all
D-
aldotetrofuranoses; (ii) the - and -anomers of all D-aldopentofuranoses; (iii) the - and anomers of all D-aldohexofuranoses; (iv) the - and -anomers of all D-ketohexofuranoses; (v) the - and -anomers of all D-ketopentofuranoses; (vi) oxolane and its mono-substituted derivatives; (vii) sucrose; (viii) disaccharide of α(1→5)-linked disaccharide of (2→6)-linked
D-fructofuranose;
D-arabinofuranose;
(x) disaccharide of (2→1)-linked
(ix) D-
fructofuranose. For compounds listed in points (i)-(v) their O-methylated derivatives containing the methoxyl group attached directly to the anomeric carbon atom were considered as well. The total number of furanose monomers that were considered was equal to 47. Figs. 1 and 5 show the chemical structures of the studied compounds.
ACS Paragon Plus Environment
18
Page 19 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Fig. 5. The furanose monomers that were simulated during the parameterization and/or validation procedures. Compound numbering is subsequently used throughout the article, including e.g. the results given in Figs. 8 and 9. Additionally, the study involved some furanose dimers, described in the text and shown in Fig. 1.
All MD simulations were carried out with the GROMACS 5.0 and GROMACS 2016.440 packages, in combination with PLUMED 2.1.5.41 The force field parameters used for simulations are described in the previous sections (2.3 and 2.4) and summarized in Tabs. 1-3. Molecular systems under consideration consisted of one saccharide molecule and approximately 1000 simple point charge (SPC) water molecules39 placed in a cubic (31 × 31 × 31 Å3) simulation box. In the case of oxolane, the additional simulations were performed with either CHCl342 or CH3CN43 as the solvent and increased box size (51 × 51 × 51 Å3 or 43 × 43 × 43 Å3, respectively) to accommodate approximately the same number of solvent molecules. The MD simulations were carried out under periodic boundary conditions and in the isothermal-isobaric ensemble (unbiased MD) or at fixed volume (metadynamics simulations). The temperature was maintained close to its reference value (298 K) by applying the V-rescale thermostat,44 whereas for the constant pressure (1 bar, semi-isotropic coordinate scaling) the Parrinello-Rahman barostat45 was used with a relaxation time of 0.4 ps. The equations of motion were integrated with a timestep of 2 fs using the leap-frog scheme.46 The solute bond lengths were constrained by application of the LINCS procedure47,48 with a relative geometric tolerance of 10-4. The full rigidity of the water molecules was enforced by application of the SETTLE procedure.49 The translational centerof-mass motion was removed every timestep separately for the solute and the solvent. The non-bonded interactions were calculated using a twin-range scheme,50 with short- and long-
ACS Paragon Plus Environment
19
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 20 of 47
range cutoff distances set to 0.8 and 1.4 nm, respectively, and an update frequency of 5 timesteps for the short-range pairlist and intermediate-range interactions. The reaction-field correction37,38 was applied to account for the mean effect of the electrostatic interactions beyond the long-range cut off distance, using a relative dielectric permittivity of 61 as appropriate for the SPC water model51 and of 4.8 or 37.5 for CHCl3 or CH3CN, respectively. All the systems were 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 production simulations. After equilibration, all unbiased simulations were carried out for 20 ns (monomers containing only methoxyl or hydroxyl exocyclic groups) or 100 ns (all remaining compounds) and the trajectory was saved every 1 ps. Apart from the regular GROMACS trajectory, the values of selected parameters were saved every 0.1 ps. These parameters included: the pseudorotation phase (P) and amplitude (m), defined according to Altona-Sundaralingam notation1 as well as the , , and dihedral angles. The probability distributions of all these parameters were calculated on the basis of unbiased MD simulations and the corresponding free energy profiles were obtained by applying the Boltzmann inversion. The 3JHH and 3JHC coupling constants were calculated by using either the Cremer-Wu (only oxolane),33 Houseknecht et al.52 (O1-methylated aldopentofuranoses) and AltonaHaasnoot53,54 (all compounds) relations and the unbiased MD trajectories generated for those compounds for which there exist the corresponding experimental data. The anomeric ratios for unfunctionalized aldoses were calculated according to the scheme described in ref.55 According to that scheme, some of the 'improper dihedral' terms inherent in the GROMOS force field are replaced by the double-well, stereoconfigurationindependent potential in order to describe the chirality of the selected carbon atoms. In this case the chirality of the anomeric carbon atom (C1, expressed by the CH1r type) is associated with the anomery of the given compound. The simulations were carried out by applying the well-tempered metadynamics approach and the resulting free energy profiles were based on the value of the γ dihedral angle defined by the following quadruplet of atoms: C1-C2-O4-O1. The anomeric populations were calculated from the differences between the two minima of free energy profile (located in the regions of the positive and negative values of γ and corresponding to the α and β anomers, respectively). At this stage, the parallel tempering scheme was combined56 with the well-tempered metadynamics approach57 as implemented in PLUMED 2.1.5. The well-tempered metadynamics,57 relied on Gaussian local functions of widths 2.86°, an initial deposition rate of 0.01 kJ mol−1·ps−1 and a temperature parameter ΔT ACS Paragon Plus Environment
20
Page 21 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
(see eq. (2) in ref.57) of 1788 K. The parallel-tempering58 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. The anomeric ratios for ketofuranoses were calculated by using the path-based welltempered metadynamics, implemented in PLUMED 2.1.541 and applied by us previously in the analogous simulations of ketohexopyranoses.59 The two variables are defined by the path (i.e. the collection of molecular structures) connecting the two distinct states (here: the α- and β-anomers). The 20 different, path-inherent structures contained only the ring atoms and the atoms neighboring to the anomeric carbon. The calculated 2D free energy landscapes relied on the two coordinates (s and z) that describe the position of the system along the path (s) and the deviation from this path (z). The well-tempered metadynamics,32,56,57 relied on Gaussian local functions of widths 0.4 (s) and 3·10−4 (z), an initial deposition rate of 2 kJ mol−1·ps−1 and a temperature parameter ΔT of 14602 K. The value of the λ parameter was set as 100. The anomeric populations were calculated from the differences between the two minima of the converged free energy landscape. The QM calculations were performed within Gaussian0960 at the DFT/B3LYP/6311++G** level of theory.61–65 The calculations included the energy scans with the value of single dihedral angle as the coordinate or the energy scans along the pseudorotation phase with the ring shape frozen in each of the 20 canonical shapes (Fig. 1).
3. Results and Discussion 3.1. Pseudorotation of oxolane Due to the GROMOS force field build-up rules and the accepted atomic partial charges (see refs.15), the oxolane (tetrahydrofuran) molecule is described by a model containing five atoms (one ring oxygen atom and four united atoms of the CH2r type) that exhibit no mutual nonbonded interactions. This significantly reduces the possibilities to adjust the force field parameters accurately describing the conformational properties of oxolane. As the result, only the parameters associated with the torsional dihedral angles (eq. (6)) can be adjusted. Another complication is connected with the controversies about the pseudorotational energy profiles calculated for oxolane by means of quantum-chemical approaches of various types. Namely, depending on the accepted method, one can obtain profiles varying not only quantitatively but also qualitatively. More details on this subject can be found in ref.
66
Note
that all works mentioned therein relied on the ab initio methods of high accuracy, probably
ACS Paragon Plus Environment
21
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 22 of 47
even higher than those usually accepted during force field parameterization procedures. One of the possible reason for such diverse results is that the magnitude of energy changes along the pseudorotational coordinate is very low; according to the data compilation in ref.
66
this
quantity can vary within the range of 0-0.8 kJ/mol, depending on the system. Very low inaccuracies (i.e. of order of 1 kJ/mol) inherent in the given method of its implementation may result in qualitative differences between final profiles.
Fig. 6. Comparison of the 1D free energy profiles calculated at the classical (the presently proposed parameters) and QM/MM (the three different ab initio potentials: BP86/def2-SVP, MP2/UFH/631G** and MP2/UHF/cc-pVTZ) levels of accuracy for the oxolane molecule in explicit water. The QM/MM data were taken from ref.67 The free energy profiles correspond to the two descriptors of the five-membered ring shape, i.e. the pseudorotation phase (panel (A)) and amplitude (panel (B)). The force field-derived values are given as red lines.
The independent calculations performed by us confirm the high dependence of the oxolane conformational energies on the computational setup.67 Interestingly, even when using the most accurate method (i.e. MP2/UHF/cc-pVTZ) recommended in ref.,66 we were not able to quantitatively recover the results reported there. However, the authors of the original work emphasized the exceptionally high sensitivity of the oxolane pseudorotational energy profiles on the existing constraints or other aspects of computational setup. This, in practice, significantly limits the applicability of the routine procedures of the QM energy scans based on applying one (or more) constrains and optimizing the remaining degrees of freedom.
ACS Paragon Plus Environment
22
Page 23 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Finally, our recent studies67 have indicated that due to the solvent-induced charge flux within the furanose molecules, the properties calculated in the absence of explicitly-defined solvent may drastically differ from those corresponding to the condensed-phase systems. In view of these difficulties, we decided to use the ab initio-derived free energy profiles corresponding to pseudorotation of the oxolane molecule in the aqueous solvent which have been calculated in our recent, independent study.67 The adjustment of the classical force field parameters relied on recovering these two-dimensional free energy profiles with respect to their characteristic features (e.g. location and depth of the free energy minima, location and height of the free energy barriers, average puckering amplitude, etc.). Deriving the force-field parameters was based on the MD simulations of oxolane in aqueous solvent whereas their further validation relied on the simulations of the systems containing either CHCl3 or CH3CN as the solvent (see the “Methods” section for details). The final set of parameters is capable to correctly predict the types and populations of both favorable and unfavorable ring conformers and the puckering amplitude. Moreover, the theoretically calculated NMR 3JHH coupling constants agree well with the experimentally-inferred values68 (see Fig. 9 and Tab. S1 in Supporting Information for details). Therefore, as a final set, accurately describing the intrinsic flexibility of unsubstituted five-membered furanose ring, we accepted the combination of the five-fold potentials, listed in Tab. 3. Although applying such unusual multiplicity in the torsional potentials seems to be counter-intuitive and against the symmetry of the bonds within the ring, only potentials of this type are capable of reproducing not only the average puckering amplitude value but also the full free energy profile related to this quantity. Interestingly, in contrary to the recent results obtained within the GLYCAM force field dedicated to furanoses,24 our study indicates that the most populated conformations of the oxolane ring are located at the north and south of the pseudorotation wheel (Fig. 1) rather than on its east and west. That former behavior is compatible with the two-state model for conformation of the five-membered rings, assuming that only north and south regions of the pseudorotation itinerary are non-negligibly populated, whereas the energy barriers exist at the east and west. However, during testing a large number of alternative combinations of force field parameters (rejected at the further stages), we found that the reasonable agreement with the experimental data (3JHH coupling constants) can be obtained even if the pseudorotation behavior was contradictory to that presented in Fig. 6 (the average deviation not exceeding 1.5 Hz). The necessary condition for that seems to be the existence of large populations of both east/west and south/north-like conformers. Thus, we can conclude that even if the two-
ACS Paragon Plus Environment
23
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 24 of 47
state-like behavior is characteristic for the conformational equilibrium of oxolane, the related model would not be a reliable approximation, due to the large population of the ring shapes located on the free energy barriers for pseudorotation. Finally, it should be emphasized that the above-mentioned issues are associated mainly with the oxolane molecule. The substitution of the oxolane ring by at least one hydroxyl group changes its pseudorotational properties and, what more important, reduces the differences between pseudorotational energy profiles calculated by means of various ab initio methods. Analogously, the other types of energy profiles (e.g. those associated with the rotation of exocyclic dihedral angles) are not very sensitive for the choice of ab initio approach when considering the relative differences between energy minima. As the result, the remaining ab initio-based energy scans, not concerning the oxolane molecule, were performed by applying the DFT/B3LYP/6-311++G** level of theory. 3.2. Rotation of Exocyclic Groups The rotation of hydroxyl and methoxyl group attached directly to the furanose ring was studied in the context of the - and -anomers of
D-arabinofuranose
and
D-ribofuranose.
Additionally, both their unfunctionalized and O1-methylated forms were considered. The results are given in Fig. 7 as the probability distribution plots related to the value of the dihedral angle. Although, to our best knowledge, there is no experimental data that would provide the quantitative populations of the rotamers of the lactol group, the corresponding conformational properties can be checked against the preferences resulting from the presence of the exo-anomeric effect.69 Regarding the D-furanose molecule, due to this effect, the most populated rotamers are expected to be g- for the β-anomers (i.e. ≈ -60 deg) and g+ for the αanomers (i.e. ≈ +60 deg), whereas the t rotamer (i.e. ≈ ±180 deg) should always be lowpopulated. The results given in Fig. 7 show that the presently proposed parameter set reflect the expected behavior. The methylation at O1 results in even larger population of the dominating conformer (i.e. g+ for the α- and g- for β-anomers) over the remaining two in comparison to the unfunctionalized compounds. The trend is analogous to that found for pyranoses in a series of separate studies,15,17,59,70 and can be ordered as follows: g+ > t > g- (anomers); g- > t > g+ (unfunctionalized β-anomers) and g- > t ~ g+ (O1-methylated β-anomers).
ACS Paragon Plus Environment
24
Page 25 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Fig. 7. The probability distribution plots related to the values of the dihedral angle that describes the conformation of the lactol group. All the plots correspond to the unbiased MD simulations in the explicit water solvent.
3.3. Rotation of Hydroxymethyl Group The results concerning the hydroxymethyl group rotation for furanoses in water are shown in Fig. 8 in the context of the O1-methylated - and -anomers of aldopentofuranoses, including 2-deoxyribose. Fig. 8 contains the comparison of the calculated populations of hydroxymethyl group conformers with corresponding experimental estimates from the literature.71,72 The values illustrated graphically in Fig. 8 are also given in Tab. S2 (Supporting Information). For most of the considered compounds the trend in populations of the particular conformers is the same for both the experiment and results of MD simulations. A relatively small divergence (up to 20%) between theory and experiment that changes the order of rotamer populations can be observed for -D-2-deoxyribose and -D-xylose. Nevertheless, the level of agreement between the simulation results and the available experimental data is satisfactory and confirms that the force field predictions are accurate. The average deviation from the experimental data is only 8% per one compound, with the maximum value of 29% in the case of -D-xylose (note that in this case such large deviation does not change the type of dominating conformer). The predicted trend in populations of the three staggered conformers is similar for most of the considered compounds and indicates that the t conformer is usually the most populated one (for 7 out of 9 compounds). In the case of -D-ribose, the lack of any exocyclic groups at the same side of the ring with respect to the hydroxymethyl group causes the increase of the g+ conformer population which becomes the highest one. The same trend is observed for structurally similar -D-2-deoxyribose where the relatively balanced populations of the g+ and t conformers (41 vs. 48 %) are observed.
ACS Paragon Plus Environment
25
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 26 of 47
Contrary to the pyranoses, in furanoses the anomery of the given compound seems to be a relevant factor that influences the population of the hydroxymethyl group rotamers. More precisely, the -anomers always exhibit larger population of the t conformer in comparison to the corresponding -anomers. Moreover, the populations of the remaining two conformers can also be significantly altered (up to 25%) upon the inversion of anomeric center and this effect may lead to reordering in the population trends (e.g. for arabinose). As the distance between the hydroxymethyl and lactol moieties is of similar magnitude in comparison to that of pyranoses, it is most likely that this effect is caused by the large flexibility of the furanose ring which, by adopting different geometries, permits for more intensive interactions between lactol and hydroxymethyl groups. It is worth noting that in analogy to pyranoses, the g+ conformer, which would potentially allow for an intramolecular hydrogen bonding between the hydroxymethyl group and the hydroxyl group at C3 is not dominating for any of the compounds exhibiting the suitable orientation of the hydroxyl group. This is another example70,73,74 of limited significance of intramolecular hydrogen bonding for the conformational preferences of carbohydrates in polar solvents.
Fig. 8. The experimentally-measured (transparent bars) and calculated (solid bars) populations of the three staggered conformers of the hydroxymethyl group (in [%]). The calculations relied on the unbiased MD simulations within the presently proposed force field in the presence of explicit water. Numbering of compounds is in accordance with Fig. 5.
3.4. Ring Conformation The partial results concerning relative free energies of ring conformers for various groups of furanoses (see Methods for details) in water are shown in Figs. S1 (as free energy profiles) and 9 (as the values of NMR-inferred coupling constants, related to the ring geometry). Contrary to pyranoses, the five-membered furanose rings are much more flexible, adopting various conformations separated by small or no free energy barriers on the pseudorotation path. This seriously complicates the quantitative description of the conformational equilibrium within the furanose rings. There exists a simple two-state model,
ACS Paragon Plus Environment
26
Page 27 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
based on empirical observation,1,24 assuming that the pseudorotation itinerary can be divided into two hemispheres (north vs. south) which correspond to the two relatively narrow distributions of the possible ring conformers. However, there are number of growing evidences that this popular model may be oversimplified for some groups of furanoses.24,27– 30,67.
Thus, we limit ourselves in the analysis of the ring-distortion properties by providing
only basic findings and presenting the complete force field validation based on the NMR data. The comparison of the experimentally measured J-coupling constant (3JHH and 3JCH) values24,52,68,75–81 with their theoretically-calculated counterparts is given in Tab. S1 (numerical values) and Fig. 9 (graphical illustration). The comparison was performed separately for each group of the considered furanoses, namely: (A) oxolane and its monosubstituted derivatives; (B) O1-methylated and unfunctionalized (C) O1-methylated
D-aldopentofuranoses;
(D) O1-methylated
D-aldotetrofuranoses;
D-aldohexofuranoses;
(E) O2-
methylated D-fructofuranose; (F) other compounds (45-47) that do not belong to any of the above groups. For groups (B)-(E), both the - and -anomers were considered. The simulation results for the simplest compounds (group (A)) show a good agreement with the experimental 3JHH values. The average deviation between the experiment and theory for this group is equal to 1 Hz. Especially good match was obtained for oxolane, where the Wu-Cremer equation33, dedicated to oxolane was used instead of the more general Haasnoot relstionship.53,54 Relatively larger deviations (> 2 Hz) appear in the case of 3JHH measured for some of the H3-H4 proton pairs for monosubstituted oxolane derivatives. However, the 3JHH values related to the same proton pairs are reproduced very well in the remaining groups of compounds (see further paragraphs in the present section). For group (B) also a very good agreement has been obtained, with average deviations equal to 0.9 and 1.1 Hz for unfunctionalized and O1-methylated tetroses, respectively. The maximal deviation does not exceed 1.9 Hz. Interestingly, this is the only group of furanosides for which the experimental data for both unfunctionalized and methylated compounds are available. The insight into the measured NMR data shows that the influence of O1methylation on the ring-distortion properties is very limited and leads to very minor shifts in the measured 3JHH values. Namely, the average shift equals to 0.38 Hz per data point for all group B, whereas the maximal deviation does not exceed 0.8 Hz. This is in contrast with the results of MD simulations concerning aldohexopyranoses and their ring-inversion properties which seem to be highly-dependent on the functionalization at O1.16 There exists a rationalization of these results on the ground of enhancement of the anomeric effect upon the presence of substituent at O116 but the magnitude of such effect is differently reported in
ACS Paragon Plus Environment
27
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 28 of 47
refs.16,70,82. However, there are no related studies focused on furanosides which would confirm or refute the existence of the same effect and describe it in terms of its possible influence on the conformational equilibrium within the ring. Therefore, in view of this fact, supported by the above mentioned minor impact of the methylation on the 3JHH values, we decided to change the corresponding description of the force field terms in comparison to the GROMOS 56A6CARBO_R.16 More precisely, we introduced a uniform pattern for the LJ interactions between ring atoms and the ring-attached oxygen atoms. Such combination provides the reasonably well agreement with the experimental data for both unfunctionalized and methylated tetrofuranoses. Although there is no NMR data measured for unfunctionalized furanoses of other considered groups, we have kept the same parameters also for them as dictated by the uniform force field built-up rules. The main effect of this assignment is that, in accordance with the experimental observations, the O-methylation at the anomeric carbon does not significantly change the conformational equilibrium within the furanose ring. The only exception from that rule is the case of -psicofuranose, where the depths of the corresponding local minima of free energy profiles are significantly altered upon O2methylation (see Fig. S1). The agreement with experiment within group (C) is also satisfactory; when considering the 3JHH coupling constants, the average deviation is 1.3 and 0.8 Hz for - and anomers, respectively. However, contrary to group (B), one can observe one case of a very large deviation (> 4 Hz) exhibited by -D-arabinose. After examining this case, we concluded that the exceptionally high experimental value of 3JHH is beyond the range of applicability of the used variation of the Karplus equation. Actually, to achieve the deviation of range of 1 Hz, one would assume that the ring is nearly completely conformationally locked; this is against the results obtained for the remaining two pairs of protons. Moreover, the 3JCH values, originating from a competitive study are reproduced very well by the MD data combined with the variation of the Karplus equation dedicated for furanosides (average deviation equal to 0.2 Hz). The same type deviation is also observed in the case of GLYCAM force field.24 Finally, the comparison of the pseudorotation free energy profile calculated for -D-arabinose (see Fig. S1(B)) with the compiled data extracted from the X-ray structures (Fig. 9 in ref.83) revealed the agreement between the most favorable ring conformers (i.e. shapes close to 1T2 and 2E). Simulations concerning group (D) represent a similar level of agreement with experiment as those of group (C). The average deviation calculated for the whole group is
ACS Paragon Plus Environment
28
Page 29 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
equal to 0.8 Hz, whereas the only significant deviation (equal to 3.5 Hz) is observed for the 3J H1H2
value measured for -D-allofuranose. The MD results for D-fructofuranose (group (E)) reproduce the experimental data very
well for the -anomer (deviation equals to 0.4 Hz) but the agreement for the -anomer is much poorer (deviation exceeds 3 Hz). Again, after examining the limits of applicability of the Karplus equation, we concluded that to achieve the perfect agreement with the experiment, the assumption of the ring rigidity is necessary. Furthermore, the obtained free energy profile of pseudorotation of -D-fructofuranose is in a very good agreement with the compilation of the related X-ray data (see Fig. 10 in ref.83), predicting the most favorable ring conformers to be located around 2E. The results related to the most diverse group (F) of other furanosides will not be discussed in details but the good agreement of the simulation results with the experimental values of 3JHH (average deviation equals to 0.6 Hz, whereas the maximal deviation does not exceed 0.9 Hz) is worth noting. Although the X-ray diffraction (XRD) structures cannot be treated in terms of direct evidence for the conformational behavior of furanosides in solutions, they still may provide valuable qualitative information about the conformational preferences of the ring which includes the insight into the most favorable geometries or conformationally-restricted regions. We performed the comparison based on the data compilation gathered in ref.
83
(see Figs. 9
and 10 therein) which concerns D-aldopentofuranoses and D-fructofuranose. Due to the scarce amount of data, we did not perform the analogous comparison for the remaining ketohexofuranoses. In all the cases, except of -D-xylofuranose, a very good qualitative agreement between the free energy profiles shown in Fig. S1 and the corresponding XRD data can be seen. Interestingly, the match between the experimentally- and MD simulation-derived 3J
values for -D-xylofuranose is among one of the best found (deviation per data point equals
to 0.8 Hz) which indicates that the ring conformational properties in solution are accurately reproduced by the force field. Apart from the pseudorotation phase, one can additionally focus on the pseudorotation amplitude. Also in the case of this quantity, the good agreement between the XRD data and the free energy profiles can be observed; namely, the largest number of structures is located around m = 40 deg with the typical deviation of order of 10 deg. This corresponds to the location and relative depths of the free energy minima on the calculated profiles (Fig. S1).
ACS Paragon Plus Environment
29
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 30 of 47
In conclusion, the reproducibility of the ring-distortion-related NMR experimental data by the presently proposed force field is satisfactory. The best agreement has been obtained when using the variations of the Karplus equation that have been developed for special group of furanosides (e.g. equation proposed in ref. pseudorotation phase-based equation in ref.
33
52
for aldopentoses or the
for oxolane). In the remaining cases some
larger deviations are observed but the agreement for the whole group of compounds is still acceptable. This is additionally confirmed by the qualitative insight into the XRD data collected for various furanosides. The analysis of the free energy maps provides another example of inapplicability of the two-state model of conformational equilibrium in furanose rings. The two-state model assumes a relatively narrow distribution of ring shapes into the north and south hemispheres of the pseudorotation itinerary. Such behavior has been proposed in 1972 on the basis of analysis of the XRD structures.1 However, since then there appeared many examples of its poor applicability, particularly in the case of -D-arabinofuranose (see ref.
24
and references
therein) which complicates the analysis of the NMR data and the description of the fivemembered ring distortion, in general. Recently, we have shown that the good applicability of the two-state model may be explained on the ground of the solvent-induced flux of the atomic charges within the furanose ring.67 However, the height of the north vs. south energy barrier generated in this way is rather low (2.5-6 kJ/mol, depending on the considered compound). Thus, the interactions between ring substituents become more important for the ring conformation. Our simulations confirm that no north vs. south free energy barrier exists in the case of compounds exhibiting the orientation of the hydroxyl ring substituents characteristic for -D-arabinofuranose, for instance: -D-fructofuranose, -D-altrofuranose and -Dthreofuranose (see Fig. 10). However, both quantitative and qualitative patterns of deviations from the two-state model seem to be much more complicated to be solely ascribed to this specific combination of substituent dispositions. This is shown by examples of free energy profiles exhibited by -D-psicofuranose, -D-threofuranose, -D-ribofuranose (no north vs. south free energy barrier) or -D-galactofuranose (existence of the east vs. west free energy barrier). We postpone the detailed description of the influence of the substituent type and its orientation on the five-membered ring conformation to our future work, presently limiting ourselves to only brief comments on the general features of the calculated free energy maps. In the case of all compounds the free energy minima do not correspond to the single canonical
ACS Paragon Plus Environment
30
Page 31 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
conformers (see Fig 1). They are rather spanned over several of them, forming large basins of favorable ring shapes with no separating energy barriers. The typical number of ring conformers located in the most favorable free energy minimum (energy no larger than ~RT, i.e. ~2.5 kJ/mol with respect to the minimal value) varies between 5-7 (tetroses), 2-5 (aldohexoses), 3-6 (aldopentoses), 3-7 (ketohexoses) and 3-4 (ketopentoses). In general, the compounds with smaller number of ring substituents (e.g. tetroses) seem to exhibit a larger degree of ring flexibility, comparable to that of unsubstituted ring of oxolane. The puckering amplitude typically varies within the range of 20-60 deg whereas the lower and higher values of m are conformationally inaccessible. The average amplitude is close to 40 deg, however, the most favorable values of m seem to be dependent on the actual P value which is in agreement with the results of the recent QM/MM study.67 The large conformational heterogeneity covering a number of different ring shapes within single free energy minimum as well as the described above deviations from the two-state models suggest to look for a new description pattern for conformational features of the five-membered rings that would simplify the calculated free energy profiles more accurately.
ACS Paragon Plus Environment
31
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 32 of 47
Fig. 9. J-coupling constant (3JHH and 3JCH) for the group of various furanoses in water, measured experimentally (transparent bars) or calculated by MD simulation using the presently proposed force field (solid bars). The number above each bar corresponds to the given compound, according to the notation in Fig. 5. The red and green bars represent the - and -anomers, respectively. The experimental values were taken from refs.52,68,75–78,80,81 The values from the MD simulations were calculated by application of the Karplus equation by Haasnoot et al.,53,54 Wu and Cremer33 and Houseknecht et al.52
ACS Paragon Plus Environment
32
Page 33 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Fig. 10. The exemplary two-dimensional free energy maps calculated for compounds that do not obey the two-state-like behavior. The coordinates are: pseudorotation phase (expressed as the ring geometry type) and the pseudorotation amplitude (distance from the wheel center; linear scale changing from 0 to 60 deg). (A,B) Compounds that exhibit the disposition of substituents characteristic of arabinofuranose. (C-F) Compounds that exhibit the disposition of substituents different than that characteristic of -arabinofuranose. (F) Free energy map exhibiting the barrier along the north-south axis. (C,D) The comparison of maps calculated for unfunctionalized and O1-methylated compounds reveals the minor influence of the O1-substituent. The energy scale is in [kJ/mol].
ACS Paragon Plus Environment
33
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 34 of 47
3.5. Anomeric Equilibria of Furanoses The results concerning the anomeric equilibrium of aldohexofuranoses, ketopentofuranoses and ketohexofuranoses of the
D
series in water are gathered in Tab. 4. The force field
reproduces reasonably well the experimental anomeric population values with an average deviation of 14%. The most significant errors correspond to glucofuranose (49%), talofuranose (32%) and allofuranose (29%), i.e. the compounds belonging to the group which tautomeric equilibrium is strongly shifted toward pyranose form. Additionally, for the same three compounds, these deviations lead to an inversion of the α:β ratio. Although such differences seem to be exceptionally large, they represent relatively small errors in the free energies corresponding to the inversion of the anomeric center. When converting the anomeric populations into corresponding free energy values, the maximal deviation between the simulation and experiment does not exceed ~5.9 kJ/mol (glucofuranose) with the average deviation equal to 2.1 kJ/mol. When neglecting the above-mentioned three compounds, the deviation drops to only 1.1 kJ/mol per compound. The magnitude of discrepancy observed in the case of the present study is in line with those
observed
previously
for
aldohexopyranoses,15,70
ketohexopyranoses,
aldopentofuranoses, deoxy- and dideoxyaldopentopyranoses59 as well as for the indirect results concerning acyclic forms of aldo- and ketohexoses.84 In conclusion, the present force field can be used to study the anomeric equilibrium of furanoses, however, some additional refinements (of character and magnitude dependent on the compound of interest) may be necessary to increase the accuracy of calculations. Tab. 4. The anomeric populations (in [%]) calculated for selected furanosides of the D series. The experimental data were taken from refs.85–87
compound allose altrose galactose glucose gulose idose mannose talose
aldohexoses experiment α β 36 58 38 28 24 43 73 62
64 42 62 72 76 57 27 38
simulations α β 65 61 39 77 27 24 70 30
35 39 61 23 73 76 30 70
compound fructose psicose sorbose tagatose compound ribulose xylulose
ketohexoses experiment α β 18 82 72 28 100 0 36 64 ketopentoses experiment α β 75 25 23 77
simulations α β 5 95 64 36 90 10 21 79 simulations α β 73 27 3 97
ACS Paragon Plus Environment
34
Page 35 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
3.6. Conformation of Glycosidic Linkages The results concerning the conformation of the glycosidic linkages are given in Tab. 5 and in Fig. 11. Considering the NMR-inferred experimental data measured for furanoses in solution, the force field reproduces very well the 3JCH and 3JHH values recorded for sucrose and Me-O1di--D-arabinofuranose, respectively (Tab. 5). The maximal deviation equals to 0.9 Hz whereas the average deviation does not exceed 0.4 Hz. Appreciating the agreement with the experiment one has to note that our comparison concerns only the two types of furanoseinvolving glycosidic linkage out of many possible. As the analogous NMR data for glycosidic linkages of other topologies are not available, we performed the additional, qualitative comparison of the results of MD simulations with the available structural data found in the PDB database. The results are given in Fig. 11. The location of the conformationallypreferred and restricted regions is satisfactory reproduced for the three types of considered linkages. The most significant deviations (i.e. the database-extracted coordinate values located at the edge of MD-most intensively-sampled region) are observed for the compounds containing the (2→6) type of linkage. However, it is worth to note that for these compounds only very limited number of structures (4) is available. Interestingly, the conformational equilibrium of inulin in water seems to be shifted in the direction characterized by the dihedral angle value close to +60 deg rather than to -60 deg, as reported in the crystallographic studies (see Fig. 11). Such behavior is expected on the basis of the results reported in Section 3.3 and characteristic of all furanoses that contain the hydroxymethyl moiety in their structures. The XRD-like conformation appears as the secondary minimum of the free energy. Currently, no experimental data related to the inulin, 1-kestose or nystose in water exist to validate this observation. We plan to revisit this issue by performing the QM/MM-MD simulations concerning the conformation of the corresponding glycosidic linkage and/or concerning the larger chains, better reflecting the molecular topology of inulin.
ACS Paragon Plus Environment
35
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 36 of 47
Tab. 5. 3JHH coupling constant for the furanose-containing disaccharides in water, measured experimentally or calculated by MD simulation using the presently proposed force field. The experimental values were taken from refs.88,89 The values from the MD simulations rely on the Karplus equation53,54 and unbiased simulations in water. Coupling type
Me-O1-di--D-arabinofuranose (48) reducing end
non-reducing end
experiment simulation experiment simulation
3J H4H5R
5.5
5.9
5.5
4.6
3J H4H5S
3.4
3.4
3.4
3.6
ap
Coupling typea
3J p f H1 C2
sucrose (49) experiment simulation 4.2
4.8
– pyranose ring, f – furanose ring
Fig. 11. The values of the and dihedral angles, describing the orientation of the glycosidic linkages in furanose-containing di- and oligosaccharides, according to the XRD data deposited in the PDB database (squares). The black dots represent the results of the unbiased MD simulations concerning disaccharides exploiting the corresponding type of linkage in the explicit water.
ACS Paragon Plus Environment
36
Page 37 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
4. Conclusions The main aim of this study was to propose and validate the new parameter set dedicated for the simulation of furanose-based carbohydrates. The developed united-atoms GROMOS force field was designed to reproduce the conformational properties of wide class of furanosides in the presence of aqueous solution. This includes the conformation of exocyclic hydroxyl, methoxyl and hydroxymethyl groups, conformational equilibrium within the ring and the conformation of the furanose-involving glycosidic linkages. Additionally, the proposed set is capable to reasonably well reproduce the anomeric equilibria of furanoses. The proposed parameterization scheme uses the built-up rules introduced in the GROMOS 56A6CARBO set with a series of modifications which concern mainly the parameters involving the furanose ring. The proposed force field can be combined with the existing GROMOS 56A6CARBO/CARBO_R set to formulate the model description for more complex carbohydrates, containing not only furanose but also pyranose units and accounting for the different glycosidic linkage topologies, possibility of functionalization (alkylation) as well as glycosylation (including branched chains). In particular, the new set of parameters can be applied to study a series of biologically-relevant di-, oligo- and polysaccharides containing furanose units, such as: sucrose, nystose, raffinose, 1-kestose, 6-kestose, inulin, levan, pectins, etc. Finally, the results gathered during the validation procedure provide several important findings regarding the conformational features of furanose carbohydrates: 1. The intrinsic conformational properties of unsubstituted furanose ring (oxolane) indicate that the most favorable are the 3T2 and 2T3 conformers with the free energy barriers (~2 kJ/mol) located around the OE and EO shapes. 2. The conformational equilibrium in the furanose rings often deviates from the twostate model, assuming the north vs. south division of the pseudorotation itinerary. These deviations are not limited to the compounds exhibiting arabino-like topology of the substituents. Moreover, the two-state model is rather poorly applicable for the elastic rings exhibiting non-negligible conformer populations at the east and/or west regions even if those regions correspond to the expected free energy barrier (e.g. oxolane). 3. Contrary to pyranoses, in the case of furanoses, the O-methylation at the anomeric carbon atom usually does not lead to any significant alterations in the ring-distortion properties.
ACS Paragon Plus Environment
37
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 38 of 47
4. The trends in population of the staggered conformers of both the lactol and hydroxymethyl groups reflect the findings reported earlier for pyranoses. 5. The conformation of the (2→6) glycosidic linkage in water seems to systematically deviate from that known from the crystal structures. Instead, the most favorable conformation of the exocyclic dihedral angle corresponds to that characteristic of aldopento- and ketohexofuranose monomers (i.e. the t conformation); the conformer indicated by the crystallographic data corresponds to the secondary minimum of free energy. A more systematic studies aimed at issues mentioned above are planned to be the subject of our future works.
Supporting Information The Supporting Information is available free of charge on the ACS Publications website. Additional tables and graphs concerning the numerical values for JHH and JCH coupling constants, ring-distortion free energy profiles and the exemplary force field topologies (for oxolane, sucrose and monomers of: -D-threofuranose, Me-O1--D-threofuranose, -Darabinofuranose, Me-O1--D-arabinofuranose, -D-altrofuranose, Me-O1--D-altrofuranose, -D-xylulofuranose,
Me-O2--D-xylulofuranose,
-D-fructofuranose,
Me-O2--D-
fructofuranose) in the GROMACS format.
Acknowledgments The authors acknowledge the financial support of the Polish National Science Centre (contract financed in 2016–2020 under Project No. 2015/18/E/ST4/00234 SONATA BIS).
References (1)
Altona, C.; Sundaralingam, M. Conformational Analysis of the Sugar Ring in Nucleosides and Nucleotides. New Description Using the Concept of Pseudorotation. J. Am. Chem. Soc. 1972, 94 (23), 8205–8212.
(2)
Langridge, R.; Marvin, D. A.; Seeds, W. E.; Wilson, H. R.; Hooper, C. W.; Wilkins, M. H. F.; Hamilton, L. D. The Molecular Configuration of Deoxyribonucleic Acid: II. Molecular Models and Their Fourier Transforms. J. Mol. Biol. 1960, 2 (1), 38-IN12.
ACS Paragon Plus Environment
38
Page 39 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
(3)
Levitt, M.; Warshel, A. Extreme Conformational Flexibility of the Furanose Ring in DNA and RNA. J. Am. Chem. Soc. 1978, 100 (9), 2607–2613.
(4)
Meutermans, W.; Le, G. T.; Becker, B. Carbohydrates as Scaffolds in Drug Discovery. ChemMedChem 2006, 1 (11), 1164–1194.
(5)
Golgher, D. B.; Colli, W.; Souto-Padrón, T.; Zingales, B. Galactofuranose-Containing Glycoconjugates of Epimastigote and Trypomastigote Forms of Trypanosoma Cruzi. Mol. Biochem. Parasitol. 1993, 60 (2), 249–264.
(6)
Carapito, R.; Imberty, A.; Jeltsch, J.-M.; Byrns, S. C.; Tam, P.-H.; Lowary, T. L.; Varrot, A.; Phalip, V. Molecular Basis of Arabinobio-Hydrolase Activity in Phytopathogenic Fungi Crystal Structure and Catalytic Mechanism of Fusarium Graminearum GH93 Exo-α-L-Arabinanase. J. Biol. Chem. 2009, 284 (18), 12285– 12296.
(7)
Clemens, R. A.; Jones, J. M.; Kern, M.; Lee, S.-Y.; Mayhew, E. J.; Slavin, J. L.; Zivanovic, S. Functionality of Sugars in Foods and Health. Compr. Rev. Food Sci. Food Saf. 2016, 15 (3), 433–470.
(8)
Cooper, P. D.; Petrovsky, N. Delta Inulin: A Novel, Immunologically Active, Stable Packing Structure Comprising β-D-[2 → 1] Poly(fructo-Furanosyl) α-D-Glucose Polymers. Glycobiology 2011, 21 (5), 595–606.
(9)
Behrends, R.; Kaatze, U. Molecular Dynamics and Conformational Kinetics of Monoand Disaccharides in Aqueous Solution. ChemPhysChem 2005, 6 (6), 1133–1145.
(10) Immel, S.; Lichtenthaler, F. W. Molecular Modeling of Saccharides, 7. The Conformation of Sucrose in Water: A Molecular Dynamics Approach. Liebigs Ann. 1995, 1995 (11), 1925–1937. (11) Xia, J.; Case, D. A. Sucrose in Aqueous Solution Revisited, Part 1: Molecular Dynamics Simulations and Direct and Indirect Dipolar Coupling Analysis. Biopolymers 2012, 97 (5), 276–288. (12) Rehm, B. H. A., Ed. Microbial Production of Biopolymers and Polymer Precursors: Applications and Perspectives; Caister Academic Press: Norfolk: U.K., 2009. (13) Lins, R. D.; Hünenberger, P. H. A New GROMOS Force Field for HexopyranoseBased Carbohydrates. J. Comput. Chem. 2005, 26 (13), 1400–1412. (14) Pol-Fachin, L.; Rusu, V. H.; Verli, H.; Lins, R. D. GROMOS 53A6GLYC, an Improved GROMOS Force Field for Hexopyranose-Based Carbohydrates. J. Chem. Theory Comput. 2012, 8 (11), 4681–4690.
ACS Paragon Plus Environment
39
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 40 of 47
(15) Hansen, H. S.; Hünenberger, P. H. A Reoptimized GROMOS Force Field for Hexopyranose-Based Carbohydrates Accounting for the Relative Free Energies of Ring Conformers, Anomers, Epimers, Hydroxymethyl Rotamers, and Glycosidic Linkage Conformers. J. Comput. Chem. 2011, 32 (6), 998–1032. (16) Plazinski, W.; Lonardi, A.; Hünenberger, P. H. Revision of the GROMOS 56A6CARBO Force Field: Improving the Description of Ring-Conformational Equilibria in Hexopyranose-Based Carbohydrates Chains. J. Comput. Chem. 2016, 37 (3), 354–365. (17) Panczyk, K.; Gaweda, K.; Drach, M.; Plazinski, W. Extension of the GROMOS 56a6CARBO/CARBO_R Force Field for Charged, Protonated, and Esterified Uronates. J. Phys. Chem. B 2018, 122 (14), 3696–3710. (18) Guvench, O.; Greene, S. N.; Kamath, G.; Brady, J. W.; Venable, R. M.; Pastor, R. W.; Mackerell, A. D. Additive Empirical Force Field for Hexopyranose Monosaccharides. J. Comput. Chem. 2008, 29 (15), 2543–2564. (19) Guvench, O.; Mallajosyula, S. S.; Raman, E. P.; Hatcher, E.; Vanommeslaeghe, K.; Foster, T. J.; Jamison, F. W.; MacKerell, A. D. CHARMM Additive All-Atom Force Field for Carbohydrate Derivatives and Its Utility in Polysaccharide and Carbohydrate– Protein Modeling. J. Chem. Theory Comput. 2011, 7 (10), 3162–3180. (20) Damm, W.; Frontera, A.; Tirado–Rives, J.; Jorgensen, W. L. OPLS All-Atom Force Field for Carbohydrates. J. Comput. Chem. 1997, 18 (16), 1955–1970. (21) Kirschner, K. N.; Yongye, A. B.; Tschampel, S. M.; González-Outeiriño, J.; Daniels, C. R.; Foley, B. L.; Woods, R. J. GLYCAM06: A Generalizable Biomolecular Force Field. Carbohydrates. J. Comput. Chem. 2008, 29 (4), 622–655. (22) Hatcher, E.; Guvench, O.; MacKerell, A. D. CHARMM Additive All-Atom Force Field for Aldopentofuranoses, Methyl-Aldopentofuranosides, and Fructofuranose. J. Phys. Chem. B 2009, 113 (37), 12466–12476. (23) Raman, E. P.; Guvench, O.; MacKerell, A. D. CHARMM Additive All-Atom Force Field for Glycosidic Linkages in Carbohydrates Involving Furanoses. J. Phys. Chem. B 2010, 114 (40), 12981–12994. (24) Wang, X.; Woods, R. J. Insights into Furanose Solution Conformations: Beyond the Two-State Model. J. Biomol. NMR 2016, 64 (4), 291–305. (25) de Leeuw, F. A. A. M.; Altona, C. Computer-Assisted Pseudorotation Analysis of FiveMembered Rings by Means of Proton Spin–spin Coupling Constants: Program PSEUROT. J. Comput. Chem. 1983, 4, 428.
ACS Paragon Plus Environment
40
Page 41 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
(26) Altona, C.; Sundaralingam, M. Conformational Analysis of the Sugar Ring in Nucleosides and Nucleotides.
Improved Method for the Interpretation of Proton
Magnetic Resonance Coupling Constants. J. Am. Chem. Soc. 1973, 95 (7), 2333–2344. (27) Hendrickx, P. M. S.; Corzana, F.; Depraetere, S.; Tourwé, D. A.; Augustyns, K.; Martins, J. C. The use of time-averaged 3JHH restrained molecular dynamics (tar-MD) simulations for the conformational analysis of five-membered ring systems: Methodology and applications. J. Comput. Chem. 2010, 31 (3), 561–572. (28) Taha, H. A.; Castillo, N.; Sears, D. N.; Wasylishen, R. E.; Lowary, T. L.; Roy, P.-N. Conformational Analysis of Arabinofuranosides: Prediction of 3JH,H Using MD Simulations with DFT-Derived Spin−Spin Coupling Profiles. J. Chem. Theory Comput. 2010, 6 (1), 212–222. (29) Taha, H.; Roy, P.-N.; Lowary, T. Theoretical Investigations on the Conformation of the Beta-D-Arabinofuranoside Ring. J. Chem. Theory Comput. 2010, 7, 420–432. (30) Seo, M.; Castillo, N.; Ganzynkowicz, R.; Daniels, C. R.; Woods, R. J.; Lowary, T. L.; Roy, P.-N. Approach for the Simulation and Modeling of Flexible Rings: Application to the α-D-Arabinofuranoside Ring, a Key Constituent of Polysaccharides from Mycobacterium Tuberculosis. J. Chem. Theory Comput. 2008, 4 (1), 184–191. (31) Jana, M.; MacKerell, A. D. CHARMM Drude Polarizable Force Field for Aldopentofuranoses and Methyl-Aldopentofuranosides. J. Phys. Chem. B 2015, 119 (25), 7846–7859. (32) Cremer, D.; Pople, J. A. General Definition of Ring Puckering Coordinates. J. Am. Chem. Soc. 1975, 97 (6), 1354–1358. (33) Wu, A.; Cremer, D. Extension of the Karplus Relationship for NMR Spin-Spin Coupling Constants to Nonplanar Ring Systems: Pseudorotation of Tetrahydrofuran. Int. J. Mol. Sci. 2003, 4 (4), 158–192. (34) Oostenbrink, C.; Villa, A.; Mark, A. E.; van Gunsteren, W. F. A Biomolecular Force Field Based on the Free Enthalpy of Hydration and Solvation: The GROMOS ForceField Parameter Sets 53A5 and 53A6. J. Comput. Chem. 2004, 25 (13), 1656–1676. (35) Gasteiger, J.; Marsili, M. Iterative Partial Equalization of Orbital Electronegativity—a Rapid Access to Atomic Charges. Tetrahedron 1980, 36 (22), 3219–3228. (36) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. DREIDING: A Generic Force Field for Molecular Simulations. J. Phys. Chem. 1990, 94 (26), 8897–8909. (37) Barker, J. A.; Watts, R. O. Monte Carlo Studies of the Dielectric Properties of Waterlike Models. Mol. Phys. 1973, 26 (3), 789–792.
ACS Paragon Plus Environment
41
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 42 of 47
(38) Tironi, I. G.; Sperb, R.; Smith, P. E.; van Gunsteren, W. F. A Generalized Reaction Field Method for Molecular Dynamics Simulations. J. Chem. Phys. 1995, 102 (13), 5451–5459. (39) Berendsen H. J. C., Postma J. P. M., van Gunsteren, WF, Hermans J. Handbook of Intermolecular Forces; Israel: Springer-Science+Business Media BV, 1981. (40) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS:
High
Performance
Molecular
Simulations
through
Multi-Level
Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1–2, 19–25. (41) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New Feathers for an Old Bird. Comput. Phys. Commun. 2014, 185 (2), 604–613. (42) Tironi, I. G.; Gunsteren, W. F. V. A Molecular Dynamics Simulation Study of Chloroform. Mol. Phys. 1994, 83 (2), 381–403. (43) Gee, P. J.; van Gunsteren, W. F. Acetonitrile Revisited: A Molecular Dynamics Study of the Liquid Phase. Mol. Phys. 2006, 104 (3), 477–483. (44) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126 (1), 014101. (45) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52 (12), 7182–7190. (46) Hockney, R. W. Potential Calculation and Some Applications. Methods Comput. Phys. 1970, 9, 135–211. (47) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18 (12), 1463– 1472. (48) Hess, B. P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. J. Chem. Theory Comput. 2008, 4 (1), 116–122. (49) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. In Intermolecular Forces; Pullman, B., Ed.; The Jerusalem Symposia on Quantum Chemistry and Biochemistry; Springer Netherlands, 1981; pp 331–342. (50) Berendsen, H. J. C.; van Gunsteren, W. F.; Zwinderman, H. R. J.; Geurtsen, R. G. Simulations of Proteins in Water. Ann. N. Y. Acad. Sci. 1986, 482 (1), 269–286. (51) Heinz, T. N.; van Gunsteren, W. F.; Hünenberger, P. H. Comparison of Four Methods to Compute the Dielectric Permittivity of Liquids from Molecular Dynamics Simulations. J. Chem. Phys. 2001, 115 (3), 1125–1136.
ACS Paragon Plus Environment
42
Page 43 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
(52) Houseknecht, J. B.; Lowary, T. L.; Hadad, C. M. Improved Karplus Equations for 3JC1,H4 in Aldopentofuranosides: Application to the Conformational Preferences of the Methyl Aldopentofuranosides. J. Phys. Chem. A 2003, 107 (3), 372–378. (53) Haasnoot, C. A. G.; de Leeuw, F. A. A. M.; Altona, C. The Relationship between Proton-Proton NMR Coupling Constants and Substituent electronegativities—I. Tetrahedron 1980, 36 (19), 2783–2792. (54) Altona, C.; Francke, R.; de Haan, R.; Ippel, J. H.; Daalmans, G. J.; Hoekzema, A. J. A. W.; van Wijk, J. Empirical Group Electronegativities for Vicinal NMR Proton–proton Couplings along a C-C Bond: Solvent Effects and Reparameterization of the Haasnoot Equation. Magn. Reson. Chem. 1994, 32 (11), 670–678. (55) Plazinska, A.; Plazinski, W.; Jozwiak, K. Fast, Metadynamics-Based Method for Prediction of the Stereochemistry-Dependent Relative Free Energies of Ligand– receptor Interactions. J. Comput. Chem. 2014, 35 (11), 876–882. (56) Bussi, G.; Gervasio, F. L.; Laio, A.; Parrinello, M. Free-Energy Landscape for β Hairpin Folding from Combined Parallel Tempering and Metadynamics. J. Am. Chem. Soc. 2006, 128 (41), 13435–13441. (57) Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100 (2), 020603. (58) Sugita, Y.; Okamoto, Y. Replica-Exchange Molecular Dynamics Method for Protein Folding. Chem. Phys. Lett. 1999, 314 (1–2), 141–151. (59) Panczyk, K.; Plazinski, W. Pyranose Ring Puckering in Aldopentoses, Ketohexoses and Deoxyaldohexoses. A Molecular Dynamics Study. Carbohydr. Res. 2018, 455, 62–70. (60) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.; et al. Gaussian 09, Revision A.02; Gaussian, Inc.: Wallingford, CT, 2016. (61) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37 (2), 785–789. (62) Vosko, S. H.; Wilk, L.; Nusair, M. Accurate Spin-Dependent Electron Liquid Correlation Energies for Local Spin Density Calculations: A Critical Analysis. Can. J. Phys. 1980, 58 (8), 1200–1211. (63) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98 (45), 11623–11627.
ACS Paragon Plus Environment
43
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 44 of 47
(64) Becke, A. D. Density-functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98 (7), 5648–5652. (65) Hehre, W. J.; Radom, L.; Schleyer, P. von R.; Pople, J. AB INITIO Molecular Orbital Theory, 1 edition.; Wiley-Interscience: New York, 1986. (66) Rayón, V. M.; Sordo, J. A. Pseudorotation Motion in Tetrahydrofuran: An Ab Initio Study. J. Chem. Phys. 2005, 122 (20), 204303. (67) Gaweda, K.; Plazinski, W. The systematic influence of solvent on the conformational features of furanosides. Submitted. (68) Chertkov, A. V.; Pokrovskiy, O. I.; Shestakova, A. K.; Chertkov V. A., V. A. Conformational Analysis of Tetrahydrofuran and Tetrahydrothiopene Using 1H NMR Spectroscopy and Ab Initio Calculations. Chem. Heterocycl. Compd. 2008, 44 (5), 621–623. (69) Lemieux, R. U.; Koto, S.; Voisin, D. The Exo-Anomeric Effect. In Anomeric Effect; ACS Symposium Series; American Chemical Society, 1979; Vol. 87, pp 17–29. (70) Gaweda, K.; Plazinski, W. Pyranose Ring Conformations in Mono- and Oligosaccharides: A Combined MD and DFT Approach. Phys. Chem. Chem. Phys. 2017, 19 (31), 20760–20772. (71) Stenutz, R.; Carmichael, I.; Widmalm, G.; Serianni, A. S. Hydroxymethyl Group Conformation in Saccharides: Structural Dependencies of (2)J(HH), (3)J(HH), and (1)J(CH) Spin-Spin Coupling Constants. J. Org. Chem. 2002, 67 (3), 949–958. (72) Raap, J.; van Boom, J. H.; van Lieshout, H. C.; Haasnoot, A. G. Conformations of Methyl 2’-Deoxy-a-D-Ribofuranoside and Methyl 2’-Deoxy-B-D-Ribofuranoside. A Proton Magnetic Resonance Spectroscopy and Molecular Mechanics Study. J. Am. Chem. Soc. 1988, 110, 2736–2743. (73) Wang, D.; Ámundadóttir, M. L.; van Gunsteren, W. F. van; Hünenberger, P. H. Intramolecular Hydrogen-Bonding in Aqueous Carbohydrates as a Cause or Consequence of Conformational Preferences: A Molecular Dynamics Study of Cellobiose Stereoisomers. Eur. Biophys. J. 2013, 42 (7), 521–537. (74) Plazinski, W.; Drach, M. The Influence of the Hexopyranose Ring Geometry on the Conformation of Glycosidic Linkages Investigated Using Molecular Dynamics Simulations. Carbohydr. Res. 2015, 415, 17–27. (75) S. Serianni, A.; Barker, R. [13C]-Enriched Tetroses and Tetrofuranosides: An Evaluation of the Relationship between NMR Parameters and Furanosyl Ring Conformation. J. Org. Chem. 1984, 49.
ACS Paragon Plus Environment
44
Page 45 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
(76) Church, T. J.; Carmichael, I.; Serianni, A. S. 13C−1H and 13C−13C Spin-Coupling Constants
in
Methyl
β-D-Ribofuranoside
and
Methyl
2-Deoxy-β-D-Erythro-
Pentofuranoside: Correlations with Molecular Structure and Conformation. J. Am. Chem. Soc. 1997, 119 (38), 8946–8964. (77) Cyr, N.; Perlin, A. S. The Conformations of Furanosides. A 13C Nuclear Magnetic Resonance Study. Can. J. Chem. 1979, 57 (18), 2504–2511. (78) AIST:
Spectral
Database
for
Organic
Compounds,
SDBS
https://sdbs.db.aist.go.jp/sdbs/cgi-bin/direct_frame_top.cgi (accessed Jul 31, 2018). (79) Napolitano, J. G.; Gavín, J. A.; García, C.; Norte, M.; Fernández, J. J.; Hernández Daranas, A. On the Configuration of Five-Membered Rings: A Spin–Spin Coupling Constant Approach. Chem. – Eur. J. 2011, 17 (23), 6338–6347. (80) Waterhouse, A. L.; Calub, T. M.; French, A. D. Conformational Analysis of 1-Kestose by Molecular Mechanics and by N.M.R. Spectroscopy. Carbohydr. Res. 1991, 217, 29– 42. (81) Duker, J. M.; Serianni, A. S. (13C)-Substituted Sucrose: 13C-1H and 13C-13C Spin Coupling Constants to Assess Furanose Ring and Glycosidic Bond Conformations in Aqueous Solution. Carbohydr. Res. 1993, 249 (2), 281–303. (82) Weldon, A. J.; Vickrey, T. L.; Tschumper, G. S. Intrinsic Conformational Preferences of Substituted Cyclohexanes and Tetrahydropyrans Evaluated at the CCSD(T) Complete Basis Set Limit: Implications for the Anomeric Effect. J. Phys. Chem. A 2005, 109 (48), 11073–11079. (83) Taha, H. A.; Richards, M. R.; Lowary, T. L. Conformational Analysis of FuranosideContaining Mono- and Oligosaccharides. Chem. Rev. 2013, 113 (3), 1851–1876. (84) Plazinski, W.; Plazinska, A.; Drach, M. Acyclic Forms of Aldohexoses and Ketohexoses in Aqueous and DMSO Solutions: Conformational Features Studied Using Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2016, 18 (14), 9626–9635. (85) Zhu, Y.; Zajicek, J.; Serianni, A. S. Acyclic Forms of [1-13C]Aldohexoses in Aqueous Solution: Quantitation by 13C NMR and Deuterium Isotope Effects on Tautomeric Equilibria. J. Org. Chem. 2001, 66 (19), 6244–6251. (86) Lawrence Que; Jr.; Gary R. Gray. 13C Nuclear Magnetic Resonance Spectra and the Tautomeric Equilibria of Ketohexoses in Solution. Biochemistry 1974, 13 (1). (87) Köpper, S.; Freimund, S. The Composition of Keto Aldoses in Aqueous Solution as Determined by NMR Spectroscopy. Helv. Chim. Acta 2003, 86 (3), 827–843.
ACS Paragon Plus Environment
45
Journal of Chemical Theory and Computation 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 46 of 47
(88) Islam, S. M.; Roy, P.-N. Performance of the SCC-DFTB Model for Description of Five-Membered Ring Carbohydrate Conformations: Comparison to Force Fields, HighLevel Electronic Structure Methods, and Experiment. J. Chem. Theory Comput. 2012, 8 (7), 2412–2423. (89) du Penhoat, C. H.; Imberty, A.; Roques, N.; Michon, V.; Mentech, J.; Descotes, G.; Perez, S. Conformational Behavior of Sucrose and Its Deoxy Analog in Water as Determined by NMR and Molecular Modeling. J. Am. Chem. Soc. 1991, 113 (10), 3720–3727.
ACS Paragon Plus Environment
46
Page 47 of 47 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Graphical abstract
ACS Paragon Plus Environment
47