The impact of Trichoderma reesei Cel7A carbohydrate binding domain mutations on its binding to a cellulose surface: a molecular dynamics free energy study

A critical role of the Family 7 cellobiohydrolase (Cel7A) carbohydrate binding domain (CBD) is to bind to a cellulose surface and increase the enzyme concentration on the surface. Several residues of Trichoderma reesei Cel7A CBD, including Y5, N29, Y31, Y32 and Q34, contribute to cellulose binding, as revealed by early experimental studies. To investigate the interactions between these important residues and cellulose, we applied a thermodynamic integration method to calculate the cellulose–Cel7A CBD binding free energy changes caused by Y5A, N29A, Y31A, Y32A and Q34A mutations. The experimental binding trend was successfully predicted, proving the effectiveness of the complex model. For the two polar residue mutants N29A and Q34A, the changes in the electrostatics are comparable to those of van der Waals, while for three Y to A mutants, the free energy differences mainly come from van der Waals interactions. However, in both cases, the electrostatics dominates the interactions between individual residues and cellulose. The side chains of these residues are rigidified after the complex is formed. The binding free energy changes for the two mutants Y5W and Y31W were also determined, and for these the van der Waals interaction was strengthened but the electrostatics was weakened.


Introduction
Biofuels produced from biomass are a promising and environmentally friendly alternative to fossil fuels, but the process of producing biofuels faces substantial challenges, mainly due to its high cost compared to traditional energy production. One approach is to hydrolyze cellulose-the main polymer in biomass-to sugars by cellulases, which are then fermented to ethanol. Cellulose is the most abundant bioresource on Earth and the major structural component of plant cell walls. In nature, cellulose molecular chains are linear polymers composed of β-1,4 linked D-glucose units.
The enzymes involved in the hydrolysis of cellulose consist of endoglucanases (EG), exoglucanases or cellobiohydrolases (CBH), and β-glucosidases. Among them, endoglucanases are responsible for breaking down the intramolecular β-1,4 glucosidic bonds of cellulose chains to produce new chain ends; exoglucanases hydrolyze cellulose chains from their ends and release soluble cellobiose; and βglucosidases catalyze cellobiose to glucose [1,2]. The enzymes secreted by the filamentous fungus Trichoderma reesei (T. reesei), which are highly effective against native and crystalline cellulose [3], are well studied. T. reesei produces a mixture of cellulases consisting of two CBHs (CBH I or Cel7A and CBH II or Cel6A) and several EGs [4,5]. The enzyme Cel7A is composed of a large catalytic domain (CD) and a small cellulose-binding domain (CBD) connected via a glycosylated peptide linker [6]. The CD of Cel7A, containing a 50 Å tunnel-shaped active site, binds a cellulose chain and hydrolyzes the glycosidic bond to release cellobiose [3,[7][8][9].
The three-dimensional solution structure of the Cel7A CBD from T. reesei has been determined by nuclear magnetic resonance (NMR) spectroscopy [10], which Electronic supplementary material The online version of this article (doi:10.1007/s00894-011-1167-4) contains supplementary material, which is available to authorized users. revealed that several aromatic residues, Y5, Y31 and Y32, form a hydrophobic surface. However, there is no experimental structure of the Cel7A-or CBD-cellulose complex due to the insoluble nature of cellulose. The experimental study by Linder and coworkers [4] has shown that these three tyrosine residues are critical for binding to crystalline cellulose, while N29 and Q34 also contribute but are less important. Recently, the CBD cellulose complex has been modeled in order to investigate the binding mechanism by molecular dynamics (MD) simulations. It has been suggested that CBD may assist with the recognition of the cellulose chain end and facilitate the processivity of Cel7A [11][12][13]. Beckham et al. predicted that electrostatic interactions dominate the energy changes in the Cel7A CBD processivity [12]. However, it is very important to validate the CBD-cellulose complex model against experimental data.
Free energy methods provide a rigorous way to obtain important and readily measurable thermodynamic observables, such as binding free energy, from MD simulations [14,15]. Chemical alchemy is an important and accurate class of free energy calculations. This approach involves simulating a transformation of one chemical system to a different one along a nonphysical path to compute the free energy change associated with the transformation [14]. Among the alchemical methods, thermodynamic integration (TI) [15][16][17][18][19] has grown in popularity over the last few years. A parameter 1, ranging from zero to one, is introduced to define the intermediate states, with 1 =0 (1) corresponding to the initial (final) state. A series of simulations are performed for all the states, and then the partial derivative of the Hamiltonian over 1 is estimated for each state and integrated over 1 to yield the free energy change of the transformation. Although the computational demand is still large, TI has come of age with increased computer power and the wide use of parallel computers. One important application of the TI method is the study of the consequences of protein sidechain mutations, which has shown promising results [16,[20][21][22].
In this work, we calculated the relative binding free energy differences between the mutants (N29A, Q34A, Y5A, Y31A and Y32A) and wild-type CBD that bind to cellulose. Qualitative agreement between the experimental and computed free energy differences was observed, validating the effectiveness of the CBD cellulose-binding model. The absolute binding free energy calculation of CBD involves the annihilation of CBD on the cellulose surface and in water, which is a formidable task, so it was not pursued in this work. On the other hand, the absolute binding free energy was also difficult to determine (if not impossible) experimentally because of the insolubility of cellulose.

Preparation of the system
The starting coordinates for the CBD of Cel7A were obtained from the NMR structure [10] (PDB code 1CBH). The cellulose slab was modeled based on the cellulose Iβ crystal structure presented by Nishiyama et al. [23]. As shown in Fig. 1, the cellulose is three layers deep with five chains in each layer, and each chain includes 12 glucose units. Consistent with early studies [12,13], the CBD is aligned on a cellulose surface with the aromatic residues Y5, Y31 and Y32 within 3 Å of the cellulose surface simultaneously, and the Y5, Y31 and Y32 rings are stacked over the glucose rings of cellulose. The complex model is stable during the MD equilibration (see the section "Results and discussion").

Molecular dynamics simulations
The simulations were performed using the SANDER module in the AMBER 10 program [24], with the AMBER99sb [25] force field used for the protein and the GLYCAM_06 force field [26] [27] for the cellulose. The complex was soaked in a rectangular box with 12 Å TIP3P water molecules on each side using the LEAP module in AMBER. Na + ions were added to neutralize the system. These yielded about 35,000 atoms for the simulation system. The 1-4 electrostatic (SCEE) and nonbonded (SCNB) scaling factors were set to their default values, with SCEE=1.2 and SCNB=2.0. It was realized later that the SCEE and SCNB scaling factors of cellulose should be set to one to yield the correct rotamer population, which is not done properly in AMBER 10. However, in our work, cellulose is restrained by using a harmonic potential with a force constant of 2.0 kcal mol −1 A −2 on the C1 and C4 atoms, essentially removing the dihedral fluctuations along the C1-O4 and C4-O4 bonds.
Energy minimization was performed to remove bad contacts, and then the system was heated from 0 K to Fig. 1 Structural model of the CBD and cellulose complex. The side chains of residues Y5, N29, Y31, Y32 and Q34 are shown in stick representation, with the carbon atoms colored in cyan and the oxygen atoms in red. The cellulose was modeled by three layers of cellulose chains 300 K in a 50 ps constant volume MD simulation, after which a 1 ns NPT simulation was run to equilibrate the water so that the water density reached 1.0 kg/L. The particle mesh Ewald (PME) method [28] was applied to calculate long-range electrostatics interactions. The SHAKE method [29] was applied to constrain all the covalent bonds involving hydrogen atoms. Periodic boundary conditions were applied. A nonbonded pair list cutoff of 12.0 Å was used. The time step used to integrate the equations of motion was 0.002 ps. Coordinates were saved every 0.4 ps along the 20 ns production run, and the trajectories were analyzed with the PTRAJ module of AMBER.

Free energy calculation by thermodynamics integration
According to the thermodynamic cycle, four states should be considered during the calculation of the relative binding free energy. The relative binding free energy difference ΔΔG bind between the wild-type CBD cellulose complex (WT-CBD-Cel) and a mutant CBD cellulose complex (MU-CBD-Cel) could be obtained from ΔG B − ΔG A (Fig. 2). However, as discussed above, directly calculating the binding free energies ΔG A and ΔG B is a formidable task, which can be avoided by designing a thermodynamic cycle as shown in Fig. 2. We can get ΔΔG bind by calculating the difference in the free energy changes for the free (ΔG 1 ) and bound (ΔG 2 ) states.
ΔG 1 (ΔG 2 ) corresponds to the free energy change in the alchemical process of converting the wild type to the mutant residue in the absence (presence) of the cellulose. It can be computed using the following equation: where the subscripts W and M denote the wild type and mutant state, V is the potential energy, and 1 is a parameter that ranges from 0 (wild type) to 1 (mutant). The binding affinities of the five mutants, including Y5A, N29A, Y31A, Y32A and Q34A, to cellulose were determined experimentally. The free energy calculations for the mutation of each of these five residues to an alanine were performed by dividing the process into three steps, as commonly employed in TI simulations. The first step is to remove the partial charges of the side chain atoms of the residue subject to mutation, except for the C β atom and the two H β hydrogens that exist in the alanine residue. In step two, the van der Waals (VDW) interactions with the atoms whose charges are removed in step one are eliminated (this is described in more detail below). In step three, the third H β hydrogen atom attached to the C β of the alanine residue gets the partial charge back. During the VDW transformation in step two, a linear mixing soft core potential is applied in which nonbonded VDW interactions are represented by a 1 -dependent modified LJ equation [14], where ε is the well depth and σ is the collision diameter, α is an adjustable constant which is set to 0.5 in this study, and r ij is the distance between atoms i and j. Soft core scaling provides smooth 1 curves that are best suited for commonly used numerical integration schemes [14]. Equation 3 is the potential function used for atoms that disappear during a VDW transformation. When 1 =0, it has the same form as the normal LJ equation; for 1 values close to 1, it provides a smooth interaction function that allows other atoms to be situated on top of a vanishing atom, with a finite energy penalty that goes to zero at 1 =1 [14]. The last configuration of the CBD-cellulose complex from the 20 ns simulation was selected as the starting structure for the alchemical simulations to calculate ΔG 2 . A CBD system without cellulose was prepared with a size of 11000 atoms for the ΔG 1 calculations. In general, for each 1 in the three-step transformation, we carried out a 100 ps constant-volume equilibration run followed by a 400 ps production run where the energy derivative ∂V/∂1 for each 1 was collected. 1 values were set to 0.025n (n=1, 2, …, 39) in each of the three steps; therefore, there were a total of 117 windows for each mutation. The free energy change associated with each of the three steps was calculated by numerically integrating the derivative data. The integration was done using an in-house perl script. The errors were estimated by splitting the trajectories into four equal blocks and calculating the corresponding free energy differences.
Besides the five single mutations (Y5A, N29A, Y31A, Y32A and Q34A), we also performed two more alchemy simulations (Y5W and Y31W) to determine the binding affinity change after tryptophan substitution. Presumably, Fig. 2 Thermodynamic cycle for calculating the binding free energy difference between the wild-type and mutant CBD binding to cellulose the tryptophan ring prefers to be parallel to the cellulose surface to maximize the VDW interaction, similar to the tyrosine phenol ring conformations on a cellulose surface. From this reasoning, we designed two conformations with (χ 1 , χ 2 ) dihedrals equal to (−36°, 113°) and (−36°, -66°) for W5, and (−69°, 83°) and (−69°, -96°) for W31, respectively (see Fig. 3), as the starting structures, but they were not restrained during the TI simulations. In total, there were nine TI simulations with a total simulation time of~1 μs.

Side-chain conformational entropy
The side-chain conformational entropy, S conf , was calculated using T-Analyst [30]. S conf is defined as where R is the gas constant, p i is the probability distribution of angle i, which refers to the side chain (χ 1 , χ 2 ) of the mutated residues in this work. The distributions were built from 20 ns molecular dynamics trajectories of free and cellulose-bound CBD. For each residue mutation, the entropy change was estimated using the T-Analyst program, and the entropic contribution to the binding free energy was calculated with the temperature set to 300 K. The bin size for each p i was 5°for side-chain dihedrals [30].

ns MD
A putative binding position for CBD was proposed where the aromatic rings of the three tyrosine residues are nearly collinearly aligned along the cellulose fiber axis [31]. The close contact of the tyrosine and sugar rings is in concert with the fact that three tyrosine residues are critical in the CBD cellulose binding based on the experimental study. The Cel7A processivity action against insoluble highly crystalline cellulose is unambiguously toward glucose nonreducing ends from reducing ends [9], which effectively restrains the alignment of the CBD relative to cellulose. Figure 1 shows the modeled initial structure of the CBDcellulose complex for the 20 ns simulation. During the simulation, the RMSD values of all atoms compared to the starting structure reach~0.8 Å after several ps and are maintained at around 1 Å without large deviations (Fig. 4a).
In addition, the distance between the center of mass (COM) of the CBD and the COM of cellulose remains stable (Fig. 4b). Figure 4c shows the projection of the CBD COM on the cellulose surface during the 20 ns MD simulation, which is narrowly distributed. The results indicate that the system is equilibrated after the 20 ns simulation and can serve as the input geometry for the TI simulations.

Free energy calculations
The mutations of functionally important amino acids such as Y5, N29, Q34, Y31 and Y32 in the CBD domain to alanine have been shown experimentally to reduce the binding affinity and activity of Cel7A to crystalline cellulose [4,32,33]. To elucidate the fundamental principles governing the interaction between the CBD and cellulose, a three-step TI method was used to calculate the binding free energy differences between the WT-CBD-Cel and MU-CBD-Cel (see "Methods"). The calculations from the first and third steps are related to electrostatic effects, and those from the second step are the VDW interactions.
In the second step, soft core Lennard-Jones (LJ) potentials were used to provide smooth 1 curves that are best suited to commonly used numerical integration schemes [14]. Plots of ∂V/∂1 versus 1 for steps one and two are depicted in Fig. 5. Switching on the charge of the third hydrogen attached to C β , as done in step three, provides a significantly smaller contribution to the free energy changes for all five mutations, so the corresponding ∂V/∂1 plot is not shown here. The ∂V/∂1 values versus time for step two at 1 =0.125 are shown in Fig. S1 as a typical energy derivative plot (see the "Electronic supplementary material"), and this indicates that the systems are equilibrated during the production period. It is worth noting that in the VDW soft core alchemical process, a Y to A mutation generally displays larger ∂V/∂1 fluctuations than a Q or N to A mutation does, especially for the CBD-cellulose complex (Fig. 5). To further investigate the cause of the energy derivative fluctuations, we selected two neighboring windows with the largest ∂V/∂1 difference (13.8 kcal mol −1 ), corresponding to 1 =0.875 and 0.900 of the CBDcellulose complex Y32A VDW alchemy step. It can be seen (Fig. S2) that the side-chain χ 2 distribution is quite broad (ranging from −180°to 180) for both windows, suggesting that tyrosine ring rotation occurs in the simulations, which can cause a dramatic fluctuation in the VDW interaction between the tyrosine and a glucose ring of the cellulose chain. Furthermore, a shift of the χ 1 distribution was observed for the 1 = 0.875 window compared to that for 1 =0.900 (Fig. S2). Thus, from a structural point of view, the bulky size of tyrosine and a slightly different local environment around the mutated site are the main cause of the large ∂V/∂1 fluctuation in the Y to AVDW alchemical process.
Total free energies were obtained by summing the contributions of all three steps. Table 1 shows the calculated and experimental relative free energy differences between the WT-CBD-Cel and MU-CBD-Cel. According to the TI simulations, the calculated ΔΔG cal values for Q34A, N29A, Y5A, Y31A and Y32A are 3.42, 3.68, 5.40, 6.35 and 10.76 kcal mol −1 , respectively, suggesting that all five of these residues play roles in the binding to the cellulose surface. The corresponding standard deviations were 0.57, 0.74, 0.82, 0.59 and 0.87 kcal mol −1 , respectively.
The three aromatic residues Y5, Y31 and Y32 show larger binding free energy contributions than the residues N29 and Q34 do, which is in good agreement with the experimental results [4]. Among the five mutations, Y32A shows the largest binding free energy loss. The NMR structure of the CBD suggests that the substitution at position 32 could cause significant disruption of the flat CBD cellulose contact face due to the close proximity of the neighboring side chains, which might account for the binding importance of this residue. In general, the absolute values of the free energy changes calculated from TI simulations are higher than the experimental values extrapolated from the adsorption isotherms. This may be due to a number of factors. (1) The binding of CBD to cellulose is very complex. Lehtio and coworkers [34] showed that the primary binding site is the (110) face of cellulose, which is modeled here. However, as it is stated by Lehtio that "it cannot be ruled out that the initial, perhaps lower affinity, binding occurs in the (100) and (010) faces followed by tighter binding at the (110) plane as the hydrolysis proceeds," the binding on other faces is most likely less sensitive to the mutations. (2) The cellulose model built here has an ideal flat surface composed of polyglucan chains; this is restrained during the TI simulations to simplify the system and yield better convergence; however, in reality, the surface may twist due to molecular motions. The restraints on the C1 and C4 atoms also limit the conformational sampling of glucose units in contact with the CBD side chains. This effect was studied by removing the restraints, except for those on the terminal glucose units. We repeated the TI calculations for Q34A and Y5A. The calculated ΔΔG values decrease slightly, by 0.48 kcal mol −1 for the former and 1.08 kcal mol −1 for the latter-comparable to the uncertainties of the calculations (Table 1). This effect alone is unlikely to alter the predicted binding trend or account for the differences between experimental and computational ΔΔG values. (3) The calculated binding free energies may also depend on the cellulose binding surface size. Presumably, the smaller surface should yield looser binding. This effect is not considered in our work. (4) The orientation of CBD relative to cellulose sampled in the MD simulations is limited. The TI calculations adopted here are very time-consuming and costly, impeding us from exploring this factor. (5) Force field errors for the cellulose may also contribute to the deviations. The GLYCAM06 force field has not been widely tested for binding free energy calculations of carbohydrates binding to protein receptors. Further investigation is needed to clarify these factors and yield quantitative agreement with the experimental data.
Although free energy decomposition is only an approximate process [35][36][37], it offers useful information about the molecular mechanism of the CBD-cellulose binding. The three-step TI simulation conveniently separates the electrostatic and VDW contributions to the free energy. We can see that ΔΔG vdw for the three aromatic residue mutations Y5A, Y31A and Y32A are 4.13, 4.39 and 8.97 kcal mol −1 , respectively, corresponding to 76%, 69% and 83% of the total free energy changes, which indicates that the VDW interactions of the three residues play a more important role in the tight binding to cellulose than the electrostatic interactions. For the N29A and Q34A mutations, ΔΔG vdw accounts for 37% and 48% of the corresponding total free energy differences. Considering that tyrosine residues contribute more to the overall binding, we can conclude that VDW interactions are more important than electrostatics in the binding of CBD to cellulose.
However, if we integrate ∂V/∂1 for the electrostatics and VDW separately (Fig. 5), it becomes apparent that the electrostatic change is much larger than the VDW change in the mutation to alanine for all five residues (Fig. 5) of the CBD-cellulose complex. This result underlines the importance of electrostatics in the CBD-cellulose binding, in concert with an early theoretical study of Cel7A CBD processivity [12] as well as a more recent quantum mechanics (QM) study of a tyrosine methyl β-D-glucopyranoside (OMG) complex, which suggests that dispersive interactions are less important than the hydrogen bond between the tyrosine OH and sugar O4 [38]. This observation seems to suggest that electrostatics is more important than VDW in the binding. However, according to the thermodynamic cycle (Fig. 2), the reference state (CBD alone in solution) must also be considered. In fact, the electrostatics dominates the overall free energy change (ΔG 2 ) for the free CBD forms (Fig. 5). This electrostatics presumably arises from the change in solvation energy when mutating a polar residue (e.g., N or Y) to a nonpolar alanine. The comparable electrostatic change for a mutation in CBD and the CBD-cellulose complex largely cancels out. In another words, the solvation energy loss due to the binding of CBD to cellulose is compensated by the electrostatic interactions of CBD and cellulose. As a consequence, for ΔG 2 − ΔG 1 (the binding free energy Þwhere K mut and K wt are the partition coefficients of the mutant and wild type b The values cannot be extrapolated from the adsorption isotherms, but they should be higher than 1.74 kcal mol −1 , according to the adsorption isotherms [4] difference), the VDW contribution is larger and thus more important for binding. An early experimental isothermal titration microcalorimetry study of the binding of CBD from the β-1,4 exoglucanase Cex of Cellulomonas fimi to bacterial microcrystalline cellulose (BMCC) suggests that electrostatics does not drive the adsorption of CBDcex onto cellulose. CBDcex is larger than CBD Cel7A, but both have a flat surface composed of three aromatic residues (tyrosine or tryptophan). Our result is consistent with this experimental finding.
Our work confirms that the hydrophobic effects of the three tyrosine residues Y5, Y31 and Y32 that form the flat surface of the CBD are required for tight binding to cellulose. Experimental data show that the CBD binding free energy difference ΔG W5 − ΔG Y5 is −0.26 kcal mol −1 , suggesting that the Y5W mutant has a higher binding affinity to cellulose [32]. Our computational values are −0.24±0.55 kcal mol −1 for side-chain configuration 1 and 0.88±0.52 kcal mol −1 for configuration 2. While both configurations form a flat surface on contact with cellulose, configuration 1 is favored since it is more stable ( Table 2). The free energy decomposition suggests that while the VDW interaction becomes stronger for the Y5W mutant, presumably due to its larger hydrophobic surface, the electrostatic interaction is weaker, probably because of the loss of the polar tyrosine hydroxyl group. Compared to the wild type, the Y31W mutant exhibits decreases in the binding free energy by 0.34±0.35 and 0.62±0.40 kcal mol −1 for the two configurations. Similar to Y5W, the electrostatic and VDW changes offset each other. The Y32W mutation was not investigated, because substituting the Y32 phenol side chain by the larger indole group of tryptophan in the confined space could disturb the structure and make the surface less flat [39].
The experimental study of the impact of Humicola grisea CBHI CBD cellulose binding affinity on CBHI hydrolysis activity against cellulose has shown that the  (1) and (2) stand for free and cellulose-bound CBD states, respectively tighter the binding, the higher the catalytic activity [39].
Our work suggests that one way to increase the binding affinity is to strengthen either the VDW or the electrostatic interactions, but not to sacrifice the other.
Changes in the dynamics of side chains due to binding The protein side-chain dynamics represents a major component of protein conformational entropy [40]. Loss of side-chain conformational entropy could be an important influence that reduces a substrate's binding affinity to an enzyme. The torsional entropy of side chains is thought to provide a good approximation for the side-chain entropy. To investigate the side-chain dynamics, we monitored the changes in side-chain dihedral angles and entropy of the five CBD residues before and after binding to cellulose. Figure 6 shows scatter plots of Chi1 (χ 1 ) and Chi2 (χ 2 ) for the five residues Y5, N29, Y31, Y32 and Q34 from the 20 ns free CBD and cellulose-bound CBD simulations. Table 3 lists the averages of and fluctuations in the χ 1 and χ 2 angles of the five residues in free and bound states. In general, as expected, the χ 1 and χ 2 values of the side-chain residues in the free CBD spread over a larger area than those of the cellulose-bound CBD, which indicates that the conformations of the side chains are restrained due to the binding to cellulose. As shown in Fig. 6, after binding to cellulose, the number of the major side-chain states reduces from 5 and 3 to 1 for Tyr5 and Tyr31, respectively. The fluctuations in χ 1 and χ 2 for both Y5 and Y31 decease to ≤10°( Table 3). The considerably restricted side-chain rotations due to binding are further reflected in the notable drop in the entropic contribution of the side chain to the free energy (TΔS), with differences between the ligand-free and bound states for Y5 and Y31 of −1.71 and −1.02 kcal mol −1 , respectively (Fig. 7). Unlike Y5 and Y31, the distributions the of χ 1 and χ 2 values for Y32 are similar in the free and bound states, with only a small fluctuation of~10°( Table 3), suggesting that this side chain is confined in both states. The TΔS value between the free and bound states is only −0.11 kcal mol −1 for this residue. In comparison, the side-chain torsional entropy losses for N29 and Q34 are slightly larger: −0.75 and −0.30 kcal mol −1 , respectively (Fig. 7). It is worth noting that the entropy change considered here is very rough, and the sampling of the dihedral space with a 20 ns simulation is rather limited. However, qualitatively, it is clear that for some residues the entropy changes are relatively large and partially offset the gain in interaction energy after CBD binding.

Conclusions
In this paper, we studied the effect of a single mutation of Cel7A CBD on binding to the cellulose surface using a TI simulation method. Analyses of the computed free energy differences provided useful information for understanding the binding mechanism between CBD and cellulose. Our results suggest that residues Y5, Y31 and Y32 are more important than N29 and Q34 for the binding of CBD to cellulose. The VDW interactions between the amino acid Y5, Y31 and Y32 and the cellulose surface account for a large proportion of the binding affinity. Our calculations yielded relative binding affinities with a similar trend, but overestimated the experimental ΔΔG data. The current model yields only qualitative agreement with experiments. A more realistic cellulose model with properly weighted CBD binding modes will likely yield a more quantitative agreement, which will be explored in a further study. The changes in the binding affinities of the mutants Y5W and Y31W to cellulose were calculated. The mutant Y5W shows slightly improved affinity, which is in agreement with the experimental results. Although the VDW interaction between CBD and cellulose is enhanced for the Y to W mutations, weakening of the electrostatics largely offsets this effect. The side-chain torsional entropy decreases for Y5, N29, Y31, Y32 and Q34 due to rotamer rigidification after binding to cellulose. Fig. 7 Side-chain dihedral entropic contributions to the free energies TS (T=300 K) of the five residues in the free and cellulose-bound states. The difference between the two states for each residue is shown as a blue bar Table 3 Averages of and fluctuations in the angles χ 1 and χ 2 for the five residues in the free and bound states

Residue
Free state Bound state