Chirality and bound water in the hierarchical cellulose structure

The origins of the unique properties of natural fibres have remained largely unresolved because of the complex interrelations between structural hierarchy, chirality and bound water. In this paper, analysis of the melting endotherms for bleached hardwood pulps indicates that the amount of non-freezing bound water (0.21 g/g) is roughly half of the amount of freezing bound water (0.42 g/g). We link this result to the two smallest constitutive units, microfibrils and their bundles, using molecular dynamics simulations at both hierarchical levels. The molecular water layers found in the simulations correspond quite accurately to the measured amount of non-freezing and freezing bound water. Disorder that results from the microfibril twist and amphiphilicity prevents co-crystallisation, leaving routes for water molecules to diffuse inside the microfibril bundle. Moreover, the simulations predict correctly the magnitude of the right-handed twist at different hierarchical levels. Significant changes in hydroxymethyl group conformations are seen during twisting that compare well with existing experimental data. Our findings go beyond earlier modelling studies in predicting the twist and structure of the microfibril bundle.


Introduction
The hierarchical structure of cellulose is a fundamental basis not only for the majority of living plants, but also for several industries that take advantage of the striking features of natural fibres. The characteristic stiffness, hydrophilicity and water-induced swelling of fibres makes both strong natural and man-made materials possible at low density. In recent years, extensive experimental characterisation has gradually filled up the picture of cellulose synthesis and the subsequent build-up of microfibrils and their aggregates within macroscopic fibres (Jarvis 2018). At the same time, rapid development in the field of computational materials research has made atomistic simulations of cellulose microfibrils, and even larger structures, feasible. In fact, the smallest cellulose structures reachable with atomic force microscopy (AFM) (Ding et al. 2014;Usov et al. 2015;Zhang et al. 2016) are now being studied by all-atom molecular dynamics (MD) simulations (Langan et al. 2014;Oehme et al. 2015a;Kulasinski et al. 2017). In this paper, we extend such simulations to twisted macrofibrils, showing interesting interrelations between the underlying microfibril bundle structure and our measurements on bound water.
The basic unit of cellulose synthesis is a hexameric rosette structure from which elementary fibrils (in plant biology called microfibrils) have their origin . The rosette has been thought to lead to a six-multiple of cellulose chains in a microfibril (Cosgrove 2014). Most speculations in the literature are directed towards 18 (Nixon et al. 2016), 24 (Fernandes et al. 2011) or 36 chains (Cosgrove 2014) forming the fibrils in vascular plants. However, a partial fusion after synthesis could also give rise to other chain numbers (Jarvis 2018). According to a recent solid-state nuclear magnetic resonance (NMR) spectroscopy study on Arabidopsis and Brachypodium primary cell walls (Wang and Hong 2016), a microfibril seemed to contain two types of interior chains, referred to as surface-bound and core chains. By assuming one-chain difference between neighbouring cellulose layers, this leaves out the 18-chain model, as no core chains are possible in such a structure. On the other hand, the 36-chain model seems too large in many cases (Fernandes et al. 2011;Cosgrove 2014). By combining a range of spectroscopic methods with small-angle neutron and wide-angle X-ray scattering measurements, Fernandes et al. (2011) suggested an amphiphilic 24-chain microfibril for spruce wood. Such fibrils contain both hydrophilic surfaces, which expose the equatorial hydroxyl groups of the pyranose rings, and less hydrophilic surfaces, which expose the axial hydrogens.
Cellulose chains and microfibrils are the first two levels of structural hierarchy in natural fibres (see Fig. 1). Microfibrils have been frequently observed to bundle into macrofibrils (Ding et al. 2014;Zhang et al. 2016;Cho et al. 2017), which in turn form the lamellae of the layered cell wall structure. In principle, a bundle may contain any number of microfibrils, but smaller bundles are more frequent than large ones (Fahlén and Salmén 2003;Zhang et al. 2016). An intriguing feature of this structural hierarchy is the chiral order that seems to be repeated over several length scales (Wang et al. 2013). Simple continuum and geometric models have been used to discuss how chirality is transferred from macrofibrils to fibres (Burgert et al. 2005), or from fibres to their networks (ten Bosch 1996). The latter is seen in the twisting of macroscopic paper samples when their moisture content is varied (Yu et al. 2004). Mechanical analysis (ten Bosch 1996) suggests that the mean twist rate decays inversely with the size of the structural dimensions. This feature has also been seen in molecular dynamics simulations of single microfibrils (Hadden et al. 2013). The twisting of a natural fibre is generally associated with a nonzero microfibril angle (Burgert et al. 2005), and is generally not expected to be correlated with similar behaviour in the microfibrils (Usov et al. 2015) or in their bundles (Ding et al. 2014). The complexity of the structure at the fibril level (see Fig. 1) makes experimental characterisation extremely difficult. Contrary to some earlier claims (Zhang et al. 2016), our simulations suggest the existence of significant twisting in the macrofibrils in the scale of a few hundred nanometers. Similar helical structures have been recently observed in the S1 layer of spruce using cryo-transmission electron tomography (Reza et al. 2017).
Several molecular-level mechanisms have been proposed as the driving force behind the twisting of single microfibrils. Since the phenomenon was first observed in classical MD simulations (Yui et al. 2006;Matthews et al. 2006;Yui and Hayashi 2007), it has been attributed, solely or in part, to (1) changes in intra and interlayer hydrogen bonding patterns (Paavilainen et al. 2011) (2) attractive, interlayer, van der Waals forces (Hadden et al. 2013) (3) torsional forces due to the geometry of the trans-glycosidic hydrogen bonds, especially the HO2-O6' bond (Bu et al. 2015;Conley et al. 2016), and (4) a joint effect of several mechanisms (Kannam et al. 2017). In a related discussion (Zhao et al. 2013), it has been speculated that amorphous regions repeating along the length of the fibril would be needed to release stresses involved with the twist, and to keep the cellulose I b structure intact. However, there is experimental evidence for deviations from the cellulose I b structure in natural fibres (Fernandes et al. 2011;Wang and Hong 2016). The level of disorder, as measured by hydroxymethyl group conformations, seems to depend on the location within the microfibril cross-section, with the transgauche (tg) conformer 1 predominant in the microfibril core (Fernandes et al. 2011;Wang and Hong 2016). Most of the theoretical studies on microfibril twist have focused on the underlying molecular mechanisms. There are a few simulation studies on microfibril bundles (Langan et al. 2014;Kulasinski et al. 2015Kulasinski et al. , 2017Oehme et al. 2015a;Silveira et al. 2016), but the chiral order in these larger structural units has not been shown in any of the previous works. Fig. 1 Hierarchical levels of the cellulose fibre structure. From left to right: molecular visualisation of a layer of cellulose chains (in cellulose I b , hydrogen bond network shown in green), a microfibril, and a bundle of four microfibrils (cross-sections shown in the insets; outermost chains coloured in red, others in grey); followed by scanning electron micrographs of the S2 layer of birch (Lovikka et al. 2018) and a whole fibre. The orange and blue scale bars in the micrographs indicate 1 lm and 10 lm, respectively In this work, we bridge the gap between experiments and simulations by showing, using an atomistic model, how the chirality of cellulose microfibrils is inherited by their bundles. The twist rates predicted by our simulations are similar to those found in recent AFM and cryo-transmission electron tomography studies on microfibrils and microfibril bundles (Ding et al. 2014;Usov et al. 2015;Reza et al. 2017). According to our simulations, the microfibril twist and amphiphilicity affect the distribution of water molecules near microfibril bundles. Thus, we compare our simulations against thermoporosimetry measurements on the amount of non-freezing and freezing bound water in wood fibres at low hemicellulose content. The simulations were carried out assuming that microfibrils consist of 24 chains, as suggested for spruce wood (Fernandes et al. 2011). However, we expect that most of the qualitative conclusions are insensitive to the exact number of chains.

