Introduction

At present, nanocellulose is being widely considered as an innovative material for various purposes and has become a research area of great interest. Cellulose nanofibers, which are comprised of β-D-glucopyranose chains and produced from sustainable resources, have been identified as an environment-friendly material with excellent physical and chemical characteristics1,2,3. They have become an important material due to their specific structure in which molecules form a chain. They can possess an extensive set of functionalities from the different characteristics of their glycosidic units as well as their atomic positions.

Nanocellulose substrates can be tailored to produce electronic devices for different types of applications4,5,6,7. During the processing and utilization of a device, it is known that temperature changes can affect its entire service life8. In this context, a more detailed understanding of the changes β-D-glucopyranose chains undergo during heat exposure in terms of their durability is needed to ascertain the thermal stability of cellulose nanomaterials.

The reaction of nanocellulose to heat depends on different factors, for example, the chemical composition of the material and the presence of oxidants, which affect both the surface and core structure of the material9. The dynamics of the structural changes of cellulose chains during thermal exposure could be better understood at elevated temperatures—around 190 °C. This generally is the boundary temperature before the reaction of the nanocellulose components starts the irreversible degradation process. Although different types of conformational studies10,11,12 have been performed to understand the conformational changes of the cellulose chains, they still face limitations.

The three-dimensional structure of cellulose nanofibers is one of the important factors determining its structure and properties. To gain insight into the mechanism of the thermal exposure of β-D glucopyranose, it is necessary to evaluate the conformational changes of the units. Exploration of the conformational space involves studying the molecular geometry13 of glucopyranose chains during the heat exposure.

In this study, we explored the complex of conformational transformation molecules of glucopyranose chains during thermal exposure and supported this with analysis of hydrogen conformers. The high-temperature resistant cellulose-based materials became an important feature when applied to the fabrication of flexible substrate for electronic devices. The high-temperature analysis of the structural and chemical transformations of cellulose nanofibers (CNF) and its effects on CNF properties were investigated in this study14,15. We applied computational modeling, especially molecular dynamics (MD) simulations, which led to insights into nanostructures and dynamics. In contrast to other computational research on force-fields16,17, our study focuses on the analysis of ring transformation dynamics as an important part of glucopyranose chain conformation during thermal exposure. During the simulations, conformation of the glucopyranose rings showed extensive movement of the middle part compared to the sides of the chain. In addition, in order to better understand the effects of conformation, three different types of cellulose nanofibers (CNF) were used in this study, each of which had different chemical compositions (Supplementary Table 1), to confirm its structural transformation. The dihedral angles for OH (2), OH (3), and OH (6) groups were increased, which promoted intrachain hydrogen bond18 transformation and led to geometrical changes of the glucopyranose chain (Supplementary Fig. 1). By using the coupling constant (J), the MD results of the structural transformation were confirmed by NMR study. The relationship between the coupling constant and the conformational characteristics of the glucopyranose ring is represented by the Karplus equation19. Analysis of the NMR spectrum indicated that the behavior of the chemical shift was affected by temperature, which follows a change in the vicinal coupling constant. CNFs from different raw materials sources (cotton, Kraft pulp, and bleached chemi-thermomechanical pulp) were used to confirm the structural transformation and observed differences based on the chemical composition of each material.

Results and discussion

Simulation

Three parts of the glucopyranose chain (both ends and middle part) were considered to understand the conformational changes during the thermal exposure. Each part of the chain is represented by two glucopyranose units (Supplementary Fig. 1a) that are linked together by 1,4-glycosidic bonds. By understanding the dynamic changes of the cellulose chain with the θ and φ polar variables, the conformational changes of the six units (Supplementary Fig. 1c) were reconstructed (Cremer-Pople (CP) parameters)20,21.

