Biomacromolecules 2008, 9, 3133–3140
3133
Neutron Crystallography, Molecular Dynamics, and Quantum Mechanics Studies of the Nature of Hydrogen Bonding in Cellulose Iβ Yoshiharu Nishiyama,† Glenn P. Johnson,‡ Alfred D. French,‡ V. Trevor Forsyth,§,| and Paul Langan*,⊥ Centre de Recherches sur les Macromole´cules Ve´ge´tales of CNRS, affiliated with the Joseph Fourier University of Grenoble, BP 53, 38041 Grenoble Cedex 9, France, Southern Regional Research Center, United States Department of Agriculture, 1100 Robert E. Lee Boulevard, New Orleans, Louisiana 70124, Partnership for Structural Biology, Institute Laue Langevin, 6 Rue Jules Horowitz, F-38042 Grenoble, France, EPSAM/ISTM, Keele University, Keele, Staffordshire ST5 5BG, England, and Bioscience Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545 Received July 2, 2008; Revised Manuscript Received August 29, 2008
In the crystal structure of cellulose Iβ, disordered hydrogen bonding can be represented by the average of two mutually exclusive hydrogen bonding schemes that have been designated A and B. An unanswered question is whether A and B interconvert dynamically, or whether they are static but present in different regions of the microfibril (giving temporally or a spatially averaged structures, respectively). We have used neutron crystallographic techniques to determine the occupancies of A and B at 295 and 15 K, quantum mechanical calculations to compare the energies of A and B, and molecular dynamics calculations to look at the stability of A. Microfibrils are found to have most chains arranged in a crystalline Iβ structure with hydrogen bonding scheme A. Smaller regions of static disorder exist, perhaps at defects within or between crystalline domains in which the hydrogen bonding is complex but with certain features that are found in B.
Introduction Microfibrils of naturally occurring cellulose can consist of two crystal forms, IR and Iβ.1 Cell walls of some algae and bacteria are rich in IR, whereas cotton, wood, and ramie fibers are predominantly composed of Iβ.2,3 Using a combination of X-ray and neutron crystallographic techniques, we have determined the crystal structures and hydrogen (H) bonding arrangements in both allomorphs to atomic resolution (1 Å).4,5 IR corresponds to a triclinic P1 unit cell with one chain that has two similar but geometrically nonidentical adjacent glucose residues related by pseudo two-fold screw-axis symmetry. Iβ has a monoclinic P21 unit cell containing two conformationally distinct chains, which are referred to as the corner and center chains. Each Iβ chain has a two-fold screw axis that requires the adjacent glucosyl residues in that chain to be identical. The cellulose chains in both allomorphs are arranged parallelup6 and edge-to-edge, making flat sheets that are held together by H bonds. There are no strong H bonds between the stacked sheets. The major difference between IR and Iβ is the pattern of the stagger of these sheets in the chain direction. In both forms of cellulose I, all O6 primary alcohol groups (commonly called hydroxymethyl groups) are in tg7conformations. The O3 secondary alcohol group of each residue donates its proton to the O5 ring atom of a neighboring residue to form an intrachain H bond. The O2 secondary alcohol and O6 hydroxymethyl groups are involved in both intrachain and interchain H bonds within * To whom correspondence should be addressed. Tel: 505-665-8125. Fax: 505-665-3024. E-mail:
[email protected]. † Centre de Recherches sur les Macromole´cules Ve´ge´tales of CNRS. ‡ U.S. Department of Agriculture. § Institute Laue Langevin. | Keele University. ⊥ Los Alamos National Laboratory.
the sheets. A surprising result from our earlier diffraction studies is that the H-bonding arrangements involving O2 and O6 are disordered.4,5 That conclusion is based on the neutron scattering density maps that show two locations for each H atom covalently bound to O2 or O6 of the center chain and O6 of the corner chain. Together, the two locations for each atom have a total neutron scattering density that is expected for a proton in a single location, a disordered condition known as fractional occupancy. Two mutually exclusive H-bonding networks can be drawn on the basis of the alternative locations for the protons of the center chain, as illustrated in Figure 1. One (designated scheme A) occurs when the proton on O2 is donated to make an intrachain H bond to O6. In turn, O6 donates its proton to an interchain bifurcated H bond with O3 and O2. In the other (scheme B), O6 donates its proton to an intrachain H bond with O2, and then O2 donates to an interchain H bond with O6. Thus, the sense of the donor-acceptor-donor network is reversed in B compared with A, with the exception that the bifurcated bond to O3 in scheme A has no counterpart in scheme B. Representing the alternative locations of H atoms covalently bound to O2 and O6 as two mutually exclusive H-bonding networks is strictly meaningful only for the center chain. For the corner chain, the two locations for the H atom covalently bound to O6 can also be found in schemes A and B. However, there is only a single location, found in scheme A, for the H atom covalently bound to O2 of the corner chain. Nevertheless, we can define (average percentage) occupancies for either scheme by averaging the fractional occupancies of all H atoms in that scheme and multiplying by 100; therefore, the fractional occupancies of the H atom covalently bound to O2 of the corner chain would be 1.0 and 0.0 in A and B, respectively. There is an important question that could not be answered in the previous crystallographic studies. Do these alternative
10.1021/bm800726v CCC: $40.75 2008 American Chemical Society Published on Web 10/15/2008
3134
Biomacromolecules, Vol. 9, No. 11, 2008
Nishiyama et al.
Figure 1. Thick dotted lines indicate the proposed cooperative networks of hydrogen bonds, with arrows indicating the donor-acceptor-donor directions for the A and B schemes. Thin dotted lines indicate the O3-H · · · O5 hydrogen bonds, and for the A network, the O6-H · · · O3 linkage.
schemes (i) interconvert dynamically and perhaps cooperatively, giving a temporally averaged final crystallographic structure, or (ii) are they static and present in different regions of the microfibril, thus producing a spatially averaged structure? Arguments can be made in support of both possibilities. In support of (i), if the enthalpies of schemes A and B are similar, then there might be a large entropic advantage in having dynamic disorder. The fractionally occupied H positions in A and B lie close together in an O2-H · · · H-O6 arrangement that is reminiscent of the flip-flopping dynamic H-bonding arrangement observed in neutron crystallographic studies of β-cyclodextrin.8 In the center sheet of Iβ (Figure 1), the A scheme corresponds to one of the two opposite directions for a unidirectional, infinite network of H bonds that runs continuously along the microfibril. The B scheme involves the same O6-H and O2-H hydroxyl hydrogen atoms but points them in the other direction.9 Therefore, this disorder might, at least partially, take the form of oscillations between these two networks, which is analogous to the β-cyclodextrin example. The key to this possibility is the gain in entropy at a particular temperature being greater than the increase in enthalpy because of weaker overall H bonding. In support of (ii), the occupancy of scheme A is 70-80% (the exact value depends on the method of structure refinement), which suggests that A could correspond to the H-bonding network for this allomorph, with B (30-20% occupancy) present in regions where the crystallographic environment is different, such as near surfaces or defects. A recent quantum mechanical modeling study using the Car-Parrinello molecular dynamics (CPMD) ab initio pseudopotential method has produced results in very good agreement with the crystal structure of Iβ with scheme A.10 However, that study also suggested that competition between H-bonding energy and electronic-strain energy might limit the size of crystalline domains. There is evidence from small-angle neutron scattering of cellulose fibers containing disordered regions or defects.11 There is also evidence from solid-state NMR spectroscopy of the conformations of chains at the interior and surfaces of microfibrils being different, with
only the interior chains being truly periodic and the surface chains having reduced intrachain H-bonding.12 A substantial fraction of the total number of chains is on the surfaces of typical cellulose crystals. For example, a crystallite with nine chains in each cross-sectional direction (about 74 Å on a side, similar to cotton) would have 32 of its 81 chains on the surface, and even a crystallite with 500 chains (typical of the large microfibrils isolated from tunicin that were used to determine the crystal structure of Iβ) would have 15% surface chains. Furthermore, recent molecular dynamics simulation calculations have provided results that indicate the possible presence of microheterogeneity in Iβ with complex H-bonding disorder.13 To gain further insight into the nature of H-bonding in Iβ, we have investigated the occupancies of schemes A and B at different temperatures using neutron diffraction and electronic structure theory (quantum mechanics) calculations. In the neutron diffraction study, data were collected at 295 K (room temperature) and then at 15 K so that any dynamical effects might be frozen. In particular, if flip-flopping occurs, then we might expect the oscillating strings of unidirectional H bonds to take only the more stable form at low temperature. In the quantum mechanics calculations, we compared the total electronic energies of schemes A and B to determine whether they are similar enough to allow an overall energy gain from flipflopping. We also looked at the stability of scheme A with conventional molecular dynamics calculations.
Experimental Section Sample Preparation. Cellulose Iβ nanocrystals (about 18 × 30 nm2 in cross section) were extracted and purified from tunicate using previously described techniques.14 These nanocrystals were then reassembled into highly oriented samples for neutron fiber diffraction studies using an approach that was designed to ensure cylindrical symmetry around the fiber axis. The final sample was cylindrical in shape with a diameter of ∼2 mm and a length of ∼10 mm. In previous neutron crystallographic studies, we reconstituted the nanocrystals into films by using a shear flow process and then assembled the parallel
Hydrogen Bonding in Cellulose Iβ films into 1-mm-thick stacks.15 This method provided us with highresolution and high-accuracy data, but the presence of the preferred orientation of nanocrystals within the sheets meant that the diffraction data did not possess true cylindrical symmetry. Consequently, for accurate data measurement, the samples had to be rotated in steps around the fiber axis, and diffraction images were collected at each step. In the new approach used here, we prepared PVA-cellulose samples by uniaxially stretching a viscoelastic gel containing the whiskerlike tunicin nanocrystals. An aqueous suspension of crystals was first mixed with polyvinyl alcohol (PVA), and then sodium borate was added to cross-link the PVA physically. This involved mixing a 10% aqueous solution of PVA (nominal average molecular weight of 124-186 kg/ mol) with the cellulose suspension to obtain a solution containing 3.5% PVA and 3.5% cellulose. Saturated Borax solution (500 mL) was then added to the PVA-cellulose mixed suspension (3 g) to make a viscoelastic gel. The gel was centrifugated in an Eppendorf polypropylene tube to remove air bubbles. A fiber was spun from this gel by the use of a pair of forceps and was dried with a weight hanging on one end. The final PVA-cellulose sample contained ∼50% PVA by weight. PVA remained amorphous in the final sample. The dried fibers were cut into lengths of ∼3 cm, and approximately 50 mg was put into a quartz capillary of 3 mm inner diameter and immersed in D2O with 0.1 N NaOH. The capillary with heavy water was then put in a portable reactor and heated to 210 °C for 1 h. The fibers swelled to occupy the whole space of the capillary during the hydrothermal treatment. After the removal of the capillary, the gel was again dried with both ends clamped on a manual stretching device. We checked the samples for crystal orientation by using X-ray diffraction. Samples with the best orientation were selected and glued into a bundle for neutron crystallographic measurements. Data Collection. Neutron diffraction data were collected on the D19 diffractometer at the Institut Laue Langevin in Grenoble, France. D19 is a four-circle diffractometer that is typically used for chemical and small biomolecule neutron crystallography as well as for high-angle fiber diffraction studies of biological and synthetic polymers.16-18 The D19 diffractometer has recently been completely rebuilt, improving the efficiency of data collection by a factor of approximately 25. The new instrument incorporates a large multiwire gas detector (30 × 120° angular aperture), a versatile monochromator and takeoff-angle assembly, and flexible arrangements for a variety of sample environment options. The new detector has cylindrical geometry, and outputs a data array of size 256 (vertically) × 640 (horizontally), giving an angular resolution of 0.19° in the equatorial direction. Generic strategies for collecting fiber diffraction data from cellulose fibers on D19 have been previously described.19 For the study reported here, the sample was attached with epoxy resin to an aluminum pin, which itself was then attached to the copper head of a Displex cooling system.20 A data set was collected at 295 K under ambient conditions by the use of a neutron beam of wavelength 1.9 Å. The sample was then enclosed in a cylindrical vanadium can to protect it from vacuum, which itself was enclosed in an aluminum can to hold vacuum. The sample was cooled to 15 K at a rate of 3 K/min, and a data set was collected while the temperature was maintained at 15 K. A further data set was collected at 15 K under identical conditions but with the sample removed so that the effects of scattering from the sample container system could be partially removed. The measuring time for each data set was 3 days. Data Processing. For each sample position, the recorded data were corrected for spatial distortion using the D19 program DCD19 (Turner et al., in preparation). We applied a correction for the variation in response of the detector over its area by dividing each frame by a frame of isotropic incoherent scattering collected from a vanadium rod. We also applied a correction for scattering from the sample container by subtracting the corresponding frame that was recorded without the sample in place. Because of the sample geometry and the quantity of incoherent scattering from H arising from the sample, the absorption correction
Biomacromolecules, Vol. 9, No. 11, 2008
3135
Table 1. Measured Values of the Reciprocal Radii (Å-1), Fobsd, of Well-Resolved and Unambiguously Indexed Reflections and Their Calculated Values, Fcalcd, after Least-Squares Minimizationa,b reflns h
k
1 -1 1 1 2 1 3 0 3 1 1 1 2 1 2 -2 1 3 0 0 0 1 2 1 0 1 1 -1 1 1 0 0 2 -2 root mean square
room temp
low temp
l
Fobsd
Fcalcd
Fobsd
Fcalcd
0 0 0 0 0 1 1 1 1 2 2 2 3 3 3 4 4 error
0.16912 0.19026 0.29900 0.39260 0.41827 0.21321 0.31408 0.34880 0.41610 0.19479 0.23058 0.35800 0.31453 0.33598 0.34534 0.38505 0.51070 0.001251
0.16800 0.18858 0.29955 0.38921 0.42137 0.21181 0.31469 0.34956 0.41557 0.19286 0.22860 0.35626 0.31425 0.33453 0.34533 0.38572 0.51154
0.16912 0.19177 0.30200 0.39864 0.42431 0.21457 0.31710 0.35220 0.41790 0.19479 0.23058 0.36019 0.31529 0.33658 0.34609 0.38505 0.51266 0.001224
0.16968 0.19021 0.30313 0.39502 0.42693 0.21324 0.31809 0.35279 0.41720 0.19279 0.22873 0.35924 0.31428 0.33529 0.34613 0.38557 0.51365
a The program LATREF, which was kindly provided by Dr. Mahendrasingam, was used. b The corresponding values of a, b, c, and γ are given in Table 2.
was particularly important, and it varied in a complex way with the scattering angle. Because the scattering background was dominated by isotropic incoherent scattering from the sample, the absorption correction was approximated to be inversely proportional to the scattering background at each position on the detector. The background was modeled as a bicubic spline function and was fitted using an adaptation of Bayesian analysis. A linear least-squares method with weighting was employed in scaling down the contribution of diffraction peaks. This weighting scheme was modified at each iteration to mimic the maximization of the sum of the log (likelihood) function.21 The slowly varying background thus obtained was used as a point-by-point scaling factor to compensate for sample absorption. The intensity was then remapped into cylindrical reciprocal space with a background level of unity. In this way, diffraction data collected at different goniometer settings could be remapped into a single contiguous reciprocal space image. The merged reciprocal space image was processed using the software suite FiberFix.22 The program LSQINT (Denny, unpublished) from that suite was used to refine a 2D fit of the background and Bragg intensities in reciprocal space. The peak positions of 17 well-resolved and unambiguously indexed reflections were measured, and a refinement was carried out that minimized the difference between the corresponding observed reciprocal-space radii and those calculated from the unit cell parameters using a least-squares procedure with a, b, c, and γ as variables (Table 1). The resulting unit cell parameters were held fixed as the background, and intensities were then fitted. The background was fitted using the Roving aperture method. A Lorentzian distribution was used to model the angular profile of the spots, and the width of this distribution was refined during maximum entropy refinement. A comparison of the fit of the data is shown in Figure 2. The relevant crystal data and a summary of the data collection parameters are given in Table 2. At the end of the intensity fitting, a standard deviation, σ, was calculated for each reflection directly from the root-mean-square difference between the calculated and observed pixel intensity values over the predicted profile of the spot. All reflections with similar reciprocal radial distances were regarded as accidentally overlapping because of cylindrical averaging, and the intensities associated with those reflections were merged. Structure Refinement. Structure refinement was carried out using previously described strategies for applying SHELX-9723 to fiber diffraction data.24 Atomic starting positions were taken from the highresolution room-temperature structure of Iβ.4 This structure was determined by the use of a combination of X-ray and neutron data, the
3136
Biomacromolecules, Vol. 9, No. 11, 2008
Figure 2. Neutron fiber diffraction patterns collected from cellulose Iβ at room temperature (top left quadrant) and low temperature (bottom left quadrant). The right-hand quadrants show the fits to Bragg intensities that were carried out using the software suite FiberFix.22 The background has been removed, and the integrated intensities that were determined from the original pattern have been displayed according to the determined crystallite size and tilt distribution. Table 2. Experimental Details room temp chemical formula space group a (Å) b (Å) c (Å) γ (deg) V (Å3)
low temp
C12H20O10 Crystal Data P21 7.76(1) 8.20(1) 10.37(1) 96.62(5) 655.5
P21 7.64(1) 8.18(1) 10.37(1) 96.54(5) 643.9
Data Collection independent reflns range of h range of k range of l diffractometer λ (Å)
151 0-4 -4 to 4 -7 to 7
146 0-4 -4 to 4 -7 to 7 D19 1.9
Refinement refinement on R[F2 > 2σ(F2)] ωR(F2) ∆Fmax ∆Fmin Frms
F2 0.244 0.533 1.15 -1.04 0.29
F2 0.241 0.567 1.36 -1.36 0.33
positions of C and O atoms were determined by X-ray structure refinement, the positions of non-D atoms were identified in the X-ray analysis and were then fixed, and the positions of D atoms were determined by neutron structure refinement. Similarly, in the refinement of the 295 K neutron structure in this study, the positions of non-D atoms were fixed. With no D atoms included, the scale factor, K, and global thermal parameter, B, were refined against 150 reflections with d < 1.3 Å, which resulted in values of 0.393 and 0.759 for R and Rω. D atoms were then explicitly added to O3, O2, and O6, and their positions on O2 and O6 were also refined. This involved adding four additional degrees of freedom to the refinement. Refinement from starting positions for D atoms corresponding to scheme A for both center and corner chains resulted in values of 0.289, and 0.649 for Rω, respectively. Refinement from starting positions for D atoms corre-
Nishiyama et al. sponding to scheme B for both center and corner chains resulted in values of 0.304 and 0.665 for R and Rω, respectively, with the D on O2 of the corner chain refining to the same position as that in the refinement of scheme A (i.e., its fractional occupancy is 0.0 for scheme B). The refinement of both scheme A and scheme B positions for D atoms of both corner and center chains with fractional disorder resulted in values of 0.279 and 0.625 for R and Rω, respectively, and involved four additional parameters. The D atom on O6 of the center chain refined to unit and zero fractional occupancies for schemes A and B, respectively. The improvement in the agreement of scheme A compared with scheme B is significant at a confidence level of >90%, and the improvement in the agreement of fractionally disordered scheme A/B compared with scheme A is significant at a confidence level of >90%.25 The average occupancies of schemes A and B were 67 and 33%, respectively. A final refinement with 17 reflections removed to establish a flat analysis of variance resulted in values of 0.244 and 0.533 for R and Rω, respectively, and average occupancies of 76 and 24% for A and B, respectively. Refinement of the 15 K neutron structure was carried out in a similar manner. Again, the non-D atoms were fixed at their positions in the room-temperature structure. The reason for using the room-temperature structure was because we did not have an X-ray structure for Iβ at 15 K. The justification for using a room-temperature structure is a previous X-ray crystallographic study that reported structures for Iβ (without H atoms on the hydroxyl groups) at room temperature and 100 K; the molecular conformation was essentially unchanged by cooling to 100 K, and we might therefore expect the conformation to be similar at 15 K.26 Refinement from starting positions for D atoms corresponding to scheme A for both center and corner chains resulted in values of 0.306 and 0.700 for R and Rω, respectively. Refinement from starting positions for D atoms corresponding to scheme B for both center and corner chains resulted in values of 0.317 and 0.709 for R and Rω, respectively. Finally, refinement of both scheme A and B positions for D atoms of both center and corner chains with fractional disorder resulted in values of 0.287 and 0.676 for R and Rω, respectively, and involved four additional parameters. The D atoms on both O6 of the center chain and O2 of the origin chain refined to unit and zero fractional occupancies for schemes A and B, respectively. The improvement in the agreement of fractionally disordered network A/B compared with scheme A is significant at a confidence level of >90%.25 The average occupancies of schemes A and B were 79 and 21%, respectively. A final refinement with 18 reflections removed to establish a flat analysis of variance resulted in values of 0.241 and 0.567 for R and Rω, respectively, and average occupancies of 83 and 17% for A and B, respectively. Quantum Mechanics Calculations. Two-dimensional sheets of cellulose chains with H-bonding schemes A and B were based on the cellobiose repeat of the center sheet of Iβ. These sheets, extending 100 bohrs (∼53 Å) in the unit cell b-axis and c-axis directions, were used as starting structures for energy minimization, using the facility within Gaussian 0327 that allows for calculations on periodic structures. Calculations were carried out with HF theory using the 6-31G(d) basis set and with PBEPBE theory28 using the 6-31+G(d,p) basis set. HF/ 6-31G(d) is a commonly used, fairly economical level of theory. HF/ 6-31G(d,p) theory, which gives relative energies that are similar to those of the HF/6-31G(d) level used here, has reasonably small basis set superposition error (BSSE) (0.96 kcal/mol for the water dimer).29 PBEPBE, a pure density functional, was selected for anticipated advantages in computational speed. The use of a diffuse function (indicated by the + sign) in the 6-31+G(d,p) basis set should keep BSSE to reasonable levels when used with different density functionals. Besides fixed geometry calculations, we enforced periodicity during energy minimization with and without constraints that maintained the b- and c-unit cell dimensions. Results were similar for both of the minimizations, so only the values from the minimization without the unit-
Hydrogen Bonding in Cellulose Iβ
Biomacromolecules, Vol. 9, No. 11, 2008
3137
Table 3. Minimization of 2D Periodic Cellulose Iβ with Quantum Mechanicsa exptb relative E (optimized kcal/mol) relative E (initial kcal/mol) b axis (Å) c axis (Å) R φ O5-C1-O4C4 Ψ C1-O4-C4-C5
8.20 10.38 90.00 -89/-98 -147/-142
HF A schemeb
HF B schemeb PBEPBE A scheme PBEPBE B scheme
0 0 8.27 10.41 89.95 -94 -147
7.4 18.4 8.53 10.29 88.26 -92, -88 -152, -147
0 0 8.23 10.46 89.89 -95 -144
10.4 18.1 8.45 10.32 90.01 -88 -151
Hydrogen Bondingc O3-H · · · O5
1.75, 2.70, 162/1.97, 2.76, 137 1.94, 2.83, 152 1.90, 2.76, 149 1.84, 2.79, 160
O2-H · · · O6 intramolecular O6-H · · · O3 intermolecular O6-H · · · O2 intermolecular
1.90, 1.83, 1.85, 2.15, 2.72, 3.63,
1.76, 2.71, 159
A Scheme 2.86, 2.76, 2.67, 2.80, 3.29, 3.62,
165 159 140 122 117 81d
1.88, 2.84, 179
1.75, 2.75, 174
1.92, 2.84, 159
1.77, 2.74, 164
2.84, 3.56, 132
2.82, 3.46, 123
B Scheme O6-H · · · O2 intramolecular O2-H · · · O6 intermolecular
1.97 2.86, 152 1.88, 2.76, 150 2.37, 3.29, 156 na, 3.62, nae
2.28, 2.35, 2.10, 2.62,
3.22, 3.29, 3.03, 3.52,
171 168 163 157
2.50, 3.45, 164 1.90, 2.87, 166
a Models were based on center chains. b When two sets of values are listed for the experimental values, they are for the center and corner chains, respectively. Two sets of values for the models occur when there are substantial deviations from two-fold screw-axis symmetry. c Listed for each bond are the H · · · O distance (Å), the O · · · O distance (Å) and the O-H · · · O angle (deg). Minor hydrogen bonds that involve the glycosidic oxygen atoms are not reported. d These parameters would not indicate that a hydrogen bond had formed. e na ) not applicable because the corner chain proton on O2 did not have a B-network orientation.
Figure 3. Proposed A and B schemes before (black) and after (cyan) energy minimization under 2D periodic conditions at the PBEPBE/631+G(d,p) level of theory.
cell dimensional constraints are reported. The results are summarized in Table 3 and in Figure 3. Empirical Force-Field Molecular Dynamics. A minicrystal (or cluster) consisting of 19 cellooctaose chains arranged according to the Iβ crystal structure was constructed. It had a central sheet consisting of five chains, with sheets of four corner chains above and below the central sheet (when viewed down the chain-axis direction). Sheets of three central chains were then positioned above and below the fourchain sheets. Simulations were carried out with the AMBER 830 program and GLYCAM 0431 force field for 10 ns. There were no periodic boundary conditions in the simulation, and no equilibration was done. The simulation produced the canonical (constant temperature) ensemble using Langevin dynamics temperature regulation with a collision frequency of 1.0 and a temperature of 293 K. There was no nonbonded cutoff so all interactions were considered. Results are summarized in Table 4 and in Figure 4.
Table 4. Molecular Dynamics Simulations of 19-Chain Cellulose Iβ Model, Hydrogen Bonding Scheme A last frame of MDa
expt
Unit-Cell Dimensions a axis b axis c axis γ
7.51 ( 0.05 Å 8.11 ( 0.07 Å 10.58 ( 0.03 Å 97.7 ( 1.3°
7.78 Å 8.20 Å 10.38 Å 96.5°
Model-Chain Parameters φ corner chain O5-C1-O4-C4 φ center chain O5-C1-O4-C4 Ψ corner chain C1-O4-C4-C5 Ψ corner chain C1-O4-C4-C5 τ corner chain C1-O4-C4 τ center chain C1-O4-C4
-95.9 ( 1.1° -93.3 ( 0.6° -142.8 ( 1.1° -141.9 ( 1.4° 117.3 ( 0.4° 116.8 ( 0.4°
-98.5° -88.7° -142.3° -147.1° 115.1° 116.2°
a Values shown are averages of four measurements of the most central cell with the standard deviation.
3138
Biomacromolecules, Vol. 9, No. 11, 2008
Figure 4. Minicrystals before and after molecular dynamics simulation. The upper image shows the 19-chain model, along with the unit cell dimensions a, b, and γ before simulation. The lower image shows the model after simulation with first and eighth glucose residues on each chain removed. Hydrogen atoms are not shown, and Roman numerals refer to surface chains that were examined. (See the text and Figure 5.)
Results The neutron crystallographic results indicate that there are no large differences in the occupancies of H-bonding schemes A and B at 295 K (76 and 24%, respectively) compared with those at 15 K (83 and 17%, respectively) and that those occupancies agree with the results from the previously reported room-temperature structure (70-80%).4 At both temperatures, the H-bonding arrangement predominantly corresponds to scheme A. At both temperatures, the refined positions of fractionally occupied D atoms in scheme A did not change greatly from their starting positions. However, the same was not true for fractionally occupied D atoms in scheme B, and the refined positions for those atoms do not strictly correspond to an H-bond network. The fact that there appears to be no change in the occupancies of A and B as Iβ is slowly cooled to 15 K suggests that the H-bonding disorder is not dynamic but rather that these structures represent a statistical spatial average of a crystallographic structure in which the H bonding corresponds to scheme A but also contains defects or heterogeneity in H bonding. An argument could be made against this conclusion if the energies of schemes A and B were very similar so that one was not significantly preferred over the other at low temperature. However, the energy differences between A and B before (18.1 kcal/mol of cellobiose units) and after (10.4 kcal of cellobiose units at the PBEPBE level of theory, Table 3) minimization indicate that scheme A is clearly preferred and that the fraction of atoms with B-network positions would be insignificant. The HF calculations also clearly support the A scheme. During the quantum mechanical calculations, the structural changes resulting from energy minimization were small for A but more significant for B (Figure 3). In particular, there is good retention of the b-axis length after minimization for A but less so for B (Table 3). The calculations were carried out with translational symmetry maintained, but there was no imposition of two-fold screw symmetry. In the case of the HF calculations for the B scheme,
Nishiyama et al.
the model substantially deviated from screw-axis symmetry, especially the O2-H · · · O6 values that indicate a 0.5 Å discrepancy for adjacent residues after minimization. Most other details do not indicate substantial departure from symmetry. Generally, the PBEPBE calculations gave hydrogen bonds that were 0.1 Å shorter than the HF calculations. H · · · O distance values for both levels of theory were usually between the experimental values for the center and corner chains for the A scheme, whereas the angle values in the model were often somewhat closer to 180°. This is not surprising because the differences between the corner and center chains arise from the 3D packing, and that is not present in the model. Whereas both A and B schemes are stabilized by the O3-H · · · O5 H bond, the A scheme has two other high-quality H bonds, whereas the B scheme has no other strong bonds (H · · · O < 2.0 Å), which explains, at least partially, the preference for A. The crystallographic results show that the unit cell contracts significantly (by ∼1.6%) in only the direction of the a axis when Iβ is cooled from 295 K to 15 K. This contraction corresponds to a value for the thermal expansion coefficient in this direction, Ra, of 6(1) × 10-5 °C-1. In previous studies in which Iβ was cooled from 295 K to 100 K, a value of 4.6 × 10-5 °C-1 was determined for Ra.26 In a separate study in which Iβ was heated from 293 K to 473 K, a value of 4.3 × 10-5 °C-1 was calculated for Ra at 293 K.32 Furthermore, a value of (7 to 8)× 10-5 °C-1 for Ra over the temperature range of 122 K to ambient has been reported for cotton, which predominantly consists of Iβ.33 The similar values of Ra reported for Iβ over a wide range of temperatures below ambient have been previously explained by the possible presence of simple harmonic molecular librations by analogy to the anisotropic expansion of glycine.26,34 The value of Ra reported in this study seems to indicate that the effects of molecular librations determine the anisotropic contraction of Iβ from ambient to 15 K and that there is little variation in this behavior over this temperature range. In the crystallographic studies, we made an approximation that the positions of O and C atoms and therefore the conformations of the cellulose chains are essentially the same at 295 K and 15 K. Even if this approximation is not exact, it is unlikely to affect our conclusions about the refined positions and occupancies of the H/D atoms. However, it will affect the accuracy of the intersheet atomic distances in the structures at 15 K and 295 K, and we have therefore not examined these distances in detail. However, it is clear that there is no strong intersheet H bonding but rather a variety of potential weak C-H · · · O interactions. The intrasheet H bonds are of a strong O-H · · · O type that is characterized by a narrow distribution of relatively short bond distances with strong directionality, whereas the intersheet H bonds are of a weak C-H · · · O type that is characterized by a broad range of relatively long distances that are only weakly directional. In the CPMD calculations of Qian et al., the intersheet interaction energy is 8 times weaker than the intrasheet interaction energy when van der Waals interactions are included.10 The electrostatic component of the intersheet bond energy might be expected to fall off as r-1 compared with r-6 for the exchange and dispersion components. The intersheet interactions will therefore be softer, and this softness may allow the librational behavior, perhaps in the form of oscillations around the glycosidic linkage, to dominate the response of Iβ to variations in temperature. In the crystal structure of Iβ, the chains are strained, possibly to facilitate better H bonding.4 The CPMD calculations of Qian et al. show that this strain results in an increase in electronic energy.10 For small assemblies of chains, the increase in
Hydrogen Bonding in Cellulose Iβ
Biomacromolecules, Vol. 9, No. 11, 2008
3139
Figure 5. Surface chain I from the molecular dynamics study. (See Figure 4.) The top of the model is the crystallite surface, and the lower part of the model is adjacent to the crystallite interior. Hydrogen bonds in yellow correspond to bonds that are characteristic of scheme A, and those in magenta on the top side correspond to scheme B. At the left end of the model, there was substantial deviation from two-fold symmetry, and the O6 atom took a gg orientation. Dotted blue lines correspond to arrangements that are not part of either A or B.
intramolecular electronic energy is partly compensated for by the lower energy from intermolecular H bonding. However, there is a finite crystal size beyond which the balance between H bonding and electronic energy is no longer maintained, and this is predicted to correspond to four stacked sheets. If the increase in electronic energy is released by a breakdown in H bonding in scheme A in every fifth sheet, then this would result in the experimentally observed occupancies for A and B that are reported here. In the quantum calculations reported here, H-bonding scheme B is energetically unfavorable. During energy minimization of the crystallographic coordinates corresponding to scheme B, the Iβ structure was significantly distorted (Figure 3). In the molecular dynamics studies of Mazeau,13 minicrystals of Iβ with H-bonding scheme A were stable and had the lowest energy in the PCFF35 force field. Our own studies with AMBER/ GLYCAM gave results similar to those of Yui et al.36 despite our longer simulation, our smaller crystal size, and our absence of water. Although our minicrystal twisted during simulation, as reported by Yui et al., its central part retained a reasonably close approximation to the original A-scheme structure over the 10 ns simulation. Orientations of the O6-H and O2-H atoms on the central cellobiose unit stayed close to their original values throughout the simulation; none of the saved frames corresponded to the B scheme. The stability of the A system confirms Mazeau’s finding, but differs from the conclusions of Matthews et al. 37 The latter simulation, with the CHARMM program and CSFF force field, led to a substantial alteration of the original structure, highlighting a dependence on the exact model and force field. Mazeau’s models with H-bonding scheme B were unstable. However, many stable intermediate structures were found with a variety of different H-bonding arrangements involving O2, O6, and O3. Taken together, these results suggest, as first suggested in Mazeau’s paper,13 that network B is a crystallographic average of a number of disordered H-bonding arrangements and does not correspond in reality to an actual extended H-bonding network. Unlike the periodic modeling used herein with quantum mechanics and in Mazeau’s work, the finite minicrystal model in the present molecular dynamics studies also models the crystal surfaces. The lower after-simulation model in Figure 4 has six chains labeled with Roman numerals. Figure 5 shows the reasonably typical chain I. In the six cellooctaose surface chains labeled in Figure 4, there were eight instances of a cooperative network formed by O6-H · · · O2-H · · · O3-H · · · O5 (Figure 5). This is the same network that was found in scheme B except that it lacks the intermolecular H bond that would be formed by O2-H · · · O6. These three-bond networks were all on the exterior side of the chains, whereas the interior side mostly retained the intramo-
lecular components of scheme A. Missing scheme A hydrogen bonds to O2 atoms on the interior side had reasonable angular alignment but exceeded the generous H · · · O distance criterion of 2.8 Å. None of the exterior sides of the six chains retained a continuous scheme A network, although some instances of scheme A hydroxyl orientations remained after the 10 ns of dynamics. Both schemes have O3-H · · · O5 and O6 · · · H-O2 H bonds, so the basis for these observations must lie in the O2-H and O3-H hydroxyl groups. In B, they tend to form a weak hydrogen bond (so weak that it was not indicated in Figure 1), whereas in A, they are oriented in opposite directions and suffer destabilizing dipole-dipole interactions. Even when the exterior O6-H groups did not participate in the three-bond networks, the O2-H groups often reoriented to make the weak H bond to O3 in these vacuum simulations. In other instances, the surface O3-H groups turned away from O5 to make the weak link to O2, avoiding the dipole-dipole interaction one way or another. It is unlikely that crystal surface effects can explain all of the observed H-bond disorder. For the large tunicin nanocrystals used in this study, only ∼7% of the hydroxyl groups will be at the surface. In addition, we might expect some of those groups to have exchanged back from OD to OH because of surface exposure, which rendered them less visible in the neutron scattering density maps. The observed H-bond disorder must be, at least in part, interior to the microfibrils. The surface effects observed in our microcrystal simulations indicate that short networks found in scheme B may be possible under certain conditions interior to the microfibril if the extended interchain H bonding observed in scheme A is disturbed.
Conclusions The results presented here from neutron crystallography, empirical force field molecular dynamics, and ab initio quantum mechanics calculations suggest that at ambient temperatures and below, cellulose Iβ nanocrystals from hydrolyzed tunicin have most of their chains in a crystalline Iβ arrangement with H bonding as in scheme A. Smaller regions of static H-bond disorder exist, perhaps at defects within crystalline domains, interfaces between crystalline domains, or at the surfaces of the microfibril. The exact quantitative association of defects, interfaces, and surfaces with static H-bond disorder remains to be determined. Although H bonding in these disordered regions may be complex and may not strictly correspond to an extended H-bonding network, our results suggest that it could have certain features in common with scheme B. In particular, when the O2-H · · · O6 interchain H bond in scheme A is disrupted, the O2-H and O3-H hydroxyl groups can suffer destabilizing dipole-dipole interactions that favor the reorientation of O2-H
3140
Biomacromolecules, Vol. 9, No. 11, 2008
so that it forms a three 3-bond cooperative network that is a feature of scheme B. Although the experimental results presented here were obtained with large nanocrystals of cellulose from hydrolyzed tunicin, they may be relevant to our understanding of some of the different properties of cellulose microfibrils found in other naturally occurring systems. Depending on their origin, these microfibrils can contain smaller regions of the crystalline Iβ form in addition to regions of the crystalline IR form and also lessordered cellulose. If the static H-bonding disorder is present at the boundaries between crystalline domains of scheme A in the crystalline Iβ regions, then some of the properties of those microfibrils under strain/stress may be influenced by the interaction of crystal dislocations with those boundaries, although the nature and effects of crystal dislocations in polymeric systems are not well known. In addition, although the crystalline Iβ regions in microfibrils will be dominated by H bonding in scheme A, the presence of H-bonding features found in scheme B may be more relevant to the surface and interface properties of those crystalline Iβ regions, in particular the chemical reactivity and susceptibility to enzymatic hydrolysis. However, other types of more extensive disorder present in naturally occurring microfibrils may be more significant. Finally, although we have detected only static H-bond disorder in cellulose Iβ at ambient temperatures and below, it is likely that dynamical disorder will be present at higher temperatures. Therefore, there may be two temperature regimes for cellulose in which either static or dynamic disorder is important as well as a range of temperatures that corresponds to a transition between static and dynamic disorder.
References and Notes (1) Atalla, R. H.; VanderHart, D. L. Science 1984, 223, 283. (2) Horri, F.; Yamamoto, H.; Kitamaru, R.; Tanahashi, M.; Higuchi, T. Macromolecules 1987, 20, 2946. (3) Sugiyama, J.; Vuong, R.; Chanzy, H. Macromolecules 1991, 24, 4168. (4) Nishiyama, Y.; Langan, P.; Chanzy, H. J. Am. Chem. Soc. 2002, 124, 9074. (5) Nishiyama, Y.; Sugiyama, J.; Chanzy, H.; Langan, P. J. Am. Chem. Soc. 2003, 125, 14300. (6) The parallel-up arrangement used here is the one defined by French, A. D.; Howley, P. S. In Cellulose and Wood-Chemistry and Technology, Schuerch, C., Ed.; Wiley: New York, 1989; p 159. Attention was originally drawn to the distinction between parallel-up and paralleldown in Gardner, K. H.; Blackwell J. Biopolymers 1974, 13, 19752001. (7) The glycosidic torsion angles, Φ and Ψ, which describe the relative orientation of adjacent glycosyl residues in the same chain, are defined by (O5-C1-O1-C4) and (C1-O1-C4-C5), respectively. The glycosidic bond angle, τ, is defined by (C1-O4-C4). The conformation of the hydroxymethyl group is defined by two letters, the first referring to the torsion angle χ (O5-C5-C6-O6) and the second referring to the torsion angle χ′ (C4-C5-C6-O6). An ideal tg conformation would be defined as the set of two angles, 180 and60°. (8) Saenger, W.; Betzel, Ch.; Hingerty, B.; Brown, G. M. Nature 1982, 296, 581. (9) Jarvis, M. Nature 2003, 426, 611. (10) Qian, X.; Ding, S.-Y.; Nimlos, M. R.; Johnson, D. K.; Himmel, M. E. Macromolecules 2005, 38, 10580. (11) Nishiyama, Y.; Kim, U. J.; Kim, D. Y.; Katsumata, K. S.; May, R. M.; Langan, P. Biomacromolecules 2003, 4, 1013. (12) Vietor, R. J.; Newman, R. H.; Ha, M. A.; Apperley, D. C.; Jarvis, M. C. Plant J. 2002, 30, 721.
Nishiyama et al. (13) Mazeau, K. Cellulose 2005, 12, 339. (14) Nishiyama, Y.; Langan, P.; Chanzy, H. Fibre Diffr. ReV. 2003, 11, 75. (15) Nishiyama, Y.; Kuga, S.; Wada, M.; Okano, T. Macromolecules 1997, 30, 6395. (16) Sapede, D.; Seydel, T.; Forsyth, V. T.; Koza, M. M.; Schweins, R.; Vollrath, F.; Riekel, C. Macromolecules 2005, 38, 8447–8453. (17) Gardner, K. H.; English, A.; Forsyth, V. T. Macromolecules 2004, 37, 9654–9656. (18) Mason, S. A.; Clergeau, J. F.; Geurard, B.; Archer, J.; Fuard, S.; Allibon, J.; Howard, J. A. K.; Davidson, M.; Forsyth, V. T. J. Appl. Crystallogr. In preparation. (19) Langan, P.; Denny, R. C.; Mahendrasingam, A.; Mason, S. A.; Jaber, A. J. Appl. Crystallogr. 1996, 29, 383. (20) Archer, J.; Lehmann, M. S. J. Appl. Crystallogr. 1986, 19, 456–458. (21) David, W. I. F.; Sivia, D. S. J. Appl. Crystallogr. 2001, 34, 318–324. (22) Rajkumar, G.; AL-Khayat, H. A.; Eakins, F.; Knupp, C.; Squire, J. M. Fibre Diffr. ReV. 2006, 14, 10–14. (23) Sheldrick, G. M. SHELX-97: A Program for the Refinement of SingleCrystal Diffraction Data; University of Go¨ttingen: Go¨ttingen, Germany, 1997. (24) Langan, P.; Nishiyama, Y.; Chanzy, H. Biomacromolecules 2001, 2, 279. (25) Hamilton, W. C. Acta Crystallogr. 1965, 18, 502. (26) Langan, P.; Sukumar, N.; Nishiyama, Y.; Chanzy, H. Cellulose 2005, 12, 551. (27) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (28) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865–3868. (29) Lii, J.-H.; Ma, B.; Allinger, N. L. J. Comput. Chem. 1999, 20, 1593– 1603. (30) Case, D. A.; Darden, T. A.; Cheatham, T. E.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Merz, K. M.; Wang, B.; Pearlman, D. A.; Crowley, M.; Brozell, S.; Tsui, V.; Gohlke, H.; Mongan, J.; Hornak, V.; Cui, G.; Beroza, P.; Schafmeister, C.; Caldwell, J. W.; Ross, W. S.; Kollman, P. A. AMBER 8; University of California: San Francisco, 2004. (31) Woods, R. J.; Dwek, R. A.; Edge, C. J.; Fraser-Reid, B. J. Phys. Chem. 1995, 99, 3832. (32) Wada, M. J. Polym. Sci., Part B: Polym. Phys. 2002, 40, 1095. (33) Seitsonen, S.; Mikkonen, I. Polym. J. 1973, 5, 263. (34) Langan, P.; Mason, S. A.; Myles, D.; Schoenborn, B. P. Acta Crystallogr. 2002, B58, 728. (35) Maple, J. R.; Hwang, M. J.; Stockfisch, T. P.; Dinur, U.; Waldman, M.; Ewing, C. S.; Hagler, A. T. J. Comput. Chem. 1994, 15, 162– 182. (36) Yui, T.; Nishimura, S.; Akiba, S.; Hayashi, S. Carbohydr. Res. 2006, 341, 2521–2530. (37) Matthews, J. F.; Skopec, C. E.; Mason, P. E.; Zuccato, P.; Torget, R. W.; Sugiyama, J.; Himmel, M. E.; Brady, J. W. Carbohydr. Res. 2006, 341, 138.
BM800726V