Materials
Six commercial pulps were used in the experiments. A never-dried bleached birch kraft pulp (BKHW-ND) was collected from a Finnish pulp mill. Some of the pulp was dried at 40°C and then rewet (BKHW-40). Some of the dried pulp was further heat treated at 100°C over night (BKHW-100). The drying and heating causes the fibrils to aggregate, swelling to reduce and surface area to decrease. This process is known as hornification. A sample of commercial never-dried birch dissolving pulp was also collected from a Finnish pulp mill (BDHW-ND). Some of the pulp was also collected from the pulp dryer in the same production line (BDHW-dried). Some of this pulp was further heat treated at 100°C to increase hornification (BDHW-100).

Porosity
The porosity of the pulps in the wet state was evaluated by a two-point solute exclusion method. The fibre saturation point (FSP) is the water inside the cell wall, evaluated with a 2 9 10 6 Da dextran molecule with a Stokes diameter of 54 nm. A 5 kDa dextran with a 3.2 nm diameter was also used to evaluate inaccessible water. In this case, the probe dextran can penetrate into the cell wall, but not into the micropores within a fibril aggregate. Thus, we refer to the water measured with this approach as ''fibril swelling''.
Another way to determine pore structure is by using solvent exchange, followed by critical point drying from a supercritical CO 2 solvent. This approach somewhat preserves the pore structure of the fibres and allows evaluation with classical porosity methods. In this case, N 2 sorption with application of the Brunauer-Emmitt-Teller model (Brunauer et al. 1938) was used to determine surface area and pore volume from the adsorption isotherms.

