Computational prediction of the spectroscopic chemical shifts of nuclear magnetic resonance (NMR) can significantly improve its capability as an essential technique for the identification and characterization of complex biomolecular structures in solution (Duus et al., 2000; Wüthrich, 2003; Bifulco et al., 2007; Lodewyk et al., 2012; Tantillo, 2013; Iron, 2017; Krivdin, 2019). NMR techniques play a vital role in the structural elucidation of one of the least structurally characterized classes of biological molecules—carbohydrates—since they cannot otherwise be characterized via methods such as X-ray crystallography (Van Halbeek, 1994; Duus et al., 2000; Toukach and Ananikov, 2013a). Carbohydrates have a diverse set of building block residues that can be tied together with a range of inter-residue linkages to form complex macromolecules of biological significance in plants and animals (Duus et al., 1999). Plants harness the sugar building blocks available in nature to build large, complex, and heterogenous polysaccharide structures that constitute the plant cell wall.

Rhamnogalacturonan-II (RG-II) is a pectin molecule that exemplifies this phenomenon as the most complex polysaccharide known in nature, consisting of 12 monosaccharide units interconnected by 21 glycosidic linkages (Ndeh et al., 2017; O’Neill et al., 2020). Apiose is one of the 12 monosaccharide units in RG-II and has special significance as it is involved in dimerizing RG-II via a borate–ester cross-link (O'Neill et al., 1996; Bharadwaj et al., 2020). The formation of this borate–ester cross-link is required for normal plant growth, and most RG-II found in the plant cell wall remains cross-linked (O'Neill et al., 2004). The chemical structure and topology of RG-II (Figure 1) have only been recently determined by systematically deconstructing RG-II using a multitude of linkage-specific hydrolase enzymes (Ndeh et al., 2017). However, the three-dimensional (3D) structure of RG-II in solution in its monomeric and dimerized states still remains to be elucidated.

Characterizing the 3D structure of complex polysaccharides using NMR is challenging due to the myriad linkages between monosaccharides that often only differ from each other in stereochemistry (Duus et al., 2000). Furthermore, the diverse palette of available monosaccharides allows for structurally complex molecules, all of which share very similar chemical groups (hydroxyl groups and glycosidic linkages) resulting in NMR chemical shifts located in a much narrower spectral region than for proteins and nucleic acids (Toukach and Ananikov, 2013a), leading to overlaps in spectral data that limit the differentiation between specific residues within a polysaccharide (Bubb, 2003).

Computational approaches play an important role in overcoming these challenges. The use of empirical methods that are based on experimental chemical shift databases for spectral assignments may offer quick predictions, but are frequently inaccurate for carbohydrates (Pierens, 2014). The wide diversity of substructures and the complex 3D arrangements of polysaccharides far exceed the parameterization space employed for the development of these databases (Pierens, 2014). As a result, these methods are often unable to distinguish between diastereomers or account for interspatial coupling (Willoughby et al., 2014; Benassi, 2017), greatly hindering the elucidation of stereochemistry within a polysaccharide. On the other hand, quantum mechanical (QM) approaches to predicting NMR shifts have been more successful in aiding experimentalists to improve the reliability of NMR and correlate the spectra to the structure by predicting the chemical shift values of complex compounds (Pierens, 2014; Tarazona et al., 2018; Ahmed et al., 2020). This approach has even been used to correct previously misassigned spectra that relied on experimental databases and empirical-based predictions (Pierens, 2014).

QM-based approaches use the gauge-invariant atomic orbital (GIAO) method to calculate the isotropic magnetic shielding values from electronic structure calculations (Ditchfield, 1974; Schreckenbach and Ziegler, 1995). Early studies applied the GIAO approach to predict NMR shifts for small molecules and relied on Hartree–Fock (HF)- and post-Hartree–Fock (post-HF)-based methods (Schreckenbach and Ziegler, 1995; Helgaker et al., 1999). However, the high computational costs of the HF and post-HF-based ab initio methods for larger molecules, such as carbohydrates, have prompted investigations into using density functional theory (DFT) as a more tractable method (Taubert et al., 2005; Bagno et al., 2007; Toukach and Ananikov, 2013a). However, protocols that utilize DFT have been shown to have limitations (Teale et al., 2013; Benassi, 2017), with many previous studies having large discrepancies from experimental data, including being unable to correctly predict the ordering of chemical shifts (Taubert et al., 2005; Bagno et al., 2007), which is a vital piece of information for spectral analysis. Additionally, several studies have struggled to replicate either ¹H or ¹³C data, often with proton data showing little variation despite changes in the levels of theory while carbon data were shown to vary greatly depending on the chosen level of theory (Benassi, 2017).