The CP parameter represented by the Mercator projection shows the position of the conformers in the southern and northern hemispheres. The dense concentration of the points shows a low mobility of the 1st, 2nd, and 5th units during the thermal exposure (as shown in Fig. 1). The difference in the angles θ and φ (θ—azimuthal and φ—meridian) before and after thermal exposure confirms the changes that occurred during the rotation of each glucopyranose ring (Supplementary Fig. 2). Therefore, the compliance (conformability) of units 3 and 4 was in the two parts of the hemisphere, which made the movement amplitude of the cellulose chains in the middle part higher. Angles θ and φ resulted in rearrangement of the glucopyranose rings 3 and 4 from 4C1 to 1C4 in the CP coordinates (representation of the dynamical changes of each unit can be found in Fig. 1). The intensive amplitude of conformation was also inherent to unit 6, which showed the highest magnitude of Δθ and Δφ (Supplementary Fig. 2). From these results, it is clear that the different conformational behaviors of the glucopyranose rings provide different types of stability to each part of the cellulose chain during thermal exposure. Such behavior of the conformational units can be explained by the interaction with the neighboring units and interaction around bonds of the rings. Since the side units are less connected to neighboring units, less movement can be observed compared to the central part of the chain. Upon further conformation, increasing structural changes may lead to the twisting transformation of the glucopyranose chain (both ends).

Fig. 1: Mercator representation of parameters for the six glucopyranose rings.
figure 1

Six glucopyranose rings during the conformation process a Ring 1, b Ring 2, c Ring 3, d Ring 4, e Ring 5, f Ring 6.

Characterization of glucopyranose units by NMR in conformational analysis

During thermal exposure of nanocellulose samples, partial dehydration may have occurred, which resulted in a double bond between the C2 and C3 positions and a OH group at C3 carbon (see Supplementary Fig. 3). 1H NMR spectroscopy was used to analyze the outcome of the thermal exposure of CNF films (Fig. 2). 1H NMR spectra of thermally treated CNF samples clearly indicated the presence of aldehyde groups at 10.1 and 9.5 ppm for each sample (see Fig. 2b). Also, peak sizes of OH (2) and OH (3) protons were reduced by up to 40% compared to original samples. Peaks at 4.65 and 4.5 ppm were attributed to the =C–OH group and =CH group, respectively, which clearly indicated the presence of an enol structure in the pyranose ring.

Fig. 2: 1H-NMR spectra of CNF (DMSO-d6).
figure 2

1H-NMR spectra of CNF a before and b after thermal exposure.

Moreover, the NMR spectra were recorded before and after the thermal exposure of CNF films at 190 °C and were designed to prove the conformational changes of the β-D-glucopyranose chain found from MD simulation results (Supplementary Fig. 1c). The flexibility of the chain (units in the chain) demonstrated a change in the coupling constant based on the observed chemical shift. We measured spectra and coupling constants before (at room temperature) and after heat at 190 °C that impacted the molecular bonding of glucopyranose units. Analysis of the NMR spectra provided insights into the chemical shifts and coupling constants for the OH (2), OH (3), and OH (6) groups.

Comparing the results of the samples before and after the heat exposure shows a decrease of the coupling constant (Fig. 3, Supplementary Table 2). The coupling constants provide insight into the conformational preferences22 of the backbone in the twisted cellulose chain. Compared with OH (2) and OH (6), the OH (3) proton signals of all 3 samples had a smaller coupling constant, which was specified by a smaller upfield migration. Comparing the coupling constants of the treated samples, it was found that the coupling constant of the OH (3) group was halved for all samples (after heat exposure, CN dropped by 50%, while CTMN dropped by 42.5%). The coupling constant of OH (6) decreased for CN by 7.9%, KN by 7.8%, and CTMN by 4.1% after thermal exposure. In contrast with CN and KN, CTMN samples had the smallest change after the heat exposure in their coupling constant.

Fig. 3: Coupling constant JHOCH.
figure 3

JHOCH for OH (2), OH (3), and OH (6) before and after thermal exposure. a CN, b KN, and c CTMN samples.

The proton coupling constant was used to evaluate the conversion of cellulosic rings23. Experimental JHOCH was applied to determine a θ dihedral angle by using the Karplus Eq. (1)19,24.

$$J_{{HOCH}} = 10.4\;{\rm{cos}}^2\theta - 1.5\;{\rm{cos}}\theta + 0.2$$
(1)

