Multiple Functions of Aromatic-Carbohydrate Interactions in a Processive Cellulase Examined with Molecular Simulation*

Background: Aromatic residues line glycoside hydrolase active sites mediating ligand binding. Results: Binding affinity is significantly altered upon tryptophan to alanine mutation, although relative to the location in the active site. Conclusion: Aromatic-carbohydrate interactions are employed in a variety of functionalities within the purview of ligand binding. Significance: Understanding the functional role of aromatic residues in the active site is necessary for the rational design of new carbohydrate-active enzymes. Proteins employ aromatic residues for carbohydrate binding in a wide range of biological functions. Glycoside hydrolases, which are ubiquitous in nature, typically exhibit tunnels, clefts, or pockets lined with aromatic residues for processing carbohydrates. Mutation of these aromatic residues often results in significant activity differences on insoluble and soluble substrates. However, the thermodynamic basis and molecular level role of these aromatic residues remain unknown. Here, we calculate the relative ligand binding free energy by mutating tryptophans in the Trichoderma reesei family 6 cellulase (Cel6A) to alanine. Removal of aromatic residues near the catalytic site has little impact on the ligand binding free energy, suggesting that aromatic residues immediately upstream of the active site are not directly involved in binding, but play a role in the glucopyranose ring distortion necessary for catalysis. Removal of aromatic residues at the entrance and exit of the Cel6A tunnel, however, dramatically impacts the binding affinity, suggesting that these residues play a role in chain acquisition and product stabilization, respectively. The roles suggested from differences in binding affinity are confirmed by molecular dynamics and normal mode analysis. Surprisingly, our results illustrate that aromatic-carbohydrate interactions vary dramatically depending on the position in the enzyme tunnel. As aromatic-carbohydrate interactions are present in all carbohydrate-active enzymes, these results have implications for understanding protein structure-function relationships in carbohydrate metabolism and recognition, carbon turnover in nature, and protein engineering strategies for biomass utilization. Generally, these results suggest that nature employs aromatic-carbohydrate interactions with a wide range of binding affinities for diverse functions.

