Prediction and elucidation of cellulose solubility in ionic liquids under high pressure using all-atom molecular dynamics simulations

Abstract

This study investigated the pressure dependence of cellulose solubility in a 60 wt% ionic-liquid mixture of 1-ethyl-3-methylimidazolium acetate and dimethyl sulfoxide ([EMIm][OAc]/DMSO) using all-atom molecular dynamics simulations conducted over a wide pressure range (P = 0.1–1000 MPa) at temperature T = 500 K. The dissolution free energy obtained via the umbrella sampling method (36 windows, 200 ns each) decreased monotonically as pressure increased, indicating that solubility enhanced at elevated pressures. The underlying molecular mechanism was elucidated by performing 100 ns NPT-MD simulations of a 36-chain cellulose system in 60 wt% [EMIm][OAc]/DMSO at each pressure. Analysis of the radial distribution functions, coordination numbers, interaction energies, chain dispersibility, hydrogen bond populations, and conformational transitions revealed a dual mechanism: (1) pressure-induced compression strengthens solute–solvent interactions and weakens cellulose–cellulose contacts, promoting chain dispersion; and (2) conformational rearrangements of the C6 hydroxyl groups disrupt intramolecular hydrogen bonds and favor new intermolecular hydrogen bonding with [OAc]⁻ anions. This molecular-level insight demonstrated that pressure enhances cellulose solubility through both volumetric compression and active conformational mechanisms, thereby providing guidance for the design of high-pressure cellulose dissolution processes.

---

1 Introduction

Cellulose is the primary structural polysaccharide in plant cell walls and accounts for approximately 50% of their dry weight. It is the most abundant biomass resource on Earth and is recognized for its importance as a sustainable carbon source.¹ Owing to its non-toxicity, biocompatibility, and biodegradability, cellulose has a diverse range of potential applications, including wound dressings,^2,3 scaffolds for tissue engineering,⁴ and drug delivery systems.^5–7 In environmental applications, cellulose-based materials have been explored as filtration and separation membranes,⁸ adsorbents for aqueous organic pollutants,^9,10 and materials for heavy metal ion recovery.^11,12 Additionally, they can be used in agriculture as water-retaining and slow-release fertilizer carriers^13–15 and in consumer products as absorbent components in diapers and feminine hygiene products.¹⁶

However, natural cellulose exhibits a high degree of crystallinity (50–95%) and is stabilized by strong intra- and intermolecular hydrogen bonds, making it insoluble in common solvents like water and organic solvents.^17–20 Therefore, solvents such as N-methylmorpholine-N-oxide (NMMO), N,N-dimethylacetamide/lithium chloride (DMAc/LiCl), dimethyl sulfoxide/tetra-n-butylammonium fluoride trihydrate (DMSO/TBAF), choline acetate, and aqueous sodium hydroxide/urea solutions are required. However, many of these solvents are problematic from environmental and practical standpoints because they are often highly toxic, generate volatile organic compounds, are prone to chemical degradation (limiting their reusability), and produce difficult-to-treat liquid waste or secondary pollutants.^21–23

In recent years, ionic liquids have attracted significant attention as green solvents owing to their ability to dissolve cellulose, low volatility, low toxicity,²⁴ high thermal stability, and reusability. Among the ionic liquids that can dissolve cellulose, imidazolium-based carboxylate ionic liquids and imidazolium-based halide (chloride and bromide) ionic liquids are well known. Compared to the latter, imidazolium-based carboxylate ionic liquids exhibit superior cellulose dissolution capability, lower viscosity, and lower melting points, allowing cellulose dissolution at room temperature (≈25 °C).²⁵

Furthermore, Tomimatsu et al. demonstrated that adding DMSO as a cosolvent increased cellulose solubility, with a 60 wt% 1-alkyl-3-methylimidazolium acetate/DMSO mixture yielding the highest solubility.²⁶ In this system, DMSO does not directly contribute to cellulose dissolution, since it lacks inherent cellulose-dissolving capacity, but instead reduces solution viscosity, thereby accelerating the dissolution process. Additionally, it weakens the interactions between ionic liquid cations and anions, increasing the number of free ions, which in turn enhances solvent–cellulose interactions and further increases cellulose solubility.²⁷

Moreover, several studies have reported that thermal treatment increases cellulose solubility by enhancing both the dissolution rate and solubility.^28–30 Heating is a straightforward approach from both theoretical and practical standpoints, as it requires minimal modifications to existing processes or equipment. However, thermal treatment also presents several challenges, including cellulose degradation, which may reduce mechanical strength, and caramelization reactions that cause discoloration. Additionally, heating requires relatively high energy input.

In contrast, applying pressure offers several advantages. Unlike heat, pressure generally does not break covalent bonds, requires less energy, and applies a uniform force. Therefore, understanding the effects of pressure on cellulose solubility is of considerable importance. Chang et al. investigated the solution structure of a cellulose/1-ethyl-3-methylimidazolium acetate ([EMIm][OAc]) system under high pressure using a diamond anvil cell at pressures of up to 3 GPa. Fourier-transform infrared spectroscopy (FT-IR) analysis revealed that pressure enhanced the interactions between cellulose and the ionic liquid, particularly between cellulose and cations.³¹ In general, stronger solute–solvent interactions are known to improve solubility. Based on this study, the application of pressure is expected to increase the solubility of cellulose in ionic liquids.

In dissolution experiments, including those involving cellulose, the dissolution process must be monitored, and the solution must be stirred. However, specialized equipment is typically required to conduct these procedures under high pressures, making such experiments technically challenging. Consequently, experimental studies on the dissolution behavior of cellulose under high-pressure conditions remain extremely limited. Chang et al. provided valuable insights into the interactions between cellulose and ionic liquids in high-pressure environments. However, the specific details of cellulose solubility in ionic liquids under high pressure remain largely unexplored.

Therefore, in this study, we used all-atom molecular dynamics (AA-MD) simulations to predict the solubility of cellulose in a 60 wt% [EMIm][OAc]/DMSO solvent system across a wide pressure range (P = 0.1–1000 MPa), which is difficult to achieve experimentally. Furthermore, we aimed to elucidate the molecular mechanisms by which pressure affects cellulose solubility.

In Section 3.1, we evaluate the validity of the force field parameters used in this study by calculating the densities of the ionic liquid, DMSO, and the ionic liquid/DMSO mixture at various temperatures and pressures through MD simulations and comparing them with experimental data. In Section 3.2, we calculate the dissolution free energy using the umbrella sampling method and use these results to predict the pressure dependence of cellulose solubility. In Section 3.3, we describe MD simulations of the cellulose/[EMIm][OAc]/DMSO system. We analyzed the radial distribution functions, coordination numbers, interaction energies, cellulose chain dispersibility, hydrogen bond populations, and conformational changes to elucidate the mechanism by which pressure affects solubility at the molecular level.

2 Computational methods