Bound water
Interfacial water has the outstanding property that it does not freeze, even at very low temperatures. Nonfreezing water (NFW) is one measure of bound water and gives an indirect way to probe the amount of hydrated surface area. Beyond the first monolayer, several layers of water will have a reduced melting temperature compared to bulk water. This gives rise to a second fraction called freezing bound water (FBW). The NFW and FBW can be determined from water saturated pulp samples by careful analysis of the melting endotherms under controlled conditions. The NFW and FBW were measured in a carefully calibrated Mettler DSC3?. Samples were 2-3 mg with a moisture content of 0.8-1.2 g water/g solids. The FBW is defined as the water in the sample, which melts at a temperature less than -0.2°C. Further details of the method can be found in (Maloney and Paulapuro 1999;Maloney 2015).

Fibre twisting
In order to reflect the behaviour on the macrofibril level against the macroscopic scale, the twisting behaviour of pulp fibres was also considered. BKHW-ND was dried in a turbulent air stream to prevent the formation of fibre-fibre bonds and to allow the free shrinkage of individual fibres. Under these drying conditions, the fibres experience a torsional strain and twist about their longitudinal axis. To quantify the twisting, dried fibres were investigated by scanning electron microscopy (SEM, Zeiss Sigma 123 VP). Prior to the imaging, the fibres were placed on carbon tape at low concentrations. The procedure was repeated for four replicates, which were simultaneously sputtered with a 4 nm layer of gold. For each replicate, a minimum of three images was collected. The area covered by the images was randomly selected. Magnification was adjusted to ca. 100X. The images were analysed using the ImageJ Fiji software. For each image, the number of 360°twists of ca. 5-15 random fibres was counted and the fibre length measured. The number of twists per mm was then calculated.
In this work, MD simulations of the twisting behaviour and water interactions of cellulose microfibrils are extended to their bundles-the next level of structural hierarchy in a cellulose fibre (see Fig. 1). The simulations were carried out in both vacuum and water environments to determine the qualitative response of the microfibrils to surrounding water. An individual microfibril was first simulated in vacuum until it adopted an equilibrium twist around its longitudinal axis (see Fig. 2). Four replicas of the fibril were then aggregated and again simulated in vacuum until equilibrium. After this, both the single fibril and the bundle were simulated in a water environment. The twisting behaviour, along with conformational changes of the hydroxymethyl groups, was quantified from the level of individual chains up to the twisting of fibrils around the bundle axis. Water penetration into the fibril bundle and related structural changes were also followed.
The MD simulations were performed using the GROMACS software package (Pronk et al. 2013). A variant of the OPLS-AA force field for carbohydrates (Jorgensen et al. 1996;Damm et al. 1997;Price et al. 2001;Kony et al. 2002), adapted for cellulose by Paavilainen et al. (Paavilainen et al. 2011), was used to describe chemical and physical bonding (see Supporting Information for details). The OPLS-AA force field is optimised for condensed-phase simulations, and thus the results obtained in water are considered more reliable than those in vacuum. However, the comparison with vacuum simulations gives interesting insights to the contribution of water molecules in conformational and chirality changes. A broad description of the simulations is given in what follows. Further details of the force field and its parameters, the simulation set-up and the post-processing steps are given in the Supporting Information.