Protein-carbohydrate recognition is a ubiquitous phenomenon in biology that is important for a vast range of biochemical functions (1)(2)(3)(4)(5)(6)(7). A common means of selective protein binding to carbohydrates is accomplished by the aromatic amino acids, tryptophan, tyrosine, and phenylalanine, wherein aromatic rings can stack with planar faces of carbohydrate rings via carbohydrateinteractions (8 -22). To date, many studies have been conducted to determine the effect of aromatic-carbohydrate interactions in model (9 -15) and real systems (16 -22), with the aim of quantifying the strength and types of interactions for this important biological recognition pattern. A particularly dramatic example of the prevalence of aromatic-carbohydrate interactions is the ubiquity of these motifs observed throughout the 125 families of glycoside hydrolases (GHs), 2 which catalyze hydrolysis of the glycosidic linkages in carbohydrates, and the related 63 families of carbohydrate-binding modules, which are responsible for increasing the binding affinity of GH enzymes to carbohydrate substrates by many orders of magnitude (8,23). In each of the structurally characterized GH families, aromatic residues line the enzyme tunnels and clefts for binding, processing, and stabilization of carbohydrate ligands. In carbohydrate-binding modules, aromatic residues are commonly found on the substrate-binding faces, which form favorable electrostatic and hydrophobic interactions with carbohydrates (24,25).
In nature, processive GHs provide the majority of direct hydrolytic potential in the conversion of insoluble carbohydrate polymers, such as cellulose and chitin, to soluble sugars for utilization as food sources (26 -30). From crystal structures of processive chitinases and cellulases, it is known that these enzymes typically exhibit long tunnels or deep clefts lined with aromatic and polar residues that can thread single polymer chains to their active sites for hydrolysis (28 -30). Processive GHs are generally hypothesized to bind to their insoluble, polymeric substrates, find a free chain end, thread the chain into the enzyme tunnel to form the catalytically active complex, conduct hydrolysis, expel the soluble product, and then repeat this catalytic cycle processively (31,32) until the enzyme reaches the end of the polymer chain, becomes inhibited by obstacles on the surface, or detaches from the polymer chain (33). It should be noted, however, that processive cellulases have also been observed to exhibit endo-initiation activity (34,35). The primary hypothesis in the literature to date is that aromatic residues in processive enzyme tunnels facilitate carbohydrate processivity, but the thermodynamic and kinetics of processivity have not been elucidated directly at the molecular level (28 -30).
Several studies have examined the role of aromatic residues in cellulase and chitinase tunnels with respect to processivity and enzyme activity on substrates of varying accessibility (36 -42). Watanabe et al. (41) conducted many of these original studies related to various chitinases. More recently, Horn et al. (36) mutated two tryptophans in the Serratia marcescens family 18 chitinase B (ChiB). They demonstrated that a single mutation (W97A) directly upstream of the catalytic site in ChiB significantly reduced the ability of the enzyme to degrade crystalline chitin, but increased the conversion rate by 29-fold on soluble chitosan. The authors also measured the degree of processivity between the wild type ChiB and the W97A mutant and demonstrated that the W97A mutation converted ChiB to a nonprocessive enzyme, even on soluble substrates (36). The same group later demonstrated similar trends for another S. marcescens chitinase, ChiA (38). Overall, this work suggests that the ability for an enzyme to be processive on insoluble substrates comes at a significant detriment to its substrate turnover frequency.
For aromatic-carbohydrate interactions in cellulases, a tryptophan residue (Trp-272) at the entrance to the tunnel of the processive family 6 cellulase, Cel6A, from Trichoderma reesei (an anamorph of Hypocrea jecorina), was mutated to alanine and aspartic acid (39). In both cases, the activity of the mutant on crystalline cellulose was significantly reduced, whereas the activity on amorphous cellulose was not affected (39). The authors interpret these results to suggest that the entrance tryptophan residue in Cel6A is crucial in binding to cellodextrin chains in crystalline substrates but plays little functional role when the substrate is more accessible. Similar observations have been made for Cel7A (43,44). Additionally, Vuong and Wilson (37) mutated a tryptophan (Trp-464) in Thermobifida fusca Cel6B directly upstream of the catalytic site to alanine and tyrosine. For both mutations, they observed reduced activity on crystalline cellulose and increased activity on amorphous and soluble substrates, although the extent of processivity on filter paper was similar to the wild type. The authors concluded that Trp-464 helps to bind the cellulose chain to the active site tunnel, but based on increased activity on various accessible substrates, this function is not necessary for noncrystalline substrates. It is worthwhile to note that processivity has been measured in distinctly different ways for the cellulases and chitinases discussed above and thus may not be directly comparable.
We recently combined observations from the work of Horn et al. (36) on chitinases with a computational analysis of cellu-lose and chitin decrystallization to propose a microkinetic model of processive enzyme action (45)(46)(47). We hypothesize that the work that an enzyme must conduct to decrystallize a carbohydrate polymer from a crystalline substrate is directly offset by a favorable ligand binding free energy (31,32,45,46). By removing aromatic residues in processive enzyme tunnels, the ligand binding free energy may be reduced sufficiently to switch the function from a processive enzyme to a nonprocessive enzyme; however, in the case of modified substrates that are either more accessible (48) or easier to decrystallize (45,49), it is a direct extension from this work that engineered enzymes and enzyme mixtures significantly different from those found in nature may offer significant gains in conversion rates for industrial processes using modified cellulose. This, in turn, could lower production costs for the burgeoning renewable fuels and chemicals industry. Quantitative understanding of the effects of the aromatic-carbohydrate interactions in cellulase tunnels is thus a key factor in the rational design of new cellulases and other carbohydrate-active enzymes. Understanding the role that aromatic-carbohydrate interactions play in the structure and function of carbohydrate-active enzymes is similarly important. This remains an open question, for which molecular level simulations can offer unique insights.

