Cellulose Iβ behaviors in non-solvent liquid media: Molecular dynamic simulations

Kong, Y., Fu, S., Yang, X., Leu, S-Y., and Hu, C. (2023). “Cellulose Iβ behaviors in non-solvent liquid media: Molecular dynamic simulations,” BioResources 18(4), 8223-8248.

Abstract

The structural changes of cellulose in non-solvent liquid media can provide insights into the high-value utilization of cellulose. This study includes molecular dynamics simulations of 36-chain cellulose Iβ microfibril model (Iβ-MF) behavior in 16 non-solvent liquids with different polarities at room temperature using two carbohydrate force fields (CHARMM36, GLYCAM06). Iβ-MF in CHARMM36 retains more than 70% of the tg conformation in 16 liquids, and the retention of the tg conformation increased with decreasing liquid polarity. Liquid polarity can affect the hydroxymethyl conformation of cellulose, which is only an appearance, and the real driving force behind is the electrostatic interaction between liquid molecules and cellulose. Furthermore, changing the 1,4 electrostatic scaling factor of GLYCAM06 can effectively affect the structural convergence of Iβ-MF. The Iβ-MF forms an alternating layer structure in the gg/gt conformation in a medium to high polarity non-solvent liquid, while the model undergoes untwisting. Model untwisting is inextricably linked to the degree of alternate layer structure formation. This paper provides a theoretical basis for the molecular study of nanocellulose structures from an energy-structure-property perspective.

Download PDF

Full Article

Cellulose Iβ Behaviors in Non-solvent Liquid Media: Molecular Dynamic Simulations

Yi Kong,^a Shiyu Fu,^a,b,* Xuedi Yang,^b Shao-Yuan Leu,^c and Chuanshuang Hu ^d

DOI: 10.15376/biores.18.4.8223-8248

Keywords: Cellulose; Liquid polarity; Molecular dynamics; Structural transformation

Contact information: a: State Key Laboratory of Pulp and Paper Engineering, South China University of Technology, Guangzhou 510640, PR China; b: South China University of Technology-Zhuhai Institute of Modern Industrial Innovation, Zhuhai 519175, PR China; c: Department of Civil and Environmental Engineering, The Hong Kong Polytechnic University, Hong Kong SAR 999077, PR China; d: College of Materials and Energy, South China Agricultural University, Guangzhou 510640, PR China;

* Corresponding author: shyfu@scut.edu.cn

INTRODUCTION

Cellulose, the main component of plant cell walls, is nature’s most widely distributed and abundant homoglycan (Matthews et al. 2006; Zhao et al. 2013; Zhang et al. 2019; Zhou et al. 2021). Cellulose consists of repeating β-D-glucopyranose units linked by β (1-4) glycosidic bonds. These pyranose rings are in a chair-like conformation with the hydroxyl group in the equatorial position (Kamide 2005). In nature, wood cellulose chains’ polymerization (DP) is about 10,000 pyranose units, and cotton cellulose is about 15,000. Natural cellulose is defined as cellulose I, which is known from ¹³C CP/MAs NMR spectra to exist in two different forms called cellulose Iα and Iβ (Horii et al. 1987). The main difference between cellulose Iα and Iβ is the stacked arrangement of the hydrogen-bonded layers. The layers in Iα (P1 symmetry) are permanently displaced +c/4 along the c-axis, while the layers in cellulose Iβ (P2₁ symmetry) are alternately displaced at +c/4 and -c/4 (center and origin chains) (Gardner and Blackwell 1974). Cellulose I from naturally occurring lower plants (e.g., algae and bacteria) is rich in cellulose Iα. In contrast, cellulose Iβ, mainly from higher plants (e.g., cotton and wood), is the most studied and widely used type of cellulose (Imai and Sugiyama 1998; Nishiyama et al. 2008).

Hydrogen bonding is an essential component of the cellulose crystal structure. For cellulose I, hydrogen bonding within the cellulose chains allows a linear arrangement of cellulose, and two adjacent cellulose chains are bound together mainly by hydrogen bonding, forming cellulose sheets. Stacked sheets, conversely, are thought to have no hydrogen bonding interactions between them but rather form cellulose microfibril crystals through van der Waals interactions (Zhou et al. 2021). Notably, the orientation of the C6 hydroxymethyl group will highly affect the hydrogen bonding and chain conformation. This pyranose ring substituent (C6 hydroxymethyl) has three possible minimum energy orientations: trans–gauche (tg), gauche–trans (gt), and gauche–gauche (gg) (Shefter and Trueblood 1965). The cellulose models for each of the three orientations (tg, gt, and gg) were compared with the X-ray data to obtain the best fit, described as the model with the lowest reliability factor (R). For cellulose I, the tg direction gave the best fit with R values of 0.242, 0.292, and 0.349 for the tg, gt, and gg models, respectively (Gardner and Blackwell 1974). The majority view in the literature is that cellulose I has a tg conformation throughout the chain (Sarko et al. 1976; Nishiyama et al. 2008).

