University of Birmingham Effects of Moisture on the Mechanical Properties of Microcrystalline Cellulose and the Mobility of the Water Molecules as Studied by the Hybrid Molecular Mechanics–Molecular Dynamics Simulation Method

: A hybrid molecular mechanics – molecular dynamics simulation method has been performed to study the effects of moisture content on the mechanical properties of microcrystalline cellulose (MCC) and the mobility of the water molecules. The speci ﬁ c volume and diffusion coef ﬁ cient of the water increase with increasing moisture content in the range studied of 1.8 – 25.5 w/w %, while the Young ’ s modulus decreases. The simulation results are in close agreement with the published experimental data. Both the bound scission and free-volume mechanisms contribute to the plasticization of MCC by water. The Voronoi volume increases with increasing moisture content. It is related to the free volume and the increase enhances the mobility of the water molecules and thus increases the coef ﬁ cient of diffusion of the water. Moreover, with increasing moisture content, the hydrogen bonding per water molecule between MCC – water molecules decreases, thus increasing the water mobility and number of free water molecules. © 2019 The Authors. Journal of Polymer Science Part B: Polymer Physics published by Wiley Periodicals, Inc. J. Polym. Sci., Part B: Polym. Phys. 2019 , 57 , 454 – 464

INTRODUCTION Microcrystalline cellulose (MCC) is widely used in the pharmaceutical industries as an excipient in combination with other active pharmaceutical ingredients. MCC has special properties such as inactivity, the absence of toxicity, and considerable hygroscopicity. 1 It is essentially alpha cellulose derived mainly from wood pulp, which has been processed to produce particles by acid hydrolysis. Cellulose itself is a linear homopolymer composed of anhydroglucose units connected by 1,4-βglucosidic bonds. There are four major polymorphs of cellulose, namely, Cellulose I, II, III, and IV, where the most abundant native crystalline form is Cellulose I. 2 There are two types of Cellulose I, Cellulose Iα, which is metastable, and Cellulose Iβ that is more stable. 3 To improve the performance of the oral dosage form of pharmaceutical drugs and develop the processing, an understanding of the physical properties of the active and excipient components is crucial. 4 Twin-screw granulation has been considered as a potential processing route for the size enlargement of the feed powders for pharmaceutical tablets. It is a continuous wet granulation method that employs an open-channel twin-screw extruder to form the granules. When MCC is processed using any wet granulation technique, which involves liquid-solid interactions, the effects of water addition on the properties of MCC will influence the final quality of the oral dosage product. For example, tablets made from nongranulated MCC have a greater strength than those made from wet granulated MCC. 5,6 A wet granulation process also reduces the compactibility of MCC to various degrees. 7 Moreover, the Young's modulus and tensile strength of MCC have been found to decrease as the moisture content increases. 8 MCC granules formed by a roller compactor under constant operational settings and mass feed rate show a decrease in tensile strength and Young's modulus with increasing moisture content. 9 The tensile strength, elongation, and elasticity of cotton fibers vary with the water content in cellulose up to about 20% water content. 10 At 50 and 80 C, the dynamic modulus of regenerated cellulose has been shown to decrease rapidly at about 90-95% relative humidity, which may be due to the onset of micro-Brownian motion in the plasticized chain molecules. 11 In addition to experimental studies, molecular dynamics (MD) simulations have been conducted to examine the interaction of water and cellulose. This technique supports experimental studies with more insights about the mechanisms governing the interaction at the atomistic/molecular levels. MD simulations also provide an alternative approach for predicting water diffusion in MCC, which is difficult to be measured experimentally, for example, using the pulsed-field-gradient stimulated-echo nuclear magnetic resonance (NMR) technique. 12, 13 Heiner et al. 14 performed MD simulations of the interface between water and four surfaces of crystalline cellulose. They found that only the first surface layer of the crystalline cellulose was structurally affected by the water while the 1(−1)0 and 100 surfaces were more hydrophilic than the 110 and 010 surfaces. Matthews et al. 15 reported MD simulations of MCC in water. They showed that the surfaces of MCC highly structured the water solution in contact with them primarily because of direct hydrogen bonding with the first layer of the cellulose. They suggested that the 110, 1(−1)0, and 010 surfaces might be more susceptible to hydrolysis than the 100 surface. Yui et al. 16 presented MD simulations of the swelling behavior of the cellulose crystal. Its crystal twisted rapidly to form a slightly right-handed shape and then remained in that form as a steady swollen state for the remainder of the simulation. Mazeau and Rivet 17 studied the wetting of cellulose surfaces by water using MD simulations. They found that the wetting of the 100 surface was less than that of the 110 surface.
The objective of the current study is to investigate the influence of absorbed water on the physical properties of MCC. For validation purposes, the specific volume and elastic modulus are compared with the available experimental values corresponding to relatively low moisture contents (>10 w/w%). The model then is employed to calculate the specific volume and elastic modulus at higher moisture contents. It also is employed to calculate the water diffusion coefficient in MCC, which to the authors' knowledge, published data are not available. The accuracy of these data is validated by the calculating the activation energy for a water molecule in MCC, which is compared to the experimental value. Moreover, the aim is to evaluate the contributions of the free volume and bond scission mechanisms of plasticization by water based on the mobility of the water molecules in MCC. In this context, the roles of hydrogen bonds and mobility of the water molecules are discussed. This work also serves as a further validation for our recently developed hybrid molecular mechanics (MM)-MD simulation method 18 for the calculation of the Young's modulus of organic polymers. The method overcomes the limitation of conventional MD, that is, extremely high strain rates are involved to calculate the Young's modulus.

