Theoretical Study of the Structural Stability of Molecular Chain Sheet

Jul 22, 2014 - Kevin Conley , Louis Godbout , M.A. (Tony) Whitehead , Theo G.M. van de Ven. Carbohydrate Polymers 2016 135, 285-299 ...
0 downloads 0 Views 3MB Size
Download PDF
Article pubs.acs.org/JPCB

Theoretical Study of the Structural Stability of Molecular Chain Sheet Models of Cellulose Crystal Allomorphs Takuya Uto, Sho Mawatari, and Toshifumi Yui* Department of Applied Chemistry, Faculty of Engineering, University of Miyazaki, Nishi 1-1 Gakuen-kibanadai, Miyazaki 889-2191, Japan S Supporting Information *

ABSTRACT: The structural stabilities of the molecular chain sheets constituting the crystal structures of the cellulose allomorphs Iα, Iβ, II, and IIII were investigated by density functional theory (DFT) optimization of the isolated chain sheet models with finite dimensions. The DFToptimized chain sheet models of the two native cellulose crystals developed a right-handed twist with a similar amount of twisting. The DFT-optimized cellulose II (010) and (020) models twisted in opposite directions with right- and left-handed chirality, respectively. The cellulose IIII (1−10) model retained the initial flat structure after the DFT-optimization. The structural features of the DFT-optimized chain sheet models were reflected in the structures of the parent crystal models observed in solvated molecular dynamics (MD) simulations. The minor conformations of the hydroxymethyl groups proposed in the real crystal structures were detected in the MD crystal models and the DFT-optimized (010) model of cellulose II. The crystal chain packing and crystal conversions are interpreted in terms of principal chain sheet stacking.

■

INTRODUCTION

In the native cellulose crystal structures, molecular chains form flat chain sheet structures connected by O3···O6 intermolecular hydrogen bonds.13−15 The polar functional groups are hydrogen bonded, causing hydrophobic interactions to dominate between the chain sheets. The main difference between cellulose Iα and Iβ is the relative displacement of the chain sheets along the chain axis. The structure of cellulose IIII, however, exhibits a two-dimensional O2···O6 intermolecular hydrogen-bonded network with respect to the ab projection, giving rise to two types of corrugated chain sheets along either the (100) or (1−10) lattice planes.10 Conversely, two parallel and one antiparallel chain sheets can be identified in cellulose II, with the parallel chain sheets along the (010) and (020) lattice planes and the antiparallel chain sheet along the (110) lattice plane.11−13 Microscopic observations of a native cellulose microfibril using atomic force and transmission electron microscopy often reveal the formation of a right-handed twist.19,20 For example, narrow twisted ribbons of cellulose were observed during biosynthesis of bacterial cellulose,21,22 and the twisted cellulose microfibrils in spruce wood may restrict the distance along the fiber axis with their coalescence at the crystal faces.23 We have reported molecular dynamics (MD) simulations of finite crystals of the different allomorphs in aqueous solution using the AMBER force field coupled with Glycam04 parameters, and the right-handed fiber twists were reproduced in the solvated cellulose Iα and Iβ crystal models with the Iα crystal models showing greater twisting.24,25 Just prior to our reports,

