A Molecular Dynamics Examination into How Temperature Influences the Dissolution of Cellulose in an Aqueous Cuprammonium Hydroxide Solution

1. Introduction

Cellulose stands as one of the most frequent biological materials on Earth, it is a material employed for thousands of years across diverse realms of human endeavors including the manufacture of paper, Film Production (Cellophane), Packaging, textile, and more [1,2,3]. Additionally, it serves as the primary raw material for the industrial production of liquid fuels using biomass conversion technologies and is a subject of considerable current interest [4]. Cellulose represents a syndiotactic biopolymer comprising anhydroglucose units (AGU) (Figure 1) linked by (1,4) β- glycosidic bonds [5]. It possesses fascinated physicochemical characteristics and appealing properties, including biocompatibility, biodegradability, and non-toxicity. [6,7]. This substance can form at least six distinct crystalline polymorphs and exists in less organized structures, commonly termed amorphous or paracrystalline cellulose [8]. It is categorized as a semi-crystalline polymer, incorporating disordered (amorphous) and well-ordered (crystalline) regions.

Inside the crystalline regions, cellulose chains form robust and complex intra- and inter-molecular hydrogen bonds acting as the cohesion among cellulose molecules, while in the amorphous domains, distribution occurs uniformly along fascinatingthe microfibers. [9,10]. The stability and compactness of the cellulose fiber arise from intra- and interchain hydrogen bonds, forming a network that significantly contributes to cellulose’s insolubility in both organic and non-organic solvents. This network also plays a crucial role in the resistance to natural microbial and enzymatic cellulose degradation [11].

Cellulose causes challenges in terms of dissolution, as it remains insoluble in both water and organic solvents, with limited solubility observed in selected classes of solvents. Available results from the literature show that there are some very important factors influencing the solubility of cellulose, such as the kinetic and thermodynamic control of the dissolution [12].

From a thermodynamic stand point, the dissolution process will occur spontaneously when the free energy change is negative [13]. That is, the process of a polymer dissolving in a solvent, is governed by the free energy of mixing where the entropy component can play a predominant role. This is particularly true in the case of semi-crystalline polymers such as cellulose where the translational and rotational entropy of polymer chains increase by dissolution of the polymer. Thus, cellulose in solution may have a much greater translational and rotational entropy than when it is enclosed in crystal form [4]. As a consequence, the dissolution of cellulose should be favored by increasing the temperature because of the increase of the entropic component of the free energy. However, several studies have shown that lowering the temperature favors the dissolution of cellulose, which suggests that the temperature vs volumetric fraction (ϕ₂) profile of cellulose dissolution correspond to a lower critical solution temperature (LCST) [14]. The mechanism of this unusual behavior is still controversial, but there is strong evidence for the link between conformational changes and temperature. Indeed, segments with conformational flexibility around the C-C bond in O-CH₂CH₂-O can change their conformation depending on temperature. At higher temperatures, they adopt a less polar state around the C-C bond, while at lower temperatures, they move more polar states [15]. As a result, when the temperature decreases and polar states become prevalent and favorable interactions with the polar solvent are promoted by facilitating the dissolution of cellulose [12,16].

Another factor generally involved in the insolubility of cellulose is its amphiphilic property (Figure 2), which arises from the hydrophilic nature of the three hydroxyl groups situated at the equatorial positions of the ring and the hydrophobic nature of the C-H bonds positioned at the axial positions. Therefore, many studies have hypothesized that disturbing the equilibrium between the hydrophilic and hydrophobic nature of cellulose can significantly impact its dissolution [17,18]. Among the works approving such a hypothesis, Bergenstrahle et al [4] recently employed molecular dynamics simulations to calculate mean force potentials for the dissociation of cellulose oligomers in an aqueous solution.

The free energy calculation showed a significant hydrophobic pairing energy, at about 2.0 kcal/mol/residue, which favors the association of cellulose oligomer chains in a manner comparable to that found in the different crystal structures of cellulose.

It is noteworthy to mention that, according to Lindman et al [12] hydrophobic cellulose interactions are often disregarded in favor of the more robust hydrophilic hydrogen bonds, which are considered to be more crucial. From this point of view, the cellulose intra- and inter chain hydrogen bonds are considered to be a key factor for a better understanding of various physico-chemical properties associated with the solubility of cellulose and cellulose-based materials [20,21,22]. Therefore, tracking alterations in these hydrogen bonds may offer a viable approach for studying the cellulose dissolution process.