AA-MD simulations were conducted using the GROMACS 2023.4.³² The applied force field parameters included the Chemistry at HARvard Macromolecular Mechanics (CHARMM) Force Field^33,34 for cellulose (degree of polymerization, DP = 8), Sambasivarao's parameter set³⁵ for [EMIm][OAc], and AM1-BCC-derived charges³⁶ for DMSO. Regarding the cellulose model, a chain with a DP of 8 was selected due to the prohibitive computational cost of simulating high-molecular-weight polymers. While this model effectively captures the essential local solute–solvent interactions governing dissolution, it represents an idealized system that excludes long-range polymer entanglement and crystalline cooperativity. Consequently, it should be noted that the calculated solubility may be overestimated compared to experimental values for high-DP cellulose. Intermolecular interactions were modeled using the Optimized Potentials for Liquid Simulations-All-Atom (OPLS-AA/L) force field,³⁷ whereas intramolecular interactions were described using the General Amber Force Field (GAFF).³⁸ Temperature regulation was implemented using a velocity-rescaling thermostat,³⁹ and pressure control was managed with a Berendsen barostat⁴⁰ during equilibration and a Parrinello–Rahman barostat^41,42 for production runs. The electrostatic interactions were calculated using the Smooth Particle Mesh Ewald (SPME) method.⁴³ A cutoff radius of 1.3 nm was applied to both the van der Waals interactions and the real-space component of the SPME calculation. A time step of 2 fs was considered, and all bond distances involving hydrogen atoms were constrained using a LINear Constraint Solver (LINCS) algorithm.⁴⁴ Initial packing of the solvent molecules was performed using Packmol.⁴⁵ Molecular and ionic structures were visualized using visual molecular dynamics (VMD).⁴⁶

2.1 Force field evaluation

A cubic simulation box (6 nm × 6 nm × 6 nm) with periodic boundary conditions (PBCs) was constructed, and the molecules and ions were randomly distributed according to the compositions listed in Table 1. The system underwent energy minimization for 50 000 steps, followed by a 10 ns NVT (number of particles, volume, and temperature) equilibration at T = 500 K. Subsequently, a 10 ns NPT (number of particles, pressure, and temperature) equilibration and a 50 ns NPT production run were conducted under the temperature and pressure conditions specified in Table 1. System densities were calculated from the production trajectories, and the validity of the force field used in this study was evaluated based on the relative error with respect to the experimental density values.

Table 1 Calculation conditions for each system [number of pairs or molecules (N), temperature (T), and pressure (P)]

	[EMIm][OAc]	DMSO	60wt% [EMIm][OAc]/DMSO
Number of pairs (or molecules)	400	667	400/667
T (K)	293, 323, 373 (P = 0.1 MPa)	293, 333, 373 (P = 0.1 MPa)	293, 323, 373 (P = 0.1 MPa)
P (MPa)	0.1, 50, 100 (T = 323 K)	0.1, 15, 35 (T = 333 K)	0.1, 50, 100 (T = 323 K)

---

2.2 Umbrella sampling

2.2.1 Preparation of the solvated system. The potential of mean force (PMF), which represents the free energy profile for the dissolution of a single cellulose chain from its crystal surface, was calculated using the umbrella sampling method combined with the Weighted Histogram Analysis Method (WHAM).^47,48

The initial setup of the simulation system is illustrated in Fig. 1. An initial cellulose Iβ crystal, composed of 36 chains each with a DP of 8, was generated using cellulose-builder.⁴⁹ The system used for the simulation was excised from this crystal and consisted of a crystal slab containing 11 cellulose chains (DP = 8) and a single isolated cellulose chain (DP = 8). These components were placed in a 7 nm × 7 nm × 10 nm simulation box with PBCs. The box was then solvated with 1150 [EMIm][OAc] ion pairs and 1670 DMSO molecules to obtain a 60 wt% [EMIm][OAc]/DMSO solvent mixture.

Fig. 1 Simulation setup for umbrella sampling. The figure shows an overlay of all 36 simulation windows, representing the path of the single cellulose chain as it moved away from the crystal surface. For clarity, the solvent molecules are not shown.

The reaction coordinate for the dissolution process is defined as the distance between the center-of-mass (COM) of the crystal slab and the single chain. This coordinate was chosen because it intuitively and directly describes the physical process of dissolution, specifically, the dissociation of the solute molecule from the solid surface into the bulk solvent. Sarmad et al.⁵⁰ used a similar approach in recent studies.⁵⁰ As shown in Fig. 1 and Table S1, this distance was sampled across 36 windows. The windows were spaced at 0.05 nm intervals from 0.6 to 1.55 nm (20 windows) and at 0.2 nm intervals from 1.6 to 4.6 nm (16 windows). In each window, the COM distance was constrained by the harmonic biasing potential. A strong force constant of 10 000 kJ mol⁻¹ nm⁻² was applied for the first 10 windows (0.6–1.05 nm), where strong intermolecular forces were expected. For the remaining 26 windows (1.1–4.6 nm), a weaker force constant of 1000 kJ mol⁻¹ nm⁻² was used.

All molecular dynamics simulations in this study were performed at an elevated temperature (T = 500 K). This relatively high temperature was chosen because the ionic liquid used as the solvent exhibits extremely high viscosity, particularly under high-pressure conditions, where both translational and conformational relaxations are significantly slowed. At ambient temperatures (e.g., T = 298 K), the system displayed sluggish dynamics, and even submicrosecond simulations failed to yield sufficient conformational sampling for reliable statistical analysis. Under such conditions, the free energy profiles obtained from umbrella sampling exhibited poor convergence and large statistical uncertainties. Therefore, simulations were performed at an elevated temperature to accelerate molecular relaxation, enhance conformational sampling, and obtain statistically meaningful results at a feasible computational cost. The validity and possible limitations related to the elevated temperature are discussed in detail in Section 3.2.3.

2.2.2 System equilibration and production run. For each window, the system was energy-minimized for 50 000 steps, followed by a 1 ns equilibration in the NVT ensemble at T = 500 K. Subsequently, a 200 ns NPT simulation was conducted at the same temperature (T = 500 K) and at each target pressure (P = 0.1, 200, 400, 600, 800, and 1000 MPa), where the initial 20 ns was discarded for equilibration, and the remaining 180 ns was used for production sampling. The structural integrity of the crystal slab throughout the simulation was maintained by applying a harmonic potential of 100 kJ mol⁻¹ nm⁻² to restrain the positions of the atoms.

2.2.3 Analysis of the sampling data. Finally, the PMF was constructed from the biased COM distance distributions collected from all the windows using WHAM. The statistical uncertainty was estimated by dividing the last 180 ns NPT of the data (from 20 to 200 ns) into six independent 30 ns blocks. The final PMF profile and its error bars represent the mean and standard deviation of the results for the six blocks, respectively.

For enthalpy analysis, the average system enthalpies were calculated from the equilibrium trajectory for each of the 36 windows. The resulting plot of the enthalpy values as a function of the reaction coordinates constitutes the dissolution enthalpy profile.

2.3 Structural analysis of the cellulose solution (cellulose/[EMIm][OAc]/DMSO system)

The simulations for the structural analysis were performed under the same temperature and pressure conditions (T = 500 K; P = 0.1–1000 MPa) as the umbrella sampling calculations to ensure consistent sampling conditions and to elucidate the underlying mechanisms under comparable thermodynamic states.

2.3.1 System preparation and simulation protocol. The structural properties of cellulose in the solution were investigated by preparing an initial system. The crystalline cellulose (DP = 8, 36 chains, Iβ-type crystal) constructed using cellulose-builder⁴⁹ was placed at the center of a simulation box and solvated with 1500 [EMIm][OAc] ion pairs and 2180 DMSO molecules. This resulted in a system containing 10 wt% cellulose in a 60 wt% [EMIm][OAc]/DMSO solvent mixture.