where θ—dihedral angle. Equation (1) provides estimation of the dihedral angle of the HOCH moiety in the conformation of OH (2), OH (3), and OH (6) groups. The absolute data of the calculated dihedral angle are shown in Supplementary Table 3. The decrease of the coupling constant with the thermal exposure increased the dihedral angle, which led to the structural changes of the cellulose chain. Observation of the temperature dependence of the OH group showed similar dihedral changes that follow a linear trend, as can be clearly observed in Fig. 4. Therefore, the strong correlation between the dihedral angles is mostly attributed to the active participation of OH groups in the conformation process during the thermal exposure. Dihedral angles plotted separately for each sample (Supplementary Fig. 4a–c) clearly indicate similar behavior for CN and KN samples. Interestingly, compared with CTMN and CN samples, the dihedral angles of KN were least affected by heat. The greatest changes were observed for OH (2) dihedral angles for CTMN samples, which were attributed to higher reactivity25 and more intensive conformation during the thermal exposure.

Fig. 4: Temperature dependence of dihedral angles.
figure 4

Temperature dependence of ψ and φ for OH (2), OH (3), and OH (6) based on coupling constants from experimental data for CTMN—1, CN—3, and KN—5 before and CTMNtr—2, CNtr—4, and KNtr—6 after thermal exposure.

The dihedral angles (ψ and φ) depend on temperature for OH (2), OH (3), and OH (6) and identified the hydrogen bond area according to different conformational directions. The correlation between coupling constants and dihedral angles is described by the parameterized Karplus curve (see Fig. 5). Figure 5 shows the distribution of the dihedral angles and provides the values for both hydrogen orientation directions during the conformations (Supplementary Fig. 5). For example, the dihedrals of OH (2) group for the CTMN sample, based on a coupling constant, span a range from ± 30° to ±135° before and ±49° to ±120° after the heating. This indicates a broad range of conformations during thermal exposure.

Fig. 5: Graph of the coupling constant mapped as a function of dihedral angles for CN, KN, and CTMN samples.
figure 5

CN, KN, CTMN indicate before and CNtr, KNtr, CTMNtr indicate after thermal exposure. Dashed line—clockwise direction.

Raman spectroscopy

The thermal reactivity of the CNF films is evidenced by the Raman spectra for the samples before and after the thermal exposure (Fig. 6). Based on the spectra, the medium peak intensity at 1710 cm−1 characterizes the C=O vibrations that were expected during the thermo-oxidative reaction and ketonic group formation (Fig. 6)26. Moreover, during the reactivity of KN and CTMN components, C=O stretching vibrations occurred between 1660 and 1683 cm−1 that connected to the lignin27,28,29. The signal intensity from 1470 cm−1 increased and belongs to COH mode assignment30. CH2 dominates as the contributing mode of the 1400–1500 cm−1 region. Moreover, the bands of a scissor mode at 1476 cm−1 are assigned to tg conformers and 1461 cm−1 to gg conformers of the CH2OH group31,32. At the same time, stretching vibrations between 1275 and 1346 cm−1 are assigned to the CH2 twisting that emphasizes the conformational changes of the CH2OH group during the thermal exposure26.

Fig. 6: Raman spectra of CN, KN, and CTMN CNF films.
figure 6

Raman spectra before and after the thermal exposure of CN, KN, and CTMN CNF films.

Figure 7a shows the TGA curves and derivative (DTG) curves of the samples at a heating rate of 10 °C min−1. Compared with lignin, a large amount of hemicellulose reduces the initial decomposition temperature of CTMN. However, the presence of lignin and its rigid nature hindered the thermal degradation since the well-packed and more entangled chains prevented heat diffusion through the sample and increased the thermal stability of the lignin-containing samples (Supplementary Fig. 6).

Fig. 7: Kinetic analysis of CNF films using the KAS method.
figure 7

a TG and DTG profiles of CN, KN, and CTMN CNF; b activation energy KAS.

As shown in Table 1, the mean activation energies calculated using the Kissinger method for CN, KN, and CTMN were 119.96, 144.21, and 164.55, respectively. The higher the lignin content of the samples at the same conversion rate, the greater the activation energy values. This indicated that the higher the lignin content, the larger the energy that was required for thermal degradation, and consequently the better the thermal stability, which is in agreement with the results of TGA. Figure 7b shows that at higher conversions (>0.80), the thermal degradation of the samples followed multiple reaction mechanisms due to the complex decomposition of lignin33. This complexity can be observed by the significant standard errors and poor correlation factors in this region.

Table 1 Effective activation energy for CN, KN, and CTMN obtained by the KAS isoconversional method.