It has been shown that changing the specified density functional and basis set (FBS) combination can greatly affect the level of accuracy of a computational prediction (Pierens, 2014). Another important factor is accounting for the effect of solvent, which can be modeled using either implicit (Cramer and Truhlar, 1999) or explicit (Taubert et al., 2005; Bagno et al., 2007) solvation models in QM calculations, molecular dynamics (MD), and molecular mechanics (MM) simulations. Implicit solvation models approximate the presence of a solvent using a dielectric constant, thereby imparting a screening effect. While explicit solvation models may be considered to be intuitively more accurate, Bagno et al. demonstrated that the presence of explicit water molecules did not result in improved ¹³C and ¹H shift predictions for glucose (Bagno et al., 2007). This may be attributed to the fact that, in explicit solvation models, the shielding constant calculations are affected by the conformational distribution of the hydroxyl and hydroxymethylene groups in glucose, which in turn are impacted by the presence of explicit solvent molecules, the number of solvation shells considered, and the inherent effects of the force fields considered for water and glucose. Hence, considering the additional computational costs and the lower accuracy of explicit solvation models for NMR shift calculations, implicit solvation models were considered in this study.

Here, we applied the DFT approach and evaluated the accuracy of various FBS combinations and solvation models in predicting the chemical shifts of ¹³C NMR for four apiosyl compounds (Figure 2). The ability to accurately predict how the chemical shifts of apiose change upon dimerization with boric acid will aid experimentalists in using NMR to probe the dimerization process of RG-II. Therefore, this study is focused on apiose and its borate esters and is not a generalized study on the efficacy of these methods for all carbohydrates. Another reason for this focus is the fact that the experimental NMR shift data from Ishii and Ono (Ishii and Yanagisawa, 1998; Ishii and Ono, 1999) for apiose and its borate esters in D₂O using methanol as the ¹³C NMR reference compound are the only data available for benchmarking the computational predictions presented here.

DFT calculations for chemical shift predictions involve two steps, a geometry optimization followed by an NMR shielding constant calculation, each of which may be performed using a specific FBS combination and a solvation model. Recent literature on NMR shift predictions in carbohydrates often has not accounted for solvation during geometry optimizations and has occasionally done so only during NMR shielding constant calculations (Bagno et al., 2007; Toukach and Ananikov, 2013b; Benassi, 2017). We evaluated 23 FBS combinations based on the levels of theory recommended in the literature for NMR chemical shift predictions and applied them using implicit solvation models and a reference correction to study apio–borate compounds. We demonstrate that considering implicit solvation during both the geometry optimization and the shielding constant calculation results in lower mean absolute errors (MAEs) for computational predictions with respect to the experimental data for apio–borate compounds. Additionally, we found that while considering solvation is important during both geometry optimizations and shielding constant calculations, there was a greater positive impact from using implicit solvation during geometry optimization for both charged and charge-neutral compounds. Across all 23 FBS combinations tested, we list the top-performing FBS combinations that had lower MAEs for ¹³C data than those reported in the literature for NMR predictions of carbohydrates (Taubert et al., 2005; Toukach and Ananikov, 2013b; Benassi, 2017), while also replicating the correct ordering of chemical shifts.

RG-II has numerous exotic sugar moieties that have not been extensively characterized through experiments or computational models, and the protocols recommended in the literature for NMR chemical shift predictions for carbohydrates remain to be evaluated for these sugars. Here, we have shown that common techniques, such as gas-phase geometry optimizations and the use of scaling factors to convert chemical shifts, lacked accurate predictive capabilities for the NMR chemical shifts of the apiosyl compounds of RG-II. Modeling of this system is complicated by the diastereomers, charged states, dimerized configurations, solvent interactions, and the effects of boron.