A multi-step protocol was used to equilibrate the dissolved states. First, the system underwent energy minimization over 50 000 steps, followed by a 10 ns simulation in the NVT ensemble at T = 500 K to relax the solvent molecules, and then a 10 ns simulation in the NPT ensemble was conducted at T = 500 K and P = 0.1 MPa. Subsequently, a 400 ns NPT simulation was performed under the same conditions to ensure complete dissolution, which yielded a fully equilibrated cellulose solution at ambient pressure.

Using this equilibrated structure as a starting point, separate systems were prepared for higher-pressure conditions. Each system was re-equilibrated for 20 ns in the NPT ensemble at T = 500 K and the target pressures of P = 200, 400, 600, 800, or 1000 MPa. Finally, 100 ns NPT production runs were conducted for each pressure condition (including P = 0.1 MPa) to generate trajectories for analysis.

2.3.2 Data analysis. Trajectories from the production runs were used to analyze the structural and energetic properties of the cellulose solution. Subsequently, the interaction energies, coordination numbers, radial distribution functions (RDFs), cellulose chain dispersibility, hydrogen bond populations, and cellulose dihedral angle distributions were analyzed.

The interaction energies, coordination numbers, and RDFs were calculated for cation–anion, cellulose–solvent (anions, cations, and DMSO), and cellulose–cellulose pairs (atom labels are shown in Fig. 2). In the cation–anion interactions, the RDF and coordination numbers were computed between the H2 and C1 atoms; in the cellulose–anion interactions, between the H atom of the OH group (H2, H3, and H6) and C1; in the cellulose–cation interactions, between the O atoms of the OH group (O2, O3, and O6) and H2; in the cellulose–DMSO interactions, between the H atom of the OH (H2, H3, and H6) and the O atoms; and in the cellulose–cellulose interactions, between the carbon atoms (C1, C2, C3, C4, C5, and C6) of the different molecules (ions).

Fig. 2 Crystalline cellulose, amorphous cellulose, cellulose chain, ionic liquid ions, and DMSO structures with their corresponding atom type labels modeled in this study.

The interaction energy (U_nonbonded) between the single cellulose chain and crystalline cellulose at each pressure condition was calculated as the sum of the coulombic potential energy (U_Coulomb) and van der Waals potential energy (U_vdW):

U_nonbonded = U_Coulomb − U_vdW (1)

The cellulose distribution inhomogeneity index (σ_dist) was calculated as the spatial standard deviation of the time-averaged cellulose atom counts over all the mesh cells. That is,

where 〈n_i〉_t is the time-averaged number of cellulose atoms in mesh cell i, 〈n〉_t,mesh is the average of 〈n_i〉_t over all meshes, and N_mesh = 1000 is the total number of mesh cells. A smaller σ_dist indicates a more spatially homogeneous (uniformly dispersed) cellulose distribution.

Hydrogen bonds were identified when the distance between the H atom and the acceptor atom was less than 0.35 nm, and the angle between the donor, hydrogen, and acceptor atoms exceeded 150°. In this study, donor atoms were defined as the hydroxyl O atoms in cellulose (O2, O3, and O6), whereas acceptor atoms included the hydroxyl O atoms (O2, O3, and O6), ring O atoms (O5), and glycosidic O atoms (O4) within the cellulose structure. The atom labels are presented in Fig. 2.

3 Results and discussion

3.1 Evaluation of the force field

The validity of the force field used in this study was evaluated using MD simulations to calculate the densities of [EMIm][OAc], DMSO, and 60 wt% [EMIm][OAc]/DMSO mixtures across a broad range of temperatures and pressures reported in the literature.^51–53 The relative errors between the simulated and experimental values were calculated and are presented in Tables S2–S4. The relative errors ranged from −0.0765% to −0.9332% under all conditions, demonstrating excellent agreement with the experimental data. These results confirmed the accuracy and reliability of the force field used in this study.

3.2 Prediction of the pressure dependence of solubility

Kamlet–Taft (KT) parameters—the hydrogen bond acidity (α), basicity (β), and polarizability (π*)—have been widely used to describe the characteristics of solvents for cellulose dissolution. In particular, the β value, representing the hydrogen bond accepting ability of the solvent, is known to have a strong correlation with the cellulose solubility.^54–61 Similarly, theoretical studies using MD simulations and density functional theory (DFT) calculations have frequently reported that the interaction energy between cellulose and the solvent, especially the strength of the hydrogen bonds, is the dominant factor governing solubility.^62–65 These approaches, which primarily focus on the enthalpic contributions of solute–solvent interactions, are considered effective metrics for understanding the dissolution mechanism, as they show good agreement with experimental findings.

However, these discussions are centered on interaction energies and do not explicitly account for other critical thermodynamic contributions such as changes in entropy and volume that accompany dissolution. In particular, pressure is a key thermodynamic parameter that can dramatically affect the solubility by influencing the volume change (ΔV) and solvation structure of the system, which in turn affects the entropy (ΔS), as described by the fundamental relation dG = VdP − SdT. Therefore, a rigorous evaluation of the pressure dependence of the solubility necessitates a discussion based on the Gibbs free energy, which incorporates both enthalpic and entropic terms. In this study, we addressed this issue by using the umbrella sampling method to calculate the free energy profile along a reaction coordinate defined by pulling a single cellulose chain from its crystal. The change in this profile yielded the free energy of dissolution.

The free energy change calculated by this method, ΔG_dis, corresponds to the standard Gibbs free energy of dissolution: the free energy change for the transfer of cellulose from its solid (crystalline) state to a standard state (e.g., 1 mol L⁻¹) in an infinitely dilute solution. The relationship between this value and the saturated solubility is given by the phase equilibrium condition, where the chemical potential of the solid is equal to that of the solute in a saturated solution:

ΔG_dis(P) = −RT ln[α_sat(P)] = −RT ln[γ_sat(P)X_sat(P)] (3)

Here, α_sat is the activity of the solute (cellulose) in the saturated solution, X_sat is the solubility expressed as a mole fraction, and γ_sat is the activity coefficient at saturation. For an ideal dilute solution where solute–solute interactions are negligible (γ_sat ≈ 1), this equation would allow for the direct quantitative determination of solubility X_sat from ΔG_dis.

Experimentally, the system in this study—cellulose dissolved in the [EMIm][OAc]/DMSO system at concentrations up to approximately 20 wt%—represents a concentrated polymer solution. In such systems, strong steric repulsion and intermolecular interactions among the polymer chains lead to significant nonideality. Consequently, the activity coefficient is expected to deviate substantially from unity (γ_sat ≫ 1), making a quantitative prediction of X_sat from ΔG_dis using eqn (3) extremely challenging.

Although ΔG_dis(P) provides valuable qualitative insight into the effect of the pressure on the overall thermodynamic driving force of dissolution, it should be interpreted as an indicator of the intrinsic pressure dependence of the dissolution process rather than a direct measure of absolute solubility.

3.2.1 Features of the PMF and their physical meanings. The dissolution free energy profile (i.e., PMF) was first analyzed for convergence to ensure the statistical reliability of the calculations. The PMF was calculated over progressively longer sampling time ranges (20–50 ns to 20–200 ns), and the profile converged well with the 20–200 ns dataset (Fig. S1). Furthermore, sufficient overlap between the histograms of adjacent sampling windows confirmed that the statistical sampling was adequate for accurate PMF reconstruction (Fig. S2).