Supplementary Fig. 7 shows the application results of the KAS methods with conversions from 0.01 to 0.99. These plots were used to provide the apparent activation energy values and revealed thermal degradation mechanisms for different stages of solid-phase reactions at different values of conversion. The parallelism of the lines is associated with only one type of dominant reaction mechanism in the degradation process. In this context, we focus on the conversion in the pseudo cellulose region (from approximately 0.1 to 0.80–0.84) instead of the entire process, since this range can provide important information for thermal decomposition kinetic modeling of CNF-based films. The fitted lines at different conversions for CN, KN, and CTMN in this mentioned range are almost parallel with each other, suggesting one reaction mechanism for pseudo cellulose. The reaction mechanism changed at higher conversion values since the lines were not parallel. This is an indication of the complexity of the decomposition mechanism, including parallel and competitive reactions. Spaced lines at lower conversion degrees (<0.07) were also observed, which were assigned to the dehydration process. Interestingly, noticeable lines were observed from 0.08 and 0.11 for CTMN because of hemicellulose decomposition.

For KN and CTMN, due to the presence of lignin, its structural complexity, and its interactions with the other biomass components, the decomposition mechanism was found to be different when compared with KN. For the KN sample, Fig. 8 shows that random scission kinetic models (L2) better approximate the experimental curves than the nucleation models34 obtained for the cellulose compound. Master plots in Fig. 8 show that the thermal decomposition of CN is in agreement with a random scission which is related to the arbitrary break of polymer chains into smaller segments. On the other hand, for KN, the experimental curve is close to the F1 master curve, indicating the validity of the mechanism of random nucleation with one nucleus on the individual particles. This observation is similar to the observation for wood pyrolysis35.

Fig. 8: Master-plot functions.
figure 8

Theoretical and experimental master-plot functions for different stages at 10 °C min−1 according to the 11 most common models, where a CN, b KN, c CTMN stage I, and d CTMN stage II.

Based on these observations it is possible to infer that the degradation initiates from random points acting as growth centers for the propagation of degradation reactions. After the removal of inherent moisture, the degradation of KN proceeds through the rupture of bonds of long-chain cellulose and hemicellulose molecules which result in low molecular weight and act as sites for random nucleation and growth. For CTMN film, kinetic analysis of the hemicellulose and cellulose pseudo components was also performed in this study. As is known, hemicellulose is predominantly amorphous which may cause better transfer of heat to the interior of the particles resulting in the formation of nucleation sites36. The decomposition stage (stage 1) associated with hemicellulose degradation suggested degradation behavior corresponds to random nucleation with three nuclei on an individual particle proceeded through the two-dimensional diffusion mechanism during cellulose degradation (stage 2). This behavior may be associated with the highest degradation of the extractives, and the hemicellulose content can lead to higher volatility of the main wood components at relatively lower temperatures than cellulose alone, which leads to the degradation of the cellulose by diffusion. The D2 mechanism (stage II) indicates that thermolysis of CTMN is a diffusion-limited process owing to its compact structure37.

These results indicate that the thermal stability of CNF will expand the application area of nanocellulose composite materials to advanced electronics and medical devices. High-temperature stability can optimize the manufacturing of CNF film substrates during device processing. As well, the flexible capability of the glucopyranose chain makes it possible to use CNF in user-interactive electronic skin (flexible sensors) applications and as a green reinforcement for lightweight multifunctional materials. The lightweight multifunctional materials fabrication is now more possible based on the thermal stability of CNF and the application of different types of techniques, including hot inks deposition and 3D printing.

In summary, structural changes of the glucopyranose chain and the chemical behavior response of nanocellulose during thermal exposure were described in this study. Three types of CNF films that differ by chemical composition demonstrated conformational reactivity of OH groups (OH (2), OH (3), and OH (6)) during thermal exposure.

