Insights into the Electronic and Structural Properties of Cellulose and Amylose: A Comparative Force Field Study
Esmat Mohammadi, Justin A. Lemkul
TL;DR
This study compares the structural and electronic properties of amylose and cellulose in water using different simulation methods.
Contribution
The paper introduces a comparative analysis of polarizable and nonpolarizable force fields for amylose and cellulose.
Findings
CHARMM simulations showed stable hydrogen bonding and rigid structures in amylose.
Drude simulations revealed dynamic polarization and flexible conformations in amylose.
Cellulose showed similar structural behavior with both simulation methods.
Abstract
Amylose and cellulose are important biopolymers with diverse applications in biotechnology and materials science. Understanding their structural, dynamic, and solvation properties at the molecular level is critical for harnessing their potential. This study investigates the electronic and structural properties of single-chain cellulose and single- and double-chain amylose in aqueous solution using molecular dynamics simulations with both nonpolarizable (CHARMM) and polarizable (Drude) force fields. CHARMM simulations show stable hydrogen bonding between amylose and water, higher glucose ring dipole moments, increased rigidity, adoption of chair conformations, and less variation in dihedral angles. In contrast, Drude simulations captured dynamic electronic polarization, enhanced conformational flexibility, and resulted in heterogeneous inter- and intramolecular hydrogen bonds. For…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15| Label | Simulated chain | Number of water molecules | Simulation box size (Å) | Simulation time (ns) |
|---|---|---|---|---|
|
| Single-stranded amylose without any restraints | 10723 | 68 | 1000 |
|
| Single-stranded amylose with restrained dihedrals (Φ = 103° and Ψ = 115°) | 1924 | 39 | 200 |
|
| Single-stranded amylose extracted from double helix system with restrained dihedrals (Φ = 85,92° and Ψ = −145°, –153°) | 8965 | 64 | 200 |
|
| Double-stranded amylose | 8886 | 64 | 1000 |
|
| Single-stranded cellulose | 18153 | 81 | 1000 |
- —National Institute of General Medical Sciences10.13039/100000057
- —Division of Materials Research10.13039/100000078
- —National Institute of Food and Agriculture10.13039/100005825
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsAdvanced Cellulose Research Studies · Electrostatics and Colloid Interactions · Surfactants and Colloidal Systems
Introduction
Carbohydrates comprise structurally diverse molecules with relevance to chemistry, biology, materials science, and related fields. ?,? In biological studies, carbohydrates, along with nucleic acids and proteins, play critical functional roles and influence numerous physiological processes.? Consequently, the study of carbohydrates has garnered considerable attention due to their relevance to fundamental biochemical mechanisms and potential developments in the pharmaceutical field. ?,? Two prototypical polysaccharides, amylose and cellulose, serve as a key example of the need to study carbohydrates due to their different structures and functions. Amylose, with its helical shape and role in plant energy storage, consists of linear glucose units linked by α-1,4 glycosidic linkages. ?−? ? Cellulose, which provides structural support in plant cell walls, consists of linear glucose chains linked by β-1,4 glycosidic linkages that form crystalline fibrils. ?,? These polysaccharides differ only in the stereochemistry of their glycosidic linkages, which confers these different properties. Beyond structural differences, amylose is digestible by human enzymes, while cellulose is not due to its β-linkages. ?,?
The behavior of carbohydrates and polysaccharides in water is of specific significance because the dynamic structure and properties of carbohydrates are strongly dependent on their interactions with water. ?−? ? These interactions play a crucial role in many biological processes, such as enzyme recognition and molecular transport. ?,? Furthermore, for industrial use, knowledge of carbohydrate behavior in water is essential for the design of biocompatible materials, such as hydrogels for drug delivery and wound healing. ?,? Overall, elucidation of the intricate interactions of carbohydrates and polysaccharides with water highlights the role of carbohydrates in biological systems and enables the development of new materials with diverse applications.
Molecular dynamics (MD) simulations have long been an important tool in exploring the behavior and interactions of biomolecules, including carbohydrates, at the atomic level. ?,? In carbohydrate science, MD simulations have been particularly valuable for understanding carbohydrate dynamics in water, probing their conformational flexibility, and analyzing their interactions with other biomolecules, including proteins and nucleic acids.? These insights are invaluable to applications ranging from drug design to biomaterials engineering. ?,?
The accuracy and reliability of MD simulations depend strongly on the quality of force field (FF) parameters, which define the interactions between atoms and molecules within the system being simulated. ?,? Specifically, FFs define the potential energy surface responsible for driving the molecular motion and interactions being simulated in MD simulations.? Nonpolarizable FFs such as CHARMM,? GROMOS, ?,? OPLS-AA,? and GLYCAM06? are commonly employed in molecular simulations of carbohydrates and polysaccharides. However, their reliance on fixed partial charges limits their ability to accurately model the intricate environments found in biological systems. For example, quantum mechanical (QM) studies ?,? have demonstrated that MD simulations utilizing nonpolarizable FFs struggle to adequately balance both inter- and intramolecular phenomena such as solvent effects, counterion interactions, hydrogen bonding, and molecular dipole variations, which involve complex charge redistributions. In response, the adoption of polarizable FFs capable of incorporating polarization serves to bridge the gap between QM investigations and classical MD simulations.?
The Drude polarizable FF, based on the classical Drude oscillator model, introduces additional auxiliary particles with negative charge to model electronic degrees of freedom and account for polarization effects, enabling explicit modeling of induced dipoles. ?−? ? This convention proves particularly useful in simulating complex systems with strong electrostatic interactions. However, the parametrization of the Drude FF is challenging, requiring care to avoid overpolarization and extensive testing because of its increased complexity.? Parametrizing carbohydrates is particularly challenging due to dynamic ring puckering, the interdependence of neighboring hydroxyl groups, and the sensitivity of glycosidic linkages to electronic effects. Significant efforts have gone into developing the current Drude polarizable FF for carbohydrates. For instance, laying a foundation for carbohydrate modeling using polarizable FFs, Drude parameters for hexopyranoses were derived and validated against quantum mechanical data. These parameters were tested against dipole moments, vibrational spectra, and torsional profiles to accurately reproduce conformational energetics and molecular polarizability.? Additional work expanded the Drude polarizable FF to enable modeling of polysaccharides containing other monosaccharides such as pyranose and furanose. Parametrization of these monosaccharides was based on QM potential energy scans, dipole moments, and molecular polarizabilities, resulting in a model that agrees well with both QM and experimental data, and improves upon the additive CHARMM36 carbohydrate FF.? Further development extended the parameters to include both N- and O-linked glycoproteins, where the dihedral terms were optimized using solution NMR J-coupling data. This model was validated against several glycopeptides and glycoproteins, enabling compatibility between carbohydrate and protein components in polarizable MD simulations.? Altogether, these developments, along with other studies on the Drude FF development for carbohydrates, ?,? have significantly advanced the ability to simulate carbohydrate systems with greater accuracy, particularly in capturing the delicate balance of electronic effects that govern their structural and dynamic properties.
Cellulose and amylose, which are key biopolymers with wide applications in material science and biotechnology, must be understood at the molecular level to harness their potential. Computational studies have been instrumental in elucidating their structural, conformational, and thermodynamic properties. For amylose, MD simulations have revealed its flexibility and its ability to adopt a helical conformation when complexed with fatty acids, such as linoleic acid, providing insights into its solubility and interactions in aqueous environments.? Other studies have examined the flexibility of single-chain amylose in different solvents, revealing solvent-dependent helical stability and chain extension, which are critical for understanding its behavior in processing and biological environments. ?,?,? In the case of cellulose, atomistic simulations have characterized both its crystalline and amorphous phases, highlighting differences in chain conformations, torsional flexibility, and ring puckering.? The crystalline phase of cellulose exhibits well-defined hydrogen bonding networks and structural parameters consistent with experimental data, while amorphous regions display greater conformational variability and distinct hydrogen bonding patterns that influence mechanical and solvation properties.? Additionally, MD simulations have explored the twisting behavior of cellulose microfibrils, revealing that this phenomenon is driven by van der Waals interactions and modulated by intrachain hydrogen bonding and solvent effects.? Other investigations have modeled cellulose hydration and swelling, demonstrating that water content affects both chain mobility and the stability of hydrogen bond networks.? For instance, studies have shown that water molecules can disrupt the hydrogen bonding between cellulose chains, leading to increased molecular mobility and swelling.? These additive FF studies have advanced our understanding of amylose and cellulose but may be limited by their fixed-charge treatment of electrostatics, which overlooks polarization effects that may be important for modeling the delicate balance between intramolecular hydrogen bonding and interactions with water. The Drude polarizable FF overcomes this limitation by explicitly modeling induced dipoles, offering a more accurate depiction of solvation, hydrogen bonding, and molecular interactions. Thus, it is expected that such a convention will lead to better predictions of structure and dynamics of these carbohydrates.
Here, we have investigated the behavior of amylose and cellulose in water through MD simulations with nonpolarizable (CHARMM36) and polarizable (Drude) FFs. Solvation has been studied through the analysis of parameters like radial distribution function and hydrogen bonding patterns. Comparing results from both FFs allowed us to highlight the influence of electrostatic interactions on solvation and helical structure. Structural analysis, e.g., radius of gyration and dihedral angle analysis, provided descriptions of amylose and cellulose conformations with atomistic resolution. The use of the Drude polarizable FF allowed for a quantitative description of electric fields around cellulose and amylose, revealing details of their electrostatic interactions with solvating water molecules. This approach offers a combined picture of the dynamic behavior and solvation properties of these biopolymers.
Methods and Computational Details
Systems Studied
In this study, we employed both CHARMM36 and Drude FFs to explore the behavior of single-stranded and double-stranded amylose [α(1→4)-d-glucopyranan] and a single-stranded cellulose [β(1→4)-d-glucopyranan] chain. Both polysaccharides contained 12 units of α-d-glucopyranose or β-d-glucopyranose, respectively. All the systems of amylose and cellulose chains were solvated in a box of water. The choice of 12-unit oligomers is consistent with previous MD studies of amylose and cellulose, which have employed chain lengths ranging from 9 to 30 sugar units to probe intrinsic conformational tendencies in dilute aqueous solution. ?,?,?
To generate the initial geometry of a single amylose chain, we used the CHARMM internal coordinate builder and dihedral values of QM-optimized disaccharides. Additionally, we focused on investigating ideal amylose helices, known to form upon interaction with various molecules such as fatty acids.? To this end, we investigated two distinct amylose structures with restrained dihedral angles (Φ, O5_ i -C1 i -O4 i–1 -C4 i–1 , and Ψ, C1 i -O4 i–1 -C4 i–1 -C5 i–1 _). Our goal was to create amylose chains with channel-like central cavities capable of accommodating other molecules. This approach allowed us to compare the properties of water molecules both confined within the helical structure and in the surrounding environment. It also provided insights into the electronic properties, such as the dipole moment and electric field, induced by the amylose chain in its ideal helical conformation when interacting with other molecules. We constructed these initial coordinates using the CHARMM internal coordinate builder with manually specified dihedral values for Φ and Ψ. For the first restrained amylose system, the dihedral angles were taken from an experimental crystal structure of cycloamylose, with values shown in Table.? For the second restrained amylose system, we extracted a single amylose chain from an experimental crystal structure of a double-helix amylose system.? The Φ and Ψ dihedral angles were restrained at the experimental values throughout the simulation.?
1: Simulations Descriptions for Amylose and Cellulose Systems
Furthermore, we performed simulations on a double-helix crystal structure composed of two amylose chains.? This structural arrangement was derived from an experimental study detailing the three-dimensional configuration of the crystalline segment of A-starch.? Each chain consisted of 12 glucopyranose residues distributed across two left-handed, parallel-stranded double helices. The details of the MD simulations conducted for the amylose and cellulose systems are outlined in Table. Renderings of the initial structures of the amylose and cellulose systems are shown in Figure.
Snapshots of the side and front views of different amylose and cellulose systems studied.
General MD Protocol
The initial setup applied the CHARMM36 FF to generate the coordinates and prepare the systems.? Each amylose or cellulose chain was centered in a cubic box of TIP3P water molecules. ?−? ? Energy minimization was carried out using CHARMM36 FF,? utilizing 500 steps of steepest descent minimization followed by 500 steps of adopted-basis Newton–Raphson minimization. Following minimization, a 1 ns equilibration was conducted under an NPT ensemble using OpenMM 7.7.0.? During this period, position restraints were applied to all non-hydrogen atoms with a force constant of 500 kJ mol^–1^ nm^2^. Temperature was maintained at 298 K using the Langevin thermostat method and pressure set to 1 atm employing the Monte Carlo barostat method. ?,? Periodic boundary conditions were imposed in all dimensions, and the short-range van der Waals forces were smoothly switched to zero over a distance of 10–12 Å. Electrostatic interactions were computed using the particle mesh Ewald (PME) method with a real-space cutoff of 12 Å.? Bonds to hydrogen atoms were constrained using SHAKE algorithm,? allowing for a time step of 2 fs in CHARMM simulations.
After the CHARMM equilibration phase, the systems were transformed to the Drude polarizable model. This process involved adding Drude oscillators and lone pairs into the equilibrated coordinates. Simultaneously, the TIP3P water molecules were replaced with the polarizable SWM4-HLJ model.? The relaxation of Drude oscillators was achieved through 1000 steps of steepest descent minimization, followed by 500 steps of adopted-basis Newton–Raphson energy minimization. The Drude systems were then equilibrated under an NPT ensemble, with temperature and pressure maintained at 298 K and 1 atm, respectively, using the same algorithms noted above for the CHARMM systems. The short-range Lennard-Jones potential was smoothly switched to zero from 10–12 Å, and electrostatic forces were calculated using PME.? The polarizable systems were subjected to the same harmonic position restraints and bond constraints as the additive systems, although the integration time step was set to 1 fs to account for the high-frequency oscillations of Drude-atom bonds. A 1 ns equilibration period was employed for the polarizable systems.
Upon completion of equilibration, the production simulations were performed using OpenMM, maintaining the NPT ensemble with the previously described thermostat and barostat settings. The durations of the production simulations are provided in Table. For each system, three independent replicate simulations were produced by generating different, random velocities at the outset of equilibration.
Clustering Analysis
To investigate the sampled conformational diversity across the simulations, we performed clustering analysis employing the MDANCE toolkit, which is used to analyze MD trajectories using stable, ensemble-based methods for clustering.? For this work, we applied k-means N-Ary Natural Initiation (NANI) algorithm to identify dominant structural states. NANI is a centroid initialization method that enhances standard k-means clustering by emphasizing high-density regions of conformational space and selecting diverse initial cluster centers. Contrary to traditional stochastic approaches, NANI is deterministic and provides reproducible results when used with the same data, enhancing clustering accuracy and convergence efficiency.
To determine the appropriate number of clusters, we used the Davies-Bouldin (DB) index, which is a metric for evaluating clustering algorithms. To choose the appropriate number of clusters, two criteria are supported in MDANCE: (1) the lowest DB value and (2) the point corresponding to the maximum second derivative of the DB curve. In our case, we chose the number of clusters based on the lowest DB value, which helps to capture the most clearly defined conformational groupings. Prior to clustering, each of the trajectories was also aligned to a reference frame to remove global translational and rotational motion. Structural similarity was quantified in terms of root-mean-square deviation (RMSD) of heavy atoms, and the representative conformations were defined as the structures nearest the centroid of each cluster in RMSD space.
Results and Discussion
The principal aims of our work were to characterize amylose and cellulose behavior at the molecular level and to compare the properties produced by the CHARMM and Drude FFs to quantify the inclusion of electronic polarization on dynamics and molecular interactions. Our approach explored single-chained cellulose and different conformations of amylose to systematically analyze the structural, hydration, and electronic characteristics within and surrounding these polysaccharide species.
Structural Characterization of Single-Stranded Amylose and Cellulose
We first assessed the structural change in each of the single-stranded systems, in terms of the RMSD, i.e., the difference from their reference configurations. The results of this analysis are shown in the Supporting Information, Figure S1. This analysis revealed that the amylose chains behaved differently using the CHARMM36 and Drude FFs. The CHARMM simulations displayed relatively stable RMSD values throughout the 1-μs trajectories, but frequent transient spikes indicated some conformational variability in the amylose chains. On the other hand, the Drude simulations showed more dynamic amylose behavior, indicating increased conformational dynamics and enhanced flexibility relative to the additive FF. Moreover, there was considerable variability across the Drude replicate simulations, emphasizing the importance of running multiple independent simulations to ensure robust sampling of the conformational space and to assess the variability in the observed dynamics.
To further characterize the conformational variability suggested by the RMSD analysis, clustering was performed on the trajectories using MDANCE, an efficient clustering package for identifying dominant conformational states based on RMSD metrics.? The CHARMM simulations of amylose resulted in an optimal result of 13 clusters based on the DB index (Supporting Information, Figure S2), with the top four clusters each accounting for approximately 10–12% of the total frames (Supporting Information, Figure S3). The top four clusters shared a common feature: extended helical conformations with a longer helical pitch compared to the initial structure (FigureA and Supporting Information, Figure S4). In contrast, the Drude simulation produced 28 clusters, with the largest cluster containing only ∼8% of the total frames, and 25 of these clusters contributing less than 5% (Supporting Information, Figure S3). This distribution from Drude simulations of amylose reflects a more diverse conformational ensemble. Only the most populated Drude cluster retained a semihelical structure (FigureB), whereas the second and third clusters with the highest population adopted compact conformations, in which the ends of the amylose chain were closer together (Supporting Information, Figure S5). These findings are consistent with the broader conformational sampling observed in the RMSD profiles. Overall, the RMSD analysis and the subsequent clustering suggest that the inclusion of polarization effects leads to a broader exploration of the conformational landscape of amylose compared to the CHARMM FF. This increased plasticity may be related to the explicit treatment of electronic polarization in the Drude model, which could allow for more sensitive interactions within the system.
Representative structures of all clusters identified for amylose and cellulose simulations. (A) Amylose with the CHARMM FF, (B) amylose with the Drude FF, (C) cellulose with the CHARMM FF, and (D) cellulose with the Drude FF.
In the case of cellulose, the RMSD values produced by both the CHARMM and Drude FFs were more consistent than in the case of amylose (Supporting Information, Figure S6). Across all three replicates, the cellulose systems, whether simulated with CHARMM or Drude, maintained relatively low RMSD values throughout the 1-μs simulations, with only transient increases. For both CHARMM and Drude simulations of cellulose, the average RMSD was 4 ± 1 Å. In contrast, the CHARMM simulations of amylose yielded an average RMSD of 7 ± 2 Å and the Drude simulations showed even higher value of 11 ± 3 Å. This result suggests that cellulose, irrespective of the FF used, retains a narrower conformational ensemble than amylose. The higher RMSD values observed for amylose reflect its greater conformational variability in dilute aqueous solution. That is, the α(1→4) linkages introduce flexibility that allows the chain to adopt diverse coil-like geometries. In contrast, the β(1→4) linkages in cellulose impose geometric constraints and promote stabilizing intrachain hydrogen bonding, limiting structural drift. Thus, even when modeled as single short oligomers, amylose exhibits the expected flexibility while cellulose retains its characteristic rigidity. This comparison provides a controlled benchmark for evaluating how CHARMM and Drude FF capture intrinsic conformational behavior. Consistent with the RMSD analysis, clustering of the cellulose trajectories with MDANCE revealed a narrower conformational landscape for both FFs. The clustering detected two distinct clusters for each FF, with similar occupancies. The top cluster from the CHARMM simulations contained 69% of total frames and that of the Drude simulations contained 65% of the frames. The representative structures from the CHARMM and Drude simulations are shown in the Supporting Information, Figures S7 and S8, respectively, and showed elongated conformations (FigureC,D). These findings further support the conclusion that cellulose maintains a relatively narrow conformational ensemble, with minimal differences between the CHARMM and Drude FFs. This outcome agrees with earlier simulations of 9-unit cellulose oligomers, which likewise reported elongated conformations throughout the trajectory.?
Root-mean-square fluctuation (RMSF) measures the flexibility of individual atoms or functional groups, offering valuable insights into the local dynamics of the molecule.? RMSF plots for unrestrained amylose in the simulations using CHARMM and Drude FFs are shown in FigureA, showing fluctuation levels of all 12 glucose residues. These values were calculated by averaging across the three replicate simulations. There is a clear difference in flexibility between the two FFs, with the Drude FF exhibiting systematically greater fluctuation along the entire amylose chain than the CHARMM FF. Both CHARMM and Drude systems showed the expected increase in fluctuation at the terminal residues, indicative of increased flexibility at the termini of the polysaccharide chains. Conversely, residues near the middle of the amylose chain were more rigid. However, even among these central residues, the Drude simulations consistently produced slightly greater RMSF values than the CHARMM simulations. Also, the standard deviations shown by error bars on the Drude points are larger than those for CHARMM, reflecting a greater heterogeneity in the conformational sampling with the polarizable Drude model. These observations suggest that the Drude FF allows for a greater range of conformational flexibility in amylose than CHARMM, likely due to the explicit treatment of polarization effects.
Structural analyses and their comparison between CHARMM and Drude FFs for unrestrained amylose. (A) Per-residue RMSF of the amylose chain. (B) An extended conformation of amylose chain in CHARMM and a more compact conformation of amylose in Drude system. Probability densities of the (C) end-to-end distance and (D) the radius of gyration, respectively.
End-to-end distance (R ee) is calculated as the separation between the two termini of a polymer chain.? We calculated the evolution of the average R ee over three replicates of CHARMM and Drude simulations. The CHARMM FF yielded systematically greater average R ee values compared to those of the Drude FF over the simulation time (Supporting Information, Figure S9A). Additionally, the R ee probability density (FigureC) shows that although simulations using CHARMM and Drude FFs sampled the same range of R ee, the CHARMM FF favored values around 35 Å, whereas the Drude FF featured a dominant peak near 10 Å. This difference indicates that the Drude FF samples more collapsed states of amylose, compared to the more extended conformations produced by the CHARMM FF (FigureB). This result further contextualizes the clustering results discussed above (Supporting Information, Figures S4 and S5), in that the collapsed structures produced by the Drude FF were disordered with proximal termini rather than being coiled in helical-like structures. The imperfect helices and collapsed states we observed align with previous MD analyses demonstrating that amylose helices are destabilized by band-flips and kinks in aqueous solution, yielding irregular hydrogen-bonding networks rather than perfectly ordered structures.?
The radius of gyration (R g) assesses the three-dimensional compactness of a molecule, representing the root-mean-square distance of its atoms from the center-of-mass.? As with our findings on R ee for the unrestrained amylose, the evolution of R g (Supporting Information, Figure S9B) and its probability density (FigureD) show the same range of R g values for CHARMM and Drude systems. However, their probability densities demonstrate that the Drude model exhibits a sharper and higher peak at a lower R g value compared to the CHARMM model. This result implies that the Drude model predominantly adopted more compact conformations. In contrast, the distribution of R g in CHARMM system is more diffuse with lower peak in the probability density graph (FigureD), reflecting a wider range of amylose conformations throughout the simulations. The similarity in the probability density graphs of the R g and R ee suggests that the observed compaction of amylose in the Drude simulations is primarily due to the two terminal ends of the chain approaching each other, rather than compaction of inner residues. We performed the same analysis on the cellulose simulations to assess the consistency of conformational behavior across FFs. As shown in FigureA, RMSF calculations of the single cellulose chain in the CHARMM and Drude simulations reveals nearly identical flexibility profiles in both CHARMM and Drude simulations. The terminal residues of the cellulose chain were more flexible than the residues in the middle of the chain, which is an expected outcome due to their exposure and lack of neighboring interactions. The R ee of the cellulose chain in CHARMM and Drude systems were essentially indistinguishable (FigureB and Supporting Information, Figure S10A) indicating that both FFs capture a comparable degree of extension along the chain. The same holds true for the R g, which also showed minimal differences between the two simulations (FigureC and Supporting Information, Figure S10B).
Structural analyses and their comparison between CHARMM and Drude FFs for cellulose. (A) Per-residue RMSF of the cellulose chain. Probability densities of the (B) R ee and (C) the R g, respectively.
Amylose and cellulose likely respond very differently in the Drude polarizable FF because of their distinct chemistries. Amylose, an α(1→4) glucan, is relatively flexible and forms helix-like coils in water with only “imperfect” intrachain hydrogen bonding.? Its helical cavity is largely hydrophobic and encourages helix formation and internal packing,? so when electronic polarization is included, induced dipoles strengthen the effective attractions among the glucose units. As a result, our Drude simulations showed amylose folding into tighter conformationsmanifesting as a lower R g and shorter R ee than in the corresponding CHARMM simulations. In contrast, cellulose is a β(1→4) glucan whose chain is intrinsically stiff and equipped with equatorial −OH groups that form a regular pattern of intrachain hydrogen bonds.? As such, even with explicit electronic polarization, the balance of forces does not change substantially. Indeed, the probability densities of RMSD, R g, and R ee for cellulose are essentially the same with the Drude and CHARMM FFs, and clustering yielded similar (extended) conformations with both FFs. In short, the intrinsic flexibility and hydrophobic character of amylose make its conformational ensemble more sensitive to the induced-dipole forces in the Drude model and thus it collapses to a more compact form. In contrast, the rigidity from intrachain hydrogen bonding in cellulose restricts its shape so that polarization effects make little difference in its observed properties. Therefore, use of a polarizable model when simulating amylose and its binding to hydrophobic compounds may be more appropriate.
The overall conformation of an oligosaccharide is primarily influenced by the Φ and Ψ orientations of its torsion angles between glycosyl residues (FigureA). ?,? Specifically, the conformational sampling of the acetal linkages are described by the rotation around the C1–O (Φ angle) and O–C4 (Ψ angle) bonds.? Characterizing these glycosidic torsion angles is essential for describing the structural flexibility and functional properties of oligosaccharides in biological systems. We calculated the Φ and Ψ values for the linkages between all 12 sugar units in unrestrained amylose and cellulose chains under both FFs (Supporting Information, Figure S11A–D). In the case of amylose, the probability density of Φ and Ψ angles reveal distinct conformational preferences in each FF. Both CHARMM and Drude models show a dominant peak around Φ = 100°. In contrast, the Ψ distribution in CHARMM is centered at −142°, whereas Drude samples a wider range around −135°. In addition, both force fields sampled secondary Ψ states, with CHARMM showing a peak near 80° and Drude near 75°. The broader range of Ψ values observed in the Drude simulations suggests greater flexibility around the glycosidic linkages and, consequently, increased conformational flexibility of amylose under the Drude FF. This observation is consistent with our clustering analysis, which revealed greater conformational variability in amylose with the Drude FF compared to CHARMM. To further illustrate these differences, we constructed two-dimensional free-energy surfaces of the Φ and Ψ torsions for both force fields (FigureB–E). These plots reveal two major basins in each force field, with one basin occupying a larger area and corresponding to lower free-energy values. In the Drude model, both basins are noticeably broader than CHARMM, indicating enhanced conformational sampling and increased flexibility. The dominant basin in CHARMM is defined near (100°, −142°), whereas in Drude it spans a wider region centered near (100°, −135°), with an additional shallow minimum near (100°, 75°). These differences underscore the impact of explicit polarization in the Drude model, which allows amylose to access a more diverse set of low-energy conformations.
Characterization of glycosidic torsion angles in amylose and cellulose. (A) Schematic of the Φ and Ψ angles around glycosidic bonds for a polysaccharide. Free energy surfaces of of Φ and Ψ for unrestrained amylose in (B) CHARMM and (C) Drude. Free energy surfaces of Φ and Ψ for cellulose in (D) CHARMM and (E) Drude.
For cellulose, the probability density plots of dihedral angles (Supporting Information, Figure S11C–D) show that the Φ distributions are closely aligned with both force fields exhibiting a dominant peak around −80°. For the Ψ angle, the distributions for CHARMM and Drude are largely centered around −125°; however, CHARMM shows a more sharply defined peak at this value and a secondary peak around ∼90°, whereas in the Drude simulation this secondary population is shifted toward values closer to 60°. The two-dimensional Φ–Ψ free-energy surfaces (FigureD–E) further illustrate the similarity between the CHARMM and Drude force fields for cellulose. Both maps show the same major low-energy basin centered around Φ ≈ −80° and Ψ ≈ −125°, consistent with the dominant populations observed in the probability distributions. A secondary, higher-energy basin appears at positive Ψ values for both models, though CHARMM displays a slightly more localized minimum compared to the broader, more diffuse basin observed with Drude. Overall, the location and shape of the free-energy minima are highly comparable between CHARMM and Drude, reinforcing that both force fields stabilize similar glycosidic conformational ensembles for cellulose.
Another interesting aspect of the dynamics of these polysaccharides is sugar puckering. We used the Cremer-Pople convention? to analyze and quantify the deviations from planarity in each of the glucopyranose rings. In the Cremer-Pople method, the puckering geometry with respect to a mean plane is described by θ and φ parameters. The θ parameter is the polar angle within a spherical coordinate system and indicates the general type of puckering conformation such as chair or boat. The φ parameter serves as the phase angle, and it defines the specific type of puckering deformation that occurs within the ring. It captures the nature of the puckering, such as whether the puckering forms a boat conformation, twist-boat conformation, or other typical conformations observed within cyclic compounds. The important ^4^ C 1 or ^1^ C 4 chair conformations are indicated by θ = 0° and θ = 180°, respectively. When θ = 90°, the ring samples boat (φ = 0, 60, 120, 180, 240, or 300°) or twist-boat conformations (φ = 30, 90, 150, 210, 270, or 330°).
With the CHARMM FF, all rings in amylose exhibited θ values near 0°, indicating a preference for the ^4^ C 1 chair conformation (FigureA). With the Drude FF, while predominantly sampling the ^4^ C 1 chair conformation, amylose also sampled the ^1^ C 4 chair and boat conformations (FigureA). The sampling of both ^4^ C 1 and ^1^ C 4 conformations in the Drude simulations reveals greater flexibility of ring structures compared to the CHARMM simulations. Experimentally, it has been demonstrated that the ^4^ C 1 chair-puckering conformation is favored.? This conclusion was also reached in another comparative study by Chytra et al. on the CHARMM and Drude FFs in the context of α-d-glucose monosaccharides.? Their work highlighted differential conformational preferences using CHARMM and Drude FFs. Specifically, the analysis determined that simulations with the CHARMM FF samples predominantly the ^4^ C 1 conformation but simulations with the Drude FF exhibited greater exploration of midconformations and the ^1^ C 4 conformation.? Our results confirm that these puckering dynamics extend to polysaccharide systems and emphasize the conformational flexibility that exists with the Drude FF, which may arise due to lower energy barriers to structural interconversion. Importantly, experimental evidence also suggests that pyranoses such as glucose and GlcNAc are not exclusively rigid ^4^ C 1 chairs, but can transiently access nonchair geometries. Atomic force microscopy studies of glucose polymers inferred the presence of non-^4^ C 1 conformations. ?−? ? Thus, the broader sampling observed with the Drude FF is consistent with experimentally supported flexibility, even though the ^4^ C 1 chair remains the predominant conformation.
Probability densities of θ values using the Cremer-Pople convention for (A) amylose and (B) cellulose chains.
We also calculated the θ values for the cellulose chain in the CHARMM and Drude systems (FigureB). The probability densities of θ values in both FFs exhibited a prominent peak near 0° and the Drude FF produced a very small population at 90°. Thus, the conformational plasticity with the Drude FF persists for cellulose, indicated by the presence of boat and twist-boat conformations, but to a smaller degree than in the case of amylose. Nevertheless, both the CHARMM and Drude FFs consistently favor the ^4^ C 1 chair conformation, the expected geometry for this system.
Hydration of Single-Stranded Amylose and Cellulose
The radial distribution function (RDF) provides valuable insights into the spatial distribution of particles around a reference particle, shedding light on liquid structure, and thus the intermolecular interactions within a system. In aqueous solution, hydrogen bonds, are critically important for determining the structure of water and its interactions with solutes. We calculated the RDF for water around each of the oxygen atoms within sugar units in both CHARMM and Drude FF systems and analyzed the total number of intramolecular hydrogen bonds within the polysaccharide chain as well as those formed between the chain and water molecules for all unrestrained and restrained amylose systems.
The resulting RDF plots for unrestrained amylose reveal distinctive hydration characteristics between the CHARMM and Drude systems. Notably, in both FFs, hydroxyl oxygen atoms O2, O3, and O6 were more hydrated than O4 and O5, the atoms involved in the glycosidic linkage and functioning solely as hydrogen bond acceptors (FigureA,B). Consistently, these hydroxyl groups also formed a greater number of hydrogen bonds with water molecules (FigureC,D) as a result of their solvent accessibility and their dual role as hydrogen bond donors and acceptors. Our findings align well with the study by Khatami et al.,? who used MD simulations to investigate hydrogen bonding patterns of single amylose chains in water. They found that amylose chains formed a greater number of hydrogen bonds with water molecules at the O2, O3, and O6 positions. Additionally, their study highlights the dynamic nature of hydrogen bond formation and breakage, aligning with the fluctuating behavior we observed in different FFs.?
Radial distribution functions for unrestrained amylose and cellulose. Amylose RDFs are shown for the CHARMM and Drude FFs in panels (A) and (B), respectively. Cellulose RDFs with the CHARMM and Drude FFs are shown in panels (C) and (D), respectively.
Hydrogen bonding in CHARMM and Drude systems for unrestrained amylose. Total number of (A) amylose-water (intermolecular) and (B) amylose-amylose (intramolecular) hydrogen bonds with both CHARMM and Drude FFs. Panels (C) and (D) show the number of hydrogen bonds formed between each of the oxygen atoms in the sugar rings and water in the CHARMM and Drude simulations, respectively. Panels (E) and (F) show the number of polysaccharide intramolecular hydrogen bonds in the CHARMM and Drude simulations, respectively.
Within the CHARMM system, the oxygen atoms O2, O3, and O6 exhibited much greater hydration (FigureA,B) and a greater number of hydrogen bonds with water (FigureC,D) compared to their counterparts in the Drude system, suggesting a stronger attraction of water molecules to these hydroxyl groups. These hydrogen bonds were relatively consistent throughout the simulation, reflecting stable hydration of hydroxyl groups. In contrast, the Drude FF produced more intramolecular hydrogen bonds (FigureB), and these interactions fluctuated more frequently during the simulations, indicating dynamic exchange events. FigureE,F illustrate the intramolecular hydrogen bonding for individual oxygens within the amylose chain using CHARMM and Drude models. These results indicate that the CHARMM model tends to maintain stable intramolecular hydrogen bonds, contributing to a more rigid structure. In contrast, the Drude model’s dynamic behavior highlights greater structural flexibility, allowing the amylose chain to explore a wider range of conformations.
It is worth noting that the O1 atom, situated at the reducing end of the chain, was prominently hydrated and formed hydrogen bonds with water molecules, consistent with its exposed position allowing interaction without steric hindrance from neighboring glucose units. As expected, O4 consistently exhibited the lowest hydration peak in the RDF and the fewest hydrogen bonds with water, reflecting its role as the glycosidic ether oxygen with limited hydrogen-bonding capacity. Hydration and hydrogen-bond behavior were similar in the restrained systems of amylose (Supporting Information, Figure S12 and S13), with O2, O3, and O6 maintaining higher RDF peaks and a greater frequency of hydrogen bonds than O4 and O5. The reduced interaction of O4 and O5 with water is consistent with their ether character, in contrast to the hydroxyl oxygens that readily form hydrogen bonds with solvent.
We also calculated RDFs and hydrogen bonding patterns for individual oxygen atoms in cellulose chains (FigureC,D and ?). The RDF profiles and H-bond analysis yielded similar results as in the amylose chain. In CHARMM, hydroxyl oxygen atoms O2, O3, and O6 were more hydrated and formed more hydrogen bonds with water compared to O4 and O5, indicating a greater propensity for solvent interaction. For intrachain hydrogen bonds, O3 and O5 were especially prominent in cellulose, differing from amylose, where O2, O3, and O6 dominated intramolecular interactions. O4 again exhibited the lowest hydration and fewest hydrogen bonds. These results are consistent with prior studies by Shen and Gnanakaran demonstrating the stability of O3H3···O5 intrachain hydrogen bonds in cellulose.?
Hydrogen bonding in CHARMM and Drude systems for cellulose. Total number of (A) cellulose-water (intermolecular) and (B) cellulose-cellulose (intramolecular) hydrogen bonds with both CHARMM and Drude FFs. Panels (C) and (D) show the number of hydrogen bonds formed between each of the oxygen atoms in the sugar rings and water in the CHARMM and Drude simulations, respectively. Panels (E) and (F) show the number of intramolecular polysaccharide hydrogen bonds in the CHARMM and Drude simulations, respectively.
Taken together, the combined RDF and hydrogen-bond analysis suggests that while the hydration properties of amylose and cellulose are broadly similar, there are subtle differences as a function of inclusion of electronic polarization in the simulation. CHARMM simulations indicated a greater number of water molecules around hydroxyl groups and more stable hydrogen bonds with water, whereas Drude simulations favored intramolecular hydrogen bonding and exhibited more dynamic exchange. This difference may be explained by the anisotropic treatment of electronic polarization and the inclusion of “lone pair” virtual sites on oxygen atoms in the Drude FF, which reflect the directionality of electrostatics intrinsic to these functional groups and manifest as slightly reduced hydration compared to the additive CHARMM FF.
Electrostatic Properties of Single-Stranded Amylose and Cellulose
and Hydrating Waters
Figure illustrates the probability densities of dipole moments for the sugar rings in amylose and cellulose chains using CHARMM and Drude polarizable FFs. Given the flexibility of the sugar rings (puckering) and rotation of hydroxyl groups, these dipole moments reflect contributions from both nuclear geometry and electronic plasticity. For both amylose and cellulose, the dipole moment distributions with the CHARMM FF are generally unimodal and have a slightly wider range compared to the values produced by the Drude FF. For amylose, the average dipole moment was 5 ± 2 D in both the CHARMM and Drude simulations. For cellulose, the CHARMM simulations yielded an average dipole moment of 4 ± 2 D, while the Drude simulations produced a value of 3 ± 1 D. This outcome indicates that the sugar rings in the CHARMM system generally possess higher dipole moments. This property can be attributed to the intrinsic mean-field charge assignment in additive FFs, which is designed to reflect average properties in an aqueous medium. Conversely, the dipole moment distributions from the Drude FF simulations are slightly narrower with sharper peaks, indicating a different, and responsive representation of dipole moments. The polarizable nature of the Drude FF captures real-time electronic fluctuations, in response to the local electrostatic environment and resulting in a dynamic hydrogen bonding network.
Probability densities of sugar unit dipole moments in (A) amylose and (B) cellulose.
In addition to investigating the dipole moments of constituent glucose units, it is also of interest to calculate dipole moments of water molecules that directly interact with the polysaccharides. By assessing variations in water dipole moments, we sought to quantify the influence of induced polarization by carbohydrates on the surrounding solvent structures and solvation dynamics.? Our analysis focused on calculating the dipole moments of water molecules confined within the restricted chain of amylose (amylose restrained 1, which set Φ = 103° and Ψ = 115°). We selected this structure due to the presence of clearly defined cavity within its helical configuration. Prior to this analysis, we aligned the helical axis along the x-axis to ensure a consistent reference frame for both dipole moment and electric field vector analysis (FigureA). We investigated the components of the molecular dipole moments because they allow for insight into whether or not there is a preferential alignment of these molecules in the cavity.
Internal hydration of amylose in a helical configuration. (A) Side and front views of the amylose helical structure aligned with x-axis. Blue circles show the water molecules confined inside the helix. Probability densities of water molecule dipole moment components for (B) CHARMM and (C) Drude simulations.
In simulations using the CHARMM FF, we observed a distinct peak at ∼2 D in the probability densities of water dipole moments, indicating a preference for the alignment of water dipole moments along the + x-axis (FigureB). However, the distribution displayed a near-flattened profile along the y- and z-axis, suggesting no preferential orientations along these axes. It is important to note that given the rigid geometry of the TIP3P water model, any apparent change in dipole moment here simply reflects rotational motion of the water such that it is oriented differently along the Cartesian axes. In contrast, simulations employing the Drude FF produced a bimodal distribution of water dipole moments, indicating that the waters were aligned in both directions, with some disorder in the middle (FigureC). In the Drude simulations, the SWM4-HLJ water model? has a variable dipole moment given that it is explicitly polarizable. Across all three axes, the dipole vector components of water molecules confined within the amylose helix exhibited a range between −3 and 3 D in the Drude system, wider than that observed in CHARMM simulations. Whereas the water dipole moments from the CHARMM simulations had a hard truncation at ±2.34 D (the fixed value of TIP3P, FigureB), the gradual decay of these values in the Drude simulations with SWM4-HLJ (FigureC) reflects the impact of explicit polarization. Given the bulk dipole moment value of ∼2.58 D for SWM4-HLJ,? our dipole vector component analysis suggests that along the x-axis, the water molecules are slightly polarized at the terminal ends of the amylose helix. The distinct peaks at the ends of the amylose helix reveal another difference between the two FFs compared here, in that the CHARMM FF (in conjunction with the TIP3P model) leads to unidirectional alignment of water molecules, while the Drude FF with the SWM4-HLJ model leads to bidirectional alignment of water. The impact of this predicted phenomenon will require additional experimental and theoretical work to validate and understand.
Properties of the First Hydration Shell around Single-Stranded
Amylose and Cellulose
The first hydration shell around a biomolecule often has different structural, dynamic, and thermodynamic properties from bulk water. In MD simulations, therefore polarizable FFs may serve as a better model of these properties.? Recent work on monosaccharides supports this view, as the polarizable AMOEBA force field produces hydration shells that are denser compared to additive models, and, importantly, produces solution behavior and thermodynamic properties in closer agreement with experiment.? Moreover, AMOEBA eliminates spurious sugar aggregation commonly observed with fixed-charge force fields, suggesting a better balance of intra- and intermolecular forces. We examined the dipole moments of water molecules in the first hydration shell of amylose to determine the impact of electronic polarization and to understand if the observed behavior of the AMOEBA FF is also produced by the Drude model. As in the analysis of cavity-bound water molecules, the amylose chain was aligned along the x-axis, and dipole moments were decomposed into their x-, y-, and z-components.
In CHARMM simulations, the distribution of the x-axis dipole moment component shows a distinct peak, indicating preferential water molecule alignment along the x-axis in which the amylose chain is oriented (FigureA). Given that the distributions along the y- and z-axes are rather flat, the waters exhibited no preferred orientation in these dimensions. Additionally, as discussed above, in the CHARMM simulations, water molecules participated in a greater number of hydrogen bonds with amylose. As a result, the water molecules in the CHARMM simulations are likely somewhat restricted because they are engaged in relatively strong hydrogen bonds with surface-accessible hydroxyl groups on amylose.
Probability densities of the dipole moments of water molecules in the first hydration shell of amylose chain in (A) CHARMM and (B) Drude simulations, respectively.
The dipole properties of first-shell water molecules in the Drude simulations were similar to those of the CHARMM simulations (FigureB). The y- and z-axis components of the water dipole moments were relatively flat, and along the x-axis there was a slight alignment effect along the +x-axis. As with the internal water molecules, the SWM4-HLJ water molecules manifested a wider distribution of dipole moments (from −3 to 3 D) along the three axes in the Drude simulations. Thus, we conclude that there is a subtle, but potentially important, difference in water alignment of internal water molecules in the amylose helix with the different FFs, but the first-shell water molecules behaved similarly.
Electric Fields inside the Amylose Helix
We calculated the electric field (E⃗) exerted by the amylose chain along its axis within the helical structure using the TUPÃ software.? E⃗ calculations were performed at 24 equidistant points along the axis of the amylose chain (FigureA). In this convention, point 1 is closest to the reducing end of the chain and point 24 is closest to the nonreducing end. E⃗ was calculated according to eq:
where E⃗(x,y,z) is the electric field vector at the point (x,y,z), N is the number of atoms considered in the calculation, Q _ i _ is the charge of atom i, r is the distance between the point of interest and atom i, r̂_(x,y,z)_ is the unit vector pointing from atom i to the point (x,y,z), and ε_0_ is the vacuum permittivity constant.?
Electric field calculations along the axes of the “restrained 1” and “restrained 2” amylose structures. (A) A schematic representation of helical amylose chain and the points along the axis where E⃗ values were calculated. E⃗ values for (B) “restrained 1” and (C) “restrained 2” amylose structures. Error bars correspond to the 95% confidence interval.
The calculations were performed on the trajectories obtained from both restrained amylose systems, allowing for a comparison of the electric fields acting within a helical structure as a function of its solvent accessibility. For the single-chain system, the entire structure was divided into two blocks, each containing 6 glucose units. The dynamic coordinates of each of 24 points along the amylose chain axis were determined by calculating the minimum and maximum values of each set of 6 residues. This approach ensured that the coordinates were dynamically adjusted for each point along the helix. E⃗ values were averaged across all frames in the simulation.
E⃗ values calculated from both CHARMM and Drude “restrained 1” and “restrained 2” amylose simulations are shown in FigureB, C. A pronounced oscillatory pattern emerged in the E⃗ values for both FFs, signifying periodic variations in electrostatic interactions at different points along the axis. In the case of the “restrained 1” structure, which has a wider central cavity, E⃗ values with both FFs varied between 30–60 MV/cm (FigureB) and steadily declined to ∼20 MV/cm toward the nonreducing end of the amylose chain. Periodicity and symmetry of the E⃗ values in the CHARMM simulation of the “restrained 2” structure, which has a very narrow central cavity, were more pronounced and E⃗ values were higher than in the “restrained 1” simulation, 50–150 MV/cm (FigureC). Thus, it appears that compaction of the amylose structure leads to larger internal E⃗ values. This outcome is sensible given the distance-dependence of electric fields; the strength of the field decays linearly as a function of distance between the charged atom and the point in space (eq). In the Drude simulation of the “restrained 2” structure, some periodicity was apparent but E⃗ values were higher toward the nonreducing end of the chain than in the case of CHARMM, reaching ∼200 MV/cm before reducing to ∼50 MV/cm at the final point along the axis.
E⃗ values differed in their variation between the two FFs. The E⃗ values from the CHARMM simulation fluctuated less, leading to smaller error bars, and indicating of a less variable electrostatic environment along the amylose chain. Conversely, the E⃗ values from the Drude simulation exhibited larger error bars and greater variability, reflecting the dynamic and flexible nature of the polarizable Drude FF given the ability of induced dipoles to respond to changes in the local environment.
Double-Stranded Amylose
In addition to investigating single-stranded amylose chains, we extended our study to explore the behavior of double-helical structures of amylose. In all the three simulation replicates using CHARMM and Drude FFs, we observed remarkable robustness in the double-stranded helical structure, which remained largely intact throughout the 1-μs simulation duration. Although transient unzipping events occurred, the helical structure promptly reformed, indicating a high degree of stability. FigureA illustrates the RMSD evolution of the double-helical amylose employing the CHARMM and Drude FFs. The observed stability of the double-helical structure is consistent with experimental studies reported previously. Techniques like NMR spectroscopy have confirmed the presence of a double helix structure of amylose in water, demonstrating that these duplexes are thermodynamically stable entities at sufficient chain lengths (greater than 12 glucose units) and concentrations, with stability increasing linearly with chain length and driven primarily by favorable enthalpic interactions.?
Structural characterization of the double-helical amylose structure. (A) Average RMSD over three replicates, with lighter shading indicating the standard deviation across the replicates. Probability distributions of (B) R ee, (C) R g, (D) Φ, (E) Ψ, and (F) the θ puckering parameter for CHARMM and Drude simulations. Experimental values are indicated with triangles in panels D and E.
We also computed the R ee and R g for the double-helical amylose for both the CHARMM and Drude FFs (FigureB,C). The R ee distribution reveals that the CHARMM model exhibits a sharp peak around 43 Å, indicating a consistent and rigid double helix structure. In contrast, the Drude model shows a broader distribution, with peaks around 19 Å and 29 Å, suggesting greater variability and flexibility in the helical structure. Similarly, the R g distribution (FigureC) highlights a sharp peak around 14 Å for the CHARMM model, emphasizing a stable conformation. The Drude model, however, demonstrates a broader distribution, reflecting increased flexibility and dynamic behavior. These differences underscore the CHARMM FF’s propensity for rigidity, while the Drude FF predicts a more flexible structure. The transient partial unzipping observed in Drude simulations suggests that amylose duplexes of chain length 12 are at the threshold of stability and can partially open at their termini. This observation aligns with experimental findings showing that at least 12 glucose units are required for duplex stability, which further improves with longer chains beyond a degree of polymerization of 12.?
We analyzed the Φ and Ψ dihedral angles for both chains of the double-helical amylose structure to assess whether the observed global structural stability was maintained at the level of glycosidic geometries. In the CHARMM simulations, we observed small fluctuations in the Φ and Ψ values throughout the simulations. Consequently, the corresponding probability densities exhibited sharper peaks compared to the Drude simulations, indicating a greater degree of structural rigidity (FigureD,E). The Φ sampling with the CHARMM FF aligns closely with the experimental value (FigureD), but for Ψ sampling, the CHARMM FF produced two peaks (FigureE). The dominant peak coincides with the experimental value but the second peak deviates substantially. Conversely, in the Drude simulations, both Φ and Ψ agreed well with the experimental values, though the Φ sampling deviated slightly from the experimental value.
Additionally, employing the Cremer-Pople convention,? we calculated θ values to evaluate ring puckering. In the CHARMM system, θ values clustered near zero, indicating predominantly ^4^ C 1 chair conformations for the constituent glucose rings. The ^4^ C 1 chair conformation was also dominant in the Drude simulations, although with the polarizable model, some ^1^ C 4 chair conformations were sampled (FigureF). This subtle shift in the distribution of θ values in the Drude simulations underscores a slight increase flexibility and susceptibility to chair conformational transitions in the Drude-modeled amylose, similar to what we observed in the case of single-chain amylose. The double-helical nature of the duplex structure likely restricts conformational fluctuations to some extent, as the emergence of ^1^ C 4 chair conformations in the single-stranded amylose simulations was more pronounced (Figure).
Finally, we examined the hydrogen bonding between the two parallel chains as well as within each individual chain in the double helix. We observed that a greater number of hydrogen bonds formed between the two chains than within each individual chain using the CHARMM FF (FigureA). This observation suggests that interchain hydrogen bonding stabilizes the double-helix structure. In the Drude system, the opposite was true, such that the hydrogen bonds in a single chain were greater than those between the chains (FigureB). This behavior suggests weaker interchain interaction in Drude system, which can lead to the greater flexibility of the double helix that we have described above. We also compared the number of hydrogen bonds between each chain and water. The CHARMM FF produced a greater number of amylose-water hydrogen bonds than the Drude FF (FigureC,D). Each of these values was systematically lower than in the case of single-chain amylose with water (FigureA), likely due to some hydroxyl groups being engaged in interchain hydrogen bonds that stabilize the double helix. Finally, we counted the hydrogen bonds formed between the amylose chain and water molecules as a function of each hydroxyl group (Supporting Information, Figure S14). With the CHARMM FF, O2 and O3 formed a greater number of hydrogen bonds with water molecules than with the Drude FF, whereas O6 formed an equal number of hydrogen bonds with both the CHARMM and Drude FFs. These findings reflect the impact of FF selection on the balance of intermolecular interactions and stability of biomolecular structures. As in the single-chain amylose systems, the properties of double-helical amylose appear to be subtly dependent on softer interactions with the Drude FF compared to the additive CHARMM FF.
Intermolecular and intramolecular hydrogen bonding in the double-helical amylose system. The number of intermolecular and intramolecular hydrogen bonds in (A) CHARMM and (B) Drude simulations, shown as an average over three replicates. The lighter shading indicates the standard deviation across the three replicates. Hydrogen bonds formed between each amylose chain and water are shown for (C) CHARMM and (D) Drude simulations.
Conclusions
In this study, we investigated the structural, hydration, and electronic properties of amylose and cellulose through MD simulations using CHARMM and Drude FFs. Our findings highlight the importance of the FF selection in modeling the conformational dynamics, hydration behavior, and electrostatic properties of the amylose chains in aqueous solution. Structural analyses revealed that amylose is more flexible than cellulose, and that the Drude FF samples a larger conformational space, likely due to its directional electrostatics (modeled via induced dipoles) and lower intrinsic energy barriers. End-to-end distance, radius of gyration, glycosidic torsion angles, and ring puckering corroborated that amylose adopts a more flexible conformation using the Drude FF, but cellulose adopts a relatively rigid conformation with both FFs.
Hydration analysis revealed distinct differences in water interactions between the two FFs. The CHARMM FF produced more pronounced hydration of hydroxyl groups and stronger hydrogen bonding networks with water molecules, whereas the Drude FF preferred more intramolecular hydrogen bonding in amylose. These distinctions emphasize the importance of polarization effects in governing solvation and hydrogen bond dynamics in carbohydrate systems. Additional theoretical and experimental work should be conducted to test these predictions.
Electronic property calculations provided insight into the role played by polarization in carbohydrate-water interactions. The Drude model captured variable dipole moments, giving rise to a more dynamic electrostatic environment. Comparing dipole moments of water confined inside the amylose helix to water present in the first hydration shell indicated stronger directionality with CHARMM, but with Drude, water molecules were more dynamic and heterogeneous in their orientation. Calculation of the electric field along the amylose helical axis also revealed periodic variation in electrostatic forces.
Our extension of simulations of double-helical amylose structures confirmed the stability of the helical structure, with some FF-dependent variations in transient disorder. CHARMM simulations produced more interchain hydrogen bonding, therefore favoring structural stability, whereas Drude simulations favored increased intrachain hydrogen bonding, suggesting greater flexibility. The structural stability of the double helix was generally preserved for both FFs, consistent with experimental findings.
Overall, this comparative FF study provides a systematic, molecular-level understanding of amylose and cellulose in water. The results emphasize the significance of selecting an appropriate FF based on the target properties in carbohydrate research. Future studies using enhanced sampling techniques may provide additional insights into these complex biomolecules and potentially enable improved FF parametrization and use in biomaterials engineering, drug delivery, and carbohydrate-based therapeutics.
Supplementary Material
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Horton, D. The development of carbohydrate chemistry and biology. In Carbohydrate chemistry, biology and medical applications; Elsevier, 2008; pp. 1–28.
- 2Lepenies B.Yin J.Seeberger P. H.Applications of synthetic carbohydrates to chemical biology Curr. Opin. Chem. Biol.201014340441110.1016/j.cbpa.2010.02.01620227905 · doi ↗ · pubmed ↗
- 3Jelinek R.Kolusheva S.Carbohydrate biosensors Chem. Rev.2004104125987601610.1021/cr 030028415584694 · doi ↗ · pubmed ↗
- 4Werz D. B.Seeberger P. H.Carbohydrates as the next frontier in pharmaceutical research Chem. Eur. J.200511113194320610.1002/chem.20050002515798968 · doi ↗ · pubmed ↗
- 5Obiro W. C.Sinha Ray S.Emmambux M. N.V-amylose structural characteristics, methods of preparation, significance, and potential applications Food Rev. Int.201228441243810.1080/87559129.2012.660718 · doi ↗
- 6Li Y.Liu P.Ma C.Zhang N.Shang X.Wang L.Xie F.Structural disorganization and chain aggregation of high-amylose starch in different chloride salt solutions ACS Sustainable Chem. Eng.20208124838484710.1021/acssuschemeng.9b 07726 · doi ↗
- 7Li W.Corke H.Beta T.Kinetics of hydrolysis and changes in amylose content during preparation of microcrystalline starch from high-amylose maize starches Carbohydr. Polym.200769239840510.1016/j.carbpol.2006.12.022 · doi ↗
- 8Cosgrove D. J.Growth of the plant cell wall Nat. Rev. Mol. Cell Biol.200561185086110.1038/nrm 174616261190 · doi ↗ · pubmed ↗