MODEL AND SIMULATION
A 5 × 5 × 5 unit cell of an MCC model containing 12,810 atoms was created using the crystalline cellulose-Atomistic Toolkit. 19 The 5 × 5 × 5 unit cells of an MCC consist of 11 cellulose layers each has 10 chains of 10 cellulose units (C 6 H 10 0 5 ). This is slightly larger than the model employed Heiner et al., 14 that is, six cellulose layers, each consisting of six chains of three cellulose units, to study the interface between water and native crystalline cellulose by MD simulations. Moreover, periodic boundary conditions were applied in the MD simulations to achieve a bulk representative of MCC. The atomic coordinates of cellulose Iβ were those from ref. 20 , which were determined by synchrotron X-ray and neutron diffraction studies. Figure 1 shows a unit cell of this cellulose model as visualized using Jmol. 21 The ReaxFF force field was based on the parameter set developed for glycine conformers and glycine-water complexes. 22 It has been reported that the ReaxFF force field is able to predict some properties of cellulose with reasonable accuracies: lattice constants and angles (deviations from experimental values between <1 and~6%), elastic constants (within the range of experimental values) and thermal expansion (within the range of experimental values for one direction but not for two other directions). 23 The Young's moduli of crystalline cellulose predicted using the ReaxFF force field are about 190 and 30 GPa for the directions parallel and perpendicular to the cellulose axial direction, respectively. 23 These values are within the range of the corresponding experimental values, which are about 120-220 and 8-50 GPa. [24][25][26][27] Therefore, in terms of the lattice parameters and elastic constants, this force field can be considered as a reasonable model to represent the properties of MCC for the current study. Although the force field is not able to predict all thermal expansion constants accurately, they are not the current focus, which is to understand more quantitatively the origin of the attenuation in the Young's modulus resulting from plasticization. A particular advantage of the ReaxFF force field is that it incorporates an explicit hydrogen bonds energy function, thus it can provide explicit intrachain and interchain hydrogen bonding information that is required in this study. Moreover, the water parameters in this force field have been validated for bulk water by comparison against experimental data on the cohesive energy, diffusion constant, and structure of water. 22,28 Consequently, it is expected to be able to predict the interaction between MCC and water molecules investigated here. All the calculations and simulations described in this section were performed using LAMMPS 29 with periodic boundary conditions. The Polak-Ribiere version of the conjugate gradient algorithm 30 was used to minimize the energy of the initial structure with a periodic boundary condition, using a force stopping criteria of 10 −9 . The structure was equilibrated for 0.45 ns at 300 K and 1 atm pressure with a timestep of 0.5 fs. The Young's modulus of the MCC (dry and with water) was calculated using the hybrid MM-MD method developed in ref. 18 . This method overcomes the limitation of conventional MD, that is, extremely high strain rates are involved to calculate the Young's modulus. The combination of MM and MD drives the system to the stable state faster than conventional MD, thus the Young's modulus can be calculated at a strain rate that corresponds to typical experimental quasi-static loading rates. For example, the method was able to reproduce accurately the effect of temperature on the Young's modulus of poly (methyl methacrylate) up to values greater than the glass transition temperature (T g ), 18 which is not possible with traditional MD. 31 The Young's modulus was calculated perpendicularly to the polymer chains and compared to values obtained from three-point bending test on compacted MCC. 8,32 Due to the compression in preparing the compact samples, microfibrils are expected to be orientated in a horizontal direction. 8 The MCC model described above was used to study the interaction of water with cellulose. Water molecules were randomly created using Packmol 33 and inserted into the cellulose bulk model. The inserted water molecules will modify the crystallinity of the cellulose and result in a more amorphous structure. Moreover, during the MD simulations, the thermal energy also perturbs the structure resulting in local disorder. The regions with greater local disorder show an ability to absorb more water as in the case of amorphous solids. 34 Therefore, the system imitates the condition where water molecules have already penetrated in the amorphous or disordered regions of MCC. Initially, the water molecules were relaxed using the conjugate gradient algorithm method before the water-cellulose system was further relaxed using the same method. The numbers of water molecule used in the simulations were 100, 400, 800, and 1400, which is equivalent to water contents of 1.8, 7.3, 14.6, and 25.5 w/w%, respectively. Subsequently, thermal energy was introduced in order to increase the temperature from 0.1 to 300 K within 50 ps. The system was equilibrated at 300 K and 1 atm pressure for 1 ns with a timestep of 0.5 fs. In all the MD simulations, the temperature and pressure were controlled by using Nose-Hoover thermostat and barostat, respectively. The computations were performed on the GPU Supercomputing System Mole-8.5 at the Chinese Academy of Sciences.
The climbing-image nudged elastic band (CI-NEB) method 35 implemented in LAMMPS was used to find the minimum energy path and the energy barrier of water when diffusing in the MCC from one energy basin to another. In the CI-NEB, the fast inertial relaxation engine optimization method 36 was used with an energy stopping criteria of 2 × 10 −8 . The number of images was 5 and a spring constant of 10 was applied in the CI-NEB calculation. The starting configuration was built on the relaxed 5 × 5 × 5 unit cells of MCC (minimized MCC structure at 0 K without performing MD simulations) as described previously. One water molecule was then inserted randomly in the void space between the polymer chains and was perturbed by thermal energy at 300 K to find the local minima; the MCC was not perturbed during the process to obtain a clear potential energy surface. The resulting water-MCC configuration with the minimum energy level was used as the initial and final points of transition path in the CI-NEB calculation.

