FULL PAPER

WWW.C-CHEM.ORG

Solvents Effects on the Mechanism of Cellulose Hydrolysis: A QM/MM Study Claudia Loerbroks,[a] Andreas Heimermann,[b] and Walter Thiel*[a] This article reports a combined quantum mechanics/molecular mechanics (QM/MM) investigation on the acid hydrolysis of cellulose in water using two different models, cellobiose and a 40-unit cellulose chain. The explicitly treated solvent molecules strongly influence the conformations, intramolecular hydrogen bonds, and exoanomeric effects in these models. As these features are largely responsible for the barrier to cellulose hydrolysis, the present QM/MM results for the pathways and reaction intermediates in water are expected to be more real-

istic than those from a former density functional theory (DFT) study with implicit solvent (CPCM). However, in a qualitative sense, there is reasonable agreement between the DFT/CPCM and QM/MM predictions for the reaction mechanism. Differences arise mainly from specific solute–solvent hydrogen bonds that are only captured by QM/MM and not by DFT/CPCM. C 2015 Wiley Periodicals, Inc. V

Introduction

gen bonding situation, as HO(Y)AC(1) interactions can increase the acceptor strength of r*(O(Y)AC(1)) and OO(X) interactions can increase the donor strength of n(O(X)). Several previous theoretical studies have addressed the mechanism of cellulose hydrolysis. They made use of different computational methods (ab initio,[34–38] density functional theory (DFT),[39,40] and molecular dynamics (MD)[25,40–45]) and considered different models (dimethoxymethane,[34] 2-methoxytetrahydropyran,[35] 2-oxanol,[37] 2-methoxyoxane,[38] and cellobiose[33,42]) and solvents (water,[33] alcohols,[39] and ionic liquids[43,46,47]). Computational and experimental kinetic studies proposed a stepwise A1 mechanism via a carbocation,[33,48–50] which involves an initial protonation, followed by a conformational change of the glucose unit, the cleavage of the glycosidic linkage, and the addition of water to the anomeric carbon (Scheme 2). In this mechanism, the puckering of the nonreducing glucose ring attached to the C(1)AO(1) bond, the presence of inter- or intramolecular hydrogen bonds, and the anomeric effect were deemed important.[5,33,51,52] There are only few previous theoretical studies that address the influence of the solvent on the hydrolysis mechanism.[40,46,53] For acid hydrolysis, water and ionic liquids are the two most important solvents. Conformational analyses suggest that

