CARBO_R Force Field for

Article Cite This: J. Phys. Chem. B 2018, 122, 3696−3710

pubs.acs.org/JPCB

Extension of the GROMOS 56a6CARBO/CARBO_R Force Field for Charged, Protonated, and Esterified Uronates Karina Panczyk,† Karolina Gaweda,† Mateusz Drach,‡ and Wojciech Plazinski*,† †

Jerzy Haber Institute of Catalysis and Surface Chemistry, Polish Academy of Sciences, ul. Niezapominajek 8, 30-239 Cracow, Poland Department of Theoretical Chemistry, Faculty of Chemistry, M. Curie-Sklodowska University, pl. M. Curie-Sklodowskiej 3, 20-031 Lublin, Poland

‡

Downloaded via DURHAM UNIV on September 2, 2018 at 07:58:05 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: An extension of the GROMOS 56a6CARBO/CARBO_R force field for hexopyranose-based carbohydrates is presented. The additional parameters describe the conformational properties of uronate residues. The three distinct chemical states of the carboxyl group are considered: deprotonated (negatively charged), protonated (neutral), and esterified (neutral). The parametrization procedure was based on quantum-chemical calculations, and the resulting parameters were tested in the context of (i) flexibility of the pyranose rings under different pH conditions, (ii) conformation of the glycosidic linkage of the (1 → 4)-type for uronates with different chemical states of carboxyl moieties, (iii) conformation of the exocyclic (i.e., carboxylate and lactol) moieties, and (iv) structure of the Ca2+-linked chain−chain complexes of uronates. The presently proposed parameters in combination with the 56a6CARBO/CARBO_R set can be used to describe the naturally occurring polyuronates, composed either of homogeneous (e.g., glucuronans) or heterogeneous (e.g., pectins, alginates) segments. The results of simulations relying on the new set of parameters indicate that the conformation of glycosidic linkage is nearly unaffected by the chemical state of the carboxyl group, in contrary to the ring conformational equilibria. The calculations for the poly(α-D-galacturonate)−Ca2+ and poly(α-L-guluronate)−Ca2+ complexes show that both parallel and anitiparallel arrangements of uronate chains are possible but differ in several structural aspects. Finally, glucuronans13,14 are linear, homogeneous polymers of (1 → 4)-linked β-D-glucuronate (β-D-GlcU) residues. The most relevant processes involving polyuronates occur in the aqueous environment, therefore the experimental determination of structural and conformational properties of these compounds is significantly hindered. In view of this fact, the methods of molecular modeling are a promising alternative offering insight into desired information at the molecular scale.15−21 However, the accuracy of such methods relies primarily on the quality of the molecular models that represent real physical systems. Although modeling the structure/ conformation of uronates as well as their interactions with metal ions has involved a wide spectrum of theoretical methods from rigid docking22 to QM-MM-based simulations,18 the large dimensions of the typical, polyuronate-containing system make the molecular dynamics (MD) simulations performed within the classical force fields the most promising choice. The non-GROMOS carbohydrate-dedicated force fields applied to simulate the uronate molecules include:

1. INTRODUCTION Polyuronates are linear polymers of uronic acids that contain (1 → 4)-linkage (see Figure 1).1 As their molecules bear the carboxyl groups, they are polyelectrolytes able to change the ionization state upon different pH conditions. The three most important natural derivatives of polyuronates are pectins, alginates and glucuronans. They have diverse biological functions in plants, related mainly to the preservation of the structure, texture and flexibility.2−4 Furthermore, they are able to form hydrophilic gels,5−7 typically in the presence of divalent metal cations. Pectins4−9 exhibit a complex sequence involving mainly homogeneous segments of (1 → 4)-linked α-D-galacturonate (α-D-GalU) residues but may also comprise neutral sugar molecules (e.g., rhamnose, galactose, and arabinose) and include branched segments. Furthermore, the carboxyl groups of pectins may be methoxy-esterified to different extents which strongly influences their physicochemical properties.10,11 Alginates3,5−7,12 are linear copolymers of β-D-mannuronate (β-D-ManU) and α-L-guluronate (α-L-GulU) residues in (1 → 4)-linkage that predominantly involve homogeneous segments of either types of residues and regular repeats of their dimer. © 2018 American Chemical Society

Received: November 23, 2017 Revised: March 16, 2018 Published: March 20, 2018 3696

DOI: 10.1021/acs.jpcb.7b11548 J. Phys. Chem. B 2018, 122, 3696−3710

Article

The Journal of Physical Chemistry B

Figure 1. (A) Atom numbering used in the article for the three different chemical states of carboxyl group of uronates. (B) Structural formulas of the uronate residues comprising the natural polyuronates (alginates, pectins and glucuronans) and the monomer of α-L-iduronate. Aliphatic hydrogen atoms (treated implicitly in the GROMOS force field) are not shown.

AMBER14,23 GLYCAM0624 (with the extensions proposed by Dentini et al.25 or Hackbusch et al.26), CHARMM,27 OPLSAA,2 CVFF,28 and (CHARMM-based) CSFF.29 The main field of their application is connected mainly with investigating the structural and conformational properties of naturally occurring mono- (e.g., iduronic acid and its derivatives) and polyuronates (e.g., alginates or pectins) as well as their interactions with divalent metal ions.4,21,30−36 The most recent version of the carbohydrate-dedicated GROMOS force field is the 56a6CARBO set37 which in combination with the correction named 56a6CARBO_R38 describes the nonfunctionalized or O-alkylated pyranoses as well as di-, oligo- and polymers composed of (1 → 1)-, (1 → 2)-, (1 → 3)-, (1 → 4)-, and (1 → 6)-linked residues. So far, no other parameters extending the applicability of this force field in the context of the functionalized carbohydrates have been proposed. Although the GROMOS force field exhibits a large degree of transferability, 38−42 the direct adaptation of parameters for various functional groups from the existing sets may be problematic, especially in reference to the ringinversion properties which are extremely sensitive to alterations in the substituent-substituent and ring-substituent interactions.38 The approach relying on the direct adaptation of the GROMOS-inherent parameters for carboxylate residues has been tested in the past34,35 in the case of the GROMOS 45a4 force field41 but without extensive, quantitative validation against the ab initio or experimental data. In the present article we introduce the parameters that extend the functionality of the GROMOS 56a6CARBO and 56a6CARBO_R sets for the general case of uronates, i.e. aldohexopyranoses substituted by the carboxyl group at position C5 (see Figure 1 for atom numbering). The derived parameters take into account the possible different chemical states of the carboxyl moiety which can be deprotonated, protonated or esterified. Although the parametrization was focused mainly on the most relevant uronates (e.g., α-D-GalU, β-D-GlcU, and α-L-Gul) the final set can be applied to any other uronate residue. Furthermore, because of adaptation of the already existing parameters, the current set can be used to study the longer uronate chains as well as the chains of nonuniform composition, containing either different uronate residues (e.g., alginates) or both uronate and nonuronate residues (e.g., pectins). It is worth mentioning that the direct combination of the existing GROMOS parameters (the 53a6 set to describe the interactions involving carboxyl moieties and the 56a6CARBO/CARBO_R set to describe the rest of carbohydrate

molecule) may not be a viable procedure when considering several different conformational descriptors. More precisely, a preliminary study has shown that such combination leads to inaccurate energy profiles associated with the rotation of a carboxyl group (see Figure S1 in the Supporting Information). Furthermore, the ring-inversion free energies calculated for esterified and protonated uronates are often in disagreement with the available structural data. The newly developed and validated set of parameters is used in order to investigate several issues related to the structure/ conformation of uronate mono- and oligomers as well as to the interactions of natural uronate polymers (alginates and pectins) with calcium ions.

2. METHODS 2.1. Parameterization Strategy. The parameters describing the interactions that do not involve carboxyl moiety were transferred from the GROMOS 56a6CARBO force field37,38 without any changes. This includes the 56a6 CARBO_R correction38 which is relevant for nonmonomeric residues as well as for O1-alkylated residues (see Figure 1 for atom numbering). The parameters present in the original force field and involving the hydroxymethyl moieties were removed and replaced by the new parameters, related to the carboxyl moiety. For the clarity, we do not describe in details every attempt made in order to develop the parameter set but limit ourselves to the final procedure that has led to successful results. The parametrization strategy relied on the following consecutive steps: (i) assigning atomic charges for atoms present in the carboxyl (−COO−), protonated carboxyl (−COOH) of esterified carboxyl (−COOCH3) groups; (ii) assigning the nonbonded (Lennard-Jones) interaction parameters for atoms present in −COO−, −COOH, and −COOCH3; (iii) assigning the bonded (bond stretching, bond angle bending, dihedral angle motion, improper angle bending) interaction parameters involving − COO−, − COOH and − COOCH3; (iv) validation procedure performed by comparison of the energy profiles with the ab initio-originating data and by the analysis of the MD simulation results. Partial atomic charges were assigned on the basis of quantum-mechanical (QM) calculations at the DFT/B3LYP/ 6-311++G** level of theory43−47 and the subsequent RESP48 fitting of the resulting electrostatic potentials. The charge calculations involved α-D-galacturonate and β-D-glucuronate with carboxylate groups in three different chemical states (see above). The RESP procedure involved the constraints on the allowed atomic charges which were introduced in order to (i) 3697