Atomistic model and simulations of single fibrils
The atomistic model of a cellulose microfibril was constructed using crystallographic unit cell data given for cellulose I b by Nishiyama et al. (2002). The fibril consists of 24 cellulose molecules, each 40 glucose residues long. This crystallite size, containing both surface-bound and core chains, has been suggested by Fernandes et al. (2011) for wood, and is supported by the solid-state NMR studies of Wang and Hong (2016). The stacking of cellulose molecules is such that the fibril's cross-section has an elongated hexagonal shape (see Fig. 2, upper left corner). This results in two kinds of exterior surfaces: ones that expose the hydroxymethyl groups, corresponding to the ð110Þ and ð1 10Þ crystallographic planes, and ones that expose the non-polar hydrogens, corresponding to the ð100Þ and ð200Þ planes. The height of the resulting fibril (2.8 nm; see Fig. 2, upper left corner) agrees well with the wide-angle X-ray scattering results by Fernandes et al. (2011) for the microfibril thickness (2.9 nm) normal to the ð200Þ plane. Moreover, in the fibril model, the average thickness normal to the ð110Þ and ð1 10Þ planes is 2.5 nm, which is close to the experimental mean of 2.6 nm in the same directions.
The cellulose is initially in the I b crystalline form, with hydrogen bonding corresponding to the network A structure (Nishiyama et al. 2008). Consequently, all the hydroxymethyl groups are in the trans-gauche (tg) conformation. The length of the fibril, 21 nm, was chosen so that the twist in both single fibrils and their bundles could be simulated with sufficient accuracy while keeping the system size feasible for simulations over several hundred nanoseconds. The length is roughly twice the distance over which crystallographic coherence is maintained in experiments (Fernandes et al. 2011). In the current simulations, no amorphous regions along the fibrils were assumed. Instead, disorder was introduced by adding loose cellulose chains to the fibril bundle model. The main motivation was to crudely describe the effect of hemicellulose, but similar disorder could also be caused by chain ends or mechanical damage, as suggested by Jarvis (2018).
The MD simulation set-up for the individual fibrils was the following. A single fibril was embedded in a periodic simulation domain of (10.2 9 12.6 9 29.2) nm 3 . Potential energy minimisation was then performed to bring the system to a local minimum. Next, the system was equilibrated in a 500 ns MD simulation in the canonical (NVT) ensemble, at 293 K temperature. After the vacuum equilibration, the simulation domain was filled with water molecules, yielding a total density of 0.99 g/cm 3 . After the addition of water, a potential energy minimisation was again performed. The system was then equilibrated in a 500 ns MD simulation in the isothermal-isobaric (NPT) ensemble, at 293 K temperature and 1 atm pressure. Two parallel simulations were performed at each stage to assess how sensitive the system was to the choice of initial velocity vector. The microfibril twist rate (°/nm) was defined as the average twist rate of its constituent cellulose chains around the fibril axis. Fig. 2 The simulation procedure began from the upper-left 24-chain microfibril model (geometric mean diameter 3.5 nm) with surface (yellow), sub-surface (cyan) and core (indigo) cellulose chains. After equilibration in vacuum, the twisted fibril was studied either in water as such (surface water molecules shown in orange), or four copies of the fibril were used to form a bundle. The bundle model contained also added disordered chains (shown in dark red). The bottom-middle image shows the bundle after aggregation and twisting in vacuum. Finally, the response of the bundle to a water environment was simulated (bottom right)