Figure 6 illustrates the impact of the consideration of solvation for geometry optimizations and NMR shielding constant calculations (with each quadrant representing a different solvation combination) for the four compounds of interest (indicated by the four different colors). For each NMR FBS combination, the vertical and horizontal arrows are quantitative representations of MAE improvements upon solvation for the geometry optimizations and shielding constant calculations, respectively. The arrows point to the protocol that enables the highest improvements in MAE for computational predictions. The vertical arrows, which are much larger, consistently point to the lower quadrants, showing the large and consistent benefit of using implicit solvation during geometry optimizations. Modeling implicit solvation benefits all compounds, particularly charged compounds 2–4. The shorter horizontal arrows indicate a lower benefit of applying implicit solvation during shielding constant calculation.

It may be noted that based on the MAE value alone, NMR FBS 23 and 16 appeared to perform best in the lower left quadrant with gas-phase NMR shielding constant calculations. However, a closer look at the gas-phase NMR predictions (Supplementary Table S8) revealed an incorrect ordering of the shifts, with C2/C3 either being in the wrong order or being too similar to distinguish. When considering chemical shift ordering, the r ² values, and SDs in addition to the MAE values across all four quadrants (Supplementary Table S9), the lower right quadrant (both implicit solvation) was the most accurate solvation combination for all NMR FBS.

While there are examples of NMR predictions of carbohydrates in the literature, these primarily studied uncharged compounds without solvent interactions. Our protocol, when compared with other studies in the literature, more accurately reproduced shift values for both charged and uncharged carbohydrates (Taubert et al., 2005; Ishii and Ono, 1999; Toukach and Ananikov, 2013b; Chełmecka et al., 2007; Boese and Martin, 2004). Taubert et al. utilized B3LYP/TZVP (NMR FBS 17) without solvation to predict the chemical shifts of a few saccharides (MAE = 2.28 ppm, r ² = 0.97) (Taubert et al., 2005), but our top-performing NMR FBS coupled with implicit solvation for our compounds produced better MAE (1.50 ppm) and r ² value (0.99) than those in their study. Additionally, Taubert’s top-performing method (NMR FBS 17 in this study) was not adept at modeling compounds 1–4, as seen in Figure 5 (MAE = 2.89 ppm) (Taubert et al., 2005).

The accuracy of our top-performing NMR FBS 1 paired with implicit solvation can also be analyzed to report a root mean squared error (RMSD) of 1.82 ppm, which is equivalent to the MAE value of 1.51 ppm (Supplementary Table S6). This is an improvement over the use of B3LYP/6-31G(d) (NMR FBS 1), as done by Chelmecka et al., for studying glucopyranoside derivatives as their non-solvated gas-phase geometry optimizations and NMR shielding constant calculations yielded larger RMSEs (3–6 ppm) (Toukach and Ananikov, 2013b; Chełmecka et al., 2007). On the other hand, Kupka et al. reported a better RMSE value of 1.75 ppm during their study of glucose derivatives; they did not consider solvent effects when they used B3PW91/6-31+G* [comparable to B3PW91/6-31+G(d), NMR FBS 19] (Toukach and Ananikov, 2013b; Kupka et al., 1999). Applying NMR FBS 19 to our protocol showed that implicit solvation was required to correctly predict shift order (Supplementary Table S8). Belyakov et al. studied saccharides, but their use of PBE/TZ2p (NMR FBS 21) in gas-phase conditions without any solvation considerations yielded a higher RMSE of 2.43 ppm, and our study of NMR FBS 21 showed that it was not adequate at modeling compounds 1–4 (MAE = 2.04 ppm) (Belyakov and Ananikov, 2011; Toukach and Ananikov, 2013b).