The interaction between cellulose microfibrils and solvents or the solubilization process of cellulose microfibrils is a research hotspot (Zhang et al. 2019; Zhou et al. 2021). Still, studies have yet to investigate cellulose’s behavior and conformational changes in non-solvent liquids. It is crucial to observe the effect of polarity on the internal and external structure of crystals in polar or nonpolar liquid media and to explore the mechanism behind it. Related works have provided a substantial basis for this study, including cellulose’s room- and high-temperature behavior in a single solvent (Matthews et al. 2006, 2011, 2012; Zhang et al. 2011) and cellulose phase transitions (Bellesia et al. 2011). These works have demonstrated a helpful tool – computer simulation. Molecular dynamics (MD) computational methods have gained much attention in cellulose structure studies because of their ability to obtain valuable information regarding cellulose structure and properties to have sufficient sampling to explore cellulose conformation. MD simulations were used to model the wetting of cellulose surfaces by water. The theoretical wetting limit of cellulose was solved by considering the surfaces (110) and (100) of the Iβ variant, respectively (Mazeau and Rivet 2008). Gross and Chu (2010) used MD simulations to observe cellulose microfibrils in water and concluded that inter-sheet interactions of cellulose microfibrils were the most robust component. The following year, the team used the same method to study the state of cellulose in two liquids (water and BmimCl), showing that the insolubility of cellulose in water results mainly from a reduction in liquid entropy (Gross et al. 2011). In addition, a considerable number of MD simulations focused on the distortion phenomenon of cellulose microfibrils. Matthews et al. (2006) first reported the tendency of right-handed twisting of cellulose Iβ microfibrils under Charmm force field in a short time (< 200 ps). In the same year, Yui et al. (2006) reported similar behavior of different microfibril models under Glycam force fields. Hadden et al. (2013) investigated the driving forces behind the twisting behavior, including the role of the liquid, the effect of non-bonding force field parameters, and the use of explicitly modeled oxygen lone pairs in the solute and liquid. The results showed that the twisting of microfibrils is influenced by van der Waals interactions and is counteracted by intra-chain hydrogen bonding at the microfibril surface and liquid effects.

MD calculations strongly depend on molecular force fields. The appropriate force fields can match experimental phenomena well or even be used to accurately calculate the system’s internal state to predict the experimental phenomena. MD simulations of cellulose are usually based on the molecular force fields of carbohydrates, and the most commonly used molecular force fields for cellulose are CHARMM36 and GLYCAM06 force fields (C36 and G06) (MacKerell et al. 1998; Kirschner et al. 2008; Huang et al. 2017).

In this work, two common carbohydrate force fields (C36 and G06) were used to study the behavior of cellulose Iβ microfiber in a large number of polar or nonpolar liquids, focusing on the crystal parameters, conformational changes, and fiber twisting of cellulose Iβ microfiber under the influence of liquid polarity. The effect of liquid polarity on model conformation was analyzed from an energy perspective. The reason why these liquid media cannot dissolve cellulose may be that the liquid media cannot penetrate into the crystal, cannot change the internal hydroxymethyl structure, and cannot destroy the hydrogen bond network and van der Waals forces in the crystal structure. It was noteworthy that a short time might exist to not easily discern the simulation results in liquids of similar polarity. Therefore, G06 with an electrostatic scale factor of 1,4 was chosen for the simulation (G06 is fully compatible with the Amber14 force field with a default factor of 5/6). The purpose was to “speed up” the simulation process and to observe more simulation results. The “C36 model” and “G06 model” in this paper refer to the cellulose Iβ microfiber model (named Iβ-MF) in the C36 and G06 force fields, respectively.

EXPERIMENTAL

Cellulose Iβ Microfiber Modeling

A crystal structure file of cellulose Iβ was constructed through the Cellulose-Builder website (Gomes and Skaf 2012), and the model contained 36 glucan chains with 15 cellobiose per chain (Fig. 1). The structure was a monoclinic P2₁ space group with lattice parameters of a=7.784 Å, b=8.201 Å, c=10.380 Å, and y=96.5°. Current studies generally agreed that the cross-sectional shape of cellulose microfibrils was hexagonal. The number of cellulose microfibrils containing cellulose chains was controversial, and both 18-24 and 36 chains models had been reported (Cosgrove 2014; Ding et al. 2014). The 36-chain Iβ-MF was chosen because a more significant number of chains would allow more structural information and data to be obtained for better analysis of the behavior of Iβ-MF in liquids.

Liquid Polarity Calculation

The inability to disrupt the crystal structure was a prerequisite for studying the behavior of Iβ-MF in liquid media, so none of the selected liquids could dissolve cellulose. The preferred liquids were water, N,N-dimethylformamide (DMF), formic acid, acetic acid, propionic acid, methanol, ethanol, propanol, pyridine, carbon tetrachloride (CCl₄), trichloromethane (HCCl₃), dichloromethane (H₂CCl₂), benzene, carbon disulfide, n-hexane, and cyclohexane. The molecular polarity index (MPI) was an important parameter to express the polarity of a molecule and can measure the overall polarity of the molecule. MPI was defined as shown in Eq. 1 (Liu et al. 2021),

(1)

where V is the molecular electrostatic potential, A is the molecular surface area, and the function is integrated over the molecular surface S. The system’s charge distribution determines the polarity’s magnitude. Uneven charge distribution led to differences in the molecular surface’s electrostatic potential, causing a polarity change. The more inhomogeneous the charge distribution, the greater the polarity and MPI.

The computational software optimized all molecular structures at the B3LYP-D3(BJ)/def2-SVP level. The single point energy was calculated at B3LYP-D3(BJ)/def2-TZVP, followed by the quantitative analysis of the molecular surface module of Multiwfn (Lu and Chen 2012) to obtain the MPI.

Fig. 1. Iβ-MF and color scheme of hydroxymethyl group conformation. a) Iβ-MF cross-section (ab surface) and chain numbering scheme. The layers were numbered from top to bottom, left to right, e.g., 11, 12, and 13 for the first layer, 21, 22, 23, and 24 for the second layer, and so on. The blue line was a dividing line between the inner and outer cellulose chains. b) The ac surface of cellulose Iβ microfibrils, each chain contained 15 cellobiose. c) The three monomers in the single chain were shown from two orientations, TG in yellow, GG in purple, and GT in green.

