- Email: [email protected]
Potential of mean force for separation of the repeating units in cellulose and hemicellulose
Carbohydrate Research 346 (2011) 867–871
Contents lists available at ScienceDirect
Carbohydrate Research journal homepage: www.elsevier.com/locate/carres
Note
Potential of mean force for separation of the repeating units in cellulose and hemicellulose Suma Peri, M. Nazmul Karim, Rajesh Khare ⇑ Department of Chemical Engineering, Texas Tech University, Box 43121, Lubbock, TX 79409-3121, USA
a r t i c l e
i n f o
Article history: Received 3 June 2010 Received in revised form 5 January 2011 Accepted 7 January 2011 Available online 15 January 2011 Keywords: Cellobiose Xylobiose Disaccharides Potential of mean force Umbrella sampling Molecular dynamics simulations
a b s t r a c t As a first step toward understanding the energetics of removal of cello-oligomers from the cellulose surface, we have performed umbrella sampling calculations to determine the free energy required for separation of repeating units of cellulose and hemicellulose from each other. Molecular dynamics (MD) simulations were performed for both the stacked and non-stacked arrangements of the cellobiose pair system and the xylobiose pair system. These stacked and non-stacked arrangements were taken as representative systems for the crystalline and amorphous domains in cellulose and hemicellulose. In addition, similar calculations were also carried out to determine the energetics involved in the separation of the cellobiose–xylobiose molecule pair in the non-stacked arrangement. The potential of mean force profiles exhibit a single minimum in all cases and are qualitatively similar. Our results show that the location of the minimum as well as the depth of the well can be correlated with the size of the disaccharide molecules. Ó 2011 Elsevier Ltd. All rights reserved.
Lignocellulosic biomass serves as a potential renewable source of energy in the production of biofuels. The main constituents of lignocellulosic biomass are cellulose, hemicellulose, and lignin. Cellulose is the major constituent of biomass (35–50% by dry weight), and is formed by repeating units of cellobiose, a disaccharide molecule with a b-(1?4)-glucosidic linkage. Hemicellulose mainly consists of five carbon sugars; it constitutes around 20– 35% (dry weight) of biomass. Lignin is a hydrophobic network of phenyl propanoid unit, which binds cellulose and hemicellulose chains together.1,2 An important step in the production of ethanol from biomass consists of depolymerization of the cellulose and hemicellulose components by cellulase and xylanase enzymes, respectively. The action of these enzymes on their substrates results in the formation of several oligosaccharides which are the reaction intermediates. These intermediates are eventually broken down to the product molecules, namely, glucose and xylose.3 Among these cello-oligosaccharides, only the short ones (degree of polymerization less than seven) are soluble; their solubility decreases with an increase in the chain length. The longer cello-oligosaccharide intermediates thus accumulate during the reaction and are known to adsorb back onto the cellulosic surfaces.3,4 This phenomenon hinders the release of these chain fragments into the solvent, and thus inhibits further action of the enzyme thereby effectively slowing down the rate of the hydrolysis reaction.5,6 Thus, the knowledge of the interactions within and between the
⇑ Corresponding author. Tel.: +1 806 742 0449; fax: +1 806 742 3552. E-mail address: [email protected] (R. Khare). 0008-6215/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.carres.2011.01.008
cellulose and hemicellulose derived oligosaccharides will be crucial for developing strategies for increasing the rate of cellulose hydrolysis reaction. As a first step toward this goal, in this note, we have quantified these interactions for cellobiose and xylobiose molecules using molecular simulations. Molecular simulations provide the ability for direct quantification of the effects of specific molecular interactions. A number of previous simulations have focused on the sugar conformations, intrachain hydrogen bonding and solvent structuring around the mono- and disaccharide molecules.7–15 Simulations have also been reported to characterize the interactions of tetramers of glucose, xylose, and mannose with the cellulose crystal.16 These results showed that the binding affinity as measured by the potential energy was the largest for cellotetraose followed by mannotetraose and xylotetraose, respectively. The aggregation behavior of cello-oligosaccharides in aqueous solutions was also studied by MD simulations;17 these simulations showed that two strands of cellohexaose remained aggregated in a helical conformation over a time scale of 1 ns, whereas two strands of cellotetraose separated from each other over the same time scale. Umbrella sampling calculations were performed to determine the free energy for a glucose and benzene system which was considered to be a prototype for an enzyme–substrate system.18 Similarly, potential of mean force (PMF) profile for the process of adsorption of D-glucose on a continuum graphite surface was also determined.19 In both cases, the PMF for the binding of the molecules in aqueous solutions was shown to have an oscillatory behavior where the peak positions corresponded to the water layers. The PMF for the separation of cello-oligomers with a degree
868
S. Peri et al. / Carbohydrate Research 346 (2011) 867–871
of polymerization ranging from 1 to 4 (i.e., glucose, cellobiose, cellotriose and cellotetraose) was reported in a recent study.20 For both the cases studied in that work (a configuration in which the molecule planes face each other and another in which the molecules were placed adjacent to each other in a co-planar fashion), the PMF showed weak oscillatory behavior for cellobiose and the two longer cello-oligosaccharides. The focus of this paper is on the determination of the free energy involved in the separation of the disaccharides (which are the shortest oligomers of cellulose and hemicellulose polymers). Simulation systems containing cellobiose and xylobiose pairs in water were prepared to study the free energy profiles for two types of arrangements: (1) Non-stacked arrangement (shown in Fig. 1a and c), where the disaccharide molecules were allowed to have complete orientational freedom with respect to the other molecule; this is considered to represent the amorphous regions of the substrate, and (2) Stacked arrangement (shown in Fig. 1d and e), where the disaccharide molecules were maintained parallel to each other in the system at all times during the simulation; this is considered to represent the crystalline regions in the substrate. As is the case for crystalline cellulose, both sets of disaccharide molecules were maintained in an arrangement where they are off-set from the ‘face-to-face’ configuration (see Fig. 1d and e). We note that the cellobiose pair system was also studied in Ref. , however, the arrangements studied in that work were different from the ones studied here. Specifically, the two arrangements studied in that work consisted of the stacked arrangement in which the cellobiose molecules were in the face-to-face (not off20
set as in a crystal) configuration and a co-planar arrangement of the cellobiose molecules. In addition, the interactions between the repeating units of cellulose and hemicellulose were studied by considering the cellobiose–xylobiose-water system (Fig. 1b). General Amber force field (GAFF) was used to describe the carbohydrate models in this study.21 Further details of the procedure used for the calculation of the PMF curves are provided in the Methods section. The PMF profiles for the process of separation of a cellobiose pair, a pair of cellobiose–xylobiose molecules, and a xylobiose pair in the non-stacked arrangement in the presence of water at room temperature are shown in Figure 2a. The constructed PMFs are the averages of three sets of umbrella sampling simulations, with the uncertainties being determined as the standard deviations (shown for only one point in the figure for the sake of clarity). All three profiles are qualitatively similar and exhibit a single shallow minimum; furthermore, as the separation distance between the molecules increases, the profiles flatten around a distance of 7.5 Å, indicating the absence of any significant interactions between the molecules at longer distances. A comparison of the positions of the minima and the well depth for the three systems indicates a monotonic relationship that can be correlated with the size of the molecules, where we note that xylobiose is smaller in size than cellobiose. The location of the minimum for the cellobiose pair is 4.75 Å while for the cellobiose–xylobiose separation, this location decreased to 4.5 Å; the value further decreased to 4.2 Å for the xylobiose pair. The separation distance at which the minimum occurs is slightly larger than the Lennard–Jones (LJ) diameter for the atoms (e.g., the LJ diameter for carbon and oxygen atom pair is 3.19 Å), we attribute this to the axial projection of the protons on both sides of the solute molecules. The well depth for the cellobiose pair is found to be 1.84 kcal/mol, that is nearly three times the value of RT. This value gets progressively smaller as one
Figure 1. Schematic representation of the disaccharide molecules in the non-stacked arrangement: (a) cellobiose pair, (b) xylobiose–cellobiose pair and (c) xylobiose pair; and in the stacked arrangement: (d) cellobiose pair and (e) xylobiose pair. The carbon atoms are shown in gray, oxygen in black, and hydrogen in white colors.
S. Peri et al. / Carbohydrate Research 346 (2011) 867–871
869
Figure 2. Potential of mean force (PMF) profiles for (a) cellobiose pair, xylobiose pair and cellobiose–xylobiose in the non-stacked arrangement and (b) cellobiose pair and xylobiose pair in the stacked arrangement. Profiles represent the average of three sets of simulations and lines are shown as a guide to eye.
moves to the cellobiose–xylobiose (1.2 kcal/mol) and the xylobiose pair (1.0 kcal/mol) systems. The well depth of the free energy profile for the separation of the xylobiose pair is approximately half of that for the separation of the cellobiose pair. In these simulations with non-stacked arrangement of disaccharide molecules, the molecules were allowed complete orientational freedom. It was observed that these molecules tend to align in such a way that the distance between the hydrophilic end of one molecule and the hydrophilic face of the other molecule is minimized (similar to the configuration shown in Fig. 1a), while still maintaining the distance constraint on their centers of masses. Hence, the favorable orientation is the configuration in which the hydrophilic groups of each molecule face toward one another. In such a configuration, as seen from Figure 2a, the free energy profile for the separation of these disaccharide molecules exhibits only one minimum, similar to that observed for the separation of two hydrophilic solute molecules in water.22
Umbrella sampling simulations were also carried out to determine the PMF for the process of separation of disaccharides that are arranged in the stacked arrangement (Fig. 1d and e) which is representative of the structural arrangement in crystalline cellulosic substrates. Figure 2b presents the PMF profiles for the separation of a pair of cellobiose molecules and a pair of xylobiose molecules that are in the stacked configuration in aqueous solutions. The profiles show the average from three sets of umbrella sampling simulations. As can be seen from the figure, the PMF profiles show only one minimum at a separation distance of about 4.25 Å with a well depth value of around 3.2 kcal/mol for xylobiose pair and 4.0 kcal/mol for cellobiose pair; the profiles flatten out at distances greater than 6 Å. A schematic representation of the configuration of the two cellobiose molecules that are in the stacked arrangement is shown in Figure 3. It can be seen that some of the oxygen atoms of the hydroxyl groups which are capable of forming hydrogen bonds with water molecules are found on the in-
870
S. Peri et al. / Carbohydrate Research 346 (2011) 867–871
ner surface. Hence, even in the stacked arrangement, some of the hydrophilic ends of the cellobiose molecules face each other. As a result, the PMF profile is qualitatively similar to that for hydrophilic solutes, as was also the case for the non-stacked arrangement. However, the well depth for the stacked arrangement is much higher than that for the non-stacked case. This can be explained by considering that while the configuration in the free energy minimum state for the stacked configuration (at 4.25 Å) is energetically favorable (since it is similar to that in the crystal), the stacked configuration at larger separations is of a higher free energy due to the lack of orientational freedom; the net result is a larger well depth. The position of the first minimum obtained by us (4.25 Å) matches quite well with the first minimum location (4.6 Å) recently reported in the PMF profile for the stacked configuration of cellobiose;20 also, our well depth (4 kcal/mol) is close to the value (3.3 kcal/mol) reported in that work. We note that the stacked arrangement of the cellobiose pair in our system (off-set from face-to-face) is different from that in the literature study20 where the cellobiose molecules interact with each other in a face-to-face configuration during the simulations. Finally, we also note that our result for the stacked configuration should be considered to be an upper bound on the free energy required for separation of cellobiose in the crystalline domain since although the cellobiose molecule under consideration will be in the ‘stacked’ configuration when it is near a crystal surface, it will gain orientational freedom as it moves away from the surface. The calculations reported here have focused on the smallest units of the cellulose and hemicellulose substrates. For the enzymatic hydrolysis process, the main goal is to elucidate the energetics involved in the separation of cello-oligosaccharides from the crystalline cellulose surface; further calculations are currently underway to achieve this. 1. Methods 1.1. System preparation and molecular dynamics simulation details The structures of the repeating unit of cellulose and hemicellulose were created with the user interface program xleap in AMBER
Figure 3. Schematic representation of two cellobiose molecules in the stacked arrangement. The oval shows the atoms present in the area encompassing the planes of the two cellobiose molecules. Carbon, oxygen, and hydrogen atoms are shown in gray, black, and white colors, respectively.
8.0.23 The pair of disaccharide molecules was solvated with 2050 TIP3P water24 molecules using xleap in a cubic periodic box. The force field parameters for water and carbohydrates were taken from ff03 and GAFF force fields, respectively.21,24,25 All the simulations reported in this work were performed using AMBER 8.0 molecular dynamics simulation package.23 Conjugate gradient energy minimization was performed to relax the systems. A sequence of constant NVT (constant number of molecules, constant volume and constant temperature) and constant NPT (constant number of molecules, constant pressure and constant temperature) simulations employing the Berendsen weak coupling algorithm were carried out on the systems to equilibrate them at a temperature of 300 K and a pressure of 1 atm.26 These resulted in final density values that are close to 1 g/cc for all the five systems. In these simulations, the time-step for integrating the MD trajectories was set to 1 fs and the OH bond length for the solvent in the systems was constrained using the SHAKE approach.27 The non-bonded and electrostatic interactions were truncated at a cut-off distance of 10 Å. Particle Mesh Ewald (PME) method was employed to account for the long-range electrostatic interactions.28 1.2. Umbrella sampling simulations Straightforward MD simulations, in general, tend to avoid the high energy regions of the configuration space; this aspect leads to inadequate sampling of the phase space. Hence to adequately sample these high energy regions in free energy calculations, an additional biasing harmonic potential (umbrella potential) is employed which is later subtracted from the total free energy to give the net or unbiased free energy change. Umbrella sampling method as developed originally as well as its modifications have been extensively used to estimate the free energy difference between two states at equilibrium.18,20,29–32 In this work, the reaction coordinate for the umbrella sampling simulations for the systems described above (see Fig. 1) was defined as the center of mass separation distance between the molecule pair. In a series of simulations, this distance was varied to represent the process of separation of the disaccharide molecules; these simulations resulted in individual windows in the umbrella sampling procedure. A series of 48 overlapping umbrella sampling windows were considered which spanned the center of mass separation distance between the disaccharide pair from 15 Å to 3.25 Å in increments of 0.25 Å. The separation distance between the molecule centers of mass in each of these windows was maintained near the specified reaction coordinate value by applying a harmonic biasing potential with a force constant of 25 kcal/mol Å2 on the center of mass separation distance. A 500 ps long constant NVT simulation was then performed for equilibrating the system at each separation distance (umbrella sampling window). One set of 1 ns long and two sets of 500 ps long simulation for the cellobiose pair and three sets of 1 ns long production runs for cellobiose– xylobiose and xylobiose pair (initial velocities were reassigned in each case) were then carried out and the trajectories were written every 100 steps. PMF for each set was computed from the generated trajectory data for all of the windows using the Weighted Histogram Analysis Method (WHAM).33 The same procedure was employed to determine the PMF profiles in the stacked arrangement with the exception that, in this case, additional potentials were required to maintain the stacked configuration. Thus, in addition to the umbrella potential on the center of mass distance, harmonic restraining potentials were also applied on the distances between the chosen carbon atoms on the cellobiose molecules (atoms C1, C3, C5, C7, C9, and C11) and on the xylobiose molecules (atoms C1, C3, C5, C7, C9, and C10). These potentials ensured parallel arrangement of the molecules. A 500 ps long constant NVT simulation was performed for equilibrat-
S. Peri et al. / Carbohydrate Research 346 (2011) 867–871
ing the system at each umbrella sampling window. One set of 1 ns long and two sets of 500 ps long simulations for each of the umbrella sampling windows for the cellobiose pair and three sets of 1 ns long simulations for the xylobiose pair were performed. Acknowledgment We gratefully acknowledge the partial financial support of this work by the National Science Foundation (Grant number: CBET 0854463). References 1. Lynd, L. R.; Weimer, P. J.; van Zyl, W. H.; Pretorius, I. S. Microbiol. Mol. Biol. Rev. 2002, 66, 506–577. 2. Zhang, Y. H. P.; Lynd, L. R. Biotechnol. Bioeng. 2004, 88, 797–824. 3. Wyman, C. E.; Dale, B. E.; Elander, R. T.; Holtzapple, M.; Ladisch, M. R.; Lee, Y. Y. Bioresour. Technol. 2005, 96, 1959–1966. 4. Peri, S.; Karra, S.; Lee, Y. Y.; Karim, M. N. Biotechnol. Prog. 2007, 23, 626–637. 5. KlemanLeyer, K. M.; SiikaAho, M.; Teeri, T. T.; Kirk, T. K. Appl. Environ. Microbiol. 1996, 62, 2883–2887. 6. 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–152. 7. Brady, J. W. Carbohydr. Res. 1987, 165, 306–312. 8. French, A. D.; Johnson, G. P. Cellulose 2004, 11, 449–462. 9. French, A. D.; Johnson, G. P. Can. J. Chem./Revue Can. Chim. 2006, 84, 603–612. 10. French, A. D.; Johnson, G. P. Cellulose 2009, 16, 959–973. 11. Hardy, B. J.; Sarko, A. J. Comput. Chem. 1993, 14, 848–857. 12. Liu, Q.; Brady, J. W. J. Phys. Chem. B 1997, 101, 1317–1321. 13. Shen, T. Y.; Langan, P.; French, A. D.; Johnson, G. P.; Gnanakaran, S. J. Am. Chem. Soc. 2009, 131, 14786–14794. 14. Stortz, C. A.; French, A. D. Mol. Simul. 2008, 34, 373–389.
871
15. Stortz, C. A.; Johnson, G. P.; French, A. D.; Csonka, G. I. Carbohydr. Res. 2009, 344, 2217–2228. 16. Leeflang, B. R.; van Kuik, J. A.; Kroon-Bratenburg, L. M. J. In NMR Spectroscopy and Computer Modeling of Carbohydrates: Recent Advances. ACS Symposium Series; 2006; Vol. 930, pp 133–155. 17. Umemura, M.; Yuguchi, Y.; Hirotsu, T. J. Phys. Chem. A 2004, 108, 7063–7070. 18. Palma, R.; Himmel, M. E.; Brady, J. W. J. Phys. Chem. B 2000, 104, 7228–7234. 19. Glennon, T. M.; Merz, K. M. J. Mol. Struct. (Theochem) 1997, 395–396, 157–171. 20. Bergenstrahle, M.; Wohlert, J.; Himmel, M. E.; Brady, J. W. Carbohydr. Res. 2010, 345, 2060–2066. 21. Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem 2004, 25, 1157–1174. 22. Durell, S. R.; Brooks, B. R.; Ben-Naim, A. J. Phys. Chem. 1994, 98, 2198–2202. 23. Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Merz, K. E.; 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, 2004. 24. Mark, P.; Nilsson, L. J. Phys. Chem. A 2001, 105, 9954–9960. 25. Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G. M.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J. M.; Kollman, P. A. J. Comput. Chem. 2003, 24, 1999–2012. 26. Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. 27. Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327– 341. 28. Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089–10092. 29. Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187–199. 30. Trzesniak, D.; Kunz, A. P. E.; van Gunsteren, W. F. Chemphyschem 2007, 8, 162– 169. 31. Chipot, C.; Jaffe, R.; Maigret, B.; Pearlman, D. A.; Kollman, P. A. J. Am. Chem. Soc. 1996, 118, 11217–11224. 32. Hansen, H. S.; Hunenberger, P. H. J. Comput. Chem. 2010, 31, 1–23. 33. Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, R. M. J. Comput. Chem. 1992, 13, 1011–1021.
Contents lists available at ScienceDirect
Carbohydrate Research journal homepage: www.elsevier.com/locate/carres
Note
Potential of mean force for separation of the repeating units in cellulose and hemicellulose Suma Peri, M. Nazmul Karim, Rajesh Khare ⇑ Department of Chemical Engineering, Texas Tech University, Box 43121, Lubbock, TX 79409-3121, USA
a r t i c l e
i n f o
Article history: Received 3 June 2010 Received in revised form 5 January 2011 Accepted 7 January 2011 Available online 15 January 2011 Keywords: Cellobiose Xylobiose Disaccharides Potential of mean force Umbrella sampling Molecular dynamics simulations
a b s t r a c t As a first step toward understanding the energetics of removal of cello-oligomers from the cellulose surface, we have performed umbrella sampling calculations to determine the free energy required for separation of repeating units of cellulose and hemicellulose from each other. Molecular dynamics (MD) simulations were performed for both the stacked and non-stacked arrangements of the cellobiose pair system and the xylobiose pair system. These stacked and non-stacked arrangements were taken as representative systems for the crystalline and amorphous domains in cellulose and hemicellulose. In addition, similar calculations were also carried out to determine the energetics involved in the separation of the cellobiose–xylobiose molecule pair in the non-stacked arrangement. The potential of mean force profiles exhibit a single minimum in all cases and are qualitatively similar. Our results show that the location of the minimum as well as the depth of the well can be correlated with the size of the disaccharide molecules. Ó 2011 Elsevier Ltd. All rights reserved.
Lignocellulosic biomass serves as a potential renewable source of energy in the production of biofuels. The main constituents of lignocellulosic biomass are cellulose, hemicellulose, and lignin. Cellulose is the major constituent of biomass (35–50% by dry weight), and is formed by repeating units of cellobiose, a disaccharide molecule with a b-(1?4)-glucosidic linkage. Hemicellulose mainly consists of five carbon sugars; it constitutes around 20– 35% (dry weight) of biomass. Lignin is a hydrophobic network of phenyl propanoid unit, which binds cellulose and hemicellulose chains together.1,2 An important step in the production of ethanol from biomass consists of depolymerization of the cellulose and hemicellulose components by cellulase and xylanase enzymes, respectively. The action of these enzymes on their substrates results in the formation of several oligosaccharides which are the reaction intermediates. These intermediates are eventually broken down to the product molecules, namely, glucose and xylose.3 Among these cello-oligosaccharides, only the short ones (degree of polymerization less than seven) are soluble; their solubility decreases with an increase in the chain length. The longer cello-oligosaccharide intermediates thus accumulate during the reaction and are known to adsorb back onto the cellulosic surfaces.3,4 This phenomenon hinders the release of these chain fragments into the solvent, and thus inhibits further action of the enzyme thereby effectively slowing down the rate of the hydrolysis reaction.5,6 Thus, the knowledge of the interactions within and between the
⇑ Corresponding author. Tel.: +1 806 742 0449; fax: +1 806 742 3552. E-mail address: [email protected] (R. Khare). 0008-6215/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.carres.2011.01.008
cellulose and hemicellulose derived oligosaccharides will be crucial for developing strategies for increasing the rate of cellulose hydrolysis reaction. As a first step toward this goal, in this note, we have quantified these interactions for cellobiose and xylobiose molecules using molecular simulations. Molecular simulations provide the ability for direct quantification of the effects of specific molecular interactions. A number of previous simulations have focused on the sugar conformations, intrachain hydrogen bonding and solvent structuring around the mono- and disaccharide molecules.7–15 Simulations have also been reported to characterize the interactions of tetramers of glucose, xylose, and mannose with the cellulose crystal.16 These results showed that the binding affinity as measured by the potential energy was the largest for cellotetraose followed by mannotetraose and xylotetraose, respectively. The aggregation behavior of cello-oligosaccharides in aqueous solutions was also studied by MD simulations;17 these simulations showed that two strands of cellohexaose remained aggregated in a helical conformation over a time scale of 1 ns, whereas two strands of cellotetraose separated from each other over the same time scale. Umbrella sampling calculations were performed to determine the free energy for a glucose and benzene system which was considered to be a prototype for an enzyme–substrate system.18 Similarly, potential of mean force (PMF) profile for the process of adsorption of D-glucose on a continuum graphite surface was also determined.19 In both cases, the PMF for the binding of the molecules in aqueous solutions was shown to have an oscillatory behavior where the peak positions corresponded to the water layers. The PMF for the separation of cello-oligomers with a degree
868
S. Peri et al. / Carbohydrate Research 346 (2011) 867–871
of polymerization ranging from 1 to 4 (i.e., glucose, cellobiose, cellotriose and cellotetraose) was reported in a recent study.20 For both the cases studied in that work (a configuration in which the molecule planes face each other and another in which the molecules were placed adjacent to each other in a co-planar fashion), the PMF showed weak oscillatory behavior for cellobiose and the two longer cello-oligosaccharides. The focus of this paper is on the determination of the free energy involved in the separation of the disaccharides (which are the shortest oligomers of cellulose and hemicellulose polymers). Simulation systems containing cellobiose and xylobiose pairs in water were prepared to study the free energy profiles for two types of arrangements: (1) Non-stacked arrangement (shown in Fig. 1a and c), where the disaccharide molecules were allowed to have complete orientational freedom with respect to the other molecule; this is considered to represent the amorphous regions of the substrate, and (2) Stacked arrangement (shown in Fig. 1d and e), where the disaccharide molecules were maintained parallel to each other in the system at all times during the simulation; this is considered to represent the crystalline regions in the substrate. As is the case for crystalline cellulose, both sets of disaccharide molecules were maintained in an arrangement where they are off-set from the ‘face-to-face’ configuration (see Fig. 1d and e). We note that the cellobiose pair system was also studied in Ref. , however, the arrangements studied in that work were different from the ones studied here. Specifically, the two arrangements studied in that work consisted of the stacked arrangement in which the cellobiose molecules were in the face-to-face (not off20
set as in a crystal) configuration and a co-planar arrangement of the cellobiose molecules. In addition, the interactions between the repeating units of cellulose and hemicellulose were studied by considering the cellobiose–xylobiose-water system (Fig. 1b). General Amber force field (GAFF) was used to describe the carbohydrate models in this study.21 Further details of the procedure used for the calculation of the PMF curves are provided in the Methods section. The PMF profiles for the process of separation of a cellobiose pair, a pair of cellobiose–xylobiose molecules, and a xylobiose pair in the non-stacked arrangement in the presence of water at room temperature are shown in Figure 2a. The constructed PMFs are the averages of three sets of umbrella sampling simulations, with the uncertainties being determined as the standard deviations (shown for only one point in the figure for the sake of clarity). All three profiles are qualitatively similar and exhibit a single shallow minimum; furthermore, as the separation distance between the molecules increases, the profiles flatten around a distance of 7.5 Å, indicating the absence of any significant interactions between the molecules at longer distances. A comparison of the positions of the minima and the well depth for the three systems indicates a monotonic relationship that can be correlated with the size of the molecules, where we note that xylobiose is smaller in size than cellobiose. The location of the minimum for the cellobiose pair is 4.75 Å while for the cellobiose–xylobiose separation, this location decreased to 4.5 Å; the value further decreased to 4.2 Å for the xylobiose pair. The separation distance at which the minimum occurs is slightly larger than the Lennard–Jones (LJ) diameter for the atoms (e.g., the LJ diameter for carbon and oxygen atom pair is 3.19 Å), we attribute this to the axial projection of the protons on both sides of the solute molecules. The well depth for the cellobiose pair is found to be 1.84 kcal/mol, that is nearly three times the value of RT. This value gets progressively smaller as one
Figure 1. Schematic representation of the disaccharide molecules in the non-stacked arrangement: (a) cellobiose pair, (b) xylobiose–cellobiose pair and (c) xylobiose pair; and in the stacked arrangement: (d) cellobiose pair and (e) xylobiose pair. The carbon atoms are shown in gray, oxygen in black, and hydrogen in white colors.
S. Peri et al. / Carbohydrate Research 346 (2011) 867–871
869
Figure 2. Potential of mean force (PMF) profiles for (a) cellobiose pair, xylobiose pair and cellobiose–xylobiose in the non-stacked arrangement and (b) cellobiose pair and xylobiose pair in the stacked arrangement. Profiles represent the average of three sets of simulations and lines are shown as a guide to eye.
moves to the cellobiose–xylobiose (1.2 kcal/mol) and the xylobiose pair (1.0 kcal/mol) systems. The well depth of the free energy profile for the separation of the xylobiose pair is approximately half of that for the separation of the cellobiose pair. In these simulations with non-stacked arrangement of disaccharide molecules, the molecules were allowed complete orientational freedom. It was observed that these molecules tend to align in such a way that the distance between the hydrophilic end of one molecule and the hydrophilic face of the other molecule is minimized (similar to the configuration shown in Fig. 1a), while still maintaining the distance constraint on their centers of masses. Hence, the favorable orientation is the configuration in which the hydrophilic groups of each molecule face toward one another. In such a configuration, as seen from Figure 2a, the free energy profile for the separation of these disaccharide molecules exhibits only one minimum, similar to that observed for the separation of two hydrophilic solute molecules in water.22
Umbrella sampling simulations were also carried out to determine the PMF for the process of separation of disaccharides that are arranged in the stacked arrangement (Fig. 1d and e) which is representative of the structural arrangement in crystalline cellulosic substrates. Figure 2b presents the PMF profiles for the separation of a pair of cellobiose molecules and a pair of xylobiose molecules that are in the stacked configuration in aqueous solutions. The profiles show the average from three sets of umbrella sampling simulations. As can be seen from the figure, the PMF profiles show only one minimum at a separation distance of about 4.25 Å with a well depth value of around 3.2 kcal/mol for xylobiose pair and 4.0 kcal/mol for cellobiose pair; the profiles flatten out at distances greater than 6 Å. A schematic representation of the configuration of the two cellobiose molecules that are in the stacked arrangement is shown in Figure 3. It can be seen that some of the oxygen atoms of the hydroxyl groups which are capable of forming hydrogen bonds with water molecules are found on the in-
870
S. Peri et al. / Carbohydrate Research 346 (2011) 867–871
ner surface. Hence, even in the stacked arrangement, some of the hydrophilic ends of the cellobiose molecules face each other. As a result, the PMF profile is qualitatively similar to that for hydrophilic solutes, as was also the case for the non-stacked arrangement. However, the well depth for the stacked arrangement is much higher than that for the non-stacked case. This can be explained by considering that while the configuration in the free energy minimum state for the stacked configuration (at 4.25 Å) is energetically favorable (since it is similar to that in the crystal), the stacked configuration at larger separations is of a higher free energy due to the lack of orientational freedom; the net result is a larger well depth. The position of the first minimum obtained by us (4.25 Å) matches quite well with the first minimum location (4.6 Å) recently reported in the PMF profile for the stacked configuration of cellobiose;20 also, our well depth (4 kcal/mol) is close to the value (3.3 kcal/mol) reported in that work. We note that the stacked arrangement of the cellobiose pair in our system (off-set from face-to-face) is different from that in the literature study20 where the cellobiose molecules interact with each other in a face-to-face configuration during the simulations. Finally, we also note that our result for the stacked configuration should be considered to be an upper bound on the free energy required for separation of cellobiose in the crystalline domain since although the cellobiose molecule under consideration will be in the ‘stacked’ configuration when it is near a crystal surface, it will gain orientational freedom as it moves away from the surface. The calculations reported here have focused on the smallest units of the cellulose and hemicellulose substrates. For the enzymatic hydrolysis process, the main goal is to elucidate the energetics involved in the separation of cello-oligosaccharides from the crystalline cellulose surface; further calculations are currently underway to achieve this. 1. Methods 1.1. System preparation and molecular dynamics simulation details The structures of the repeating unit of cellulose and hemicellulose were created with the user interface program xleap in AMBER
Figure 3. Schematic representation of two cellobiose molecules in the stacked arrangement. The oval shows the atoms present in the area encompassing the planes of the two cellobiose molecules. Carbon, oxygen, and hydrogen atoms are shown in gray, black, and white colors, respectively.
8.0.23 The pair of disaccharide molecules was solvated with 2050 TIP3P water24 molecules using xleap in a cubic periodic box. The force field parameters for water and carbohydrates were taken from ff03 and GAFF force fields, respectively.21,24,25 All the simulations reported in this work were performed using AMBER 8.0 molecular dynamics simulation package.23 Conjugate gradient energy minimization was performed to relax the systems. A sequence of constant NVT (constant number of molecules, constant volume and constant temperature) and constant NPT (constant number of molecules, constant pressure and constant temperature) simulations employing the Berendsen weak coupling algorithm were carried out on the systems to equilibrate them at a temperature of 300 K and a pressure of 1 atm.26 These resulted in final density values that are close to 1 g/cc for all the five systems. In these simulations, the time-step for integrating the MD trajectories was set to 1 fs and the OH bond length for the solvent in the systems was constrained using the SHAKE approach.27 The non-bonded and electrostatic interactions were truncated at a cut-off distance of 10 Å. Particle Mesh Ewald (PME) method was employed to account for the long-range electrostatic interactions.28 1.2. Umbrella sampling simulations Straightforward MD simulations, in general, tend to avoid the high energy regions of the configuration space; this aspect leads to inadequate sampling of the phase space. Hence to adequately sample these high energy regions in free energy calculations, an additional biasing harmonic potential (umbrella potential) is employed which is later subtracted from the total free energy to give the net or unbiased free energy change. Umbrella sampling method as developed originally as well as its modifications have been extensively used to estimate the free energy difference between two states at equilibrium.18,20,29–32 In this work, the reaction coordinate for the umbrella sampling simulations for the systems described above (see Fig. 1) was defined as the center of mass separation distance between the molecule pair. In a series of simulations, this distance was varied to represent the process of separation of the disaccharide molecules; these simulations resulted in individual windows in the umbrella sampling procedure. A series of 48 overlapping umbrella sampling windows were considered which spanned the center of mass separation distance between the disaccharide pair from 15 Å to 3.25 Å in increments of 0.25 Å. The separation distance between the molecule centers of mass in each of these windows was maintained near the specified reaction coordinate value by applying a harmonic biasing potential with a force constant of 25 kcal/mol Å2 on the center of mass separation distance. A 500 ps long constant NVT simulation was then performed for equilibrating the system at each separation distance (umbrella sampling window). One set of 1 ns long and two sets of 500 ps long simulation for the cellobiose pair and three sets of 1 ns long production runs for cellobiose– xylobiose and xylobiose pair (initial velocities were reassigned in each case) were then carried out and the trajectories were written every 100 steps. PMF for each set was computed from the generated trajectory data for all of the windows using the Weighted Histogram Analysis Method (WHAM).33 The same procedure was employed to determine the PMF profiles in the stacked arrangement with the exception that, in this case, additional potentials were required to maintain the stacked configuration. Thus, in addition to the umbrella potential on the center of mass distance, harmonic restraining potentials were also applied on the distances between the chosen carbon atoms on the cellobiose molecules (atoms C1, C3, C5, C7, C9, and C11) and on the xylobiose molecules (atoms C1, C3, C5, C7, C9, and C10). These potentials ensured parallel arrangement of the molecules. A 500 ps long constant NVT simulation was performed for equilibrat-
S. Peri et al. / Carbohydrate Research 346 (2011) 867–871
ing the system at each umbrella sampling window. One set of 1 ns long and two sets of 500 ps long simulations for each of the umbrella sampling windows for the cellobiose pair and three sets of 1 ns long simulations for the xylobiose pair were performed. Acknowledgment We gratefully acknowledge the partial financial support of this work by the National Science Foundation (Grant number: CBET 0854463). References 1. Lynd, L. R.; Weimer, P. J.; van Zyl, W. H.; Pretorius, I. S. Microbiol. Mol. Biol. Rev. 2002, 66, 506–577. 2. Zhang, Y. H. P.; Lynd, L. R. Biotechnol. Bioeng. 2004, 88, 797–824. 3. Wyman, C. E.; Dale, B. E.; Elander, R. T.; Holtzapple, M.; Ladisch, M. R.; Lee, Y. Y. Bioresour. Technol. 2005, 96, 1959–1966. 4. Peri, S.; Karra, S.; Lee, Y. Y.; Karim, M. N. Biotechnol. Prog. 2007, 23, 626–637. 5. KlemanLeyer, K. M.; SiikaAho, M.; Teeri, T. T.; Kirk, T. K. Appl. Environ. Microbiol. 1996, 62, 2883–2887. 6. 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–152. 7. Brady, J. W. Carbohydr. Res. 1987, 165, 306–312. 8. French, A. D.; Johnson, G. P. Cellulose 2004, 11, 449–462. 9. French, A. D.; Johnson, G. P. Can. J. Chem./Revue Can. Chim. 2006, 84, 603–612. 10. French, A. D.; Johnson, G. P. Cellulose 2009, 16, 959–973. 11. Hardy, B. J.; Sarko, A. J. Comput. Chem. 1993, 14, 848–857. 12. Liu, Q.; Brady, J. W. J. Phys. Chem. B 1997, 101, 1317–1321. 13. Shen, T. Y.; Langan, P.; French, A. D.; Johnson, G. P.; Gnanakaran, S. J. Am. Chem. Soc. 2009, 131, 14786–14794. 14. Stortz, C. A.; French, A. D. Mol. Simul. 2008, 34, 373–389.
871
15. Stortz, C. A.; Johnson, G. P.; French, A. D.; Csonka, G. I. Carbohydr. Res. 2009, 344, 2217–2228. 16. Leeflang, B. R.; van Kuik, J. A.; Kroon-Bratenburg, L. M. J. In NMR Spectroscopy and Computer Modeling of Carbohydrates: Recent Advances. ACS Symposium Series; 2006; Vol. 930, pp 133–155. 17. Umemura, M.; Yuguchi, Y.; Hirotsu, T. J. Phys. Chem. A 2004, 108, 7063–7070. 18. Palma, R.; Himmel, M. E.; Brady, J. W. J. Phys. Chem. B 2000, 104, 7228–7234. 19. Glennon, T. M.; Merz, K. M. J. Mol. Struct. (Theochem) 1997, 395–396, 157–171. 20. Bergenstrahle, M.; Wohlert, J.; Himmel, M. E.; Brady, J. W. Carbohydr. Res. 2010, 345, 2060–2066. 21. Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem 2004, 25, 1157–1174. 22. Durell, S. R.; Brooks, B. R.; Ben-Naim, A. J. Phys. Chem. 1994, 98, 2198–2202. 23. Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Merz, K. E.; 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, 2004. 24. Mark, P.; Nilsson, L. J. Phys. Chem. A 2001, 105, 9954–9960. 25. Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G. M.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J. M.; Kollman, P. A. J. Comput. Chem. 2003, 24, 1999–2012. 26. Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. 27. Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327– 341. 28. Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089–10092. 29. Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187–199. 30. Trzesniak, D.; Kunz, A. P. E.; van Gunsteren, W. F. Chemphyschem 2007, 8, 162– 169. 31. Chipot, C.; Jaffe, R.; Maigret, B.; Pearlman, D. A.; Kollman, P. A. J. Am. Chem. Soc. 1996, 118, 11217–11224. 32. Hansen, H. S.; Hunenberger, P. H. J. Comput. Chem. 2010, 31, 1–23. 33. Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, R. M. J. Comput. Chem. 1992, 13, 1011–1021.