# Advancing plant cell wall modelling: Atomistic insights into cellulose, disordered cellulose, and hemicelluloses – A review

Ali Khodayari, Ulrich Hirn, Stefan Spirk, Yu Ogawa, David Seveno, Wim Thielemans

► To cite this version:

Ali Khodayari, Ulrich Hirn, Stefan Spirk, Yu Ogawa, David Seveno, et al., Advancing plant cell wall modelling: Atomistic insights into cellulose, disordered cellulose, and hemicelluloses – A review. Carbohydrate Polymers, 2024, 343, pp.122415. 10.1016/j.carbpol.2024.122415. hal-04877525

HAL Id: hal-04877525

https://hal.science/hal-04877525v1

Submitted on 9 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Advancing plant cell wall modelling: Atomistic insights into cellulose, disordered cellulose, and hemicelluloses – A review

Ali Khodayari<sup>a, *</sup>, Ulrich Hirn<sup>b</sup>, Stefan Spirk<sup>b</sup>, Yu Ogawa<sup>c</sup>, David Seveno<sup>a</sup>, Wim Thielemans<sup>d</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

### Keywords:

Molecular modelling

Cellulose

Hemicellulose

Non-crystalline cellulose

Interactions

## ABSTRACT

The complexity of plant cell walls on different hierarchical levels still impedes the detailed understanding of biosynthetic pathways, interferes with processing in industry and finally limits applicability of cellulose materials. While there exist many challenges to readily accessing these hierarchies at (sub-) angström resolution, the development of advanced computational methods has the potential to unravel important questions in this field. Here, we summarize the contributions of molecular dynamics simulations in advancing the understanding of the physico-chemical properties of natural fibres. We aim to present a comprehensive view of the advancements and insights gained from molecular dynamics simulations in the field of carbohydrate polymers research. The review holds immense value as a vital reference for researchers seeking to undertake atomistic simulations of plant cell wall constituents. Its significance extends beyond the realm of molecular modeling and chemistry, as it offers a pathway to develop a more profound comprehension of plant cell wall chemistry, interactions, and behavior. By delving into these fundamental aspects, the review provides invaluable insights into future perspectives for exploration. Researchers within the molecular modeling and carbohydrates community can greatly benefit from this resource, enabling them to make significant strides in unraveling the intricacies of plant cell wall dynamics.

## 1. Introduction

In recent years, there has been a tremendous increase in using computational methods to enrich our view on cellulose and hemicellulose properties at the nanoscale. The overarching goal is understanding the structure and the interplay of the constituents in the plant cell wall on different hierarchical levels (Fadda, 2022; Pieczywek, Chibrikov, & Zdunek, 2023). However, even basic questions regarding the nanostructure of the cell wall features remain unresolved. Research has been focusing on providing proper responses to such questions, but some have still remained intact. As an example, the ultrastructure of the cellulose microfibrils in lignocellulosic plant cell walls is still not well-understood. Microfibrils, often referred to as elementary fibrils, are of tremendous interest in materials research as they form the building blocks for

nanocellulose, one of the few naturally occurring nanomaterials. Elucidating the cross section of nanocellulose, defined by the number of chains and their configuration, is of utmost importance as it defines surface interactions and mechanical characteristics of the nanocrystal.

Another example is the twist of nanocellulose, which is slowly being unraveled. The twisting of cellulose chains and the nanoparticles they constitute has many implications for their self-assembly into chiral nematic materials (Hu & Abidi, 2016; Majoinen, Kontturi, Ikkala, & Gray, 2012), their use in optical materials (Kose, Tran, Lewis, Hamad, & MacLachlan, 2019), and in catalysis (Kaushik et al., 2015) where nanocellulose can be applied as an auxiliary chiral reagent, avoiding the need for chiral catalysts. The twist may also be involved in the generation of longitudinal dislocations in nanocellulose, i.e., small areas covering a few nanometers where the cellulose crystallinity in

<footnote><sup>a</sup> Department of Materials Engineering, KU Leuven, Kasteelpark Arenberg 44, Leuven 3001, Belgium</footnote>

<footnote><sup>b</sup> Institute of Bioproducts and Paper Technology, TU Graz, Inffeldgasse 23, Graz 8010, Austria</footnote>

<footnote><sup>c</sup> Centre de recherches sur les macromolécules végétales, CERMAV-CNRS, CS40700, 38041 Grenoble cedex 9, France</footnote>

<footnote><sup>d</sup> Sustainable Materials Lab, Department of Chemical Engineering, KU Leuven, Campus Kulak Kortrijk, Etienne Sabbelaan 53, 8500 Kortrijk, Belgium</footnote>

<footnote>* Corresponding author.
E-mail address: ali.khodayari@kuleuven.be (A. Khodayari).</footnote>
elementary fibrils is broken.

Computational studies could, therefore, provide an alternative view, added to the pre-existing experimental studies, to better understand similar phenomena, as the above-mentioned examples stated.<fnref n="1" /> Attempts to find an interplay between computational studies and experimental investigations have been become popular in recent years (Penttilä, Paajanen, & Ketoja, 2021). One example would be finding the best agreement between several models of the secondary cell walls and results obtained through solid-state Nuclear Magnetic Resonance (ssNMR) (Addison et al., 2024). Same concepts have also been valorized through combination of quantum mechanics and ssNMR (Kirui et al., 2022).

Larger scale simulations, e.g., coarse-grained molecular dynamics, have made studying such an interplay between the computational and experimental realm much easier (Zhang, Yu, et al., 2021).

In this work we are reviewing the state of molecular dynamics (MD) research on some of the plant cell wall constituents. The importance of understanding the properties of cellulose and hemicelluloses and their interactions on a molecular level cannot be overstated. We are, hence, compiling a bottom-up perspective on the topic, starting with the smaller size scales and moving up from there. The work ultimately aims to present a comprehensive and critical physico-chemical view of the contribution of molecular dynamics simulations to the field of studying the plant cell wall, assessing the state of the research and the current limitations.

It is worth noting that we have avoided including a standalone section on lignin in the review, although we do consider its interactions with other substances in the plant cell wall. The reason is mainly that lignin is a highly complex and heterogeneous macromolecule, and its structural and chemical features are still not fully understood. This poses a significant challenge for molecular modeling studies, as the accuracy of the results strongly depends on the quality of the force field parameters and the experimental data used as a reference. The limitations posed by such complications have slowed progress on modeling lignin considerably. Additionally, lignin interactions with cellulose and hemicellulose are also highly dependent on the plant species, tissue type, and environmental conditions, making it difficult to draw general conclusions. Furthermore, chemical treatments performed on lignin during extraction alter its structure from how it is within the plant cell wall. Therefore, while lignin is undoubtedly a crucial component of the plant cell wall, we believe that a comprehensive review of its molecular dynamics simulations warrants a separate and focused effort.

To aid better discussions, some terminologies are clarified below, as they are regularly used in this review.

**Cellulose nanocrystals (CNC)**: Smallest crystalline segments isolated from cellulose fibres, typically with a thickness of 3–4 nm. The length of CNCs can extend up to 150 nm, although this scale is primarily influenced by whether one considers individual monocrystals or an aggregate of them.

**Cellulose microfibrils (CMF)**: Connection of crystalline segments via disordered cellulose regions. Microfibrils mainly have the same diameter as the CNCs while having lengths that can grow up to µm ranges.

**Dislocations**: Disordered regions in the CMFs, along the length of the fibrils.

**Cellulose macrofibrils**: Aggregated CMFs can build up large clusters with diameters in the range of tens of nm. Cellulose macrofibrils connected by a hemicellulose-lignin matrix are comprising the fibre cell walls.

**Plant cell walls**: The primary reason of plant growth, made of distinct layers, including the primary wall, secondary walls (S1, S2, S3 with S2 being the predominant layer), and a central lumen. Depending on the plant source, the number of layers, the percentage of the constituents, and their type could differ (Bourmaud, Beaugrand, Shah, Placet, & Bailey, 2018). Plant cell walls are believed to be the source of strength of the plant, intercellular communication, and water movement in the plants (Cosgrove, 2005; Ishii, Thomas, Darvill, & Albersheim, 1989).

**Plant fibres**: A complex hierarchical structure made of cellulose, hemicelluloses, lignin, pectin, water (depending on the source), formed by aggregates of several macrofibrils. These fibrous materials are mainly obtained from stems, leaves, or seeds of plants.

## 2. Force fields

Various computational methods contribute to our understanding of materials, including Density Functional Theory (DFT), tight-binding, and less detailed methods such as rigid body simulations or finite element methods. However, each method has its inherent limitations. Conventional ab-initio methods, for instance, are usually constrained by specific size scales, making them less suitable for larger simulations that demand molecular precision, with the exception of recent developments in the quantum chemistry software, such as TeraChem. In contrast, molecular dynamics offer a middle ground, providing a broader length-time scale that is more accurate than continuum scale computations. Nevertheless, each method introduces its own advantages and disadvantages. An example would be non-realistic (ballistic-like) strain rates usually observed in MD simulations, which is not a matter of concern in continuum scale simulations.

Since the primary implementation of molecular dynamics solvers involves solving Newton's equations of motion for all particles in the system, the distinctions in simulation works typically arise from the specific parameter sets used to define interactions between atoms.

Force fields, comprising a set of parameters for modeling interactions among all atoms in a simulated system, play a crucial role. Atomistic simulations, in particular, categorize atoms (in pairs or more) into bonded and non-bonded modes. Interactions are characterized by various equations expressing distinct types of interatomic potentials, such as van der Waals (vdW) and Coulombic non-bonded terms. Force fields encompass parameters for defining bonds, angular potentials among sets of three atoms, and interactions involving two planes, such as dihedral angles.

Several parameter sets have been developed throughout the years to model cellulose, and mainly carbohydrates. A summary of the force fields and simulation packages used in the literature is brought in Table 1. The aim of preparing such a comprehensive list of constitutes, force fields, and packages is to provide the readers with a database including all possible combinations of packages and force fields able to model specific constituents. Obviously, the compatibility of different packages and force fields would always be modulated through continuous development of the molecular dynamics resources. Noteworthy is that the table shows certain packages/force fields are more favoured in modeling specific constituents.

Force fields, in general, are mainly divided into two categories: reactive and non-reactive. Reactive force fields are primarily able to model the formation and breakage of the covalent bonds (Aktulga, Pandit, Van Duin, & Grama, 2012), whereas, in the non-reactive force fields, only non-bonded terms can form or break, and the covalent bonds are modelled as springs able to be stretched as long as they are asked to. This is in principle because the bonded terms in non-reactive force fields are defined as, for example:

<footnote n="1">¹ To better understand the behavior of cellulose nanocrystals, we can consider a specific case. While the precise ultrastructure of these nanocrystals is still unknown, experimental studies could provide insight into their hydration properties. For example, the sorption isotherms (which measure the amount of water absorbed at different humidity levels) of cellulose nanocrystals could be experimentally measured. Computational modeling could then be performed on various proposed ultrastructures, and the simulated sorption isotherms compared to the experimental data. By relating the computational and experimental findings in this way, we may be able to infer the most likely arrangement of the nanocrystals.</footnote>
<table><caption>**Table 1**
The plant cell wall constituents modelled with different packages and force fields.</caption><thead><tr><th>Constituent type</th><th>Force field</th><th>Package</th><th>Reference</th></tr></thead><tbody><tr><td>Crystalline cellulose</td><td>Cerius2/ Discover</td><td>CVFF91</td><td>(Mazeau & Vergelati, 2002)</td></tr><tr><td>Crystalline cellulose</td><td>CHARMM</td><td>CHARMM</td><td>(Besombes & Mazeau, 2005a, 2005b; Bu, Himmel, & Crowley, 2015; Ciesielski et al., 2019; Matthews et al., 2006)</td></tr><tr><td>Crystalline cellulose</td><td>CHARMM36</td><td>GROMACS</td><td>(Chae et al., 2023; Fujisawa, Daicho, Yurtsever, Fukuma, & Saito, 2023; Ma et al., 2021; Nascimento et al., 2022; O'Neill et al., 2017; Sridhar, Berglund, & Wohlert, 2023; Thu, Moreira, Weber, & Poma, 2022; Tolonen, Bergensträhle-Wohlert, Sixta, & Wohlert, 2015)</td></tr><tr><td>Crystalline cellulose</td><td>CHARMM36</td><td>NAMD</td><td>(Izumi, Saitoh, Sato, Takuma, & Takahashi, 2023; Malaspina & Faraudo, 2019; Oehme et al., 2015; Oehme, Yang, & Kubicki, 2018)</td></tr><tr><td>Crystalline cellulose</td><td>CHARMM36</td><td>NAMD 2.9</td><td>(Kannam et al., 2017; Oehme et al., 2015)</td></tr><tr><td>Crystalline cellulose</td><td>CHARMM36</td><td>LAMMPS</td><td>(Bregado, Tavares, Secchi, & Segtovich, 2021; He, Wu, Xia, Hou, & Zhu, 2023; Ito & Matsumoto, 2022)</td></tr><tr><td>Crystalline cellulose</td><td>COMPASS</td><td>Materials Studio</td><td>(Deng, Wu, & Wang, 2020; Eichhorn & Davies, 2006; Eichhorn, Young, & Davies, 2005)</td></tr><tr><td>Crystalline cellulose</td><td>DREIDING</td><td>LAMMPS</td><td>(Säenz Ezquerro, Crespo Minana, Izquierdo, & Laspalas, 2019)</td></tr><tr><td>Crystalline cellulose</td><td>GLYCAM</td><td>AMBER</td><td>(Yui, Nishimura, Akiba, & Hayashi, 2006)</td></tr><tr><td>Crystalline cellulose</td><td>GLYCAM04</td><td>AMBER8</td><td>(Nishiyama, Johnson, French, Forsyth, & Langan, 2008; Yui & Hayashi, 2007)</td></tr><tr><td>Crystalline cellulose</td><td>GLYCAM 04/06</td><td>AMBER 9/10</td><td>(Nishiyama, Johnson, & French, 2012)</td></tr><tr><td>Crystalline cellulose</td><td>GLYCAM06</td><td>AMBER</td><td>(Chen, Jiang, Wu, & Zheng, 2020; Jiang, Chen, Yuan, & Zheng, 2019)</td></tr><tr><td>Crystalline cellulose</td><td>GLYCAM06</td><td>AMBER12</td><td>(Hadden, French, & Woods, 2013)</td></tr></tbody></table>

<table><caption>**Table 1 (continued)**</caption><thead><tr><th>Constituent type</th><th>Force field</th><th>Package</th><th>Reference</th></tr></thead><tbody><tr><td>Crystalline cellulose</td><td>GLYCAM06</td><td>LAMMPS</td><td>(Bering, Torstensen, Lervik, & de Wijn, 2022)</td></tr><tr><td>Crystalline cellulose</td><td>GLYCAM06</td><td>DL_POLY</td><td>(Kong et al., 2023; Maurer, Sax, & Ribitsch, 2013)</td></tr><tr><td>Crystalline cellulose</td><td>GLYCAM06</td><td>GROMACS</td><td>(Chen, Lo Re, Berglund, & Wohlert, 2020; Chen, Wohlert, Berglund, & Furó, 2022; Djahedi, Bergensträhle-Wohlert, Berglund, & Wohlert, 2016; Hadden, French, & Woods, 2014; Ishida, 2020; Karna, Wohlert, Lidén, Mattsson, & Theliander, 2021; Lindh, Bergensträhle-Wohlert, Terenzi, Salmén, & Furó, 2016; Pang, Mehandzhiyski, & Zozoulenko, 2023; Wohlert, Chen, Berglund, & Lo Re, 2023; Zhang, Brumer, Ågren, & Tu, 2011; Zhang, Bulone, Ågren, & Tu, 2011)</td></tr><tr><td>Crystalline cellulose</td><td>GLYCAM06</td><td>NAMD</td><td>(Wada et al., 2011)</td></tr><tr><td>Crystalline cellulose</td><td>GLYCAM06<br>CHARMM35</td><td>AMBER</td><td>(Matthews et al., 2011)</td></tr><tr><td>Crystalline cellulose</td><td>GLYCAM06<br>COMPASS<br>ReaxFF</td><td>Materials Studio</td><td>(Dri, Wu, Moon, Martini, & Zavattieri, 2015)</td></tr><tr><td>Crystalline cellulose</td><td>GROMOS 45a4</td><td>GROMACS</td><td>(Bergensträhle, Berglund, & Mazeau, 2007; Bergensträhle, Thormann, Nordgren, & Berglund, 2009; Wohlert, Bergensträhle-Wohlert, & Berglund, 2012)</td></tr><tr><td>Crystalline cellulose</td><td>GROMOS 53a6</td><td>GROMACS</td><td>(Chen, Nishiyama, & Mazeau, 2012, 2014; Zhang, Keten, Derome, & Carmeliet, 2021)</td></tr><tr><td>Crystalline cellulose</td><td>GROMOS 54a7</td><td>GROMACS</td><td>(Lombardo et al., 2018)</td></tr><tr><td>Crystalline cellulose</td><td>GROMOS<br>56A<sub>carbo</sub></td><td>GROMACS</td><td>(Chen et al., 2014; Molnár et al., 2018; Ogawa, Nishiyama, & Mazeau, 2020)</td></tr><tr><td>Crystalline cellulose</td><td>GROMOS87</td><td>GROMOS87</td><td>(Heiner, Kuutti, & Teleman, 1998; Heiner, Sugiyama, & Teleman, 1995; Heiner & Teleman, 1996, 1997)</td></tr><tr><td>Crystalline cellulose</td><td>OPLS</td><td>NAMD</td><td>(Sánchez-Badillo, Gallo, Rutiaga-Quinones, &amp; Sánchez-Badillo, 2015)</td></tr></tbody></table>

(continued on next page)
<table><caption>**Table 1 (continued)**</caption><thead><tr><th>Constituent type</th><th>Force field</th><th>Package</th><th>Reference</th></tr></thead><tbody><tr><td>Crystalline cellulose</td><td>OPLS-AA</td><td>GROMACS</td><td>López-Albarrán, 2021)<br>(Garg et al., 2021; Garg, Linares, & Zozoulenko, 2020; Muthoka et al., 2020; Muthoka, Panicker, & Kim, 2022; Muthoka, Shishir, Kim, Kim, & Kim, 2018; Paajanen, Ceccherini, Maloney, & Ketoja, 2019; Paavilainen et al., 2012; Paavilainen, Rög, & Vattulainen, 2011; Pang, Mehandzhiyski, & Zozoulenko, 2022)<br>(Hardy & Sarko, 1996)</td></tr><tr><td>Crystalline cellulose</td><td>PARM20</td><td>CHARMM</td><td>(Mazeau & Rivet, 2008)</td></tr><tr><td>Crystalline cellulose</td><td>PCFF</td><td>Cerius</td><td>(Zhu, Gu, & Zhu, 2020)</td></tr><tr><td>Crystalline cellulose</td><td>ReaxFF</td><td>Materials Studio LAMMPS</td><td>(Diaz, Wu, Martini, Youngblood, & Moon, 2013; Geng et al., 2023; Hou et al., 2021; Wu, Moon, & Martini, 2013, 2014)<br>(Xia, Qin, Zhang, Sinko, & Keten, 2018)</td></tr><tr><td>Crystalline cellulose, Amorphous cellulose</td><td>CHARMM</td><td>NAMD</td><td>(Kulasinski, Keten, Churakov, Derome, & Carmeliet, 2014)<br>(Mazeau & Heux, 2003)</td></tr><tr><td>Crystalline cellulose, Amorphous cellulose</td><td>PCFF</td><td>–</td><td>(Khodayari, Hirn, Spirk, Van Vuure, & Seveno, 2021; Khodayari, Hirn, Van Vuure, & Seveno, 2020; Khodayari, Van Vuure, Hirn, & Seveno, 2020)<br>(Chen, Zhang, & Ke, 2023)</td></tr><tr><td>Crystalline cellulose, Dislocated cellulose</td><td>GLYCAM06</td><td>GROMACS</td><td>(Youssefian & Rahbar, 2015)</td></tr><tr><td>Crystalline cellulose, Dislocated cellulose</td><td>PCFF</td><td>LAMMPS</td><td>(Zhao, Crespi, Kubicki, Cosgrove, & Zhong, 2014)<br>(Jaafar et al., 2019)</td></tr><tr><td>Crystalline cellulose, Amorphous cellulose, Hemicellulose, Lignin</td><td>CHARMM36</td><td>CHARMM</td><td>(Heinonen, Henriksson, Lindström, Vilaplana, & Wohlert, 2022)<br>(Mazeau & Charlier, 2012)</td></tr><tr><td>Crystalline cellulose, Xyloglucan</td><td>Universal force field</td><td>Materials Studio</td><td>(Busse-Wicher et al., 2014; Kong, Li, & Fu, 2022)</td></tr><tr><td>Crystalline cellulose, Xylan</td><td>Universal force field</td><td>Cerius2</td><td></td></tr><tr><td>Crystalline cellulose, Xylan</td><td>CHARMM</td><td>NAMD</td><td></td></tr></tbody></table>