Article


and ζ0, respectively)] are shown in Table 2. These parameters were obtained on the basis of the experimental spectroscopic

maintain the charges on noncarboxyl atoms unchanged with respect to the original set and (ii) keep the charge groups neutral or of integer charge. The final atomic charges are reported in Table 1. The comparison with the charges inherent

Table 2. Bond Stretching, Bond-Angle Bending, and Improper-Dihedral Deformation Parameters for the Carboxyl Group in Uronatesa

Table 1. Atom Types and Atomic Partial Charges for the Carboxyl Group in Uronatesa atom

atom type

bond type

partial atomic charge [e−]

C5−C6

−COO− group C6 O61 O62

C OM OM

C6−O61, C6−O62

+0.284 −0.642 −0.642

C6−O61 C6−O62 O62−H6

−COOH group C6 O61 O62 H62 C6 O61 O62 C7

C OM OA H −COOCH3 group C O OE CH3

+0.874 −0.642 −0.642 +0.410

C6−O61 C6−O62 O62−C7 bond-angle type

+0.874 −0.642 −0.464 +0.232

kb [106 kJ/mol/nm4] 7.15 −COO− group 13.4 −COOH group 16.6 10.2 15.7 −COOCH3 group 16.6 10.2 8.18 kθ [kJ/mol]

−COO− group C5−C6−O61, C5−C6−O62 635.0 O61−C6−O62 770.0 −COOH group C5−C6−O61 685 C5−C6−O62 610 O61−C6−O62 730 C6−O62−H6 450 −COOCH3 group C5−C6−O61 685 C5−C6−O62 610 O61−C6−O62 730 C6−O62−C7 380 improper dihedral type kζ [kJ/mol/deg2]

a

The charge group is created by all atoms belonging to any of the above moieties. See Figure 1 for atom numbering