Details of MD

The restrained electrostatic potential (RESP) charge (Bayly et al. 1993) was first calculated for all molecular monomers. Their RESP charges for liquid molecules could be calculated from the file in the previous step. For cellulose chains, the excessive number of atoms caused difficulties in structural optimization. To obtain the cellulose chain’s RESP charge, an oligomer was constructed containing five glucose units. The structure was geometrically optimized at the B3LYP-D3(BJ)/def2-SVP level. The single-point task was computed at the B3LYP-D3(BJ)/def2-TZVP level, and the resulting wave function file was imported into the Multiwfn program to compute the RESP charge. The RESP charge was calculated three times to constrain the total charge of the head, central, and tail units to 0, respectively. Finally, the RESP charge of the oligomer was extended to the cellulose chain.

MD simulations were performed using GROMACS 2020.6 (Kutzner et al. 2019) software and two force fields of carbohydrates (G06 and C36). Iβ-MF was placed in the middle of a 60 × 80 × 200 box and wrapped with approximately 16 molecules in each direction. When describing Iβ-MF using G06, all liquid molecules used the Gaff force field to construct molecular topology files. When Iβ-MF were defined using C36, all liquid molecules were required to generate the molecular topology files using the Glycan Reader and Modeler in CHARMM-GUI (Lee et al. 2016). Periodic boundary conditions (PBC) were applied to the three directions of the box to ensure that the number of particles inside the box was constant. The Langevin thermostat and the Nose-Hoover Langevin pressure regulator stabilized the temperature and pressure at 298.15 K and 1 atm, respectively. The TIP3P water model was used with a 2-fs time step for the dynamics. The Particle-Mesh Ewald (PME) method was used to handle the long-range electrostatic forces, and a non-bonding cutoff distance of 10 Å was applied. In the C36 and G06 simulations, the steepest descent scheme was used to minimize the system energy, requiring a maximum force of less than 100 kJ mol^-1 nm^-1, followed by gradual heating to the prescribed temperature. After reaching the desired temperature, the energy-minimized system did a restricted MD on the cellulose microfibers by the NVT ensemble to ensure that the cellulose microfibers would not move before liquid relaxation. Finally, a long-time MD simulation with a duration of 50 ns was performed. The trajectory data of the system were analysed using the Visual Molecular Dynamics (VMD) package (Humphrey et al. 1996).

Interaction Energies

The interaction energy is expected to be a powerful tool for studying weak interactions between molecules, as shown in Eq. 2,

(2)

where E_elec was the electrostatic interaction energy, acting as an attraction (negative) or repulsion (positive); E_vdw was the dispersive interaction energy, corresponding to the long-range Coulomb correlation of electrons, acting as an attractive force. In medium-strength hydrogen bonds and dihydrogen bonds, electrostatic interactions were dominant, complemented by dispersive interactions. Both van der Waals (vdW) interactions and π-π stacking interactions are dispersive.

RESULTS AND DISCUSSION

Liquid Polarity Data

The polarity data for all liquids are shown in Table 1, sorted by MPI from highest to lowest. Water, formic acid, and DMF molecules (MPI>15 kcal mol-1) are strongly polar, acetic acid, methanol, propionic acid, ethanol, pyridine, propanol, and H₂CCl₂ molecules (8<MPI<15 kcal mol^-1) are moderately polar, benzene and HCCl₃ molecules (4<MPI<8 kcal mol^-1) are weakly polar, and CCl₄, CS₂, n-hexane, and cyclohexane molecules (MPI<4 kcal mol^-1) could be considered as non-polar. The pA_polar content, V_min, and V_max could be used as auxiliary indicators for the polarity of liquid molecules. Liquid molecules with a high percentage of pA_polar or large absolute values of V_min and V_max could be considered a prerequisite for high polarity.

Table 1. Polarity Data for 16 Liquid Molecules

^a) The pA_polar was the molecular polar surface area percentage over the total surface area.

^b)The minimum and maximum molecular surface electrostatic potentials were defined as V_min and V_max.

Cellulose Iβ Microfiber Crystal Structure

The RMSD of Iβ-MF (Fig. S1; see Appendix) showed that all systems reached equilibrium within 50 ns. In C36, Iβ-MF showed the most minor structural change in the formic acid liquid, considered an anomaly and discussed in detail in the subsequent subsection.

Table 2. Crystal Parameters of Crystalline Cellulose Iβ

^a) Experimental data at room temperature (Nishiyama et al. 2003)

^b) Lattice parameter values from the simulation were averaged over the final 10 ns. Standard deviations were shown in parentheses.

^c) Lattice parameter values of the center chains

^d) Lattice parameter values of the origin chains

Table 2 shows the simulated lattice parameters of Iβ-MF and the lattice parameters of Iβ determined by synchrotron radiation and neutron diffraction experiments. The lattice parameters of the Iβ-MF were generally similar in the G06 and C36. Combined with Table 2 and Fig. S1, the crystal structure of Iβ-MF did not change substantially in the liquid and could maintain the hexagonal cross-section well. Detailed data on the lattice parameters of the Iβ-MF for all liquid systems are listed in Table S1. For the G06 model, the c-lattice parameter (cellobiose length) was 0.31 Å larger than the experimental value, which could be a force field defect. Table S1 exhibited the pattern of a-lattice parameters with liquid polarity, and the Iβ-MF interlayer distance was positively correlated with polarity. There was a tendency for the Iβ-MF to swell in liquids higher than benzene polarity, while the opposite was true in liquids lower than benzene polarity. In addition, the perturbation of the α- and β-angles of the center chain by the liquid was higher than that of the origin chain. For the C36 model, the standard deviation indicated that the lattice parameters of the model were largely unaffected by the liquid.