Acid hydrolysis of cellulose is a process to utilize biomass in chemical industry, for example for the production of platform molecules like sorbitol and hydroxymethylfurfural (Scheme 1). It has therefore received much attention in recent research.[1–7] Cellulose hydrolysis is not a facile reaction:[1,2,8–15] in water it has to overcome a reaction barrier of 30–40 kcal mol 21 and thus needs to be carried out at high temperatures with strong acids (pKa < 23).[16–22] Since the first usage of biomass hydrolysis, several improvements of reaction conditions were made, for example, by changing the solvent from water to ionic liquids,[23–25] by changing the catalyst from strong acids like sulphuric acid to solid phase and metal catalysts,[12,18] and by introducing pretreatments like ball milling.[26,27] The necessity for harsh reaction conditions is commonly attributed to the rigid structure of cellulose. In this work, we focus on an abundant stock of biomass: cellulose Ib (Fig. 1). It consists of 1,4-b-connected glucose chains, held together by hydrogen bonds and van der Waals forces.[28–31] The glycosidic linkage C(1)AO(1) that is broken in acid hydrolysis is surrounded by intramolecular hydrogen bonds [O(30 )HO(5), O(2)HO(60 ), and O(2)HO(60 )] as well as by the intermolecular hydrogen bond O(6)HO(300 ) that connects the chains.[28–31] In this three-dimensional structure electronic effects also contribute to the rigidity of cellulose, in particular the anomeric effect in the acetal moiety O(5)AC(1)AO(1) (Fig. 2). In general, the anomeric effect describes the delocalisation of electronic density from a lone-pair orbital n of atom X to an antibonding orbital r* of the bond Y-Z (n(X) ! r*(Y-Z)), which is aligned antiperiplanar to the lone-pair. In cellulose there are two types of such interactions: the endoanomeric effect [n(O(5)] ! r*[C(1)AO(1)] that elongates the C(1)AO(1) bond, and the exoanomeric effect [n(O(1)]!r*[C(1)AO(5)] that shortens the C(1)AO(1) bond and therefore hinders cellulose hydrolysis. The exoanomeric effect is dominant in cellulose.[32,33] These anomeric effects are regulated not only by the antiperiplanar alignment, but also by the hydro1114

Journal of Computational Chemistry 2015, 36, 1114–1123

DOI: 10.1002/jcc.23898

[a] Claudia Loerbroks, W. Thiel Max-Planck-Institut f€ ur Kohlenforschung, Kaiser-Wilhelm-Platz 1, 45470 M€ ulheim an der Ruhr, Germany E-mail: [email protected] [b] A. Heimermann Theoretische Chemie, Technische Universit€ at Kaiserslautern, Erwin-Schr€ odinger-Str. 52, 67663 Kaiserslautern, Germany Contract grant sponsor: Research Cluster “Sustainable Chemical Synthesis—A Systems Approach (SusChemSys)”; Contract grant sponsor: European Regional Development Fund (ERDF) and the state of North Rhine-Westphalia, Germany, under the Operational Program “Regional Competitiveness and Employment” 2007–2013 (SusChemSys); Contract grant sponsor: Deutsche Forschungsgemeinschaft; Contract grant number: Cluster of Excellence RESOLV (EXC 1069) C 2015 Wiley Periodicals, Inc. V

WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

Scheme 1. Different pathways of biomass utilization.

minimum conformations in water are determined by the intramolecular hydrogen bonds of the model, especially O(2)O(60 ) and O(30 )HO(5), and their competition with intermolecular hydrogen bonds to solvent molecules.[54] The O(30 )HO(5) hydrogen bond is persistent and keeps the dihedral O(5)AC(1)AO(1)AC(4)0 in a minimum conformation between 270 to 280 , while the O(2)H and O(6)0 H groups are often interacting with water molecules rather than with each other.[52,54–57] In water puckered ring structures like boat conformations (Fig. 3)[58] are seldom found, but they play an important role in the acid hydrolysis of cellulose and in enzyme-catalysed depolymerisations of cellulose.[45,51,59,60] Glycosides behave differently in ionic liquids. Puckered conformations become more frequent and negatively charged ions are positioned close to the hydroxyl groups, disrupting hydrogen bonds, according to experimental and computational studies.[3,43,47] This is consistent with the higher polarity of ionic liquids, which is believed to affect the strength of anomeric effects and can thus enhance or diminish the conformational rigidity of glycosides.[61–67] Experimental studies indicate that acid hydrolysis in ionic liquids requires slightly less activation than in water and allows for milder reaction conditions regarding the choice of acid strength and the dissolution of cellulose.[68] In this work, we investigate the influence of the solvent on the mechanism and on the structural and electronic barriers to cellulose hydrolysis in atomistic detail. Rather than studying cellobiose in implicit solvent as before,[33] we treat the solvent

Figure 1. Part of the cellullose Ib structure with O(2)HO(60 ) (top) and O(2)HO(60 ) (bottom) hydrogen bonds. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

Figure 2. Orbital interactions involved in the most important exoanomeric and endoanomeric effects in cellobiose. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

explicitly and also investigate a larger cellulose chain model. We discuss the importance of explicit solvent and model size for the mechanism.

Methods We used the NAMD program suite for the MD runs. All simula˚ for the nonbonded tions employed a cutoff radius of 12 A interactions and the particle mesh Ewald summation method for calculating the electrostatic potential with periodic boundary conditions.[69] We used the following procedure to prepare

Scheme 2. A1 mechanism for cellulose hydrolysis via a cyclic carbocation.

Journal of Computational Chemistry 2015, 36, 1114–1123

1115

FULL PAPER

WWW.C-CHEM.ORG

Figure 3. Possible conformations of the nonreducing ring of cellulose; the notations of spiroaminal conformers in this paper follow the recommendations for comparable sugar compounds.[58] Hydrogens and hydroxyl groups are omitted for clarity.

the MD production runs: (a) minimising the energy; (b) heating to 300 K (water) with Langevin temperature control; (c) equilibration for 20 ns in the NPT ensemble with Langevin piston pressure control (1 atm); and (d) production run of 100 ns for each model. The steps from (a) to (c) were first carried out for the solvent box and were then repeated with the glycoside model inside the box (without restraints). The time step in production runs was 2 fs, and energy, volume, and pressure were saved every 2 ps as a snapshot. The QM/MM calculations were done with the package ChemShell,[70] which calls the programs Gaussian09[71] and TURBOMOLE[72] (BP86[73–77] and BB1K[77]) for computing QM energies and gradients, and DL_POLY[78] (GLYCAM06,[79] and TIP3P) for computing MM energies and gradients. Hydrogen link atoms and the charge-shift model were used to treat the QM/MM boundary, and the QM/MM interactions were handled ˚ was applied by electronic embedding.[80,81] A cutoff of 1000 A for the nonbonding interactions. Six snapshots for cellobiose and for all three regions of the cellulose chain were selected from the MD production runs. For each of them, a proton was added close to O(1) and a region of 35 A˚ around the O(1) atom of the glycosidic bond was cut out. In the resulting structures, the inner part (radius ˚ ) was taken as active region, while the outer part was of 25 A kept fixed in the subsequent QM/MM calculations. In water, the QM region contained 55 atoms (small; two glucose units, two water molecules, one hydronium ion) or 79 atoms (large; two glucose units, ten water molecules, one hydronium ion). The charge of the QM region was 11. The chosen QM region and the hydrogen link atom scheme were also advised in the literature.[82] The selected snapshots were first optimised and then scanned along possible reaction pathways with a step size of 1116

Journal of Computational Chemistry 2015, 36, 1114–1123

0.05 Bohr. The optimizations were carried out with the HDLCopt[83] and the DL-FIND[84] modules, employing the limited-memory quasi-Newton algorithm (L-BFGS) for energy minimisations and the P-RFO algorithm for transition state searches. The scans were performed at the QM(BP86/SVP)/ GLYCAM06, TIP3P level (with resolution-of-identity approximation for BP86), while stationary points and frequencies were calculated at the QM(BB1K/TZVP)/GLYCAM06, TIP3P level. The proper connection between transition states (one imaginary mode) and minima was verified by careful scans. Intrinsic reaction coordinates could not be used due to low imaginary frequencies. This problem had also been encountered in other studies.[33] Natural charges and electron delocalisation were determined by natural bonding orbital (NBO) analysis of the QM region using Gaussian09.[85] The cellulose models, cellobiose and the 40-unit glucose chain, were generated with the cellulose builder.[86] The solvent box had a volume of 503 A˚3 for cellobiose and of ˚ 3 for the cellulose chain. The sugar molecules 50*50*265 A were represented by the GLYCAM06 force field, which is known to give good results in comparison with other force fields.[79] To establish a common energy scale, QM-only single-point calculations (BB1K/TZVP, CPCM(water)) of the QM region with hydrogen link atoms were performed for the stationary points INT1 (reactant) and INT3 (product). One single snapshot was chosen as an energy reference for all other snapshots, and the relative energies of the relevant intermediates were calculated. This procedure yielded two reaction pathways for each snapshot, one each with reference to reactant and to product. As all pathways refer to the same reference, the relative energies of a given stationary point obtained from the two pathways can be averaged for each snapshot, resulting in pathways which are comparable, even though the active MM region around the QM region may be different.[87] In the following section, we compare the present QM/MM results with previously reported results from a DFT/CPCM study,[33] which made use of the BB1K functional and the 6– 3111G** basis in the DFT calculations and of the CPCM implicit solvent model for describing water.

Results and Discussion Hydrolysis in water In our QM/MM study on the mechanism of acid hydrolysis we employed two cellulose models: cellobiose and a 40-unit glucose cellulose chain. For the chain we were especially interested in which part is cleaved first, so we split the chain into three regions (Fig. 4): region N (nonreducing end), region R (reducing end), and region C (center of the chain). We took six starting structures for cellobiose and six structures for each region in the cellulose chain (named a to f) and performed QM/MM calculations with two QM regions (Fig. 5): a small region containing two connected glucose units with one hydronium ion and the minimum number of water molecules (2 H2O, H3O1) and a larger region additionally including the WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

Figure 4. Different regions [nonreducing end (N), reducing end (R), and centre of chain (C)] in one cellulose chain.

Figure 5. QM regions for cellobiose (snapshot b) in water: small (left) and large (right).

water molecules connected to those of the minimal QM region around the glycosidic bond (10 H2O, H3O1). For all these setups, the reaction pathway of acid hydrolysis was computed and an NBO analysis was carried out for the starting structure. We now discuss the influence of the explicit solvent on (a) the pathways and geometries, (b) the energies, and (c) the relation between structures and reaction barriers, and compare with the previous DFT results.[33] Pathways and geometries A comparison between the results obtained with implicit and explicit solvation is useful to assess the merits and shortcom-

Scheme 3. Mechanism of acid hydrolysis of cellobiose with implicit (left) and explicit solvent (right). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

ings of these solvation models and to analyse the role of the solvent for the reaction. In Scheme 3 and Table 1 we display the mechanistic features obtained with these two approaches. Detailed numerical results are given in Tables 2-6. In the DFT/CPCM and QM/MM calculations in water, we find both A1 and A2 pathways: A1 is a stepwise mechanism, which features a cyclic carbocation (INT2) as a stationary point.A2 is a concerted mechanism, which directly connects reactant INT1 and product INT3. Fewer stationary points are present when using explicit rather than implicit solvation. At the QM/MM level, protonated cellobiose reactants are rarely stable and nonchair conformers without a broken C(1)AO(1) bond cannot be found. The nonchair conformations located for transition states and carbocations in explicit solvent are not only similar to those in implicit solvent (e.g., 2,5B), but also include novel conformers like E3. The optimized structures are similar for both solvent models (Table 1). Although the QM/MM optimized distances and dihedrals cover a greater range of values (i.e., larger distances and more flexible dihedrals), the DFT values always lie within this range. The same hydrogen bonds are normally found in DFT and QM/MM, although in the QM/MM calculations with the larger QM region intramolecular hydrogen bonds are easily replaced by intermolecular hydrogen bonds with water. For example, in Figure 5 the O(30 )HO(5) hydrogen bond (top, left) is broken when increasing the number of QM water molecules, with the proton being transferred away from O(1) (top, right). The only major discrepancies involve barrierless pathways (B) and cases of no reaction (N) that are sometimes Journal of Computational Chemistry 2015, 36, 1114–1123

1117

FULL PAPER

WWW.C-CHEM.ORG

Table 1. Hydrolysis of cellobiose in implicit (DFT/CPCM) and explicit (QM/MM) solvent: similarities and differences. QM/MM DFT/CPCM Intermediates Types Conformers C(1)-O(1) (INT1)[a] E(2) exo (INT1)[b] u (INT1)[c] C(1)-O(1) (TS1)[a] O-C(1) (TS1)[a]

Small

prot., non-chair, INT2 A1, A2 4 H3, 2,5B, B2,5; 4E 1.38 18 292 1.87–1.96 2.99

prot., INT2 A1, A2 2,5 B, E3 1.38–1.50 6–17 2134 to 259 1.74–2.60 1.99–4.02

Large

Exp. [[48,49]

INT2 A1[48,49]

1.40[88,89]

1.39–1.44 7–14 2107 to 258 1.88–2.57 1.90–3.11

273 to 280[52,56,57]

˚ . [b] kcal mol21. [c]  . [a] A

encountered in QM/MM and not in DFT/CPCM. At the QM/MM level, a barrierless transition from cellobiose to two glucose units may happen for nonchair conformers with a large endoanomeric effect (e.g., cellobiose-large-e, C-d, R-c), and there may be no reaction if the proton is located at O(5) or too far away from O(1) for a protonation to occur (e.g., R-d, R-b-large, R-dlarge). In the QM/MM calculations, the major structural features are similar for cellobiose and for the cellulose chain in explicit solvent (Tables 5 and 6). To summarise the main differences between QM/MM and DFT/CPCM results: in explicit solvent there are fewer intermediates and the flexibility of the existing conformers is higher. To identify the reasons for these differences, we checked the hydrogen bonds and the interactions of cellulose with explicit solvent (Tables 5 and 6). Protonated cellulose is in most cases not a stable minimum in QM/MM, as the QM water molecules surrounding O(1) are more basic than O(1) and abstract the proton in a barrierless process.[18,33] This cannot be modeled in DFT/CPCM calculations with implicit solvent. At the QM/MM level, an intermediate that is protonated at O(1) can occur only if no water molecule is close enough to abstract it (cellobiose-d, cellobiose-e, C-d, cellobiose-a-large, Nb-large). Moreover, structures of the flipped nonreducing glucose unit with an unbroken C(1)AO(1) bond are not stable in explicit solvent. Already in our previous DFT study,[33] such nonchair minima were very shallow, as the energy differences to the next transition state were below 1 kcal mol21. The QM/ MM results are consistent with the literature, as flipped conformers are known to be high in energy.[90,91] For the cyclic carbocation intermediate, we find similar conformers in DFT and QM/MM (2,5B), but mostly those with C(3) below the ring (E3); the E3 conformer is highly populated because C(3) carries no large group but instead a hydroxyl group that is easy to flip, and because the C(6) hydroxyl group in E3 is not in the unfavorable axial position.

Table 2. Hydrolysis of cellobiose in explicit solvent: QM/MM energies for small and large QM region and experimental data (kcal mol21).

Reaction barrier (range) Reaction energy (range)

1118

Small

Large

24 (14–31) 10 (3–14)

21 (15–26) 2 (24 to 6)

Journal of Computational Chemistry 2015, 36, 1114–1123

Exp.[[92–94] 30–40 23

We also examined the starting structures with regard to conformational freedom (Tables 5 and 6). We find that intra- and intermolecular hydrogen bonds impact different distances and dihedrals in the reactant. O(1)H2O hydrogen bonds are associated with long C(1)AO(1) bonds and weak exoanomeric effects (e.g., in cellobiose-e), while in the absence of O(1)H2O interactions, O(30 )HO(5) hydrogen bonds lead to short C(1)AO(1) bonds and large exoanomeric effects (e.g., in cellobiose-c). The strong influence of solute–solvent hydrogen bonds cannot be captured in implicit solvent.[33] The variation of O(5)AC(1)AO(1)AC(40 ) dihedrals also seems related to the stabilization by explicit solvent. The less intramolecular hydrogen bonds are present, the more flexible is the dihedral (u5 2160 to 260 , e.g., N-d), while O(30 )HO(5) hydrogen bonds reduce the conformational freedom (u5 2100 to 280 , e.g., N-e). Explicit solvent has a strong impact on which pathway is chosen. Whether an A1 or an A2 pathway is taken is mostly determined by the position of the water molecules in the starting structures. If water molecules are positioned within ˚ under C(1), an A2 mechanism is more likely than A1 3.8 A (Tables 5 and 6). In the literature, structures with water molecules under the ring, and therefore close to C(1), are considered less likely, as the ring is hydrophobic and the water molecules are more attracted by the hydrophilic hydroxyl groups.[2] This explains why the A1 mechanism is found to be preferred for cellulose hydrolysis in kinetic experiments. Energies In our previous DFT study on cellobiose, the initial protonation was the energetically most demanding step (uphill by 28 kcal Table 3. Overall QM/MM reaction barriers in explicit solvent: DE(TS1), average values (kcal mol21) for small and large QM region (see text).

All A1 A2 Cellobiose Region N Region C Region R

Small

Large

28 29 28 26 28 30 29

30 30 30 23 32 33 31

WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

Table 4. QM/MM energies (kcal mol21) relative to QM single-point calculations for INT1 and INT3 in water. Water-small Structure Cellobiose

Region N

Region C

Region R

INT1 a b c d e f a b c d e f a b c d e f a b c d e f

15.8 13.1 23.3 15.3 18.0 6.7 5.4 16.0 12.8 12.7 12.1 9.6 19.8 25.4 13.4 44.0 15.8 15.9 13.3 11.0 17.1 10.3 2.3 10.0

TS1

INT2

28.6 22.6 30.5 25.3 20.8

28.5

24.1 19.8

21.7 31.4 31.1 29.4 26.7 30.0 33.1 38.5

Water-large TS2

INT3

INT1

TS1

INT2

14.0

6.2 13.9 3.3 13.9 9.6 12.9

1.7 6.3 24.2

25.5 22.8

20.4

3.8 26.4

21.8

6.4

18.9 18.1 18.2 14.3 4.2 7.3 19.8 17.0 22.1 23.3 9.2 16.3 9.9 33.8

19.0 21.8 8.2 23.1 4.6 14.1 0.7 12.7 41.0 4.3 5.9 21.4 1.4

25.1 33.9 39.0 30.1 34.1 36.4

9.2 10.4 25.5 4.5 14.0

22.9

24.8 28.1 42.9

29.1 20.2 29.2 22.4

26.9 30.4

INT3

14.8

3.6 23.8

30.6 38.8 42.6

38.0 45.5

15.6 34.1 13.7 0.0 9.0

28.9 24.8 28.9

12.6 11.7 4.5

18.7 7.7

TS2

22.9

40.2

11.2

Empty entries indicate that no results are available.

mol21). A1 was the preferred mechanism with an overall free energy barrier of 31 kcal mol21 (A2 38 kcal mol21) and a reaction free energy of 23 kcal mol21.

Tables 2–4 list the computed QM/MM energies for cellobiose and the cellulose chain in explicit solvent. The QM/MM barriers for the different snapshots range from 14 to 31 kcal mol21

Table 5. Geometries, energies in kcal mol21, and hydrogen bonding (HBB) situation for all QM/MM models (small QM region) in water. E(2)

[a]

TS1

DE(TS1-INT1)[a] exo endo C(1)-O(1) Cellobiose a b c d e f Region N a b c d e f Region C a b c d e f Region R a b c d e f

– 15.2 29.4 10.3 8.5 19.1 – 12.9 17.3 17.3 18.0 21.2 15.2 10.7 9.6 – 8.1 3.6 24.7 – – – 24.8 21.5

7 11 17 8 8 6 15 11 3 3 7 7 5 11 8 4 6 4 8 6 7 21 17 9

4 4 2 4 30 4 3 4 4 4 4 4 4 4 5 48 4 5 3 2 31 1 2 4

– 1.79 1.88 1.74 2.00 1.93 – 1.81 1.94 1.93 1.95 1.98 1.93 1.71 1.91 – 2.02 1.89 2.18 – – – 1.97 2.03

[b]

INT1 O- (C1) – 3.27 2.72 2.98 4.02 3.29 – 3.15 2.97 2.91 3.11 2.79 3.71 2.83 4.32 – 2.83 3.29 2.47 – – – 2.66 2.58

[b]

C(1)-O(1) 1.41 1.41 1.38 1.45 1.50 1.41 1.38 1.42 1.43 1.43 1.42 1.42 1.43 1.43 1.48 1.46 1.41 1.43 1.41 1.39 1.51 1.35 1.39 1.41

[b]

O- (C1) 3.62 3.80 3.37 3.53 4.49 4.12 4.57 3.75 3.33 3.33 3.51 3.84 4.20 3.11 4.55 4.30 3.33 3.54 3.27 3.51 3.32 4.94 3.12 3.42

H bonds in INT1 [b]

u/



2134 275 288 260 285 259 272 276 2127 2142 299 286 154 282 297 285 2118 258 277 2147 263 2117 2100 263

Intra O(1) / O(5) 3 1 5 5 5 5 5 5 3 3 1 5 2 5 3 5 5 4 1 5 5 3 5 5

3 3 5 4 4 3 2 3 3 3 3 3 3 3 3 4 3 3 3 5 5 5 5 3

Nonchair Type of conformation mechanism TS1/INT2 2 1 2 2 1 1 4 2 2 2 2 2 1 2 1 3 2 2 2 2 3 4 2 2

2,5

B E3 E3 E3 2,5 B flip, E3 O S2 (INT1) E3 O C3 E3 E3 E3 flip, E3 4 E E3 E3 (INT1) E3 E3 E3 E3 2,5 B 4 C1 (INT1) E3 E3

Number codes: intramolecular HBB (intra) 1 5 O(30 )HO(5), 2 5 O(2)O(30 ), 3 5 O(2)O(60 ), 4 5 O(30 )O(6), 5 5 none; HBB to O(1) or O(5) 1 5 O(1)H2O, ˚. 2 5 none, 3 5 O(1)H3O1, 4 5 O(1)-H1, 5 5 O(5)-H1; type of mechanism (mech) 1 5 A1, 2 5 A2, 3 5 B, 4 5 N. [a] kcal mol21. [b] A

Journal of Computational Chemistry 2015, 36, 1114–1123

1119

FULL PAPER

WWW.C-CHEM.ORG

Table 6. Geometries, energies in kcal mol21, and hydrogen bonding (HBB) situation for all QM/MM models (large QM region) in water. E(2)

[a]

TS1

DE(TS1-INT1)[a] exo endo C(1)-O(1) Cellobiose a b c d e f Region N a b c d e f Region C a b c d e f Region R a b c d e f

[b]

INT1 O- (C1)

[b]

C(1)-O(1)

[b]

O- (C1)

H bonds in INT1 [b]

u/



intra O(1) / O(5)

Nonchair Type of conformation mechanism TS1/INT2 2,5

B E3 4 C1

– 24.5 26.8

7 14 17

16 3 2

– 2.05 1.88

– 3.11 2.73

1.44 1.40 1.39

3.73 4.25 3.35

2107 272 288

3 5 5

4 1 2

2 1 2

– 25.0

14 11

19 4

– 1.96

– 2.72

1.41 1.40

5.25 3.40

288 258

5 5

3 3

3 2

2,5

5.4 40.5 18.4 30.6 21.1 25.7 – 12.1 2.0 11.5 – 32.4 29.1

8 13 5 14 3 5 16 14 2 9 8 16 5

5 4 3 3 4 4 3 3 6 4 5 3 4

1.82 1.83 1.97 1.93 1.92 1.89 – – – 2.02 – 2.14 2.01

3.24 2.55 3.07 3.36 3.17 4.9 – – – 2.83 – 2.66 2.80

1.46 1.41 1.41 1.39 1.41 1.41 1.38 1.4 1.50 1.41 1.41 1.38 1.41

3.24 3.40 3.54 3.59 1.39 4.46 1.40 4.42 3.67 3.62 3.72 3.47 3.56

269 282 2150 2100 288 163 283 296 2107 2117 261 283 2154

5 5 3 1 5 2 1 3 5 5 4 1 5

4 3 1 1 3 3 2 3 4 3 3 2 3

2 2 2 2 2 1 4 1 1 2 3 2 1

E3 E3 E3 O C3 E3 E3 2 C5 (INT1)

– – 35.9

9 17 9

3 2 3

– – 2.18

– – 2.77

1.39 1.38 1.40

4.94 4.40 3.43

2135 298 257

3 5 5

5 5 1

4 4 2

B (INT1) flip, E3

4 E E3 4 C1 (INT1) E3 E3

4

4 C1 C1 (INT1) E3

Number codes: intramolecular HBB (intra) 1 5 O(30 )HO(5), 2 5 O(2)O(30 ), 3 5 O(2)O(60 ), 4 5 O(30 )O(6), 5 5 none; HBB to O(1) or O(5) 1 5 O(1)H2O, 2 5 none, 3 5 O(1)H3O1, 4 5 O(1)-H1, 5 5 O(5)-H1; type of mechanism (mech) 1 5 A1, 2 5 A2, 3 5 B, 4 5 N. [a] kcal mol21 [b] A˚.

and are thus lower than the experimental values. Unlike DFT we find no clear preference for the A1 or A2 mechanism. While this seems contrary to the experimental results (Table 1) we have to remember that our QM/MM calculations only yield electronic energies DE. Our previous DFT study has already shown that entropy has a large impact: the rate-determining A1 and A2 transition states are nearly isoenergetic without considering entropy [DDH(A1 versus A2) 5 0.6 kcal mol21] but differ strongly in free energy [DDG(A1 versus A2) 5 7 kcal mol21].[33] At the QM/MM level (Scheme 3), protonation is not a separate step, and the concerted transformation from INT1 to INT2/INT3 via TS1, which involves protonation, conformational change, and C(1)AO(1) bond cleavage, is ratedetermining in explicit solvent. Concerning reactivities (Table 3), the QM/MM calculations with the large QM region indicate that cellobiose is by 8 kcal mol21 easier to hydrolyse than the cellulose chain, while there are no striking differences between the different regions of the chain. This is consistent with our previous metadynamics simulations for the same cellulose models, which gave no conformational differences between the regions of the chain and cellobiose in water.[95] Relation between geometry and energy Finally, we correlate the computed energies with the reactant and transition state structures, to determine the main 1120

Journal of Computational Chemistry 2015, 36, 1114–1123

obstacles to hydrolysis (in comparison with DFT/CPCM). In our preceding DFT study,[33] we had identified several such obstacles: the exoanomeric effect, hydrogen bonds, protonation, and conformational changes. As before, we gauge the strength of the anomeric effect from the associated secondorder perturbation energy lowering [E(2) exo] obtained from NBO analysis.[33] Figure 6 (top) shows a plot of DE(TS1-INT1) versus E(2) exo. In general, a larger exoanomeric effect leads to a higher energy barrier, as can best be seen for cellobiose (small and large, in blue). However, the corresponding QM/MM data points are widely scattered, and not for all models the trend is true. Some of these irregularities can be explained by unusual conformations that increase or reduce the barrier by effects unrelated to the E(2) exo value. For example, flipped conformers (Fig. 3, cellobiose-f, cellobiose-f-large, C-a, C-a-large) have higher barriers than expected from the strength of their anomeric effect. Our DFT calculations in implicit solvent indicated that a value of u around 290 leads to the highest possible exoanomeric effect.[33] This is confirmed by the QM/MM calculations with explicit solvent. However, Figure 6 (middle) also shows that u is not the main criterion, as dihedrals around 290 can have the highest E(2) exo values, but not necessarily so. Indeed, hydrogen bonds influence the exoanomeric effect even more: For example, contrary to our expectations regarding the correlation of u with E(2) exo, R-d has a high WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

hydrogen bond present in INT1. Most structures do not possess such intramolecular hydrogen bonds and vary in DE between 5 and 40 kcal mol21. In the subset of all hydrogenbonded structures, O(30 )HO(5) hydrogen bonds are associated with slightly higher barriers, probably because they increase the exoanomeric effect. This effect is stronger for the large QM region. For example, in N-e and R-a, which both have O(30 )HO(5) hydrogen bonds, we find for the small QM region DE(TS1-INT1) 5 18 or 25 kcal mol21 and E(2) exo 5 7 or 8 kcal mol21, respectively, while for the large region both the barrier and the exoanomeric effect increase (DE(TS1-INT1) 5 31 or 32 kcal mol21 and E(2) exo 5 14 or 16 kcal mol21). In our previous DFT/CPCM study, the O(2)O(60 ) intramolecular hydrogen bond was found to be important due to the high basicity of O(60 ); here, it is not associated with a high barrier, and hence a hindrance of protonation by O(2)O(60 ) is not confirmed by the present calculations.

Conclusions

Figure 6. Top: Correlation of the energy difference DE between TS1 and INT1 with the exoanomeric effect, E(2) exo. Middle: Correlation of the exoanomeric effect in INT1, E(2) exo, with the dihedral u (O(5)-C(1)-O(1)-C(40 )). Bottom: Correlation of the energy difference DE between TS1 and INT1 with the intramolecular hydrogen bonds in INT1.

exoanomeric effect at a dihedral of 2117 due to a strong O(30 )HO(5) hydrogen bond, and C-d (u 5 285 ) has a low exoanomeric effect because O(1) is protonated in INT1. Therefore, we further investigated the influence of the hydrogen bonds on the reaction barrier. In Figure 6 (bottom) we correlate DE(TS1-INT1) with the type of intramolecular

We investigated the acid hydrolysis of cellulose using the QM/ MM method and two different models: cellobiose and a 40unit cellulose chain. In the cellulose chain we examined three regions, the two ends of the chain and the center of the chain. Comparison to our earlier DFT/CPCM study shows that explicit solvation by water has a distinct influence on the reaction mechanism: energetically the A1 and A2 pathways now have similar overall QM/MM barriers, and only the reactant, the carbocation (in A1), and the product are stationary points. Whether an A1 or A2 pathway is taken depends on the solvent environment: short C(1)-water distances favor the A2 mechanism. As the ring containing C(1) is hydrophobic, such short C(1)-water distances are less likely and A1 is the preferred path. Hydrogen bonds strongly influence the mechanism, as they can stabilize different conformations and can reduce or increase the exoanomeric effect and hence the reaction barrier. In water we find no significant energetical preferences between hydrolysis at the end and in the middle of the cellulose chain. The QM/MM calculations with explicit solvent indicate that the obstacles to cellulose hydrolysis in water (hydrogen bonding, conformations, anomeric effects) cannot be regarded separately, as they are all interconnected. To improve acid cellulose hydrolysis, it is thus not sufficient to focus on only one of these factors but instead one has to address them all. The QM/MM mechanism confirms some of the important qualitative insights gained from the previous DFT study (pathways, conformations, energies) but provides a more diverse and presumably more realistic picture (stationary points, conformational freedom, A1 or A2 preference) due to the explicit treatment of the solvent molecules and their interactions with the solute, especially by hydrogen bonds. Keywords: cellulose  solvation  combined quantum mechanics molecular mechanics method  reaction mechanisms Journal of Computational Chemistry 2015, 36, 1114–1123

1121

FULL PAPER

WWW.C-CHEM.ORG

How to cite this article: C. Loerbroks, A. Heimermann, W. Thiel J. Comput. Chem.. 2015, 36, 1114–1123. DOI: 10.1002/jcc.23898

[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44]

1122

R. Rinaldi, F. Sch€ uth, Energy Environ. Sci. 2009, 2, 610. R. Rinaldi, F. Sch€ uth, ChemSusChem 2009, 2, 1096. H. M. Cho, A. S. Gross, J.-W. Chu, J. Am. Chem. Soc. 2011, 133, 14033. K. Kulasinski, S. Keten, S. Churakov, D. Derome, J. Carmeliet, Cellulose 2014, 21, 1103. M. Bergenstra˚hle, E. Thormann, N. Nordgren, L. A. Berglund, Langmuir 2009, 25, 4635. A. Pinkert, K. N. Marsh, S. Pang, Ind. Eng. Chem. Res. 2010, 49, 11121. W.-H. Hsu, Y.-Y. Lee, W.-H. Peng, K. C. W. Wu, Catal. Today 2011, 174, 65. G. T. Beckham, J. F. Matthews, B. Peters, Y. J. Bomble, M. E. Himmel, M. F. Crowley, J. Phys. Chem. B 2011, 115, 4118. C. M. Payne, M. E. Himmel, M. F. Crowley, G. T. Beckham, J. Phys. Chem. Lett. 2011, 2, 1546. S. P. S. Chundawat, G. T. Beckham, M. E. Himmel, B. E. Dale, Annu. Rev. Chem. Biomol. Eng. 2011, 2, 121. P. Dhepe, A. Fukuoka, ChemSusChem 2008, 1, 969. H. Kobayashi, T. Komanoya, S. K. Guha, K. Hara, A. Fukuoka, Appl. Catal. A: Gen. 2011, 409, 13. L. R. Lynd, C. E. Wyman, T. U. Gerngross, Biotechnol. Progr. 1999, 15, 777. S. Van de Vyver, J. Geboers, P. A. Jacobs, B. F. Sels, ChemCatChem 2011, 3, 82. J. Verendel, T. Church, P. Andersson, Synthesis 2011, 2011, 1649. C. Li, Q. Wang, Z. K. Zhao, Green Chem. 2008, 10, 177. C. Li, Z. K. Zhao, Adv. Synth. Catal. 2007, 349, 1847. R. Rinaldi, N. Meine, J. vom Stein, R. Palkovits, F. Sch€ uth, ChemSusChem 2010, 3, 266. R. Rinaldi, R. Palkovits, F. Sch€ uth, Angew. Chem. Int. Ed. 2008, 47, 8047. N. Meine, R. Rinaldi, F. Sch€ uth, ChemSusChem 2012, 5, 1449. L. Vanoye, M. Fanselow, J. D. Holbrey, M. P. Atkins, K. R. Seddon, Green Chem. 2009, 11, 390. J. B. Binder, R. T. Raines, Proc. Natl. Acad. Sci. USA 2010, 107, 4516. R. P. Swatloski, S. K. Spear, J. D. Holbrey, R. D. Rogers, J. Am. Chem. Soc. 2002, 124, 4974. R. Rinaldi, Chem. Commun. 2011, 47, 511. H. Liu, K. L. Sale, B. A. Simmons, S. Singh, J. Phys. Chem. B 2011, 115, 10251. H. Zhao, J. H. Kwak, Y. Wang, J. A. Franz, J. M. White, J. E. Holladay, Energy Fuels 2005, 20, 807. F. Sch€ uth, R. Rinaldi, N. Meine, M. K€aldstr€ om, J. Hilgert, M. D. K. Rechulski, Catal Today 2014, 234, 24. A. C. O’Sullivan, Cellulose 1997, 4, 173. Y. Nishiyama, P. Langan, H. Chanzy, J. Am. Chem. Soc. 2002, 124, 9074. Y. Nishiyama, J. Sugiyama, H. Chanzy, P. Langan, J. Am. Chem. Soc. 2003, 125, 14300. M. Jarvis, Nature 2003, 426, 611. G. A. Jeffrey, J. A. Pople, L. Radom, Carbohyd. Res. 1972, 25, 117. C. Loerbroks, R. Rinaldi, W. Thiel, Chem. Eur. J. 2013, 19, 16282. C. W. Andrews, J. P. Bowen, B. Fraser-Reid, J. Chem. Soc., Chem. Commun. 1989, 1913. C. W. Andrews, B. Fraser-Reid, J. P. Bowen, J. Am. Chem. Soc. 1991, 113, 8293. P. Deslongchamps, S. Li, Y. L. Dory, Org. Lett. 2004, 6, 505. B. J. Smith, J. Am. Chem. Soc. 1997, 119, 2699. B. J. Smith, J. Phys. Chem. A 1998, 102, 4728. J. M. Beck, S. M. Miller, M. W. Peczuh, C. M. Hadad, J. Org. Chem. 2012, 77, 4242. H. Satoh, H. S. Hansen, S. Manabe, W. F. van Gunsteren, P. H. H€ unenberger, J. Chem. Theory Comput. 2010, 6, 1783. H. Dong, M. R. Nimlos, M. E. Himmel, D. K. Johnson, X. Qian, J. Phys. Chem. A 2009, 113, 8577. X. Liang, A. Montoya, B. S. Haynes, J. Phys. Chem. B 2011, 115, 10682. H. Liu, K. L. Sale, B. M. Holmes, B. A. Simmons, S. Singh, J. Phys. Chem. B 2010, 114, 4293. F. A. Momany, U. Schnupf, Carbohyd. Res. 2011, 346, 619.

Journal of Computational Chemistry 2015, 36, 1114–1123

[45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58]

[59] [60] [61] [62] [63]

[64] [65] [66] [67] [68] [69] [70]

[71]

[72] [73] [74] [75] [76] [77] [78] [79] [80]

W. Plazinski, M. Drach, RSC Adv. 2014, 4, 25028. Y. Nishimura, D. Yokogawa, S. Irle, Chem. Phy. Lett. 2014, 603, 7. Z. Jarin, J. Pfaendtner, J. Chem. Theory Comput. 2014, 10, 507. T. E. Timell, Can. J. Chem. 1964, 42, 1456. S. R. Wann, M. M. Kreevoy, J. Org. Chem. 1981, 46, 419. J. N. Bemiller In Acid-Catalyzed Hydrolysis of Glycosides, Vol. 22; L. W. Melville, R. S. Tipson, Eds., Academic Press, 1967, pp. 25–108. H. B. Mayes, L. J. Broadbelt, G. T. Beckham, J. Am. Chem. Soc. 2013, 136, 1008. L. Peric-Hassler, H. S. Hansen, R. Baron, P. H. H€ unenberger, Carbohyd. Res. 2010, 345, 1781. K. L. Fleming, J. Pfaendtner, J. Phys. Chem. A 2013, 117, 14200. E. Hatcher, E. S€awen, G. Widmalm, A. D. MacKerell, J. Phys. Chem. B 2011, 115, 597. N. J. Christensen, P. I. Hansen, F. H. Larsen, T. Folkerman, M. S. Motawia, S. B. Engelsen, Carbohyd. Res. 2010, 345, 474. A. D. French, G. P. Johnson, Can. J. Chem. 2006, 84, 603. C. S. Pereira, D. Kony, R. Baron, M. M€ uller, W. F. van Gunsteren, P. H. H€ unenberger, Biophys. J. 2006, 90, 4337. (a) IUPAC-IUB Joint Commission on Biochemical Nomenclature (JCBN), Conformational nomenclature for five, six-membered ring forms of monosaccharides, their derivatives, Recommendations 1980 Arch. Biochem. Biophys. 1981, 207, 469; (b) Eur. J. Biochem. 1980, 111, 295; (c) Pure Appl. Chem. 1981, 53, 1901. C. B. Barnett, K. J. Naidoo, Mol. Phys. 2009, 107, 1243. B. M. Sattelle, S. U. Hansen, J. Gardiner, A. Almond, J. Am. Chem. Soc. 2010, 132, 13132. J. P. Praly, R. U. Lemieux, Can. J. Chem. 1987, 65, 213. R. U. Lemieux, Pure Appl. Chem. 1971, 25, 527. T. G. A. Youngs, J. D. Holbrey, C. L. Mullan, S. E. Norman, M. C. Lagunas, C. D’Agostino, M. D. Mantle, L. F. Gladden, D. T. Bowron, C. Hardacre, Chem. Sci. 2011, 2, 1594. Y. Zhao, X. Liu, J. Wang, S. Zhang, Carbohydr. Polym. 2013, 94, 723. J. Zhang, H. Zhang, J. Wu, J. Zhang, J. He, J. Xiang, Phys. Chem. Chem. Phys. 2010, 12, 1941. B. D. Rabideau, A. Agarwal, A. E. Ismail, J. Phys. Chem. B 2013, 117, 3469. Y. Zhao, X. Liu, J. Wang, S. Zhang, Chem. Phys. Chem. 2012, 13, 3126. L. Vanoye, M. Fanselow, J. D. Holbrey, M. P. Atkins, K. R. Seddon, Green Chem. 2009, 11, 390. T. Darden, D. York, L. Pedersen, J. Chem. Phys. 1993, 98, 10089. P. Sherwood, A. H. de Vries, M. F. Guest, G. Schreckenbach, C. R. A. Catlow, S. A. French, A. A. Sokol, S. T. Bromley, W. Thiel, A. J. Turner, S. Billeter, F. Terstegen, S. Thiel, J. Kendrick, S. C. Rogers, J. Casci, M. Watson, F. King, E. Karlsen, M. Sjøvoll, A. Fahmi, A. Sch€afer, C. Lennartz, J. Mol. Struc.-Theochem 2003, 632, 1. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. € Farkas, J. B. Foresman, J. V. Dannenberg, S. Dapprich, A. D. Daniels, O. Ortiz, J. Cioslowski, D. J. Fox, Gaussian 09, Revision A.1, Gaussian: Wallingford CT, 2009. R. Ahlrichs, M. B€ar, M. H€aser, H. Horn, C. K€ olmel, Chem. Phys. Lett. 1989, 162, 165. A. D. Becke, Phys. Rev. A 1988, 38, 3098. J. Slater, Phys. Rev. 1951, 81, 385. S. H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 1980, 58, 1200. J. Perdew, Phys. Rev. B 1986, 33, 8822. Y. Zhao, B. J. Lynch, D. G. Truhlar, J. Phys. Chem. A 2004, 108, 2715. W. Smith, T. R. Forester, J. Mol. Graphics 1996, 14, 136. V. Spiwok, B. Kralova, I. Tvaroska, Carbohyd. Res. 2010, 345, 530. D. Bakowies, W. Thiel, J. Phys. Chem. 1996, 100, 10580.

WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

[81] P. Sherwood, A. H. de Vries, S. J. Collins, S. P. Greatbanks, N. A. Burton, M. A. Vincent, I. H. Hillier, Faraday Discuss. 1997, 106, 79. [82] W. Crous, M. J. Field, K. J. Naidoo, J. Chem. Theory Comput. 2014. [83] S. R. Billeter, A. J. Turner, W. Thiel, Phys. Chem. Chem. Phys. 2000, 2, 2177. [84] J. K€astner, J. M. Carr, T. W. Keal, W. Thiel, A. Wander, P. Sherwood, J. Phys. Chem. A 2009, 113, 11856. [85] E. D. Glendening, A. E. Reed, J. E. Carpenter, F. Weinhold, NBO Version 3.1. [86] T. C. F. Gomes, M. S. Skaf, J. Comput. Chem. 2012, 33, 1338. [87] S. Metz, W. Thiel, J. Am. Chem. Soc. 2009, 131, 14885. [88] J. T. Ham, D. G. Williams, Acta Cryst. B 1970, 26, 1373. [89] A. Rencurosi, J. R€ ohrling, J. Pauli, A. Potthast, C. J€ager, S. P erez, P. Kosma, A. Imberty, Angew. Chem. Int. Ed. 2002, 41, 4277.

[90] X. Biarn es, A. Arde`vol, A. Planas, C. Rovira, A. Laio, M. Parrinello, J. Am. Chem. Soc. 2007, 129, 10686. [91] F. A. Momany, M. Appell, G. Strati, J. L. Willett, Carbohyd. Res. 2004, 339, 553. [92] J.-H. Q. Pinto, S. Kaliaguine, AIChE J. 1991, 37, 905. [93] Z.-E. Dadach, S. Kaliaguine, Can. J. Chem. Eng. 1993, 71, 880. [94] Y. B. Tewari, R. N. Goldberg, J. Biol. Chem. 1989, 264, 3966. [95] C. Loerbroks, E. Boulanger, W. Thiel, Chem.-Eur. J. 2015, DOI:10.1002/ chem.201405507.

Received: 23 December 2014 Revised: 21 February 2015 Accepted: 25 February 2015 Published online on 21 March 2015

Journal of Computational Chemistry 2015, 36, 1114–1123

1123

Copyright of Journal of Computational Chemistry is the property of John Wiley & Sons, Inc. and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.