Cellulose, the most abundant linear poly(1 → 4)-β-D-glucan biopolymer, has many current applications as well as potential new uses, such as nanomaterials and biofuel. Cellulose is produced naturally in highly crystalline microfibrils. It is wellknown that cellulose forms various crystalline allomorphs derived from the two native crystalline phases: cellulose Iα and Iβ.1−3 One important physical aspect of these allomorphs is that the Iα phase is metastable and can be converted into the Iβ phase by hydrothermal treatment4 or high-temperature treatment in organic solvents or even in helium gas.5 Cellulose IIII can be obtained by treatment of the cellulose Iβ form with liquid ammonia or amines and is readily converted back into cellulose Iβ by hot-water treatment.6 In the crystal structures of the native and IIII allomorphs, all of the molecular chains are aligned with the same polarity, which is known as parallel chain packing.7 In contrast, the unit cell of cellulose II has one chain oriented in the opposite direction to the other,8,9 although it has a similar ab projection structure to that of cellulose IIII.10 The native crystal forms are irreversibly converted into cellulose II by either mercerization or regeneration processes, which involve a complete change in the cellulose chain polarities from parallel to antiparallel.11−13 Another difference in the structures of the cellulose allomorphs is the orientation of the hydroxymethyl groups: the native forms adopt a rare trans−gauche (tg) conformation,14,15 whereas the cellulose IIII10 and II11−13 crystals, along with the related β-D-glucose, cellobiose, and cellotetraose crystals,16,17 all adopt a gauche− trans (gt) conformation. Some of the small carbohydrate crystals involve diversity in a hydroxymethyl group conformation, including the gauche−gauche (gg) conformation.18 © 2014 American Chemical Society

Received: April 10, 2014 Revised: July 21, 2014 Published: July 22, 2014 9313

dx.doi.org/10.1021/jp503535d | J. Phys. Chem. B 2014, 118, 9313−9321

The Journal of Physical Chemistry B

Article

(110) chain sheet, antiparallel chains (see Figure 1). The set of atomic coordinates proposed by X-ray and neutron diffraction

Matthews et al. also reported the right-handed twist of the cellulose Iβ crystal in a similar MD study using the CHARMM/ CSFF force field.26 The fiber twist of the native cellulose crystals has been further examined using various force fields and carbohydrate parameters, including the atomistic force fields AMBER with Glycam06,27,28 CHARMM C35,28 GROMOS 45a4,28 and OPLS29 and the coarse-grained force fields M3B,30 REACH,31 and MARTINI.32,33 In MD simulations with the coarse-grained MARTINI force field, the crystal models showed either a right-32 or left-handed33 twist depending on the number of constituent cellulose chains. Matthews et al. reported that the twisted crystal model untwists after long simulations, along with other internal changes such as rotation of the hydroxylmethyl groups. These changes occurred at both high and low temperatures with the CHARMM C35 and AMBER/Glycam06 force fields.26 Hadden et al. suggested that microfibril twisting was favored by van der Waals interactions and counteracted by both the intramolecular hydrogen bonds and solvent effects at the microfibril surface from MD simulations using the Glycam06 parameters.34 In addition to the native crystal forms, we have performed similar crystal model studies of the cellulose IIII form. It was found that the solvated cellulose IIII crystal models retained their initial structure during the MD simulations.25 However, at a temperature of 370 K, irreversible rotations occurred in the hydroxymethyl group and the conformation converted to tg.35,36 This was accompanied by a change in the intermolecular hydrogen bonding scheme of the O3···O6 bond in the (1−10) chain sheet, which was converted into a cellulose I-like flat chain sheet along with dissipation of the cellulose IIII (100) chain sheet,36 which is consistent with the crystalline conversion scheme from cellulose IIII to Iβ proposed by Wada based on in-situ X-ray diffraction measurements.37 As the structure conversion proceeded, the chain sheet developed a right-handed twist.36 This scheme of chain sheet conversion was further examined by density functional theory (DFT) calculations of the isolated chain sheet models, where the cellulose Iβ (100) and cellulose IIII (100) and (1−10) models with finite dimensions were optimized.38 The cellulose Iβ (100) model formed a slight right-handed twist and the cellulose IIII (100) model significantly deviated from the flat sheet structure, whereas the (1−10) model retained the initial sheet structure. Here, we report systematic DFT optimization calculations of the isolated molecular chain sheets composing the cellulose Iα and Iβ, II, and IIII allomorphs combined with MD simulations of their crystal models in aqueous solution. In contrast to our previous studies, the present study focuses on discussing the structural changes of the DFT-optimized chain sheet models using more accurate DFT calculations. Although it is unrealistic for a chain sheet to be isolated, the DFT optimization was performed to examine sheet deformation not to reproduce a sheet structure. The effect of chain sheet deviation on the parent crystal structure is assessed by MD simulations of the corresponding solvated crystal model.