<table><caption>**Table 1 (continued)**</caption><thead><tr><th>Constituent type</th><th>Force field</th><th>Package</th><th>Reference</th></tr></thead><tbody><tr><td>Crystalline cellulose, Xylan</td><td>CHARMM</td><td>GROMACS</td><td>(Gupta, Rawal, Dupree, Smith, & Petridis, 2021)</td></tr><tr><td>Crystalline cellulose, Arabinoxylan</td><td>CHARMM36</td><td>GROMACS, NAMD</td><td>(Shrestha et al., 2019)</td></tr><tr><td>Crystalline cellulose, Glucuronoarabinoxylan</td><td>CHARMM</td><td>NAMD</td><td>(Busse-Wicher et al., 2016)</td></tr><tr><td>Crystalline cellulose, Galactoglucomannan, Glucuronoarabinoxylan</td><td>GLYCAM06</td><td>GROMACS</td><td>(Paajanen, Zitting, Rautkari, Ketoja, & Penttilä, 2022; Zitting, Paajanen, & Penttilä, 2023; Zitting, Paajanen, Rautkari, & Penttilä, 2021)</td></tr><tr><td>Crystalline cellulose, Galactoglucomannan, Glucuronoarabinoxylan, Xyloglucan</td><td>GLYCAM06</td><td>GROMACS</td><td>(Khodayari, Thielemans, Hirn, Vuure, & Seveno, 2021)</td></tr><tr><td>Crystalline cellulose, Galactoglucomannan</td><td>CHARMM36</td><td>NAMD</td><td>(Yu et al., 2018)</td></tr><tr><td>Crystalline cellulose, Hemicelluloses</td><td>CHARMM</td><td>GROMACS</td><td>(Kumar et al., 2018)</td></tr><tr><td>Crystalline cellulose, Hemicelluloses</td><td>GLYCAM06</td><td>AMBER</td><td>(Zhang, Li, Xiong, Hong, & Chen, 2015)</td></tr><tr><td>Crystalline cellulose, Hemicelluloses</td><td>GROMOS 53a6</td><td>GROMACS</td><td>(Kulasinski, Guyer, Derome, & Carmeliet, 2015a)<br>(Hanus & Mazeau, 2006)</td></tr><tr><td>Crystalline cellulose, Hemicelluloses</td><td>PCFF</td><td>–</td><td>(Pereira, Silveira, Dupree, & Skaf, 2017; Sarkar et al., 2023)</td></tr><tr><td>Crystalline cellulose, Hemicelluloses</td><td>CHARMM36</td><td>NAMD</td><td>(Addison et al., 2024; Cresswell et al., 2021)</td></tr><tr><td>Crystalline cellulose, Hemicelluloses, Lignin</td><td>CHARMM36</td><td>NAMD</td><td>(Charlier & Mazeau, 2012)<br>(Hao, Ho Tam, Lu, & Lau, 2018)</td></tr><tr><td>Crystalline celluloses, Hemicelluloses, Lignin</td><td>GROMOS 53a6</td><td>GROMACS</td><td>(Kulasinski, Derome, & Carmeliet, 2017)</td></tr><tr><td>Crystalline cellulose, hemicelluloses, lignin</td><td>COMPASS</td><td>LAMMPS</td><td>(Aguilera-Segura, Di Renzo, & Mineva, 2020)</td></tr><tr><td>Crystalline cellulose, Hemicelluloses, Lignin</td><td>GROMOS 53a6</td><td>GROMACS</td><td>(Lindner, Petridis, Schulz, & Smith, 2013; Zhang et al., 2022)</td></tr><tr><td>Crystalline cellulose, Lignin</td><td>CHARMM36</td><td>GROMACS</td><td>(Bregado, Tavares, Secchi, & Segtovich, 2020)</td></tr><tr><td>Crystalline cellulose, Lignin</td><td>CHARMM</td><td>GROMACS</td><td>(Bregado, Tavares, Secchi, & Segtovich, 2020)<br>(Chen, Lickfield, & Yang, 2004)</td></tr><tr><td>Amorphous cellulose</td><td>PCFF, COMPASS</td><td>LAMMPS</td><td>(Zhu et al., 2016)</td></tr><tr><td>Amorphous cellulose</td><td>PCFF</td><td>–</td><td>(Ren et al., 2021)</td></tr><tr><td>Amorphous cellulose</td><td>PCFF</td><td>Materials Studio</td><td>(Bregado, Secchi, Tavares, de Sousa Rodrigues, & Gambetta, 2019)</td></tr><tr><td>Amorphous cellulose</td><td>CHARMM36</td><td>LAMMPS</td><td>(Kulasinski et al., 2014)<br>(Zhang, Tschopp, Horstemeyer, Shi, & Cao, 2013)</td></tr><tr><td>Amorphous cellulose, Hemicellulose</td><td>GROMOS 53a6</td><td>GROMACS</td><td>(Kulasinski, Guyer, Keten, Derome, & Carmeliet, 2015)</td></tr><tr><td>TEMPO-oxidized cellulose, Xyloglucan</td><td>GLYCAM06</td><td>GROMACS</td><td>(Kishani, Benselfelt)</td></tr></tbody></table>

(continued on next page)
<table><caption>**Table 1 (continued)**</caption><thead><tr><th>Constituent type</th><th>Force field</th><th>Package</th><th>Reference</th></tr></thead><tbody><tr><td>TEMPO-oxidized cellulose</td><td>OPLS-AA</td><td>GROMACS</td><td>Wågberg, & Wohlert, 2021)<br>(Paajanan et al., 2016)</td></tr><tr><td>TEMPO-oxidized cellulose</td><td>INTERFACE</td><td>-</td><td>(Asgarpour Khansary et al., 2020)<br>(Houtman & Atalla, 1995)</td></tr><tr><td>Cellulose surface</td><td>CVFF</td><td>INSIGHTII/< DISCOVER</td><td>(Gupta, Dupree, Petridis, & Smith, 2023)</td></tr><tr><td>Cellulose surface</td><td>CHARMM</td><td>GROMACS</td><td>(Mudedla et al., 2021)</td></tr><tr><td>Defected cellulose</td><td>CHARMM36</td><td>GROMACS</td><td>(Heasman, Mehandzhiyski, Ghosh, & Zozoulenko, 2023)</td></tr><tr><td>Cellulose chains</td><td>OPLS-AA</td><td>GROMACS</td><td>(Paajanan & Vaari, 2017)<br>(Guo, Wang, & Jiang, 2023)</td></tr><tr><td>Cellulose chains</td><td>ReaxFF</td><td>LAMMPS</td><td>(Azimzadeh Irani, Askari, Jahanfar, Nasehi, & Hamedi, 2023)</td></tr><tr><td>Cellulose chain</td><td>GLYCAM04</td><td>AMBER</td><td>(French & Johnson, 2009)</td></tr><tr><td>Cellobiose, Xylobiose</td><td>CSFF<br>GAFF</td><td>CHARMM<br>AMBER8</td><td>(Peri, Nazmul Karim, & Khare, 2011)<br>(Mohan, Viswanath, Banerjee, & Goud, 2018)</td></tr><tr><td>Cellobiose, Xylose</td><td>CHARMM</td><td>NAMD</td><td>(Berglund et al., 2016)<br>(Martínez-Abad, Jiménez-Quero, & Wohlert, 2020)</td></tr><tr><td>Galactoglucomannan, Xyloglucan, xylan</td><td>GLYCAM06</td><td>GROMACS</td><td>(Martínez-Abad et al., 2017)<br>(Berglund et al., 2020)</td></tr><tr><td>Mannan</td><td>GLYCAM06</td><td>GROMACS</td><td>(Berglund et al., 2018)<br>(Zhang, Coasne, Guyer, Derome, & Carmeliet, 2020)</td></tr><tr><td>Lignin, Glucuronoarabinoxylan</td><td>COMPASS</td><td>-</td><td>(Youssefian, Jakes, & Rahbar, 2017)<br>(Liu, Ku, & Jin, 2022)</td></tr><tr><td>Lignin, Arabinoxylan</td><td>ReaxFF</td><td>LAMMPS</td><td>(Petridis, Schulz, & Smith, 2011; Rawal, Zahran, Dhital, Akbilgic, & Petridis, 2020)</td></tr><tr><td>Lignin</td><td>CHARMM<br>lignin</td><td>GROMACS</td><td>(Manna, Datta, & Ghosh, 2021; Petridis & Smith, 2016)<br>(Li, Xu, Jin, Luo, & Fan, 2019)</td></tr><tr><td>Lignin</td><td>ReaxFF</td><td>-</td><td>(Liu et al., 2020)</td></tr><tr><td>Lignin</td><td>ReaxFF<br>PCFF</td><td>LAMMPS</td><td>(Vu, Chaffee, & Yarovsky, 2002)</td></tr><tr><td>Lignin</td><td>COMPASS<br>DREIDING</td><td>DL_POLY</td><td>(Zhang & LeBoeuf, 2009)</td></tr></tbody></table>

$$E_{bonded} = k(d - d_0)^2 \qquad (1)$$

where $E_{bonded}$ is the potential energy (stretching energy) of a covalent bond, $k$ is the spring's constant, $d_0$ is the equilibrium distance between

the two atoms making bond, and $d$ is the current distance at which a bond is compressed or stretched. To bring some examples of such functional forms defining force fields, it can be mentioned:

$$E(d_N) = \sum_{i \in \text{bonds}} k_{b,i} (d_i - d_{0,i})^2 + \sum_{i \in \text{angles}} k_{a,i} (\Theta_i - \Theta_{0,i})^2 \\ + \sum_{i \in \text{torsions}} \sum_{n} \frac{1}{2} V_i^n [1 + \cos(n\Omega_i - \Gamma_i)] \\ + \sum_{j=1}^{N-1} \sum_{i=j+1}^{N} f_{ij} \left\{ \epsilon_{ij} \left[ \left( \frac{d_{ij}^0}{d_{ij}} \right)^{12} - 2 \left( \frac{d_{ij}^0}{d_{ij}} \right)^6 \right] + \frac{q_i q_j}{4\pi\epsilon_0 d_{ij}} \right\}, \qquad (2)$$

which shows an expression for the AMBER general force fields (Wang, Wolf, Caldwell, Kollman, & Case, 2004). In this example, the first term computes the bonded stretching energies (as an ideal spring), the second term represents the energy of the angles between sets of three atoms, the third term represents the torsional energies defined by two planes, and the last term represents the non-bonded energies calculated for all the atom pairs including 12–6 Lennard-Jones and the coulombic potentials. These expressions mainly include constants (here $k, n, \Gamma, d, q, \epsilon_0$) which are defined based on the quantum calculations and are provided by the developers of the specific force field. While many conventional force fields use the same principle, the differences are mainly in the defined constants, as well as the equilibrium distances, angles, etc. For instance, CHARMM force field (MacKerell et al., 1998) is defined as:

$$E(d) = \sum_{i \in \text{bonds}} k_{b,i} (d_i - d_{0,i})^2 + \sum_{i \in \text{angles}} k_{a,i} (\Theta_i - \Theta_{0,i})^2 \\ + \sum_{i \in \text{dihedrals}} \frac{1}{2} V_i^n [1 + \cos(n\Phi_i - \delta_i)] \\ + \sum_{i \in \text{improvers}} k_{\Omega_{m,i}} (\Omega_i - \Omega_{0,i})^2 + \sum_{\text{Urey-Bradley}} k_{\text{Uuk,i}} (u_i - u_{0,i})^2 \\ + \sum_{\text{nonbonded}} \left( \epsilon \left[ \left( \frac{R_{\text{minij}}}{d_{ij}} \right)^{12} - \left( \frac{R_{\text{ij}}}{d_{ij}} \right)^6 \right] + \frac{q_i q_j}{\epsilon d_{ij}} \right). \qquad (3)$$

According to what CHARMM is parametrized by, torsional angles are defined as two terms, i.e., dihedrals and improper dihedrals. Additionally, the Urey-Bradley term is included in this potential accounting for the angle bending energies through 1,3 non-bonded interactions. Potentials can become even more complicated when computing different energies, such as dihedrals. For instance, OPLS (Jorgensen, Maxwell, & Tirado-Rives, 1996) expresses the interatomic energies as:

$$E(d) = \sum_{i \in \text{bonds}} k_{b,i} (d_i - d_{0,i})^2 + \sum_{i \in \text{angles}} k_{a,i} (\Theta_i - \Theta_{0,i})^2 \\ + \sum_{i \in \text{dihedrals}} \left( \frac{1}{2} V_1^n [1 + \cos(\phi_1 - \phi_2)] \right. \\ + \sum_{i \in \text{dihedrals}} \left. \frac{1}{2} V_2^n [1 + \cos(2\phi_1 - \phi_2)] \right. \\ + \sum_{i \in \text{dihedrals}} \left. \frac{1}{2} V_3^n [1 + \cos(3\phi_1 - \phi_3)] \right) \\ + \sum_{i \in \text{dihedrals}} \left. \frac{1}{2} V_4^n [1 + \cos(4\phi_1 - \phi_4)] \right) \\ + \sum_{\text{nonbonded}, i \to j} f_{ij} \left( \frac{A_{ij}}{r_{ij}^6} + \frac{C_{ij}}{r_{ij}^6} + \frac{q_i q_j e^2}{4\pi\epsilon_0 r_{ij}} \right). \qquad (4)$$

On the contrary, reactive force fields, e.g., ReaxFF (Zhang & Van Duin, 2018), usually use a bond order by which breakage of the bonds is defined; when the distance between two atoms is beyond a certain point bonds are considered broken. While the calculations performed are not simply defined in words, the system energies are calculated as:
$$E_{system} = E_{bond} + E_{lp} + E_{over} + E_{under} + E_{val} + E_{pen} + E_{coa} + E_{C2} + E_{triple} + E_{tors} + E_{conj} + E_{hb} + E_{vdw} + E_{Cb} \quad (1)$$

where $E_{bond}$ is the bonded energies calculated by the bond order, $E_{lp}$ is the lone pair energy, $E_{over}$ is the overcoordination energy of an atom for a broken-up lone electron pair, $E_{under}$ is the undercoordination energy, $E_{val}$ is the angle energy, $E_{pen}$ is called the penalty energy, $E_{coa}$ is the three-body conjugation term where the latter three construct the valance angle terms, $E_{C2}$ is the correction energy for a $C_2$ pair due to a very strong triple bond, $E_{triple}$ is the triple bond energy correction for a carbon monoxide, $E_{tors}$ is the torsion rotation barrier energy, $E_{conj}$ is the four-body conjugation energy term, $E_{hb}$ is the hydrogen bond energy, and $E_{vdw}$ and $E_{Cb}$ are the non-bonded Lennard-Jones and Coulombic terms (please refer to ReaxFF documentation for an explanation of each of these terms: https://www.scm.com/doc/ReaxFF/).

Conventional force fields are limited in simulating large systems over long timescales. While united atom force fields (providing coarse-grained-like potentials) offer faster computations compared to all-atom ones, still the urge for enhanced methods can be felt. In such force fields, groups of atoms are considered as one unit, for instance, the hydrogens on aliphatic carbons of amino acids, being united with the respective carbons (Yang et al., 2006). Coarse-grained (CG) models enable faster computations by grouping atoms into single particles. However, CG force fields are usually property/system-specific. Hence, there still seem to be huge possibilities in this field of research to expand the versatility of such parameter sets. Section 3.8 briefly covers some CG models for cellulose nanocrystals.

Development of carbohydrate force fields is an ongoing attempt for decades. Early force fields for carbohydrates include the study by Neyertz et al. (2000). While recently, Charvati, Zhao, Wu, and Sun (2021) parametrized an AMBER-based force field using quantum mechanics, reproducing density and lattice parameters of cellulose I/β.

As reliability of MD simulations mostly depend on the accuracy of the force fields, ongoing assessments and comparisons have been made thoroughly in the literature. Stortz, Johnson, French, and Csonka (2009) assessed 18 force fields for carbohydrates (54 conformers). While most of the force fields found proper minima, CHARMM force field gave odd dihedral angles for the minima. Overall, MM3, GLYCAM06, and GROMOS were reasonable choices in studying disaccharides.

In a later study, Foley, Tessier, and Woods (2012) reviewed force fields including AMBER, CHARMM, GLYCAM, GROMACS and MM3 for carbohydrates. A brief review of these force fields' history, challenges, and development strategies are brought together.

Comparison of the force fields has brought a lot of added value to the community of carbohydrates and modeling, signifying the main capabilities of each force field. Examples can be provided in how various force fields could provide different structural characteristics and properties on cellulose. Matthews et al. (2012) found GLYCAM06 and CHARMM35 provide conformations of I/β cellulose similar to those observed at 500 K, unlike GROMOS 45a4 reproducing significantly different structures. The same shift in the $γ$ value (also including deviations in a unit cell parameter) reported in studies performed by Bergenstråhle et al. (2007) for cellulose was reported here for GROMOS 45a4, as well as longer c parameter also modelled by GLYCAM06. Nevertheless, GLYCAM06 seems to best reproduce cellulose crystal structures, although some lattice parameter deviations exist (Zhang, Bulone, et al., 2011).

Modifying Lennard-Jones potentials and charges in GROMOS 53a6 and 56A$_{carbo}$ is shown to improve cellulose crystal predictions (Chen et al., 2014). For instance, while both the 53a6 and 56A$_{carbo}$ versions modelled the a axis parameter in cellulose I/β with an increase of 7% and 8.3%, respectively, estimations could be significantly improved when the charges on single hydrogen bonded carbons were adjusted according to those of the CHARMM35 force field.

Dri et al. (2015) compared three ReaxFF versions with two other force fields GLYCAM and COMPASS, and found no force field ideally

predicts all cellulose I/β properties. COMPASS best predicted lattice parameters. COMPASS was also previously assessed by Eichhorn and Davies (2006) in 2006, and the results of the their minimized cellulose structures showed an increase in the c-parameter for the I/β and II models, whereas this value was lower than the experimental observations for the I/α.

Another important factor in simulating carbohydrates and respective force fields, is the choice of water models. As different water molecules have been shown to exhibit different behavior in MD simulations and extensive care must be taken when water molecules, solvated structures, and force fields are chosen, it is necessary to conduct studies to benchmark the effect of water molecule models. The choice of water models can lead to variations in computational speeds (Sauter & Grafmüller, 2015). Nevertheless, recent developments in polarizable water models for atomistic simulations have demonstrating promising performances, which could change the picture in the new future (Lambros & Paesani, 2020; Xiong, Izadi, & Onufriev, 2022).

When selecting a force field for molecular dynamics simulations of carbohydrates and natural fibres, several important parameters should be considered. One must carefully assess force field accuracy, compatibility, parameterization, electrostatics, and water interactions when selecting for carbohydrate/fibre simulations. In short, target property determines best choice. For instance, it is compulsory to consider ahead, where the usage of a non-reactive or reactive force field meets the expectations of the study. As a more detailed example, while non-reactive studies would provide insights on mechanisms and elasticity, extraction of ultimate properties and reactions are only feasible through reactive molecular dynamics. It is not easy to draw a clear conclusion which force field fits best modeling carbohydrates, but literature offers three main force fields, i.e., GLYCAM, CHARMM, and GROMOS. Examples provided above and many of the current study motivates us to choose GLYCAM06 as one of the most accurate force fields in modeling cellulose and hemicelluloses. Yet for cellulose, slight structural inaccuracies could be seen by GLYCAM. The literature offers an extended list of research performed on modeling plant cell wall constituents aiming at various characteristics. Here, we try to elaborate on the weaknesses and strengths of such simulation studies.

### 3. Cellulose

Cellulose, being the most abundant polymer on earth, has been the focus of many studies in recent decades (Richely, Bourmaud, Placet, Guessasma, & Beaugrand, 2022; Zhang et al., 2019) due to its superior mechanical properties (Jin, Qin, & Buehler, 2015), interesting hygrothermal behavior (Pietrucci, Boero, & Andreoni, 2021; Zhang et al., 2021), and sustainability motivating studying it from many experimental and computational points of view (Eichhorn et al., 2010; Habibi, Lucia, & Rojas, 2010; Moon, Martini, Nairn, Simonsen, & Youngblood, 2011). Researchers have debated cellulose's structure regarding fibril shape and size, arrangement of chins, and their placement in the plant cell wall (Cosgrove, 2005, 2014; Cosgrove & Jarvis, 2012). Molecular dynamics simulations have elucidated some not easily experimentally observable aspects. This section would first aim at introducing cellulose in as a whole, its crystal types, and structural characteristics. We will then introduce several challenges in cellulose community targeted by molecular dynamics simulations and their results in further sub-chapters.

Cellulose chains consist of anhydroglucose units (AGU, and hereafter referred to as glucose unit, GU) connected via β-1,4 glycosidic linkages. Assemblies of such chains construct the basic building block, the cellulose microfibril (CMF), often also called elementary fibril (Payne et al., 2015). The repeating unit of cellulose is occasionally referred to as cellobiose, Fig. 1a, however, studies offer the importance of anhydroglucose as the sole repeating unit of cellulose and not cellobiose (French, 2017). Each atom in cellulose is distinguished by a number, and main angles and dihedrals ($ω, φ$, and $ψ$) generally show specific
distributions in different cellulose polymorphs (see Fig. 1b).

Cellulose commonly occurs in different crystalline polymorphs I(α, β) (Nishiyama, Langan, & Chanzy, 2002; Nishiyama, Sugiyama, Chanzy, & Langan, 2003), II (Langan, Nishiyama, & Chanzy, 1999, 2001; Langan, Sukumar, Nishiyama, & Chanzy, 2005; Leng et al., 2021; Miyamoto et al., 2009), and III (Ford, Stevens, Johnson, & French, 2005; Wada, Chanzy, Nishiyama, & Langan, 2004; Yui & Hayashi, 2009). Cellulose IV (Hermans & Weidinger, 1946) is generally known as a disordered form of cellulose I/β and not a main polymorph (Wada, Heux, & Sugiyama, 2004). Cellulose chains in Iα form a triclinic lattice, while their arrangement has a monoclinic lattice structure in I/β. Fig. 1c-f shows Iα and I/β unit cells.

Cellulose II is a regenerated crystal form of cellulose with a monoclinic unit cell with P2₁ symmetry. Cellulose III can be derived from the treatment of cellulose I and II with liquid ammonia, having a unit cell with P2₁ symmetry as well (Maurer et al., 2013). Additionally, cellulose IV can be produced as a result of cellulose II at high temperatures (Hermans & Weidinger, 1946). Comparisons of cellulose I, II, and III structures are depicted in Fig. 2.

Such cellulose arrangements induce different positioning of atoms with respect to each other. This will, hence, cause alternative hydrogen bonding patterns (more on this in Section 3.2). The rotation of the hydroxymethyl groups on cellulose chains facilitates different hydrogen

bonding patterns between the sheets and within the sheets (inter- and intra-chain hydrogen bonds). These conformations of the hydroxymethyl groups can be in three main forms: tg, gg, and gt as shown in Fig. 3. As an example, in cellulose I, the orientation of the hydroxymethyl groups and the respective torsional angle ω is in a tg conformation. These torsional conformations can be indicative of the crystallinity of cellulose structures, as will be discussed later in Section 3.4.

MD has aided research to inspect interactions of various substances with cellulose with a detailed atomistic vision, providing illustrative views which might not be possible to achieve experimentally. Cellulose modeling dates back to the 90's, when its crystalline structure (Iα and I/β) were studied (Heiner et al., 1995). A 1996 study found I/β models more stable than Iα (exhibiting structures close to those at high temperature), which showed more fluctuations that could induce transitions to I/β through hydroxymethyl rotations and slip planes in the c direction (Hardy & Sarko, 1996).

MD provides valuable atomistic insights into cellulose interactions and behavior, often difficult to obtain experimentally. With enhanced computational power and methods, MD will likely continue advancing cellulose characterization. In the following sections, we review the findings of MD works on different aspects of cellulose.

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** a) One form of repeating unit in cellulose chains, a cellobiose. The reducing end would connect to the right side and the non-reducing end on the left side of the cellobiose. The degree of polymerization of the chain will, hence, become 2N + 2, where the addition of 2 refers to 2 GUs on each end, i.e. the reducing and non-reducing ones. b) Representation of part of a cellulose chain and respective atom nomenclature. Cellulose generally shows a twofold (2₁) helical screw, i.e. there is a 360° rotation of the monomers by every cellobiose (the third monomer has rotational identity to the first monomer). Hydrogen bonds are shown in black dotted lines. The main dihedral angles about the glycosidic linkage are ψ (C1-O4-C4-C5) and φ (O5-C1-O4-C4), and the hydroxymethyl dihedral angle is ω (O5-C5-C6-O6). Carbon atoms are shown in cyan, oxygen in red, and hydrogen in white spheres. Figure adapted with permission from (Khodayari, Van Vuure, et al., 2020). c) A 3D comparison of the Iα and I/β unit cells. γ represents the angle between a and b lattice vectors. Figure adapted and redrawn with permission from (Sugiyama, Vuong, & Chanzy, 1991). d) and e) depict the hydrogen bonding sheet displacement for Iα and I/β lattices, respectively. f) Representation of the chains along the main axis (c) and the a- and b-parameters. Views 1, 2, and 3 refer to the (010), (110), and (100) planes for Iα and (110), (200), and (1-10) for I/β, respectively. View 4 represents the (1-10) and (010) planes for Iα and I/β lattices, respectively. Unit cells are represented with dashed lines for Iα and with solid lines for β. Figures d-f are adapted and redrawn with permission from (Imai, Putaux, & Sugiyama, 2003). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 2. Representations of unit cells of a) cellulose I/β, b) cellulose II, c) cellulose III. Cases for cellulose I/β and III are arbitrarily chosen to represent their category.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. Conformations of the hydroxymethyl group forming three distinct χ (dihedral of C4-C5-C6-O6) values of -60° (tg), 60° (gg), and 180° (gt). H61 and H62 are two hydrogen bonds connected to the C6, and H5 is the hydrogen bond on the C5 atom. The center of the circles is where the C5 and C6 atoms are positioned. Figure drawn with inspirations from (Maurer et al., 2013).</figcaption></figure>