The conformational behavior of the glucopyranose rings provided different types of stability for each part of the cellulose chain during thermal exposure. The dynamical change trajectory for the rings in the middle part of the chain showed a transition path from 4C1 to 1C4 states (chair to invert-chair). Conforming rings on the ends of the glucopyranose chain were less dynamically active with the exception of ring 6. The reactivity of the components during thermal exposure was studied by NMR and RAMAN analysis, which confirmed the formation of an enol structure. Moreover, NMR studies demonstrated the effect of thermal exposure on the coupling constant of each CNF sample and the change of the dihedral angle of the glucopyranose unit. Increasing dihedral angles resulted in conformational changes, mostly in the middle part of the chain, where the impact was the highest. Comparison of the dihedral angles for CNF samples indicates that KN and CN samples were more conformationally stable compared to CTMN samples. A kinetic study proved a strong relationship between the chemical composition of the CNF samples and various stages of the reaction, supported by TGA data. This work demonstrated the thermal reactivity and stability of the glucopyranose chain using the structural conformational “response” of the nanocellulose to high-temperature exposure.

Methods

Film preparation

Cotton linters, softwood kraft pulp, and chemi-thermomechanical pulp were used to produce CNFs to be carried through the mechanical fibrillation process (Supplementary Table 1).

CNFs from a chosen pulp were produced by using a Supermasscolloider. A CNF suspension with 2–2.3% of fibers was diluted to 0.2% consistency and vacuum filtered by using a polyethersulfone filter. After filtration, the CNF film was formed and taken out of the funnel (Supplementary Fig. 8). Subsequently, the CNF film was sandwiched between two filtering papers to absorb water and then placed under a pressure of 344 kPa for 30–45 min. The formed CNF sheet was kept overnight at room temperature. CN CNF film from cotton pulp, KN CNF film from kraft pulp, and CTMN CNF film from BCTMP were obtained. Later, the prepared samples were heated in an oven at 190 °C for 5 h for further study.

Simulation details

Classical molecular dynamics simulations were conducted using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)38. The ReaxFF potential was used to describe three-dimensional changes during the thermal equilibration39. To integrate the equation of motion, a time step of 0.25 fs was used. The glucopyranose chain was heated from 1 to 463 K in a span of 0.125 ps with the application of an initial random Gaussian velocity distribution. These steps were followed at 463 K for 9.0 ns to achieve a thermally equilibrated state. It is necessary to mention that oxygen atoms, as a part of intrachain hydrogen bonds, were fixed to prevent breaking of the glucopyranose chain during the heating process. An isothermal-isobaric ensemble was applied to achieve structure relaxation. Periodic boundary conditions were employed along the three orthogonal directions.

Characterization

Structural characterization of CNFs before and after thermal exposure was carried out by NMR (nuclear magnetic resonance) spectroscopy. The NMR spectrum was obtained from a 500 MHz Agilent DD2 NMR Spectrometer (with 2 channel system). 1H spectra were received from CNF samples dissolved in DMSO-d6 solvent. The 1H spectrum was recorded at 25 °C; a 45-degree pulse flipping angle, 6.000 s acquisition times, and a 1.0 s relaxation delay time were used.

Raman spectra were collected on a Bruker SENTERRA Raman microscope using 785 nm diode laser excitation equipped with a camera to select the spots to be analyzed. The Bruker OPUS software program was used to process Raman data for normalization, baseline correction, and peak position determination.

Thermogravimetric analysis (TGA) was carried out using TGA Q500, from 25 to 800 °C at four heating rates of 10, 15, 20, and 30 °C/min under nitrogen flow.

Differential scanning calorimetry (DSC) was conducted using DSC Q1000, from 25 to 300 °C with a heating rate of 10 °C/min in N2.

Kinetics

The thermal stability and reaction mechanism of the samples were evaluated by using model-free and model-based kinetic methods, through the application of state-of-the-art kinetic software Kinetics NEO (NETZSCH-Gerätebau GmbH). It helps to determine the number of reaction steps and to analyze the temperature-dependent chemical process. The Kissinger-Akahira-Sunose (KAS) isoconversional and Master-plot methods were also utilized to estimate the activation energy of thermal degradation and investigate the governing degradation mechanism of a solid-state reaction. Using a non-isothermal thermogravimetric analysis (Thermal Gravimetric Analyzer – TGA, TA TGA-Q500, TA Instruments, USA), after determining the activation energy values of the reaction, the reaction mechanism was determined by plotting master plots40. The model-free results from approximately 0.1 to 0.80–0.84, instead of the entire process, were employed to obtain the master plots, usually wherever the lowest standard deviations were present. This range can provide important information about reactivity and thermal decomposition kinetic modeling of films.