The Behavior of Cellulose Microfibers in Liquids

Herein, 36-chain Iβ-MF (DP=30) behavior was reported in different liquids in C36 and G06. Figure S2 organizes snapshots of the cellulose microfibrils for each of the 100 layers at 50 ns for two force fields. Figure S3 depicts the variation of the number of gg, gt, and tg conformations over time for the two force fields, and the conformational data (number of gg, gt, and tg conformations) for the last moments are collated in Fig. 2.

Fig. 2. Average hydroxymethyl conformation of a) CHARMM36 and b) GLYCAM06 simulations at 50 ns. Chart c was an additional illustration of the data with large deviations in the three replicate simulations in Chart b, correlated with the twist data in Fig. 7 (explained in Fig. 7). The data sets H₂CCl₂_3, and CS₂_1 in Chart c were not counted in chart b because of the large difference between them and the remaining two data sets in the same group.

The initial conformation of hydroxymethyl groups in all systems was tg. Several relevant observations could be derived from the simulation results for each force field. C36 simulations showed that at most 30% of the tg conformations were transformed in 50 ns. More tg conformations were retained as the liquid polarity decreased, with a negative correlation between liquid polarity and the number of tg conformations (except for the formic acid liquid system). In liquids other than formic acid, the difference in the number of gg conformations in each system was slight, ranging from 5.6% to 9.4%. The liquid media can only change the hydroxymethyl conformation on the surface and fails to change the internal conformation. It is therefore believed that these liquid media cannot dissolve cellulose because they cannot break the hydrogen bonding network and van der Waals forces of cellulose. The formic acid system is abnormal in each of the three replicate calculations. The simulation results show that the number of tg conformations in the formic acid system fluctuated very little, and there was no conversion of tg conformation to other conformations. The version number of the computer operating system and simulation software was changed to overcome this anomaly, and five calculations were carried out using the simulation parameters; the results are shown in Table S2. Most subsequent results were consistent with the earlier results, with only one set of exceptions, with 157 tg conformations transformed. According to the conclusion above, formic acid and DMF were both highly polar molecules, and the simulation results for both liquids should be similar, but they were very different. C36 might not be well described for both formic acid and cellulose microfibrils.

The anomaly of the formic acid system was not reflected in the G06 simulations. The transformation of cellulose microfibrils is evident in G06, which is in line with an expectation that the narrowing of 1,4 electrostatic scaling factors will necessarily affect the structure of cellulose substantially. The stock of tg conformation in the low-polarity liquid was less than 40%. It could be assumed that the tg conformation of cellulose microfibrils in the low polarity liquid could be transformed entirely given enough time.

The hydroxymethyl group exposed to the liquid was susceptible to conformational transitions, including odd glucose units at the upper edge, even those at the lower edge, and those at the non-reducing ends, which could be observed in the simulation results for C36. In contrast, interlayer conformational differentiation occurred in the G06 model (alternating layer pattern), with the center layers (even layers) mostly in the gg conformation and the origin layers (odd layers) mostly in the gt conformation, a phenomenon that was particularly evident in liquid systems with high polarity. In the center layer, the hydroxymethyl group exhibiting the gg conformation was essentially perpendicular to the mean plane of the pyranose ring and therefore pointed upward and downward to the original chain in the upper and lower layers.

In this alternating layer pattern, the spatial conflict between these center chain hydroxymethyl groups and the upper and lower layers (origin layers) forced the cellulose chains to tilt significantly for the plane of their layer. The center chain (gg conformation) was tilted counterclockwise, and the origin chain (gt conformation) was tilted clockwise when viewed from the non-reducing end, as shown in Fig. 3. The tilting of the chains led to changes in the hydrogen bonding network, and Figs. 4a, b, and c showed the hydrogen bonding network before and after the simulation, respectively. The extra-ring groups could form effective O6-O2 hydrogen bonds between chains in the same layer (Figs. 4b and c) or O6_centerH6_center-O2_origin hydrogen bonds between adjacent layers. This simulated structure, in terms of chain tilt and hydroxymethyl conformation, which was very similar to the high-temperature phase structure of I-HT (Matthews et al. 2012).

Fig. 3. Schematic diagram of the alternating layer pattern of cellulose microfibers

Fig. 4. Schematic diagram of the hydrogen bonding network before and after the simulation. a) interlayer hydrogen bonding network between the center and origin layers, mainly O6_centerH6_center-O2_origin hydrogen bonds; b) interchain hydrogen bonding network diagram inside the origin layer with the center layer in the background; c) interchain hydrogen bonding network diagram inside the center layer with the origin layer in the background

Interactions Between Iβ-MF and Liquids

The interaction between Iβ-MF and liquids was described by the interaction energy (Lu and Chen 2023), as shown in Fig. 5. The E_elec between liquids and Iβ-MF showed a positive correlation with the polarity of liquids. Due to the lack of hydrogen bond donors for DMF and pyridine liquids, the hydrogen bond formation ability was poor compared to short-chain alcohols and short-chain carboxylic acids, resulting in a weaker electrostatic interaction (attraction). The weakening of this interaction was not found to be significantly different in the simulation results for C36 for a short time (Fig. 2a). However, in the G06 model (Fig. 2b), the retention of the tg conformation confirmed that electrostatic interactions between the liquid and the Iβ-MF dominated the Iβ-MF structural changes. Dispersive interactions only played an auxiliary role, and electrostatic interactions were almost non-existent in low polarity liquids. The conformational change of Iβ-MF was guided by dispersive interactions, and about 10% of tg conformations were transformed in the C36 model. In fact, the relationship between liquid polarity and conformational change is energy-driven and dominated by electrostatic interactions between the liquid and Iβ-MF. Moreover, the E_total of the G06 simulation was generally larger than that of the C35 simulation, which was caused by the large structural changes of the G06 model, and the E_total difference between C35 and G06 was smaller at the beginning of the simulation.