Fig. 3 illustrates the final PMF profiles as a function of the COM distance between the crystal and the single cellulose chain under each pressure condition (an enlarged view of the PMF in the range of 0.2–1.8 nm is shown in Fig. S3). All the profiles exhibited common characteristic features. A local minimum was observed at a COM distance (D_COM) of approximately 0.6 nm, corresponding to the thermodynamically stable state of the single chain bound to the crystal surface. Immediately following this, a local maximum appeared at approximately D_COM = 0.7 nm, representing the activation free energy for the initial separation of the chain from the crystal. As the distance increased from D_COM = 0.7 to 2.5 nm, the PMF decreased until it converged to a stable plateau between D_COM = 2.5 and 4.5 nm. This plateau region indicates that the cellulose chain had fully dissociated from the crystal and was stably solvated by the solvent.

Fig. 3 Pressure dependence of the PMF under P = 0.1, 200, 400, 600, 800, and 1000 MPa.

3.2.2 Activation free energy. Before discussing the overall thermodynamics of dissolution, we analyzed the kinetic barriers. The activation free energy of dissolution (ΔG_act) was calculated using the equation

ΔG_act = G^‡ − G_xtal (4)

where G^‡ is the value of the PMF at the local maximum (D_COM = 0.7 nm) and G_xtal is the value at the local minimum (D_COM = 0.6 nm).

The results show that ΔG_act increased monotonically with the applied pressure (Fig. 4 and Tables S5–S10). According to the Arrhenius-like relationship given by the equation

where k is the dissolution rate constant. The increase in the ΔG_act indicates that the dissolution rate decreases exponentially as pressure increases. Crucially, this trend is in excellent agreement with the pressure dependence of the dissolution rate that we previously reported for a similar solute–solvent system—a 36-chain cellulose crystal (DP = 8) in [EMIm][OAc]/DMSO based on direct AA-MD simulations,⁶⁶ whereby simulations over a pressure range of P = 0.1–1000 MPa also revealed that the dissolution rate decreased exponentially as the pressure increased. This consistency serves as a critical validation of the approach used in this study, demonstrating that the PMF calculated using the umbrella sampling method is statistically robust and provides a reliable foundation for subsequent analysis of thermodynamic solubility.

Fig. 4 Pressure dependence of the activated free energy (ΔG_act) estimated from the PMF under P = 0.1, 200, 400, 600, 800, and 1000 MPa.

3.2.3 Prediction of the pressure dependence on solubility. The overall free energy of dissolution (ΔG_dis) was determined from these profiles using the equation

ΔG_dis = G_solv − G_xtal (6)

where G_solv is defined as the average PMF value in the solvated plateau region (D_COM = 2.5–4.5 nm) and G_xtal is the local minimum PMF value in the bound state at 0.6 nm.

The calculated values reveal a clear and consistent trend (Fig. 5 and Tables S5–S10). As the applied pressure increases, ΔG_dis decreases monotonically. This is the main finding of the present study. A lower ΔG_dis corresponds to a more spontaneous and thermodynamically favorable process, indicating that an elevated pressure acted as a strong driving force promoting the dissolution of cellulose in the [EMIm][OAc]/DMSO solvent.

Fig. 5 Pressure dependence of the dissolution free energy (ΔG_dis), dissolution enthalpy (ΔH_dis), and dissolution entropy term (−TΔS_dis) under P = 0.1, 200, 400, 600, 800, and 1000 MPa. Error bars represent the standard deviation.

However, ΔG_dis itself does not allow for the quantitative determination of solubility, as the solubility X_sat(P) depends on the activity coefficient, which deviates significantly from unity in such non-ideal, concentrated polymer solutions. Consequently, although the observed decrease in ΔG_dis with increasing pressure clearly demonstrates that solubility is thermodynamically enhanced, it should be noted that relying solely on ΔG_dis may overestimate the extent of the actual solubility enhancement due to the large activity coefficients (γ_sat ≫ 1) in concentrated solutions. Therefore, the precise magnitude or functional form of the change in solubility cannot be inferred directly from the present calculations.

Moreover, all simulations in this study were conducted at T = 500 K, as discussed in Section 2.2.1. This elevated temperature was necessary to achieve sufficient conformational sampling of the cellulose and solvent molecules in a highly viscous ionic liquid medium. At lower temperatures, the molecular relaxation and diffusion of the polymer chains are extremely slow; hence, statistically converged free energy profiles cannot be obtained within feasible simulation timescales. Despite this elevated temperature, the present results capture physically reasonable dissolution behavior. For example, at ambient pressure (P = 0.1 MPa), the calculated dissolution free energy (ΔG_dis = −4.38 kJ mol⁻¹) was negative, which is consistent with the experimentally observed solubility of cellulose in [EMIm][OAc]/DMSO mixtures.

Furthermore, numerous studies have shown that, although temperature strongly affects the absolute magnitudes of the thermophysical properties of ionic liquids, such as viscosity and density, the qualitative pressure dependence of these properties remains largely invariant with temperature.^51–53 This indicates that the temperature primarily rescaled the dynamic timescales rather than altering the underlying pressure response. Accordingly, although the simulations were performed at an elevated temperature, the observed trend of decreasing ΔG_dis with increasing pressure is expected to be robust and physically meaningful. Nevertheless, we must acknowledge that the present findings strictly correspond to the thermodynamic conditions at T = 500 K and that quantitative extrapolation of the dissolution free energies to experimental temperatures should be performed with caution. Therefore, the key outcome of this study is the identification of a consistent pressure-induced enhancement of cellulose solubility, rather than the precise numerical value of ΔG_dis under ambient conditions.

3.3 Elucidation of the underlying mechanism

The umbrella sampling results presented in the previous section demonstrate that pressure thermodynamically promotes cellulose dissolution. In this section, we describe the mechanisms underlying this phenomenon. First, we decomposed the dissolution free energy into enthalpic and entropic contributions to identify the thermodynamic origin of the driving force (Section 3.3.1). Next, we discuss the “macroscale” changes in the solvation shells and chain dispersibility through solution structure analysis (Section 3.3.2). Finally, we investigate the specific molecular mechanisms at the “microscale”, including hydrogen bonding networks and conformational transitions, to explain the pressure-enhanced solubility in detail (Section 3.3.3).

3.3.1 Decomposition of ΔG_dis into enthalpic and entropic contributions. The calculated ΔG_dis described in the preceding section was decomposed to elucidate whether the pressure-induced enhancement of solubility is driven by enthalpic (ΔH) or entropic (−TΔS) contributions. The average enthalpy of the solvated state (H_solv) was taken from the D_COM = 2.5–4.5 nm plateau region, whereas the enthalpy of the crystalline state (H_xtal) was taken from the local minimum at approximately D_COM = 0.6 nm (Fig. 3). The enthalpy of dissolution (ΔH_dis) was then calculated using the equation

ΔH_dis = H_solv − H_xtal (7)

The results are shown in Fig. S4. Subsequently, the entropic contribution (−TΔS_dis) was determined by decomposing the free energy according to the equation

ΔG_dis = ΔH_dis − TΔS_dis (8)

Fig. 5 and Tables S11 and S12 illustrate the pressure dependence of these thermodynamic components. The entropic contribution (−TΔS_dis) remained largely independent of the pressure across the entire range in focus. In contrast, the enthalpic contribution (ΔH_dis) decreased monotonically as the pressure increased. This clearly indicated that the pressure-induced enhancement in the cellulose solubility in the [EMIm][OAc]/DMSO system is an enthalpy-driven process.