Figure 1. Projections of the ab base planes of the seven crystal models with the unit cell frames and lattice plane labels. The “u” and “d” notations in the cellulose II B36-1 and -2 models denote the “parallel up” and “parallel down” chain sheets composing the (010) and (020) lattice planes, respectively.

were used for model construction. Each chain sheet model consisted of four cello-oligomers with degrees of polymerization (DP) = 6 or 8. The positions of the hydrogen atoms were initially optimized with the AMBER/Glycam 0639 parameters followed by full DFT optimization of all of the atom positions. In the DFT optimization, the present study used either the B3LYP or the CAM-B3LYP40 functional coupled with the 6-31G(d) basis set. The latter functional introduced long-range corrections by the CAM functional and was only used for the 4× cello-octamer chain sheet models. The single-point energy of the optimized model was then calculated using the 6-311G(d) basis set. The binding energy (ΔEbind) between adjacent oligomers was calculated as the difference between the DFT energies of the total model (Etotal) and its separate parts (E1 and E2), as defined by the following equation: 1 ΔE bind = [Etotal − (E1 + E2)] 2 For the four-chain sheet models, E1 and E2 represent the energies of an interior oligomer and the remainder of the sheet, respectively. Since the four-chain sheet models contained two interior oligomers, two values of ΔEbind were obtained from the respective interior oligomers, and they were averaged to give E1. The stabilization energy of the chain sheet (ΔEsheet) was defined as the sum of E1 and ΔEbind:

■

EXPERIMENTAL METHODS The six chain sheet models were constructed based on the corresponding crystal structure data, which were the (110) chain sheet of cellulose Iα,14 the (100) chain sheet of cellulose Iβ,15 the (1−10) chain sheet of cellulose IIII,10 and the (010), (020) and (110) chain sheets of cellulose II.12 In the three chain sheets of cellulose II, the (010) and (020) chain sheets consist of parallel “up” and “down” chains, respectively, and the

ΔEsheet = E1 + ΔE bind

The basis set superposition error was estimated using the counterpoise method.41,42 All calculations were carried out using the Gaussian 09 program package.43 Figure 1 shows the ab base planes of the cellulose crystal models consisting of either 36 or 48 cellulose chains with DP = 40 that were investigated by MD simulations. The crystal 9314

dx.doi.org/10.1021/jp503535d | J. Phys. Chem. B 2014, 118, 9313−9321

The Journal of Physical Chemistry B

Article

Table 1. Dimensions and Number of Molecules of the Crystal Models allomorph Iα Iβ II IIII a

label of crystal model

A36 B36-1/-2 A48 B48

base plane dimensionsa 6(010) 6(110) 6(010) 6(110) 8(010) 8(010)

× × × × × ×

6(100) 6(1−10) 6(1−10) 6(1−10) 6(1−10) 6(100)

no. of chains

DP

no. of waters

36 36 36 36 48 48

40

28350 36792 40116 38125/38571 41216 30492

Number of chains (lattice plane).