There are several classifications of cellulose solvents, according to their nature (aqueous or non-aqueous) and their effect on the cellulose structure (derivatizing or non-derivatizing) [12,23]. Derivatizing solvents are systems that modify the cellulose structure with the formation of an intermediate such as an ester, ether, or acetal to facilitate the decomposition of the cellulose and its dissolution. Among the derivatizing solvents, we can mention the DMF/N₂O₄ (N, N- dimethylformamide/dinitrogen tetraoxide) system, which leads to the formation of cellulose trinitrate [20,24], while the mixture of DMSO/paraformaldehyde leads to the formation of methyl cellulose soluble in non-aqueous solvents [20,25]. Non-derivatizing solvants are systems that dissolve cellulose without any chemical modification. As an exemple, The LiCl/DMAc (dimethylacetamide) [26,27], ZnCl2/H2O [28], LiSCN/H2O [26,29], LiBr/H2O, NaOH/H2O [30,31], this class of solvents dissolve the cellulose chains by breaking the intra- and interchain hydrogen bonds without any structural modification.

On the other hand, the literature review showed that the addition of urea as co-solvents or ZnO in the NaOH/H2O system reinforces the dissolution of cellulose by decreasing the hydrophobic interaction between the hydrocarbon part of the AGU [32]. Furthermore, solution of transition metal cations such as Ni (‘Nioxam’) [33], Zn (‘Zincoxen’) [34], Cd (‘Cadoxen’) [35], constitute an important class of solvent systems used in industrial processes for dissolving cellulose. The best known are cuprammonium hydroxide (Cuam) [36,37] and cupriethylene diamine hydroxide (Cuen) [38]. In this type of solvent system, complexes of cellulose AGU units and the transition metal cation is generally observed and the complexed moeity is considered as non-derivatizing intermediates. Thus, Cellulose with a degree of polymerization not exceeding 1000 AGU unit can be dissolved with this type of solvent.

Nowadays many studies have been focused on understanding the mechanism of cellulose’s solubility using molecular dynamics (MD) in various experimental conditions of solvent systems, temperature and pressure. Cai et al [39] investigated the interplay between an individual cellulose chain and a urea-water solvent mixture under different temperature conditions. The MD simulation result revealed a preference for the cellulose chains to engage in hydrogen bond formation with urea molecules more than with water molecules and maintain their stability even at low temperatures. However, when increasing the temperature in the range of 265 to 283 K the interaction between urea and cellulose decreases suggesting that urea molecules are forming an inclusive layer around the cellulose chain, thereby reinforcing the dissolution of cellulose in the urea-containing solvent mixture by minimizing self-interactions among cellulose chains. In the same theme, Bregado et al [40], examined the effect of temperature on dissollution of a 36-chain cellulose Iβ micro-fibril crystal model at 25 MPa in the range of temperature between 298 to 660 K. The results of this MD investigation showed that the cellulose chains dissolved completely at temperatures close to 600 K. The mechanism of the cellulose dissolution starts between 560 K and 580 K where the initial hydration layer separates from the outer chains situated on the hydrophilic planes through hydrogen bond interactions, which promotes the relaxation of the crystal lattice. Both effects potentially contribute to the advantageous dissolution of cellulose in compressed water. Ramakrishnan et al [41] performed MD simulations of cellulose in the mixture of an ionic liquid 1-ethyl-3-methylimidazolium acetate [C2C1Im][OAc] with water at two distinct temperatures, 300K and 433K. The results showed that the RMSD, and the number of cellulose intra and interchain bonds during the dissolution process are primarily influenced by the temperature. Thus, at T=300K the cellulose chains remain intact whereas at 433K a disruption of the chains takes place. The decrease in the number of H-bonds at higher temperature was due to the weakening of the interactions inherent in cellulose H-bonds and the interaction accumulated by IL with cellulose chains at higher temperatures.

The present work aims to report the results obtained in a molecular dynamics study of the effect of temperature on the dissolution of cellulose in Cuam solvent system. Our primary emphasis is on the Cuam solvent system and its mechanism of solubilizing cellulose typically characterized by the formation of cuprammonium complexes with AGU units of cellulose. These complexes lead to the deprotonation of the hydroxyl groups in the coordination bond connected to the C2 and C3 hydroxyl groups of each AGU. [42] (Figure 3) leading to a disruption in the interchain hydrogen bonds network of cellulose and therefore the segregation of the close packed chains of the cristalline structure.

2. Material and Methods

Throughout this study, molecular dynamics (MD) simulations at different temperatures were performed with the Schrödinger suite software (2018.4). The cellulose structure used in these simulations was extracted from an Iβ crystal structure constructed using a cellulose builder [44]. Parameterised by the OPLS-2005 force field [45] This model structure consists of 6 chains with 8 glucose units each. The Cuam solvent system contains water molecules modeled according to the TIP3P model. OH-, SO42-, and Na+ molecules were simulated using the Material package integrated in the Schrödinger 2018.4 software with the same force field. For all simulations, a cubic box of dimensions (8.12 nm × 8.12 nm × 8.12 nm) was used with periodic boundary conditions (PBC). This box contains 15901 atoms. Each simulation was assigned a different temperature, ranging from 250K to 350K. The pressure was controlled with a Berendsen barostat P=1bar, with a time of 10ns and a time interval of 2.0 fs, under a canonical ensemble (NVT) and the simulations were performed using the Desmond package.

