Mehdi
Mobli
and
Andrew
Almond
*
Manchester Interdisciplinary Biocentre, 131 Princess Street, University of Manchester, Manchester, UK M1 7DN. E-mail: Andrew.Almond@manchester.ac.uk; Fax: +44 161 30 68918; Tel: +44 161 30 64199
First published on 31st May 2007
N-Acetylated amino sugars are essential components of living organisms, but their dynamic conformational properties are poorly understood due to a lack of suitable experimental methodologies. Nuclear magnetic resonance (NMR) is ideally suited to these conformational studies, but accurate equations relating the conformation of key substituents (e.g., the acetamido group) to NMR observables are unavailable. To address this, density functional theory (DFT) methods have been used to calculate vicinal coupling constants in N-acetylated amino sugars and derive empirical Karplus equations for 3J(HNH2) of N-acetyl-D-glucosamine (GlcNAc) and N-acetyl-D-galactosamine (GalNAc). The fitted Karplus parameters were found to be similar to those previously derived for peptide amide groups, but are consistently larger in magnitude. Local intramolecular interactions had a small effect on the calculated J-couplings and comparison with experimental data suggested that DFT slightly overestimated them. An implicit solvation model consistently lowered the magnitude of the calculated values, improving the agreement with the experimental data. However, an explicit solvent model, while having a small effect, worsened the agreement with experimental data. The largest contributor to experimentally-determined 3J(HNH2)-couplings is proposed to be librations of the amide group, which are well approximated by a Gaussian distribution about a mean dihedral angle. Exemplifying the usefulness of our derived Karplus equations, the libration of the amide group could be estimated in amino sugars from experimental data. The dynamical spread of the acetamido group in free α-GlcNAc, β-GlcNAc and α-GalNAc was estimated to be 32°, 42° and 20°, with corresponding mean dihedral angles of 160°, 180° and 146°, respectively.
The GAGs are a family of carbohydrate polymers, named according to their basic repeating disaccharide, which comprises an amino sugar and a single uronic acid sugar. For instance, the GAGs hyaluronan and heparan sulfate contain GlcNAc, while chondroitin sulfate contains N-acetyl-D-galactosamine (GalNAc), see Fig. 1. They are found in all vertebrate tissue and although it is known that GAGs partake in cell signaling, development and growth,2,6,8 many of their functions are not presently understood. One of the major reasons for this poor appreciation of GAG function is the current lack of detailed information about their dynamic 3D-organization in solution.7,9 Among experimental techniques, vibrational and nuclear magnetic resonance (NMR) spectroscopy are applicable to these investigations in aqueous solution. While vibrational spectroscopy can yield some of this information, extensive theoretical analysis is currently required to deconvolute the resultant spectra.10 A more suitable technique is NMR, which can determine 3D-structural information at atomic resolution in relatively complex molecules. Unfortunately, the repeating nature of GAGs results in little chemical shift dispersion, which has reduced the effectiveness of NMR in this area.11,12 However, recent advances in NMR instrumentation and experiments, coupled with chemical and biological methods for derivation, purification and labeling, are finally allowing the 3D solution structures of these complex carbohydrates to be investigated with unprecedented atomic accuracy.11,13–15
![]() | ||
Fig. 1 Structure and ring numbering of the α- (A) and β- (B) anomers of N-acetylated amino sugars; only the anomeric-ring hydroxyl group is shown. The two hydrogen atoms involved in the 3J(HNH2)-coupling are shown together with the corresponding angle θ. In GlcNAc (1) all ring hydroxyl groups are equatorial; GalNAc (2) differs in that the hydroxyl group at the C4-position is axial. |
Structural information is usually obtained from NMR experiments using the nuclear Overhauser effect (NOE) or residual dipolar couplings (RDC) induced by weak alignment. Furthermore, the scalar coupling (also referred to as spin–spin or J-coupling) between nuclei can be used to provide structural information if the relationship between bond geometry and coupling is known. This approach has been used to determine the conformational dependence of several saccharides,16–19 but in many cases, such as around the acetamido moiety of the important amino sugars, the geometry–coupling relationship has not been investigated. In general, this relationship can be described by the Karplus equation, an empirical relationship relating coupling constants to bond dihedral angle and a second order expansion in the cosine of this angle, which is regularly used in chemistry to obtain geometrical information. The parameters needed to define a Karplus relationship are traditionally determined using coupling constants of molecules with known geometries and are then used to analyze systems with unknown geometries.20 However, one may calculate coupling constants using theoretical quantum-chemistry methods for any arbitrary geometry. This approach has been unpopular because significant computational effort is required to produce accurate results.
With the advance of computational methods and power it is now feasible to pursue the theoretical approach more routinely. Recently, it has been shown that density function theory (DFT) can calculate scalar couplings very accurately by employing basis-sets that describe the electronic environment near the nucleus well.21,22 It has also been shown that the computational cost can be greatly reduced by calculating the dominating Fermi contact (FC) term of the scalar coupling using a larger basis-set and the computationally more demanding but smaller remaining contributions from the paramagnetic spin–orbit (PSO), diamagnetic spin–orbit (DSO) and spin dipolar (SD) terms using a smaller basis-set.21–23 Modern DFT calculations hold the promise of unprecedented accuracy in geometry–coupling estimation because all dynamics are frozen-out in calculations. This is not possible directly from experimental measurements of even the most rigid molecules, which generally report lower values than theoretical methods due to local librations.24 Furthermore, DFT calculations that include solvent effects show improved accuracy, giving results in close agreement with the experimental NMR measurements.25,26
Hydrogen-bonding interactions (intramolecular and intermolecular) are essential for defining carbohydrate 3D-structure in aqueous solution.27 The theorized role of the amide group as both a hydrogen-bond donor and acceptor in solvated GAGs makes detailed analysis of this group particularly important.28 In N-acetylated 2-deoxyamino sugars, the H2 proton is fixed within a pyranose ring and its coupling to the HN proton provides detailed information of the rotation of the entire acetamido group. Thus, we have employed the recent developments in DFT calculations to derive Karplus equations suitable for use in the amino sugars found in GAGs (GlcNAc and GalNAc). The change in the derived Karplus parameters with a change in the anomeric form of the amino sugar was investigated in GlcNAc, where both the α- and β-anomers were considered. Also, the effects of including an implicit or an explicit solvent model was explored, the latter with the aid of data from molecular dynamics (MD) simulations. Since the interaction with water has been found to be important for determining carbohydrate conformation,29 we employed classical MD simulations with explicit water molecules, rather than more complex ab initio MD methods which are generally used in vacuo.30 Information from the MD simulations was also used to model dynamical effects and an equation was derived to estimate the dynamic conformational spread of the acetamido moiety directly from a derived Karplus equation and NMR measurements.
All coupling constants were calculated using structures optimized at the B3LYP/6-31G(d,p) level of theory. The calculations were performed using the Gaussian 03 software23 rev. D (where the FC calculation can easily be calculated separately from the spin–spin coupling terms). The implicit PCM method implemented in the G03 software was used where specified.
J = Acos2(θ +φ) + Bcos(θ +φ) + C | (1) |
In this case, the angle θ was defined to be the H2–C2–NH–HN torsion (Fig. 1), which by symmetry of the electronic orbitals allows φ to be zero for the three-bond coupling. This torsion was rotated through 360° in 30° degree increments, and at each point minimization was performed, followed by spin–spin coupling constant calculation. A Karplus-type equation was derived by non-linear least-squares fitting the resultant spin–spin couplings to the general eqn (1). In this procedure, implicit solvation (when used) was applied during both minimization and spin–spin coupling calculation.
Although eqn (1) only allows calculation of spin–spin coupling for a static structure, it can be generalized to the case of harmonic motion about a mean angle θ by integration. If the motion is assumed to be Gaussian with a standard deviation of σ, then the A, B and C parameters in eqn (1) can be replaced by the modifications shown by eqn (2), as described previously.32
![]() | (2) |
Assuming that the libration (σ) is small, it is possible to simplify eqn (2) considerably assuming B′
≅
B, since the exponential term in B′ is very close to 1 for small angles. Furthermore, by substitution of these approximated equations into eqn (1), it is possible to derive a relationship between the dynamic spread, σ, the measured coupling, J, and the original Karplus coefficients from eqn (1), which is shown in eqn (3). Such an equation can be used to calculate the angular spread for librations around the cis or trans conformation (by suitable substitution of ), from an experimentally measured coupling and derived Karplus parameters. Due to the assumption made for B′, the result of eqn (3) will always be slightly larger than the correct value, but will be negligible for small angles.
![]() | (3) |
Models of the sugars that included several important explicit water molecules were constructed from a representative snapshot of the MD simulation (selected using statistical analysis of water interactions during the complete MD simulation) by deletion of non-essential water molecules and minimization using a molecular mechanics force field (MMX34 in PCModel35), which includes a hydrogen-bonding potential. The resulting structures were then minimized without any restraints using the aforementioned DFT method. This procedure was iterated until a structure was found that both satisfied the MD data and converged to a low-energy minimum according to DFT.
A | Error | B | Error | C | Error | |
---|---|---|---|---|---|---|
a Values derived using Ace–Ala–NMe. b Values derived using Ala–Ala–NH2. | ||||||
[1] 1α-Full | 9.81 | ±0.23 | −1.51 | ±0.12 | 0.62 | ±0.14 |
[2] 1α-FC | 9.56 | ±0.23 | −1.62 | ±0.11 | 0.69 | ±0.14 |
[3] 1α-FC[PCM] | 9.60 | ±0.20 | −1.51 | ±0.10 | 0.99 | ±0.12 |
[4] 1β-FC | 9.45 | ±0.26 | −2.08 | ±0.13 | 0.63 | ±0.16 |
[5] 2α-FC | 10.02 | ±0.21 | −1.79 | ±0.10 | 0.49 | ±0.13 |
Ref. 24 a | 9.44 | — | −1.53 | — | 0.07 | — |
Ref. 24 b | 9.14 | — | −2.28 | — | −0.29 | — |
![]() | ||
Fig. 2 Calculated non-FC (SD = ‘—’, PSO = ‘![]() |
![]() | ||
Fig. 3 Calculated 3J(HNH2) plotted as a function of the dihedral angle θ (note arbitrary line-fit). A) DFT calculated scalar coupling constants for α-GlcNAc. Full calculation (SD + PSO + DSO + FC) = ‘—’, FC-only = ‘![]() ![]() |
Each geometry minimization was first performed using an empirical force-field followed by DFT at fixed angle θ. In the latter high-level minimizations, only local minima could be found and hence the rotamer states of hydroxyl moieties did not change. For some values of θ, the demand that this angle should remain fixed led to deformation of the amide plane during DFT minimization (see Fig. 4), which resulted from a steric clash between the carbonyl oxygen and the O3 hydroxyl oxygen. Consequently, in successive minimizations the O3 hydroxyl group could be found either pointing towards or away from the carbonyl group. These two orientations for the O3 hydroxyl group (referred to as T and A, respectively) resulted in differing levels of amide deformation (Fig. 4). However, when the coupling constant was calculated for both the A and T orientations at θ = 180°, the calculated difference in coupling was less than 0.2 Hz (Table 2). Furthermore, although the extent of deformation was greater when θ was fixed at around 90° or 270°, the calculated vicinal coupling constants, and subsequently any errors due to deformation, were found to be very small. To investigate this further, the FC term was calculated for α-GlcNAc with θ fixed at 180°, and the out-of-plane bending of the amide nitrogen was varied between 0° and 30° in 10° increments. The result was a change of the vicinal coupling of 0.28 Hz (i.e., less than 3% of the coupling constant), further suggesting that these amide deformations have little effect on the accuracy of the scalar coupling constants.
Solvent | θ/° | OH config. | Calc. 3J(HNH2)/Hz | Pred. 3J(HNH2)/Hz | |
---|---|---|---|---|---|
α-GlcNAc | |||||
Vacuo | 180 | A | 11.9 | 11.9 | |
Vacuo | 180 | T | 12.1 | 11.9 | |
PCM | 180 | A | 11.7 | 11.9 | |
PCM | 180 | T | 12.0 | 11.9 | |
Explicit | 168 | T | 11.2 | 11.4 | |
β-GlcNAc | |||||
Vacuo | 180 | T | 12.5 | 12.2 | |
Vacuo | 180 | A | 12.5 | 12.2 | |
PCM | 180 | T | 11.9 | 12.2 | |
Explicit | 163 | A | 11.8 | 11.3 |
![]() | ||
Fig. 4 Overlay of 12 DFT minimized structures of GlcNAc (top α, bottom β), with 30° incremental rotation of the HN–H2 dihedral angle θ. The carboxyl oxygen is represented by a sphere and the NH–H2 bond by a cylinder. |
One may use the potential energy surface calculated during the geometry optimization to estimate an average J-coupling. The 12 static conformers generated were Boltzmann weighted according to eqn (4), where conformer a has internal energy Ea and a probability of observation of pa (kB is the Boltzmann constant and T is the temperature).
![]() | (4) |
The average J-coupling is then estimated from the relative populations of the 12 conformers and their calculated couplings. This was performed on GlcNAc and the resulting average values were 5.09 Hz and 1.04 Hz for the α- and β-anomers respectively (cf. experimental values of 8.9 Hz and 9.1 Hz, respectively). The very poor agreement of the calculated average couplings of the β-anomer with the experimental data is due to the potential energy surface of this anomer showing a minimum energy value around 270°, where a stabilizing intramolecular hydrogen bond is formed. Similarly for the α-anomer, the minimum energy is found in the cis orientation (near 60°), where a hydrogen bond is formed between the carbonyl group of the amide and the anomeric hydroxyl group. Similar results were found when including an implicit solvent model (data only calculated for the α-anomer). Again the stabilizing intramolecular hydrogen bond in the cis orientation (near 30°) yielded the most stable conformation giving an average J of 7.29 Hz. These results further support the need for explicit water molecules to accurately predict the populations of individual conformers.
While the implicit solvation model is computationally efficient, it ignores important electronic effects that may result from local hydrogen-bonding interactions to solvent water molecules. Therefore, in order to investigate the effect of these specific interactions on coupling constants, explicit models of water interaction were produced from 5 ns molecular dynamics (MD) simulations of GlcNAc using a molecular mechanics force-field. The mean orientation of various rotatable bonds and their standard deviations were calculated from the simulations (Table 3). From the table, it can be seen that the H2–C2–NH–HN torsion is predominantly in the trans orientation in both cases, with mean torsion angles of 161° and 180° for the α- and β-anomers, respectively. Furthermore, the MD simulations were used to identify the positions of key water molecules in GlcNAc (Table 4). Those water molecules found to be involved in hydrogen bonds and water bridges to the amide were suitably positioned in a static model, which was subsequently optimized using DFT. The two optimized structures for α- and β-GlcNAc, with attendant explicit water molecules, are shown in Fig. 5. These geometries were used as a basis for calculating the FC contribution to the three-bond scalar coupling. These results are also given in Table 2, which includes calculations of the two anomers of GlcNAc (for both the A and T orientations) excluding any solvent effects. Furthermore, the scalar couplings have also been calculated by using the Karplus relationship derived in vacuo (see above).
α-GlcNAc | β-GlcNAc | |||||
---|---|---|---|---|---|---|
Rotamer | g+ | g− | t | g+ | g− | t |
a For HN–NH–C2–H2 the cis orientation is given under ‘g−’. | ||||||
O5–C1–O1–H | ||||||
Population (%) | 94 | 2 | 4 | 29 | 64 | 7 |
![]() |
45 | 341 | 157 | 36 | 315 | 201 |
σ | 20 | 22 | 20 | 21 | 23 | 20 |
C2–C3–O3–H | ||||||
Population (%) | 13 | 57 | 29 | 9 | 52 | 39 |
![]() |
56 | 287 | 196 | 63 | 283 | 195 |
σ | 21 | 22 | 23 | 20 | 21 | 23 |
C3–C4–O4–H | ||||||
Population (%) | 78 | 13 | 9 | 76 | 12 | 12 |
![]() |
65 | 310 | 144 | 66 | 307 | 146 |
σ | 23 | 21 | 17 | 23 | 21 | 18 |
C4–C5–C6–O6 | ||||||
Population (%) | 57 | 2 | 41 | 53 | 2 | 45 |
![]() |
61 | 282 | 187 | 60 | 284 | 186 |
σ | 10 | 17 | 13 | 11 | 17 | 13 |
C5–C6–O6–H | ||||||
Population (%) | 32 | 21 | 47 | 32 | 23 | 45 |
![]() |
66 | 292 | 180 | 67 | 288 | 180 |
σ | 23 | 25 | 25 | 22 | 24 | 25 |
H–N–C2–H | ||||||
Population (%) | — | 0 | 100 | — | 13 | 87 |
![]() |
— | — | 161 | — | 3 | 180 |
σ | — | — | 29 | — | 57 | 21 |
![]() | ||
Fig. 5 Stereo view of the ab initio optimized structure of GlcNAc, including water bridges identified from MD simulations (A = α, B = β). The dashed lines show intermolecular hydrogen bonds. |
The change in the scalar coupling on inclusion of explicit water molecules does not follow a definite trend, which is expected since the solvent molecules are directly interacting with the solute and thus the change in electronic configuration is complex. It is noted that although only the α-anomer model has direct interactions between the amide proton and water molecules (Fig. 5), it is the β-anomer that shows the greatest deviation from the derived Karplus relationship defined above. Therefore, from the data in Table 2 it is concluded that the deviation from the 3J(HNH2) due to solvent effects is less than 0.6 Hz (i.e., less than 5%).
Using the H2–C2–NH–HN dihedral angle data from the MD simulations and the Karplus parameters (defined earlier), an average value for the scalar coupling constant could be determined for α/β-GlcNAc. The coupling constants were calculated by integrating the Karplus equation, as described in eqns (1) and (2), about an angle with a Gaussian distribution (standard deviation σ) The average angles and standard deviations extracted from the MD simulations were used (Table 2). This latter method was repeated for α-GalNAc, using the same
and σ parameters as for α-GlcNAc. The results are presented in Table 5 together with the experimentally-measured coupling constants.
α-GlcNAc/Hz | β-GlcNAc/Hz | α-GalNAc/Hz | |
---|---|---|---|
a Using corresponding parameters from Table 1.
b Parameters from [1] in Table 1.
c Parameters from [2] in Table 1.
d Parameters from [3] in Table 1.
e Parameters from [4] in Table 1.
f Modified parameters from [1] in Table 1; A′ = 5.9, B′ = −1.3, C′ = 0.6; ![]() ![]() ![]() ![]() ![]() |
|||
Experimental | 8.88 | 9.07 | 8.42 |
MD averagea (θ both cis and trans) | 8.93b | 10.39e | — |
8.91c | — | — | |
9.13d | — | — | |
Gaussian average (θ exclusively trans) | 9.10f | 10.91i | 9.35j |
9.07g | — | — | |
9.29h | — | — |
Based on these results it may be concluded that the calculated values for the H1–H2 scalar coupling in α-GlcNAc are slightly overestimated in the DFT calculations. One may speculate that it is the solvent effects that are poorly modeled. However, experimental results41 on a form of GlcNAc where all the hydroxyl groups are acetylated shows that there is no change in this coupling constant (the data are given to 0.1 Hz) when going from H2O to CDCl3, suggesting that solvent effects are likely to be small. The remaining sources for this overestimation are either due to the basis-set or the DFT method itself. A larger basis-set is not expected to improve the current calculations significantly, based on a previous study.22 Furthermore, the same study found that overestimation of coupling constants beyond motional dynamics is an intrinsic property of the DFT methodology.
It is also noted that there are no simple factors that can be used to correct for errors resulting from DFT calculations, suggesting that the extent of this overestimation is coupling-specific. For the trans three-bond H2–H3 coupling the error is likely to be less than 10%, as any greater decrease would lead to calculated values lower than those measured experimentally. It is therefore concluded that although this inherent systematic error is small, one must be vigilant when using Karplus parameters fitted to such calculations as they are likely to result in either overestimation of the bond librations or inaccurate mean dihedral angles.
It should be noted that the MD simulation of β-GlcNAc shows a bimodal distribution of acetamido angles (while the trans orientation of the θ dihedral is dominant, cis is also populated), which cannot be taken into account by eqn (3). Averaging only the trans conformations from the MD simulation, via fitting to eqn (3), results in a larger predicted coupling and a worse agreement with experimental data than averaging the whole simulation (Table 5). Therefore, if this dihedral angle populates both conformers in reality, the assumption that θ is exclusively trans for free GlcNAc will result in an overestimation of the spread (σ) of dihedral angles at the acetamido group. Alternatively, the MD simulation both overestimates the relative population of the cis conformer and underestimates the dynamic spread around the trans conformation.
In Table 5 the experimental and calculated coupling constants (using the various approaches) are shown. The calculated coupling constant using various Karplus equations derived for α-GlcNAc differ little. The coupling constant calculated including an implicit solvent model (PCM) is surprisingly higher than when a solvent model is excluded. This difference appears as the geometry optimization including the implicit solvent model favors the A orientation of the O3 atom, and in this orientation the coupling constant is slightly higher than for the T orientation favored in the in vacuo optimizations. It is interesting to note that local geometrical changes have a similar (but slightly larger) effect on the coupling constants than the calculated solvent effect. Similar results are obtained when using the modified Karplus parameters compared to averaging over the MD population, suggesting that the distribution of conformers is Gaussian in the MD simulation.
In α-GalNAc the experimental coupling constant is 8.42 Hz, which is in disagreement with the calculated value reported in Table 5 (9.35 Hz). This calculation is made using the Karplus parameters derived specifically for α-GalNAc, but the θ and σ values of eqn (1) and (2)are those determined for α-GlcNAc. Thus it would appear that the mean dihedral angle or the extent of libration is different in these two molecules. The crystal structure of α-GalNAc is known and the dihedral angle of interest is reported to be 146°.43 Inserting this value into eqn (3), a dynamical spread of 20° is found for this molecule. This suggests that both the dynamics and the structure of α-GalNAc are different from those of α-GlcNAc. This procedure can easily be repeated for more complex systems such as in hyaluronan where the 3J(HNH2) coupling constants for each amino sugar have been measured in the hexasaccharide.13 We find by the methods presented here that the amide group libration of the central ring (γ) in this molecule is 34°, suggesting that this group is less dynamical than the corresponding one in the free amino sugar.
In general, inclusion of the explicit solvent molecules worsens the agreement between the experimental and calculated coupling constants (see Tables 4 and 5). This may be due to the fact that the intermolecular interactions with water are very dynamic and are thus difficult to simulate using static DFT calculations. However, the models that included explicit solvent molecules are much more conformationally realistic than those without the solvent model. This is evident when the α-GlcNAc structure is extracted and superimposed with the available crystal structure (Fig. 6). The positioning of the ring hydroxyl groups that form hydrogen bonds with water molecules is nearly identical to those determined by X-ray crystallography. The only differences between the two structures are a deviation of the amide plane from planarity by 20° in the crystal structure (cf. planar in the ab initio calculation) and a deviation of the HN–H2 dihedral angle by a compensating amount in the ab initio calculated structure (cf.trans in the crystal structure). These discrepancies aside, this method does appear to be able to generate realistic 3D-structures for these highly-solvated carbohydrates in silico. Furthermore, one may argue from these data that the α-anomer is involved in more hydrogen bonds, see Fig. 6, which may explain the higher stability of this anomer in aqueous solution than predicted by in vacuo calculations (where the more sterically-relaxed β-anomer has the lower energy).
![]() | ||
Fig. 6 Overlay of the crystal structure and ab initio refined average MD-structure of α-GlcNAc. |
Calculated implicit solvent effects are small and lower the DFT-calculated coupling constants, which improves the agreement with the experimental data. The effect of local geometrical parameters, excluding the dependence on the intervening dihedral angle, is of the same order as implicit solvent effects. Due to the modest contributions from these effects it is concluded that the 3J(HNH2) and thus the Karplus parameters derived here are not sensitive to these effects. Inclusion of explicit solvent molecules did not improve the coupling constant calculations, but it was shown that these are important to maintain a realistic structure for the sugars during DFT minimization. For non-rigid bonds the largest effect on the scalar coupling constant is due to dynamical effects (bond libration). By integration of the Karplus equations, a relation is presented that can estimate the range of libration of the acetamido group from the DFT-derived Karplus parameters presented here and an experimental 3J(HNH2). Application of this method to the sugars studied here suggests that the dynamical spreads at the acetamido groups of α-GlcNAc, β-GlcNAc and α-GalNAc are 32°, 42° and 20°, respectively.
Footnote |
† Electronic supplementary information (ESI) available: Calculated energies, relative populations and calculated coupling constants of GlcNAc. See DOI: 10.1039/b705761j |
This journal is © The Royal Society of Chemistry 2007 |