Some of the protocols recommended in the literature used a different type of correction—the scaling factor correction—to convert computed isotropic magnetic shielding values to chemical shifts (Bagno et al., 2007; Pierens, 2014). This uses a computational method on a set of compounds, which are then compared to their experimental values so as to calculate scaling factors to convert computed isotropic magnetic shielding values to chemical shift data using a non-constant correction. Pierens et al. used scaling factors to predict experimental values with MAEs of <1.94 for ¹³C and 0.154 ppm for ¹H for a non-carbohydrate compound using the same geometry optimization FBS of B3LYP/6-31+G(d,p) as in this study, but with non-solvated gas-phase optimizations, combined with NMR FBS 1–10, using a mix of gas-phase and implicit solvation models. We found that Pierens’s scaling factor values were unable to reliably predict the NMR data for compounds 1–4, yielding errors of 2.5–8.8 ppm (Supplementary Table S10). Additionally, our protocol of modeling implicit solvation for geometry optimizations and shielding constant calculations followed by a constant reference correction method produced lower MAEs (30% lower for ¹³C and 15% lower for ¹H) than those calculated using the scaling factor correction reported by Pierens.

We also quantified the effect of implicit solvation on the performance of the scaling factors, showing in Supplementary Table S10 that consideration of implicit solvation during both steps greatly reduced the MAEs and improved the relative shift ordering regardless of the correction method. Additionally, considering implicit solvation improved the reference correction method regardless of the reference compound. While using purely gas-phase conditions, TMS as a reference performed rather poorly compared to using methanol. However, as seen in Supplementary Table S5, the accuracy of using any reference greatly improved when implicit solvation was used during both steps, to the extent that TMS could yield lower MAEs than methanol for certain methods. As discussed in Supplementary Information, we continued to use methanol due to other benefits besides MAE and to match the experimental reference, as recommended in the literature (Kupka et al., 1999; Taubert et al., 2005; Sarotti and Pellegrinet, 2009; Pierens, 2014). This study demonstrated that implicit solvation improved predictions regardless of the correction method or the reference used.

Another type of correction uses an internal reference. Bagno et al. referenced the predicted chemical shifts of all carbon nuclei in glucose to the known isotropic magnetic shielding value of anomeric carbon instead of using a reference compound (Bagno et al., 2007). In addition to the limited applicability of internal references, the anomeric atoms in our study tended to have higher errors (Supplementary Table S6), making this approach unreasonable. Additionally, Bagno et al. only studied one compound, so their low MAE of 1.12 ppm was for one FBS that was optimized for one neutral sugar complex, just as our MAE for NMR FBS 1 for compound 3 was only 1.11 ppm. While Bagno et al. showed promise in using MD to take a sampling of 50–100 geometry conformations and averaging their chemical shift predictions, our NMR FBS 3–4 (similar to Bagno’s use of B3LYP/cc-pVTZ to reoptimize structures following MD simulations) did not perform well when applied to compounds 1–4.

The metrics by which we measured the efficacy of a method can also guide improvements to the protocol. On the one hand, a well-performing protocol (i.e., low MAE value) that exhibited high precision (i.e., r ² close to 1) but had low accuracy (i.e., a large MSE) was likely limited by constant systematic error that still was able to replicate the relative ordering of shifts. On the other hand, a well-performing protocol (i.e., low MAE value) that exhibited high accuracy (i.e., a small MSE) but had low precision (i.e., r ² far from 1) was likely limited by a non-constant error that significantly hindered the applicability of the protocol for chemical shift predictions.

NMR FBS 16 exhibited high precision, as evidenced by the best r ² value and the lowest SD, shown in Table 3 and Supplementary Table S9. However, out of all NMR FBS tested, 16 had the second largest MSE, which indicates very low accuracy. This can be observed in Figure 7 where the red linear best-fit line of NMR FBS 16 was significantly but consistently offset from the experiment’s black line. NMR FBS 16 is an example where the protocol predicting the relative shift values and ordering to high accuracy, but is unable to perfectly replicate the magnitude of the chemical shift.

Figure 7 illustrates that the computed values of the top 6 NMR FBS combinations were typically overpredicted compared to the expected experimental values, as seen by the red line (linear fit of the computed values) being above the black line. NMR FBS 8, on the other hand, had the best (smallest) MSE values while still maintaining one of the best r ² values, suggesting high accuracy that came at the price of having a poor MAE value. Figure 7 shows that, unlike the constant offset displayed by NMR FBS 16, both NMR FBS 7 and 8 had different slopes from those of the experimental values, suggesting that they could benefit from a non-constant correction, such as the scaling correction method employed by Pierens (2014). This is consistent with Supplementary Table S10, where our application of Pierens’s non-constant scaling factors performed best for NMR FBS 7 and 8, with MAEs of just over 2.5 ppm. Conversely, the non-constant correction performed worst for NMR FBS 1, whose isotropic magnetic shielding values already showed a linear relationship with the experimental values and greatly benefited from a constant reference correction.