2.1. Root Mean Square Deviation RMSD

The root mean square deviation RMSD indicates the mean deviation of a particle from a reference position. It is an important parameter that allows the control of structural changes during simulation processes and also to follow the diffusion of particles in a given solvent [46]. Several authors have studied the average RMSD of cellulose polymer in different solvent systems; as an example Manna and Ghosh [47] showed that the dissolution of cellulose in water-[C2mim][OAc] (1-ethyl-3-methylimidazolium (Emim+)) solvent is associated with high RMSD values.

In this study, the comparison between RMSD of cellulose chains containing 8 units of glucose dissolved in an aqueous Cuam solution at different temperatures ranging from 250k to 350K showed that the average values fluctuate around (1-3.5Ǻ) at T=250K, 270K, and 280K (Figure 4A). The time-dependent variation of the RMSD indicates stability across all systems, with minor fluctuations remaining below 3Å, suggesting the absence of significant dissolution at these temperatures. On the other hand, at higher temperatures, from 320K to 350K (Figure 4B), the cellulose also faces stability of RMSD values which fluctuated around (2-4Ǻ), indicating that no dissolution occurs.

However, between 300K and 310K (Figure 4C), the RMSD curve indicates an initial increase in RMSD vs of simulation’s time. For example, at a temperature of T=300K, the RMSD value increases to (22.5 Å) before stabilizing. Therefore, a temperature of 300K is considered an effective temperature for the dissolution of cellulose in the Cuam system.

According to Manna and Ghosh [47], the high average RMSD values obtained in our study indicate a modification in the conformation of the cellulose crystal structure within the Cuam solvent system leading to the dissolution of cellulose. This observation is confirmed by the visual representations of the final structures in the simulations at 300K, 330K, 350K, 280 and 250K (Figure 5), which show an aggregated form of the cellulose chains at 280K, 250K, 330K and 350K, while a significant separation of the chains is observed at T=300K.

2.2. Radius of Gyration Rg

In addition to the RMSD results, we also studied the radius of gyration (Rg) of cellulose in Cuam solvent system at different temperatures from 250K up to 350K, to observe the effect of temperature in the dissolution processus during simulation. The radius of gyration (Rg) is a measure defined as the root mean square distance of the atoms of a molecule from their common center of mass [48]. The radius of gyration (Rg) provides a general overview of the compactness of a structure by applying the following formula [49]:

$$Rg=\sqrt{\frac{1}{N}}\sum _N^{i=0}{\left(ri-rm\right)}^2$$

The simulation of cellulose chains at low temperatures between 250 to 280K reveals that the Rg values (12-12.6 Ǻ) are very stable, with a fluctuation around 0.6 Ǻ (Figure 6A). Therefore, the stability of the values indicates that the cellulose is still compact, and no dissolution has occurred in the Cuam solvent system at low temperature. The Rg value increases to 30 Ǻ (Figure 6C) at a temperature between 300K and 310K. Due to the dissolution of cellulose chains which induces significant changes of the conformational arrangement of cellulose crystalline structure.

Above 310 K, the Rg values of the system decrease significantly (Figure 6B) and range in the same order than for temperatures below 300 K. It is noted that when the temperature is increased to values higher than T=310K, the cellulose becomes more compact due to cellulose insolubility.

2.3. Interaction Energy

The dissolution of cellulose is the process of breaking the chain bonds of cellulose into individual chains. Indeed, the interaction energy represents the bonding strength between the cellulose chains and the increase in this interaction energy, leads to a close packing of the cellulose chains [24]. Therefore, to get a better understanding of the effect of temperature in the dissolution of cellulose in the Cuam solvent system, we examined the interaction energy between the cellulose chains and solvent at different temperatures along the simulation process (Figure 7).

According to the presented curve (Figure 7A), it can be observed that the interaction energy between the Cuam solvent and cellulose chains remains stable at temperatures of 280K, 270K, and 250K. However, at temperatures ranging from 300K to 310K (Figure 7B), the interaction energy starts to decrease, from a value of -300 Kcal/mol to -350 Kcal/mol, indicating that cellulose chains have a preference for interaction with the solvent within this temperature range. In contrast, when the temperature increases from 320K to 350K (Figure 7C), the cellulose-solvent interaction energy remains stable compared to the T=300K and T=310K. This means that the binding force between cellulose chains is stronger, and therefore, cellulose is not dissolved in the solvent. The cellulose-solvent interaction energy at different temperatures is in agreement with the analysis of Rg cellulose-solvent and RMSD calculations, where cellulose interacts more with the solvent at T=300K than with the solvent at low and high temperatures.