### 3.1. Cross-section and chain number in plant microfibrils

A continuing source of debate in the cellulose and plant physiology community is how many cellulose chains build up an individual plant CMF (Cosgrove, 2014). While many studies proposed a 36-chain structure is most probable (Elazzouzi-hafraoui et al., 2008; Endler & Persson, 2011; Mutwil, Debolt, & Persson, 2008), recent studies on cellulose synthase complexes and transmission electron microscopy assert 18- or 24-chain models likely fit within six lobes best (Cosgrove, 2022; Fernandes et al., 2011; Kubicki et al., 2018; Nixon et al., 2016; Rosén et al., 2020; Vandavasi et al., 2016; Willhammar et al., 2021; Zhong, Cui, & Ye, 2019). In our opinion, at current knowledge level, we believe that the 18 chain model is the most likely model for the CMF as it can also explain in large parts results on the 36-chain model, which can be interpreted as a dimer of the 18-chain model. Apart from discussions on the reliable chain number in plant CMFs, there are also open questions related to the chain arrangement, resulting in different cross sections of the plant CMF (Hill, Hammudi, & Tien, 2014; Wang & Hong, 2016). Consequently,

cellulose models used currently are selected based on exposed faces and desired investigations. Examples of some possible cross-sections used in the literature are depicted in Fig. 4.

There are more CMF modeling computational tools available now than 20 years ago. Cellulose-builder generates Iα, Iβ, II, and III crystal-line structures of any size, with/without periodic boundary conditions (PBC) (Gomes & Skaf, 2012), helping produce structures reported by Nishiyama et al. (2002) and Nishiyama et al. (2003). Different cellulose face planes can also be efficiently modelled.

Of the earliest set of studies, the works of Heiner et al. (1995) on 18-chain models of Iα and Iβ could be mentioned. 36-chain structures also seem to be popular from the 90's (Heiner et al., 1998; Heiner & Teleman, 1997). Minor deviations from the lattice parameters, e.g. increased *c* parameter, have been widely reported with many force fields since the early years (Baird, O'Sullivan, & Banks, 1998).

To conclude on probable CMF conformations, Oehme, Downton, et al. (2015) modelled 36-, 24-, 18-chain I/β cellulose, arguing 18-, 24-chain models are more viable than the 36-chain, although deviations
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Some popular cross-sections suggested for non-relaxed untwisted CNCs. a) 18-chain (34443), b) 18-chain (234432), c) 36-chain (hexagonal), and d) 36-chain (square) models. The numbers in the parenthesis refer to the number of chains in each plane. The numbers in (c) are the miller indices of each cellulose (sur)face.</figcaption></figure>

were noticeably greater than the 36-chain from initial conformations and from any phase identified experimentally.

Later, comparisons between 18-, 24-, 36-chain Iβ and Iα cellulose, showed that the 36-chain model lattice parameters matched experiments, while smaller ones exhibited high-temperature-like structures (Oehme et al., 2018). Smaller models also had more non-tg hydroxymethyl groups than the 36-chain model. However, the model with non-tg conformations better fitted the in silico NMR spectra peak sizes, especially for the Iα allomorph. Experimental evidence obtained through ssNMR generally indicates differences in the hydroxymethyl conformations on the surface compared to those in the inner regions (Phyo, Wang, Yang, O'Neill, & Hong, 2018).

Comparisons of 36- and 32-chain Iβ cellulose resulted in unit cell angle changes nearing 90° and tilts in every second layer (hydrogen bonding plane in Fig. 1f) in both models (Matthews et al., 2006).

In a study by Nishiyama et al. (2012), models were developed for 4–169 chain Iβ cellulose to investigate X-ray diffraction (XRD) patterns with and without water, examining the impact of size and disorder. The investigation focused on three sets of structures: 1) structures without solvation, 2) structures with one layer of solvation shell, and 3) structures with two layers of solvation shell. Across all models, the 110 peak was influenced by small oscillations in the 2θ range of 10°-15°, possibly connected to a slight deterioration in the crystallinity of the models.

Low-angle scatterings were found to be at their minimum when considering one layer of solvation shell. These scatterings diminished under two conditions: the presence of water and certain structural deviations from the original coordinates. Notably, the effect of water was primarily observed on these low-angles, with no other pronounced effects detected.

Interestingly, at the 200 peak (around 2θ = 23°), a shift to higher values was observed when the number of chains increased from 36 to 144. This shift was associated with changes, such as shorter plane distances, in the intermolecular spacing (d<sub>200</sub>). The 004 peaks at 2θ = 34.5°

slightly shifted to the left with an increase in the number of chains, suggesting that chains were somewhat stretched in structures with higher chain numbers than in smaller ones.<fnref n="2" /> Similar observations were made in 36-chain simulations with GLYCAM06 (Khodayari, Hirn, et al., 2021), indicating that computational models tend to have slightly bigger c parameter than their experimental counterparts, especially for larger numbers of chains.

In summary, literature suggests that increasing the chain numbers could lead to more compact, yet stretched chains. Hence, such observations might give a clue about the number of the chains in CMFs. Lower numbers, e.g., 18, might have smaller c parameter (axial compactness), comparable with those observed in the experimental powder diffraction patterns. Since different cross-sections have varied hydrophilic/hydrophobic surface ratios with different water interactions, hydration research could help illustrating the true cross-section. In particular, water shows a two-layer shell on the surface of cellulose (more details in Section 3.7), one at about 2.2 Å where the density of water is 2.5 times the density of the bulk, and one at 5.5 Å where the density of the hydration shell reaches 1.3 times that of the bulk. The first shell above the hydrophobic faces is shown to be centered at about 3.6 Å (Matthews et al., 2006).

A summary of different works and the cellulose model used in the studies is brought in Table 2. The above studies, as well as data shown in Table 2, illustrate that there is still no consensus in the community on the accurate number of chains in the plant CMF. Atomistic studies, dating back to decades ago, and the current research have all covered a

<footnote n="2">² It must be noted that, often in the literature, the experimental lattice spacings are taken from large crystals. However, experimental values determined from smaller crystals (such as wood, pulp, or CNCs) can better suit comparisons with computational studies. Examples are larger d<sub>200</sub> lattice spacing for smaller crystals mentioned by (Nishiyama et al., 2012). Also see (Newman, 2008), (Huang et al., 2018), and (Zitting et al., 2023).</footnote>
<table><caption>**Table 2**
The CNC models used in the studies reported in the literature review and their respective cross-sections, degree of polymerization (DP).</caption><thead><tr><th>Cellulose type</th><th>Cross-section [# chains]</th><th>Length [DP]</th><th>Reference</th></tr></thead><tbody><tr><td>Iβ</td><td>7</td><td>8</td><td>(Aguilera-Segura et al., 2020)</td></tr><tr><td>Iβ</td><td>7</td><td>10</td><td>(Ishida, 2020)</td></tr><tr><td>Iβ</td><td>14, 18, 23, 28, 34, 40</td><td>PBC (34)</td><td>(Zitting et al., 2023)</td></tr><tr><td>Iβ</td><td>16</td><td>8</td><td>(Besombes & Mazeau, 2005a, 2005b)</td></tr><tr><td>Iβ</td><td>16 (square)</td><td>20</td><td>(Tolonen et al., 2015)</td></tr><tr><td>Iβ</td><td>16 (square)</td><td>20</td><td>(Djahedi et al., 2016)</td></tr><tr><td>Iβ</td><td>18</td><td>22</td><td>(Pereira et al., 2017)</td></tr><tr><td>Iβ</td><td>18</td><td>20</td><td>(Muthoka et al., 2018)</td></tr><tr><td>Iβ</td><td>18 (3 × 6), 16 (4 × 4)</td><td>PBC</td><td>(Diaz et al., 2013)</td></tr><tr><td>Iβ</td><td>18 (234432)</td><td>30, PBC (34)</td><td>(Zitting et al., 2021)</td></tr><tr><td>Iβ</td><td>18 (234432)</td><td>16</td><td>(Garg et al., 2020)</td></tr><tr><td>Iβ</td><td>18 (234432)</td><td>20</td><td>(Muthoka et al., 2022; Oehme, Doblin, et al., 2015)</td></tr><tr><td>Iβ</td><td>18 (234432)</td><td>28</td><td>(Zhang et al., 2022)</td></tr><tr><td>Iβ</td><td>18 (234432)</td><td>40</td><td>(Addison et al., 2024; Sarkar et al., 2023)</td></tr><tr><td>Iβ</td><td>18 (234432)</td><td>80</td><td>(Fujisawa et al., 2023)</td></tr><tr><td>Iβ</td><td>18 (234432)</td><td>PBC</td><td>(Heinonen et al., 2022; Jaafar et al., 2019; Lombardo et al., 2018; Paajanen et al., 2022)</td></tr><tr><td>Iβ</td><td>18 (234432)</td><td>–</td><td>(Cresswell et al., 2021)</td></tr><tr><td>Iβ</td><td>19</td><td>8</td><td>(Nishiyama et al., 2008)</td></tr><tr><td>Iβ</td><td>20</td><td>PBC</td><td>(He et al., 2023)</td></tr><tr><td>Iβ</td><td>24</td><td>20</td><td>(Busse-Wicher et al., 2014)</td></tr><tr><td>Iβ</td><td>24, 48</td><td>10, 20</td><td>(Yui et al., 2006)</td></tr><tr><td>Iβ</td><td>24</td><td>PBC</td><td>(Charlier & Mazeau, 2012)</td></tr><tr><td>Iβ</td><td>24 (345543)</td><td>40</td><td>(Paajanen et al., 2019)</td></tr><tr><td>Iβ</td><td>25</td><td>20</td><td>(Zhao et al., 2014)</td></tr><tr><td>Iβ</td><td>25</td><td>70</td><td>(Shrestha et al., 2019)</td></tr><tr><td>Iβ</td><td>25 (2344543)</td><td>8</td><td>(Hanus & Mazeau, 2006)</td></tr><tr><td>Iβ</td><td>30 (5 × 6)</td><td>24</td><td>(Mazumder & Zhang, 2023)</td></tr><tr><td>Iβ</td><td>32 (8 × 4)</td><td>10</td><td>(Kulasinski, Keten, Churakov, Derome, & Carmeliet, 2014)</td></tr><tr><td>Iβ</td><td>32 (8 × 4)</td><td>PBC</td><td>(Wu et al., 2014)</td></tr><tr><td>Iβ</td><td>32 (8 × 4)</td><td>8</td><td>(Bergenstråhle et al., 2007; Bergenstråhle et al., 2009; Lindh et al., 2016)</td></tr><tr><td>Iβ</td><td>34, 36</td><td>20</td><td>(Paavilainen et al., 2012)</td></tr><tr><td>Iβ</td><td>36, 18 (234432)</td><td>20</td><td>(Muthoka et al., 2020)</td></tr><tr><td>Iβ</td><td>36 (9 × 4)</td><td>PBC</td><td>(Wu et al., 2013)</td></tr><tr><td>Iβ</td><td>36</td><td>40</td><td>(Zhang et al., 2015)</td></tr><tr><td>Iβ</td><td>36 (square)</td><td>6</td><td>(Heiner et al., 1998; Heiner & Teleman, 1997)</td></tr><tr><td>Iβ</td><td>36 (square)</td><td>8</td><td>(Jiang et al., 2019; Zhang, Bulone, et al., 2011)</td></tr><tr><td>Iβ</td><td>36 (square)</td><td>10</td><td>(Chen, Lo Re, et al., 2020; Kong et al., 2022; Kulasinski et al., 2017; Kulasinski, Guyer, Derome, & Carmeliet, 2015a)</td></tr><tr><td>Iβ</td><td>36 (square)</td><td>10, PBC</td><td>(Wohlert et al., 2023)</td></tr><tr><td>Iβ</td><td>36 (square), 32 (8 × 4)</td><td>14</td><td>(Matthews et al., 2006)</td></tr><tr><td>Iβ</td><td>36 (square), 24 (345543), 18 (234432)</td><td>16</td><td>(Pang et al., 2022)</td></tr><tr><td>Iβ</td><td>36 (square)</td><td>30</td><td>(Khodayari, Hirn, et al., 2021; Khodayari, Thielemans, et al., 2021; Khodayari, Van Vuure, et al., 2020)</td></tr><tr><td>Iβ</td><td>36 (square)</td><td>40</td><td>(Matthews et al., 2011; Ogawa et al., 2020)</td></tr><tr><td>Iβ</td><td>36 (square), 59</td><td>20, 40</td><td>(Paavilainen et al., 2011)</td></tr><tr><td>Iβ</td><td>36 (square)</td><td>40</td><td>(Chen, Wohlert, et al., 2022; Kong et al., 2023)</td></tr><tr><td>Iβ</td><td>36 (square)</td><td>80</td><td>(O'Neill et al., 2017)</td></tr><tr><td>Iβ</td><td>36 (square)</td><td>100</td><td>(Ciesielski et al., 2019)</td></tr><tr><td>Iβ</td><td>36 (square)</td><td>119</td><td>(Khodayari, Hirn, et al., 2020)</td></tr><tr><td>Iβ</td><td>36 (hexagonal)</td><td>8</td><td>(Bering et al., 2022)</td></tr><tr><td>Iβ</td><td>36 (hexagonal)</td><td>10</td><td>(Bregado et al., 2021; Zhang, Keten, et al., 2021)</td></tr><tr><td>Iβ</td><td>36 (hexagonal)</td><td>20</td><td>(Gupta et al., 2021; Sáenz Ezquerro et al., 2019)</td></tr><tr><td>Iβ</td><td>36 (hexagonal)</td><td>40</td><td>(Thu et al., 2022)</td></tr><tr><td>Iβ</td><td>36 (hexagonal)</td><td>–</td><td>(Geng et al., 2023; Hou et al., 2021)</td></tr><tr><td>Iβ</td><td>36 (hexagonal)</td><td>80</td><td>(Chen et al., 2023)</td></tr><tr><td>Iβ</td><td>36 (hexagonal)</td><td>61, 160</td><td>(Lindner et al., 2013)</td></tr><tr><td>Iβ</td><td>36 (hexagonal)</td><td>20, 30, 40, 60, 80</td><td>(Kannam et al., 2017)</td></tr><tr><td>Iβ</td><td>36 (hexagonal), 18 (3 × 6)</td><td>24</td><td>(Yu et al., 2018)</td></tr><tr><td>Iβ</td><td>36 (hexagonal), 24 (3 × 6)</td><td>22</td><td>(Busse-Wicher et al., 2016)</td></tr><tr><td>Iβ</td><td>36 (hexagonal), 24 (345543), 18 (234432)</td><td>PBC</td><td>(Oehme, Downton, et al., 2015)</td></tr><tr><td>Iβ</td><td>39</td><td>12</td><td>(Sánchez-Badillo et al., 2021)</td></tr><tr><td>Iβ</td><td>41</td><td>10</td><td>(Izumi et al., 2023)</td></tr><tr><td>Iβ</td><td>64 (square)</td><td>8</td><td>(Mazeau & Charlier, 2012)</td></tr><tr><td>Iβ</td><td>69 (hexagonal)</td><td>40</td><td>(Nascimento et al., 2022)</td></tr><tr><td>Iβ</td><td>81 (square)</td><td>20</td><td>(Hadden et al., 2014)</td></tr><tr><td>Iβ</td><td>9–289 (square)</td><td>20</td><td>(Hadden et al., 2013)</td></tr><tr><td>Iβ</td><td>1–100 (square)</td><td>20–400</td><td>(Bu et al., 2015)</td></tr><tr><td>Iβ</td><td>4–169</td><td>20</td><td>(Nishiyama et al., 2012)</td></tr><tr><td>Iβ</td><td>16, 25, 36, 144 (square)</td><td>PBC</td><td>(Wohlert et al., 2012)</td></tr><tr><td>Iβ</td><td>14 × 4</td><td>14</td><td>(Zhang, Brumer, et al., 2011)</td></tr><tr><td>Iβ</td><td>surface</td><td>12</td><td>(Sridhar et al., 2023)</td></tr></tbody></table>

(continued on next page)
<table><caption>**Table 2 (continued)**</caption><thead><tr><th>Cellulose type</th><th>Cross-section [# chains]</th><th>Length [DP]</th><th>Reference</th></tr></thead><tbody><tr><td>I/β</td><td>surface</td><td>20</td><td>(Mazeau & Rivet, 2008)</td></tr><tr><td>I/β</td><td>surface</td><td>PBC</td><td>(Karna et al., 2021; Ma et al., 2021; Malaspina & Faraudo, 2019; Mazeau & Vergelati, 2002; Mudedla et al., 2021; Xia et al., 2018; Zhu et al., 2020)</td></tr><tr><td>I/β</td><td>–</td><td>20</td><td>(Kumar et al., 2018)</td></tr><tr><td>I/β</td><td>–</td><td>PBC</td><td>(Chen et al., 2012; Deng et al., 2020; Dri et al., 2015; Hao et al., 2018; Molnár et al., 2018)</td></tr><tr><td>I/β, Iα</td><td>16, 18</td><td>4</td><td>(Mazeau & Heux, 2003)</td></tr><tr><td>I/β, Iα</td><td>18</td><td>PBC</td><td>(Heiner et al., 1995)</td></tr><tr><td>I/β, Iα</td><td>25</td><td>–</td><td>(Hardy & Sarko, 1996)</td></tr><tr><td>I/β, II, III<sub>1</sub></td><td>36 (square)</td><td>20</td><td>(Pang et al., 2023)</td></tr><tr><td>I/β, Iα</td><td>36 (hexagonal) 24 (345543) 18 (234432) 18 (34443) 18 (6 × 3)</td><td>20</td><td>(Oehme et al., 2018)</td></tr><tr><td>I/β, Iα, II</td><td>surface</td><td>PBC</td><td>(Maurer et al., 2013)</td></tr><tr><td>I/β, Iα, III</td><td>48</td><td>PBC</td><td>(Yui & Hayashi, 2007)</td></tr><tr><td>I/β, Iα, II, III</td><td>surfaces</td><td>PBC</td><td>(Trentin et al., 2021)</td></tr><tr><td>I/β, Iα, II, III</td><td>8 × 8, 6 × 6, 8 × 5, 9 × 5</td><td>8</td><td>(Chen et al., 2014)</td></tr><tr><td>I</td><td>surface</td><td>PBC</td><td>(Houtman & Atalla, 1995)</td></tr><tr><td>I, III</td><td>9 (3 × 3)</td><td>PBC</td><td>(Wada et al., 2011)</td></tr><tr><td>I, II</td><td>chain</td><td>14</td><td>(Heasman et al., 2023)</td></tr><tr><td>II</td><td>16 (4 × 4)</td><td>8</td><td>(Ito & Matsumoto, 2022)</td></tr><tr><td>II</td><td>61</td><td>20</td><td>(Chen, Jiang, et al., 2020)</td></tr><tr><td>II</td><td>–</td><td>4</td><td>(Eichhorn et al., 2005)</td></tr><tr><td>TEMPO-oxidized</td><td>36 (square)</td><td>20</td><td>(Paajanen et al., 2016)</td></tr><tr><td>TEMPO-oxidized<br>Carboxymethylated</td><td>16 (4 × 4)</td><td>10</td><td>(Garg et al., 2021)</td></tr><tr><td>TEMPO-oxidized</td><td>surface</td><td>PBC</td><td>(Kishani et al., 2021)</td></tr><tr><td>TEMPO-oxidized</td><td>single chain</td><td>–</td><td>(Asgarpour Khansary et al., 2020)</td></tr><tr><td>chain</td><td>–</td><td>2, 3, 4</td><td>(French & Johnson, 2009)</td></tr><tr><td>chain</td><td>–</td><td>8, 16, 32, 64</td><td>(Paajanen & Vaari, 2017)</td></tr></tbody></table>

wide range of cross-sections for cellulose. While the 18-, 24-, and 36-chain models seem to be the most accepted among all, biochemical and structural biology data collectively support the biosynthesis of 18-chain elementary microfibrils (Hill et al., 2014; Purushotham, Ho, & Zimmer, 2020). Yet, it must be noted that the perspectives acquired from cellulose biosynthesis may not directly correlate with the experimentally observed final and averaged structures, e.g. the 36-chain models. Further experimental strategies should target exploring this important characteristic of cellulose. We believe that incorporating experimental/computational studies on sorption properties, employing imaging techniques, and surface modification of cellulose could eventually relate findings with those resulted from the cellulose synthesis complex.

## 3.2. Hydrogen bonds

Hydrogen bonding has shown to be an important parameter in defining the crystallinity and structure of cellulose (Jarvis, 2023). The majority of the literature studies have focused on this aspect, as transformations between cellulose types stem from changes in hydrogen bonds. (H-bond). For instance, MD is capable of showcasing how the conversion of cellulose I to cellulose III can take place by this rearrangement (Wada et al., 2011).

In 2002, a fundamental study using diffraction techniques providing key insights on I/β cellulose structure was performed (Nishiyama et al., 2002). In particular, cellulose I/β exhibits two major intra-chain H-bonds (O2-HO2...O6 and O3-HO3...O5) and two inter-chain H-bonds (O6-HO6...O2 and O6-HO6...O3), depicted in Fig. 5.

Later in 2008, Nishiyama et al. (2008) inspected two possible hydrogen bonding schemes. Scheme A has O2 acting as an acceptor for interchain O6-HO6...O2 bonds and a donor for the O2-HO2...O6 intrachain bond throughout the fibril (see Fig. 6). Scheme B conversely has O2 as a donor in the interchain bond and acceptor in the intrachain bond. While these two schemes are mutually available, simulations showed chains in I/β nanocrystals arrange in scheme A bonds at room temperature and below. Disordered regions like crystal surfaces can exhibit scheme B patterns. In particular, a disrupted scheme A can

facilitate hydroxyl group orientations leading to scheme B H-bonds.

In 2016, Djahedi et al. (2016) presented findings on the role of H-bonds in the tensile behavior of CMFs - referred to as the “leverage effect”, first described by Altaner, Thomas, Fernandes, and Jarvis (2014). The intrachain O2-HO2...O6 distance remains constant when a cellulose chain is stretched while the O3-HO3...O5 distance lengthens (see Fig. 7a). The main question here is whether H-bonds are weak interactions, and whether they really contribute to 20–60 % of tensile modulus of CMFs as reported by Tashiro and Kobayashi (1991). Using potential energy calculations, Djahedi et al. (2016) reported that H-bonds contribute only 11 % to the chain's stiffness. The leverage effect increases this contribution; without it, H-bonds would only contribute 3 % to elasticity. Assuming pyranose ring infinite rigidity (which is probably not the case in real samples), this contribution was about 29 %, aligning with values later reported by Khodayari, Van Vuure, et al. (2020).

While hydrogen bonding defines crystallinity and contributes to cellulose's elastic modulus, Wohlert et al. (2021) concluded its role is less influential in other properties. H-bonds likely play a negligible role in solubility, polymer adsorption, and paper strength. However, they indirectly impact fibril aggregation and directly affect molecular deformation. Since H-bonds are short-range, longer-range interactions like van der Waals may dominate in inter-fibril binding. Hence their role depends on conditions.

However, a recent study on I/β cellulose reported that H-bonds govern interface shear interactions regardless of moisture or fibril alignment (Zhang, Keten, et al., 2021). The periodic make-and-break of H-bonds causes a “stick-slip” frictional behavior with a periodicity of 1.06 nm (approximately, length of a cellobiose), coherent with what also observed by Wei, Sinko, Keten, and Luijten (2018), Khodayari, Thielemans, et al. (2021), and He et al. (2023). Neither moisture nor misalignments altered the fundamental role of H-bonds in controlling this mechanism, being more pronounced for hydrophilic-hydrophilic surfaces than the hydrophilic-hydrophobic ones.

Multiple works have highlighted hydrogen bonding's role in the intrinsic right-hand twist of CMFs. Comparisons between cellulose fibrils
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Possible inter- and intra-chain H-bonds in cellulose *I*β, shown as thick and thin dashed lines, respectively, for the origin (a and b) and center (c and d) sheets. Adapted and redrawn with permission from (Nishiyama et al., 2002).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** The schematics for two H-bond networks A and B in cellulose *I*β. Arrows describe the donor-hydrogen-acceptor directions. Figure adapted and redrawn with permission from (Nishiyama et al., 2008).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** a) The leverage effect observed in cellulose crystals under tensile load, explained by (Djahedi et al., 2016). Each part of the chain acts as a separate spring, i.e. the glycosidic linkage (in blue), the ring (in green), and HO3...O5 bond (in red), controlling the tensile behavior of the chains. The intrachain HO2...O6 bond maintains its distance (anchor point in the figure), while the distance in HO3...O5 increases. b) Twist dependence on the cross-section size of cellulose, with and without water. Cellulose structures show lower twist in vacuum compared to the solvated states. The green circles show the situation when the charges of HO2 and O6 are removed, hence no HO2...O6 H-bond. Figure adapted and redrawn with permission from (Bu et al., 2015). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