RESULTS AND DISCUSSION
The density of the MCC model at 300 K is 1606 kg m −3 where the mean measured value is 1540 kg m −3 for samples with about a 61-64% degree of crystallinity 37 while the density of a perfect cellulose crystal is between 1582 and 1599 kg m −3 . 38 The calculated Young's modulus of the MCC model is 13.45 GPa while experimental values are in the range of 9.2-11.1 GPa. 32,39,40 Therefore, in terms of both density and elastic modulus, the current MCC model can be considered as a reasonable numerical analogy of the real material.

Effect of Moisture on Swelling Behavior
The specific volume of MCC with various water moisture contents as calculated by the simulations and from experimental measurements is shown in Figure 2. Khan et al. reported that for MCC (Avicel PH-101; FMC Corporation) particles, the specific volume increases as a function of water content from around 0.5 up to 9.5% 8 ; the densities of these MCC particles were measured with an air pycnometer as a basis for calculating the specific volumes. Sun 38 presented a similar trend for MCC (Avicel PH-102; FMC Corporation) samples with moisture contents in the range of 3-9%; this involved measuring the density of MCC tablets at each moisture content and extrapolating the tablet density to that at a zero porosity. As presented in Figure 2, the calculated specific volumes of the MCC model are is in close agreement with the experimental values. Figure 3 presents the changes in energies of the MCC as a function of moisture content. The energies are normalized for the clarity, by normalizing the values by the respective maximum value. The energies decrease with increasing moisture content since there are more water molecules available in MCC to change its structure, as indicated by the increase of the specific volume (Fig. 2). In the case of the hydrogen bonding energy for MCC only [ Fig. 4(a)], the value decreases which are related to the breaking of some hydrogen bonds within the MCC structure due to the addition of water molecules. The normalized number of hydrogen bond in the MCC decreases as moisture content increases [ Fig. 4(b)]. The numbers of hydrogen bond are normalized by dividing the values by the maximum value. Here, the hydrogen bond is defined by the distance between the hydrogen atom and an acceptor atom being <0.245 nm based on previous studies. 41,42 The numbers of hydrogen bonds per unit MCC cell are 101.96, 96.54, 92.59, and 85.74 for the moisture content of 25.5, 14.6, 7.3, and 1.8 w/w%, respectively. Decreasing number of hydrogen bonds in the MCC contributes to the expansion of the MCC dimensions since the polymer chains are less restricted, which facilitates their motion under the action of the available thermal energy. Increasing the molecular mobility also influences the stiffness and will be discussed in the next section. At greater moisture contents, the presence of additional water molecules has a more significant effect of increasing the hydrogen bonding energy of the system (MCC + water) than that arising from bond scission within the MCC [ Fig. 4(a)].
The valence angle and, more significantly, the torsion angle energies decrease as the moisture content increases (Fig. 3). The torsional energy is related to conformational variations, which depend principally on the rotations about the MCC chains. The reduction is due to the water molecules present in the MCC that accumulates internal pressure and repel the adjoining MCC chains, and thus expanding the volume of the overall system. The changes in the hydrogen bonding and the resulting changes in the torsional energies of the system as a function of the moisture content are critical factors in the plasticization of MCC by water molecules that causes a greater molecular disorder swelling.

