Extended ensemble molecular dynamics study of cellulose I–ethylenediamine complex crystal models: atomistic picture of desorption behaviors of ethylenediamine

Cellulose I crystals swell on exposure to ethylenediamine (EDA) molecules to form a cellulose I–EDA complex, and successive extraction of EDA molecules converts the complex crystalline phase to either original cellulose I or cellulose IIII, depending on the treatment procedure. The present study reports the extended ensemble molecular dynamics (MD) simulation of the cellulose I–EDA complex models. An accelerated MD simulation allows most of the EDA molecules to desorb from the crystal model through a hydrophilic channel between the piles of cellulose chains, one at a time. Migration of a single EDA molecule along the channel is simulated by the adopted steered MD method combined with the umbrella sampling method to evaluate the potential of mean force (PMF) or free energy change on its movement. The PMF continues to increase during the migration of an EDA molecule to give a final PMF value of more than 30 kcal/mol. The PMF profiles are largely lowered by the removal of EDA molecules in the neighboring channels and by the widening of the channel. The former suggests that the EDA desorption cooperates with that in the neighboring channels, and, in the latter case, an EDA migration is efficiently promoted by solvation with water molecules in the expanded channel. We conclude that the atomistic picture of the EDA desorption behaviors observed in the crystal models is applicable to the real crystalline phase.

Abstract Cellulose I crystals swell on exposure to ethylenediamine (EDA) molecules to form a cellulose I-EDA complex, and successive extraction of EDA molecules converts the complex crystalline phase to either original cellulose I or cellulose III I , depending on the treatment procedure. The present study reports the extended ensemble molecular dynamics (MD) simulation of the cellulose I-EDA complex models. An accelerated MD simulation allows most of the EDA molecules to desorb from the crystal model through a hydrophilic channel between the piles of cellulose chains, one at a time. Migration of a single EDA molecule along the channel is simulated by the adopted steered MD method combined with the umbrella sampling method to evaluate the potential of mean force (PMF) or free energy change on its movement. The PMF continues to increase during the migration of an EDA molecule to give a final PMF value of more than 30 kcal/mol. The PMF profiles are largely lowered by the removal of EDA molecules in the neighboring channels and by the widening of the channel. The former suggests that the EDA desorption cooperates with that in the neighboring channels, and, in the latter case, an EDA migration is efficiently promoted by solvation with water molecules in the expanded channel. We conclude that the atomistic picture of the EDA desorption behaviors observed in the crystal models is applicable to the real crystalline phase.