of various cross-section and length, determined that twist is dictated by intrachain H-bonds – O2-HO2...O6 induces twist while O3-HO3...O5 restricts twist (Bu et al., 2015). Removing charges on HO2 and O6 to prevent this bond decreased twist in solvated fibrils. Larger cross-sections showed less twisting, as also observed by Hadden et al. (2013), exhibiting more pronounced effects in dry environments and for smaller cross-section structures (see Fig. 7b). However, similar analyses in another study reported neither intrachain H-bond exclusively governs twist (Kannam et al., 2017). Constraining angles/charges did not consistently increase or decrease twist in their models. In other words, the intrachain H-bonds do not drive or restrict twists. Additionally, a slight decrease in the twist was observed when longer models were

simulated, something that contradicts the results found by Bu et al. (2015).

Interestingly, hydrogen bonding patterns in the twisted fibrils (definite lengths) are shown to be equivalent to a periodically modelled cellulose crystal (infinitely long), in which twist does not exist. This confirms the validity of both models (periodic and non-periodic) and removes the stress from results obtained from fibrils/crystals with no proper or accepted twist value.

In short, owing to the extensive list of experimental and computational evidence, the hydrogen bonding pattern of cellulose is quite clear now. The versatile role of hydrogen bonding in cellulose emerges as a focal point in understanding its properties and behavior. As the

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** a) Twist of a 7-CNC bundle. Figure reprinted with permission from (Garg et al., 2020). b) Cross-section dependency of cellulose twist, modelled by Hadden et al. (2013). Larger cross-sections show a lower twist magnitude, while length does not impose a noticeable effect on the twist. Results are comparably larger than the experimentally measured values. Models shown by circles have a DP20. The model shown with a triangle has a longer DP106. The experimental estimations are taken from Hanley, Revol, Godbout, and Gray (1997). Figure adapted and redrawn with permission from (Hadden et al., 2013).</figcaption></figure>
literature unfolds from many discussions surrounding intrachain H-bonds to the diverse impacts on crystallinity, mechanical properties, and even hydrogen storage capacities, the ongoing research underscores the complexity of H-bonds in cellulose. While several mechanisms, such as the leverage effect during cellulose tension or the stick-slip behavior during cellulose shear, have been established for years now, researchers are still eagerly investigating H-bonds' role on some cellulose characteristics such as twist (see Section 3.3) or mechanical properties (see more in Section 3.5). The literature reports a wide range of 11–60 % for the contribution of H-bonds to mechanical properties, which appears to be valid only within its lower bound.

### 3.3. Cellulose twist

Twist of the cellulose crystals and microfibrils has been the center of attention for many years (Chen, Ogawa, Nishiyama, Bergenstråله-Wohlert, & Mazeau, 2015; Ogawa, 2021). CMFs have been shown to exhibit a right-hand twist when modelled in atomistic simulations. While there are several studies confirming the existence of the twist (e.g. by electron tomography, electron diffraction, and AFM), its origins are not fully understood. Drying is an important factor in inducing twist in cellulose materials on different scales and has been observed also for CMFs. The twist in CMFs, i.e. in their native environment is less studied but, also there, drying is considered to play a major role. A main issue is that resolution of electron tomography becomes an issue, making precise geometric statements challenging. Glucan chains (and all β-1,4 linked polysaccharides) develop a slight twist along their chain axes, i.e. not perfectly twofold conformations. Researchers believe that chain twist could be the source of twist in CMFs. Nevertheless, evidence suggests the nanoscale twist of CMFs causes the macroscopic twist of fibres observed experimentally (Garg et al., 2020).

Simulations of 7-CNC bundle (see Fig. 8a) demonstrated higher twist values for sulfate CNCs with Na⁺ counterions versus Ca²⁺ counterions, while both showing more twist than native cellulose (Garg et al., 2020). In all cases, the CNC bundle showed slightly less twist than individual CNCs, confirming macrofibrillar twist arises from twisting of individual CNCs. However, reported twist values for cellulose models (3–4.2°/nm) exceeded typical literature values (Shklyaev, Kubicki, Watts, & Crespi, 2014).

Other simulations mimicked CMF bundles, finding individual vacuum-equilibrated fibrils twisted 2.6–3.1°/nm versus 0.8°/nm for the bundle (Paajanen et al., 2019). Wet fibrils showed 4°/nm individually and 0.9°/nm for the bundle. Twist mismatches between bundles and individual fibrils was related to water penetrating interfaces, explaining experimental results (Hanley et al., 1997; Usov et al., 2015). Nevertheless, twist is generally observed disregarding water effects (Yui et al., 2006).

Although experiments and quantum mechanics (QM) calculations have targeted to clarify twist (Conley, Whitehead, & van de Ven, 2017; Hanley et al., 1997), the origins of CMF twist are still debated in literature. Different force fields have yielded varying twist values, decreasing over simulation time with GLYCAM06 and CHARMM35 but increasing with GROMOS 45a4 (Matthews et al., 2012). Nonetheless, a twofold screw axis (2₁) pseudo-symmetry is supported in the literature, where cellulose fibres can exhibit 2₁ symmetry without fully meeting the glycosidic torsion angle criteria $φ + ψ = -240°$ (French & Johnson, 2009).

For CMFs specifically, greater twist has been observed for surface chains versus interior chains, attributed to different H-bonding densities (Paavilainen et al., 2011). According to Paavilainen et al. (2011), simulation results showed that twist values were lower initially, converging to 2.3°/nm after 140 ns at 310 K. Increasing temperature to 350 K increased twist rates. A 36-chain model showed comparable twist around 0.8°/nm with what is usually obtained at room temperature, rising to 2.1–2.2°/nm at higher temperatures. The rapid initial twist appearance correlated with rapidly decreasing interchain H-bonds

except for intrasheet bonds. This synergistic H-bond change facilitates twist, challenging claims that twist is a force field artifact. Additional evidence like induced circular dichroism or electron tomography supports this (Conley et al., 2017; Majoinen et al., 2014). Results contradict findings that twist decreases around 500 K (Matthews et al., 2011). Comparisons of Iα and Iβ cellulose with cellulose III showed higher twisting in cellulose I than in cellulose III (Yui & Hayashi, 2007). Greater Iα twist was attributed to weaker intersheet interactions. But, temperature increased twist, also attributed to weaker intersheet interactions.

Computationally, twist alters scattering patterns like WAXS showing smaller smaller $d_{200}$ lattice spacing for twisted versus untwisted models (Zitting et al., 2023). Hadden et al. (2014) inspected twisting effects on cellulose diffraction data, modeling 81-chain Iβ cellulose fibrils as finite twisted and infinite periodic untwisted models. The twisted model showed a 1.17°/cellobiose twist rate. Both models indicated more compact structures (shorter distance between the 200 planes) than experimental data, based on 200 peaks shifting to larger 2θ values. Small-angle scatterings ($2θ = 10–13°$) were present except for non-periodic models, attributed to end distortions mimicking disordered regions. Trends were similar between models, indicating crystallographic discrepancies arise from unit cell changes rather than twist alone. As a conclusion, modest twist, observed in the modeling studies, agrees with the results from scattering techniques. DFT calculations and $^{13}$C NMR chemical shifts suggest experimental data fits twists between 2°/nm to -2°/nm (Shklyaev et al., 2014). Additionally, it must be noted that, isolated CMF twists may differ from surface-adhered CMFs (Kong et al., 2023).

CMF twist also depends on model size (Yui & Hayashi, 2007). For Iβ crystals, 9–289 chain models showed twists of ~9–0.35°/cellobiose, indicating cross-sectional twist dependency (Hadden et al., 2013), as also observed by Bu et al. (2015), with no effect from the degree of polymerization (see Fig. 8b). An 81-chain solvated model showed ~1°/cellobiose twist versus ~6°/cellobiose for vacuum, confirming water reduces twist. Twist appears to stem from vdW interactions maximizing crystalline packing efficiency, countered by interchain H-bonds and

<caption>**Table 3**</caption>

The computational twist angles observed in the literature based on the models and the degree of polymerization (Reported values are converted to °/nm to improve consistency, considering a cellobiose length of 10.38 Å).

<table><thead><tr><th>Model</th><th>Periodicity [DP]</th><th>Twist [°/nm]</th></tr></thead><tbody><tr><td>4-Chain (Bu et al., 2015)</td><td>40</td><td>19.3</td></tr><tr><td>4-Chain (Zhao et al., 2013)</td><td>20–40</td><td>9.9</td></tr><tr><td>9-Chain (Zhao et al., 2013)</td><td>20–40</td><td>5.8</td></tr><tr><td>9-Chain (Hadden et al., 2013)</td><td>20</td><td>8.7</td></tr><tr><td>9-Chain (Bu et al., 2015)</td><td>40</td><td>11.6</td></tr><tr><td>16-Chain (Bu et al., 2015)</td><td>40</td><td>4.8</td></tr><tr><td>16-Chain (Zhao et al., 2013)</td><td>20–40</td><td>3.1</td></tr><tr><td>18-Chain (Fujisawa et al., 2023)</td><td>80</td><td>4.7</td></tr><tr><td>18-Chain (Shklyaev et al., 2014)</td><td>–</td><td>0.2</td></tr><tr><td>24-Chain (Paajanen et al., 2019)</td><td>40</td><td>4</td></tr><tr><td>25-Chain (Hadden et al., 2013)</td><td>20</td><td>3.4</td></tr><tr><td>25-Chain (Bu et al., 2015)</td><td>40</td><td>2.9</td></tr><tr><td>36-Chain (Bu et al., 2015)</td><td>40</td><td>1.9</td></tr><tr><td>36-Chain (Kong et al., 2023)</td><td>40</td><td>1.8</td></tr><tr><td>36-Chain (Paavilainen et al., 2011)</td><td>20</td><td>0.8</td></tr><tr><td>36-Chain (Matthews et al., 2006)</td><td>14</td><td>1.4</td></tr><tr><td>36-Chain (Zhao et al., 2013)</td><td>20–40</td><td>1.3</td></tr><tr><td>36-Chain (Khodayari, Hirn, et al., 2021)</td><td>30</td><td>1.8</td></tr><tr><td>36-Chain (Gupta et al., 2022)</td><td>30</td><td>1.2</td></tr><tr><td>49-Chain (Hadden et al., 2013)</td><td>20</td><td>1.7</td></tr><tr><td>49-Chain (Bu et al., 2015)</td><td>40</td><td>1.5</td></tr><tr><td>59-Chain (Paavilainen et al., 2011)</td><td>40</td><td>0.6–0.8</td></tr><tr><td>64-Chain (Bu et al., 2015)</td><td>40</td><td>1.1</td></tr><tr><td>81-Chain (Hadden et al., 2014)</td><td>20</td><td>1.1</td></tr><tr><td>81-Chain (Hadden et al., 2013)</td><td>20</td><td>1.1</td></tr><tr><td>100-Chain (Bu et al., 2015)</td><td>40</td><td>1</td></tr><tr><td>289-Chain (Hadden et al., 2013)</td><td>20</td><td>0.3</td></tr><tr><td>4 CNCs (Paajanen et al., 2019)</td><td>40</td><td>0.8</td></tr><tr><td>7 CNCs (Garg et al., 2020)</td><td>40</td><td>3–4.2</td></tr></tbody></table>
solvent effects.

Table 3 provides a summary of reported degrees of twist for various models based on their length and number of chains. In general, it can be observed that the degree of twist is more closely tied to the number of chains (cross-section dependency) in the fibrils rather than the length of the models.

In conclusion, the investigation into the twist phenomenon in cellulose microfibrils (CMFs) reveals a complex interplay of factors influencing the observed twist. Atomistic simulations consistently demonstrate a right-hand twist in CMFs, but challenges arise in fully validating these observations experimentally due to scale differences in microfibril lengths. Additionally, as illustrated in Table 3, while summarizing measured twist angles for models of various sizes, force fields, and periodicities, discrepancies in the computed values persist. For instance, in a 36-chain structure, the values can range from 0.8 to 1.9°/nm, and for a 9-chain model, the values can range from 5.8 to 11.6°/nm. Hence, in reality, the decay in the twist versus size is not as smooth as depicted in Figs. 7b and 8b. Nevertheless, the study of bundles, such as the model of 7 CNCs presented by Garg et al. (2020), supports the idea that nanoscale twists contribute to the macroscopic twist observed in fibres. However, discrepancies in reported twist values, influenced by force fields and other simulation parameters, raise questions about the physical nature of CMF twist. While experiments and quantum mechanics contribute to the understanding, the debate continues regarding whether twist is an inherent characteristic of CMFs or a computational artifact. Notably, recent studies emphasize the influence of water, with simulations indicating changes in twist rates in the presence of water, further highlighting the need for a comprehensive understanding.

### 3.4. Cellulose at high temperatures