EXPERIMENTAL PROCEDURES
We use thermodynamic integration (TI) (50,51) to calculate the change in ligand binding free energy for removal of aromatic residues in the catalytic tunnels of the well characterized processive family 6 T. reesei Cel6A cellulase (29,39,52). Cel6A is a nonreducing end-specific enzyme with six binding sites numbered ϩ4 (upstream of the catalytic site) to Ϫ2 (down the tunnel from the catalytic site), as shown in Fig. 1. The Cel6A catalytic domain tunnel exhibits tryptophan residues (Trp-135, Trp-269, Trp-272, and Trp-367), which we mutate individually to alanine.
The TI calculations were performed using NAMD with dual topology methodology, wherein a hybrid residue containing atoms from both the tryptophan and alanine interacts with the system over a scaled coupling parameter during the transformation from reactant to product (50,53,55). The atoms of the hybrid residue do not interact with each other over the course of the simulation. In the TI simulations, the tryptophan is referred to as the reactant and the alanine as the product. The change in free energy was determined through numerical inte- gration of sets of simulations probing the electrostatics and van der Waals interactions separately. The free energies, G, for each state were then combined according to Equation 1 to obtain the relative ligand binding free energy, ⌬⌬G, where the states are wild type (reactant) to mutant (product).
The alchemical cycle for measuring ⌬⌬G is shown under supplemental Figs. S1 and S2 along with error analysis details. From the ⌬⌬G for each mutation, we determined the binding affinity for all point mutations as in Equation 2, where ⌬⌬G is the relative binding free energy, R is the universal gas constant, T is absolute temperature, and K mut and K WT are the partition coefficients for the mutant and wild type structures, respectively. All of the TI simulations were started using a snapshot at 18 ns from molecular dynamics (MD) simulations as described below. Simulations from which the TI data were collected were equilibrated for an additional 0.5 ns prior to collection of 14.5 ns of TI data. In addition to the TI calculations, we use long MD simulations (250 ns) to quantify the energetics and fluctuations along the catalytic tunnels of the wild type and four mutant enzymes. These MD simulations provide insights into the role of the aromatic residues as a function of location in the tunnel in ligand acquisition and stabilization. CHARMM (56) was used to build, minimize, and equilibrate the protein structure from the wild type crystal structure (Protein Data Bank code 1QK2) (57). Final system sizes were on the order of 52,000 atoms, of which roughly 46,000 atoms were water. Subsequent MD simulations were performed using NAMD in the NVT ensemble at 300 K for 250 ns. The CHARMM27 force field with the CMAP correction (56,58,59) was used to describe the protein, whereas glycans and cellulose used the CHARMM35 carbohydrate force field (60 -62). The system preparation and computational methods are described in more detail under supplemental data. Table 1 lists the predicted change in relative binding free energy and binding affinity from the TI calculations. These results reveal several interesting observations. The ligand binding free energies of the tunnel entrance tryptophan, Trp-272, and the product side tryptophan, Trp-135, are quite large at 3.8 and 9.4 kcal/mol, respectively. Surprisingly, the results for Trp-367 and Trp-269, at the ϩ1 and ϩ2 sites, respectively, show only minor changes in relative ligand binding free energy. The relative binding affinity, K wt /K mut , was calculated from the ligand binding free energy as discussed under supplemental Table S1. The binding affinities of W135A and W272A indicate a clear preference for the wild type in these positions. The W367A and W269A mutations only slightly alter the relative binding affinity in comparison to W135A and W272A.