Simulations of fibril bundles
The model of a fibril bundle was created as follows. Four replicas of the vacuum-equilibrated fibril were arranged in a cruciform setting parallel to each other (see Fig. 2, lower left corner). The two central fibrils had ð100Þ and ð200Þ planes facing each other at the reducing end, while the two peripheral fibrils had ð110Þ and ð1 10Þ planes facing each other. Due to the vacuum-equilibration, the cellulose chains had adopted a helical twist around the longitudinal axis of the fibril. Thus, the crystallographic planes facing each other would change along the length of the bundle. Six cellulose molecules, each 40 glucose residues long, were added between the fibrils to represent hemicellulose, or other disordered forms of cellulose (e.g., loose chain ends sticking out from a microfibril), which are expected to exist together with crystallites (Jarvis 2018). These extra chains amounted to 6 wt% of all cellulose, which corresponds to a typical amount of hemicellulose in dissolving pulp (Khanjani et al. 2017). Moreover, dual-axis electron tomography (Xu et al. 2007) has shown an intimate packing of hemicellulose with microfibril bundles.
The four fibrils and the free chains were embedded in a periodic simulation domain of (17.4 9 19.8 9 30.5) nm 3 . Potential energy minimisation was then performed to bring the system to a local minimum. The fibrils were aggregated and equilibrated in a 100 ns MD simulation in the canonical ensemble (NVT), at 293 K temperature. After the vacuum equilibration, the simulation domain was filled with water molecules, yielding a total density of 1.00 g/ cm 3 . After the addition of water, a potential energy minimisation was again performed. The system was then equilibrated in a 300 ns MD simulation in the isothermal-isobaric (NPT) ensemble, at 293 K temperature and 1 atm pressure. Unlike in the single fibril case, only one simulation was performed for the bundle. The bundle twist rate (°/nm) was defined as the average twist rate of its constituent microfibrils around the bundle axis.

Results and discussion
The simulations revealed strong couplings between the chirality, chain conformations, hierarchical structure and water interactions. None of these fundamental properties can thus be studied alone. Rather than focusing on a specific property, our aim is to form an overall picture of the interrelations and demonstrate the power of the simplified atomistic model when reflected against current experimental evidence.
Our simulations were able to predict several measurable parameters with reasonable accuracy, taking into account the limited statistics and the small number of systems studied. Table 1 shows this for a number of key parameters at different levels of the structural hierarchy. The predicted values are compared with both measurements carried out in this paper and earlier experimental works. The details of the simulation results and their comparison with the experiments are discussed in further below.