In summary, the results presented thus far demonstrate that within P = 0.1–1000 MPa, increasing pressure makes the dissolution of cellulose more thermodynamically favorable, primarily due to favorable changes in the enthalpy of the system. This enthalpic origin suggests that the mechanism is rooted in the modulation of intermolecular interactions. Therefore, a detailed analysis of the interaction energies, radial distribution functions, and hydrogen bonding networks between solute–solvent and solute–solute species is necessary to elucidate the molecular-level origins of this pressure effect. These aspects will be discussed in the following sections.

3.3.2 “Macroscale” analysis of the solution structure of the cellulose/[EMIm][OAc]/DMSO system. The solution structure of the dissolved cellulose was investigated at the molecular level by performing MD simulations of the cellulose/[EMIm][OAc]/DMSO system (36 cellulose chains, 1500 pairs of ionic liquids, and 2180 DMSO molecules). Visual inspection of the simulation snapshots after 100 ns confirmed that the cellulose chains remained well-dissolved and did not aggregate at any pressure (P = 0.1–1000 MPa) (Fig. S5).

Subsequently, the interaction energies, coordination numbers, and RDFs were calculated for the cation–anion, cellulose–solvent (anion, cation, and DMSO), and cellulose–cellulose pairs. The interaction energies and coordination numbers are summarized in Fig. 6(A–E), and the RDFs are summarized in Fig. 6(A′–E′). For clarity, the RDFs in Fig. 6(A′–E′) correspond to representative pressure conditions (P = 0.1 and 1000 MPa), as overlaying all pressure-dependent curves would obscure the trends. The complete set of RDFs for all pressures is shown in Fig. S6.

Fig. 6 Interaction energy (U_nonbonded) and coordination number (CN) and radial distribution function [g(r)] for (A and A′) cation–anion, (B and B′) cellulose–anion, (C and C′) cellulose–cation, (D and D′) cellulose–DMSO, and (E and E′) cellulose–cellulose under various pressures. Reference and target atoms for the RDF and coordination number are defined in Section 2.3.2. Error bars indicate the estimated errors.

3.3.2.1 Cation–anion interactions. As pressure increased, the cation–anion interactions became more favorable, strengthening by a factor of 1.10 from P = 0.1 to 1000 MPa [Fig. 6(A) and Table S13]. The anion coordination number around the cations also increased by the same factor [Fig. 6(A) and Table S14]. The RDF analysis revealed that the first and second solvation shell peaks shifted to shorter distances at higher pressures [Fig. 6(A′)]. These results indicate that pressure enhances cation–anion interactions, primarily by decreasing the interionic distance and increasing the coordination number.
3.3.2.2 Cellulose–anion interactions. The cellulose–anion interactions were significantly strengthened under pressure, becoming 1.20 times more favorable at P = 1000 MPa than at P = 0.1 MPa [Fig. 6(B) and Table S15]. This was accompanied by a 1.26-fold increase in the anion coordination number around cellulose [Fig. 6(B) and Table S16]. The RDFs exhibited a corresponding shift in the solvation peaks to shorter distances [Fig. 6(B′)]. This demonstrates that pressure strengthens cellulose–anion interactions through both closer contact and increased coordination.
3.3.2.3 Cellulose–cation interactions. Similarly, the cellulose–cation interactions were enhanced by a factor of 1.24 at P = 1000 MPa, with a 1.14-fold increase in the anion coordination number [Fig. 6(C), Tables S17 and S18]. The RDFs showed a peak shift to shorter distances and the formation of a new shoulder peak at approximately 0.4–0.5 nm at higher pressures [Fig. 6(C′)]. This shoulder suggests that the pressure induced the formation of a new, distinct solvation structure by pushing cations closer to the cellulose chain. This finding provides a basis for a molecular-level explanation of the experimental results reported by Chang et al., who observed pressure-enhanced cellulose–cation interactions using FT-IR spectroscopy.³¹ Our results reveal that this enhancement is driven by a combination of shorter intermolecular distances and the formation of new solvation shells.
3.3.2.4 Cellulose–DMSO interactions. The cellulose–DMSO interactions were also moderately enhanced, with the interaction energy becoming 1.20 times more favorable and the anion coordination number increasing by a factor of 1.08 at P = 1000 MPa [Fig. 6(D), Tables S19 and S20]. The RDFs confirmed that the cellulose–DMSO distance decreased with increasing pressure [Fig. 6(D′)].
3.3.2.5 Cellulose–cellulose interactions: a counter-intuitive finding. Regarding the cellulose–cellulose interactions, the RDFs and coordination numbers [Fig. 6(E′), (E) and Table S21] indicate that the interaction strength at P = 1000 MPa decreases by approximately 2% compared with that at P = 0.1 MPa (Fig. 6(E) and Table S22). This reduction in solute–solute association may facilitate the dispersion of cellulose chains in the solvent, and such enhanced dispersion could contribute to the thermodynamic promotion of cellulose solubility under high-pressure conditions.
3.3.2.6 Spatial distribution homogeneity of cellulose chains in cellulose solution. As previously discussed for the solute–solvent and solvent–solvent pairs in Sub-sections (a)–(d), increasing the pressure consistently enhanced the attractive interactions. This is consistent with the physical intuition that compression reduces the intermolecular distances [the system density of the cellulose solution increased (Fig. S7 and Table S23)]. In contrast, the cellulose–cellulose interactions discussed in Sub-section (e) exhibited a counterintuitive trend, and the interactions weakened even though the interatomic distance showed no significant change, as shown by the RDFs [Fig. 6(E′)].

We also analyzed the spatial distribution homogeneity (dispersibility) of cellulose chains in the solvent system at each pressure level to investigate the origin of this behavior. We defined a distribution inhomogeneity index (σ_dist), as described in Section 2.3.2. Fig. 7 shows a plot of σ_dist as a function of pressure. This shows that σ_dist decreases as the pressure increases, confirming that cellulose becomes more homogeneously dispersed throughout the solution as the pressure increases.

Fig. 7 Pressure dependence of the cellulose distribution inhomogeneity index (σ_dist).

3.3.2.7 Summary and proposed mechanism for pressure-enhanced solubility at “macroscale”. The molecular-level insights provided in Sub-sections (a)–(e) allowed us to propose a mechanism for the effect of pressure on cellulose solubility, as illustrated in Fig. 8. As the pressure increased, the interactions between the cellulose and the solvent, particularly the ionic liquid components, were significantly strengthened. This enhanced solvation of the individual cellulose chains by ions creates more robust solvent shells. These shells act as barriers, preventing cellulose chains from approaching each other and aggregating. Ultimately, this results in a more uniform and dispersed dissolution of cellulose under high-pressure conditions, which reduces unfavorable cellulose–cellulose interactions (aggregation).

Fig. 8 Schematic of the solution structures under relatively low and high pressures.

This molecular mechanism is fully consistent with the thermodynamic result from Section 3.2.3, providing a clear rationale for why the dissolution free energy (ΔG_dis) becomes more favorable under pressure and supporting our central conclusion of pressure-enhanced cellulose solubility.