While some studies have suggested that a change in the basis set does not significantly change the predicted shift values (Benassi, 2017), our findings indicated that this does not apply for methyl apioside and its borate esters. While we did see an example of this behavior when the recommended NMR FBS of B3LYP/6-31G* (Chełmecka et al., 2007) produced indistinguishable chemical shifts from NMR FBS 1 that used B3LYP/6-31G(d), the majority of our results pointed to the importance of specifying a FBS combination as opposed to merely recommending a functional for chemical shift predictions. While NMR FBS 23 [WP04/6-31G(d)] was our third top-performing method with an MAE of 1.79 ppm, changing the basis set to aug-cc-pVDZ (NMR FBS 9) increased the MAE by 0.1 ppm, as seen in Figure 5. Conversely, changing the functional to B3LYP (NMR FBS 1) reduced the MAE by 0.29 ppm. Similar trends were observed when comparing NMR FBS 1–4, 13, 17, and 20, as they all shared the same functional (B3LYP) with different basis sets and thus yield various chemical shift predictions. Additionally, while the literature suggested that very large basis sets are needed to obtain accurate NMR predictions (Taubert et al., 2005), we did not find this to be the case for our system as our top-performing NMR FBS combinations did not have particularly large basis sets. These findings, in agreement with similar observations by Kupka et al. (2010) and Ahmed et al. (2020), highlighted the fact that generalized rules of using large basis sets do not inherently lead to accurate chemical shift predictions and thereby necessitate the evaluation of the efficacy of each combination.

In this study, we evaluated the efficacy of various functional and basis set combinations and implicit solvation protocols for the accurate prediction of ¹³C chemical shifts for methyl apioside and its borate derivatives. In addition to the improved ¹³C shift predictions compared to those in recent literature on carbohydrates, we have been able to replicate the correct ordering of shifts for ¹³C data and distinguish between the apiofuranosides upon dimerization. We achieved this by considering implicit solvation during both the geometry optimization and the NMR shielding constant calculations. This protocol correctly predicted the same 6- to 7-ppm change in ¹³C chemical shift for C2/C3 that experimentalists used to distinguish between compound 1 and its borate esters (2–4) (Ishii and Ono, 1999).

Among the 23 FBS combinations that were assessed based on the SD values, r ² analyses, and correct shift ordering, we observed four combinations (FBS 1, 16, 19, and 23) that all reproduced the correct ordering and were well within 2 ppm of the experimental values. We observed that our protocol was also able to demonstrate relatively accurate ¹H predictions in most cases, with similar MAPE values to those observed for ¹³C. In some instances, the ordering of the ¹H shifts was incorrectly predicted and could be attributed to the proton shifts of apiose being highly similar, a fact that was also experienced by experimentalists who struggled to distinguish them without using two-dimensional NMR (Ishii and Ono, 1999). Specifically, our protocol did not replicate the very precise 0.003-ppm change in the ¹H chemical shift for H7 that experimentalists noted to distinguish between the borate esters (2–4).

From the RG-II context, in the future, in addition to being able to predict the ordering of ¹H shifts to greater accuracy, we would also like to extend our protocol to predict chemical shifts in other sugar moieties such as galactose. We believe that this protocol can successfully be applied to study other complex RG-II mono- and oligosaccharides with NMR spectra that are difficult to assign and will play a key role in aiding our understanding of the role of specific sugars in the RG-II dimerization process.

In general, the use of NMR shifts to characterize non-covalent interactions and structural conformations is a crucial approach toward resolving the conformational ambiguities in complex organic molecules (Adams and Lerner, 1992; Siskos et al., 2015). The potential of computational predictions in elucidating structure–spectra relationships hinges on such detailed studies that are pivotal in enabling unique insights into the function and mechanisms of various chemical and biochemical systems.