Fig. 5. Interaction energy between liquid molecules and Iβ-MF in a) CHARMM36 and b) GLYCAM06

Large changes in atomic charge would inevitably produce large biases, and misleading simulation results would exist (Bayly et al. 1993; Zhang et al. 2023). However, small perturbations in a large system might not produce detectable biases. If the atomic charge of the molecule was near 0, then small perturbations might produce simulation results with large deviations. In the HCCl₃ system, the RESP2(0.1) and RESP2(0.9) charges were calculated, when the average charges of the Cl atoms were -0.0095 and 0.004, respectively, while the RESP charges of the C atom were -0.2323 and -0.2713, respectively.

For the simulation of HCCl₃ systems with different charges, the E_elec of the C36 simulation differed by 463 kJ mol^-1 with no significant change in the tg conformation. The difference in E_elec for G06 was as high as 4040 kJ mol^-1 and the difference in tg conformation was 30%. Atomic charge variations lead to energy changes, which result in different interactions between molecules. The C36 was less sensitive to atomic charges and could withstand some small perturbations. The G06 reduced the 1,4 electrostatic scaling factor, and small changes in the values near the atomic charge critical point produced large deviations. In addition, simulations were performed for H₂CCl₂, CCl₄, and CS₂ under the same conditions with deviation values of 50 kJ mol^-1, 332 kJ mol^-1, and 71 kJ mol^-1, respectively.

Twist of Cellulose Microfiber

The hydrated cellulose microfibrils in both force fields were twisted at 1 ns moment, and the twisting direction is the same. The non-reduced end (distal end, ball-and-stick model) was kept in its original position. The chain orientation is indicated in Fig. 6 with blue and white spiral lines from the non-reduced end, ending at the reducing end. This finding was very similar to the work of Paavilainen et al. (2011). The distortion of the G06 model was significantly higher than that of the C36 model at 1 ns.

Fig. 6. Schematic diagram of the twisting of cellulose microfibrils at the 1ns. The non-reduced end was represented as a ball-and-stick model and placed at the distal end, and the blue and white spiral lines indicated the course of the cellulose chain, with the end being the reduced end.

Inspired by Matthews et al. (2006), the middle 5 cellobiose units were chosen to calculate the overall twist of this short oligosaccharide segment. The dihedral angle C1₁₁-C1₁₁-C1₂₁-C1₂₁ (indicated by the black line in Fig. S4) was used to define the twist, and the average twist values of the last 5 ns are summarized in Fig. 6. The difference in a twist between the center and origin layers of the model under the two force fields in a single simulation was slight and would not be discussed separately. The C36 model had little difference between the overall twist in different systems, except for the formic acid liquid system. In contrast, the overall twist in the G06 model was not related to the liquid polarity. There were data with a large degree of deviation in three replicate simulations (solid black squares in Fig. 7). The results suggest that the Iβ-MF would go through twisting and untwisting stages, and the C36 model was basically in the twisting stage.

A part of the G06 model had already been untwisted. For example, in Fig. S5a, the G06 model underwent rapid untwisting at about 20 ns in the presence of water molecules, and the tg conformation was fully transformed at the same moment (Fig. S3, Inside the black box), and it was clear that the conformational transition influences the untwisting process. Combining Fig. 7 and Fig. 2c revealed that the less tg conformation contained in the same system, the lower the twisting degree (except for the water system), which undoubtedly indicated that the degree of transformation of tg conformation was one of the essential factors affecting the twist degree. In the water liquid system, three repeated calculations showed the complete transformation of the tg conformation. Still, the overall twist of one of them was 4.79, which did not reach the state of untwisting because no alternate layer structure was formed in the 8^th layer of this calculation. According to the trend, this layer should be mostly in gg conformation (e.g., Fig. 3), but in fact, it was mostly in gt conformation (Fig. S6, inside the red box), resulting in a slight increase in the overall twist. Therefore, the untwisting of cellulose microfibers requires two conditions: 1) almost complete transformation of the tg conformation and 2) the formation of the alternate layer structure.

Fig. 7. Average twist values for the last 5 ns of a) CHARMM36 and b) GLYCAM06 simulations. The anomalies in the b-plot correspond precisely to the points with large deviations in Fig. 3c.

DISCUSSION

Because the atomic charge is one of the important factors affecting electrostatic interactions, the effect of the atomic charge at the 0 points on the simulation results was explored. In the HCCl₃ system, changing the Cl atomic charge (-0.0095 and 0.004) resulted in a difference of 463 kJ mol^-1 in E_elec between the two simulations (C36), with no significant change in the tg conformation. However, the difference in E_elec between the two simulations of G06 was up to 4040 kJ mol^-1 and the tg conformation differed by 30%. The change in atomic charge led to a change in energy and thus affected the intermolecular interactions. A charge change of about 0.05 in the G06 simulation could produce large deviations, although such significant deviations did not occur in other systems. Therefore, it was necessary to focus on the charge at the location of the critical point during the charge calculation to avoid obtaining results with large deviations.

