12936
J. Phys. Chem. A 2010, 114, 12936–12944
Free Energy Landscape for Glucose Condensation Reactions Dajiang Liu,† Mark R. Nimlos,‡ David K. Johnson,§ Michael E. Himmel,§ and Xianghong Qian*,† Department of Mechanical Engineering, Colorado State UniVersity, Fort Collins, Colorado 80523, United States, and National Bioenergy Center and Bioscience Center, National Renewable Energy Laboratory, Golden, Colorado 80401, United States ReceiVed: August 18, 2010; ReVised Manuscript ReceiVed: October 27, 2010
Ab initio molecular dynamics and metadynamics simulations were used to determine the free energy surfaces (FES) for the acid catalyzed β-D-glucose condensation reaction. Protonation of C1-OH on the β-D-glucose, breakage of the C1-O1 bond, and the formation of C1 carbocation is the rate-limiting step. The effects of solvent on the reaction were investigated by determining the FES both in the absence and presence of solvent water. It was found that water played a critical role in these reactions. The reaction barrier for the protoncatalyzed glucose condensation reaction is solvent induced because of proton’s high affinity for water. During these simulations, β-D-glucose conversion to R-D-glucose process via the C1 carbocation was also observed. The associated free energy change and activation barrier for this reaction were determined. I. Introduction Lignocellulosic biomass represents an abundant renewable resource for the production of biobased products and biofuels. Cellulosic biomass is mainly composed of hemicelluloses (∼15-32%), cellulose (∼30-50%), and lignin (∼15-25%). Hemicelluloses, mostly xylan, are natural polymers of β-Dxylose and other minor sugars, whereas cellulose is made of β-D-glucose. Lignin is a polymer composed of nonfermentable phenylpropane monomer units. Both xylose and glucose sugar monomers are connected via the β-1,4 glucosidic linkage to form xylan and cellulose, respectively. During the biochemical conversion of biomass to biofuels, biomass is first hydrolyzed to monomer sugars (mostly xylose and glucose) by chemical and/or biological catalysts in the liquid phase. These monomer sugars are then converted to ethanol or butanol via fermentation or converted directly to liquid alkanes via catalytic liquid processes.1 Glucose is the most abundant monosaccharide from biomass; therefore, achieving high glucose yield is critical for the cost-effective production of biofuels. Thermochemical pretreatment of biomass opens up the biomass structure and has long been recognized as a critical step to produce cellulose with acceptable enzymatic digestibility.2 Dilute acid (e.g., ∼0.5-3.0% sulfuric acid by weight) is one of the most common and cost-effective agent used in pretreatment.3-14 Typically, dilute acid pretreatment of biomass is carried out at an elevated temperature of 430-500 K. During this process, hemicelluloses are hydrolyzed to monosaccharides, mostly the five carbon sugar xylose. In addition, a small amount of glucose is also released. The remaining cellulosic biomass substrate can be converted to glucose via enzymatic hydrolysis. Alternatively, such pretreated biomass can be processed in solvents, such as ionic liquid with or without added mineral acid to solubilize and further depolymerize cellulose to monomer glucose sugars.15 * Corresponding author. E-mail:
[email protected]. † Colorado State University. ‡ National Bioenergy Center, National Renewable Energy Laboratory. § Bioscience Center, National Renewable Energy Laboratory.
Figure 1. The structures of monomer β-D-glucose (left) and the β-1,4 linked glucose dimer cellobiose (right). Cellobiose is the basic repeating unit of cellulose as well as one of the glucose reversion products.
Depending on the severity (time, temperature, and acidity) of the acid pretreatment, some xylose and glucose molecules undergo an undesirable degradation process that lowers the biomass conversion efficiency. 2-Furaldehyde (furfural)16-19 and 5-(hydroxymethyl)-2-furalde (HMF)17,18,20-23 are major degradation products from xylose and glucose, respectively, in an acidic environment. In addition, at high sugar concentrations, these xylose and glucose molecules could react with each other to form various disaccharides.24 These condensation reactions are termed reVersion reactions. Reversion reactions have been recognized as an important cause of limited sugar yields at high biomass solid loading (∼30%) and, thus, higher final sugar concentrations. Cellobiose (shown in Figure 1) is a β-1,4 linked disaccharide of glucose, one of the many reversion products. Evidence from experiment and pilot scale dilute acid pretreatment of biomass24 indicates that some of these reversion products cannot be easily hydrolyzed and are toxic to the microorganisms for fermentation. Recent work by Pilath and co-workers24 has shown that glucose reversion reactions in dilute acid solutions at elevated temperatures produce multiple disaccharides. These products include trehalose (1,1-R-D-glucopyranosyl-R-D-glucopyranoside), neotrehalose (1,1R-D-glucopyranosyl-β-D-glucopyranoside), isotrehalose (1,1-β-Dglucopyranosyl-β-D-glucopyranoside), kojibiose (1,2-O-R-D-glucopyranosyl-D-glucose), sophorose (1,2-O-β-D-glucopyranosyl-Dglucose),nigerose(1,3-O-R-D-glucopyranosyl-D-glucose),laminaribiose
10.1021/jp1078407 2010 American Chemical Society Published on Web 11/18/2010
Free Energy of Glucose Condensation Reactions (1,2-O-β-D-glucopyranosyl-D-glucose), maltose (1,4-O-R-D-glucopyranosyl-D-glucose), cellobiose (1,4-O-β-D-glucopyranosylD-glucose), isomaltose (1,6-O-R-D-glucopyranosyl-D-glucose), and gentiobiose (1,6-O-β-D-glucopyranosyl-D-glucose). All of these disaccharides are linked via the C1 of the first glucose to various positions of the second glucose molecule. In addition, the estimated reaction barriers for the formation of these disaccharides are all between ∼120 and 140 kJ mol-1 (∼28 and 34 kcal mol-1), very similar to each other. This indicates that the formation of C1 carbocation is a critical step for these reversion reactions. Protonation of the C1-OH and the subsequent breaking of the C1-O1 bond are considered the ratelimiting steps for these reactions.25 Our recent work on acid-catalyzed sugar reactions indicates that water and water structure play an important role in these reactions at elevated temperatures.18,25 High temperature water no longer acts as an inert medium but is an active participant in many chemical reactions. Water can not only compete for protons with the sugar molecules, but it can also donate and extract a proton from the reaction intermediates. Moreover, it was found that the barriers for xylooligomer hydrolysis and for xylose condensation reactions are solvent-induced due primarily to the high proton affinity of the water molecules. Protonation of the hydroxyl groups on the sugar ring is the rate-limiting step during these reactions and contributes significantly to the observed reaction barrier. Importantly, proton migration from the bulk solution to the neighborhood of the sugar molecules via the so-called “partial proton desolvation process” also contributes substantially to this barrier.25 Even though β-Dglucose and β-D-xylose have very similar structures, earlier experimental evidence suggests that xylose degradation to furfural is more much rapid than glucose degradation to HMF, and with much higher yields.18 Since glucose and xylose behave somewhat differently, the effects of solvent on glucose reactions were investigated in this study. Car-Parrinello based ab initio molecule dynamics simulations (CPMD)26 coupled with metadynamics (MTD)27 were used to investigate the effects of water on glucose reversion reactions. Free energy surface (FES) reconstructed from our CPMD-MTD simulations indicate that the barrier for the glucose reversion reaction is also solvent induced. In addition to glucose reversion products, the formation of R-D-glucose was also observed via the mutarotation process in acidic aqueous solutions. II. Computational Methods and Procedure The metadynamics (MTD) method was developed by Parrinello and co-workers27,28 for investigating chemical reactions and exploring associated free energy surfaces (FES). MTD is designed to enhance the probabilities of the energy-barrier crossing events during a chemical reaction or process.29,30 When the energy barrier is more than a few times higher than thermal energy (kBT), conventional simulation techniques are unable to sample the configurational space that is important to determine the free energy surface. The MTD method is based on the ideas of extended Lagrangian27,31-33 and coarse-grained non-Markovian dynamics,27 which allows for very efficient exploration of the FES of the reactive system. This method assumes that several collective coordinates that distinguish reactants from products are able to characterize the reaction process. These collective coordinates (e.g., distances between atoms and coordination numbers) must include the relevant modes that cannot be sampled within the typical time scale of the ab initio MD simulations. The details of the method used for investigating FES of sugar reactions can be found in our earlier work.25
J. Phys. Chem. A, Vol. 114, No. 49, 2010 12937 SCHEME 1: Protonation of C1-OH on Glucose and the Associated Two CVs
The selection of the collective variables (CVs) for the protonation of the C1-OH in glucose is shown in Scheme 1. The first CV (CV1) is the coordination number (CN) of C1 with respect to O1. The second CV (CV2) is the CN of O1 with respect to all the four H atoms (three from H3O+ and one from C1-OH) assuming H atoms are indistinguishable. The equation of CN is given by
CN(i, j) )
1 - (dij /d0)p 1 - (dij /d0)q
(1)
where dij is the distance between atoms i and j, d0 is the cutoff distance, and p and q are high-power integers used to distinguish between the coordinated and noncoordinated states. For CV1 and CV2, the values for d0 are chosen to be 2.0 and 1.5 Å, respectively. The values p ) 6 and q ) 12 were routinely chosen for both CVs. The dynamics of the CVs are controlled by the force constant k and fictitious mass m. In our simulations, k ) 2.0 au and m ) 50 amu were used for both CVs as in our earlier calculations.25 MTD sampling technique is based on the filling of the reactant and product wells with bias potentials to facilitate barrier-crossing. Once the reactant well is filled close enough to the reaction barrier, the system crosses the barrier and moves to the product well. When the product well is also filled, the FES becomes flat and the system can move back and forth without any barrier. Thus, FES can be reconstructed on the basis of the amount of bias potentials added to each well. The bias potential V(s) used in our current CPMD-MTD simulations is a commonly used Gaussian functional. The height of the Gaussian is usually a few percent of the reaction barrier.34,35 The height and width of the Gaussian bias potentials were chosen to be 0.01 and 0.05 au initially for all the simulations. When the first barrier-crossing was observed, the Gaussian height was reduced to 0.005 au and was fixed for the rest of simulations. The bias potentials were added whenever the displacements in the CVs were larger than 1.5 times the width, but no shorter than 100 MD steps. Studies have shown that this choice of parameters is efficient in filling energy wells within an uncertainty of 1-2 kcal/mol.29,35,36 All our calculations were carried out using the CPMD software package.26 The Becke, Lee, Yang, and Parr (BLYP) functional was used to describe valence and semicore electrons.37,38 The potentials exerted by the core electrons were approximated by the Goedecker pseudopotentials.39 An energy cutoff of 70 Ry was used to for the plane-wave basis sets. To effectively separate the motions of electrons from those of slow-moving nuclei, a fictitious mass of 800 amu and a time step of 0.125 fs were used. The MD simulations were carried out under NVT at 300 K with a Nose-Hoover chain thermostat.40 In the gas phase, the simulation box has a dimension of 15 × 15 × 15 Å and was decoupled from its image by using the Hockney’s
12938
J. Phys. Chem. A, Vol. 114, No. 49, 2010
Liu et al.
Figure 2. The dynamics of the two reaction coordinates CV1 (C1-O1 in black) and CV2 (O1-H in red) during the course of CPMD-MTD simulations for protonation of C1-OH on β-D-glucose and the subsequent breakage of C1-O1 bond in the gas phase.
method with an extra 4 Å added to each dimension of the simulation box.41 In solution, the periodic boundary condition was used for a box containing one glucose, one H3O+, one Clcounterion, and 76 H2O molecules. The simulation box size has a dimension of 14.4 × 14.4 × 14.4 Å and a corresponding density of 0.92 g/cm3. Ewald summation was used to integrate the long-range electrostatic interaction energies. The FES generated from CPMD and MTD simulations also yields information on the reaction barriers and the transition states of chemical reactions. The transition states and activation barriers thus obtained should be more reliable compared to the static quantum chemical methods, since the whole FES is sampled in CPMD-MTD simulations whereas only certain points are typically calculated in static calculations. III. Results and Discussions 1. Formation of β-D-Glucose C1+ in the Gas Phase. In order to fully understand the effects of solvent on glucose reversion reaction, the reaction was initially investigated in the absence of the solvent in the gas phase. As mentioned earlier, protonation of C1-OH on β-D-glucose and the subsequent breakage of the C1-O1 bond is a critical step in the formation of glucose disaccharides via condensation reactions. The CPMD-MTD simulations were carried out at 300 K to reconstruct the FES based on the bias potentials added to the reactant and product wells assuming that FES is not temperature dependent. Figure 2 shows the fluctuations of the two collective variables as a function of simulation time in MTD steps. At the beginning of the simulations, a proton is attached to the C1-OH group. During the course of the CPMD-MTD simulations, CV1 representing the C1-O1 bond is seen fluctuating between 0 (broken) and 1 (bonded) states, whereas CV2 representing the O1-H bond fluctuates between 1 (proton on water) and 2 (proton transferred to the OH group) states. At around 1300 MTD steps, the proton was transferred back to the water molecule and stays there for a period of time during which the C1 and O1 atoms are bonded. After sufficient sampling of the CV configuration space as seen in Figure 3 with about 2550 MTD steps, the FES was reconstructed on the basis of the amount of bias potentials added along the sampling pathways. Figures 4 and 5 are the 2-D free energy contour plot and 3-D FES, respectively, for protonation of C1-OH on β-D-
Figure 3. The trajectories of CV1 (C1-O1) and CV2 (O1-H) during the sampling of the reaction and production wells for protonation of C1-OH on β-D-glucose, breakage of C1-O1 bond, and formation of C1 carbocation in the gas phase.
Figure 4. The 2-D free energy contour plot for protonation of C1-OH and breakage of the C1-O1 bond in the gas phase. The structures of the reactant (unprotonated glucose), intermediate (protonated glucose), and the final product (carbocation) are shown in the insets. The energy scale is in kcal/mol.
glucose and the breakage of C1-O1 bond. It can be seen that the unprotonated β-D-glucose is situated on FES at a shallow minimum (CV1 ≈ 0.8 and CV2 ≈ 1). A small barrier