Introduction
Among the versatile properties of cellulosic materials, crystal polymorphism has been attracting much attention in terms of both scientific and technological interests. In particular, cellulose forms complex phases with various amines and liquid ammonia as an intermediate state in the crystalline conversion between the cellulose I b and III I forms. (Chanzy et al. 1987;Wada 2001;Wada et al. 2001Wada et al. , 2004Wada et al. , 2006Wada et al. , 2008Wada et al. , 2009a In comparison with the native cellulose allomorphs, I a and I b , cellulose III I is more susceptible to certain chemical treatments (Clark and Parker 1937;Davis et al. 1943;Pentoney 1966;Perez et al. 2003) and enzymatic degradation (Igarashi et al. 2007(Igarashi et al. , 2011Chundawat et al. 2011). Pretreatment of native cellulosic materials with liquid ammonia and amines, followed by the chemical and enzymatic methods, has been expected to be an efficient method for biomass conversion to produce biofuel. In particular, ammonia fiber expansion (AFEX) pretreatment effectively decreases the crystallinity of the original native cellulose material and reduces the amounts of hemicellulose and lignin (Lau and Dale 2009;Zheng et al. 2009). Information on the underlying molecular events of the crystal conversions, such as rearrangement of cellulose chains, penetration and desorption of guest amines, and exchange of the hydrogen bonding network, will improve the amine treatments for an increased effective biomass conversion. The high-resolution crystal structures of the relevant cellulose allomorphs have exclusively been reported by groups in Tokyo, Grenoble, Los Alamos, and Oak Ridge, using synchrotron fiber diffraction analysis, sometimes combined with nuclear magnetic resonance (NMR) measurement and theoretical calculations (Nishiyama et al. 2002(Nishiyama et al. , 2010Wada et al. 2004Wada et al. , 2009bWada et al. , 2011Sawada et al. 2013). These allomorphs include cellulose I a and I b , cellulose III I , and cellulose Iethylenediamine (EDA) and ammonia complexes. Among them, the crystal structure of the cellulose I-EDA complex includes four cellulose chains in parallel polarity and eight EDA molecules, which are packed in a monoclinic unit cell (a = 4.546 Å , b = 11.330 Å , c = 10.368 Å , and c = 94.017°in the space group P2 1 ) (Wada et al. 2009a). The ab projection structure reveals that cellulose chains are piled by their pyranose faces upon each other along the a-axis and that EDA molecules align along the hydrophilic channel composed by the hydrophilic edges of the piled cellulose chains (see Fig. 1). Intermolecular hydrogen bonds are formed between cellulose chains and EDA molecules (Wada et al. 2009a), which have been exchanged for the direct hydrogen bonds between cellulose chains in the parent allomorphs, cellulose I b and III I (Nishiyama et al. 2002;Wada et al. 2004). A time-resolved X-ray microscopy diffraction measurement detects that no intermediate phase appears during the crystalline conversion, which suggests the successive crystalline conversion from cellulose I b to the I-EDA complex (Nishiyama et al. 2010).
Theoretical calculations using molecular dynamics (MD) and density functional theory (DFT) have often been applied to examine the cellulose crystal structures to define the positions of the hydroxyl hydrogen atoms and, for the complex crystals, to refine the ligand structures. The former method has also revealed the dynamic features of the crystal structures such as exchange in the hydrogen bonding network coupled with rotation of hydroxyl groups. For both methods, periodic systems are mostly adopted; the crystalline phase continues beyond the boundary of a unit cell, usually consisting of several crystal unit cells, so that molecular chains of infinite length can be established. In contrast to the infinite system, French proposed a molecular modeling methodology based on the miniature crystals, or mini-crystals, for carbohydrates, which deals with isolated small crystal models with finite dimensions (French and Dowd 1993;. Intensive hydrogen bonding networks, which are characteristic of most carbohydrate crystals, enable the mini-crystal models to substantially retain the original crystal structures after energy minimization. One can work with the optimized models to define the hydrogen bonding network and thermodynamics properties. We first combined mini-crystal modeling with the MD method to examine the structure dynamics of the b-chitin structure (Yui et al. 2007). The combination was found to be more valid for an extended application to simulate a largescale structural change, which includes a fiber twist of the native cellulose crystal (Matthews et al. 2006;Yui et al. 2006), the decrystallization of crystalline fibers , and the conversion of cellulose III I (Yui and Hayashi 2008;Yui et al. 2010;Uto et al. 2013). In a similar fashion, the finite cellulose chain sheet models, which were extracted from the corresponding cellulose allomorphs, were studied by energy minimization using the DFT method to examine the deformation of the initial sheet forms (Uto et al. 2014;Uto and Yui 2018). We recently reported an MD study of the cellulose I-EDA crystal models in a solution state (Uto et al. 2019). Although our primary objective of the MD calculations was to observe the diffusion and desorption of EDA molecules, it was observed that desorption of EDA molecules was limited to the first and second EDA layers at the early simulation time and that it apparently ceased for the remainder of the 100-ns simulation time. Despite the small sizes of the crystal models, which consist of 48 cellulose chains with a degree of polymerization (DP) of 40, compared with the real crystalline phase, EDA desorption would be an excessively slow event under the ordinary simulation time. Similarly, an MD study aimed at the complexation of a native cellulose crystal with ammonia molecules was carried out only to observe their penetration into the surface cellulose layers ). In the present study, MD simulations of the cellulose I-EDA complex models have been carried out, this time by adopting the extended ensemble theories. We have successfully achieved almost complete desorption of the EDA molecules and estimated the free energy change on diffusion of an EDA molecule in the crystal model.

Accelerated molecular dynamics
Accelerated molecular dynamics (aMD) is an enhanced sampling technique that flattens the potential energy surface by introducing a non-negative boost potential, by which one can observe the longtime behaviors of simulation systems with limited computational time (Hamelberg et al. 2004(Hamelberg et al. , 2007. Figure 1a shows the ab base planes of the two crystal models, which consist of 24 and 48 cellulose chains with DP = 40, denoted as the 40mer 9 24 and 40mer 9 48 models, respectively. Both crystal models share the (010) lattice planes of the same length, whereas the length of the (100) lattice plane of the former is doubled in the latter. The 40mer 9 24 was solvated in 50,000 TIP3P water molecules (Jorgensen 1981;Jorgensen et al. 1983) containing 546 EDA molecules, and in 90,000 water molecules and 1,365 EDA molecules for the 40mer 9 48 model. The present aMD simulations adopted two acceleration levels, termed as a dual boost, which applied a boost potential to both dihedral and total potentials. The aMD parameters, the threshold boost energy, E, and the acceleration factor, a, are defined for each potential according to the Amber Tutorial (Pierce et al. 2012;Salomon-Ferrer et al. 2013b;Miao et al. 2014), following: where N residues is the number of solute residues including EDA, N atoms is the total number of atoms, and V dihedral_avrg and V total_avrg are the average dihedral and total potential energies from conventional MD calculations, respectively. As described in the next section, many EDA molecules are released from the crystal models in the early stage of the aMD calculations, which fills a solvent box with desorbed EDA molecules that should be removed from the estimation of the aMD parameters. The present aMD simulations were carried out by the multi-cycle procedure involving a solvent displacement step accompanied by re-estimation of the aMD parameters, following: 1. Equilibration step. The simulation system, the cellulose I-EDA complex crystal model with water molecules, is equilibrated by the standard procedure described below. The final part of the procedure involves the conventional MD run for 1 ns for estimation of the aMD parameter. 2. aMD calculations. The equilibrated system is subject to aMD calculations for 10 ns with the aMD parameters determined based on the conventional MD trajectory performed in the previous cycle, using the NPT ensemble at a constant temperature of 300 K or 370 K and pressure of 1 bar. The number of EDA molecules remaining in the crystal model is counted for each 1-ns aMD run. 3. Displacement of the solvent box. A cellulose crystal accommodating the interior EDA molecules and water molecules, and the water molecules attached to the crystal surface are extracted from the final aMD structure. A new solvate box system is built for the crystal model extracted, where the number of water molecules is adjusted to be same as the original number. The new simulation system is then moved to the equilibration step of the next cycle.
Adaptive steered molecular dynamics Steered molecular dynamics (SMD) uses a pseudo particle that applies a steering force to analyze the physical movement of atoms and molecules along a reaction coordinate at a particular velocity (Park et al. 2003;Park and Schulten 2004). SMD, in combination with the Jarzynski's non-equilibrium work relation (Jarzynski 1997), predicts the potential of mean force (PMF), which is derived from the exponential average of numbers of non-equilibrium works. Adaptive steered molecular dynamics (ASMD) is an improved method of SMD for effective sampling, where the reaction coordinate is divided into smaller segments in which the work distribution becomes narrower (Ozer et al. 2010(Ozer et al. , 2012. The 16 cello-octamers (8mer 9 16 model) model was adopted in the ASMD shown in Fig. 1b. From the crystal surface, four EDA molecules were removed from the central hydrophilic channel to produce a hole for an EDA migration path. The model was termed as the one-hole model and, in a similar fashion, the two-and three-hole models were composed by digging an additional hole on one side or both sides of the original hole, respectively. The one-hole model was additionally modified by expanding the space of the migration path, such that both crystal model halves (the 8mer 9 8 units) were slid outward along the b-axis for a given distance. The model was termed as the gap model. All the crystal models were solvated in 11,000 to 15,748 TIP3P water molecules (Jorgensen 1981;Jorgensen et al. 1983) and the one-, two-, and three-hole models contained 45, 40, and 35 EDA molecules, respectively. The EDA molecule at the bottom of the hole, which occupied the fifth EDA layer from the crystal surface, was assigned to be the target molecule to be pulled along the hydrophilic channel. The reaction coordinate was defined to connect the center of the target EDA with that of the opening of the hole at the crystal surface.
A harmonic potential with a spring constant of 71.5 kcal/(molÁÅ 2 ) was imposed on the center of the target EDA, which traveled at the constant pulling velocity of 0.5 Å /ns for 4 Å , which corresponded to the spacing of the cellulose layers. PMF, or the free energy difference, was determined from the ensemble average of the work, W, done by 50 SMD runs by the following Jarzynski's equality (Jarzynski 1997): T is the temperature. At the switching point of the SMD stages, a single trajectory was selected according to that whose work value was the closest to the PMF of the current SMD stage, and the structure at the final frame of the selected trajectory was used as the starting structure of the next SMD stage. The four SMD stages covered the overall distance of the migration path of 16 Å in the hydrophilic channel; meanwhile, the target EDA traveled through the four cellulose layers until it reached the crystal surface. The SMD calculations were performed using the NVT ensemble at a constant temperature of 300 K. Throughout the ASMD calculations, the positions of all backbone atoms of a cellulose chain, C1, C2, C3, C4, C5, and O5, and of the two carbon atoms of the EDA molecules, except those of the target EDA, were weakly constrained with a force constant of 0.1 kcal/(molÁÅ 2 ) to prevent the crystal model from irreversible distortion resulting from EDA movement.

Umbrella sampling method
In contrast to the SMD and ASMD methods with a non-equilibrium basis, an umbrella sampling (US) simulation is a classified equilibrium molecular dynamics method Valleau 1974, 1977), although it is also used to calculate PMF or a free energy change along a given reaction path, where biased potentials are applied to an atom or molecule at certain closely spaced intervals along the path. In the present US simulations, the PMF trajectory obtained by the SMD calculations was adopted as a reaction coordinate of which the overall length, 16 Å , was divided into 0.01 Å intervals, which gave 1,600 umbrella windows. At each interval, a harmonic potential with a force constant of 100 kcal/(molÁÅ 2 ) was applied to the target EDA molecule and the NPT production run was performed for 1 ns at a constant temperature of 300 K and pressure of 1 bar. As was performed in the present ASMD, all of the backbone atoms of cellulose and the carbon atoms of all EDA molecules, except the target EDA, were positionally constrained with a force constant of 0.1 kcal/(molÁÅ 2 ). The biased EDA distributions obtained from the production runs were unbiased and combined using the weighted histogram analysis method (WHAM) to obtain a PMF profile (Kumar et al. 1992(Kumar et al. , 1995Roux 1995).

General MD procedures and parameters
All energy minimizations and MD calculations of the solvated cellulose I-EDA complex systems were carried out using the PMEMD and PMEMD.CUDA modules (Le Grand et al. 2013;Salomon-Ferrer et al. 2013a) of the Amber 18 program package (Case et al. 2018) with an NVIDIA Kepler GPU system. The equilibration procedure of the system was composed by two-step energy minimizations: optimization of water configurations and full optimization of the whole system, and the gradual heating of the system from 100 to 300 K for 500 ps using the NVT ensemble. The density of the system was then equilibrated by the NPT run at a constant temperature of 300 K and pressure of 1 bar for, typically, 1 ns. In the first energy minimization and the heating runs, the positions of the solute molecules, cellulose chains, and EDA molecules were constrained with a force constant of 50 kcal/(molÁÅ 2 ). In the ASMD calculations, the density equilibrating NPT runs were extended to 50 ns until the holes in the hydrophilic channels had filled with waters with the positional constraint. The temperature and pressure of the system were regulated by a Langevin thermostat (Hoover et al. 1982;Evans 1983) with a collision frequency of 5 ps -1 and a Berendsen barostat (Berendsen et al. 1984), respectively. A time step of 2.0 fs, enabled by the SHAKE option (Ryckaert et al. 1977), was employed to ensure rigid hydrogen atoms. The particle mesh Ewald (PME) method (Essmann et al. 1995) was adopted for longrange interactions and the van der Waals interactions were cut off at 10 Å . The three extended ensemble MD options for aMD, SMD, and US, implemented by the AMBER package, were used. Cellulose molecules were described by the GLY-CAM06 force field (Kirschner et al. 2008) and EDA molecules were modeled using the revised GAFF parameters (Wang et al. 2004) developed by us (Uto et al. 2019).
The MD trajectories were analyzed using the CPPTRAJ module (Roe and Cheatham 2013) of the AmberTools 19 and, for a visual analysis, VMD 1.9.3 (Humphrey et al. 1996). US-PMF profiles were calculated using the WHAM code ver. 2.0.10 written by Grossfield Lab. (Grossfield 2018). The molecular graphics software, PyMOL ver. 2.5.0 (Schrödinger LLC. 2021), was also used for visualizing the system.

Results and discussion
Desorption behaviors of EDA molecules Figure 2a shows changes in the amounts of remaining EDA molecules in percent with respect to simulation time. Note that a 1 ns trajectory of every 11 ns trajectories corresponds to the conventional MD for sampling thermodynamics values, as described in the Computational methods section. Both the 40mer 9 24 and 40mer 9 48 models finally released more than 80% of the initial EDA molecules at 300 K. Both PMF profiles almost coincided with each other having a similar desorption rate of EDA. At the elevated temperature of 370 K, the desorption rate remarkably increased in the 40mer 9 24 crystal model and the PMF finally reached a constant value of approximately 10%. This was consistent with the observation that EDA molecules in the complex crystal became mobile at a temperature larger than 70°C (Wada et al. 2008). Figure 2b and c depicts the cross-views of the sliced structures of the 40mer 9 48 model and its overall structure obtained from the final simulation time, respectively. The respective EDA desorption behaviors are given in Movies S1 and S2. As one can clearly see in Movie S1, EDA molecules diffused along the hydrophilic channel between the cellulose chain piles, one at a time, and then moved away from the crystal model, rather than simultaneous diffusion involving more than one EDA molecule. The vacant space, once occupied by an EDA molecule, was quickly filled by water molecules that diffused from the bulk. As the hydrophilic channel lost most of intercalated EDA molecules, the cellulose chain piles collapsed and enclosed the remaining EDA molecules in the disintegrated crystal model. This was obviously the cause of the deacceleration of EDA desorption. As shown in Fig. 2a, another ASMD simulation for the 40mer 9 24 crystal model with the constrained cellulose chain structure exhibited almost complete desorption, because the hydrophilic channels were retained throughout the simulation. In the real crystalline phase, such a noticeable cellulose chain collapse is unlikely, owing to less mobile cellulose chains with a large DP, so that the crystalline phase can successfully convert to another cellulose crystal form after complete EDA desorption. In the present crystal model, EDA molecules most quickly disappeared in the central channel (see Movie S2), which led to dividing the crystal model into halves along the longitudinal (fiber) direction (see Movie S1). As a result, the overall desorption behavior can be represented by those of the crystal model halves, as also suggested by the similar desorption curves between the 40mer 9 24 and 40mer 9 48 crystal models. We have reported in a previous MD study that the cellulose I-EDA crystal model developed a right-handed twist as it loses the guest EDA molecules (Uto et al. 2019). Similarly, the present halves of the crystal model twisted and moderately intertwined with each other. The crystalline deformation involving fiber twisting and entangling likely takes place in the real crystal fibers, but on a much larger scale, which results in the degradation of crystallinity and, consequently, a decrease in the dimensions of the single-crystalline phase thorough the crystal conversion. In comparison with a native cellulose crystal, the enhanced accessibility and chemical reactivity of the cellulose III I crystal have been attributed to not only the physical nature of the latter form but also crystal degradation resulting from amine treatments. Salminen et al. reported the solid-state polymorphic transition of cellulose nanocrystals immobilized on a chemically modified silicon surface from cellulose I to III I through its EDA complex (Salminen et al. 2017). The width of the cellulose III I nanocrystals decreased to half the original size of the cellulose I nanocrystals, whereas their length remained unchanged. It should be noted that they also reported no change in crystal dimensions of the cellulose nanocrystals when treated in a bulk dispersion. Nishiyama et al. studied the crystalline conversion from the cellulose I to cellulose I-EDA complex phases by time-resolved X-ray diffraction measurements, which suggested that the lateral width dimensions were limited to half of the cellulose I phase as a result of penetration of EDA molecules (Nishiyama et al. 2010). On the basis of the present aMD results, we wish to propose that crystal thinning can also take place on removal of EDA molecules even in the dispersion state.

Diffusion behaviors of EDA molecule in the crystal models
The ASMD simulations were intended to simulate the diffusion of a single EDA molecule in the hydrophilic channel of the complex crystal, which is based on the discontinuous desorption of EDA molecules observed in the preceding aMD simulations. The time-course changes in PMF, with respect to EDA movement, for the one-, two-, and three-hole models are shown in Fig. 3. In the one-hole model, the PMF continued to increase with the movement of the target EDA molecule and reached an unexpected, large amount of endothermic free energy of more than 30 kcal/mol Fig. 2 a Change in amount of residual EDA in the crystal model with respect to simulation time. b Time course of EDA desorption viewed from the sliced cross-sections of the 40mer 9 48 model. c The sliced cross-section view and overall structure of the 40mer 9 48 models at the final time frame of 495 ns after the EDA molecule passed through the four cellulose chain layers. The PMF profile exhibited the four peaks being barely appreciable, where each peak corresponded to the position of a cellulose chain layer. Obviously, spontaneous EDA desorption would be impossible with the endothermic free energy change. Therefore, the PMF of the one-hole model was reevaluated by applying the US method to the PMF trajectory derived from the ASMD simulation as a reaction coordinate. Figure 3 shows that both the ASMD-PMF and US-PMF profiles, initially overlapping, deviated to each other with a difference of approximately 5 kcal/mol near a 1.2 nm movement, while the latter exhibited a smoother curve. For further validation of the PMF of the one-hole model, ASMD simulations with a slower pulling speed of 0.2 Å /ns, combined with the US method, were also carried out only to reproduce similar PMF profiles ( Figure S1).
The first US-MD trajectory was examined to analyze the spatial directions and conformations of the target EDA molecule by monitoring the changes in direction of the C-C bond vector and its rotation ( Figures S2 and S3). The C-C rotations, defined by the N-C-C-N atom sequence, mostly preferred to be around the ± gauche conformations. Therefore, the apparent rotational motions of the target EDA molecule observed in the visual inspection of the trajectory may reflect rotations of the molecule as a whole, rather than the C-C bond rotation. The directional change in the C-C bond vector with respect to the initial direction of the vector in the crystal structure somewhat diverged at the EDA positions from 0.5 nm to 1.0 nm and, otherwise, they were restricted to certain angles.
In contrast to the PMF profiles of the one-hole model, those of the two-and three-hole models, shown in Fig. 3, displayed moderate increases in the baselines with final values of a few kcal/mol, accompanied by four rounded peaks. The free energy changes indicated that EDA molecules can surmount the energy landscapes by the thermal motions at an ambient temperature. On the basis of these findings, we propose a plausible desorption scheme applicable to the real crystalline phase as follows. Once an initial EDA hole appears on the crystal surface after releasing several EDA molecules, those EDA molecules located in the adjacent hydrophilic channels can cooperatively migrate and diffuse into the bulk water. An alternative scheme would be that EDA desorption is initiated from the interfacial regions with amorphous or paracrystalline phases, or the terminals of crystal fibers where EDA molecules are loosely bound. In fact, EDA desorption apparently proceeds from both the terminals along the fiber axis in the aMD simulation shown in Movie S2.
We performed supplemental ASMD calculations for the one-hole model with the modified potential functions to confirm the above discussions. Figure S4 compares the PMF profiles obtained by the ASMD calculations adopting the modified potential functions with no torsional potential about the C-C and N-C bonds of an EDA molecule and by those implemented with the various dielectric constant values of 1.5, 2.0, and 4.0, along with the original PMF curve of the onehole model given in Fig. 3. The PMF with no torsional potential of an EDA molecule was lowered from the original one, although both exhibited similar appearances in their profiles. The amount of decrease was moderately limited despite the conformation of the EDA molecule being allowed to change freely. More notable changes were observed as the dielectric constant value was increased, which resulted in not only effectively suppressing the PMF level but also changing the shapes of the PMF profile to be similar to those of the two-and three-hole model, as represented by the rounded peaks with moderate or gentle slopes. In the original crystal structure, EDA molecules PMF profiles with respect to the pulling distance of the target EDA molecule obtained from ASMD simulations for the one-, two-, and three-hole models and combined MD simulation using the US method for the one-hole model participate in the intermolecular hydrogen bonding network with O2, O3, and O6 hydroxyl groups of the cellulose chains along the a-axis direction. It is unlikely that, in the one-hole model, the unfavorable endothermic PMF is caused by breaking the intermolecular hydrogen bonds involving the target EDA, because it can form another hydrogen bond in the next cellulose chain layer when it moves. A high frequency of breaking and formation of hydrogen bonds caused by thermal motions of the hydroxyl groups was suggested in the real cellulose I-EDA complex (Sawada et al. 2013). Instead, a certain dragging force against the movement of the target EDA may arise in the one-hole channel, such that EDA migration successively generates an endothermic free energy, although we do not have a straightforward interpretation for the origin of the dragging force at present. The present ASMD simulation constrained the positions of the pyranose backbones and carbon atoms of EDA molecules, except the target EDA molecule, while allowing the exocyclic and amino groups to rotate. The functional groups can more frequently rotate as they become liberated from the intermolecular hydrogen bonds, which is introduced either by removal of adjacent EDA molecules in the two-and three-hole models or attenuation of electrostatic interactions introduced by larger dielectric values. Consequently, the target EDA molecule can find a lower energy path with a higher probability, especially as it moves through the cellulose chain layers. This may be applicable to lowering the PMF in the simulation with no torsional potential of an EDA molecule, although the effect is relatively restrictive compared with that adapting larger dielectric values.
While the introduction of the multi-hole models was to modify the crystal models along the fiber axis, the gap models were examined for the effect of expanding the crystal models along the b-axis (see Fig. 1). Figure 4 compares the PMF profiles of the five gap models with different gap spaces and that of the original model shown in Fig. 3, where the ASMD calculations were performed only for the 1st step with a pulling distance of 0.4 nm. The PMF curves drastically decreased as the gap space increased from 0.4 nm to 0.6 nm. The final PMF value of the latter model indicated that spontaneous diffusion of EDA molecules can occur from the interior of the crystal model. Slight exothermic free energy changes were detected in the PMF curves near a pulling distance of 0.03 nm for the three gap models with a gap space of 0.1, 0.2, or 0.3 nm, where positional relaxation of the target EDA molecule occurs as a result of the developing attractive interactions with the surrounding cellulose chains. However, the target EDA molecules in the 0.4 and 0.6 nm gap models, which exhibited no relaxations, were solvated with water during the preceding equilibration process. The results suggest that water solvation of an EDA molecule plays an important role in its spontaneous departure from the original, crystalline position, which is consistent with the discontinuous diffusion observed in the aMD simulations. It is likely that diffusion of EDA molecules is also promoted by their solvation in the real crystalline phase.

Conclusions
We have reported desorption behaviors of guest EDA molecules from the cellulose I-EDA complex crystal models studied by adopting a combination of three extended ensemble MD methods. EDA molecules migrate along the hydrophilic channel between the cellulose chain piles on the (010) lattice plane. The same EDA migration path is proposed experimentally on complexation of native cellulose fibers with EDA (Nishiyama et al. 2010). The diffusion behavior of Original: the PMF profile of the original one-hole model EDA molecules appears to be discontinuous; they migrate in the channel one at a time. A vacancy in the channel emerging after EDA diffusion is quickly filled by water molecules. This obviously suggests that a departure of an EDA molecule from the crystal position is, at least partly, triggered by its solvation. Following complete desolvation and rearrangement of cellulose chains, crystalline conversion is allowed. The scheme may reasonably be consistent with the observation that the destiny of the crystalline phase depends on the solvent used to rinse the cellulose I-EDA complex specimen; cellulose I by hot-water treatment and cellulose III I by methanol or ethanol treatment (Loeb and Segal 1955;Lokhande et al. 1977;Roche and Chanzy 1981;Sugiyama et al. 1987). It has been also suggested by the final structure of the aMD crystal model that drastic desorption of EDA may lead to crystal degradation and fiber thinning in the real complex crystals. Inspecting the final aMD structures of the crystal models, it was found that residual EDA molecules dispersed over the hydrophilic channels, which prevented cellulose chains from forming an extended, hydrogen bonding chain sheet. No transformed region involving either cellulose I b or III I feature was detected. The fiber twist of the final crystal models, as shown in Fig. 3c, raised chain sliding along the fiber axis in individual cellulose chain piles to give a chain staggering from the original non-staggering of the cellulose I-EDA crystal structure (Wada et al. 2009a). All of the chain piles displayed one-way staggering except the one chain pile in the 40mer 9 24model which involved an alternative chain staggering of the cellulose I b feature (Nishiyama et al. 2002).
The present MD simulations have also suggested that while a large endothermic free energy is required for EDA diffusion under the regular crystalline environment, EDA molecules readily cooperatively migrate in the hydrophilic channel once the adjacent channels have discharged their own EDA molecules along with the loss of intermolecular hydrogen bonds. Diffusion of EDA molecules may be promoted by the increased rotational frequencies of the exocyclic groups of cellulose chains, which are free from restrictions of intermolecular hydrogen bonds. Naturally, enlarging the space of the hydrophilic channel would be a direct solution to accelerate EDA diffusion and desorption. However, the space needs to be expanded so that it is large enough to afford solvation of EDA molecules for spontaneous diffusion of them from the crystal interior.
One may be critical of applying the present simulation results to describe EDA desorption of the real crystalline phase, especially those obtained from the ASMD simulations where the crystal models essentially retain the regular crystalline framework, in contrast to the real crystalline phase that possibly involves disordered features on EDA desorption. However, disordered or para-crystal modeling would be a much more difficult issue; there is neither applicable theory nor experimental evidence to describe a representative configuration of a disordered cellulose crystal with atomistic detail. In the present MD study, the various crystal models have been assessed by combination of the three extended ensemble approaches. The results are reasonably correlated to each other and lead us to the specific conclusions described above, some of which are in accordance with the experimental results. The present simulation methodology can be applicable to similar complex systems, such as a cellulose I-ammonia complex crystal, which is currently underway.