Cellulose is shown to be a resistant material to temperature (D'Acierno, Michal, & MacLachlan, 2023), with a glass transition temperature ($T_g$) of about 500 K, significantly higher than many other polymers (Kubát & Pattyranie, 1967; Szcześniak, Rachocki, & Tritt-Goc, 2008), and a degradation temperature of above 588 K (Yang, Yan, Chen, Lee, & Zheng, 2007). Conformation of cellulose crystals has been inspected at high temperatures, e.g. elimination of twist at the glass transition temperature (Matthews et al., 2011; Yui & Hayashi, 2007), while in some works the reverse has been shown, i.e., the cellulose tends to show a faster twisting at these temperatures (Paavilainen et al., 2011). While studies have reported experimental evidences on phase transition of cellulose at high temperatures (Hori & Wada, 2005; Wada, Hori, Kim, & Sasaki, 2010) not too much attention has been dedicated to cellulose properties at elevated temperatures in the computational community of carbohydrates. The non-reactive nature of many simulation force fields has provided the opportunity of modeling disordered cellulose by increasing the temperature of the system above the $T_g$ to disrupt the crystallinity of cellulose. Thus, current works link temperature, hydrogen bonding, and crystallinity of cellulose.

Models of cellulose I/β at 500 K demonstrated hydrogen bonding transformed from two-dimensional to three-dimensional schemes (Matthews et al., 2011). Hydroxymethyl groups changed conformations in chain groups. More frequent intrachain O3-HO3...O5 bonds occurred with fewer interchain O6-HO6...O2 bonds. Additionally, interchain bonds with O6 as acceptor dominated over those with O6 as donor. Along (1–10) planes, H-bonds occurred between O6 and glycosidic O4 and ring O5. Thus, high temperature disrupted H-bonding patterns significantly. Results suggested thermal annealing may enhance aggregation by untwisting fibrils.

An interesting finding would be that thermal annealing might be an option for enhanced aggregation of cellulose fibrils, as it triggers untwisting. This will hence decrease the accessible surface area. Combining this with the results observed by Ogawa et al. (2020) regarding cellulose drying and enhanced interfibrillar bondings, extremely packed aggregates of cellulose crystals might be achieved.

Similar responses were seen in I/β nanocrystals between 298 and 550 K (Chen et al., 2012; Zhang, Bulone, et al., 2011). A phase transition at 475-500 K increased cross sectional area and dropped density. High temperatures caused backbone rotation (30° rotation around the main axis) and changed hydroxymethyl groups from tg to gg/gt conformations between sheet layers. This deteriorated intrachain O2-HO2...O6 bonds but formed interchain ones, agreeing with Matthews et al. (2011). Comparable results were previously obtained at 300-550 K (Bergensträhle et al., 2007), showing 25% lowered modulus (156 to 117 GPa) and density (1.5 to 1.37 g/cm³) by 500 K. Thermal expansion parameters also rose compared to experimental observations. Hydroxymethyl groups again changed (tg → gg/gt) conformations with 30° chain rotation between sheets at 450 K transition temperature. High temperatures removed O2-HO2...O6 bonds.

Similar $T_g$-linked phase changes near 500 K occur regardless of model details (Bregado et al., 2021; Jiang et al., 2019; Muthoka et al., 2018). Above $T_g$, density declines and outer layers detach while inner fibrils maintain crystallinity initially. Diaz et al. (2013) also found cellulose nanocrystal films expand anisotropically with temperature, predominantly transversely between crystals as opposed to within crystals. Ultimately, higher temperatures, such as 1000 K and above can fully decompose cellulose structures, even in non-reactive MD (Khodayari, Hirn, et al., 2021), as well as cleaving glycosidic bonds in reactive MD (Paajajanen & Vaari, 2017), forming pyrolysis products. Glycosidic linkages, the mechanically and chemically weakest points, cleaving first, like what Gupta et al. (2022) observed under strain.

According to Ishida (2020), disrupting the three-dimensional H-bond pattern in cellulose is difficult even in the presence of ionic liquids. In particular, it is only at temperatures above 400 K that the total number of H-bonds starts to decrease, leading to dissolution and decomposition. MD simulations show that ionic liquids could cause faster dissolution when results are compared with those performed at high-temperature water. It must also be noted that pressure is an additional factor in cellulose stability at high temperatures. Literature suggests that high pressure can help maintain cellulose stability even at temperatures as high as 673 K (Tolonen et al., 2015).

In summary, cellulose's behavior at high temperatures, often overlooked in literature, unveils intriguing dynamics. With a glass transition temperature around 500 K, cellulose exhibits resilience to temperature variations. Studies on its conformation at elevated temperatures demonstrate a complex interplay, including the elimination or

<table><caption>**Table 4**
Effects of temperature on cellulose properties. See Fig. 3 for illustrations of the hydroxymethyl conformations.</caption><thead><tr><th>Temperature [K]</th><th>Observed changes in fibrils</th><th>Hydrogen bonding patterns</th><th>Conformation of hydroxymethyl groups</th><th>Density changes</th></tr></thead><tbody><tr><td>500 (Matthews et al., 2011)</td><td>Elimination of twist</td><td>Shift from 2D to 3D scheme</td><td>tg → gg, gt at T<sub>g</sub></td><td>–</td></tr><tr><td>298–550 (Zhang, Bulone, et al., 2011)</td><td>Increase in twist, phase transition at 475–500 K</td><td>Changes in intrachain and interchain H-bonds</td><td>tg → gg, gt at T<sub>g</sub></td><td>Decrease</td></tr><tr><td>298–660 (Bregado et al., 2021)</td><td>Increase in a and b</td><td>Disruption of crystallinity at T<sub>g</sub></td><td>–</td><td>Decrease</td></tr><tr><td>300–550 (Chen et al., 2012)</td><td>Sudden increase in a at 500 K, decrease in γ at 450 K</td><td>Change due to hydroxymethyl groups</td><td>tg → gg, gt at T<sub>g</sub></td><td>–</td></tr><tr><td>298–550 (Jiang et al., 2019)</td><td>Dissociation above 500 K</td><td>Shifted</td><td>tg → gg, gt at T<sub>g</sub></td><td>Decrease</td></tr><tr><td>463–500 (Muthoka et al., 2018)</td><td>Phase transition</td><td>Disruption</td><td>–</td><td>–</td></tr></tbody></table>
acceleration of twisting. Hydrogen bonding patterns undergo notable transformations, mainly from intrachain to interchain, influencing cellulose fibril untwisting. Thermal annealing emerges as a potential strategy for enhancing fibril aggregation due to surface hydroxyl orientation changes and untwisting of the fibrils. Majority of the current force fields reveal changes such as decrease of the density, alternation of hydroxymethyl conformation, and different hydrogen bonding patterns. Insights into cellulose's anisotropic thermal expansion and untwisting offer valuable avenues for future research, for instance, paving the way into responding to the cross-section related questions by evaluating the changes in their sorption properties at elevated temperatures. Above-mentioned points are summarized in Table 4.

### 3.5. Mechanical properties

The nanoscale size of CMFs makes direct measurement of mechanical properties difficult, requiring diffraction methods (Cheng & Wang, 2008; Diddens, Murphy, Krisch, & Müller, 2008). Nonetheless, elastic properties of crystalline regions in the CMFs are reported in the range of 130–138 GPa, measured through X-ray diffraction (Nishino, Takano, & Nakamae, 1995; Sakurada, Ito, & Nakamae, 1964; Sakurada, Nukushina, & Ito, 1962), and 145–150 GPa, measured by Atomic Force Microscopy (AFM) (Iwamoto, Kai, Isogai, & Iwata, 2009) along fibre direction. The same was reported for CNCs to be 143 GPa, measured through Raman Spectroscopy along main axis (Šturcová, Davies, & Eichhorn, 2005). A transverse elastic modulus of 18–50 GPa is also reported for CNCs by AFM (Lahiji et al., 2010). From the theoretical perspective, MD simulations have helped understand CMF mechanical behavior extensively (Khodayari, Van Vuure, & Seveno, 2019).

In 2005, Eichhorn et al. (2005) modelled cellulose deformation by changing lattice spacing c-parameter and computed cellulose chain elastic modulus as 98 GPa, slightly higher than later work (Wohlert et al., 2012), but almost as what was later reported in 2014 (Wu et al., 2014). This was related to isolated chains lacking lateral interactions present in CMFs. They argued the marginally higher value was because of emphasized intrachain hydrogen bonding without amorphous regions acting in parallel/series.

They later incorporated two modulus calculation methods – increasing cell length or applying strain – gave 155, 149, 109 GPa for cellulose Iα, Iβ, and II, respectively (Eichhorn & Davies, 2006). Removing H-bonds substantially decreased modulus, confirming their

importance in defining the mechanical properties of cellulosic structures.

Tanaka and Iwata (2006) studied 4×4×10 and 1×1×10 cellulose models with different parameter sets. Elastic modulus deviations, for models simulated with different force fields, were lower for the larger cell, confirming interchain H-bonds' role in mechanical properties.

Most early work computed mechanical properties from the potential energy change during deformation, neglecting entropy changes. This could lead to neglecting the entropic contributions. Wohlert et al. (2012) showed including entropy gives significant temperature-dependent contribution. A modulus temperature-dependence of –0.05 to –0.1 GPa K-1 was found around room temperature due to entropic effects, confirming computations require considering entropy. The modulus was computed to be 124 GPa at 0 K, 6–9 % higher than room temperature, affected by lack of entropy.

Additionally, according to Wohlert et al. (2012), not more than 20 % of cellulose elasticity was attributed to H-bonds. As fibrils get wider, average H-bonds per GU increases then plateaus (see Fig. 9). Full H-bond dependency would predict monotonic increase in the modulus, which was not evident from the results, i.e., consistent results were obtained for different crystal sizes. However, CHARMM, in contrast with GROMOS 45a4 and 56Acarbo force fields, showed higher mechanical properties for larger crystals.

Cintrón, Johnson, and French (2011) used QM and MM calculations, finding 85.2 GPa modulus for a cellulose chain, reducing to 37.6 GPa without H-bonds, confirming their significant contribution. Chain length (10–40 GU) effect gave a peak modulus at 20 GUs. Aggregated chains showed little effect from interchain H-bonds on the modulus. However, entropy effects were again neglected in these computations.

Muthoka et al. (2020) studied cellulose crystals under tension and shear using steered molecular dynamics (SMD) simulations. In SMD, a form of non-equilibrium enhanced sampling method, one or a group of atoms are enforced motion along a certain direction. Due to the non-equilibrium nature of these techniques, results must be usually well-benchmarked. As an example, SMD simulation results have been shown to be extremely force-rate dependent (see for instance, (Lemkul & Bevan, 2010)). 161 GPa tensile modulus was found for 18 chain cellulose, with applied force magnitude affecting H-bonds. The H-bonds dropped by increasing the pulling force magnitude. This decrease of the H-bonds was then attributed to the change of the *tg* conformations of the hydroxymethyl groups.

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** The relationship between the cross section size of the fibrils (different surface to core ratio) and the number of intramolecular H-bonds per GU. The smaller the fibrils are, the lower the number of H-bonds. Data to generate the figure is taken from (Wohlert et al., 2012).</figcaption></figure>
Nanoindentation simulations of cellulose Iβ demonstrated transverse modulus of 5.1 GPa along the (1–10) face (Wu et al., 2013), lower than what reported for the (100) face by Geng et al. (2023), i.e. 17 ± 2 GPa. Large anisotropy gave 28.8, 7 GPa modulus along b, a directions and 139.5, 120.3 GPa with and without H-bonds, respectively, along the crystal axis. Surprisingly, H-bonds only changed modulus by 14 %, contrasting some predictions of their stronger role, as was recently discussed by Wohlert et al. (2021) in 2021. On the contrary, this is not entirely coherent with what was found by Khodayari, Hirn, et al. (2021). Lack of H-bonds distorts the crystallinity of the fibrils, which can lead to a significant decrease in the mechanical properties. It has been shown that elastic properties could drop to 50 %, should approximately 60 % of the H-bonds be lost in a CNC.

Ciesielski et al. (2019) studied bending deformation effects in cellulose. Breakage criteria were based on earlier QM bond length calculations. They observed that during bending simulations, the length of the C4-O4 bond changes from 1.41 to 1.58 Å, whereas the C1-O4 bond stretches from 1.39 to 1.88 Å. The authors concluded that this leads to the primary breakage of the C1-O4 bond. The issue with the non-reactive force fields when computing the mechanical properties, despite being extremely faster than the reactive force fields, is strength computation, as bond breakages are not modelled with such potential parameters. Several works incorporated reactive potentials like ReaxFF to enable modeling bond breaking for strength properties. Reactive MD suggests bonds fail based on synergistic changes in multiple structure parameters, not bond-lengths alone (Gupta et al., 2022). An ultimate tensile strength of 9.2 GPa at 8.5 % of strain is reported, thanks to such techniques (Gupta et al., 2022).

According to Hao et al. (2018), the glycosidic angle at failure is 180°, while the equilibrium value of the glycosidic angle is 117° (Nishiyama et al., 2002). Hao et al. (2018) suggests that the failure of cellulose is primarily due to the breakage of hydrogen bonds, slippage of the sheets on each other, and their rotation. However, it must be noted that the failure of covalent bonds in cellulose requires much higher forces than those breaking hydrogen bonds. While hydrogen bonds play a crucial role, they cannot be the main factors in the failure mechanism of cellulose. Although proposed factors may initiate and trigger the failure mechanism, the ultimate failure properties must be beyond the control of hydrogen bonds.

In another study, Wu et al. (2014) found stable modulus and Poisson's ratio over strain rates, while yield/failure properties increased with rate. Elastic modulus again showed anisotropy along different directions. Axial modulus exceeded transverse by an order of magnitude.

<table><caption>**Table 5**

Elastic modulus of cellulose models reported in the literature.</caption><thead><tr><td>Cellulose model</td><td>Modulus [GPa]</td></tr></thead><tbody><tr><td>Iβ (Wu et al., 2014)</td><td>107.8–112.9</td></tr><tr><td>Iβ (Wohlert et al., 2012)</td><td>113–140</td></tr><tr><td>Iβ (Chen et al., 2023)</td><td>119</td></tr><tr><td>Iβ (Pang et al., 2022)</td><td>120</td></tr><tr><td>Iβ (Wu et al., 2013)</td><td>120.3–139.5</td></tr><tr><td>Iβ (Tanaka & Iwata, 2006)</td><td>124–155</td></tr><tr><td>Iβ (Hao et al., 2018)</td><td>142</td></tr><tr><td>Iβ (Sturcová et al., 2005)</td><td>145</td></tr><tr><td>Iβ (Khodayari, Van Vuure, et al., 2020)</td><td>146.6</td></tr><tr><td>Iβ (Khodayari, Hirn, et al., 2020)</td><td>146.7</td></tr><tr><td>Iβ (Kulasinski, Keten, Churakov, Derome, & Carmeliet, 2014)</td><td>150</td></tr><tr><td>Iβ (Khodayari, Hirn, et al., 2021)</td><td>152</td></tr><tr><td>Iβ (Bergenstråhle et al., 2007)</td><td>156</td></tr><tr><td>Iβ (Paavilainen et al., 2012)</td><td>157</td></tr><tr><td>Iβ (Muthoka et al., 2020)</td><td>161</td></tr><tr><td>Iα, Iβ, II (Eichhorn & Davies, 2006)</td><td>155, 149, 109</td></tr><tr><td>II (Eichhorn et al., 2005)</td><td>98</td></tr><tr><td>Single-chain (Wu et al., 2014)</td><td>90–95</td></tr><tr><td>Single-chain (Cintrón et al., 2011)</td><td>85.2</td></tr><tr><td>Single-chain (Wohlert et al., 2012)</td><td>80</td></tr><tr><td>Single-chain (Paavilainen et al., 2012)</td><td>4</td></tr></tbody></table>

Single chain and crystal ultimate stress/strain matched, attributed to brittle failure of cellulose crystals. Accordingly, as soon as a first chain fails, the remaining chains fail without additional strain, therefore, the ultimate strains and ultimate stresses would match. However, this does not explain their different elasticity. In particular, a different elasticity for a single chain and that of a whole crystal was observed, while, different elasticity could not result to both equal ultimate stress and ultimate strains in brittle failures. Hence, a nonlinearity or strain hardening should have been present in either of the cases in order to explain such observation in the ultimate properties. A non-linear trend in the stress-strain curves for cellulose has been already observed by Hao et al. (2018).

Paavilainen et al. (2012) studied nanocrystals with introduced defects, finding modulus dropped proportionally to the defected chains. In particular, models with ten connecting chains showed an elastic modulus of 54±2 GPa and those with three connecting chains had a modulus of 11.9±0.2 GPa. The authors concluded that each connecting chain contributes to 4 GPa of the (apparent) stiffness. Nevertheless, the actual elastic modulus (and not the apparent) of a single cellulose chain is calculated to be about 80–95 GPa in the presence of H-bonds (Cintrón et al., 2011; Wohlert et al., 2012; Wu et al., 2014). Furthermore, as shown before, the absence of H-bonds can lead to a drop of the modulus to values of 35 GPa.

The computational mechanical properties of cellulose models reported in the literature are summarized in Table 5.

In short, the computational exploration of cellulose's mechanical properties has evolved significantly, revealing crucial insights into its strength, primarily attributed to hydrogen bonding. Early models, exemplified by Eichhorn et al. in 2005 and 2006 (Eichhorn et al., 2005; Eichhorn & Davies, 2006), laid the foundation for understanding cellulose's elastic modulus, considering the impact of hydrogen bonding. Subsequent studies highlighted the previously underestimated entropic contributions, challenging conventional methods when measuring the mechanical response of cellulose nanocrystals. Older studies underscored the substantial role of hydrogen bonds in cellulose, while almost a decade later it was elucidated that they should not be considered as the main components in control of the mechanical properties. Reactive molecular dynamics like those performed by Wu et al. (2014) and Hao et al. (2018) exemplify the importance of incorporating reactive force fields for better understanding the ultimate properties of carbohydrates, as well as their reactive behavior with environmental substances. Yet, currently, the need for rethinking, benchmarking, and developments can still be extensively sensed.

### 3.6. Cellulose-cellulose interactions

Understanding cellulose-cellulose interactions is paramount, as these mainly sustain macrofibril formation, influencing strength, stability, and integrity. Hence, delving into the specifics of these interactions could unravel an essential milestone to better understanding of cellulose. Such interactions primarily involve H-bonds and an interplay between vdW and electrostatic terms (Thu et al., 2022), which can mediate the macromolecular structures and their properties.

For fibril aggregation forces, H-bonds and electrostatic interactions were shown to be responsible (Sáenz Ezquerro et al., 2019). SMD simulations between the (100) faces of CMFs showed that the calculated bond strengths are noticeably force-rate dependent, (a common consequence of SMD simulations). The magnitude of the calculated forces in their studies was two to three orders higher than those measured experimentally. Hence, particular attention should be given to when SMD simulations are performed for evaluating the interaction potentials. In particular, the minimum free energy surface might be surpassed should the wrong set of parameters in a SMD simulation is used.

In studying paracrystalline matrix interactions, SMD and umbrella sampling as well as free energy calculations using Jarzynski's equality (Jarzynski, 1997) between non-polar cellulose surfaces showed limited
direct cellulose adsorption due to an intervening water layer (Oehme, Doblin, et al., 2015). However, similar surfaces or polar ones likely enable direct interactions through H-bonds, stabilizing the aggregation.

MD also allows for precisely observing the changes in the effect of surface modifications on the interactions energies between different components. Characterization of surface-modified cellulose has also been performed and reported in the literature. Umbrella sampling simulations between native, sulfated, and TEMPO-oxidized CNCs showed native cellulose had higher interaction magnitudes (Garg et al., 2020). Cation type (Na⁺ or Ca²⁺) also impacted interaction strengths. Importantly, pushing and pulling simulations between CNCs exhibited different potential energy characteristics, attributed to water's barrier role.

One must note that, in umbrella sampling simulations of cellulose nanocrystals, the sampling duration could drastically affect the obtained results. In particular, Matthews et al. (2012) showed that noticeable structure deviations, such as untwisting of fibrils, could take place when cellulose simulations were performed for hundreds of nanoseconds. Such structural deteriorations can significantly influence the computed interaction energies. For instance, assessment of interactions between non-functionalized and functionalized fibrils was reported based on the umbrella sampling methods, and contrary to Garg et al. (2020), showed weaker interactions between non-functionalized models compared to the surface modified ones (Paajanen et al., 2016). Longer sampling could possibly lead to changes in the inter-fibrillar interactions.

Besides moisture, temperature impacts have additionally been quantified. Absolute free energy of aggregation decreased with more water and lower temperatures, explained by deteriorated hydration shell repulsion (Silveira, Stoyanov, Kovalenko, & Skaf, 2016). Findings concur with weakening cellulose-water attractions at high temperatures enabling cellulose-cellulose interactions (Dri et al., 2015; Ogawa et al., 2020).

Other factors defining the aggregation mechanism is surface acetylation. Surface acetylation impacts cellulose-liquid and cellulose-cellulose interfaces by disrupting crystallinity and decreasing fibril work of adhesion in water (Chen, Lo Re, et al., 2020), despite increased hydrophobicity. Additionally, surface acetylation causes a disruption of the crystallinity at the interface of the aggregated fibrils by influencing the hydrogen bonding and vdW interactions. As the change in the energies between cellulose fibrils is greater than that of the hydrophobization, a significant decrease in the work of adhesion between fibrils is observed. In other words, acetylation leads to a lower tendency for the cellulose fibrils to aggregate in water. Recent research shows that acetylation on the surface of CNCs can lead to weaker CNC-CNC adhesion, more than the declined CNC-polymer work of adhesion. Hence, a better miscibility can be achieved, indirectly also controlling the wettability of the CNCs (Wohlert et al., 2023).

Shear mechanisms between CNCs/fibrils have been also studied (Molnár et al., 2018), exhibiting stick-slip tendencies and shear strength mainly from continuously rearranging H-bonds and dispersion forces (Khodayari, Thielemans, et al., 2021; Zhang, Keten, et al., 2021). The highly nonlinear and anisotropic shear mechanism dominates tensile strength of crystals (bearing lower ultimate stresses), while introducing dislocations in the middle of crystals (generating more realistic models for CMFs) could generate more plastic behavior.

These investigations unfold a panoramic view of cellulose interactions, from the influence of surface modifications on aggregation energies to the intriguing dynamics of shear mechanisms within CNCs and fibrils. Understanding the delicacy of the role of hydrogen bonding, electrostatic forces, and surface alterations on these interactions not only illuminates the fundamental science of cellulose but also holds promise for applications spanning materials science and biomimetic design. In particular, as cellulose-cellulose interactions are mainly controlled by H-bonds and vdW, surface modification of cellulosic fibres can drastically change the aggregation dynamics. Additionally, considering that long dynamics and relaxation phases incorporate structural

changes, e.g. twist rates, researchers much carefully take into account the computational time required for quantitative prediction of such aggregations.

### 3.7. Cellulose and water

Arguably one of the most critical parameters in defining plant cell wall properties is moisture content (MC), and relative interaction of the constituents of the plant cell wall with the contented water (Etale, Onyianta, Turner, & Eichhorn, 2023). Humidity drastically impacts the mechanical properties of plant-sourced fibres, induced by swelling and general plasticity in the plant cell wall, mainly caused by positioning of water molecules at interfaces, leading to dissociation and shear softening (Khodayari, Thielemans, et al., 2021). MD simulations have provided insights into water diffusion in plant fibres.

When CNC/CMF interactions are considered (in absence of other constituents), hydration mainly occurs due to hydroxyl-water interactions. Hydration affects interchain properties when chains are solvated, or induces swelling. MC of 10–15 % signifies the start of cellulose fibre aggregate swelling (Paajanen et al., 2022). Research shows MC can lead to water penetration between CNCs/CMFs providing a favorable electrostatic environment by decreasing net potential energy, compensating for water entropy loss. This adhesive water behavior can decrease material ductility, improving brittleness (Chen, Wohlert, et al., 2022).

Comparing the behavior of separated cellulose chains, in vacuum

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** a) Mean Square Displacement (MSD) of the surface and core hydroxymethyl groups (representative of the mobility of selected groups) versus temperature, for two different hydration levels of 5 % (dashed lines) and 20 % (solid lines). b) The persistence length of the chains versus temperature for the 20 % hydration level ($h = 0.2$). The trend follows the MSD of the hydroxyls' hydrogen shown above. Figures adapted and redrawn with permission from (Petridis et al., 2014).</figcaption></figure>
and aqueous environments, resulted that chains aggregated quickly via H-bonds without water (Tanaka & Fukui, 2004). When solvated, this aggregation was weaker, separated by water. A similar phenomenon is seen in CMF clusters.

CMF cluster hydration usually causes overall fibre swelling mainly
due to cellulose surface hydrophilicity and accumulated interfacial
water with slower diffusion than bulk water (Ito & Matsumoto, 2022).
Once hydrated, full cellulose dehydration is difficult due to strong
confined water-cellulose interactions. While removing non-confined
water is possible, removing confined water is extremely difficult. Sim-
ulations of a cellulose crystal indicated two water diffusion coefficients,
confirming two states (O'Neill et al., 2017). Increasing temperature
mobilizes non-confined water, gaining mobility at ~260 K. This water
usually fills CMF interspaces.

Differently hydrated cellulose structures induce various effects. 5 % and 20 % systems were modelled to study mechanical response, corresponding to dry and hydrated cellulose (Petridis et al., 2014). In the 20 % system, higher hydroxymethyl mobility was observed compared to pyranose rings, accentuated with more moisture (Fig. 10a). This increased agility increased hydroxymethyl $tg \rightarrow gt/gg$ reorientation, further enhanced by higher temperatures, especially above 240 K.

The persistence length trend overlapped the hydroxymethyl reor-
ientation, indicating possible rigidity increases via hydroxymethyl dis-
turbances, improved by moisture and temperature (Fig. 10b). Thus,
modest hydration may enhance single fibril properties. However, while
single-component properties may increase with humidity, inter-fibrillar
plasticization likely controls overall CMF composite behavior, which
can be mediated by surface modifications (Matsuba, Kubota, &

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** a) Wetting simulations of (1–10) face of cellulose with water at the nanoscale. b) The $tg$ conformation of the hydroxyl groups on cellulose's (110) and (1–10) surfaces. The orientation of hydroxyl groups on the (110) face allows for improved hydrophilicity compared to the (1–10) face. Figures reprinted with permission from (Trentin et al., 2021).</figcaption></figure>
Matubayasi, 2022).

Cellulose exhibits packed crystallinity making it water-impermeable, yet its surface hydroxyls allow wettability (Heiner & Teleman, 1997). Different cellulose faces behave variably due to hydroxyl densities, however cellulose surfaces are generally considered hydrophilic. In other words, when the term"hydrophobic surface" is assigned to one of the cellulose faces, it basically refers to less hydrophilicity compared to its other faces. Charged systems like sulfated CNCs show expanded bundle gaps affecting interfacial waters (Garg et al., 2020).

Trentin et al. (2021) modelled different cellulose faces to investigate the hydration of the cellulose faces (Fig. 11a). According to their results, the (010) and (1-10) faces of Iα, (010) and (110) faces of Iβ, as well as the (100), (1-10), (010), and (110) faces of cellulose II and III showed a highly hydrophilic behavior. The computed contact angles were all in the range of 11° or lower. The (1-10) face of Iβ cellulose showed a less hydrophilic behavior with a contact angle of about 32°.

The authors used the low-bond axisymmetric drop shape analysis (LBADSA) plugin in ImageJ to fit an arc to the boundary of equilibrated water droplet onto the cellulose surfaces. Defining surfaces in MD remains challenging; the water-air border has been shown to best be where the density of water reaches half of the bulk density (Benhassine, Saiz, Tomsia, & De Coninck, 2009; De Coninck & Blake, 2008; de Ruijter, Blake, & De Coninck, 1999; Herrera, García, Atilan, & Aparicio, 2015; Sridhar et al., 2023; Webb & Grest, 2002). The discrepancy between the results here and other simulations may arise from differences in methodologies.

Additionally, the water model used can significantly alter results. TIP4P/2005 showed a contact angle of 16°, while TIP3P showed complete wetting of cellulose (010) (Malaspina & Faraudo, 2023). This relates to different surface tensions of water models, also observed by Sridhar et al. (2023). Depending on the model, computed contact angles differed: TIP3P completely wetted; SPE/E showed 11.6°; TIP4P/2005 showed 15.6°. Angles were higher with acetylated surfaces (17.9° to 34.3°).

Water may also deteriorate crystalline cellulose surfaces by inducing amorphization (Khodayari, Van Vuure, et al., 2020). Hence, despite hydroxyl/hydroxymethyl freedom in Trentin et al. (2021)'s simulations, constraints on cellulose models could be another reason for differences. An interesting conclusion was that, unlike the general belief of hydrophobic cellulose surfaces like (100) and (200), these are just less hydrophilic, not hydrophobic.

Another observation was the difference between cellulose (110) and (1-10), the former showing more hydrophilicity (Trentin et al., 2021). This arises from dihedral angles of surface hydroxymethyl groups. The O6 on (110) of cellulose Iβ points outward (gt) unlike on (1-10) (tg), as shown in Fig. 11b. tg allows intrachain hydrogen bonding, lessening hydroxyl availability to interact with water, reducing wettability. This was also seen when comparing Iα (010) and (001).

The reorganization of the surface hydroxymethyl groups and, consequently, modification of the wettability properties has been studied thoroughly in the literature. Maurer et al. (2013) modelled cellulose Iα (100) face, Iβ (100) face, and cellulose II (010) and (120) faces, all consisting of 4 layers of cellulose chains (Maurer et al., 2013). They showed that cellulose surface slabs reconstruct based on the allomorph under study. For instance, the (100) face of Iα shows a quick irregularity after equilibrium by switching the interchain H-bonds to interchain H-bonds. This phenomenon was less observed in the other two cellulose models studied.

The authors concluded that three parameters define the reconstruction at the surface of cellulose: 1) the phase on the top of cellulose surface (vacuum/solvent), 2) the cellulose allomorph, and 3) the studied crystal face. Another important observation was that even cellulose models such as cellulose II, which are not as tightly packed as cellulose I cannot show significant solvation and swelling, at least within the mesoscopic time frame of the simulations.

Contact angles of 16° and 23° were found for the (010) and (100) planes, when Iβ cellulose surfaces were wetted, in another study (Malaspina & Faraudo, 2019). The difference with those obtained by Trentin et al. (2021) could possibly come from the difference in the number of water molecules used in both studies. Moreover, in this work as well, the authors mainly used a 2-D density processing methodology to define the border between water and vacuum space around and to calculate the contact angle. However, as also discussed before, distinguishing the correct border of the liquid could potentially alter the results.

Wetting dynamics of the (100) and (110) faces of cellulose Iβ with SPC water models were also studied in 2008 (Mazeau & Rivet, 2008). The contact angles were estimated 43° for the (110) face and 90° for the (100), noticeably higher than other studies (Malaspina & Faraudo, 2019; Trentin et al., 2021). Above-mentioned parameters, including water model, droplet size, simulation time, or the cut-offs distinguishing the droplet edges, etc. could be possible reasons behind the observed differences in these simulations.

Cellulose surface can undergo changes due to the hydrolysis, during delignification, or CNCs production. Simulations on sulfated CNCs showed that water completely spreads over the (010) face, whereas the (100) face showed a contact angle of about 9° (Ma et al., 2021). This angle is relatively lower than that observed in other studies (Malaspina & Faraudo, 2019; Trentin et al., 2021). This could be the consequence of the sulfate ester groups on the surfaces, facilitating better hydration and increasing the hydrophilicity of the surfaces. Interestingly, it was observed that the contact angle of the (100) was just smaller with hexadecane than water, attributed to the lower surface tension of hexadecane.