models were constructed according to the published atomic coordinates and lattice plane parameters of each corresponding cellulose allomorph.10,12,14,15 The shapes of the base planes of the two native cellulose crystal models, the generally accepted base plane shape for a real microfibril, involve hydrophilic faces mostly exposed to the surrounding solvent environment. Unlike the native cellulose crystals, there is no experimental information about the shapes of the cellulose II and IIII crystals. It has been observed that the behavior of finite crystal models during MD varies quantitatively depending on the base plane shape rather than the dimensions.36,38 Two types of base planes, designated either A or B, were designed for both of the cellulose II and IIII crystal models to construct the layered structures of both of the possible chain sheets. For the cellulose II crystal models, the B-type models were subclassified as B-1 and B-2, consisting of the “parallel-up” (010) and “paralleldown” (020) chain sheets, respectively, at the central diagonal chain sheet. The crystal models were placed in a rectangular periodic box containing 28000−42000 TIP3P water molecules,44 depending on the type of crystal model. The buffer size between the edges of the crystal model and the periodic box was set to 15 Å. Table 1 summarizes the components of the crystal models, including the dimensions, the number of molecules, and the labels assigned to the models. Similar procedures to those described in our previous reports24,25,35,36 were used to perform the minimization and dynamics simulations of the crystal models. The solvated crystal models were minimized by a two-step procedure consisting of partial optimization of the water configurations and then full optimization of the whole system. The water configurations were equilibrated using the NVT ensemble for 300 ps with gradual increase in the temperature from 20 to 300 K. The production NPT dynamic runs at 300 K and 1 bar were carried out for 20 ns. Throughout the 300 ps of the equilibration NVT simulations, the positions of the heavy atoms of the solute crystal model were constrained with a relatively weak force of 10 kcal/(mol Å2). All of the MD simulations were performed using the PMEMD module of the AMBER 12 package45 combined with the Glycam0639 carbohydrate parameter set. The SHAKE option46 coupled with a 2 fs time step for integration and the particle mesh Ewald method47 for longrange nonbonded interactions were used throughout the MD simulations. The 1−4 scaling factors for the electrostatic and nonbonded interactions were set to one. The MD trajectories were analyzed by the PTRAJ module of the AMBER 12 package45 and visualized using VMD 1.9.1 software.48 Molecular graphics were created using PyMOL 1.6.0.0 software (Schödinger, LLC).49 The deformation or fiber twist of the crystal models during the dynamics simulations was estimated by the degree of deviation from planarity of the central chain sheet. As shown in Figure 2, the virtual bonds connecting the gravity centers of the

Figure 2. Molecular chain sheet along with the residue position labeling. The chain labeling is arbitrary among the chain sheets. Sheettwisting torsion angles are defined by the virtual bonds connecting the centers of gravity (G) of the residue(s): θ+20 = G(a,+20) − G[(c,+20), (d,+20)] − G[(c,−20), (d,−20)] − G( f,−20), θ−20 = G(a,−20) − G[(c,−20), (d,−20)] − G[(c,+20), (d,+20)] − G(e,+20), and so on. The solid line in blue represents the virtual bonds defining θ+20 and in the solid red line represents θ−20.

residues define the torsion angle (θ) representing the twisting angle of the chain sheet along the fiber axis.24,25 In the MD crystal models, the θ value was defined as the average of the θ values observed between the positive and negative residue positions. Obviously, a proportion of θ depends on the residue position along the fiber axis. The largest θ value would be expected at either of the terminal ends (θ±20), and the angles decrease to nearly 0° on approaching the middle position or θ±1 if an overall uniform twist occurs. When the chain sheet twists in the clockwise direction to give a right-handed twist, the θ value is defined as being positive. 9315

dx.doi.org/10.1021/jp503535d | J. Phys. Chem. B 2014, 118, 9313−9321

The Journal of Physical Chemistry B

■

Article

RESULTS AND DISCUSSION DFT Calculations of Isolated Cellulose Chain Sheets. Figure 3 shows the six superimposed structures of the initial

The amount of twisting of both of the cellulose II chain sheet models was larger than that of the native cellulose chain sheet models, as is clear from their cross-section views in Figure 3. In addition to the minor differences in the hydrogen-bonding schemes, the initial arrangements of molecular chains are significantly different between the two chain sheets; viewing them from the reducing end to the nonreducing end, a molecular chain tilts clockwise for the (010) chain sheet and counterclockwise for the (020) chain sheet (see the corresponding cross-section views of the initial structures in Figure 3). This is likely to be the cause of the opposite twisting of the two chain sheet models. In contrast to the (010) and (020) models with parallel chain polarity, the DFT-optimized structure of the (110) model, consisting of the diagonally aligned antiparallel chains, completely lost the original sheet form. The DFT-optimized structure of the cellulose IIII (1−10) model essentially retained the sheet form and showed the smallest deviation from the original sheet form among the six chain sheet models investigated, which indicates that the chain sheet structure is located close to the potential energy minimum even without the intersheet interactions occurring in the crystal packing environment. We have previously reported that the isolated cellulose IIII (100) model (the other chain sheet of cellulose IIII and not examined in the present study) was unstable by DFT optimization and spontaneously rolled into a tubular form.38,50 Figure 4 shows the variation of the sheet twisting angle (θ) of the B3LYP and CAM-B3LYP DFT-optimized structures for the