It was clear that the reduction of the 1,4 electrostatic scaling factors weakened the torsional potential of the cellulose structure and made the whole system more chaotic. By linking this “acceleration” effect with heating, it was easy to understand why an alternating layer structure similar to the high-temperature phase intermediate I-HT appeared. This alternate layer structure had two characteristics: 1) the center layer was occupied by the gg conformation and the origin layer was occupied by the gt conformation; 2) the cellulose chains of the center and origin layers were inclined counterclockwise and clockwise, respectively. The hydroxymethyl group in the gg conformation was perpendicular to the mean plane of the pyran ring and acted like an “arm” to capture the cellulose chain of the origin layer through the O6_centerH6_center-O2_origin hydrogen bond, forcing the cellulose chain of the origin layer to tilt clockwise. When the origin layer was tilted to a certain extent, the origin layer’s hydroxymethyl group would have a spatial conflict with the center layer, resulting in a slight counterclockwise tilt of the center layer. The formation of the alternating layer structure appears to untwist. The four groups of systems (water, formic acid, hexane, cyclohexane) were re-simulated using a 1,4 electrostatic scaling factor of 1.0, and the simulation results are summarized in Table 3. The G06 simulation results were very similar to those with C36 at the same 1,4 electrostatic scaling factor and the formic acid system could be explained. Changing the electrostatic scaling factor of GLYCAM06 could affect the thermodynamics and kinetics, both of which have a significant impact on the study of cellulose structures with MD simulations.

Table 3. G06 Simulation Data with a 1,4 Electrostatic Scale Factor of 1

In experiments to identify cellulose I by CP/MAS ¹³C NMR spectroscopy, the spectra of cotton staple immersed in methanol in the C-4 region showed the largest deviation from the spectra of cotton staple immersed in water in the spectral region typical of the unordered form of cellulose (80 to 86 ppm). The typical spectral region for the ordered form of cellulose (86 to 92 ppm) remained relatively intact. Crystalline cellulose was not substantially altered by solvent exchange, i.e., the crystalline form was an “in-core” or “in-fiber” form that was inaccessible to the surrounding solvent (Wickholm et al. 1998).

Some experimental evidence has suggested that cellulose microfibrils are twisted. Narrow bands of twisted fibers were observed during the biosynthesis of bacterial cellulose (Benziman et al. 1980). Hanley and Gray reported the right-handed twisting of Micrasterias denticulata microfibrils using atomic force and transmission electron microscopy (Hanley et al. 1997). The authors suggested that the twisting may be an artificial product of sample preparation and that when microfibers are suspended in an aqueous solution before drying onto the substrate for study, they will exhibit a smooth, uniform twist along their length. The untwisting process of cellulose microfibril was more exciting, and although this simulation result was not necessarily correct, it should help to strengthen the general understanding of the structural behavior of I-HT.

Finally, without further experimental data, it was impossible to state which force field simulation results were “correct definitively”. Microfibrils’ size and surface morphology might significantly impact the observable results obtained through molecular dynamics simulations, so further experimental work will be needed to understand the morphology of small-diameter cellulose microfibrils to get a better comparative reference and to be of great value for more accurate cellulose simulation. The results of this study should provide insight into the behavior of cellulose microfibers in different polar liquids and suggest how MD simulations can best be applied to cellulose studies.

CONCLUSIONS

Two force fields, G06 and C36, were used to study the behavior of cellulose Iβ microfibrils in non-solvent liquids, including the conformational transition and driving forces behind the cellulose microfibrils, the interaction with each liquid, and the fiber twisting. Both force field models were twisted in the equilibrium phase at the beginning of the simulation, and the twist direction was consistent. With a twist, the hydroxymethyl conformation inside the cellulose microfibril was transformed (initial conformation: tg). The C36 model retained more than 70% of the tg conformation in 16 liquids, and the retention of the tg conformation increased with decreasing liquid polarity. In fact, the ability of liquid polarity to affect the hydroxymethyl conformation of cellulose was only an appearance, and the real driving force behind it was the electrostatic interaction between liquid molecules and cellulose. In general, the electrostatic interaction was closely related to the hydrogen bonding ability of the molecule and had some connection with the molecular polarity. Since DMF and pyridine molecules lack hydrogen bond donors, the electrostatic interaction energy (absolute value) was small compared to similarly polar liquid molecules, but the simulation data of C36 did not show conclusive evidence in a short time. Therefore, another set of simulations was designed using the G06 force field with a 5/6 1,4 electrostatic scaling factor, with the aim of reducing the barriers to tg conformational rotation and obtaining an “acceleration” effect. In this case, the tg conformation of the model in highly polar liquids (water, formic acid, methanol, etc.) underwent a complete transformation, while about 60% of the tg conformation of the low-polar liquids was transformed. The DMF and pyridine systems were amplified by electrostatic interactions and showed off-trend results, which further verified that the driving force of the conformational transformation of cellulose molecules is the electrostatic interaction between liquid molecules and cellulose.

ACKNOWLEDGMENTS

This work was supported by the National Key Research and Development Program (2021YFE0104500), the National Natural Science Foundation of China (22078114), Key Research and Development Program of Guangzhou Science and Technology Program (202103000011), and the Natural Science Foundation of Guangdong Province (2021A1515010360).

REFERENCES CITED

Bayly, C. I., Cieplak, P., Cornell, W., and Kollman, P. A. (1993). “A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model,” The Journal of Physical Chemistry 97(40), 10269-10280. DOI: 10.1021/j100142a004

Bellesia, G., Chundawat, S. P. S., Langan, P., Dale, B. E., and Gnanakaran, S. (2011). “Probing the early events associated with liquid ammonia pretreatment of native crystalline cellulose,” The Journal of Physical Chemistry B 115(32), 9782-9788. DOI: 10.1021/jp2048844