TEMPO-oxidized cellulose has negatively charged carboxylate groups on the surface instead of hydroxymethyl groups. Oxidation is considered to improve dispersion of cellulose in water at room temperature. Simulation shows TEMPO-oxidized glucan chains exhibit a threefold helical screw structure, with outward-pointing -COO⁻ groups (Asgarpour Khansary et al., 2020). These negative charges likely cause better dispersion, separability, and increased hydrophilicity.

Mudedla et al. (2021) studied oxidation defects on crystalline cellulose's (100) face. They modelled cellulose structures with defects from C1 and C4 carbon oxidation. Defects broke glycosidic bonds, generating different glucose derivatives. Chains had one or two oxidations, either close or distant. Results showed more water binding to defective surfaces, increasing hydrophilicity. Local chain protrusions occurred outward from the surface, raising flexibility and solvent-accessible surface area.

In recent work, Garg et al. (2021) modelled crystalline and amorphous states of TEMPO-oxidized and carboxymethylated cellulose. According to the results, MC (water weight% on dry basis) was increased by increasing the RH or decreasing the temperature. Additionally, moisture uptake in amorphous cellulose was shown to be higher than crystalline cellulose, regardless of the modification types.

Simulations have inspected water adsorption in more complex systems, like crystalline cellulose with hemicelluloses (Galactoglucomannan: GGM) (Kulasinski, Guyer, Derome, & Carmeliet, 2015a). Results showed water filled the cellulose-GGM interface as MC increased. More water resided in GGM and the interface, decreasing GGM density and swelling the structure. Less water showed ordered binding, decreasing with higher MCs. Swelling decreased Young's modulus, with a noticeable shear stiffness drop at 5 % MC, beyond which swelling occurred. Below 5 % MC, stiffness slightly increased, attributed to water forming H-bonds between polymers. This increase was previously also shown to be attributed to the bridges facilitated by H-bonds made between water and the polymer (Petridis et al., 2014). Thus, swelling mainly comes from free water over bound water.

Kulasinski et al. (2017) later added lignin to cellulose-GGM models to expand simulations of the plant cell walls. Less-ordered cellulose crystal surfaces acted as amorphous boundaries. At low MCs, water occupied GGM's empty sites, while strongly-bound water resided at
cellulose-amorphous interfaces due to cellulose's high order. MD brought insights on water interactions unexamined experimentally.

Breaking down the model's adsorption isotherms showed the lignin-hemicellulose interface contributed very little to total adsorption. Hemicellulose adsorbed the most moisture, while interfaces became more important at higher RHs. Over 40 % of hemicellulose's (dry) internal H-bonds were lost due to moisture, and replaced to hemicellulose-water H-bonds. However, lignin's H-bonds increased. Lignin slightly shrank due to moisture uptake and increase of pore volumes, putting a negative strain on lignin, and hence, resulting in more H-bonds. Still, H-bonds decreased at all polymer-polymer interfaces with moisture.

Elastic and shear moduli substantially dropped with higher MCs in three stages. First, small pores formed and diffusivity was low, with slight mechanical property drops but no H-bond loss. Next, swelling, H-bond loss, and noticeable property drops occurred from increased porosity. Finally, more pores formed larger amorphous phase clusters, drastically reducing mechanical properties. Such self-aggregating behavior of lignin, forming nanodomains, connecting to hydrated hemicellulose structures (xylans) has been well-studies recently through ssNMR and dynamic nuclear polarization (Kang et al., 2019).

While swelling in cellulosic materials takes place due to the moisture uptake of the surface of cellulose, drying the fibrils does not lead to full deswelling of the system. Zitting et al. (2021) simulated cellulose microfibril bundles with hemicellulose at 110 % MC to study swelling and drying effects. One non-periodic 7-fibril bundle and one 4-fibril periodic bundle were modelled. Although fibril-fibril distances decreased during drying, fibril morphology did not change. In the periodic bundle, water diffusivity increased with higher MCs, saturating above 40 % MC.

Separately, Ogawa et al. (2020) studied CNC deformation during drying to understand cellulose fibre hornification, i.e. lateral shrinkage of cellulose when dried, without recovered properties. One CNC was fixed while the other moved freely. Drying was modelled by removing isolated water molecules. Results showed CNCs deform laterally during drying due to water surface tension, explaining cellulosic fibre shrinkage and stiffening.

As discussed before, swelling with excessive increased MC reduces mechanical properties, while in low MC, enhanced properties could be observed (Kulasinski et al., 2017; Kulasinski, Guyer, Derome, & Carmeliet, 2015a). Lindh et al. (2016) provided an explanation, finding CNC elastic modulus relies mainly on intrachain H-bonds and not the interchain ones. As the intrachain H-bonds are preserved with water, properties only marginally drop at moderate hydration levels.

Similarly, Hou et al. (2021) saw moist CNC complexes strain hardened under tension, increasing failure stress and strain. Experiments confirmed strain hardening occurring between 30 and 50 % relative humidity. Though, these humidity levels are much higher than Kulasinski, Guyer, Derome, and Carmeliet (2015a)'s swelling threshold.

While cellulose shows a hydrophilic nature, it is insoluble in water due to hydrophobic stacking and hydrogen bonding between chains, as shown via molecular dynamics simulations. Single cellulose chains is shown to prefer an extended twisted conformation rather than behaving as flexible coils in water (Bergenstråhle, Wohlert, Himmel, & Brady, 2010). However, cellulose-water suspensions demonstrate an effective change in thermal properties like heat capacity and conductivity (Deng et al., 2020). Specifically, the axial thermal conductivity of crystalline cellulose is superior to lateral conductivity when dispersed in water.

Cellulose dissolution can potentially be feasible through surface modifications (Shi, Li, & Zhao, 2022). For CNCs, dissolution can take place in supercritical water (Tolonen et al., 2015). Simulations revealed thermodynamically favorable dissolution of cellulose nanocrystals in supercritical water models, while ambient water caused only minor twisting. Intercalation of ions and formation of H-bonds has also been explored as a mechanism for dissolving cellulose using ionic liquids (Sánchez-Badillo et al., 2021). Surface deterioration could be the trigger to start the dissolution of cellulose. Kang, Wang, Zhang, and Zhou

(2021) examined the dissolution of cellulose carbamate for fibre spinning in water, water-NaOH, and water-NaOH-ZnO solutions. Simulations confirmed carbamate cellulose forms metastable H-bonds with NaOH. Adding ZnO intensified the H-bond network, producing a stable spinning dope.

On the other side of solubility and chain dissociation, regeneration of cellulose has also been studied in aqueous environments, where water enhances cellulose-cellulose interactions, especially at higher temperatures. However, the regenerated cellulose models show hydroxymethyl group orientations different from that of crystalline cellulose and hence exhibit an amorphous state (Gupta, Hu, & Jiang, 2013b). Anti-solvents like acetone can also promote aggregation through hydrogen bonding when breaking cellulose-anion bonds (Gupta, Hu, & Jiang, 2013a). Comparing solvation in water and organic solvents like cyclohexane, simulations indicate dissociation of cellulose chains proceeds smoothly in water but step-wise in cyclohexane (Bergenstråhle et al., 2009).

Investigations into different moisture content in cellulosic materials uncover diverse effects on cellulose systems. Simulations of differently hydrated cellulose structures reveal the nuanced impact of moisture on the mechanical response of cellulose. Particularly at humidity range of 10–15 % significant swelling of cellulose aggregates can be observed, attributed to the changes of the hydrogen bonding on the surface of crystals/fibrils. The agility of distinct cellulose subgroups, such as exposed and buried hydroxymethyls and pyranose rings, varies with moisture levels, with minor MC (less than 5 %) inducing stiffening effects. While many studies targeted studying different cellulose surfaces hydration, achieving consensus has proven challenging due to the difference in the water models, force fields, and MD solvers. In this regard, users may take advantage of methodologies, such as those offered by Sauter and Grafmüller (2015), to adjust target properties with water models (see Section 2). Nevertheless, what MD community agrees on is that cellulose primarily shows a hydrophilic character, with different faces showing varying degrees of hydrophobicity. In a more general perspective, the intricate details of the plant cell wall interactions, particularly in the context of moisture content, offer a comprehensive understanding of the dynamic behavior of plants and plant fibres. To this end, it is essential to also seek the interactions between cellulose and other constituents, such as hemicelluloses and lignin (Section 5.3).

### 3.8. Coarse-gained models

This review primarily aims to provide a comprehensive atomistic overview of the current progress in modeling carbohydrate polymers. However, we will briefly introduce the perspective of building upon current atomistic studies by developing larger size/time scale models. Therefore, we will only superficially discuss the capabilities of such strategies without delving into much detail in this review. Coarse-grained (CG) molecular dynamics allows modeling of larger systems by grouping atoms into beads with singular centers of mass (Ramezani & Golchinfar, 2019). CG enables longer simulations of phenomena like diffusion while losing some accuracy. Considering many different methods for performing CG-MD, such as force matching (Izvekov & Voth, 2005), Iterative Boltzmann Inversion (Moore, Iacovella, & McCabe, 2014), or machined learned CG-MD (Wang et al., 2019), CG models are often parameterized for specific properties rather than general behavior. Several CG cellulose models have been developed for cellulose.

Martini force field and its usage in CG have been a popular option in the modeling community (Pang et al., 2022). Martini CG force field was used to model a 36-chain DP40 cellulose crystal compatible with proteins (Wohlert & Berglund, 2011). It gave reasonable diffusion behavior at room temperature.

Another CG Iβ cellulose model with 3 sites per glucose unit matched lattice parameters except the *a* axis (Bu et al., 2009) (Fig. 12).

Recently, a supra-CG model with implicit water, simulating assembly and properties like modulus was modelled (Mehandzhiyski et al., 2020).
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** The representation of a) the beads superimposed on the atomistic structure, b) a cellulose chain, and c) the coarse-grained equivalent of the atomistic model. Figure adapted and redrawn with permission from (Bu et al., 2009).</figcaption></figure>

The model contained bead of 20 GUs, having five cellulose chains of 4-GU long, positioned on the surface and a core of 6 chains of 4-GU long beads. The model was able to exhibit proper twisting angle, self-diffusion coefficients, and axial elastic modulus. Conventional evaporation modeling was not possible, so applying high pressure served as

proxy.

As discussed before, H-bonds are essential parameters defining the properties of cellulosic systems. Wu, Beltran-Villegas, and Jayaraman (2020)'s hydrogen bonding CG model simulated ordered aggregation of disordered chains, matching all-atom and crystallographic data. The

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Representation of the dislocations within the CMFs. Molecular dynamics simulations offer a range of 3.1–5.8 nm for such regions in 36-chain models. Figure reprinted with permission from (Khodayari, Hirn, et al., 2021).</figcaption></figure>
aggregated models then showed perfect glycosidic linkage, interchain H-bond pattern, and relative alignment. Considering mechanical properties, for example, Shishehbor and Zavattieri (2019)'s model captured CNC modulus and shear mode stick-slip between fibrils. They reported a value of 170±30 GPa for the tensile modulus, which is comparable to an all-atom simulation.

While there exist many other studies aiming to represent a CG model for cellulose (Fan & Maranas, 2015; Glass, Moritsugu, Cheng, & Smith, 2012; Hynninen, Matthews, Beckham, Crowley, & Nimlos, 2011; López et al., 2015; Markutsya et al., 2012; Poma, Chwastyk, & Cieplak, 2016; Ramezani & Golchinfar, 2019; Srinivas, Cheng, & Smith, 2011, 2014; Wu, Zhan, Wang, & Ju, 2012), current research seem to be progressing in the correct direction, which can eventually help in generating and studying a complete model of the plant cell wall of considerable size. Steps to efficiently valorize such progresses include offering bottom-up approaches to computations of plant fibre mechanics. Examples to these steps could be what was performed by Adler and Buehler (2013) in generating models which are able to identify mechanical behavior at various microfibrillar angles, represent comparable stress-strain relationships with experimental measurements, and to incorporate other cell wall components, such as hemicellulose alongside cellulose. Current state-of-the-art offers models for different cell walls able to predict various mechanical aspects such as elasticity, stiffening, and plasticity beyond yield point of plant fibres (Zhang, Yu, et al., 2021). Notably here, only a few examples were selected to convey our message in capabilities of CG strategies in that regards. Further steps should consider the incorporation of other cell wall components, such as lignin, as well as hydration effects. Experimental evidence should inspire the generation of these plant cell wall models.

## 4. Disordered cellulose

As discussed in the previous sections, several parameters could lead to the generation of non-crystalline segments. In particular, non-crystalline cellulose exists in two main different forms, namely amorphous cellulose and dislocated cellulose. Both of these disordered media exhibit different properties and shapes. Amorphous cellulose refers to isolated or entangled cellulose chains either man-made or found on crystal surfaces. Though, there are limited evidence that such form of cellulose actually exists in plant cell walls. On the other hand, CMFs themselves are series of crystalline cellulose segments connected via lower crystallinity regions, also referred to as strain-induced regions or dislocations. Dislocations refer to regions a couple of nanometers in length where the crystallinity of the segments is disrupted, probably due to the induced strain in the process of plant growth (Fig. 13) (Eichhorn et al., 2005; Habibi et al., 2010; Hidayat, Felby, Johansen, & Thygesen, 2012). Another form of amorphous cellulose can be also obtained by spin-coating cellulose derivatives such as trimethylsilyl cellulose followed by regeneration via the gas phase (Kontturi et al., 2011; Kontturi & Spirk, 2019). The resulting thin films are highly amorphous and feature a 2D short range order over 30–60 Å (Jones et al., 2020). They have been also employed as models to explore water–cellulose interactions for both swelling and drying of cellulose (Ehmann et al., 2015; Kittle et al., 2011; Reishofer et al., 2022). While recent research has focused more on modeling amorphous cellulose, dislocated cellulose in cellulose microfibrils warrants further study. Here, some works on disordered cellulose, its generation procedure, and properties are summarized.

### 4.1. Techniques and properties

The generation of disordered regions in MD simulations of cellulose has mainly been performed by either insertion of random chains within simulation boxes, or disrupting crystallinity through high temperatures, above $T_g$ or even oxidation temperatures of cellulose, without breaking bonds. The latter fully changes hydroxymethyl group conformations and

hydrogen bonding patterns, resulting in random structures that still show preferential group orientations (Khodayari, Hirn, et al., 2020).

In 2003, amorphous cellulose was built by randomly inserting DP10/20/40 chains equilibrated at room temperature (Mazeau & Heux, 2003). Hydroxymethyl groups showed *gg* and *gt* conformations unlike the *tg* conformation in crystals. Glycosidic bond torsional angles, however, showed values according to the defined potential parameters in the force fields, restricted from any change. The calculated glass transition temperature was 650 K, higher than experiments and that of crystalline cellulose.

Paracrystalline (dislocated) structures have been modelled as disordered states of crystal models (Kulasinski, Keten, Churakov, Derome, & Carmeliet, 2014). The models were made by tempering the crystalline cellulose to 700 K and quenching to 300 K. Amorphous cellulose was also modelled by inserting DP10–40 chains into simulation box and performing same tempering/quenching procedures. The achieved amorphous density was 14 % lower than crystals. Temperling decreased $\gamma$, increased $\alpha$, while decreasing $c$ lattice parameters, as also seen by Mazeau and Heux (2003). Transition occurred around 450–550 K. Crystalline structures retained $tg$ conformations while others showed *gg* and *gt*. $\phi$ and $\psi$ probability density distributions drastically changes in non-crystalline models. Crystalline modulus values were anisotropic unlike the isotropic amorphous structure.

Applying similar tempering/quenching techniques resulted that small segments, i.e., DP4 models, could fully amorphize but longer chains exhibited paracrystalline cores with amorphized surface chains due to restricting interchain bonds (Bregado et al., 2019). This is also confirmed later in 2021 (Khodayari, Hirn, et al., 2021) that when models shorter than 6 GU were bound between two crystalline segments, breakage of the H-bonds becomes challenging. Dislocations as long as 25 GU in CMFs with mechanically graded zones have been modelled recently, while, overall, such models might not be realistic (Chen et al., 2023). More experimental investigations are needed to clarify size distribution of the dislocations in CMFs. Two mechanisms were proposed by Bregado et al. (2019) during these processes: changes along and across cellulose chains, with back-and-forth slippage and with rotations around the main axis, respectively. Cross-sectional changes were faster for $\alpha$ structures, while $\beta$ models showed more pronounced longitudinal shifts. Spectroscopy, radial distribution functions, and XRD results showed inner chains were shielded with a slightly paracrystallized core in high DP models.

Disorder depletes mechanical properties. It was shown in Section 3.2 that intrachain H-bonds are crucial parameters when mechanical properties of cellulosic materials are considered. Amorphous cellulose modulus is less than 10 % of crystals, with a Poisson's ratio of 0.23 (Chen et al., 2004), with a density of 1.39 g/cm$^3$. Glass transition temperature, in contrast with previous mentions, was not different from that of crystalline cellulose, being in the range of 450–500 K. According to Chen et al. (2004), stress-strain curves showed yielding around 8 % strain alongside dropping intermolecular H-bonds. Comparisons to experiments are challenging at this scale.

Nevertheless, dislocations and amorphous cellulose seems to follow a slightly different stress-strain trend compared to crystals of cellulose. In particular, cellulose crystals show a brittle behavior at failure (Gupta et al., 2022), while disordered cellulose chains allow for gradual changes in the structure and hence, a more ductile behavior.

Not only the mechanical properties of disordered cellulose has shown to expectedly be different, but amorphous cellulose also shows different thermophysical properties (Bregado et al., 2020). Glass transition temperature increased slightly with chain length in disordered cellulose, compared to crystalline structures. $T_g$ also slightly decreased with decreasing the molar mass. Temperature affected density, decreasing it regardless of force field, and increased thermal expansion and isothermal compressibility, especially for shorter chains. The authors concluded that while CHARMM36 shows to be capable of reproducing some of the thermophysical properties of disordered cellulose,
estimations of the heat capacity are not as accurate as they must be; hence, a re-parametrization of the force field was suggested for specific purposes. Overall, the authors claimed that the properties studied in their work with CHARMM36 better represent generic experimental data than those obtained by GLYCAM06 and COMPASS.

In summary, methodologies for generating disordered cellulose involve disrupting crystallinity through techniques like high-temperature Molecular Dynamics simulations with non-reactive force fields. These methodologies are mainly employed to generate dislocated segments in CMFs. The second set of techniques include insertion of chains into simulation boxes, to generate fully random structures, known as amorphous cellulose. Regardless of the exercised methods, studies reveal preferential orientations of hydroxymethyl groups towards gg and gt and disruptions in the lattice parameters (if a crystalline structure is annealed). Investigations emphasize that the transitions in crystallinity mainly take place around the glass temperature of cellulose. Moreover, simulations suggest that short crystalline segments are more difficult to disrupt (Khodayari, Hirn, et al., 2021), while others showed that these short segments have lower T<sub>g</sub> (Bregado et al., 2020). Experimental investigations are required to validate either of these observations. Lastly, mechanical properties, as demonstrated in various studies, highlight lower elastic moduli in disordered cellulose compared to crystalline cellulose.

## 4.2. Interactions with water

One of the other essential features of the disordered phase of cellulose is its interaction mode with water. Many works have considered different properties of amorphous cellulose, which can easier be altered with water due to the permeability of this state of cellulose in contrast with impermeable crystalline cellulose (Kulasinski, Guyer, Derome, & Carmeliet, 2015b). Simulation results show that at low MCs, the water stays in a bound state with no diffusivity, while at saturation point, the diffusivity increases by two orders of magnitude compared to the dry state.

(Kulasinski, Keten, Churakov, Guyer, et al., 2014) modelled absorption of water in amorphous cellulose, and investigated properties of the moist material. To model hydration, water molecules were randomly added to an amorphous structure up to an MC of 50 %. Noticeable changes in the properties were observed at 10 % MC due to merging of water clusters and a substantial decrease in elastic modulus accompanied by increased water diffusion. The change in diffusivity imposed pressure on the amorphous cellulose leading to excessive swelling.

In another work, Kulasinski, Guyer, Keten, et al. (2015) used the same simulation settings to model amorphous cellulose and hemi-celluloses and inspect hydration effects on mechanical properties. They observed a linear relationship between porosity and swelling and concluded that the main reason controlling mechanical properties is hydrogen bond failure, rather than chemical composition. The exponential decay of the bulk and shear stiffness, for instance, was shown to follow the same trend observed for the H-bond failure.

The investigations on non-crystalline cellulose reveal its intricate interactions with water, showcasing permeability and moisture-induced phenomena, such as remarkable increase in water diffusivity. The literature uncovers well the influence of water content on mechanical properties, with a shift happening at 10 % of moisture content. It is worth noting that despite structural changes in disordered cellulose, and specifically on its surface, in composite materials, cellulose surface disorder enhances adhesion energy, emphasizing its potential benefits (Ren et al., 2021; Xia et al., 2018). This can, therefore, be considered as a potential strategy to process composites with enhanced properties.

# 5. Hemicelluloses

Hemicelluloses, abundant substances in the plant cell walls (Rao, Lv, Chen, & Peng, 2023), encompasses diverse types with distinct backbone

structures and side-chain substitutions, varying significantly depending on the plant species, cell type, and cell wall composition. Hemicelluloses are heterogeneous polysaccharides with β-1,4 linked backbones in mannans, xylans, and xyloglucans, and consisted of β-1,4 and β-1,3 linked backbones in mixed-linkage glucans. Hemicelluloses are mainly shown to be synthesized through glycosyltransferases in the Golgi membranes (Scheller & Ulvskov, 2010).

Mannans, prominent in the secondary cell wall of softwoods (Sorieul, Dickson, Hill, & Pearson, 2016), feature β-1,4 linked d-mannopyranosyl and d-glucopyranosyl units within the backbone (glucomannan), occasionally substituted with acetyl groups at C2 or C3 (Hannuksela & Du Penhoat, 2004). Galactoglucomannan (GGM), an extension of this family, includes β-d-galactopyranosyl substitutions on the C6 position of mannose units (Eichinger, Rahkila, Willför, & Xu, 2019) (Fig. 14a).

Xylans, another hemicellulose family, abundant in the cell wall of grasses and hardwood (Sorieul et al., 2016), consist of β-1,4 linked d-xylopyranosyl units in the backbone. Occasional substitutions of d-glucuronic acids on the O2 of xylose units form glucuronoxylans (GX) (de Carvalho et al., 2019; Heinen, Betini, & Polizeli, 2019), often acetylated at C2 and/or C3 (Fig. 14b). Glucuronoarabinoxylans (GAX) have fewer acetyl substitutions, with arabinose substitutions occurring on C2 and/or C3 (Harris, 2006; Kozlova, Mikshina, & Gorshkova, 2012; Vogel, 2008) (Fig. 14c). Ferulic acid esters are also observed as substitutions on the C5 of arabinofuranosyl units in grass cell walls (Scheller & Ulvskov, 2010).