RESULTS
To further elucidate the role of the various aromatic residues in the Cel6A tunnel, we conducted five separate 250-ns MD simulations of the individual mutants and the wild type enzyme. From these simulations, we calculated the root mean square fluctuation (r.m.s. fluctuation) of the protein (Fig. 2) and of the ligand on a binding site basis (Fig. 3). The interaction energies of the protein with the ligand and the hydrogen bonds between the ligand and protein were also calculated and are provided under supplemental Fig. S3. Analyses of the glycosylation effects on the enzyme structure, which are minimal, are provided under supplemental Figs. S11-13. These analyses included comparisons of the r.m.s. fluctuation and r.m.s. deviation of the entire enzyme along with comparisons of ligand hydrogen bonding and interaction energies on a binding site basis.
Evaluating the r.m.s. fluctuation of the protein over the course of the 250 ns MD simulations reveals interesting behavior in the loops forming the active site tunnel in some cases. The W135A mutant ( Fig. 2A) exhibits increased fluctuations in both the 170 -190 and 400 -410 residue regions. Each of these residue segments correspond to tunnel-forming active site loops. The former loop, residues 172-182, has been previously linked to a conformational change as a function of active site mutations, increasing active site solvent accessibility although having no effect on ligand conformation (52). The W272A mutant (Fig. 2C) similarly exhibits increased fluctuation of the 400 -410 active site loop region with minimal differences in the other loop region. This is likely a result of the increased interaction of the protein with the ligand as a result of the increased ligand flexibility in the ϩ4 binding site, as will be shown below in Fig.  3. Only moderate effects on the fluctuation of the protein are observed in the W269A and W367A mutants with respect to wild type (Fig. 2, B and D).
The r.m.s. fluctuation of the ligand over the course of the MD simulations also lends insight into the molecular level aspects contributing to the relative ligand binding free energy. Comparing the MD simulations of the W272A mutant to that of the wild type Cel6A enzyme reveals that mutation of tryptophan to alanine in this binding site substantially alters the cellodextrin fluctuations, as shown in Fig. 3, which suggests that this residue is crucial for stabilizing the ligand at the entrance of the tunnel, as hypothesized by experimental studies (39). Fig. 3 also illustrates destabilization of the entire ligand in all six binding sites for the W269A and W367A mutants. It is particularly interesting that for these two mutants, the fluctuations in neighboring binding sites (e.g. the Ϫ1 and ϩ2 sites for the W367A mutant) are affected, whereas for W135A and W272A the mutations to alanine have less of an impact on neighboring site fluctuations. The interaction energies of either wild type tryptophans or mutant alanines also exhibit this trend as illustrated under supplemental Fig. S4. The higher fluctuation of the Ϫ1 binding site in W269A and W367A than in all other mutations likely contributes entropically to the decreased relative binding free energy of these two mutants relative to W135A and W272A. All mutations result in increased fluctuations at the Ϫ2 binding site, and all reactant side (upstream of active site) mutations increase fluctuations at the ϩ4 binding site. The change in ligand binding free energy is most dramatic for W135A. Because it is located on the product side of the tunnel, it is directly involved in cellobiose binding, which is related both to stabilization of the product during catalysis and product inhibition (i.e. through competitive binding) (63). Although evaluation of the interaction energy of the entire protein with the individual binding sites (supplemental Fig. S3) yields little additional insight in terms of the relative ligand binding free energy, it does provide further evidence for the importance of the binding affinity of the product side residues, where the interaction energy of glucose in the Ϫ2 binding site for the wild type and the four mutants is substantially greater than any other site along the tunnel. This is the same general result that Mertz et al. (64) obtained using computational docking techniques. Reduction of this evaluation to either tryptophan or alanine with the glucose unit, as shown under supplemental Fig. S4, demonstrates that removal of any of the tryptophans in the tunnel has a significant local effect on interaction energy.
As mentioned above, the relative ligand binding free energies for W269A and W367A in Cel6A are small in comparison to  W135A and W272A. These results suggest that tryptophans at the ϩ1 and ϩ2 sites in Cel6A are less important for direct binding to the ligand. Alternatively, these residues may create the necessary physical shape of the tunnel in the catalytically active complex, and play a role in the ring distortion at the Ϫ1 site wherein the glucose exhibits the 2 S O conformation instead of the chair conformation, which is energetically favorable in solution (57). To examine this hypothesis, we calculated the Cremer-Pople ring amplitude (65) of the glucose residue in the Ϫ1 site from MD simulations for the wild type and four mutants as shown in Fig. 4. The 2 S O conformation, corresponding to that of glucose in the Ϫ1 site of Cel6A, is defined by an amplitude of 0.73 Å. The chair conformation, 1 C 4 or 4 C 1 , is defined by an amplitude of 0.57 Å (66). As the information presented in Fig. 4 was obtained from MD trajectories, the ring pucker amplitude fluctuates around these conformational definitions over time. In all cases except the W272A mutant, the conformation of the glucose ring in the Ϫ1 subsite relaxes to the chair conformation. Interestingly, the wild type Ϫ1 subsite ring cycles between chair and skew conformations, with skew as the preferential conformation. This suggests that the tryptophans local to the catalytic site, Trp-135, Trp-269, and Trp-367, are crucial for maintaining the ring conformation needed for catalysis, and that Trp-272 does not contribute to this function. Given the similarity of the tunnel region surrounding the catalytic site in W272A compared with wild type, we hypothesize that W272A MD simulation would likely exhibit similar Ϫ1 subsite conformational fluctuations as observed in the wild type upon extended sampling time.
Detailed analyses of the hydrogen bonding of active site residues and general active site behavior can also provide insights into the results obtained from the relative ligand binding free energy and the long MD simulations. In the wild type MD simulation in the Ϫ1 binding site, we observed hydrogen bonding of residue Asp-221, the catalytic residue, to Asp-175 over the course of the simulation (supplemental Fig. S7). The hydrogen bonding of Asp-221 to another carboxylic acid residue allows for the sharing of a proton between the two residues. When this hydrogen bond is broken, Asp-221 is able to act as a catalytic acid pushing a proton onto a free electron pair on the glycosidic oxygen. Hydrogen bonding of these two residues does not occur in the W135A, W269A, and W367A MD simulations. For the wild type, Asp-221/Asp-175 hydrogen bonding is disrupted concomitantly with the conformational change of the Ϫ1 glucose ring. This crucial hydrogen bond is relatively stable over the course of the W272A MD simulation. In the wild type simulation, the distance of the carboxyl oxygen of Asp-221 from the glycosidic oxygen connecting the cellulose chain between the ϩ1/Ϫ1 binding sites indicates hydrogen bonding between Asp-221 and Asp-175 ceases as Asp-221 moves away from Asp-175 and closer to the glycosidic oxygen. When the Ϫ1 glucopyranose ring changes from the skew to chair conformation in the wild type, this distance exhibits an average value of ϳ4 Å. In the W135A, W269A, and W367A simulations, the Asp-221 residue quickly moves away from the protein toward the ligand preceding the glucopyranose ring conformational change, likely as a result of the increased volume of the active site tunnel near the catalytic acid. This movement of Asp-221 toward the glycosidic oxygen is consistent with the observations of Koivula et al. (52). We also observe hydrogen bond formation between the catalytic acid, Asp-221, and the primary alcohol group of the Ϫ1 site glucose following a conformational change in this ring. Hydrogen bonding of the glucose primary alcohol to the Tyr-169 residue is consistent regardless of glucose conformation.
Related to the increased fluctuations observed in the 170 -190 residue region in Fig. 2A Fig. S5B), we did not observe a similar conformational change over the course of the 250-ns MD simulations, likely a function of the limited time scale for conformational sampling. We suspect, however, based on the r.m.s. fluctuation in Fig. 2A, that W135A is undergoing a conformational change of the 172-182 region. Decreased hydrogen bonding between Tyr-169 and Arg-174 is observed in the W135A simulation in comparison to wild type, although all of the mutants exhibit this to some extent (supplemental Fig. S9). Opposite to that observed by Koivula et al. (52), we find increased hydrogen bonding between Arg-174 and Asp-175 for W135A (supplemental Fig. S10). It is likely that this behavior, which we attribute to large conformational fluctuations, is a function of the mutation.
A final analysis of the wild type residue cross-correlation obtained from normal mode analysis (Fig. 5) shows strong positive correlations of residues Trp-135, Trp-269, and Trp-367 with the active site region, whereas Trp-272 is much less correlated. The positive correlation of Trp-135 with the active region, residues 163 to 178, includes a portion of the active site loop associated with increased fluctuations, as discussed above. Residue Trp-269 is correlated with section 221-231, and residue Trp-367 is correlated with section 399 -403. The correla- tion value of each of these three residues, Trp-135, Trp-269, and Trp-367, is generally in the range of 0.6. Notably, the positive correlation of residue Trp-272 with any region of the active site is less than 0.4. Finally, whereas the MD simulation results indicate the dynamics of the 400 -410 loop region are affected by the W135A and W272A mutations, normal mode analysis does not show a cross-correlation. This is likely a shortcoming of the capabilities of normal mode analysis, as it is well accepted that normal mode analysis may miss some loop activity, particularly in small sections such as the 10-residue loop here (67). The results from normal mode analysis further verify the roles of the aromatic residues as suggested by TI and MD simulations.