Benziman, M., Haigler, C. H., Brown, R. M., White, A. R., and Cooper, K. M. (1980). “Cellulose biogenesis: Polymerization and crystallization are coupled processes in Acetobacter xylinum,” Proceedings of the National Academy of Sciences 77(11), 6678-6682. DOI: 10.1073/pnas.77.11.6678

Cosgrove, D. J. (2014). “Re-constructing our models of cellulose and primary cell wall assembly,” Current Opinion in Plant Biology 22, 122-131. DOI: 10.1016/j.pbi.2014.11.001

Ding, S.-Y., Zhao, S., and Zeng, Y. (2014). “Size, shape, and arrangement of native cellulose fibrils in maize cell walls,” Cellulose 21(2), 863-871. DOI: 10.1007/s10570-013-0147-5

Gardner, K. H., and Blackwell, J. (1974). “The structure of native cellulose,” Biopolymers 13(10), 1975-2001. DOI: 10.1002/bip.1974.360131005

Gomes, T. C. F., and Skaf, M. S. (2012). “Cellulose-Builder: A toolkit for building crystalline structures of cellulose,” Journal of Computational Chemistry 33(14), 1338-1346. DOI: 10.1002/jcc.22959

Gross, A. S., Bell, A. T., and Chu, J.-W. (2011). “Thermodynamics of cellulose solvation in water and the ionic liquid 1-butyl-3-methylimidazolim chloride,” The Journal of Physical Chemistry B 115(46), 13433-13440. DOI: 10.1021/jp202415v

Gross, A. S., and Chu, J.-W. (2010). “On the molecular origins of biomass recalcitrance: The interaction network and solvation structures of cellulose microfibrils,” The Journal of Physical Chemistry B 114(42), 13333-13341. DOI: 10.1021/jp106452m

Hadden, J. A., French, A. D., and Woods, R. J. (2013). “Unraveling cellulose microfibrils: A twisted tale,” Biopolymers 99(10), 746-756. DOI: 10.1002/bip.22279

Hanley, S. J., Revol, J.-F., Godbout, L., and Gray, D. G. (1997). “Atomic force microscopy and transmission electron microscopy of cellulose from Micrasterias denticulata; evidence for a chiral helical microfibril twist,” Cellulose 4(3), 209-220. DOI: 10.1023/A:1018483722417

Horii, F., Hirai, A., and Kitamaru, R. (1987). “CP/MAS carbon-13 NMR spectra of the crystalline components of native celluloses,” Macromolecules 20(9), 2117-2120. DOI: 10.1021/ma00175a012

Huang, J., Rauscher, S., Nawrocki, G., Ran, T., Feig, M., de Groot, B. L., Grubmüller, H., and MacKerell, A. D. (2017). “CHARMM36: An improved force field for folded and intrinsically disordered proteins,” Biophysical Journal 112(3, Supplement 1), 175a-176a. DOI: 10.1016/j.bpj.2016.11.971

Humphrey, W., Dalke, A., and Schulten, K. (1996). “VMD: Visual molecular dynamics,” Journal of Molecular Graphics 14(1), 33-38. DOI: /10.1016/0263-7855(96)00018-5

Imai, T., and Sugiyama, J. (1998). “Nanodomains of Iα and Iβ cellulose in algal microfibrils,” Macromolecules 31(18), 6275-6279. DOI: 10.1021/ma980664h

Kamide, K. (2005). “1 – Introduction,” in: Cellulose and Cellulose Derivatives, K. Kamide (ed.), Elsevier, Amsterdam, pp. 1-23.

Kirschner, K. N., Yongye, A. B., Tschampel, S. M., González-Outeiriño, J., Daniels, C. R., Foley, B. L., and Woods, R. J. (2008). “GLYCAM06: A generalizable biomolecular force field. Carbohydrates,” Journal of Computational Chemistry 29(4), 622-655. DOI: /10.1002/jcc.20820

Kutzner, C., Pall, S., Fechner, M., Eszternnann, A., de Groot, B. L., and Grubmuller, H. (2019). “More bang for your buck: Improved use of GPU nodes for GROMACS 2018,” Journal of Computational Chemistry 40(27), 2418-2431. DOI: 10.1002/jcc.26011

Lee, J., Cheng, X., Swails, J. M., Yeom, M. S., Eastman, P. K., Lemkul, J. A., Wei, S., Buckner, J., Jeong, J. C., Qi, Y., Jo, S., Pande, V. S., Case, D. A., Brooks, C. L., III, MacKerell, A. D., Jr., Klauda, J. B., and Im, W. (2016). “CHARMM-GUI Input Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM Simulations Using the CHARMM36 Additive Force Field,” Journal of Chemical Theory and Computation 12(1), 405-413. DOI: 10.1021/acs.jctc.5b00935

Liu, Z., Lu, T., and Chen, Q. (2021). “Intermolecular interaction characteristics of the all-carboatomic ring, cyclo[18]carbon: Focusing on molecular adsorption and stacking,” Carbon 171, 514-523. DOI: /10.1016/j.carbon.2020.09.048

Lu, T., and Chen, F. W. (2012). “Multiwfn: A multifunctional wavefunction analyzer,” Journal of Computational Chemistry 33(5), 580-592. 10.1002/jcc.22885

Lu, T., and Chen, Q. (2023). “Visualization analysis of weak interactions in chemical systems,” in: Reference Module in Chemistry, Molecular Sciences and Chemical Engineering, Elsevier. DOI: 10.1016/B978-0-12-821978-2.00076-3