Effect of Moisture on the Young's Modulus
The effect of moisture on the Young's modulus of compressed MCC (Avicel PH-101; FMC Corporation) tablets has been reported previously, 8,32 where the modulus was measured using three point bending test of compacted MCC particles at various  moisture contents. The Young's modulus decreased with increasing water content as presented in Figure 5. The calculated values from the simulations, also depicted in Figure 5, are in a close agreement with the experimental values. The Young's modulus from the simulations was calculated from the gradient of linear fit for stress-strain curve. Further detail on procedure to calculate the Young's modulus can be found elsewhere. 18 Swelling caused by moisture, as discussed in the previous section, may transform polymers from brittle solids to soft and rubbery ones. Each polymer chain will have water molecules in its vicinity that reduces the intermolecular forces between the chains and increases the molecular mobility. As presented in Figure 4(a), number of hydrogen bonds in MCC decreases with moisture content, thus increases the molecular mobility of polymer chains. As a result, the mean displacement of MCC molecules during a time period of 100 ps is computed to be 0.07, 0.11, 0.19, and 0.29 nm for the moisture content of 1.8, 7.3, 14.6, and 25.5 w/w%, respectively. The increasing displacements indicate increasing of the molecular mobility with moisture content that leads to a softening or plasticization of the MCC. Figure 6 presents the longest displacements of carbon atoms in the polymer chains during the observed time period for each moisture content. Thus, increasing water content softens cellulose in a similar way to increasing temperature 11 but involves a different mechanism. Increasing water content decreases the number of hydrogen bond [ Fig. 4(b)] thus increases molecular mobility leading to reducing the stiffness. Shishehbor et al. 43 found that hydrogen bonds contribute significantly to the stiffness of cellulose nanocrystals calculated perpendicularly to the polymer chains.
In addition, the reduction in the torsional energy (Fig. 3) indicates that the moisture facilitates the rotation of the polymer chains that allow them to align with the direction of the applied stress. Consequently, the chains can be strained more easily thus reducing the Young's modulus. Moreover, increasing the specific volume of the MCC (Fig. 2) indicates that the free volume also increases with increasing moisture content, which also facilitates the molecular mobility leading to a reduction in the Young's modulus. 44 The free volume also affects the mobility of water molecules that will be discussed in the next section. Thus, it may be concluded that there are both bond scission and free volume contributions to the plasticization of MCC by water.