3.3.3 “Microscale” analysis of the cellulose structure in the [EMIm][OAc]/DMSO system. We analyzed the hydrogen bonds between the cellulose chains and solvent species, together with the dihedral angles of cellulose, to characterize the pressure-induced conformational changes and gain deeper molecular insights into pressure effects.
3.3.3.1 Hydrogen bond network analysis. First, the pressure dependence of the hydrogen bond network was investigated by analyzing the intramolecular (within the cellulose) and intermolecular (between different molecules) hydrogen bonds. Although the number of intermolecular hydrogen bonds showed negligible pressure dependence, the number of intramolecular hydrogen bonds decreased by approximately 12% as the pressure increased to P = 1000 MPa [Fig. 9(A) and Table S24]. A more detailed analysis revealed that the O6^H⋯O3 intramolecular hydrogen bond, which is generally recognized as the most dominant hydrogen-bonding motif within cellulose chains, showed the largest decrease with the increasing pressure [Fig. 10(A) and Table S25].

Fig. 9 (A) Pressure dependence of the number of intra-hydrogen bonds (N_intra-HB) between the cellulose and the number of inter-hydrogen bonds (N_inter-HB) within the cellulose at P = 0.1, 200, 400, 600, 800, and 1000 MPa. (B) Schematic showing the structure of the intra-hydrogen bonds between the cellulose and inter-hydrogen bonds within the cellulose.

Fig. 10 (A) Pressure dependence of the number of intra-hydrogen bonds (N_intra-HB) within cellulose under P = 0.1, 200, 400, 600, 800, and 1000 MPa. (B) Schematic showing the structure of the intra-hydrogen bonds within the cellulose.

In contrast, the number of intermolecular hydrogen bonds between cellulose and the solvent increased with increasing pressure, with the most pronounced increase observed for cellulose–anion interactions. These were not only the most abundant but also exhibited the strongest pressure dependence [Fig. 11(A) and Table S26]. Specifically, within the cellulose–anion hydrogen bonds, the O6^H⋯O_anion hydrogen bond demonstrated the strongest pressure sensitivity [Fig. 11(B) and Table S27].

Fig. 11 Pressure dependence of the number of inter-hydrogen bonds (A) between the cellulose and the solvent (anion, cation, and DMSO) and (B) between the cellulose and the anion under P = 0.1, 200, 400, 600, 800, and 1000 MPa.

In summary, increasing pressure primarily reduced the number of O6^H⋯O3 intramolecular hydrogen bonds while promoting the formation of O6^H⋯O_anion intermolecular hydrogen bonds. This suggests that the C6 hydroxymethyl (CH₂OH) group undergoes a pressure-induced transition from forming intramolecular hydrogen bonds within the cellulose chain to forming intermolecular hydrogen bonds with the anions. This behavior may be attributed either to a pressure-induced conformational change of the CH₂OH group or to enhanced interactions between the cellulose and solvent species due to system compression.

3.3.3.2 Pressure-induced conformational changes. Conformational analysis was conducted to clarify this point. Specifically, the dihedral angles related to the C6 position—φ6 (O5–C5–C6–O6) and τ6 (C5–C6–O6–H6)—were analyzed (Fig. S8). The left panel of Fig. 12 shows a two-dimensional probability difference map, which illustrates how the joint distribution of the two dihedral angles (φ6 and τ6) changed from P = 0.1 to 1000 MPa. The blue regions [(A) and (B)] indicate conformations that become less populated as pressure increased, whereas the red regions [(C) and (D)] indicate conformations that become more populated as pressure increased. The right panel depicts representative conformations associated with each region, where the arrows indicate the orientation vector of the hydroxyl group at the C6 position. As the pressure increases, the C6–OH vector, which initially points “inward” toward the cellulose backbone—an orientation favorable for forming intramolecular hydrogen bonds—gradually rotates “outward”. This outward orientation is favorable for establishing intermolecular hydrogen bonds with solvent molecules. Consequently, the population of the anti-conformer increases under high-pressure conditions, leading to a more solvent-exposed hydroxyl group configuration.

Fig. 12 (Left) Two-dimensional probability difference map illustrating the variation in the joint distribution of the dihedral angles φ6 (O5–C5–C6–O6) and τ6 (C5–C6–O6–H6) between P = 0.1 and 1000 MPa. The blue and red regions indicate decreased and increased conformer populations, respectively. (Right) Representative cellulose conformations labeled (A)–(D), corresponding to the specific regions marked in the left panel. The arrows show the orientation vector of the C6 hydroxyl group.

3.3.3.3 Summary and proposed mechanism for pressure-enhanced solubility at “microscale”. By integrating the results of the hydrogen bonding and conformational analyses of cellulose, the following microscopic molecular mechanism can be proposed to explain the pressure-enhanced solubility of cellulose (Fig. 13).

Fig. 13 Mechanism of pressure-induced enhancement of cellulose solubility. Hydrogen bonds switch from cellulose intramolecular hydrogen bonds (intra-Hbonds) to intermolecular hydrogen bonds (inter-Hbonds) between the cellulose and the anions.

At ambient pressure, cellulose molecules form intramolecular O6^H⋯O3 hydrogen bonds with a certain probability. However, under high pressure, the system equilibrium shifts toward states with smaller partial molar volumes in accordance with Le Chatelier's principle.

The conformational rotation of the hydroxymethyl (C6–OH) group, coupled with the penetration of [OAc]⁻ anions, leads to a more compact packing of the solvation shell, thereby reducing the total system volume (ΔV < 0). This volume-driven conformational change is analogous to pressure-induced protein denaturation, in which structural unfolding is favored to minimize the specific volume through enhanced hydration.

Crucially, this process is not driven by volume reduction alone but by a synergistic effect: the thermodynamic drive to reduce the volume brings the anion and cellulose into close contact, while the specific high hydrogen bond affinity of the anion energetically stabilizes the rotated conformation. As illustrated in Fig. 13, this leads to a concerted mechanism in which the disruption of the intramolecular bond and the formation of a stable intermolecular O6^H⋯O_anion bond occur simultaneously.

Therefore, the enhanced solubility of cellulose observed under high-pressure conditions is not merely a passive consequence of system compression, as described in Section 3.3.2, but can also be attributed to an active, conformationally driven process in the cellulose molecule that reorganizes the hydrogen-bonding network to favor solvation.

4 Conclusion

In this study, AA MD simulations were used to predict the solubility of cellulose in a 60 wt% [EMIm][OAc]/DMSO solvent over a wide pressure range (P = 0.1–1000 MPa). Comprehensive analyses were conducted to elucidate the molecular mechanisms by which pressure influences cellulose solubility.

First, umbrella sampling simulations were performed to calculate the dissolution free energy (ΔG_dis) and to predict the solubility of cellulose at higher pressures. The results revealed that ΔG_dis decreases monotonically as the pressure increases, which provides clear thermodynamic evidence that higher pressures enhance the solubility of cellulose in the [EMIm][OAc]/DMSO system.

The MD simulations of the dissolved cellulose system were analyzed in detail to elucidate the molecular origin of this behavior. The interaction energies, coordination numbers, and radial distribution functions showed that increasing the pressure strengthened cellulose–solvent interactions and weakened cellulose–cellulose interactions. Furthermore, analysis of the spatial distribution homogeneity demonstrated that the cellulose chains became more uniformly dispersed at higher pressures, forming a more homogeneous and well-solvated polymer solution. Together, these results indicate that high pressures reinforce ion–cellulose interactions, stabilize the solvation shells around each polymer chain, and simultaneously suppress cellulose–cellulose aggregation, thereby promoting dissolution under compression.