DISCUSSION
Overall, this work presents several roles for aromatic-carbohydrate interactions in cellulase tunnels. First, residues near the entrance of the tunnels in Cel6A and Cel7A have been shown to be crucial for hydrolyzing insoluble substrates experimentally (39,43,44). Here we show that removal of this residue in Cel6A, Trp-272, leads to significantly reduced binding affinity and greater fluctuations of the oligomer chain in the tunnel. Further MD simulations, including that of the entire Cel6A complex on a crystalline cellulose microfibril, similar to work done on Cel7A (68), are warranted to examine how mutation of this tunnel entrance tryptophan affects the behavior and stability of the catalytically active complex on an insoluble substrate. Second, we have demonstrated that the aromatic residue on the product side of the catalytic domain tunnel in Cel6A plays a significant role in binding, as the removal of this residue dramatically affects the ligand binding free energy. Specifically, our MD simulations of the Cel6A wild type indicate that the interactions at the Ϫ2 site are much higher than anywhere else along the cellodextrin chain, suggesting that this aromatic residue plays a significant role in stabilization of the product state dur-ing catalysis and likely contributes to product inhibition via association of the product to the Ϫ1 and Ϫ2 sites.
Last, and most surprisingly, our results for the ϩ1 and ϩ2 site residues directly preceding the hydrolysis site suggest functionality not directly associated with tight binding. Rather, we propose the roles of Trp-367 and Trp-269 include structuring the active site tunnel appropriately for ideal positioning of the ligand for catalysis. Mutation of these residues, specifically Trp-367, would almost certainly result in a decrease in activity on crystalline cellulose substrates. Without the requisite distortion at the Ϫ1 binding site, the glycosidic oxygen is further from the catalytic residue, shown in Fig. 1A, and the anomeric carbon is in a less favorable position for nucleophilic attack, and therefore not appropriately preactivated for catalysis (69). Interestingly, an experimental study of T. fusca Cel6B, a family 6 GH similar to T. reesei Cel6A, showed that a mutation corresponding to W367A (W464A in T. fusca Cel6B) resulted in a significant decrease in both activity and binding to bacterial microcrystalline cellulose, with a more moderate reduction in activity on phosphoric acid-swollen cellulose (37). Processivity of the mutated T. fusca Cel6B remained similar to wild type. The T. fusca Cel6B experimental results combined with our results suggest that mutations of tryptophans in the ϩ1 and ϩ2 binding sites of family 6 cellulases may have little effect on processivity, which in this GH family is linked to low ligand binding free energies of aromatic residues near the catalytic site. This observation can be contrasted to the observed processivity differences for GH family 18 chitinases, further suggesting that, in general, the results from Horn et al. (36) are linked closely to higher ligand binding affinities. This will be examined in future work.
From the simulation results obtained in this study and previous experimental results (36,37), we hypothesize that the thermodynamic nature of aromatic-carbohydrate interactions in glycoside hydrolases is a function of multiple factors, including the ligand conformation (both the glycosidic linkage conformations and the ring conformations), the ligand type, the hydrolytic mechanism (which dictates the position of the ring distortion), the nature of the binding site (i.e. a tunnel or a cleft with surrounding binding site residues), and the GH protein fold (i.e. the GH family). This hypothesis may suggest discontinuity in the strength of aromatic-carbohydrate interactions in the equivalent ligand binding sites across GH families. For example, evaluation of the potential energies of the glycosidic bond between the Ϫ1 and ϩ1 sites from the ␤-1,4 / angles in the crystal structures of Cel6A and the family 18 chitinase ChiB indicates the potential energy of the ChiB ligand conformation is 1-2 kcal/mol higher than the cellodextrin chain in the Cel6A crystal structure, which may manifest itself as a higher ligand binding free energy in ChiB than in Cel6A (70,71). Additionally, ChiB hydrolyzes the glycosidic linkage through a substrate-assisted retaining mechanism, wherein the substrate acts as a nucleophile (72). In contrast, Cel6A uses a single displacement inverting mechanism, where the nucleophilic attack is made by a water molecule and a catalytic base is required to abstract the proton from the water, increasing its nucleophilicity (73,74). Each of these differences likely directly contributes to the residue functionality in active sites of GHs. It is possible that aromatic-carbohydrate functional roles may only be specific to GHs of sufficiently similar fold and hydrolysis mechanisms, which we will examine in future work.
There are several direct extensions to this work and related work (36,38,45,46). In terms of converting processive enzymes to nonprocessive enzymes, as shown by Horn et al. (36), it would be of significant interest to determine the relative extents of processivity of a processive enzyme with the appropriate aromatic mutations (e.g. W97A for ChiB and possibly W367A for Cel6A) relative to a corresponding natural nonprocessive enzyme (i.e. an endoglucanase). Nonprocessive enzymes, such as endoglucanase Cel6B from Humicola insolens (75), would be ideal for such comparisons within this family. In a future study, we will examine the role of aromatic residues in corresponding processive and nonprocessive enzymes via a similar TI approach. Moreover, experimental measurements using techniques such as isothermal titration calorimetry (22,76,77), could potentially validate the trends observed in the changes in ligand binding free energy calculated here with TI. However, this would require substrates with modified linkages (e.g. with thioether bonds) such as those used in structural studies so that the catalytically active complex can be formed without hydrolytic activity. Calculation of the entire range of torsional parameters for glycosidic linkages of chitobiose and cellobiose, and similarly for ring distortion, using a consistent and sufficiently accurate model would also be another step toward identifying the energetic contributions to ligand binding free energies across the many GH families.