Xyloglucan (XyG), most abundant hemicellulose in the primary cell walls of grasses, hardwoods, and softwoods, features backbones mainly composed of d-glucopyranosylglucose units with α-D-xylopyranosylse substitutions (α-1,6 on C6 of glucose backbone) and bare unsubstituted ones (Fig. 14d). Xylose substitutions occasionally include α-L-fucopyranosyl-substituted β-D-galactopyranosyl units at the C2 (β-1,2 linked) and unsubstituted galactose units, or α-L-arabinofuranosyl-substitutions (Lerouxel, Cavalier, Liepman, & Keegstra, 2006; Pauly et al., 2001; Picard et al., 2000). Galactose units are occasionally acetylated (Pauly et al., 2001). More in depth studies on the variation of the side chains have been performed in the literature, pointing out 17 different substitutions/branches, composed of D-galactose, D-xylose, L-arabinopyranose, L-arabinofuranose, D-galacturonic acid, L-fucose, and L-galactose, which will not be discussed further in this review (Schultink, Liu, Zhu, & Pauly, 2014). Atomistic studies on this category of hemi-celluloses dates back to the 90's (Levy, Maclachlan, & Staehelin, 1997).

Another major hemicellulose type in the primary cell walls of grasses is mixed-linkage glucans (MLG). MLGs (Fig. 14e), with β-1,4 and β-1,3 linked glucosyl residue backbones, are abundant in Equisetum and Poales (Fry, Nesselrode, Miller, & Mewburn, 2008; Sørensen et al., 2008). MLGs are linear, unbranched, and unsubstituted (Purushotham et al., 2022), as shown in Poalean, for instance (Fry et al., 2008).

There are reports on the biological functions of the hemicelluloses in the plant cell walls, yet, a complete model illustrating the structure and organization of the components in the plant cell walls is yet to be presented. Nevertheless, proposals for such arrangements are presented in the last two decades (Cosgrove, 2005, 2014, 2022). XyG binds to cellulose through hydrogen bonds, or entrapped between microfibrils during synthesis. The binding is believed to happen in many conformations including extended or coiled (Cosgrove, 2022). The interactions between XyG and cellulose assert that XyGs are not the major load bearing elements of the plant cell walls (Thompson, 2005). On the contrary, in grasses for instance, xylans have an vital role in strengthening the secondary cell walls (Scheller & Ulvskov, 2010), while recent research undermines the concept of considering them as the major load bearing tether between CMFs (Cosgrove, 2022).

## 5.1. Hemicelluloses composition and flexibility

Considerable attention is devoted to studying hemicelluloses flexibility, particularly in relation to the ratio and position of substitutions on
<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Examples for a) Galactoglucomannan (GGM), b) Glucuronoxylan (GX), c) Glucuronoarabinoxylan (GAX), d) Xyloglucan (XyG), and e) mixed-linkage glucan (MLG). It must be noted that the structures depicted above are examples of possible chemical compositions, and the arrangement of the backbone, type and number of substitutions, as well as the degree of polymerization differs in different plants. In this figure β-D-glucopyranosyl units are shown in army green, β-D-mannopyranosyl units in black, β-D-galactopyranosyl units in maroon red, β and α-D-xylopyranosyl units in cyan, α-D-glucuronic acid units in orange, α-L-fucopyranosyl units in deep pink, α-L-arabinofuranosyl units in blue, acetyl units in red, and methyl units in purple. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
their backbone, which governs their binding mechanism to cellulose (Chen, Cathala, & Lahaye, 2022). Hemicelluloses are considered rigid, if the polymeric chains show a twofold (2₁) helical screw (360° of twist per two residues), and flexible, if they show a threefold (3₁) helical screw (360° of twist per three residues). From a geometrical point of view, the summation of the two dihedral angles about the glycosidic angle (ψ and φ) becomes ~120° in a 2₁ conformation, and ~50° or ~190° in a 3₁ conformation (Addison et al., 2024; Khodayari, Thielemans, et al., 2021). Considering the 2₁ conformation of cellulose chains in CMFs, hemicellulose models with a rigid backbone can better adsorb onto cellulose. On the contrary, adsorption of grass xylans, for instance, on the disordered cellulose regions has been shown to occur as weak interactions dictated by their 31 conformation through vdW and occasional H-bonds (Gao, Lipton, Wittmer, Murray, & Mortimer, 2020). The chemical composition of the backbone and substitutions type of the hemicelluloses is of utmost importance in defining their flexibility (Berglund et al., 2016).

### 5.1.1. Mannans

Mannans with different acetyl substitutions on C3, C2, or both carbons are studied extensively in the literature. Structure exhibit greater flexibility when substitution occurred on C3, or both C2 and C3 on the same sugar. Flexibility diminished with substitution only on C2. This stiffer conformation improved adsorption onto cellulose, correlated well with their experiments.

In another study, Yu et al. (2018) modelled parallel and antiparallel modes,<fnref n="3" /> of glucomannan structured of alternating glucose and mannose backbones with galactose substitutions on the mannose units. Hemicellulose models displayed 2₁ conformations despite hydrophilicity of the surface. Binding between cellulose and hemicellulose models was mentioned to occur with the C2 of mannose units pointing away from cellulose's hydrophobic surface. Patterned mannosyl substitutions diminishes a steric hinderance, acting as an advantage in binding to the hydrophobic face of cellulose.

Martínez-Abad et al. (2020) also modelled a CNC with various hemicelluloses including GGGG, MMMM, MGMG, MMGM and MLMM (where G stands for Glc, M for Man, and L for galactose-substituted Man) to analyze flexibility effects. The GGGG model was rigid with a 2₁ conformation, while the MMMM model was more flexible. The MGMG model showed similar stiffness to GGGG. Umbrella sampling showed GGGG and MGMG models had twice the interaction energies with cellulose versus the MMMM model. MMGM and MLMM also had higher energies than MMMM but not as strong to GGGG and MGMG. The conclusion is that backbone composition significantly affects cellulose affinity (as also seen for branching composition in xylans).

The axial C2 hydroxyl configuration in mannans hinders glucose-like shapes and rigidity near cellulose. Hence, when GGGG and MMMM models were compared, less rigidity was observed for the MMMM model. This issue is recovered by even mannose spacing with C2 hydroxyls pointing outward (as in MGMG). Consecutive mannosyl units drastically decrease adsorption energies. Substitutions on C2 likely increase mannan and xylan stiffness. Despite lower affinity than glucans, xylose lacking C6 has less impact on interaction energies than mannose's C2 hydroxyl.

Berglund et al. (2018) also simulated glucomannan models differing in substitutions and glucose content in alkaline conditions. Linkages were slightly stiffer with glucose in the backbone, and even stiffer by galactose substitutions. The authors also check the interactions between the C2 hydroxyl and the alkaline, as hydrolysis is assumed to happen due to the deprotonation of this group. OH- ions were added to the simulation boxes to model alkalinity. Glucose increased OH- coordination

around OH2 and OH3 versus mannose, attributed to the equatorial versus axial hydroxyl orientations in glucose and mannose, respectively, hinting at faster hydrolysis.

In general, computational strategies have confirmed the role of Mannans' C2-substitution on the rigidity of the polymeric chains, owing to the orientation of the C2 hydroxyls. Despite all, while flexibility plays a role in cellulose binding, interaction energies might depend more on the chemical composition (Khodayari, Thielemans, et al., 2021). In particular, it is shown that while the flexibility of hemicelluloses plays a role in its binding to cellulose, the magnitude of interaction energies is defined based on the chemical composition of the structures, and hence, more attention must be dedicated to that aspect of the hemicelluloses.

### 5.1.2. Xylans

Xylans are shown to primarily bind to cellulose in a twofold helical screw manner (Simmons et al., 2016). As an example, ssNMR results confirm that architecture of softwood includes xylan which exists in 2₁ conformation, indicating a bound state to cellulose surface (Terrett et al., 2019). The role of GAX side-chains on adsorption to cellulose was inspected in 2017 (Pereira et al., 2017), concluding that position rather than composition determines interactions, with α-1,2 links stabilizing the GAX on cellulose, similar to the C2 substitutions for mannans (Berglund et al., 2020).

Other work modelled cellulose Iβ and two GAX variations (Shrestha et al., 2019). The highly C3-substituted model had a lower radius of gyration, potentially facilitating self-aggregation versus the stiffer, less substituted model, which exhibited better cellulose affinity. Though highly substituted GAX had more overall cellulose contacts, weaker interactions were concluded due to higher mobility. Importantly, models differed in substitution position, affecting stiffness comparisons, as pointed out by Berglund et al. (2020) and Pereira et al. (2017). The stiffness of the low substituted model could stem from the substitution on the C2, and the flexibility of the highly substituted one from C3 substitution. Moreover, interaction magnitudes require energetic measurements.

Backbone substitution effects on xylan cellulose affinity were studied using xylohexaose models (Martínez-Abad et al., 2017). Free energy calculations showed substitutions increased interaction versus unsubstituted xylan. A doubly substituted model with glucuronic acids showed highest interaction energy via side chain cellulose surface contacts. Interactions were higher on the (200) face compared to (1-10), likely due to fewer water hindering interactions. Detachment progressed from ends distant to substitutions, confirming their positive effect.

Acetylation patterns on xylan folding and cellulose affinity were investigated using cellulose and xylan models (Busse-Wicher et al., 2014). Unsubstituted and evenly glucuronic acid/acetyl substituted xylans maintained 2₁ conformation and similar interaction energies (10-15 kcal/mol per xylose unit) on hydrophilic faces. Acetyls decreased (100) face interactions while glucuronic acid-substituted models showed improved interactions on the (200) face. Acetylations had minor water interaction effects versus glucuronic acid substitutions. Spread out substitutions maintained cellulose hydrogen bonding, should the substitutions point outwards with respect to the cellulose hydroxyl groups. These observation are similar to what also observed by Gupta et al. (2023). C2-acetylation promoted 2₁ conformation beneficial for cellulose versus lignin affinity needing 3₁ conformations from absence of acetylation, C3 acetylation, or full substitution.

Additional work modelled cellulose crystals with xylan variations, finding acetylation rigidified xylan (Gupta et al., 2021). The C2 acetylation specifically led to stiffness near cellulose over 3₁ conformations in solvents. Surrounding C2 acetylations on a specific residue provided further rigidity. Usubstituted and glucuronic acid substituted xylans maintained 3₁ conformations near cellulose. C3 acetylation showed partial 2₁ conformations but overall 3₁ structure, confirming less affinity for cellulose. Nevertheless, these studies suggest it is the patterned substitutions (evenly spaced, and mainly on C2) that are essential for

<footnote n="3">³ Parallel: when both structures showed reducing ends on the same side; antiparallel: when reducing end of hemicellulose faced the non-reducing end of cellulose</footnote>
providing a 2₁ conformation of xylans, facilitating better docking onto the hydrophilic faces of cellulose (Grantham et al., 2017). These patterned structures dictate the acetyl groups to point outwards when xylan is adsorbed onto cellulose (Addison et al., 2024).

Other modeling found (010) and (020) cellulose grooves stabilized GAX 2₁ conformations on these faces (Busse-Wicher et al., 2016). Single GAX binding fluctuations exhibited temporary threefold structures, confirming cellulose interactions also stabilize GAX to 2₁ conformations. Presence of multiple bound GAX structures showed greater dihedral angle stability and outward substitutions. However, side chain effects in confined spaces of two adjacent CMFs needs further inspection.

The studies above suggest that the flexibility of hemicellulose models is influenced by their chemical composition and the position of backbone substitutions. Notably, substitutions at C2 significantly impact the rigidity of hemicellulose structures, particularly in those with full mannose backbones. Given mannose's epimeric relationship with glucose, the orientation of the C2 hydroxyl imparts greater flexibility, resulting in lower affinity for cellulose. It remains essential to quantitatively measure and compare the effects of substitution position, ratio, type, and backbone chemical composition on the interaction energies of hemicellulose structures with cellulose for robust conclusions.

## 5.2. Hemicelluloses and water

Interactions of hemicelluloses and water are important in defining the mechanical properties of plant fibres, given hemicellulose's role as a hydrophilic cell wall component. Its flexibility can significantly change with moisture content and composition. Moreover, H-bonds at cellulose-hemicellulose interfaces can quickly deteriorate, directly affecting mechanical properties.

A recent study simulated GAX models to inspect moisture effects (Zhang et al., 2020). Heat of adsorption, heat capacity, thermal expansion, and elastic modulus shifted trends at 30 % moisture content, despite GAX-GAX and GAX-solvent H-bond saturation at 18 %. This shows water forms two layers on GAX-bound and unbound. The bound layer, starts to form from the beginning and saturates at 30 % moisture, while the unbound layer only starts emerging at 18 % of moisture, continuing to grow even after 30 %. No further GAX-water contacts occurred after 30 %, indicating added water attaches to the first layer rather than onto GAX. The property changes were thus attributed to the two water layers. Further confirming simulations, like varying GAX degree of polymerization and molar mass, could be insightful.

Moisture also distinctly affects hydrophobic lignin versus hydrophilic hemicelluloses in their coexistence within fibres. Lignin and GAX structures were modelled to study moisture effects on system elastic moduli (Youssefian et al., 2017). Experiments on bamboo fibres showed increased longitudinal modulus up to 3 % moisture, then a decrease of both transverse and longitudinal moduli at higher MC. Simulations assessed this - lignin's elastic modulus increased at low moisture then decreased at higher levels, while GAX's depleted with moisture. Hence, transverse modulus was deemed governed by GAX, while longitudinal by lignin. Mechanisms were explained as such; at low moisture, hydrogen bridges between lignin chains decrease free volume, increasing modulus. For GAX, water causes swelling and free volume gains, decreasing modulus. At higher moistures, nano-droplets form in both, lowering elastic moduli.

Additionally, cellulose-hemicellulose-lignin complexes were simulated to determine hydration effects (Cresswell et al., 2021). Hemicelluloses (including GAX and GGM) adsorbed onto cellulose's hydrophilic face assumed a 2₁ conformation when hydrated. Drying collapsed the complex into a dense, shrunk structure, while hemicelluloses interacted more strongly with cellulose, drawing surrounding lignin closer. Rehydration increased cellulose-hemicelluloses and lignin distances, though interfacial water did not fully recover at the interface of cellulose and hemicelluloses. Additionally, it was shown that GAX interacts stronger with cellulose than acetylated-GGM upon

rehydration, indicating their differing cell wall functionalities.

The interaction between hemicelluloses and water significantly influences the mechanical properties of plant fibres. The decrease in the mechanical properties of the cell wall with moisture, can partly be explained by positioning of the water molecules at the interface of cellulose and hemicellulose, thanks to the hydrophilic nature of the two structures. Hemicellulose flexibility, responsive to moisture and composition, undergoes shifts under varying conditions. Studies on hemicellulose-lignin complexes have shown diverse responses to moisture, being an increase in the mechanical properties at low moisture contents followed by abrupt decreases from 3 % onwards up to a saturation point at 30 %. Further investigations into cellulose-hemicellulose-lignin complexes are necessary to elucidate the dynamic structural changes in plant fibres during hydration cycles, as, so far, not many studies considering a full model of the plant cell wall, i.e. including cellulose, hemicelluloses, lignin, pectin, and water is investigated in the literature.

## 5.3. Cellulose and hemicelluloses interactions

Interactions between cellulose and hemicelluloses have been a major research focus when modeling plant fibres at the nanoscale. As the second most abundant component in plant cell walls after cellulose, hemicelluloses act as a hydrophilic binding material whose properties are greatly affected by moisture. Understanding the binding mechanisms and energies between cellulose and hemicelluloses, as well as the mediating role of water, is thus important. As a starting point, studying the adsorption energies of cellulose and hemicellulose building blocks can provide insight into the affinity between these substances.

Umbrella sampling simulations have been used to calculate interaction energies between pairs of cellobiose molecules, xylobiose molecules, and cellobiose-xylobiose pairs (Peri et al., 2011). The pairs were simulated either freely interacting or constrained in a defined stacked position. All models showed a single potential well, deepest for the cellobiose pair (1.84 kcal/mol) followed by the cellobiose-xylobiose pair (1.2 kcal/mol) and xylobiose pair (1 kcal/mol). This suggests cellulose building blocks form more robust interactions with each other than with hemicelluloses, although cellulose-hemicellulose interactions are still greater than intermolecular interactions within the hemicellulose matrices. Constraining the molecules into crystalline stacked series increased interaction energies and decreased separation distances, indicating crystallization strengthens these interactions.

It is believed hydrogen bonding is of the controlling interactions between cellulose and hemicellulose (Khodayari, Thielemans, et al., 2021), though the specific binding modes in plant cell walls are still unclear. Different bonding modes for GX-cellulose complexes, i.e. - bridging, looping, or random, were proposed and studied (Zhang et al., 2015). Simulation results suggested that the bridging mode required the highest force to shear the CNCs apart, while the random alignment was weakest. However, the singular shear rate applied may have biased these findings. The contact area and corresponding hydrogen bonding were deemed responsible for the system's shear strength, with water acting as a lubricant to soften the GX matrix.

Different hemicellulose types can have different interaction energies with cellulose (Khodayari, Thielemans, et al., 2021). It is known that hemicellulose adsorption restricts the dynamics of cellulose surface glucose units. Xylans tend to align both parallel and antiparallel, with a higher probability for the latter, when adhering to cellulose fibrils (Heinonen et al., 2022), though this orientation may depend on polymer chain lengths, decoration patterns, moisture content, and cellulose face. Therefore, studies observing antiparallel orientation of xylans would not reject parallel adsorption (Kong et al., 2022).

Some studies have shown xylan acetylation patterns influence the supramolecular organization and hydration of xylan-cellulose complexes, with deacetylated segments more hydrated and acetylated fragments more rigid (Jaafar et al., 2019). This confirms that
acetylations could potentially lead to a more rigid conformation, but less
hydrated structure, mediating cellulose-hemicellulose interactions. In
contrast, Mazeau and Charlier (2012) modelled bare and substituted 5-
unit xylans interacting with crystalline cellulose. All xylans showed the
same affinity for cellulose, leading them to conclude substitutions have a
negligible effect. However, bare xylans perfectly maintained alignment
along the cellulose axis compared to the substituted variants.

Simulations by *Zhao et al.* (2014) found xyloglucan side-chain sub-
stitutions and cellulose surface topography and hydrophobicity signifi-
cantly impact xyloglucan’s conformational and binding behavior. For
example, xyloglucan binds stronger to hydrophobic cellulose faces in a
flat conformation compared to the weaker random coil structures it
forms on hydrophilic faces. Side-chains on the xyloglucan deteriorated
its flat conformation on the hydrophobic surface of cellulose, leading to
less binding strength. However, xyloglucan affinity for the hydrophilic
surface of cellulose remained intact with/without substitutions.

In 2011, Zhang, Brumer, et al. (2011) modelled the (1-10) face of cellulose I/β alongside three 9-monomer xyloglucans with different side-chains. The interaction energies was computed based on the work of adhesion (Jarzynski, 1997), showing minor effects of the galactose and fucose substitutions compared to unsubstituted xyloglucan (3.4-3.8 kcal/mol/glucose). The adsorption free energies were higher than those recently obtained by Khodayari, Thielemans, et al. (2021) for 8-monomer xyloglucans.

Hanus and Mazeau (2006) studied adsorption of different glucose-, galactose-, or fucogalactose-substitutions xyloglucan models onto multiple cellulose faces. Results showed no orientation preference (parallel, anti-parallel, or perpendicular) for xyloglucan adsorption, unlike xylans (Mazeau & Charlier, 2012). vdW interactions were weaker than electrostatic interactions. Considering H-bond energy of 5–15 kcal/mol, the contribution of the H-bonds was 10–30 %. Longer side-chains increased interaction strength but reduced bonding efficiency due to higher mobility that can trigger desorption. Maximum interaction energy per residue was 23 kcal/mol for a DP10 xyloglucan.

Kishani et al. (2021) studied cellulose and xyloglucan interactions with native and TEMPO-oxidized cellulose using umbrella sampling. They found that surface charges on cellulose lowered interaction energies with xyloglucan compared to native cellulose. Interestingly, they showed that at very low temperature of 100 K, the value of the δG is positive, which means that spontaneous desorption takes place. This indicates an endothermic sorption process driven by favorable entropy due to release of bound water at the interface. Reported interaction energies were consistent with values by Khodayari, Thielemans, et al. (2021) but lower than estimates by Mazeau and Charlier (2012). According to Mazeau and Charlier (2012), the interaction energies between xyloglucan and cellulose is estimated at 10 kcal/mol/sugar, which should provide an energy of ~ 330 kJ/mol for Kishani et al. (2021)'s models. Nevertheless, the cellulose models in this study were TEMPO-oxidized cellulose which could possibly affect the magnitude of the interaction energies.

Kumar et al. (2018) studied xylan and cellulose interactions at different temperatures, finding more contacts at high temperatures due to changes in water solvation shells depleting around cellulose, facilitating cellulose-xylan interactions. The effects were largely maintained when the structures were cooled down to room temperature.

While hemicellulose content is known as a binding substance be-
tween CMFs, lignin fills out the voids, acting as a binding filler in the
plant cell wall. Charlier and Mazeau (2012) modelled a plant cell wall
section with cellulose, xylan, lignin and water. Xylans with a density of
1.5 g/cm³ laid flat on cellulose, penetrating into lignin, facilitated by
chemical bonds. Lignin showed varying alignment to cellulose and
densities from 0.5 to 2 g/cm³. Despite these differences, all lignin
models showed an amorphous structure and accommodated water
molecules equal to those in the xylan content. However, the diffusion
coefficients computed for the water molecules in the lignin content were
50 % higher than in the xylan fraction. Hence, lignin seems to be holding

more free water compared to other constituents, while more bound
water seems to locate at the interface of cellulose and hemicelluloses.

In 2015, Youssefian and Rahbar (2015) modelled bamboo fibrils containing cellulose, amorphous cellulose, GAX and lignin. The mechanical properties of the fibrils were mainly attributed to the GAX content, while lignin was shown to have the effect of enhancing the adhesion between the amorphous matrix and cellulose fibrils. GAX hydrogen bonding and lignin vdW interactions enhanced matrix adhesion and mechanical properties. The amorphous cellulose region was the weakest link determining composite strength. A more complete model for Populus secondary cell wall was investigated through a combined ssNMR-MD study (Addison et al., 2024). Results of several different models indicated that an aggregated structure of 8 CMFs where xylan was partially trapped between the fibrils, covering the whole exterior surface of the aggregate, and mediating the direct lignin-cellulose interactions was the best fit to the ssNMR results.

In summary, the exploration of cellulose-hemicellulose interactions
reveals crucial insights into the intricate dynamics shaping plant cell
wall architecture. Many studies underscore the nuanced effects of side
chains in hemicelluloses, exemplified by various types such as xylans.
The investigation into flexibility, particularly in xyloglucans, empha-
sizes the role of surface topography in dictating structural and binding
mechanisms. Temperature emerges as a critical factor, showcasing
amplified cellulose-hemicellulose contacts at elevated temperatures,
attributed to changes in water solvation shell properties. Interestingly,
comprehensive modeling incorporating lignin unveils its role in
enhancing adhesion between the amorphous matrix and cellulose fibrils.

6. Conclusions