A deeper investigation into hydrogen bonding and conformational dynamics further revealed a more intricate active mechanism. Increasing the pressure-induced conformational transitions of the CH₂OH groups, which disrupted the critical intramolecular hydrogen bond (O6^H⋯O3) and facilitated the formation of new, favorable intermolecular hydrogen bonds between the cellulose and solvent anions.

Thus, pressure enhances cellulose solubility through passive densification and enhanced solvation, and active pressure-induced conformational remodeling of the hydrogen bond network thermodynamically favors dissolution.

From a practical perspective, improved solubility under pressure implies reduced solvent consumption, leading to potential economic and environmental benefits for cellulose processing. Overall, this study provides new fundamental insights into how pressure—a key thermodynamic parameter—governs solubility. The molecular-level understanding obtained in this study offers valuable guidance for designing next-generation biorefinery processes and sustainable material technologies. Specifically, the finding that cellulose–anion interactions exhibit the highest-pressure sensitivity suggests a clear strategy for solvent design: for high-pressure applications, ionic liquids containing anions with smaller steric bulk and higher hydrogen bond basicity, such as formate ([HCOO]⁻) or chloride ([Cl]⁻), can be selected to maximize the pressure-induced solvation enhancement.

Author contributions

Kodai Kikuchi: conceptualization, methodology, investigation, data curation, writing–original draft, funding acquisition. Kazushi Fujimoto: methodology, supervision, writing-reviewing. Kazuyoshi Kaneko: methodology, supervision, writing-reviewing, funding acquisition. Akio Shimizu: conceptualization, funding acquisition. Tatsushi Matsuyama: methodology. Junichi Ida: methodology, writing-reviewing and editing, supervision, funding acquisition.

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that may have influenced the work reported in this paper.

Data availability

The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5ra08753h.

Acknowledgements