Twisting dynamics and the hierarchical chiral order
The simulations used a non-twisted, perfectly crystalline, microfibril as the starting point. However, the equilibrium state in both vacuum and water is a twisted structure where the initial tg conformations of the hydroxymethyl groups are replaced by a mixture of tg, gt and gg conformations. The conformational disorder increases from the microfibril core towards the surface, with the tg conformer dominant in core chains, and the gt conformer in surface chains. Similar disorder has been observed in solid-state NMR experiments (Fernandes et al. 2011;Jarvis 2018) and earlier MD simulation studies on cellulose microfibrils (Matthews et al. 2012a;Oehme et al. 2015bOehme et al. , 2018. The dynamics of twisting and conformational changes are strongly coupled so that the overall conversion from tg to other conformations is linearly related to the fibril twist rate after the rapid relaxation that takes place during the first nanosecond (see Fig. 3a). Moreover, the rate of the changes is similar in vacuum and in water. The conversion begins from the surface chains, and is followed by more gradual changes in the sub-surface and core chains (see Fig. 3b). The final proportions of tg/gt/gg conformers are 57%/37%/6%. However, one should notice that this result may be characteristic of the chosen parametrisation of the force field (see Supporting Information) and its sensitivity to possible alternative parametrisations should be studied further.
The twisting and conformational changes are associated with a gradual decrease in the force field term for electrostatic interactions, and a corresponding increase in the bonded (i.e. covalent) and van der Waals terms. This observation is in agreement with previous simulation studies that have associated microfibril twist with changes in the intrachain, intra-layer, and inter-layer hydrogen bonding patterns (Paavilainen et al. 2011;Bu et al. 2015;Conley et al. 2016;Kannam et al. 2017).
Interactions with water lead to increasing conformational disorder and an increased twist rate, even though there is no penetration of water molecules into the fibrils. In the water environment, the final proportions of tg/gt/gg conformers are 26%/69%/5%. The equilibrium twist rate saturates at 4.2°/nm, which agrees with the result of Kannam et al. (2017). The difference to the rate in vacuum, 3.6°/nm, is significant and should be observable in experiments. Bundle 20 Fibre 28-54 (Fernandes et al. 2011;Väisänen et al. 2018) Twist rate (°/nm)  Second and third water monolayers, see (Niinivaara et al. 2015) c Simulated bundle of four microfibrils with geometric mean diameter of 10.1 nm. The experimental bundles have varied diameters and numbers of fibrils d Extrapolated from the bundle using 1=D scaling, see Table 2 123 The simulated aggregation and vacuum equilibration results in a twist rate of 0.8°/nm for the bundle, and slightly reduced twist rates for the constituent fibrils (2.6-3.1°/nm). The resulting bundle (see Fig. 4, right) has a ribbon-like shape, which was found to be typical in earlier microscopic imaging (Ding et al. 2014;Usov et al. 2015). All changes took place rapidly during the first 3 ns, and the structures remained practically unchanged for the rest of the simulation. However, when putting the vacuum-equilibrated Fig. 3 Twisting and conformational changes in the single fibril case. a Twisting dynamics of a single fibril in vacuum (red curve) and after immersing the twisted fibril into water (blue curve). Similar characteristic times in these two cases suggest that they have the same driving mechanisms. The black curves show linear fittings of the conformational changes (1 À N tg =N all ) to the simulated twist rates. b The amount of water-induced conformational changes for the surface (gold), sub-surface (cyan) and core chains (indigo) during twisting (cf. blue curve in a). Similar changes and dynamics were also observed in vacuum. The figures on the right show the final hydroxymethyl group conformations in the upper-half layers (L1-3) of the fibril in water. The glucose units are coloured according to the conformation (see colour coding at the bottom) bundle in a water environment, a similar gradual change in both hydroxymethyl group conformations and twist rate were observed, as in the case of individual fibrils (see Fig. 2). The final proportions of tg/gt/gg conformers are 22%/71%/6%. The bundle twist rate settled at a level of 0.9°/nm, whereas the twist rate of the four constituent fibrils reached an average level of 4°/nm, similar to the case of a single fibril in water.
As shown in Table 1, the twist rate of both single fibrils and the bundle agree very well with earlier experimental characterisations. For bundles, the predicted level 0.9 8/nm corresponds well to the experimental value of 1.0 8/nm obtained by Reza et al. (2017) using cryo-transmission electron tomography. However, here the comparison is only between the orders of magnitude as the exact number of fibrils in the experimental bundle is not known. The other experimental values are based on re-analysis of previously published AFM and transmission electron microscopy (TEM) micrographs (Usov et al. 2015), where the span of repeating geometrical features had been indicated but not interpreted as twisting. The somewhat larger experimental values 1.1-1.9 8/nm in these cases are explained by the corresponding smaller bundle diameters suggesting a smaller number of fibrils (i.e. two or three) in a bundle. The alternation in fibril orientation and surface interactions in the constructed fibril bundle is similar to the length over which crystallographic coherence is maintained in experiments (Fernandes et al. 2011). Moreover, the proportions of hydroxymethyl group conformations can be compared with those estimated based on 13 C spin-lattice relaxation times (Fernandes et al. 2011). It appears that the experimental values fall roughly in between the simulation results for water and vacuum environments. This is reasonable, since the simulated conditions reflect the two extremes of moisture content. In comparison to the single fibril simulations, the overall proportion of tg conformations is slightly reduced. Oehme et al. (2015a) observed similar behaviour in their simulations of microfibril aggregation.
The lower twist rate in a bundle compared to a single fibril is explained by its larger effective diameter D. The 1=D scaling predicted by continuum models on the chiral elasticity (ten Bosch 1996) is crudely obeyed by the first two levels of structural hierarchy (fibril, bundle). The 1=D scaling was found earlier through atomistic simulations for a single microfibril with varied diameter (Zhao et al. 2013;Hadden et al. 2013). However, it is curious that the scaling also holds for the twist rate of the bundle (see Table 2), despite the much faster twist rates of the Fibre (experiment) 2500 3.5Á10 -3 9 The simulated scaling for the fibril and bundle extends to the measured twist in hardwood kraft fibres. Da gives effectively the fibril (or chain) angle Fig. 4 Twist rate of the microfibrils around the bundle axis in water environment (dashed lines), together with the twist rates of the four constituent fibrils around their own axes (solid lines). Colour legend for the constituent fibrils is shown on the right. Black lines indicate the average 123 individual fibrils inside the bundle. Moreover, the 1=D scaling is found experimentally to extend to two orders of magnitude larger hardwood kraft fibres, as shown in Table 2. The twist rate a of a wood fibre is dominated by its thick S2 layer. The preferred orientation l of fibril bundles in that layer (Barnett and Bonham 2004) is given by l % tanl ¼ Da (in radians) assuming that l follows the fibre twist (Burgert et al. 2005). Interestingly, the estimate Da obtained for the experimental fibril angle for hardwood kraft fibres is on a similar level to the value obtained for the simulated bundle even though one would not expect direct correlation between these two scales.
Bound water in the bundle The mismatch of the twist rates for the bundle and its fibrils led to a longitudinal variation in the bundle structure, which prevented co-crystallisation and opened up routes for water molecules to diffuse inside the bundle. The changing types of inter-fibril surface contacts are illustrated in Fig. 5 together with the distribution of water molecules after 300 ns diffusion time. The twisting behaviour of the peripheral fibrils (F1 and F4) seems to be coupled, and the same holds for the central ones (F2 and F3, see Fig. 4). Moreover, different crystallographic planes are exposed to water at varied locations along the length of the bundle. The most compact form (Fig. 5f)  (Fig. 5a). Significant diffusion of water molecules into the inner parts of the bundle is visible especially in Fig. 5c, d with steep crevices. Such inner bound water has been suggested by earlier experimental studies (Fernandes et al. 2011).
Even though some water molecules penetrate into the bundle, most of the bound water is located at outer bundle/fibril surfaces. In the single fibril case, the average number of directly hydrogen-bonded water molecules corresponds to 0.13 g/g of cellulose, while the average number of water molecules within the potential range of hydrogen bonding (3.5 Å ) sums up to 0.26 g/g. Aggregation eliminates some of the fibrils' exterior surface, reducing the number of accessible hydrogen bonding sites. Consequently, the number of directly hydrogen-bonded water molecules corresponds to 0.10 g/g cellulose, whereas 0.19 g/g is within the hydrogen bonding range. The reduction of the bound water mass fraction is proportional to the reduction of specific surface area by 25% from 870 m 2 /g for the single fibril to 520-630 m 2 /g (probe diameter 0.6-3.6 nm) for the bundle. The specific surface area of the bundle is quite close to the value, 470 m 2 /g, obtained by extrapolating the values obtained with nitrogen sorption experiments to the fibre saturation point of never-dried hardwood (see Fig. 6).
Vapour adsorption experiments have shown that the hydration of cellulose nanocrystals (CNCs) involves the adsorption of a circa 1 nm thick layer of water (three monolayers of water molecules) onto the CNC surfaces (Niinivaara et al. 2015). In our simulation of the microfibril bundle, water molecules within the hydrogen bonding range correspond to roughly one monolayer (see Fig. 5). The mass fraction of such water, 0.19 g/g cellulose, is remarkably close to the non-freezing bound water content of 0.22 g/g observed in the thermoporosimetry experiment on never-dried bleached hardwood pulp (see Table 3). Moreover, the combined mass fraction of the second and third water layers amounts to 0.49 g/g, which agrees well with the freezing bound water content of 0.42 g/g observed in the same experiment. The distribution of bound water is affected by amphiphilicity of the fibril surfaces. This was seen in an additional vapour adsorption simulation of the single fibril case with reduced water content. The initially evenly distributed water molecules (0.2 g/g cellulose) gathered into hydrophilic areas, leaving parts of the hydrophobic surfaces unoccupied (see Fig. 2, upper right corner).

Conclusions
Our experiments on bleached hardwood pulps indicate that the amount of non-freezing bound water is roughly half of the amount of freezing bound water. In order to see the reasons for this behaviour, we built a molecular model which goes beyond the individual microfibril. The large hierarchical system achieved in our simulations opens an avenue for probing various physical, chemical and structural properties, including chirality, crystalline order, fibril aggregation, amphiphilicity, and water interactions at varied moisture content.
In this work, we discussed the predictions for both individual microfibrils and their bundles, considering state-of-the-art experimental evidence for cellulose. Our relatively simple model correctly describes a surprising number of experimental features of cellulose, including fibril/bundle twist, proportions of hydroxymethyl group conformations, specific surface areas, and the mass-fractions of bound water. Moreover, we show that these properties are linked to one another. In particular, the conformational changes are Fig. 6 Specific surface area versus pore volume in the nitrogen sorption experiments. The linear least squares fit is used to extrapolate the specific surface area data of Table 3 to the fibre saturation point of never-dried bleached hardwood pulps (0.99 g water/g cellulose) 123 linearly related to the twist rate, which itself affects the packing of individual fibrils into a bundle. This packing further determines the surface area and the distribution of bound water in the bundle. In contrast to some earlier modelling observations (Langan et al. 2014), the co-crystallisation of microfibrils is prevented by a mismatch of the twist rates at different hierarchical levels. This creates repeating openings to the internal bundle structure that in principle could become sites for defect initiation in chemical and mechanical processing. However, the current simulations are not sufficient to make this conclusion. Moreover, we did not study kinks or branching that appear in long bundle structures. The small amount of disordered chains included in our bundle model had a minor role in the overall structure. Repeating patterns in earlier AFM micrographs (Ding et al. 2014;Usov et al. 2015) could be explained quantitatively by characteristic fibril and bundle twists.
The conformational disorder that increases from the microfibril core towards the surface is also in accordance with previous studies. A similar structural feature has been observed experimentally, as reported by Fernandes et al. (Fernandes et al. 2011) and Wang and Hong (Wang and Hong 2016), as well as in MD simulations, as reported by Matthews et al. (Matthews et al. 2012a) and Oehme et al. (Oehme et al. 2015b(Oehme et al. , 2018. Comparison between the simulations in vacuum and water environments suggests that water interactions introduce further disorder in the surface chains, from where it spreads towards the microfibril core. The current predictions are most relevant for chemical wood fibres with low hemicellulose content, in which other matrix polymers play a minor role. More findings can be expected to emerge in the future when more statistics on the effects of structural variation (e.g. the number of fibrils in a bundle) is accumulated. When combined with tailored experiments, the simulations may shed light on the exact chain number or fibril morphology in different cases. We hope our results will stimulate further experimental characterisations, which, together with hierarchical simulation models, will sharpen the picture of the fundamentals of cellulose.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits The never-dried pulps are considered to have the closest correspondence to the atomistic models developed for this work. The other cases are shown to indicate typical variation in properties coming from drying conditions. In particular, the extended data allows extrapolation of the specific surface area of the nitrogen sorption experiments to the relative pore volume obtained with the solute exclusion technique, see Fig. 6 unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.