CONCLUSIONS
Free energy calculations, MD simulations, and normal mode analysis have been used to examine the functional roles of aromatic residues in the active site of the well characterized T. reesei processive cellulase, Cel6A. Residues associated with ligand acquisition and product stabilization/inhibition have the greatest influence on ligand binding free energy and thus are likely to have a large impact on enzyme activity and processivity. Unexpectedly, the ϩ1 and ϩ2 site tryptophan residues appear not to participate in the tight binding to the ligand, but rather in stabilization of the ligand for the formation of the catalytically active substrate conformation. Comparison of the results here for Cel6A to experimental results for family 18 chitinases as a function of binding site indicates the role of aromatic-carbohydrate interactions in processivity and ligand binding may not be reduced to a model simply based on tunnel position relative to the location of catalysis. As a result, protein engineering strategies focused on improved biomass utilization may not be as generalized as previously expected and may require tailoring specific to the enzyme family. Based on the results of this study, it is apparent that many factors contribute to the thermodynamic nature of aromatic-carbohydrate interactions in GH families, which suggests additional studies on the strength of aromatic-carbohydrate interactions in the context of the protein folds, hydrolytic mechanisms, and the substrate types. More generally, this study suggests that nature employs aromatic-carbohydrate interactions over a very large range of binding affinities for multiple functions.