Kodai Kikuchi is a Research Fellow of the Japan Society for the Promotion of Science (JSPS Fellow). This work was supported by JSPS KAKENHI, Grant Numbers: JP24KJ2119 and JP21K17910. We also thank the Research Center for Computational Science, Okazaki, Japan (Project: 24-IMS-C607, 25-IMS-C147, 24-IMS-C604, 24-IMS-C502, 25-IMS-C118, and 25-IMS-C502) for the use of supercomputers. We thank Editage (https://www.editage.com) for the English language editing.

References

1. D. B. Sulis, N. Lavoine, H. Sederoff, X. Jiang, B. M. Marques, K. Lan, C. Cofre-Vega, R. Barrangou and J. P. Wang, Nat. Commun., 2025, 16, 1244,  DOI:10.1038/s41467-025-56472-y.
2. Y. Yuan, Z. Li, L. Wu, X. Cheng, C. Deng, Y. Yu, Q. Wang and P. Wang, Biomaterials, 2025, 322, 123360.
3. N. Rameli, B.-Y. Lim, P.-Y. Leong, C.-C. Lim and S.-F. Ng, Macromol. Res., 2024, 32, 861–872,  DOI:10.1007/s13233-024-00266-6.
4. S. Dacrory, L. M. A. Ali, S. Ouahrani-Bettache, M. Daurat, M. El-Sakhawy, P. Hesemann, N. Bettache and S. Kamel, BMC Chem., 2025, 19, 52,  DOI:10.1186/s13065-025-01408-2.
5. P. Pawaskar, S. Lekurwale, M. Kumari, S. Roy and S. Banerjee, Ind. Crop. Prod., 2025,(229), 120998,  DOI:10.1016/j.indcrop.2025.120998.
6. K. Wang, Y. Wu, J. Remón and W. Ding, Int. J. Biol. Macromol., 2025, 314, 144433,  DOI:10.1016/j.ijbiomac.2025.144433.
7. H. Dong, J. Fa, M. Yan, Y. Tan, W. Cheng, J. Huang and Z. Ling, J. Colloid Interface Sci., 2025, 686, 701–710,  DOI:10.1016/j.jcis.2025.02.005.
8. M. H. Talari, N. S. Tabrizi, V. Babaeipour and F. Halek, J. Sol-Gel Sci. Technol., 2022, 101, 345–355,  DOI:10.1007/s10971-021-05708-6.
9. W. U. Khan, S. U. Khan, S. U. Khan, A. Gul, H. Zaman, W. U. Khan, N. Ali, A. Aziz, M. I. Khan, D. Khan and M. M. Gouda, Colloids Surf., A, 2025, 718, 136919,  DOI:10.1016/j.colsurfa.2025.136919.
10. S. Yi, X. Lv, X. Zhang, J. Chen, D. Fang, H. Jia, Y. Ma, W. Xu and X. Yang, Int. J. Biol. Macromol., 2025, 307(2), 141689,  DOI:10.1016/j.ijbiomac.2025.141689.
11. C. Yi, H.-Y. Niu, L. Sui, J.-J. Zhu, Y. Tian, C.-G. Niu, Z.-L. Chen, H. Wei and D.-W. Huang, Sep. Purif. Technol., 2025, 358(A), 130213,  DOI:10.1016/j.seppur.2024.130213.
12. J. Chen, X. Wang, X. Huang, Z. Tong, J. Zhou, Y. Shen and C. Hao, Int. J. Biol. Macromol., 2025, 306(1), 141398,  DOI:10.1016/j.ijbiomac.2025.141398.
13. P. Laddawan, S. Auecharoenrat and C. Tongsook, ACS Omega, 2025, 10(17), 17296–17309,  DOI:10.1021/acsomega.4c09928.
14. M. Salimi, A. El Idrissi, B.-e. Channab, Y. Essamlali, A. G. Firouzabadi, M. Beygi, M. Zahouily and E. Motamedi, Cellulose, 2024, 31, 10679–10726,  DOI:10.1007/s10570-024-06273-1.
15. J. Fu, J. X. Yap, C. P. Leo, C. K. Chang and P. L. Show, Sep. Purif. Rev., 2023, 53(3), 276–288,  DOI:10.1080/15422119.2023.2269395.
16. A. Bashari, A. R. Shirvan and M. Shakeri, Polym. Adv. Technol., 2018, 29, 2853–2867,  DOI:10.1002/pat.4290.
17. R. L. Wu, X. L. Wang, F. Li, H. Z. Li and Y. Z. Wang, Bioresour. Technol., 2009, 100, 2569–2574,  DOI:10.1016/j.biortech.2008.11.044.
18. H. Takagi and A. Asano, Composites, Part A, 2008, 39, 685–689,  DOI:10.1016/j.compositesa.2007.08.019.
19. W. Y. Li, A. X. Jin, C. F. Liu, R. C. Sun, A. P. Zhang and J. F. Kennedy, Carbohydr. Polym., 2009, 78, 389–395,  DOI:10.1016/j.carbpol.2009.04.028.
20. M. Kimura, Y. Shinohara, J. Takizawa, S. Ren, K. Sagisaka, Y. Lin, Y. Hattori and J. P. Hinestroza, Sci. Rep., 2015, 5, 16266,  DOI:10.1038/srep16266.
21. T. Schalk, S. Schubert, A. Rollberg, D. Freitag-Stechl, A. Schubert, A. X. Elena, C. Koch and P. Krebs, Water, 2023, 15(8), 1463,  DOI:10.3390/w15081463.
22. M. Sawiak, B. A. Souto, L. Lawson, J. Lo and P. I. Dolez, Fibers, 2025, 13(1), 3,  DOI:10.3390/fib13010003.
23. P. Marzullo, M. Gruttadauria and F. D'Anna, Biomolecules, 2024, 14(8), 957,  DOI:10.3390/biom14080957.
24. M. Sakamoto, Y. Ohama, S. Aoki, K. Fukushi, T. Mori, Y. Yoshimura and A. Shimizu, Aust. J. Chem., 2018, 71, 492–496,  DOI:10.1071/CH18117.
25. H. Wang, G. Gurau and R. D. Rogers, Chem. Soc. Rev., 2012, 41, 1519–1537,  10.1039/c2cs15311d.
26. Y. Tomimatsu, Y. Yoshimura and A. Shimizu, Aust. J. Chem., 2019, 72, 669–673,  DOI:10.1071/CH19047.
27. Y. Zhao, X. Liu, J. Wang and S. Zhang, J. Phys. Chem. B, 2013, 117(30), 9042–9049,  DOI:10.1021/jp4038039.
28. M. Moniruzzaman, H. Mahmood, M. F. Ibrahim, S. Yusup and Y. Uemura, Adv. Mater. Res., 2016, 1133, 588–592,  DOI:10.4028/www.scientific.net/amr.1133.588.
29. X. Wei, Y. Wang, J. Li, F. Wang, G. Chang, T. Fu and W. Zhou, Carbohydr. Polym., 2018, 201, 387–391,  DOI:10.1016/j.carbpol.2018.08.031.
30. M. T. Dignani, T. A. Bioni, T. R. L. C. Paixão and O. A. El Seoud, Molecules, 2020, 25(24), 5975,  DOI:10.3390/molecules25245975.
31. H.-C. Chang, R.-L. Zhang and D.-T. Hsu, Phys. Chem. Chem. Phys., 2015, 17, 27573–27578,  10.1039/C5CP04607F.
32. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25,  DOI:10.1016/j.softx.2015.06.001.
33. S. Boonstra, P. R. Onck and E. van der Giessen, J. Phys. Chem. B, 2016, 120, 3692–3698,  DOI:10.1021/acs.jpcb.6b01316.
34. S. Lee, A. Tran, M. Allsopp, J. B. Lim, J. Hénin and J. B. Klauda, J. Phys. Chem. B, 2014, 118, 547–556,  DOI:10.1021/jp410344g.
35. S. V. Sambasivarao and O. Acevedo, J. Chem. Theory Comput., 2009, 5, 1038–1050,  DOI:10.1021/ct900009a.
36. A. Jakalian, D. B. Jack and C. I. Bayly, J. Comput. Chem., 2002, 23, 1623–1641,  DOI:10.1002/jcc.10128.
37. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236,  DOI:10.1021/ja9621760.
38. J. Wang, W. Wang, P. A. Kollman and D. A. Case, J. Mol. Graph. Model., 2006, 25, 247–260,  DOI:10.1016/j.jmgm.2005.12.005.
39. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101,  DOI:10.1063/1.2408420.
40. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. Dinola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690,  DOI:10.1063/1.448118.
41. S. Nose and M. L. Klein, Mol. Phys., 1983, 50, 1055–1076,  DOI:10.1080/00268978300102851.
42. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190,  DOI:10.1063/1.328693.
43. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092,  DOI:10.1063/1.464397.
44. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, LINCS: A Linear Constraint Solver for Molecular Simulations, John Wiley & Sons, Inc, 1997.
45. J. M. Martinez and L. Martinez, J. Comput. Chem., 2003, 24, 819–825,  DOI:10.1002/jcc.10216.
46. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graph., 1996, 14(1), 33–38,  DOI:10.1016/0263-7855(96)00018-5.
47. B. Roux, Comput. Phys. Commun., 1995, 91, 275–282,  DOI:10.1016/0010-4655(95)00053-I.
48. J. S. Hub, J. Chem. Phys., 2015, 6, 3713–3720,  DOI:10.1021/ct100494z.
49. T. C. F. Gomes and M. S. Skaf, J. Comput. Chem., 2012, 33, 1338–1346,  DOI:10.1002/jcc.22959.
50. S. Rizvi and H. M. Gade, Biomass Bioenergy, 2025, 196, 107758,  DOI:10.1016/j.biombioe.2025.107758.
51. L. Nisa, M. Mugang, M. D. Bermejo, J. J. S. Puras and E. Pérez, J. Chem. Eng. Data, 2019, 64, 2923–2928,  DOI:10.1021/acs.jced.8b01048.
52. L. H. G. de Medeiros, A. A. A. Alves, F. X. Feitosa and H. B. de Santana, J. Chem. Eng. Data, 2022, 67, 1078–1088,  DOI:10.1021/acs.jced.2c00036.
53. L. de Pablo Nisa, J. J. Segovia, Á. Martín, M. C. Martín and M. D. Bermejo, J. Supercrit. Fluids, 2017, 130, 76–83,  DOI:10.1016/j.supflu.2017.07.040.
54. A. Brandt, M. J. Ray, T. Q. To, D. J. Leak, R. J. Murphy and T. Welton, Green Chem., 2011, 13, 2489–2499,  10.1039/C1GC15374A.
55. R. P. Swatloski, S. K. Spear, J. D. Holbrey and R. D. Rogers, J. Am. Chem. Soc., 2002, 124(18), 4974–4975,  DOI:10.1021/ja025790m.
56. Y. Fukaya, K. Hayashi, M. Wada and H. Ohno, Green Chem., 2008, 10, 44–46,  10.1039/B713289A.
57. H. Zhao, C. L. Jones, G. A. Baker, S. Xia, O. Olubajo and V. N. Person, J. Biotechnol., 2009, 139(1), 47–54,  DOI:10.1016/j.jbiotec.2008.08.009.
58. A. Xu, J. Wang and H. Wang, Green Chem., 2010, 12, 268–275,  10.1039/B916882F.
59. S. Padmanabhan, M. Kim, H. W. Blanch and J. M. Prausnitz, Fluid Phase Equilib., 2011, 309(1), 89–96,  DOI:10.1016/j.fluid.2011.06.034.
60. M. Zavrel, D. Bross, M. Funke, J. Büchs and A. C. Spiess, Bioresour. Technol., 2009, 100(9), 2580–2587,  DOI:10.1016/j.biortech.2008.11.052.
61. T. V. Doherty, M. Mora-Pale, S. E. Foley, R. J. Linhardt and J. S. Dordick, Green Chem., 2010, 12, 1967–1975,  10.1039/C0GC00206B.
62. M. Brehm, M. Pulst, J. Kressler and D. Sebastiani, Triazolium-based ionic liquids: A novel class of cellulose solvents, J. Phys. Chem. B, 2019, 123(18), 3994–4003,  DOI:10.1021/acs.jpcb.8b12082.
63. T. Uto, K. Yamamoto and J. Kadokawa, J. Phys. Chem. B, 2018, 122(1), 258–266,  DOI:10.1021/acs.jpcb.7b09525.
64. J. Wei, Y. Long, T. Li, H. Gao and Y. Nie, Int. J. Biol. Macromol., 2024, 267(1), 131204,  DOI:10.1016/j.ijbiomac.2024.131204.
65. T. G. A. Youngs, C. Hardacre and J. D. Holbrey, J. Phys. Chem. B, 2007, 111(49), 13765–13774,  DOI:10.1021/jp076728k.
66. K. Kikuchi, K. Fujimoto, A. Shimizu, K. Kaneko, T. Matsuyama and J. Ida, J. Mol. Liq., 2025, 437(A), 128344,  DOI:10.1016/j.molliq.2025.128344.

---