to the 53a6 set39 reveals that the differences in the case of deprotonated carboxyl groups are rather minor (−0.635 vs −0.642 e and +0.270 vs −0.284 e for the O6 and C6 atoms, respectively). However, both protonation and methylation of the O62 atom lead to a significant charge flux concerning mainly the C6 atom and large divergences from the default parameters present in the 53a6 set.39 For instance, the charges on C6 characteristic of the 53a6 set vary between +0.253 e (methylated carboxyl) and +0.33 e (protonated carboxyl). The current, uronate-oriented calculations predict much higher values (+0.874 e), closer to those inherent to the 53a6OXY set.49 This divergence has a large impact on the predicted energies associated with the rotation of the carboxyl group around the −C5−C6− bond (see Figure S1 in the Supporting Information). The Lenard-Jones (LJ) parameters were unchanged with respect to the GROMOS 53a6/56a6CARBO/56a6CARBO_R sets of parameters. The LJ constants for atoms present in the newly parametrized groups were taken from the 53a6 set and involve OM, OA, OE, CH3, and C atom types. Note that the original GROMOS 56a6CARBO force field contains the “exceptional” LJ pairs which may involve the 1−4 and 1−5 neighbors. In the case of atoms not belonging to the carboxyl group, these type of interactions were kept unchanged with respect to the original set.37,38 The interactions involving C6 were also kept unmodified (the Hassel−Ottar repulsion terms between the C6−O1 and C6−O3 atom pairs as well as substituent-to-ring repulsion terms between the C6−C1 and C6−C3 atom pairs) with respect to 56a6CARBO. Note that this involves the incorporation of the additional atom type in the topological patterns N3 and N4 defined in ref 37 (see Table 5 therein). The potentials for bond stretching, bond-angle bending and improper dihedral deformation terms were retrieved directly from the GROMOS 53a6/56a6CARBO sets. The functional forms of these potentials are given elsewhere37,39−41 and the corresponding parameters [force constants (kb, kθ and kζ, respectively), optimal values of bond lengths and angles (b0, θ0

C6−O5−C4−C5 C6−O61−O62−C5

0.051 0.102

b0 [nm] 0.153 0.125 0.123 0.136 0.1 0.123 0.136 0.143 θ0 [deg] 117.0 126.0 121 115 124 109.5 121 115 124 109.5 ζ0 [deg] 35.264 0

a Functional forms of the corresponding potentials are given in ref37 by eqs 1−3. See Figure 1 for atom numbering.

(force constants) and X-ray diffraction (optimal distance and angle values) data for small molecules. The optimal bondlengths and bond-angle values agree well with the results of our QM calculations (geometry optimization followed by the RESP charge fitting; see above). The torsional potentials, associated with the rotation of the particular elements of −COO−, −COOH and −COOCH3 were adjusted on the basis of the QM-derived values obtained due to relaxed scans of energy in the function of changing values of dihedral angles defined by the following quadruplets of atoms: C4−C5−C6−O61 (anionic β-D-guluronate, α-D-galacturonate, and α-L-guluronate), C5− C6−O62−H62 (protonated α-D-galacturonate and β-D-guluronate) and C5−C6−O62−C7 (esterified α-D-galacturonate). The functional form of the torsional potential is given by the following equation: V (ϕ) = kϕ[1 + cos(mϕ − ϕ0)]

(1)

where ϕ is the dihedral angle value, m the multiplicity of the term, ϕ0 the associated phase shift, and kϕ the corresponding force constant. The values of m, ϕ0, and kϕ were adjusted on the basis of the QM data as shown in Figure 2; the final parameters are given in Table 3. It is worth noting that a given dihedral angle may be involved in more than one torsional potential energy term with different multiplicities and/or phase 3698


Article


Figure 2. Comparison of energy profiles calculated at the classical (the presently proposed parameters) and QM (DFT/B3LYP/6-311++G**) levels of accuracy. Panels A−C present the energy profiles for the rotation of the anionic carboxyl groups of β-D-glucuronate, α-D-galacturonate, and α-Lguluronate, respectively. Panels D and E present the profiles calculated for the protonated carboxyl group, whereas panels F and G for the esterified carboxyl group of α-D-galacturonate, respectively. Panel H shows the profiles for the protonated carboxyl group of β-D-glucuronate. Quadruplets of atoms that define considered dihedral angle are given in each of the panels. The QM- and force field-derived values are given in red and green, respectively.

adjustment of the QM data related to the rotation of carboxyl group is achieved when not assuming the single, 6-fold multiplicity but either the single, 2-fold potential or the combination of the two (1-fold and 2-fold) potentials. The performance of the direct combination of the charges and torsional potentials taken from the 53a6 set with those of 56a6CARBO/CARBO_R was tested in the context of energy profiles associated with the rotation of the carboxyl groups and their agreement with the QM data (analogously to the data shown in Figure 2). The results are shown in Figure S1 (Supporting Information). The QM calculations were performed in Gaussian0950 at the DFT/B3LYP/6-311++G** level of theory. The starting structures were drawn in Avogadro 1.1.151 and preoptimized at the UFF level of theory;52 they reflected the ring conformation characteristic of the series (i.e., 4C1 for Dsaccharides and 1C4 for L-guluronate). The orientation of the carboxyl group was adjusted according to the most favorable

Table 3. Torsional Dihedral Parameters for the Carboxyl Group in Uronatesa torsional dihedral type C4−C5−C6−O61 C4−C5−C6−O61 C4−C5−C6−O61 C5−C6−O62−H62 C4−C5−C6−O61 C4−C5−C6−O61 C5−C6−O62−C7

kϕ [kJ/mol] −COO− group 3.90 −COOH group 3.90 3.50 16.7 −COOCH3 group 3.90 9.50 16.7

ϕ0 [deg]

m

120

2

120 −60 180

2 1 2

120 −60 180

2 1 2

a

The functional form of the corresponding potential is given by eq 1. See Figure 1 for atom numbering.

shifts. Interestingly, and contrary to the default GROMOS parameters for carboxyl group (see, e.g., ref 34) the best 3699


Article


Table 4. Ring-Inversion (4C1 → 1C4) Free Energies (ΔFi) in kJ/mol for the Monomeric Uronates and the Uronate Residues in a Homooctameric Chain, Calculated Using the Presently Proposed Extension to the GROMOS 56a6CARBO Force-Fielda monomers ΔFi compound D-AllU D-AltU D-GlcU D-ManU D-GulU D-IdoU D-GalU D-TalU

O1-Me-L-IdoU O1-Me-D-GlcU a

residue in a chain

α-anomers (−COO−)

β-anomers (−COO−)

α-anomers (−COOH)

β-anomers (−COOH)

−2.7 −13.5 9.6 0.5 5.7 −4.9 16.0 14.1 −4.1 −

18.0 8.1 25.5 16.6 25.3 14.5 32.3 21.9 − 17.1

17.5 8.9 18.7 8.4 17.9 8.6 20.3 12.7 − -

20.2 14.5 22.3 19.5 24.1 16.1 28.3 19.9 − −

compound α-D-GalU α-D-GalUH β-D-GlcU β-D-GlcUH β-D-ManU β-D-ManUH α-L-GulU α-L-GulUH −

ΔFi 18.1 13.2 7.6 1.2 1.4 −1.5 −11.4 −22.7

The calculated values rely on the Cremer−Pople coordinates68 for the state assignment. The part of the data is illustrated graphically in Fig. 3.

The unbiased 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,59 whereas for the constant pressure (1 bar, semi-isotropic coordinate scaling) the Parrinello−Rahman barostat60 was used with a relaxation time of 0.4 ps. The equations of motion were integrated with a time step of 2 fs using the leapfrog scheme.61 The solute bond lengths were constrained by application of the LINCS procedure62,63 with a relative geometric tolerance of 10−4. The full rigidity of the water molecules was enforced by application of the SETTLE procedure.58 The translational center-of-mass motion was removed every time step separately for the solute and the solvent. The nonbonded interactions were calculated using a twin-range scheme,64 with short- and long-range cutoff distances set to 0.8 and 1.4 nm, respectively, and an update frequency of five timesteps for the short-range pairlist and intermediate-range interactions. The reaction-field correction65,66 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 model.67 All the systems were preoptimized by a 0.5−1.5 ns constantpressure 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 200 ns and the trajectory was saved every 10 ps. Apart from the regular GROMACS trajectory, the values of selected parameters were saved every 0.2 ps. These parameters included: the γ dihedral angle, describing the rotation of carboxyl group (defined by the O61−C6−C5−C4 quadruplet of atoms) and the ϕ dihedral angle, describing the rotation of lactol group (defined by the O5−C1−O1−H1 quadruplet of atoms). The distributions of dihedral angles were calculated on the basis of unbiased MD simulations initiated from the structures representing the energetically favorable chair conformer. To study the free energy profiles (FEPs) as a function of the Cremer−Pople θ puckering coordinate68 the parallel tempering69 scheme was combined70 with the well-tempered metadynamics approach71 as implemented in PLUMED 1.3. The well-tempered metadynamics relied on Gaussian local functions of widths 2.86 deg, an initial deposition rate of 0.01

rotamer found during the energy scans. The orientation of remaining groups (exocylic hydroxyl moieties) was adjusted in the positions that correspond to the lack of intramolecular hydrogen bonding (according to the findings described in ref 53). 2.2. Simulation Details. The compounds under consideration included: (i) monomers of all uronates of the D series; (ii) monomers of O1-methylated α-L-iduronate and O1methylated β-D-glucuronate in their deprotonated forms; (iii) homooctamers of (1 → 4)-linked residues of β-D-glucuronate (β-D-GlcU), β-D-mannuronate (β-D-ManU), α-L-guluronate (αL-GulU) and α-D-galacturonate (α-D-GalU); (iv) parallel- and antiparallel-oriented molecular complexes composed of the two oligomers (14 units) of α-L-guluronate (in the anionic state) associated with the contribution of seven calcium ions; (v) parallel- and antiparallel-oriented molecular complexes composed of the two oligomers (14 units) of α-D-galacturonate (in the anionic state) associated with contribution of seven calcium ions. For compounds listed in points i and iii, both their anionic and protonated (the latter denoted by subscript “H”) forms were taken into account; additionally, for α-D-GalU-containing compounds, the esterified forms were considered (denoted by subscript “Me”). Figure 1 shows the chemical structures of the studied compounds or their monomeric units. The initial structures of complexes mentioned in points iv and v were based on the results from our previous studies15 and on ref 22. All MD simulations were carried out with the GROMACS 4.5.5 package 54 in combination with PLUMED 1.3 55 (metadynamics) or GROMACS 5.056 with PLUMED 2.1.557 (unbiased MD). The parameters used for simulations are described in the previous section. Molecular systems under consideration consisted of one carbohydrate molecule and approximately 900−1600 simple point charge (SPC) water molecules58 placed in a cubic (30 × 30 × 30 Å3 or 36 × 36 × 36 Å3) simulation box. The systems of smaller or larger dimensions corresponded to simulations of monomers or homooctamers, respectively. Additionally, we studied the complexes composed of the two uronate chains and calcium ions. The initial geometry of such complexes was taken from our previous studies15 and the number of water molecules (about 19160) as well as the edges of simulation cell (97 × 97 × 97 Å3) were increased to account for the larger dimensions of related systems. 3700


Article


IdoU, α-D-AllU, α-D-ManU and α-D-GulU confirm the nonnegligible populations of both 1C4 and 4C1 conformers. The protonation of the carboxyl moiety leads to a notable increase of the ring-inversion energies; this increase is the largest for those α-anomers which exhibited the lowest values of ΔFi in the form of deprotonated compounds. This effect results in positive ΔFi values for all neutral uronates and significant reduction of the population of non-4C1 conformers for the whole group of uronates as the lowest ΔFi value equals now +8.4 kJ/mol. Therefore, our results indicate that pH-dependent protonation state may have an influence on the ring-inversion properties of uronates but this effect is limited to the α-anomers of compounds which exhibit most flexible rings. However, this group contains several among the most biologically relevant and naturally abundant uronates (e.g., iduronate, mannuronate and guluronate). Interestingly, the population of the nonchair conformers (i.e., boat and skew-boat ring geometries) is always much smaller in comparison to both the 4C1 and 1C4 ones and such effect is independent of the protonation state. This means that such ring geometries do not contribute much into the overall conformational space of uronate monomers and the ring conformational equilibria in some very flexible uronates (e.g., IdoU and its derivatives)75−77 can be described by taking into account only the 4C1 and 1C4 shapes. The similar findings have been reported in refs 78 and 79 for unfunctionalized hexopyranoses and hexopyranose residues in a chain. When comparing the ΔFi values calculated for uronate monomers with the analogous quantities obtained for corresponding unfunctionalized aldohexopyranose monomers (data taken from ref.38) one obtains a clear correlation between these two sets of data, especially in the case of protonated uronates. This correlation is slightly less for deprotonated uronates which (mostly in the case of the α-anomers) exhibit the lower ΔFi values. The existence of such correlation between chair−chair inversion energies can be expected on rational grounds, mainly due to the additivity of the interactions between ring substituents. This effect has been confirmed experimentally and became a foundation for proposing semiempirical schemes capable of predicting the ΔFi values for various pyranoses.80−82 A similar additivity effect, resulting in a correlation between ΔFi values, has been confirmed to occur for various groups of pyranoses.83 Additionally, we have calculated the ring-inversion properties for residues in a chain for selected types of homooctamers. The results are given in Table 4. Although in the majority of the cases the favorable ring shape is the same as for the corresponding monomer, some of residues in a chain are very prone to ring inversions (e.g., mannuronates and protonated glucuronate). This is also the case of some of unfunctionalized hexopyranoses in a chain (see Table 3 in ref 38) and seems to be an inherent feature of the GROMOS 56a6CARBO_R force field but is not necessarily in disagreement with the real conformational behavior of the polyuronates in aqueous solutions (note that there is no experimental data regarding the ring shapes in polyuronate molecules in solutions). Interestingly, the influence of the protonation state seems to be opposite to that observed in the case of monomeric uronates (Figure 3) but due to the fact of limited number of considered oligomers, this may be ascribed to the decorrelation of the ring-inversion properties of residues in a chain relative to monomers. Similarly to the results described elsewhere,38,53 also here we have not been able to

kJ/mol/ps and a temperature parameter ΔT (see eq 2 in ref 71) of 1788 K. The parallel-tempering69 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. For the FEP analyses along θ, the intervals (0−60 deg) and (120−180 deg) were ascribed to regular chair (4C1) and inverted chair (1C4) conformations, respectively. The freeenergy difference associated with the chair−chair (4C1 → 1C4) distortion was defined as the difference between the values at the minima of FEP in the two corresponding regions and denoted as ΔFi. The convergence of the final ΔFi values was checked by the hand-written bash script and the PLUMED inbuild sum_hills tool. The convergence plots are given in Figure S2 (Supporting Information). The conformation of glycosidic linkages was studied in terms of 2D free energy maps (FEMs) calculated by using the combination of parallel tempering and well-tempered metadynamics, analogously to the procedure described above. The φ and ψ dihedral angles defining glycosidic linkages are created by the O5−C1−O1−C′4 and C1−O1−C′4−C′3 quadruplets of atoms, respectively. The simulations relied on Gaussian local functions of widths 18 deg and an initial deposition rate of 0.01 kJ/mol/ps. The 3JH,H coupling constants for vicinal protons were calculated by using the Altona−Haasnoot relation72−74 and the unbiased MD trajectories initiated from both the 1C4 and 4 C1 conformers (Me−O1-α-L-IdoU and α-L-IdoU) or only from the 4C1 one (Me−O1-β-D-GlcU). The average 3JH,H value for Me−O1-α-L-IdoU and α-L-IdoU was calculated as weighted average of both 4C1 and 1C4 contributions with the weight assigned on the basis of the separately calculated ΔFi value.

3. RESULTS AND DISCUSSION 3.1. Flexibility of the Pyranose Rings. The chair−chair ring-inversion free energies (ΔFi) are given in Table 4 while some of their values calculated for monomers are graphically illustrated in Figure 3.

Figure 3. Ring-inversion free energies (ΔFi) calculated within the current set of parameters for both deprotonated and protonated uronates, compared to the corresponding values for unfunctionalized aldohexopyranoses (data taken from ref 38).

Regarding the ring-inversion properties of the deprotonated uronate monomers, the sign of the calculated ΔFi values shows that the dominant conformation for 13 (out of 16) compounds is a regular chair. The remaining three (α-D-AllU, α-D-IdoU, and α-D-AltU) exhibit the inverted chair geometry of the ring. The close-to-zero values of ΔFi (|ΔFi| < 5 kJ/mol) for α-D3701


Article


theoretical proportion of the 4C1 and 1C4 conformers that matches the experimental results best. The minimal deviation between the theoretically- (via. Haasnoot−Altona relation72) and experimentally derived 3JH,H values can be achieved by assuming the ΔFi value equal to about +1.1 kJ/mol which corresponds to the 39%:61% (1C4: 4C1) equilibrium mixture of chair conformers. This is in a fair agreement with the predictions of the current set of parameters (ΔFi = −4.1 kJ/ mol for O1-methylated and +4.9 kJ/mol for unfunctionalized αL-IdoU) in spite that the signs of the calculated ring-inversion energies may be opposite. However, the deviations from the experimental data are relatively minor (up to 8 Hz for the whole set of data) if accepting any ΔFi value from the range (−4.3 kJ/mol; +5 kJ/mol). This speaks for the non-negligible populations of both regular and inverted chair conformers and significant conformational heterogeneity of the α-L-IdoU molecule, in contrary to the β-D-GlcU one. 3.2. Orientation of the Exocyclic Groups. There is no direct experimental data related to the orientation of the lactol groups in uronate saccharides. However, the conformational distributions can be checked for agreement with the trends expected considering the exo-anomeric effect.84 Regarding the D-pyranose molecule in the regular chair conformation, 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 be essentially absent. These expectations are reflected very well by the results of the MD simulations shown in Figure 5 when considering GlcU and GalU monomers in

distinguish any single interaction pattern responsible for this effect and its magnitude. To our knowledge, no quantitative data exist on the ringinversion properties of any uronate monomers except of O1methylated β-D-GlcU and α-L-IdoU.75 Ring three-bond H−H vicinal coupling constant (3JH,H) measurements revealed significant differences between the ring geometries of these two compounds. The larger values of 3JH,H for β-D-GlcU speak for the strongly dominating 4C1 conformation whereas the lower ones for α-L-IdoU suggest the significant contribution of non-4C1 shapes. In ref 75, the latter finding has been interpreted as the existence of the equilibrium mixture of chair and skew-boat conformers. The results of our simulations offer an alternative interpretation based rather on the dynamic equilibrium between 1C4 and 4C1 conformers with the negligible contributions of any boat or skew-boat conformers. Figure 4 shows the experimental data compared with the results

Figure 4. Experimental (yellow bars) and calculated vicinal coupling constants for O1-methylated α-L-iduronate and β-D-glucuronate (red bars) as well as for unfunctionalized α-L-iduronate (blue bars). The experimental data were taken from ref 75 whereas the Haasnoot− Altona equation72 was used to obtain the theoretical values.

of our MD simulations. The ΔFi value calculated for Me−O1-αL-IdoU is equal to −4.1 kJ/mol which suggests the domination of the 1C4 conformation with a non-negligible (∼16%) contribution of the 4C1 conformer. The chair-boat/skew-boat inversion is associated with much higher free energy (>10 kJ/ mol) and the corresponding contribution is neglected. The agreement with the experimental data is almost ideal in the case of β-D-GlcU and satisfactory in the case of α-L-IdoU. In the later case, the agreement is of comparable quality in comparison with the predictions of the GLYCAM force field24 according to Figure 2 in ref.75 According to the predictions of the GROMOS 56a6CARBO_R force field, methylation at O1 greatly influences the ringinversion properties of pyranose monomers in comparison to the unfunctionalized ones. More precisely, the Me−O1 substitution results in a systematic shift (∼7.5 kJ/mol) of ring-inversion free energies (4C1 → 1C4) which is positive for the α and negative for β-anomers of D-hexopyranoses (the trend is opposite for L-hexopyranoses). However, recent studies53 indicate that the magnitude of such shift may be smaller than previously estimated. Because of this observation, we additionally calculated the 3JH,H constants for unfunctionalized α-L-IdoU which exhibits the ΔFi value equal to ∼5 kJ/ mol. The results given in Figure 4 are also in a good agreement with the experimental data (although the direction of deviations is opposite in comparison to the results obtained for Me−O1-αL-IdoU). On the basis of the above experimental 3JH,H values related to α-L-IdoU and the two-state model one can calculate the

Figure 5. Probability distribution plots related to the value of the ϕ dihedral angle that describes the conformation of the lactol group (defined by the O5−C1−O1−H1 quadruplet of atoms). All the plots correspond to the unbiased MD simulations in the explicit water solvent.

different protonation states. The agreement with the analogous data collected for unfunctionalized aldohexopyranoses (see, e.g., Figure 11 in ref 37) indicates that the influence of the carboxyl moiety on the orientation of the lactol group is very limited; the same is true in the case of the influence of the orientation of the hydroxyl group at C4 (which results in minor differences between results obtained for GlcU and GalU). This is additionally confirmed by the fact that also the protonation state is not relevant for the conformation of the lactol moiety. The deprotonated carboxyl groups of both GlcU and GalU exhibit a very high degree of conformational flexibility (Figure 6). The anomery of the saccharide is irrelevant for the conformation of the −COO− group, in contrast to the orientation of the substituent at C4, which, when oriented axially, reduces the rotational freedom of −COO−. However, in both possible cases, changing the value of the γ angle is associated with overcoming very low free energy barriers: ∼5 3702


Article


Figure 6. Probability distribution plots related to the value of the γ dihedral angle that describes the conformation of the carboxyl group (defined by the C4−C5−O6−O61 quadruplet of atoms). All the plots correspond to the unbiased MD simulations in the explicit water solvent.

Figure 7. Values of the γ-dihedral angle, describing the orientation of the carboxyl group in uronates, according to the XRD data deposited in the PDB database.

Although the X-ray diffraction (XRD) structures cannot be treated in terms of a direct evidence for the conformational behavior of saccharide molecules in solutions, they still may provide valuable qualitative information about the conformational preferences of interest which includes the insight into the most favorable geometries or conformationally restricted regions. Figure 7 contains the XRD data extracted from the PDB database and represented by the value of the γ dihedral angle. The comparison with the results of MD simulations (Figure 6) reveals a good agreement for a series of naturally abundant D-uronates, i.e. α-GalU, β-GlcU, and β-ManU, as well as for α-L-GulU. For compounds with an equatorially oriented hydroxyl group at C4 (with respect to the ring in the 4C1 conformation, i.e. β-GlcU and β-Man) the available conformational space covers all the possible values of the γ dihedral angle. Upon changing the mutual orientation of the mentioned hydroxyl group and carboxyl group (i.e., for α-L-GulU and α-DGalU), a much clearer preference for conformations around γ = −60/+120 deg is visible whereas some of the conformations are excluded (i.e., around γ = 0/ ± 180 deg). This remains in agreement with the course of the probability distribution plots shown in Figure 6. To our knowledge, there is no structural data on acidic uronates containing the explicitly defined location of carboxylate-bound proton, to perform an analogous comparison with the MD results. Note that the use of the default, carboxylate-related GROMOS parameters results in predicting a notably different rotational behavior of the carboxyl moiety; see Figure 3 in ref 34 for the quantitative results. In particular, the free energy barriers separating particular rotamers are much higher than those predicted in the present study. Moreover, in contrary to the results presented in Figure 6, the compounds with a possible syn orientation of both O4 and O61/O62 (i.e., GalU and GulU) exhibit the higher rotational flexibility of the carboxyl group in comparison to e.g. GlcU and ManU.34

kJ/mol (for GalU) and 2 kJ/mol (for GlcU). In the case of GalU, the least favorable conformations (γ = ∼60 deg and γ = ∼−120 deg) are associated with the close proximity of any of the carboxyl oxygen atoms and O4 (∼0.26 nm) when they are in the syn position. On the contrary, the most favorable conformations (γ = ∼145 deg and γ = ∼−35 deg) correspond to perpendicular orientation of the carboxyl group plane with respect to the C5−H covalent bond. In the case of GlcU, the most favorable orientations (γ = ∼30 deg and γ = ∼−150 deg) are represented by the perpendicular orientation of the carboxyl group plane with respect to the C5−O5 covalent bond whereas the least favorable ones by the perpendicular orientation of the −COO− plane with respect to the C5−C4 covalent bond. The protonation of the carboxyl group is associated with a significant reduction of the available conformational space. For each of the compounds only two orientations exhibit the nonnegligible populations, i.e., the γ = ∼70 deg and γ = ∼130 deg values which correspond to primary rotamers for GlcU and GalU, respectively, whereas γ = ∼−140 deg and γ = ∼−40 deg represent the secondary conformers. Again, the anomery is not relevant, in contrary to the orientation of the hydroxyl group at C4. The γ values characterizing the most populated rotamers are very close to those observed for deprotonated moieties of the corresponding compounds. However, due to the loss of internal symmetry of the carboxyl group (i.e., the attachment of proton to O62), the resulting asymmetry in distribution of rotamers is observed. Interestingly, the favorable orientations of the protonated carboxyl group are those that permit the formation of hydrogen bonding with hydroxyl group at C4 (which can be either axially or equatorially oriented as reflected by different positions of dominating peaks in Figure 6). Although this effect is not necessarily imposed by the presence of such interactions,85,86 the behavior opposite to that typical for unfunctionalized hexopyranoses and hydroxymethyl group is worth emphasizing. 3703


Article


Figure 8. Free-energy maps in the space of the glycosidic φ and ψ dihedral angles for nine different octamers composed on uronate residues exhibiting different protonation/esterification states, calculated using the presently proposed extension of the GROMOS 56a6CARBO force field in the presence of explicit water. The contours are plotted every 20 kJ/mol.

3.3. Orientation of the Glycosidic Linkages. The conformation of longer, oligo- or polymeric saccharide chains is controlled mainly by the orientation of glycosidic linkages between residues. According to the results described elsewhere,19 in the case of uronate residues, the conformation of glycosidic linkage characteristic of dimers can be slightly different than that of longer chains. Therefore, in the present study we focused on the longer, octameric chains which better reflect the conformational behavior of natural polyuronate molecules. Figure 8 shows the calculated free energy maps relying on the glycosidic φ and ψ dihedral angles. The results indicate that the chemical state of the carboxyl group is not relevant in the context of the conformation of longer saccharide chains. Both the location of the main minima of the free energy and the dimension of the available conformational phase space remain nearly unchanged, independently if the carboxyl group is protonated, deprotonated or esterified. Moreover, the favorable conformations are very close to those characteristic of the corresponding unfunctionalized hexopyranoses37 which suggests that the

type of exocyclic substituent is not crucial for the dynamic equilibrium related to the chain shape. Likeliness of free energy maps calculated for ManU and GlcU shows than the primary factor controlling the orientation of glycosidic linkage is the disposition of glycosidic oxygen atoms with respect to the ring of the given residue (i.e., axial vs equatorial). The favorable conformation of nonmonomeric saccharide connected by the given type of glycosidic linkage always corresponds to the same region of φ vs. ψ phase space. This is in agreement with the results of other studies, performed for unfunctionalized oligohexopyranoses78,87 and disaccharides.86 Furthermore, this effect is not disturbed much upon the inversion of the ring from 4 C1 to 1C4 as indicated by extremely similar maps obtained for α-L-GulU (diaxially linked 1C4 rings) and β-D-GlcU (diequatorially linked 4C1 rings). This observation is also common for unfunctionalized aldohexoses.87 Figure S3 (Supporting Information) contains the glycosidic φ and ψ dihedral angle values extracted from the XRD data measured for pectins and mannuronate-rich alginates. The comparison with the free-energy maps shown in Figure 8 3704


Article


Figure 9. Probability distribution plots related to the values of the dihedral or regular angles describing: (i) the conformation of the glycosidic linkage (φ and ψ); (ii) the helical twist of the single chain either involved or not in formation of the complex (χ); (iii) the conformation of the carboxyl group (γ). Additionally, the regular α angle (defined by the C2−Ca2+−C2 triplet where the C2 atoms belong to the two residues connected by Ca2+ ion) describing the mutual orientation of the two chains that create the complex, is considered. The schematic plot shows the differences between both possible chain−chain arrangements. All the plots correspond to the unbiased MD simulations of the α-L-GulU−Ca2+ complexes in the explicit water solvent.

Figure 10. Probability distribution plots related to the values of the dihedral or regular angles related to (i) the glycosidic linkage (ϕ and ψ), (ii) the helical twist of the uronate chain (χ), (iii) the conformation of the carboxyl group (γ), and (iv) the mutual orientation of the two chains creating the complex (α). All the plots correspond to the unbiased MD simulations of the α-D-GalU−Ca2+ complexes in the explicit water solvent. Other details as in Figure 9.

introduction of the Project’s correction90 in the default LJ parameters for the CA2+ and OM atom types is necessary. Otherwise, only the “loose” binding, contradictory to the experimental observations5,6,10 can be observed.35 The analogous correction is required also in the case of other biomolecular force fields.90 The ability for binding calcium ions and creating the stable complexes composed of two polyuronate chains and calcium

reveals an excellent agreement with respect to the location of the most favorable conformations for both β-D-ManU and α-DGalU following the trends expected on the basis of the stereoelectronic effects governing the orientation of glycosidic linkages.88,89 3.4. Chain−Chain Complexes. In order to obtain the correct Ca2+−carboxylate binding mode (i.e., the so-called “tight” binding) and stable calcium-uronate complexes, the 3705


Article

The Journal of Physical Chemistry B ions was validated in the context of stability of complexes created by the two chains of α-L-GulU and calcium ions. The geometry of such complexes is close to the so-called “egg-box” structural model5,6,10,32,91 and is described in details elsewhere.15−17 The simulations performed with use of the current version of parameters indicate that (i) the Ca2+-guluronate complex is stable during the 200 ns-long MD run; (ii) all the basic structural parameters (i.e., coordination of the Ca2+ ions, conformation of the glycosidic linkages, helical twist of the whole complex, etc.) are very close to those reported in refs 15 and 16 and are based on the application of the GROMOS 45a4 force field. This is understandable in view of the close similarities between the 45a4 and 56a6CARBO sets related to the treatment of the glycosidic linkages; these similarities are described in ref 37. The overall structure is different from the classical egg-box model91,92 and has been confirmed by the more accurate, QM-MM/MD simulations.18 As the present paper contains one of the first descriptions of the MD-based structure of the pectin−Ca2+ complexes,32 let us summarize briefly the most important structural parameters found in our investigation. At the same time, let us note that a more detailed study devoted to the conformational properties of pectins, including structural heterogeneity effects (i.e., esterification and/or presence of non-α-D-GalU units) is postponed to our future work. In the case of both parallel and antiparallel arrangements of the two Ca2+-linked α-D-GalU chains, we observed the interaction pattern similar to that observed in the Ca2+− guluronate complexes.15 The calcium ions are coordinated by eight oxygen atoms; the four of them belong to the two carboxyl moieties of the two paired chain whereas the remaining four to the water molecules. No other atoms contribute to the direct contact with the calcium ions. The antiparallel arrangement of the α-D-GalU chains results in the close proximity of the C2- and C3-linked hydroxyl groups belonging to two different chains. This pattern involves both Ca2+-binding residues and the remaining ones. Some of these proximities are associated with interchain hydrogen bonding. On the contrary, the parallel arrangement does not favor any direct contacts between chains and the whole structure is more loose than in the case of the antiparallel arrangement. In addition to the results described in ref 15, let us summarize some further findings originating from the current simulations. Some of them are graphically shown in Figures 9 and 10, whereas Figure 11 presents the simulation snapshots illustrating typical structural properties of the investigated complexes. Regarding α- L -GulU-containing complexes, although both the parallel and antiparallel binding modes appear to be stable over the time scale of hundreds of nanoseconds, there exist several differences between them. First, the conformational behavior considered in the context of single oligomeric chain may change upon complexation by the second chain and calcium ions; the degree of such change is dependent on the type of chain−chain pairing. The antiparallel arrangements of both α-L-GulU and α-D-GalU result in larger changes of the following conformational descriptors: (i) the magnitude of the helical twist of single chain (χ, defined by the C6−C2−C′2−C′6 quadruplets of atoms); (ii) the rotation of the carboxyl moieties engaged in calcium binding. Although the average twist exhibited by the α-L-GulU chains is roughly the same, independently of their molecular environment (χ ≈ −30 deg), the flexibility of the chain is significantly reduced upon anitiparallel binding. On the contrary, parallel arrangement

Figure 11. Compilation of snapshots extracted from the unbiased MD simulations performed within the current set of force field parameters. Both view along the complex axis (right) and the central fragments of complexes (lef t) are shown. Calcium ions are represented as yellow balls. Water molecules (present during simulations) are not shown for clarity.

does not significantly influence the helical twist. The same observation is also valid for the systems containing α-D-GalU chains. The only difference are the average twists which differ depending if either single chain (χ ≈ −30 deg) or parallel (χ ≈ 0 deg) and antiparallel (χ ≈ 0 deg) complexes are considered. Interestingly, the changes in the favorable χ value are not associated with the corresponding divergences in the glycosidic dihedral angle values which is in agreement with the results reported in ref 19. Note that the complexation of both α-LGulU and α-D-GalU chains may result in a small reduction of their flexibilities expressed in terms of the φ and ψ values. Furthermore, the antiparallel binding mode results in slightly modified pattern of rotational preferences of those carboxyl groups which contribute to the Ca2+-mediated chain−chain binding. Namely, the preferred conformation is shifted toward higher values of the γ angle (which results in the favorable eclipsed conformation with respect to the C4 and O6 atoms). Additionally, in the case of α-D-GalU, the rotational freedom of the carboxyl moiety is slightly reduced. The diverse results obtained for galactouronate- and guluronate-containing complexes in the context of their helical twist may contribute to the mechanistic interpretation of different gelling behaviors observed experimentally for pectins and alginates. A more linear form of chain−chain complexes typical for galactouronates may favor the association of further chains (also linked by the contribution of the Ca2+ ion) and 3706


Article


parts of uronates. The validation tests proved the agreement with the external data in the context of glycosidic linkage conformation, rotameric freedom of carboxyl group, ringinversion properties and the stability of the complexes formed by the Ca2+ ions and polyuronate chains. The results obtained for α-L-iduronate suggest that the conformational equilibria of the ring shapes is dominated by the chair conformers and the two-state model indicates that the composition of the equilibrium mixture is close to 39%:61% (1C4:4C1). The protonation state of the carboxyl moiety may influence the flexibility of the pyranose ring; the anionic monomers are usually more prone to the ring inversion in comparison to their protonated analogues and corresponding unfunctionalized aldohexopyranoses. The anionic carboxyl group is very prone to rotation, with the corresponding free energy barriers varying from ∼2 to 5 kJ/mol. The main factor controlling its conformational preferences is the disposition of the hydroxyl group at the C4 atom. The protonation of the carboxyl group results in a significant reduction of the rotational freedom. Conversely, the chemical state of the carboxyl group has no influence on the conformation of the glycosidic linkages in uronate oligomers. Regarding the Ca2+−uronate complexes, both α-L-GulU and α-D-GalU are able to create stable structures according to either parallel or antiparallel arrangements of the chains. The pattern of interactions with calcium ions is the same as that reported in ref 15; i.e., only the calcium carboxylate contacts contribute to the tight ion binding. The antiparallel arrangements of both α-L-GulU and α-D-GalU chains result in more compact structures with closer contacts between chains; on the contrary, parallel arrangements are associated with more loose structures and less intensive chain− chain contacts. Differences between guluronate and galacturonate complexes suggest the more facile contacts between the already paired chains in the latter case.

stabilize the junction zone. The larger helical twist observed in the case of guluronates results in a relatively smaller area of possible contacts between two (or more) complexes.16 Nevertheless, due to a more heterogeneous structure of naturally occurring pectins, a more detailed study is required in this matter. The differences between parallel and antiparallel arrangements of chains can be ascribed to the molecular details of the interactions between chains. Independently of the considered compounds, the antiparallelly oriented chains exhibit much closer contact with each other in comparison to the parallel arrangements. This is expressed by the value of the α angle which describes the orientation between the planes of the two chains in the manner depicted symbolically in Figures 9 and 10 (the regular α angle is defined by the C2−Ca2+−C2 triplet of atoms where the C2 atoms belong to the two residues connected by the same Ca2+ ion). Strictly distinct distributions of this angle value speaks for much tighter steric arrangement observed in antiparallel complexes which is correlated with a more intensive network of chain−chain hydrogen bonds and close distance between those edges of the chains which are not involved in the Ca2+-mediated chain−chain binding. These stronger chain−chain interactions result in alterations of the conformational degrees of freedom, discussed above. On the contrary, the parallel arrangement results in a more loose structure of the whole complex, demonstrated by the values of conformation-related parameters close to those characteristic of single chains. Apart from using the distributions of the α angle values (Figures 9 and 10), the same effect can be illustrated by calculating the radial distribution functions (RDFs) for the distance between the two residues bound by the same calcium ion. The results of such calculations are graphically presented in Figure S4 (Supporting Information). Finally, let us note that a more loose structure of parallelly arranged complexes is reflected by the location of the bound calcium ions. In the case of antiparallel chain−chain orientation, all bound Ca2+ ions are located in one straight row, i.e. the angle between any consecutive three of them is close to 180 deg. One the contrary, in the case of parallel complexes, the conformational rearrangements, occurring in a relatively large time scale (i.e., tens to hundreds of nanoseconds) can be seen, which correspond to local deformations of the complex and shift some of the calcium ions out of the ideal straight line. Such behavior was found to be independent of the initial structure and occurred multiple times during the single 200-ns-long simulation; thus, we concluded that it is an inherent feature of the parallel arrangement described in terms of the current force field parameters.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.7b11548. Plots showing the performance of the parameters adopted from the 53a6 set in comparison with the QM data, convergence plots related to the ring-inversion free energies, RDFs showing different chain−chain contact patterns for both parallel and antiparallel arrangements of chain−chain complexes, and the PDB-extracted data showing distribution of glycosidic dihedral angles for various uronates. (PDF)

■

4. CONCLUSIONS The article presents an extension of the GROMOS 56a6CARBO/CARBO_R force field aimed at description of the uronate mono-, di-, oligo-, and polysaccharides. The derived set of parameters takes into account the different chemical states of carboxylate groups (anionic vs protonated vs esterified). In combination with the previously proposed sets, the present extension can be used to study the conformational features of all naturally occurring uronate molecules including those composed of heterogeneous uronate residues (e.g., alginates) or incorporating nonuronate components (e.g., pectins). The parametrization procedure was based both on the quantummechanical calculations (in the context of carboxyl groups) and incorporating the already existing parameters for noncarboxyl

AUTHOR INFORMATION

Corresponding Author

*(W.P.) Telephone: +48815375685. Fax: +48815375685 Email: [email protected]; [email protected]. ORCID

Wojciech Plazinski: 0000-0003-1427-8188 Notes

The authors declare no competing financial interest.

■

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). 3707


Article


(22) Braccini, I.; Pérez, S. Molecular Basis of Ca2+-Induced Gelation in Alginates and Pectins: The Egg-Box Model Revisited. Biomacromolecules 2001, 2, 1089−1096. (23) The Amber Molecular Dynamics Package. http://ambermd. org/ (accessed Feb 6, 2018). (24) Kirschner, K. N.; Yongye, A. B.; Tschampel, S. M.; GonzálezOuteiriño, J.; Daniels, C. R.; Foley, B. L.; Woods, R. J. GLYCAM06: A Generalizable Biomolecular Force Field. Carbohydrates. J. Comput. Chem. 2008, 29, 622−655. (25) Dentini, M.; Rinaldi, G.; Barbetta, A.; Risica, D.; Anselmi, C.; Skjåk-Bræk, G. Ionic Gel Formation of a (Pseudo)alginate Characterised by an Alternating MG Sequence Produced by Epimerising Mannuronan with AlgE4. Carbohydr. Polym. 2007, 67, 465−473. (26) Arrate, M.; Aurrecoechea, J. M. Synthesis of Bicyclic Alcohols by Palladium-Catalyzed Et2Zn-Mediated Intramolecular Carbonylpropargylation. ARKIVOC 2017, 2017, 257−267. (27) Guvench, O.; Mallajosyula, S. S.; Raman, E. P.; Hatcher, E.; Vanommeslaeghe, K.; Foster, T. J.; Jamison, F. W., II; MacKerell, A. D., Jr. CHARMM Additive All-Atom Force Field for Carbohydrate Derivatives and Its Utility in Polysaccharide and Carbohydrate− Protein Modeling. J. Chem. Theory Comput. 2011, 7, 3162−3180. (28) Dauber-Osguthorpe, P.; Roberts, V. A.; Osguthorpe, D. J.; Wolff, J.; Genest, M.; Hagler, A. T. Structure and Energetics of Ligand Binding to Proteins: Escherichia Coli Dihydrofolate ReductaseTrimethoprim, a Drug-Receptor System. Proteins: Struct., Funct., Genet. 1988, 4, 31−47. (29) Kuttel, M.; Brady, J. W.; Naidoo, K. J. Carbohydrate Solution Simulations: Producing a Force Field with Experimentally Consistent Primary Alcohol Rotational Frequencies and Populations. J. Comput. Chem. 2002, 23, 1236−1243. (30) Liu, Y.; Chipot, C.; Shao, X.; Cai, W. Solubilizing Carbon Nanotubes through Noncovalent Functionalization. Insight from the Reversible Wrapping of Alginic Acid around a Single-Walled Carbon Nanotube. J. Phys. Chem. B 2010, 114, 5783−5789. (31) Babin, V.; Sagui, C. Conformational Free Energies of Methyl-αL-Iduronic and Methyl-β-D-Glucuronic Acids in Water. J. Chem. Phys. 2010, 132, 104108. (32) Assifaoui, A.; Lerbret, A.; Uyen, H. T. D.; Neiers, F.; Chambin, O.; Loupiac, C.; Cousin, F. Structural Behaviour Differences in Low Methoxy Pectin Solutions in the Presence of Divalent Cations (Ca(2+) and Zn(2+)): A Process Driven by the Binding Mechanism of the Cation with the Galacturonate Unit. Soft Matter 2015, 11, 551− 560. (33) Labus, K.; Radosiński, Ł.; Kuchta, B.; Trusek-Hołownia, A. Temperature Effect on Properties of β-Galactosidase Entrapped in Alginate Matrix: An Experimental Research Supported by Molecular Modeling. Inż. Apar. Chem. 2015, 3, 98−100. (34) Perić, L.; Pereira, C. S.; Pérez, S.; Hünenberger, P. H. Conformation, Dynamics and Ion-Binding Properties of Single-Chain Polyuronates: A Molecular Dynamics Study. Mol. Simul. 2008, 34, 421−446. (35) Perić-Hassler, L.; Hünenberger, P. H. Interaction of Alginate Single-Chain Polyguluronate Segments with Mono- and Divalent Metal Cations: A Comparative Molecular Dynamics Study. Mol. Simul. 2010, 36, 778−795. (36) Anderson, R. L.; Greenwell, H. C.; Suter, J. L.; Coveney, P. V.; Thyveetil, M.-A. Determining Materials Properties of Natural Composites Using Molecular Simulation. J. Mater. Chem. 2009, 19, 7251−7262. (37) Hansen, H. S.; Hü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, 998−1032. (38) Plazinski, W.; Lonardi, A.; Hünenberger, P. H. Revision of the GROMOS 56a6CARBO Force Field: Improving the Description of RingConformational Equilibria in Hexopyranose-Based Carbohydrates Chains. J. Comput. Chem. 2016, 37, 354−365.

The corresponding author acknowledges partial support from the Polish National Science Centre (Project No. 2015/17/B/ NZ9/03589 OPUS 9).

■

REFERENCES

(1) Ovodov, Y. S. Structural Chemistry of Plant Glycuronoglycans. Pure Appl. Chem. 1975, 42, 351−369. (2) Ramesh, H. P. F.; Tharanathan, R. N. Carbohydrates-the Renewable Raw Materials of High Biotechnological Value. Crit. Rev. Biotechnol. 2003, 23, 149−173. (3) Smidsrod, O.; Draget, K. I. Chemistry and Physical Properties of Alginates. Carbohydr. Eur. 1996, 14, 6−13. (4) Ridley, B. L.; O’Neill, M. A.; Mohnen, D. Pectins: Structure, Biosynthesis, and Oligogalacturonide-Related Signaling. Phytochemistry 2001, 57, 929−967. (5) Rees, D. A.; Welsh, E. J. Secondary and Tertiary Structure of Polysaccharides in Solutions and Gels. Angew. Chem., Int. Ed. Engl. 1977, 16, 214−224. (6) Morris, E. R. Molecular Interactions in Polysaccharide Gelation. Br. Polym. J. 1986, 18, 14−21. (7) Rees, D. A. Polysaccharide Shapes and Their Interactions - Some Recent Advances. Pure Appl. Chem. 1981, 53, 1−14. (8) Willats, W. G.; McCartney, L.; Mackie, W.; Knox, J. P. Pectin: Cell Biology and Prospects for Functional Analysis. Plant Mol. Biol. 2001, 47, 9−27. (9) Pérez, S.; Mazeau, K.; Hervé du Penhoat, C. The ThreeDimensional Structures of the Pectic Polysaccharides. Plant Physiol. Biochem. 2000, 38, 37−55. (10) Rees, D. A. Polysaccharide Conformation in Solutions and Gels - Recent Results on Pectins. Carbohydr. Polym. 1982, 2, 254−263. (11) Willats, W. G.; Orfila, C.; Limberg, G.; Buchholt, H. C.; van Alebeek, G. J.; Voragen, A. G.; Marcus, S. E.; Christensen, T. M.; Mikkelsen, J. D.; Murray, B. S.; et al. Modulation of the Degree and Pattern of Methyl-Esterification of Pectic Homogalacturonan in Plant Cell Walls. Implications for Pectin Methyl Esterase Action, Matrix Properties, and Cell Adhesion. J. Biol. Chem. 2001, 276, 19404−19413. (12) Draget, K. I.; Smidsrød, O.; Skjåk-Bræk, G. Alginates from Algae. In Polysaccharides and Polyamides in the Food Industry. Properties, Production, and Patents.; WILEY-VCH Verlag GmbH & Co. KGaA: Weinheim, Germany, 2005. (13) Heyraud, A.; Dantas, L.; Courtois, J.; Courtois, B.; Helbert, W.; Chanzy, H. Crystallographic Data on Bacterial (1 → 4)-β-DGlucuronan. Carbohydr. Res. 1994, 258, 275−279. (14) Braccini, I.; Heyraud, A.; Pérez, S. Three-Dimensional Features of the Bacterial Polysaccharide (1 → 4)-β-D-Glucuronan: A Molecular Modeling Study. Biopolymers 1998, 45, 165−175. (15) Plazinski, W. Molecular Basis of Calcium Binding by Polyguluronate Chains. Revising the Egg-Box Model. J. Comput. Chem. 2011, 32, 2988−2995. (16) Plazinski, W.; Rudzinski, W. Molecular Modeling of Ca2+Oligo(α-L-Guluronate) Complexes: Toward the Understanding of the Junction Zone Structure in Calcium Alginate Gels. Struct. Chem. 2012, 23, 1409−1415. (17) Plazinski, W.; Drach, M. The Dynamics of the Calcium-Induced Chain-Chain Association in the Polyuronate Systems. J. Comput. Chem. 2012, 33, 1709−1715. (18) Plazinski, W.; Drach, M. Calcium-α-L-Guluronate Complexes: Ca2+ Binding Modes from DFT-MD Simulations. J. Phys. Chem. B 2013, 117, 12105−12112. (19) Plazinski, W. Conformational Properties of Acidic Oligo- and Disaccharides and Their Ability to Bind Calcium: A Molecular Modeling Study. Carbohydr. Res. 2012, 357, 111−117. (20) Stewart, M. B.; Gray, S. R.; Vasiljevic, T.; Orbell, J. D. Exploring the Molecular Basis for the Metal-Mediated Assembly of Alginate Gels. Carbohydr. Polym. 2014, 102, 246−253. (21) Wolnik, A.; Albertin, L.; Charlier, L.; Mazeau, K. Probing the Helical Forms of Ca2+-Guluronan Junction Zones in Alginate Gels by Molecular Dynamics 1: Duplexes. Biopolymers 2013, 99, 562−571. 3708


Article

The Journal of Physical Chemistry B (39) 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 Force-Field Parameter Sets 53A5 and 53A6. J. Comput. Chem. 2004, 25, 1656−1676. (40) Oostenbrink, C.; Soares, T. A.; van der Vegt, N. F. A.; van Gunsteren, W. F. Validation of the 53A6 GROMOS Force Field. Eur. Biophys. J. 2005, 34, 273−284. (41) Lins, R. D.; Hünenberger, P. H. A New GROMOS Force Field for Hexopyranose-Based Carbohydrates. J. Comput. Chem. 2005, 26, 1400−1412. (42) Reif, M. M.; Winger, M.; Oostenbrink, C. Testing of the GROMOS Force-Field Parameter Set 54A8: Structural Properties of Electrolyte Solutions, Lipid Bilayers, and Proteins. J. Chem. Theory Comput. 2013, 9, 1247−1264. (43) 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: Condens. Matter Mater. Phys. 1988, 37, 785−789. (44) 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, 1200−1211. (45) 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, 11623−11627. (46) Becke, A. D. Density-functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (47) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. AB INITIO Molecular Orbital Theory; John Willey & Sons, Inc.: New York, 1986. (48) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. A WellBehaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97, 10269−10280. (49) Horta, B. A. C.; Fuchs, P. F. J.; van Gunsteren, W. F.; Hü nenberger, P. H. New Interaction Parameters for Oxygen Compounds in the GROMOS Force Field: Improved Pure-Liquid and Solvation Properties for Alcohols, Ethers, Aldehydes, Ketones, Carboxylic Acids, and Esters. J. Chem. Theory Comput. 2011, 7, 1016− 1031. (50) 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., Li, X., et al.. Gaussian 09, Revision A.02; Gaussian, Inc.: Wallingford, CT, 2016. (51) Hanwell, M. D.; Curtis, D. E.; Lonie, D. C.; Vandermeersch, T.; Zurek, E.; Hutchison, G. R. Avogadro: An Advanced Semantic Chemical Editor, Visualization, and Analysis Platform. J. Cheminf. 2012, 4, 17. (52) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. UFF, a Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J. Am. Chem. Soc. 1992, 114, 10024−10035. (53) Gaweda, K.; Plazinski, W. Pyranose Ring Conformations in Mono- and Oligosaccharides: A Combined MD and DFT Approach. Phys. Chem. Chem. Phys. 2017, 19, 20760−20772. (54) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; et al. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845−854. (55) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; et al. PLUMED: A Portable Plugin for Free-Energy Calculations with Molecular Dynamics. Comput. Phys. Commun. 2009, 180, 1961−1972. (56) Abraham, M. J.; Murtola, T.; Schulz, R.; Pá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.

(57) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New Feathers for an Old Bird. Comput. Phys. Commun. 2014, 185, 604−613. (58) 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, E., Ed.; Reidel: Dordrecht, The Netherlands, 1981; pp 331−342. (59) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (60) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (61) Hockney, R. W. Potential Calculation and Some Applications. Methods Comput. Phys. 1970, 9, 136−211. (62) 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, 1463−1472. (63) Hess, B. P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 116−122. (64) 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, 269−286. (65) Barker, J. A.; Watts, R. O. Monte Carlo Studies of the Dielectric Properties of Water-like Models. Mol. Phys. 1973, 26, 789−792. (66) 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, 5451−5459. (67) Heinz, T. N.; van Gunsteren, W. F.; Hünenberger, P. H. Comparison of Four Methods to Compute the Dielectric Permittivity of Liquids from Molecular Dynamics Simulations. J. Chem. Phys. 2001, 115, 1125−1136. (68) Cremer, D.; Pople, J. A. General Definition of Ring Puckering Coordinates. J. Am. Chem. Soc. 1975, 97, 1354−1358. (69) Sugita, Y.; Okamoto, Y. Replica-Exchange Molecular Dynamics Method for Protein Folding. Chem. Phys. Lett. 1999, 314, 141−151. (70) 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, 13435−13441. (71) Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100, 020603. (72) 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, 2783−2792. (73) 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, 670−678. (74) Huggins, M. L. Bond Energies and Polarities1. J. Am. Chem. Soc. 1953, 75, 4123−4126. (75) Sattelle, B. M.; Hansen, S. U.; Gardiner, J.; Almond, A. Free Energy Landscapes of Iduronic Acid and Related Monosaccharides. J. Am. Chem. Soc. 2010, 132, 13132−13134. (76) Sattelle, B. M.; Bose-Basu, B.; Tessier, M.; Woods, R. J.; Serianni, A. S.; Almond, A. Dependence of Pyranose Ring Puckering on Anomeric Configuration: Methyl Idopyranosides. J. Phys. Chem. B 2012, 116, 6380−6386. (77) Mayes, H. B.; Broadbelt, L. J.; Beckham, G. T. How Sugars Pucker: Electronic Structure Calculations Map the Kinetic Landscape of Five Biologically Paramount Monosaccharides and Their Implications for Enzymatic Catalysis. J. Am. Chem. Soc. 2014, 136, 1008−1022. (78) Plazinski, W.; Drach, M.; Plazinska, A. Ring Inversion Properties of 1→2, 1→3 and 1→6-Linked Hexopyranoses and Their Correlation with the Conformation of Glycosidic Linkages. Carbohydr. Res. 2016, 423, 43−48. 3709


Article

The Journal of Physical Chemistry B (79) Plazinski, W.; Plazinska, A. Molecular Dynamics Simulations of Hexopyranose Ring Distortion in Different Force Fields. Pure Appl. Chem. 2017, 89, 1283−1294. (80) Rao, V. S. R. Conformation of Carbohydrates; Harwood Academic Publishers: Amsterdam, The Netherlands, 1998. (81) Angyal, S. J. Conformational Analysis in Carbohydrate Chemistry. Aust. J. Chem. 1968, 21, 2737−2746. (82) Angyal, S. J. The Composition and Conformation of Sugars in Solution. Angew. Chem., Int. Ed. Engl. 1969, 8, 157−166. (83) Panczyk, K.; Plazinski, W. Pyranose Ring Puckering in Aldopentoses, Ketohexoses and Deoxyaldohexoses. A Molecular Dynamics Study. Carbohydr. Res. 2018, 455, 62−70. (84) Lemieux, R. U.; Koto, S.; Voisin, D. The Exo-Anomeric Effect. In Anomeric Effect. ACS Symp. Ser. 1979, 87, 17−29. (85) Lonardi, A.; Oborský, P.; Hünenberger, P. H. SolventModulated Influence of Intramolecular Hydrogen-Bonding on the Conformational Properties of the Hydroxymethyl Group in Glucose and Galactose: A Molecular Dynamics Simulation Study. Helv. Chim. Acta 2017, 100, e1600158. (86) Wang, D.; Á mundadóttir, M. L.; van Gunsteren, W. F.; Hü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, 521−537. (87) 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. (88) Tvaroška, I.; Bleha, T. Anomeric and Exo-Anomeric Effects in Carbohydrate Chemistry. Adv. Carbohydr. Chem. Biochem. 1989, 47, 45−123. (89) Miljkovic, M. Electrostatic and Stereoelectronic Effects in Carbohydrate Chemistry; Springer US: Boston, MA, 2014. (90) Project, E.; Nachliel, E.; Gutman, M. Parameterization of Ca+2Protein Interactions for Molecular Dynamics Simulations. J. Comput. Chem. 2008, 29, 1163−1169. (91) Morris, E. R.; Rees, D. A.; Thom, D.; Boyd, J. Chiroptical and Stoichiometric Evidence of a Specific, Primary Dimerisation Process in Alginate Gelation. Carbohydr. Res. 1978, 66, 145−154. (92) Grant, G. T.; Morris, E. R.; Rees, D. A.; Smith, P. J. C.; Thom, D. Biological Interactions between Polysaccharides and Divalent Cations: The Egg-Box Model. FEBS Lett. 1973, 32, 195−198.

3710


CARBO_R Force Field for

Recommend Documents