MacKerell, A. D., Bashford, D., Bellott, M., Dunbrack, R. L., Evanseck, J. D., Field, M. J., Fischer, S., Gao, J., Guo, H., Ha, S., Joseph-McCarthy, D., Kuchnir, L., Kuczera, K., Lau, F. T., Mattos, C., Michnick, S., Ngo, T., Nguyen, D. T., Prodhom, B., Reiher, W. E., Roux, B., Schlenkrich, M., Smith, J. C., Stote, R., Straub, J., Watanabe, M., Wiórkiewicz-Kuczera, J., Yin, D., and Karplus, M. (1998). “All-atom empirical potential for molecular modeling and dynamics studies of proteins,” J Phys Chem B 102(18), 3586-3616. DOI: 10.1021/jp973084f

Matthews, J. F., Beckham, G. T., Bergenstråhle-Wohlert, M., Brady, J. W., Himmel, M. E., and Crowley, M. F. (2012). “Comparison of cellulose Iβ simulations with three carbohydrate force fields,” Journal of Chemical Theory and Computation 8(2), 735-748. DOI: 10.1021/ct2007692

Matthews, J. F., Bergenstråhle, M., Beckham, G. T., Himmel, M. E., Nimlos, M. R., Brady, J. W., and Crowley, M. F. (2011). “High-temperature behavior of cellulose I,” The Journal of Physical Chemistry B 115(10), 2155-2166. DOI: 10.1021/jp1106839

Matthews, J. F., Skopec, C. E., Mason, P. E., Zuccato, P., Torget, R. W., Sugiyama, J., Himmel, M. E., and Brady, J. W. (2006). “Computer simulation studies of microcrystalline cellulose Iβ,” Carbohydrate Research 341(1), 138-152. DOI: 10.1016/j.carres.2005.09.028

Matthews, J. F., Himmel, M. E., and Crowley, M. F. (2012). “Conversion of cellulose Iα to Iβ via a high temperature intermediate (I-HT) and other cellulose phase transformations,” Cellulose 19(1), DOI: 297-306. 10.1007/s10570-011-9608-x

Mazeau, K., and Rivet, A. (2008). “Wetting the (110) and (100) surfaces of Iβ cellulose studied by molecular dynamics,” Biomacromolecules 9(4), 1352-1354. DOI: 10.1021/bm7013872

Nishiyama, Y., Johnson, G. P., French, A. D., Forsyth, V. T., and Langan, P. (2008). “Neutron crystallography, molecular dynamics, and quantum mechanics studies of the nature of hydrogen bonding in cellulose Iβ,” Biomacromolecules 9(11), 3133-3140. DOI: 10.1021/bm800726v

Paavilainen, S., Róg, T., and Vattulainen, I. (2011). “Analysis of twisting of cellulose nanofibrils in atomistic molecular dynamics simulations,” The Journal of Physical Chemistry B 115(14), 3747-3755. DOI: 10.1021/jp111459b

Sarko, A., Southwick, J., and Hayashi, J. (1976). “Packing analysis of carbohydrates and polysaccharides. 7. Crystal structure of cellulose IIII and its relationship to other cellulose polymorphs,” Macromolecules 9(5), 857-863. DOI: 10.1021/ma60053a028

Shefter, E., and Trueblood, K. N. (1965). “The crystal and molecular structure of D(+)-barium uridine5′-phosphate,” Acta Crystallogr 18, 1067-1077. DOI: 10.1107/s0365110x65002530

Wickholm, K., Larsson, P. T., and Iversen, T. (1998). “Assignment of non-crystalline forms in cellulose I by CP/MAS 13C NMR spectroscopy,” Carbohydrate Research 312(3), 123-129. DOI: 10.1016/S0008-6215(98)00236-5

Yui, T., Nishimura, S., Akiba, S., and Hayashi, S. (2006). “Swelling behavior of the cellulose Iβ crystal models by molecular dynamics,” Carbohydrate Research 341(15), 2521-2530. DOI: 10.1016/j.carres.2006.04.051

Zhang, Q., Bulone, V., Ågren, H., and Tu, Y. (2011). “A molecular dynamics study of the thermal response of crystalline cellulose Iβ,” Cellulose 18(2), 207-221. DOI: 10.1007/s10570-010-9491-x

Zhang, Q.-K., Zhang, X.-Q., Wan, J., Yao, N., Song, T.-L., Xie, J., Hou, L.-P., Zhou, M.-Y., Chen, X., Li, B.-Q., Wen, R., Peng, H.-J., Zhang, Q., and Huang, J.-Q. (2023). “Homogeneous and mechanically stable solid–electrolyte interphase enabled by trioxane-modulated electrolytes for lithium metal batteries,” Nature Energy 8(7), 725-735. DOI: 10.1038/s41560-023-01275-y

Zhang, Y., He, H., Liu, Y., Wang, Y., Huo, F., Fan, M., Adidharma, H., Li, X., and Zhang, S. (2019). “Recent progress in theoretical and computational studies on the utilization of lignocellulosic materials,” Green Chemistry 21(1), 9-35. DOI: 10.1039/C8GC02059K

Zhao, Z., Shklyaev, O. E., Nili, A., Mohamed, M. N. A., Kubicki, J. D., Crespi, V. H., and Zhong, L. (2013). “Cellulose microfibril twist, mechanics, and implication for cellulose biosynthesis,” The Journal of Physical Chemistry A 117(12), 2580-2589. DOI: 10.1021/jp3089929

Zhou, S., Jin, K., and Buehler, M. J. (2021). “Understanding plant biomass via computational modeling,” Advanced Materials 33(28), 2003206. DOI: 10.1002/adma.202003206

Article submitted: August 3, 2023; Peer review completed: September 23, 2023; Revisions received: October 7, 2023; Revisions accepted: October 8, 2023; Published: October 19, 2023.

DOI: 10.15376/biores.18.4.8223-8248