2.4. Hydrogen Bonds

Intermolecular hydrogen bonds, along with numerous C-H...O inter-shells and van der Waals interactions, contribute to the stabilization of cellulose crystal Iβ [12,20]. Thus, during the dissolution process, the initial step involves breaking these hydrogen bonds to initiate the dissolution process. This investigation delves into the variation in the average number of hydrogen bonds focusing on two key intrachain hydrogen bonding patterns and two primary interchain patterns between two cellulose chains (native and derivatized one and between two natives’ ones).

In our study, for temperatures below 300 K (i.e. 250K and 270K) and above 310 K (i.e. 330K and 350K), we observed that the number of interchain hydrogen bonds remained constant (see Figure 8A–C), as the intrachain hydrogen bonds remain approximately stable arround the mean value of 8 H-bonds, whether between two native chains or between a native chain and a derivatized chain (see Figure 8D–F). This stability is attributed to the absence of breakage of H-Bond bonds at these temperatures, which is due to insolubility of cellulose in the Cuam solution.

On the contrary, at T = 300 K, a significant change in the intra-chain and inter-chain hydrogen bonding patterns between a native chain and a derivatized chain during simulation was observed. From (Figure 9A), which represents the variation in the number of interchain bonds throughout the simulation time, we can see a decrease in the number of H bonds over the simulation time. This decrease can be attributed to the breaking of interchain bonds. Analysis of the intrachain hydrogen bonds (Figure 9B) also shows a decrease in the number of H bonds. An increase in the disruption of interchain and intrachain hydrogen bonds leads to greater accessibility for hydroxyl groups to interact with the solvent, thus facilitating the dissolution process.

To better understand the changes in H bonds at T = 300K, we measured the number of interchain hydrogen bonds. By presenting the frequency of each occurrence in histograms (Figure 11A,B) it is possible to get an insight of the behavior of the cellulose crystals. For interchain bonds between two native chains, the histogram shows values ranging from 1 to 16, and between the native and the derivative, the values are lower, ranging from 1 to 12 in (Figure 11B).

Both frequency distributions exhibit a symmetrical shape, demonstrating a normal distribution with distinct peaks. In (Figure 11A), the peak is centered at 9 H-bonds, while in (Figure 11B), it is centered at 6 H-bonds. As per existing literature, the stable conformation of cellulose involves intermolecular hydrogen bonds (C-O---HO) in a parallel arrangement. The dissolution of cellulose results from the solvent’s ability to disrupt these intermolecular hydrogen bonds (C-O---HO) among cellulose molecules and eliminate hydrophobic interactions (C-H--- C). In our study at T=300K, the organization of native/derived cellulose chains underwent significant changes, making them less inclined to align in parallel due to the replacement of hydrogen with Cu(NH₃)₂. As explained in our prior article [36], the new Cuam complex can permeate native cellulose because of increased chain rotation. Consequently, new, albeit weaker, hydrogen bonds are formed, primarily (N-H----O-CH) and (CH₂-CH-O NH), explaining the ease of breaking these bonds and subsequently dissolving cellulose. This work is the continuation of others ones in our group [50,51,52].

3. Conclusions

In this study, molecular dynamics simulation was conducted using cellulose (Iβ) in the Cuam solvent at different temperatures to gain a molecular-level understanding of the role of temperature in the efficient dissolution of cellulose. Cellulose dissolution occurred at temperatures of 300K and 310K, however, insolubility of cellulose was noted when the temperature was increased or decreased beyond the temperature range between 300K and 310K. Comparison of the RMSD and Rg values of the system at different temperatures suggests that deconstruction of cellulose chains occurs at 300K, while at temperatures between 250K and 280K, as well as between 320K and 350K, the values remain stable, indicating no conformational changes in these temperature ranges. The analysis of the interaction energy between the solvent and cellulose chains is also consistent with the RMSD and Rg results. The interaction energy between cellulose and the Cuam solvent decreases at T=300K-310K, indicating that cellulose chains prefer to interact with the solvent. Counting the number of intrachain and interchain hydrogen bonds during the simulation at T=300K-310K revealed a decrease in the number of interchain hydrogen bonds due to complexation of hydroxyl groups (linked to the C2 and C3 carbons of each anhydroglucose) by Cu(NH₃)₂, leading to chain deformation and disruption in the interchain hydrogen bond network of the cellulose crystal.