Figure 3. Superimposed structures of the original crystal (blue) and CAM-B3LYP DFT-optimized (red) 4× cello-octamers of the cellulose chain sheet models.

and CAM-B3LYP DFT-optimized 4× cello-octamer chain sheet models of the cellulose Iα (110), Iβ (100), II (010), II(020), II (110), and IIII (1−10) models. The DFT-optimized structures of the 4× hexamer chain models were similar to those of the corresponding 4× cello-octamer chain sheet models. All of the chain sheet models except for the cellulose II (110) model retained a sheet form after the DFT optimizations, with only slight deviations of the glycosidic linkage conformations of ϕ (O5−C1−O4−C4) and ψ (C1−O4−C4−C5) from the initial conformations (5°−10°). In the five chain sheet models that retained the sheet form (all but the cellulose II (110) model), the orientations of hydroxymethyl groups (ω, O5−C5−C6−O6) did not significantly rotate from the initial conformation of tg (ω ≈ 180°) for the two native cellulose crystal structures and gt (ω ≈ 60°) for the cellulose II and IIII crystal structures. However, the two hydroxymethyl groups of the reducing end residues in both of the CAM-B3LYP DFT-optimized 4× cello-hexamer and cellooctamer structures of the cellulose II (010) models rotated to a near tg conformation. The ϕ, ψ, and ω values of the CAMB3LYP and B3LYP DFT-optimized structures of the five chain sheet models are given in the Supporting Information (Tables S1 and S2). Both of the DFT-optimized native cellulose (cellulose Iα and Iβ) chain sheet models developed a right-handed twist, as was observed in the solvated crystal models,24−26 with a similar amount of twisting. We previously reported that the cellulose Iβ chain sheet model formed the right-handed twist using DFT calculations with the B3LYP functional.38 Conversely, the cellulose II (010) and (020) models twisted in opposite directions with right- and left-handed chirality, respectively.

Figure 4. Variation of the sheet twist angle (θ) of the B3LYP and CAM-B3LYP DFT-optimized structures of the five 4× cello-hexamer chain sheet models with respect to residue position.

five 4× cello-hexamer chain sheet models retaining the sheet form with respect to the residue positions. The two sets of the θ curves for the different DFT functionals are similar for the native and IIII allomorphs, but there are relatively large differences for the two cellulose II chain sheet models. The cellulose Iα and Iβ chain sheet models exhibit the most symmetrical and smooth θ curves, indicating uniform twisting. More importantly, the θ curves of the two cellulose I chain sheet models in each graph are similar. This suggests that the structures of the isolated cellulose Iα and Iβ chain sheets are located near the same optimized structure of the twisted form on the potential energy surface. The θ curve of the cellulose II (010) model shows the largest twist for the optimized structure as well as asymmetrical deformation. Figure 5 shows the θ variations of the CAM-B3LYP DFToptimized 4× cello-octamer structures of the five chain sheet models. The same general features are observed in the θ curves as those of the 4× cello-hexamer structures in Figure 4, but there are minor differences between the cello-hexamer and 9316

dx.doi.org/10.1021/jp503535d | J. Phys. Chem. B 2014, 118, 9313−9321

The Journal of Physical Chemistry B

Article

sheet form. The two cellulose I chain sheet models have similar isolated chain (ΔE1) and binding (ΔEbind) energies and, consequently, chain sheet (ΔEsheet) energies as a result of optimizing into similar twisted structures. These ΔEsheet values for the two cellulose I structures are much lower than the rest of the chain sheet models (by >5 kcal/mol per residue). The ΔEsheet value of the IIII (1−10) model is the third lowest, followed by the energies of the cellulose II (020) and (010) models. The ΔE1 values are virtually negligible (