The current literature extensively explores the advances made in the
community of computational carbohydrate polymers and materials
science, in particular molecular dynamics simulations, to model the
plant cell wall components, investigate their properties, and their in-
teractions. The potentials and limitations of these methods are exten-
sively discussed and compared. While the plant cell walls consist of
complex interplay between different components, including cellulose (in
crystalline and disordered morphologies), hemicelluloses, pectin, lignin,
and water, this review maintained focus on only cellulose and some
hemicelluloses, with minor references to lignin. It must be noted that
this review does not provide a complete understanding of the plant cell
wall, but rather aims at advancing towards it. To gain such a compre-
hensive perspective would require further study of other key compo-
nents like lignin and pectin, as well as computational and experimental
analysis of their synergistic roles in the plant cell walls. An in-depth
examination of these factors was beyond the scope of the current
work. However, elucidating their roles is essential for fully under-
standing the plant cell walls.

It must also be considered that full models of the plant cell walls
would require greater simulation scales than classical molecular dy-
namics can currently offer. Computational research earnestly seeks to
bridge the gap between atomistic studies and large-scale experiments,
while experimental techniques eagerly advance to higher precision
findings, pointing towards a future of enhanced understanding.

While research has had great advances through proposing a complete
model for the plant cell walls, literature still lacks a rigid understanding.
In computational studies, building the initial geometries of the models
for the plant cell walls is of utmost importance. Recent studies show that
conformation of the polysaccharides, as well as their distribution and
adsorption procedure can be heavily affected by the initial conditions of
the simulations. Geometries deviating from accurate plant cell wall
models could possibly lead to inaccurate conclusions derived from
atomistic simulations.

Over the years, force fields have been developed, enabling precise
modeling of carbohydrates and allowing the study of these materials in
combination with biomaterials, proteins, etc. Comparisons between the
force fields and their capabilities were discussed in this review. Current progress in providing accurate force fields to model carbohydrates is satisfactory, yet, the community desires faster and precise reactive force fields, as well as the development of coarse-grained potentials able to consider many components of the plant cell walls and their properties.

Our current understanding of cellulose, its structure and function in the plant cell walls seems to be well established. Hydrogen bonding patterns and cellulose twist have been extensively studied, and many findings have been reviewed in this work. Molecular dynamics simulations have facilitated inspection of cellulose different forms, in terms of mechanical and thermal properties with great details, elucidating the mechanisms and transformations at various hygrothermal conditions. Joint computational and experimental realms have contributed a lot to this advancement, yet, consensus on the number and arrangement of cellulose chains in CMFs is missing. This has caused utilization of several different models for the CMFs in the MD simulations in the last three decades, and required much more attention. While the employed models are versatile, comparisons between the findings illustrate that the obtained outcomes share a lot of similarities. Nevertheless, current research seems to mainly agree on the 18-chain models, and the possibility of generating 36-chain structures through aggregation of these 18-chain building blocks.

The review further explores the studies of cellulose-cellulose interactions, demonstrating the prominent role of vdW and H-bonds in controlling the dynamics. The impact of water on cellulose is subsequently reviewed in great detail. It focuses on many aspects such as the hydration of different faces of cellulose, different regimes of confined and non-confined water, comparisons between water models in studying cellulose, and hydration effect on surface-modified cellulose. The choice of water model, the employed force field, simulation parameters, etc. have shown to alter the outcome of the findings.

The review covers the studies on disordered cellulose. The methodologies to model amorphous and dislocated cellulose have been introduced in-depth, and comparisons were made where relevant. Apart from the techniques and their characteristics, molecular dynamics simulations have enabled extensive inspection of the properties of disordered cellulose. However, additional experimental inspections are required to clarify remaining aspects of these substances, such as dislocations' length within the CMFs.

Hemicelluloses have been studied thoroughly in the literature from both experimental and computational perspectives. These joint studies have significantly enhanced the understanding within the community of carbohydrates. Much attention has been dedicated to the flexibility of the hemicelluloses, a crucial parameter in defining their adsorption capacities to cellulose, and their hydration properties. Current research is well established on the effect of substitutions and their positions on hemicelluloses flexibility, such as the periodic substitutions in xylans, or C2-substitutions on mannose units in GGMs. Furthermore, computational studies have focused on calculating the interaction energies between cellulose and hemicelluloses, although this area requires further attention.

Lastly, the review concludes by covering studies on the interactions between cellulose, hemicelluloses, and lignin. Molecular dynamics simulations have made significant advances in providing a detailed view of the mechanisms governing these interactions. We believe that the development of more accurate force fields and models, combined with advances in computational power, will undoubtedly lead to a better understanding of these complex systems. Atomistic modeling provides a powerful tool for investigating the properties of natural fibres, and its continued development will further enhance our understanding of these materials, opening up new avenues for their use in various applications.

## CRediT authorship contribution statement

**Ali Khodayari:** Writing – review & editing, Writing – original draft, Methodology, Investigation, Data curation, Conceptualization.

**Ulrich Hirn:** Writing – review & editing, Writing – original draft. **Stefan Spirk:** Writing – review & editing, Writing – original draft. **Yu Ogawa:** Writing – review & editing, Writing – original draft. **David Seveno:** Writing – review & editing, Writing – original draft. **Wim Thielemans:** Writing – review & editing, Writing – original draft.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

All data that support this study are kept available at KU Leuven and can be accessed upon request from the corresponding author.

## Acknowledgement

The authors would like to express their sincere gratitude to Dr. Muhammad Farooq for his invaluable contribution and insightful inputs. The authors would like to acknowledge the European Union's Horizon 2020 research and innovation programme and project FibreNet, under Marie Skłodowska-Curie grant agreement No 764713. The authors would also like to acknowledge project M-Era Net 16531229 BioCel3D and project Core 16224491 BIOINTER funded by FNR Luxembourg.

## References

Addison, B., Bu, L., Bharadwaj, V., Crowley, M. F., Harman-Ware, A. E., Crowley, M. F., ... Cieslinski, P. N. (2024). Atomistic, macromolecular model of the Populus secondary cell wall informed by solid-state NMR. *Science Advances*, 10(1), https://doi.org/10.1126/sciadv.adi7965

Adler, D. C., & Buehler, M. J. (2013). Mesoscale mechanics of wood cell walls under axial strain. *Soft Matter*, 9(29), 7138–7144. https://doi.org/10.1039/c3sm50183c

Aguilera-Segura, S. M., Di Renzo, F., & Mineva, T. (2020). Molecular insight into the cosolvent effect on lignin-cellulose adhesion. *Langmuir*, 36(47), 14403–14416. https://doi.org/10.1021/acs.langmuir.0c02794

Aktulga, H. M., Pandit, S. A., Van Duin, A. C., & Grama, A. Y. (2012). Reactive molecular dynamics: Numerical methods and algorithmic techniques. *SIAM Journal on Scientific Computing*, 34(1), 1–23.

Altaner, C. M., Thomas, L. H., Fernandes, A. N., & Jarvis, M. C. (2014). How cellulose stretches: Synergism between covalent and hydrogen bonding. *Biomacromolecules*, 15(3), 791–798.

Asgarpour Khansary, M., Pouresmaeel-Selakjani, P., Aroon, M. A., Hallajisani, A., Cookman, J., & Shirazian, S. (2020). A molecular scale analysis of TEMPO-oxidation of native cellulose molecules. *Heliyon*, 6(12), Article e05776. https://doi.org/10.1016/j.heliyon.2020.e05776

Azimzadeh Irani, M., Askari, H., Jahanfar, M., Nasehi, M., & Hamedi, A. (2023). Molecular arrangement of cellulose bio-nanofibers in formation of higherorder assemblies. *Journal of Polymer Research*, 30(6), 1–15. https://doi.org/10.1007/s10965-023-03573-1

Baird, M. S., O'Sullivan, A. C., & Banks, W. B. (1998). A native cellulose microfibril model. *Cellulose*, 5(2), 89–111. https://doi.org/10.1023/A:1009220830595

Benhassine, M., Saiz, E., Tomsia, A. P., & De Coninck, J. (2009). Nonreactive spreading at high-temperature revisited for metal systems via molecular dynamics. *Langmuir*, 25 (19), 11450–11458. https://doi.org/10.1021/la902958k

Bergenstråhle, M., Berglund, L. A., & Mazeau, K. (2007). Thermal response in crystalline *l/f* cellulose: A molecular dynamics study. *Journal of Physical Chemistry B*, 111(30), 9138–9145. https://doi.org/10.1021/jp072258i

Bergenstråhle, M., Thormann, E., Nordgren, N., & Berglund, L. A. (2009). Force pulling of single cellulose chains at the crystalline cellulose - Liquid interface: A molecular dynamics study. *Langmuir*, 25(8), 4635–4642. https://doi.org/10.1021/la803915c

Bergenstråhle, M., Wohlert, J., Himmel, M. E., & Brady, J. W. (2010). Simulation studies of the insolubility of cellulose. *Carbohydrate Research*, 345(14), 2060–2066. https://doi.org/10.1016/j.carres.2010.06.017

Berglund, J., Angles d'Ortoli, T., Vilaplana, F., Widmalm, G., Bergenstråhle-Wohlert, M., Lawoko, M., ... Wohlert, J. (2016). A molecular dynamics study of the effect of glycosidic linkage type in the hemicellulose backbone on the molecular chain flexibility. *Plant Journal*, 88(1), 56–70. https://doi.org/10.1111/tpj.13259

Berglund, J., Azhar, S., Lawoko, M., Lindström, M., Vilaplana, F., Wohlert, J., & Henriksson, G. (2018). The structure of galactoglucomannan impacts the degradation under alkaline conditions. *Cellulose*, 1–21. https://doi.org/10.1007/s10570-018-1737-z

Berglund, J., Kishani, S., Morais De Carvalho, D., Lawoko, M., Wohlert, J., Henriksson, G., ... Vilaplana, F. (2020). Acetylation and sugar composition influence the (in)solubility of plant β-mannans and their interaction with cellulose surfaces.
ACS Sustainable Chemistry and Engineering, 8(27), 10027–10040. https://doi.org/10.1021/acssuschemeng.0c01716

Bering, E., Torstensen, J., Lervik, A., & de Wijn, A. S. (2022). Computational study of the dissolution of cellulose into single chains: The role of the solvent and agitation. *Cellulose*, **29**(3), 1365–1380. https://doi.org/10.1007/s10570-021-04382-9

Besombes, S., & Mazeau, K. (2005a). The cellulose/lignin assembly assessed by molecular modeling. Part 1: Adsorption of a threo guaiacol β/β-0 dimer onto a 1β/1 cellulose whisker. *Plant Physiology and Biochemistry*, **43**(3), 299–308. https://doi.org/10.1016/j.plaphy.2005.02.005

Besombes, S., & Mazeau, K. (2005b). The cellulose/lignin assembly assessed by molecular modeling. Part 2: Seeking for evidence of organization of lignin molecules at the interface with cellulose. *Plant Physiology and Biochemistry*, **43**(3), 277–286. https://doi.org/10.1016/j.plaphy.2005.02.004

Bourmaud, A., Beaugrand, J., Shah, D. U., Placet, V., & Bailey, C. (2018). Towards the design of high-performance plant fibre composites. *Progress in Materials Science*, **97** (May), 347–408. https://doi.org/10.1016/j.pmatsci.2018.05.005

Bregado, J. L., Secchi, A. R., Tavares, F. W., de Sousa Rodrigues, D., & Gambetta, R. (2019). Amorphous paracrystalline structures from native crystalline cellulose: A molecular dynamics protocol. *Fluid Phase Equilibria*, **491**, 56–76. https://doi.org/10.1016/j.fluid.2019.03.011

Bregado, J. L., Tavares, F. W., Secchi, A. R., & Segtovich, I. S. V. (2020). Thermophysical properties of amorphous-paracrystalline celluloses by molecular dynamics. *Macromolecular Theory and Simulations*, **29**(4), 1–21. https://doi.org/10.1002/mats.202000007

Bregado, J. L., Tavares, F. W., Secchi, A. R., & Segtovich, I. S. V. (2021). Molecular dynamics of dissolution of a 36-chain cellulose 1β microfibril at different temperatures above the critical pressure of water. *Journal of Molecular Liquids*, **336**, Article 116271. https://doi.org/10.1016/j.molliq.2021.116271

Bu, L., Beckham, G. T., Crowley, M. F., Chang, C. H., Matthews, J. F., Bomble, Y. J., ... Nimlos, M. R. (2009). The energy landscape for the interaction of the family 1 carbohydrate-binding module and the cellulose surface is altered by hydrolyzed glycosidic bonds. *Journal of Physical Chemistry B*, **113**(31), 10994–11002. https://doi.org/10.1021/jp904003z

Bu, L., Himmel, M. E., & Crowley, M. F. (2015). The molecular origins of twist in cellulose I-beta. *Carbohydrate Polymers*, **125**, 146–152. https://doi.org/10.1308/rcsann.2016.0321

Busse-Wicher, M., Gomes, T. C., Tryfona, T., Nikolovski, N., Stott, K., Grantham, N. J., ... Dupree, P. (2014). The pattern of xylan acetylation suggests xylan may interact with cellulose microfibrils as a twofold helical screw in the secondary plant cell wall of *Arabidopsis thaliana*. *Plant Journal*, **79**(3), 492–506. https://doi.org/10.1111/tpj.12575

Busse-Wicher, M., Li, A., Silveira, R. L., Pereira, C. S., Tryfona, T., Gomes, T. C. F., ... Dupree, P. (2016). Evolution of xylan substitution patterns in gymnosperms and angiosperms: implications for xylan interaction with cellulose. *Plant Physiology*, **171** (August). https://doi.org/10.1104/pp.16.00539. pp.00539.2016.

Chae, I., Paniagua-Guerra, L. E., Pitcher, M. L., Koshani, R., Yuan, M., Lin, Y. T., ... Kim, S. H. (2023). Relaxation dynamics of water in the vicinity of cellulose nanocrystals. *Cellulose*, **30**(13), 8051–8061. https://doi.org/10.1007/s10570-023-05361-y

Charlier, L., & Mazeau, K. (2012). Molecular modeling of the structural and dynamical properties of secondary plant cell walls: Influence of lignin chemistry. *Journal of Physical Chemistry B*, **116**(14), 4163–4174. https://doi.org/10.1021/jp300395k

Charvati, E., Zhao, L., Wu, L., & Sun, H. (2021). A new parameterization of an all-atom force field for cellulose. *JOM*, **73**(8), 2335–2346. https://doi.org/10.1007/s11837-021-04732-9

Chen, M., Cathala, B., & Lahaye, M. (2022). Adsorption of apple xyloglucan on cellulose nanofiber depends on molecular weight, concentration and building blocks. *Carbohydrate Polymers*, **296**(August), Article 119994. https://doi.org/10.1016/j.carbpol.2022.119994

Chen, M., Zhang, C., & Ke, L. L. (2023). The biomechanical role of the functionally graded microfibrils in the wood cell wall. *Journal of the Mechanics and Physics of Solids*, **175**(March). https://doi.org/10.1016/j.jmps.2023.105296

Chen, P., Lo Re, G., Berglund, L. A., & Wohlert, J. (2020). Surface modification effects on nanocellulose-molecular dynamics simulations using umbrella sampling and computational alchemy. *Journal of Materials Chemistry A*, **8**(44), 23617–23627. https://doi.org/10.1039/d9ta09105g

Chen, P., Nishiyama, Y., & Mazeau, K. (2012). Torsional entropy at the origin of the reversible temperature-induced phase transition of cellulose. *Macromolecules*, **45**(1), 362–368. https://doi.org/10.1021/ma201954z

Chen, P., Nishiyama, Y., & Mazeau, K. (2014). Atomic partial charges and one Lennard-Jones parameter crucial to model cellulose allomorphs. *Cellulose*, **21**(4), 2207–2217. https://doi.org/10.1007/s10570-014-0279-2

Chen, P., Ogawa, Y., Nishiyama, Y., Bergensträhle-Wohlert, M., & Mazeau, K. (2015). Alternative hydrogen bond models of cellulose II and IIII based on molecular force-fields and density functional theory. *Cellulose*, **22**(3), 1485–1493. https://doi.org/10.1007/s10570-015-0589-z

Chen, P., Wohlert, J., Berglund, L., & Furó, I. (2022). Water as an intrinsic structural element in cellulose fibril aggregates. *Journal of Physical Chemistry Letters*, **13**(24), 5424–5430. https://doi.org/10.1021/ja190154z

Chen, W., Lickfield, G. C., & Yang, C. Q. (2004). Molecular modeling of cellulose in amorphous state. Part I: Model building and plastic deformation study. *Polymer*, **45** (3), 1063–1071. https://doi.org/10.1016/j.polymer.2003.11.020

Chen, Y., Jiang, X., Wu, H., & Zheng, L. (2020). Thermal behavior of complex model with the cellulose II and amorphous chain. *Journal of Theoretical and Computational Chemistry*, **19**(3), 1–11. https://doi.org/10.1142/S0219633620400040

Cheng, Q., & Wang, S. (2008). A method for testing the elastic modulus of single cellulose fibrils via atomic force microscopy. *Composites Part A: Applied Science and Manufacturing*, **39**(12), 1838–1843. https://doi.org/10.1016/j.compositesa.2008.09.007

Ciesielski, P. N., Wagner, R., Bharadwaj, V. S., Killgore, J., Mittal, A., Beckham, G. T., ... Crowley, M. F. (2019). Nanomechanics of cellulose deformation reveal molecular defects that facilitate natural deconstruction. *Proceedings of the National Academy of Sciences of the United States of America*, **116**(20), 9825–9830. https://doi.org/10.1073/pnas.1900161116

Cintrón, M. S., Johnson, G. P., & French, A. D. (2011). Young's modulus calculations for cellulose 1β by MM3 and quantum mechanics. *Cellulose*, **18**(3), 505–516. https://doi.org/10.1007/s10570-011-9507-1

Conley, K., Whitehead, M. A., & van de Ven, T. G. (2017). Probing the structural chirality of crystalline cellulose with induced circular dichroism. *Cellulose*, **24**(2), 479–486. https://doi.org/10.1007/s10570-016-1130-8

Cosgrove, D. J. (2005). Growth of the plant cell wall. *Nature Reviews Molecular Cell Biology*, **6**(11), 850–861. https://doi.org/10.1038/nrm1746

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

Cosgrove, D. J. (2022). Building an extensible cell wall. *Plant Physiology*, **189**(3), 1246–1277. https://doi.org/10.1093/plphys/kiac184

Cosgrove, D. J., & Jarvis, M. C. (2012). Comparative structure and biomechanics of plant primary and secondary cell walls. *Frontiers in Plant Science*, **3**(AUG), 1–6. https://doi.org/10.3389/fpls.2012.00204

Cresswell, R., Dupree, R., Brown, S. P., Pereira, C. S., Skaf, M. S., Sorieul, M., ... Hill, S. (2021). Importance of water in maintaining softwood secondary cell wall nanostructure. *Biomacromolecules*, **22**(11), 4669–4680. https://doi.org/10.1021/acs.biomac.1c00937

D'Acierno, F., Michal, C. A., & MacLachlan, M. J. (2023). Thermal stability of cellulose nanomaterials. *Chemical Reviews*, **123**(11), 7295–7325. https://doi.org/10.1021/acs.chemrev.2c00816

de Carvalho, D. M., Berglund, J., Marchand, C., Lindström, M. E., Vilaplana, F., & Sevastyanova, O. (2019). Improving the thermal stability of different types of xylan by acetylation. *Carbohydrate Polymers*, **220**(May), 132–140. https://doi.org/10.1016/j.carbpol.2019.05.063

De Coninck, J., & Blake, T. D. (2008). Wetting and molecular dynamics simulations of simple liquids. *Annual Review of Materials Research*, **38**, 1–22. https://doi.org/10.1146/annurev.matsci.38.060407.130339

de Ruijter, M. J., Blake, T. D., & De Coninck, J. (1999). Dynamic wetting studied by molecular modeling simulations of droplet spreading. *Langmuir*, **15**(22), 7836–7847. https://doi.org/10.1021/la990171l

Deng, X., Wu, Z., & Wang, G. (2020). Thermal conductivity and heat capacity of water, 1β/β cellulose nanofluids: A molecular dynamics study. *International Journal of Modern Physics B*, **34**(27), 1–13. https://doi.org/10.1142/S0217979220502483

Diaz, J. A., Wu, X., Martini, A., Youngblood, J. P., & Moon, R. J. (2013). Thermal expansion of self-organized and shear-oriented cellulose nanofiber films. *Biomacromolecules*, **14**(8), 2900–2908. https://doi.org/10.1021/bm400794e

Diddens, I., Murphy, B., Krisch, M., & Müller, M. (2008). Anisotropic elastic properties of cellulose measured using inelastic X-ray scattering. *Macromolecules*, **41**(24), 9755–9759. https://doi.org/10.1021/ma801796u

Jahedhi, C., Bergensträhle-Wohlert, M., Berglund, L. A., & Wohlert, J. (2016). Role of hydrogen bonding in cellulose deformation: the leverage effect analyzed by molecular modeling. *Cellulose*, **23**(4), 2315–2323. https://doi.org/10.1007/s10570-016-0968-0

Dri, F. L., Wu, X., Moon, R. J., Martini, A., & Zavattieri, P. D. (2015). Evaluation of reactive force fields for prediction of the thermo-mechanical properties of cellulose 1β. *Computational Materials Science*, **109**, 330–340. https://doi.org/10.1016/j.comatsci.2015.06.040

Ehmann, H. M., Werzer, O., Pachmajer, S., Mohan, T., Amenitsch, H., Resel, R., ... Spirk, S. (2015). Surface-sensitive approach to interpreting supramolecular rearrangements in cellulose by synchrotron grazing incidence small-angle X-ray scattering. *ACS Macro Letters*, **4**(7), 713–716. https://doi.org/10.1021/acsmacrolett.5b00306

Eichhorn, S. J., & Davies, G. R. (2006). Modelling the crystalline deformation of native and regenerated cellulose. *Cellulose*, **13**(3), 291–307. https://doi.org/10.1007/s10570-006-9046-3

Eichhorn, S. J., Dufresne, A., Aranguren, M., Marcovich, N. E., Capadona, R. J., Rowan, S. J., ... Peijs, T. (2010). Review: Current international research into cellulose nanofibres and nanocomposites. *Journal of Materials Science*, **45**(1), 1–33. https://doi.org/10.1007/s10853-009-3874-0

Eichhorn, S. J., Young, R. J., & Davies, G. R. (2005). Modeling crystal and molecular deformation in regenerated cellulose fibers. *Biomacromolecules*, **6**(1), 507–513. https://doi.org/10.1021/bm049409x

Eichinger, T., Rahkila, J., Willför, S., & Xu, C. (2019). TEMPO-oxidized O-acetyl galactoglucomannan oligomers: Isolation and comprehensive structural elucidation. *Wood Science and Technology*, **53**(1), 71–85. https://doi.org/10.1007/s00226-018-1068-0

Elazzouzi-hafraoui, S., Nishiyama, Y., Putaux, J.-L., Heux, L., Dubreuil, F., & Rochas, C. (2008). The shape and size distribution of crystalline nanoparticles prepared by acid hydrolysis of native cellulose. *Biomacromolecules*, **9**, 57–65.

Endler, A., & Persson, S. (2011). Cellulose synthases and synthesis in arabidopsis. *Molecular Plant*, **4**(2), 199–211. https://doi.org/10.1093/mp/sgq079

Etale, A., Onyianta, A. J., Turner, S. R., & Eichhorn, S. J. (2023). Cellulose: A review of water interactions, applications in composites, and water treatment. *Chemical Reviews*, **123**(5), 2016–2048. https://doi.org/10.1021/acs.chemrev.2c00477