Mobility of Water in MCC
The mean square displacement (MSD) of water molecules in MCC was calculated in the MD simulations to determine the diffusion coefficient, as follows: where r(t) is the position of a molecule at time t and r(0) is the initial position of a molecule. The diffusion coefficient, D, was calculated by a linear fit of the MSD as a function of time, often called the Einstein's relation, as follows: The value of D increases with increasing moisture content (Fig. 7), which may be attributed to the increase in the free volume. To the authors' knowledge, there are not any published experimental diffusion data for water in MCC. There is one, however, for water in kraft pulp fibers, which consist of almost pure cellulose fibers. The calculated values from NMR experiments have been presented in ref. 45 ; they are of the order of 10 −10 m 2 s −1 , depending on the moisture content and diffusion time. Figure 7 shows that the diffusion coefficients that were calculated by the simulations are in a close agreement with those  measured. Moreover, the diffusion coefficient for "bound" water (water interacting with the component cell wall polymers) in bacterial cellulose, which is chemically identical with plant cellulose, produced from Gluconacetobacter sp. bacteria with 22% moisture content as calculated from quasi-elastic neutron scattering has been found to be 0.86 × 10 −10 m 2 s −1 at 250 K and increases to 1.77 × 10 −10 m 2 s −1 at 265 K, 46 which is in a similar range to those of kraft pulp fibers and the current simulations. In addition, the self-diffusion coefficient of water obtained from pulsed field gradient spin-echo NMR at 298.15 K is 2.30 × 10 −9 m 2 s −1 , 47 showing that, as expected, the diffusion coefficient of water in MCC is significantly smaller than that for bulk water.
One of the theoretical explanations of diffusion is the so-called free-volume theory, which states that a particle moves by hopping from one void or hole to another. 48 The holes are created by the Brownian movement and thermal fluctuations such that the free volume in the system is being continuously redistributed. The term free volume does not always have the same definition 49 and the value cannot be determined by methods that are independent of viscosity and self-diffusion measurements. 48 One of the definitions relates the free volume to the Voronoi volume. 50 Particles in a liquid may be assumed to exist in their enclosures, which are the Voronoi polyhedral. Adjacent Voronoi polyhedral can enlarge and contract to various degrees while maintaining their total volume. The volume that can be swapped between the adjacent Voronoi polyhedral, where their total energy is conserved, has been defined as the free volume. When the Voronoi volume increases, the free volume available for the particle also increases. 51 The Voronoi volumes were calculated using the Voronoi analysis function in OVITO, 52 it computes the Voronoi cell for each particle individually, as a polyhedron surrounding the particle. The Water in the neighborhood of larger structures can be classified according to the mobility (low, intermediate, and high) of their molecules investigated using the NMR method. [55][56][57] In the case of water in wood, where cellulose is the main part of the wood fibers, it is difficult to differentiate water with low and intermediate mobilities since there may not be a significant difference between their spin-spin relaxation times (T 2 ). 58 From the T 2 distributions, water in wood is classified as free water with long T 2 and bound water with short T 2 . Moreover, for water in cellulose ethers, the water molecules can be classified as bound and free water. 59,60 The water molecules interact with the cellulose by hydrogen bonding, reducing their fluctuations and the T 2 decreases as the concentration of water reduces. Examples of trajectories for bound water and free water within the MCC for the current simulations at a moisture content of 7.3 w/w% are depicted in Figure 9. The bound water molecule shows restricted motion and thus indicating that it is tightly bonded to the cellulose, while the free water has more freedom to move from one unit cell of the MCC to another.
The correlation time for water in complex biological systems as measured by NMR can be assumed to be continuously distributed. 61 The correlation time is the characteristic timescale of the molecular motion and is related to the frequency of the randomly tumbling molecules. For example, water molecules in the hydration layer of the solid matter in the muscle are treated as spherical molecules undergoing translational and rotational motions governed by a distribution of the correlation time. 62 The motion of the molecules may fluctuate the magnetic fields and at the resonance frequency can cause relaxation of the nuclei. Thus, the relaxation time is related to the correlation time. The relaxation time T 2 for three wood types with various moisture contents have been shown to be distributed and may have several peaks depending on the type and moisture content. 58 The gamma distribution was used to fit the histogram of the displacement of the water molecules calculated in a period of 100 ps, where the probability density function is as follows: where Γ(.) is the gamma function, a is the shape parameter, and b is the scale parameter. It has been employed to fit the histogram of asymmetric random walk simulations for modeling systems with a preferred spatial direction due to an external force. 63 The density function plot for the fitted distribution is presented in Figure 10. The density function of the Gaussian distribution for the displacement of 100 water molecules in the bulk water simulation is also depicted as a comparison. The mean, variance, and skewness of the distribution are calculated using the following equations: Variance = ab 2 ð5Þ Their values as a function of the normalized moisture content are presented in Figure 10.
Increasing the moisture content increases the mean and variance but results in a reduction in the skewness of the distribution as presented in Figure 11. The mean, variance, and skewness values were normalized by dividing the values by the respective maximum values. Skewness represents the measure of asymmetry of the distribution about the mean. For the fluctuations of a system, the symmetry of the distribution may be broken when an external field introduces a special direction and thus the distribution is non-Gaussian. [64][65][66] The interactions between water molecules and MCC molecules can be considered as an external forces that change the symmetry of the distribution. The water molecules are connected to the MCC molecules by various bonded and nonbonded interactions with different levels of freedom to move and thus have different displacement lengths during the simulations. However, the distribution of displacements in the bulk water simulation is well represented by the Gaussian distribution, which is symmetric and thus has zero skewness. This represents the fluctuations of water molecules without external forces. Increasing the moisture content in the MCC to 100 w/w% is expected to reduce the skewness to almost zero, in which most water molecules will be free water and have more freedom to move like bulk water molecules.
The normalized mean number of hydrogen bonds per water molecule as a function of moisture content for MCC-water and water-water interactions is presented in Figure 12. The normalized values were obtained by dividing the values by the respective maximum values. The mean number of MCCwater bonds per water molecule decreases up to the moisture content of 14.6 w/w% without a significant further change at 25.5%. The decreasing mean number of hydrogen bonds between MCC and water molecules per water molecule with increasing of moisture content contributes to increasing the mobility of the water molecules. The mean hydrogen bonding per water molecule for the interaction between water and epoxy-amine molecules has been modeled and found to decrease with increasing moisture content up to 8 w/w%. 67 It has been observed experimentally that the mobility of water molecules in cellulose 46,68 and in maltose 69 also increases with temperature. Therefore, as mentioned previously, increasing moisture content has a similar effect on water mobility as increasing temperature but by different mechanisms. Increasing temperature increases the kinetic energy of molecules and thus leads to a weakening of the hydrogen bonds or a reduction of the number of hydrogen bonds per water molecule, 70,71 therefore weaken the interaction between water-polymer matrix molecules and increasing mobility. However, increasing moisture content does not increase the kinetic energy. Instead, additional water molecules are linked more strongly to the existing water molecules than to the MCC polymer chains and thus reducing the molecular interactions between the water and MCC molecules that leads to plasticization of the MCC. In this simulations, water molecules that do not form hydrogen bond with MCC are defined as free water. The percentage of free water increases with the moisture content as follows: 2, 5, 8, and 12% for the moisture contents of 1.8, 7.3, 14.6, and 25.5 w/w%, respectively.
The mean number of hydrogen bonds per water molecule between water-water molecules increases with increasing moisture content as presented in Figure 12. By increasing the moisture content, the additional water molecules will more likely be attracted to the existing absorbed water molecules to form a more similar structure to bulk water. The calculated number of hydrogen bonds per water molecule between water-water molecules at the highest moisture content (25.5 w/w%) is 2.2. For lower moisture contents (14.6-1.8 w/w%), the number of hydrogen bonds are 2.1, 1.9, and 1.6, respectively. This value is, as expected, less than the estimated value from the experimental neutron diffraction data of 3.58 for bulk water. 72 Clearly, with increasing moisture content, their molecular arrangements will asymptotically approach that of bulk water. The increase of the mobility of water molecules in cellulose with increasing moisture content has been established using NMR measurements and it was concluded that this is due to extensive water-water contacts. 73 Other simulation studies also predicted that water molecules are more mobile if they are surrounded by other water molecules than if they are in contact with polymer chains, for example, poly(vinyl alcohol) 74 and poly(vinyl pyrrolidone). 75 The calculated total number of hydrogen bonds between MCC-water molecules increases with the moisture content as follows: 343, 1113, 1888, and 3254 for the moisture contents of 1.8, 7.3, 14.6, and 25.5 w/w%, respectively.   Figure 13 presents the displacements of the oxygen atoms of two different water molecules at various moisture contents. The selected oxygen atoms correspond to the water molecules that have the shortest [ Fig. 13(a)] and longest displacements [ Fig. 13 (b)] during the observed time period for each moisture content. The shortest displacements are for the water molecules trapped at binding sites; the displacements are <0.245 nm (hydrogen bond length) as depicted in Figure 13(a). The longest displacements are for the water molecules that are free to move diffusively in the MCC. The displacements continue to increase with time as depicted in Figure 13(b). The displacements are longer with increasing moisture content. As presented in Figure 12, increasing the moisture content results in a decrease in the mean number of hydrogen bonds between MCC and water per water molecule, thus there is an increase in the mean mobility of the water molecules (Fig. 11). Figure 14 presents the trajectory of water molecule moving from one energy basin to another calculated using the CI-NEB method. The related activation energy for water molecule is calculated to be 8.7 kcal mol −1 . The activation energy is equivalent to the amount of thermal energy required to move the water molecule from one binding site to another. An apparent activation energy of MCC with water moisture content of 20.5% was determined to be 7.4 AE 0.3 kcal mol −1 using the NMR technique. 76 Considering that the calculated energy is for one water molecule, the value is reasonable compared to the experimental value. Moreover, the experimental value is an average of activation energy of water molecules at various binding sites which are unlikely to be in identical molecular environments. In addition, for water in cellulose gel (synthesized from sodium cellulose xanthate) with a moisture content of 670% (6.70 g water/g dry gel), an apparent activation energy was determined to be 5.57 kcal mol −1 (23.3 kJ mol −1 ) using the NMR technique. 77 This suggests that increasing moisture content reduces the activation energy as might be expected. While increasing the molecular mobility of MCC by increasing the moisture content ( Fig. 6) results in a decrease in the mean number hydrogen bonds between MCC and water per water molecule (Fig. 12), it leads to a reduction in the activation energy. As a comparison, activation energies for water in the cell walls of Sitka spruce (Picea sitchensis) were calculated to reduce from about 10.74 to 6.13 kcal mol −1 as the moisture content increased from 0 to 30%. 78  Trajectory of a water molecule moving from one energy basin to another calculated using the CI-NEB method. MCC atoms are depicted as gray = carbon, red = oxygen, white = hydrogen. Water atoms are depicted as blue = oxygen, green = hydrogen 1, yellow = hydrogen 2. Blue, green, and yellow lines are the trajectories of oxygen, hydrogen atom 1, and hydrogen atom 2 of water molecule, respectively. [Color figure can be viewed at wileyonlinelibrary.com]

CONCLUSIONS
The hybrid MM-MD simulation method has been performed in order to study the effects of water on the properties of MCC. The specific volume of MCC and the diffusion coefficient of water in MCC increase with increasing moisture content while the Young's modulus decreases. The simulations results are in close agreement with the published experimental data. Water molecules repel the adjoining MCC chains and disrupt the hydrogen bonding, thus increasing the molecular mobility, causing swelling and a reduction in the elastic modulus.
The Voronoi volume, which is related to the free volume, increases with increasing moisture content and will contribute to the increase of the diffusion coefficient as predicted by the freevolume theory of diffusion. The interactions between water and MCC molecules act as external forces that skew the symmetry of the displacement distribution of the water molecules. Moreover, increasing the moisture content increases the mobility of water molecules in the MCC and the number of free water molecules. With additional water molecules, there is relatively less hydrogen bonding per water molecule between MCC and water molecules and thus the mean mobility of the water molecules increases. The accuracy of the simulations is verified further by the calculating the activation energy for a water molecule in MCC that is found to be in a good agreement with the experimental values.