Systematic Docking Study of the Carbohydrate Binding Module

Publication Date (Web): November 24, 2009 .... The impact of Trichoderma reesei Cel7A carbohydrate binding domain mutations on its binding to a cellul...
2 downloads 0 Views 1MB Size
Download PDF
J. Phys. Chem. B 2010, 114, 49–58

49

Systematic Docking Study of the Carbohydrate Binding Module Protein of Cel7A with the Cellulose Ir Crystal Model Toshifumi Yui,*,† Hirohide Shiiba,† Yuya Tsutsumi,† Sachio Hayashi,† Tatsuhiko Miyata,‡ and Fumio Hirata†,‡ Department of Applied Chemistry, Faculty of Engineering, UniVersity of Miyazaki, 1-1 Nishi, Gakuen Kibanadai, Miyazaki 889-2192, Japan, and Department of Theoretical and Computational Molecular Science, Institute for Molecular Science, 38 Nishigo-Naka, Myodaiji, Okazaki, Aichi 444-8585, Japan ReceiVed: August 26, 2009; ReVised Manuscript ReceiVed: October 27, 2009

A computer docking study has been carried out on the crystal surfaces of cellulose IR crystal models for the carbohydrate binding module (CBM) protein of the cellobiohydrolase Cel7A produced by Trichoderma reesei. Binding free energy maps between the CBM and the crystal surface were obtained by calculating the noncovalent interactions and the solvation free energy at grid points covering the area of the unit cell dimensions at the crystal surface. The potential maps obtained from grid searches of the hydrophobic (110) crystal surface exhibited two distinct potential wells. These reflected the 2-fold helical symmetry of the cellulose chain and had lower binding energies at the minimum positions than those for the hydrophilic (100) and (010) crystal surfaces. The CBM-cellulose crystal complex models derived from the minimum positions were then subjected to molecular dynamics (MD) simulation under an explicit solvent system. The (110) complex models exhibited larger affinities at the interface than the (100) and (010) ones. The CBM was more stably bound to the (110) surface when it was placed in an antiparallel orientation with respect to the cellulose fiber axis. In the solvated dynamics state, the curved (110) surface resulting from the fiber twist somewhat assisted a complementary fit with the CBM at the interface. In addition to the conventional Generalized Born (GB) method, the threedimensional reference interaction site model (3D-RISM) theory was adopted to assess a solvent effect for the solvated MD trajectories. Large exothermic values for the noncovalent interactions appeared correlated to and were mostly compensated by endothermic values for the solvation free energy. These gave total binding free energies of -13 to -28 kcal/mol. Results also suggested that the hydrogen bonding scheme was not essential for substrate specificity. Introduction Cellulose is a highly crystalline biopolymer that is the most abundant bioresource on Earth. It is generally present in the form of a single crystalline microfibril and functions as the major component of the cell walls of plants and marine algae and the extracellular material of some bacteria. Due to its abundance, cellulose is a potentially enormous source of renewable liquid fuel, namely, as bioethanol. However, its crystalline nature presents difficulty in degrading it into component glucose residues, probably through a swollen crystalline state. Trichoderma reesi is an asexually reproducing filamentous fungus and is considered to be an asexual, clonal line derived from a population of the tropical ascomycete Hypocrea jecorina.1 T. reesei produces two cellobiohydrolases (CBHs), called Cel7A and Cel6A, along with several endoglucanases. Each CBH consists of a large catalytic module and a small carbohydrate binding module (CBM) linked by an O-glycosylated linker peptide. CBHs hydrolyze crystalline cellulose from either the reducing (Cel7A) or nonreducing (Cel6A) terminal to release a cellobiose unit.2 It has been suggested that the CBM plays a significant role in the enzymatic action on crystalline cellulose by selectively binding to the surface of the crystalline substrate. CBMs are divided into 58 defined families on the basis of the * Author to whom correspondence should be addressed. Phone: +81985-7319. Fax: +81-985-7323. E-mail: [email protected]. † University of Miyazaki. ‡ Institute for Molecular Science.

amino acid sequence similarity (see http://www.cazy.org/ index.html).3 Both the CBMs of Cel7A and Cel6A belong to the family 1 so that they would be called Cel7ACBM1 and Cel6ACBM1, respectively.4 The three-dimensional (3D) solution structure of Cel7ACBM1 has been already determined by nuclear magnetic resonance (NMR) spectroscopy,5 which revealed a wedge-like irregular triple-stranded antiparallel β-sheet. Three tyrosine residues, Tyr5, Tyr31, and Tyr32, also characterize the 3D structure. These are aligned nearly colinearly along its processing direction. The aromatic side groups of the tyrosine residues constitute a hydrophobic patch on the flat surface of the CBM, which matches a similar hydrophobic surface on the native cellulose crystal. Experimentally, the CBM proteins have been observed to exclusively bind on the hydrophobic (110) plane on the Valonia cellulose crystals by transmission electron microscopy.6 This structural knowledge of CBM substrate recognition facilitates a computer docking study of CBM and its cellulose complex system. Such a study is expected to provide structural information at the atomic level. However, the predicted, essential features should be validated experimentally. Recently, Mulakala and Reilly attempted a docking study of a cellooligomer chain toward both the catalytic cleft of the catalytic module and the surfaces of the CBM.7 Of the later docking models, the AutoDock8 search based on the Lamarckian genetic algorithm predicted an additional docking position referred to as a “rough face” and involving the opposite side of the CBM to the original

10.1021/jp908249r  2010 American Chemical Society Published on Web 11/24/2009

50

J. Phys. Chem. B, Vol. 114, No. 1, 2010

hydrophobic, flat face of the three tyrosine residues. This suggested a novel binding model of the CBM with it wedged between the cellulose crystal surface and the cellulose chain detached from the surface. Nimlos et al. attempted to simulate the CBM-cellulose docking study9 according to the hypothesis that the CBM initially rode on the cellulose surface near an extracted cellulose chain.2 The molecular dynamics trajectories started with the CBM in different positions with respect to the broken cellulose chain. One trajectory illustrated three tyrosine residues of the CBM aligned with the cellulose chain, while a forth tyrosine (Tyr13) moved to develop van der Waals contact with the cellulose crystal surface. This conformation change was described as a surface-induced fit of the CBM. It was also observed in the simulations that the CBM moved randomly in both perpendicular and parallel directions relative to the cellulose chain, suggesting that the entire enzyme complex should be modeled to simulate progressive movement. This type of molecular dynamics simulation was reported by Zhong et al.10 The entire CBH complex, consisting of the CBM and the catalytic module connected by the linker module, was docked onto the intact surface of the cellulose crystal model. During the 1 ns simulation, it was observed that the orientations of the CBM and the catalytic module slightly changed, and the polypeptide backbone of the linker module displayed significant flexibility. In these studies, the docking searches to construct an enzyme-substrate complex model were carried out by a random sampling procedure partly coupled with an assumed basic configuration for the complex. Additionally, quite a few studies have been carried out on the adsorption behavior of an isolated CBM protein, which allows performance of comparative modeling and MD studies on the complex system. The present study has focused on elucidating atomic details of the very initial stage of complex formation where the CBM should be bound to an intact cellulose crystal surface involving the lattice periodicity. The exhaustive search for the distinct configuration of the CBM with respect to the cellulose crystal surface has been enabled by the systematic two-dimensional scan of the CBM over the surface area within the unit cell dimensions.

Yui et al.

Figure 1. Atom labels for a glucose residue (top) and a tyrosine residue (bottom) with the definitions of the torsion angles of primary interest. The atom sequences to define these torsion angles are: ω ) O5-C5-C6-O6, χ1 ) C-CR-Cβ-Cγ, and χ2 ) CR-Cβ-Cγ-Cδ.

Computational Methods Crystal Model Building. Figures 1a and 1b show the atomlabels and the torsion angle parameters of primary interest for a glucose residue and a side chain of a tyrosine residue, respectively. The three representative conformers of the hydroxymethyl group and their notations are also described in Figure 1a. Three types of cellulose IR crystal models with infinite dimensions were constructed, each consisting of different crystallographic faces. Figure 2 illustrates the ab projections of these cellulose crystal models. The atomic coordinates of cellooligomers and their packing arrangements are based on the high-resolution crystal structure data of cellulose IR presented by Nishiyama et al.11 The diagonal and square models, originally named by Matthews et al.,12 consist of 36 and 30 cellooligomer chains, respectively. The diagonal model exposes the hydrophilic (100) and (010) surfaces where the polar functional groups are readily accessible to solvent. Actual cellulose microfibrils are in fact mainly composed of this type of cross section with various dimensions. The square model, corresponding to a rather unrealistic microfibril, was constructed to emphasize the hydrophobic (110) surfaces with dimensions wide enough to perform the docking search without an edge effect. The third type of crystal model, referred to as the hexagonal model, had the largest dimensions and consisted of 42 cellulose chains. Narrow (110) surfaces composed of three cellooligomer chains

Figure 2. ab projections of the three types of cellulose crystal models: the square (a), diagonal (b), and hexagonal (c) models. The Miller indexes of the constituent crystal planes are noted.

were introduced as a target face for the CBM. This model represented a partially degraded diagonal crystal. The cellooligomer chain with 30 monomer units in length was commonly adopted for all three crystal models to construct the complex with the CBM for the solvated MD calculations. Shorter diagonal and square crystal models consisting of 12-mer cellooligomers were used in the grid docking searches. Program and Parameters for Minimization and Molecular Dynamics. The minimization and MD calculations were performed using the SANDER module of the AMBER 8

Carbohydrate Binding Module Protein of Cel7A

J. Phys. Chem. B, Vol. 114, No. 1, 2010 51

TABLE 1: Dimensions and Number of Molecules Constituting the Complex Models periodic box dimensions (Å)a

no. of molecules crystal model

docking face

CBM

cellooligomer (no. of residues)

square diagonal diagonal hexagonal

(110) (100) (010) (110)

1 1 1 2

30(900) 36(1080) 36(1080) 42(1260)

water 23 20 21 46

a

397-23 404 768/20 774 309/21 325 702-47 166

xb

y

z

191 192 192 197

80 64 68 94

69 83 79 103-105

a Values slightly change among the complex models with different minimum positions of the CBM. direction of the cellulose crystal model.

package,13 combined with the ff99 force field14 for a protein system. The GLYCAM04 parameter set15-17 was adopted for the cellulose crystal system. It should be noted that the parameter set has been developed without scaling 1-4 nonbonded interactions so that the 1-4 electrostatic (SCEE) and nonbonded (SCNB) scaling factors are set to unity, particularly for correctly reproducing the populations of the hydroxymethyl group orientations (ω). In the present study, the scaling factors were set to the default values, i.e., SCEE ) 1.2 and SCNB ) 2.0, according to the protein force field. Grid Docking Search. The initial structural model of Cel7ACBM1 was obtained from the Protein Data Bank (PDB). The PDB data (1CBH) comprised 27 NMR-derived conformations, and the coordinate set of the first entry was arbitrarily selected. The CBM was docked against either the (100) or (010) crystal surface of the diagonal crystal model or the (110) surface of the square model. We designed a systematic, as well as logical, docking procedure appropriate to the present complex system, as follows. First, average planes were defined for both the CBM and the cellulose crystal model. For CBM this comprised 18 aromatic ring carbon atoms belonging to the three tyrosine residues, Tyr5, Tyr31, and Tyr32, which play an essential role in recognition of the cellulose crystal surface. In the cellulose crystal model, the carbon and oxygen atoms of the surface cellooligomers were selected as reference points to establish the average plane. These included 12 pyranose ring atoms, C1, C2, C3, C4, C5, and O6, from two successive glucose residues along the cellooligomer chain in the square crystal model and eight O6 atoms pointing outward in the diagonal crystal model. Two sets of local Cartesian coordinates were defined such that each average plane of the CBM and the cellulose crystal model coincided with the x-y plane of each set of the coordinates. The molecular axis of the CBM, running from the center of the average plane to that of the side chain atoms of Tyr31, and the cellulose chain axis pointing to its reducing terminal were set as the x-axis. The positive direction of the z-axis was outward from the cellulose crystal surface. Next, the complex model system of the reference configuration was established by fitting the two coordinate systems of each component and then translating the CBM along the z-axis up to the distance that van der Waals contacts occur at the interface. The latter adjustment was manually achieved by inspecting the molecular graphics representations. When the z-axis of the CBM was rotated 0° from the reference configuration, the orientation was referred to as “parallel” (P) and “antiparallel” (AP) when the angle was 180°. In a parallel orientation, the CBM is headed for the reducing terminal of a cellulose chain, corresponding to the direction that the enzyme processes a cellulose chain from its nonreducing to reducing terminus. Finally, from the reference position, the CBM was translated along both the x (fiber) and y directions in 0.5 Å steps over the range of 11 × 9 Å for the (110) surface of the square crystal model and 11 × 7 Å for the (100) and (010) surfaces of the diagonal model. These x-y ranges covered the slightly larger areas of the unit cell

b

Corresponding to the fiber axis

dimensions of the corresponding crystal face. At each grid point, the entire complex structure was optimized under the implicit solvent system of the Generalized Born (GB) algorithm.18,19 The positions of hydrogen atoms were first optimized by 200 energy minimization steps while constraining the remaining heavy atoms of the cellooligomers and the backbone atoms of the CBM with a positional constraint force of 50 kcal/(mol Å2). The full optimization of the entire structure was accomplished over 500 steps. The binding energy between the CBM and the crystal surface was calculated by decomposing the optimized complex structure. The details to evaluate the binding energy for the complex model system will be described later. A complete survey of all grid points provided the potential energy surface with respect to the x-y grid positions. At each crystal surface, the two series of the grid surveys were carried out with the CBM oriented in either the parallel or antiparallel orientation. In the searches of the (110) surface, another four series of the grid survey with the CBM rotated by (60° from both the parallel and antiparallel orientations were carried out. As described above, the x-axis was arbitrarily defined so that the CBM of the standard orientation, i.e., z ) 0°, was not necessarily strictly parallel with the actual processing direction. To refine the z-rotation of the CBM, the representative minimum areas were further examined by a combined search of z-rotations with 2° steps in the ranges of (6° and the x-y translations. Molecular Dynamics Calculations. The minimum complex models obtained from the grid docking search and each of the isolated components, i.e., the cellulose crystal model and the CBM, were placed in a rectangular periodic box filled with TIP3P water models.20 Table 1 lists the components of the complex models and the dimensions of the periodic boxes. The water configurations were initially relaxed by 2000 minimization steps while constraining the atomic positions of the CBMcellulose crystal complexes with a force of 500 kcal/(mol Å2), followed by 5000 full minimization steps. Constant volume, NTV, dynamics runs were carried out for 300 ps, accompanied by gradually increasing the temperature up to 300 K at a constant heating rate of 1 K/ps. Finally, the production runs of NTP dynamics were performed for 2-5 ns at 1 bar and 300 K. Newton’s equations of atomic motion were integrated by a Verlet algorithm21 with a 2 fs time step combined with the SHAKE algorithm22 to fix bond stretching of the valence bonds involving hydrogen atoms. Nonbonding interactions were cut off at 10 Å, and the particle-mesh Ewald (PME) method23 was adopted to determine long-range electrostatic interactions. 3D-RISM Solvation Free Energy. The 3D reference interaction site model (3D-RISM) theory is a 3D statistical mechanical theory of molecular liquids based on the integral equation theory. The 3D-RISM theory describes the 3D distribution functions of the solvent around the solute, which can be related to thermodynamic quantities such as the solvation free energy by a spatial integration of the relevant correlation functions with a proper form of the integrands. The details of the 3D-RISM

52

J. Phys. Chem. B, Vol. 114, No. 1, 2010

Yui et al.

theory have been described elsewhere.24 Here, we shall present a brief outline of the theory. The 3D-RISM equation reads

∑ cuvν (r)*(wνγvv(r) + Fvνhνγvv(r))

hγuv(r) )

(1)

ν

hγuv(r)

cγuv(r)

where and are the total correlation functions and the direct correlation functions of the solvent site γ around the vv v solute molecules, respectively, and wvv νγ(r), hνγ(r), and Fν are the site-site intramolecular correlation function of the solvent, the site-site total correlation function of the solvent, and the number density of the solvent site ν, respectively. “u” and “v” stand for solute and solvent, respectively. The asterisk, *, represents the convolution integral. wvv νγ(r) reflects the molecular vv (r) is the correlation function of the shape of the solvent. hνγ bulk solvent and is obtained from the dielectrically consistent RISM (DRISM) theory25,26 prior to the application of eq 1. Equation 1 is to be solved for hγuv(r) and cγuv(r). To this end, we employed another relationship between hγuv(r) and cγuv(r) to couple with eq 1, which is referred to as the 3D-KH closure of the form

1 + hγuv(r) )

{

exp(dγuv(r))

for dγuv(r) e 0,

1+

for dγuv(r) > 0,

dγuv(r)

(2)

dγuv(r) ) -βuγuv(r) + hγuv(r) - cγuv(r) Here, uγuv(r) is the interaction potential energy. β is defined as β ) 1/kBT, where kB and T are the Boltzmann constant and temperature, respectively. The simultaneous solution of eqs 1 and 2 gives the 3D distribution function gγuv(r) through the following relation

gγuv(r) ) 1 + hγuv(r)

(3)

It has been shown that under eqs 1 and 2 the solvation free energy ∆Gsolv(RISM) is expressed by the following analytical formula24,27

∆Gsolv(RISM) ) -kBT

∑ ∫ dr[ Fγv

cγuv(r)

-

γ

1 1 Θ(-hγuv(r))hγuv(r)2 + hγuv(r)cγuv(r) 2 2

]

(4)

where Θ(x) is the Heaviside step function. Equation 4 indicates that the solvation free energy can be calculated by spatially integrating the correlation functions within the framework of the 3D-RISM theory. Computational Details of the Binding Free Energy. Variations in binding free energy (∆Gbind) between the CBM and the cellulose crystal surface were calculated based on the MD trajectories by evaluating the noncovalent interaction terms of the AMBER potential energy28 and the solvation free energy, as follows

∆Gbind ≈ ∆Ebind - T∆S + ∆∆Gsolv,bind

(5)

where ∆Ebind, ∆S,and ∆∆Gsolv,bind represent the differences in the noncovalent interactions of the potential energy, conformational entropy, and solvation free energy, respectively. The gas phase ∆Ebind was calculated using the following equation

∆Ebind ) Ecomp - (ECBM + Ecryst)

(6)

where Ecomp, ECBM, and Ecryst represent the AMBER potential energy of the whole complex model, the CBM, and the cellulose crystal model, respectively. Similarly, ∆∆Gsolv,bind was defined by

∆∆Gsolv,bind ) ∆Gsolv,comp - (∆Gsolv,CBM + ∆Gsolv,cryst) (7) Assuming a constant ∆S contribution for all time steps, we evaluated ∆Gbind* ) ∆Gbind + T∆S to assess the binding free energy of complex formation in solution. It is likely that the conformational entropy difference ∆S is insignificant irrespective of the binding state, i.e., ∆S ≈ 0, so that ∆Gbind* ≈ ∆Gbind. Throughout the MD trajectories, ∆Gsolv was estimated by adopting a combination of the Generalized Born (GB) method18,19 and the surface area (SA) option29 and referred to as a GB/SA energy. The latter option evaluates surfactant energy with a default value of surfactant tension at 0.005 kcal/(mol Å2). As an alternative approach to assess the solvent effect, the 3DRISM method24 (eqs 1 and 2) was employed for selected snap shots. In the 3D-RSIM calculation, the TIP3P model20 is adopted for water with a correction upon the Lennard-Jones parameters of the hydrogen atom (σ ) 0.4 Å and ε ) 0.05 kcal/mol, where σ and ε have the standard meanings concerning the LennardJones potential). The Lorentz-Berthelot mixing rule was applied to the water-cellulose and water-CBM Lennard-Jones parameters. The temperature (T ) 300 K), the number density of water (F ) 0.033309 Å-3), and the dielectric constant of water (ε ) 78.4) are employed. Furthermore, the number of spatial grid points and the size of cubic cell for solving eqs 1 and 2 are set as 400 × 400 × 400 and 200 × 200 × 200 Å3, respectively. The convergence of the 3D-RISM calculation is accelerated by using the MDIIS method.24,30 Programs for Trajectory Analysis. The MD trajectories were analyzed by using the PTRAJ module of the AMBER 8 package.13 The 3D-RISM calculations were performed using an original program developed by one of the present authors (F.H.) and his co-workers. VMD 1.8.3 software31 was used for molecular visualization and animation of the trajectory data. Results and Discussion Binding Free Energy Maps Calculated by Grid Docking Search. Binding free energy (∆Gbind-GB) maps for both the parallel and antiparallel orientations of the CBM with respect to the cellulose fiber axis were derived from the grid docking search of the CBM over the (110) crystal surface (Figure 3). Both the maps revealed two energy minima, which obviously reflect the 2-fold helical symmetry of the cellulose chain. The lowest (A) and second lowest (B) minima were separated by a distance of about 4-4.5 Å along the x (fiber) axis. They were not strictly aligned along the x-axis and appreciably deviated from each other by a few Å along the y-axis. The potential energy surface of the antiparallel map exhibited more distinct and steeper potential wells, resulting in lower energy minima than the rest of the maps. A further four (110) grid searches, adopting the CBM z-rotated by (60° from either a parallel or an antiparallel orientation, converted the energy surfaces to ones characterized by shallower and complicated topologies and involving several local minima. These results suggested that the CBM could bind at the (110) surface only when the aromatic rings of the three tyrosine residues were nearly colinearly aligned along the cellulose fiber axis. The additional x-y grid searches combined with the z-rotation over a limited area covering each of the minimum wells slightly improved the minimum energies by no more than about 1 kcal/mol and were accompanied by a z-rotation of the CBM a few degrees from the starting orientation. Grid docking searches were also carried out for the hydrophilic (100) and (010) crystal surfaces. As expected, the potential energy maps derived from the hydrophilic crystal

Carbohydrate Binding Module Protein of Cel7A

J. Phys. Chem. B, Vol. 114, No. 1, 2010 53 TABLE 3: Average GB/SA Binding Energies and rmsd Values of the Minimum Complex Models average energy (rmsd) (kcal/mol)a parallel positional constraint

A

B

A

B

(110)

none constrd none none

-59(4) -54(5) -42(7) -53(5)

-53(4) -71(5)

-68(4) -51(4) -37(5) -48(5)

-66(5) -72(5)

(100) (010) a

Figure 3. GB binding free energy (∆Gbind-GB) maps with respect to the translational positions of the CBM in Å over the (110) crystal surface. The potential energy surfaces are indicated by contour lines with an interval of 2 kcal/mol. The CBM is oriented in either the parallel (P) or antiparallel (AP) orientation with respect to the fiber axis of the cellulose crystal model.

TABLE 2: GB Binding Free Energies of the Minimum Complex Models minimum energy (kcal/mol)b docking crystal surface

orientation of CBMa

A

B

(110)

P P+60 P-60 AP AP+60 AP-60 P AP P AP

-52(-52) -53c -52c -56(-57) -52c -50c -38 -38 -42 -37

-51(-52)

(100) (010)

-53(-53)

-36 -35

a P, parallel; AP, antiparallel. b Energy value of the minimum complex model obtained from an additional z-rotation search is given in parentheses. c Value of the minimum complex model among the multiple minima.

surfaces gave shallower less distinctive minima. The (110) maps for the (60° z-rotated CBM and the (100) and (010) maps can be obtained from the Supporting Information. Table 2 lists the binding free energy of the minimum complexes obtained from the grid searches for the three cellulose crystal surfaces. The CBM formed more stable complexes on the (110) crystal surface than on the two hydrophilic surfaces, and the global minimum was found in the (110) complex model with the antiparallel CBM. On inspection of the four minimum complex structures, both the parallel and antiparallel A-minimum models were found to have the two aromatic rings of the Tyr31 and Tyr32 residues simultaneously partially stacked on the two pyranose rings separated by the fiber repeat in the same cellulose chain. It was also found that the aromatic rings of Tyr31 and Tyr32 were virtually stacked off the pyranose rings in the parallel Bminimum model, as was Tyr32 in the antiprallel B-minimum model. The two antiparallel complex models were preferred over the two parallel counterparts by significant differences of 2-4

antiparallel

docking crystal surface

Average values from the final 1 ns trajectory.

kcal/mol. Interestingly, the orientation with CBM directed toward the nonreducing terminal of the cellulose chain was consistent with the proposed processing direction of Cel7A.32 Assuming that the CBM makes first contact with the cellulose crystal surface in advance of the catalytic module, the directional specificity of the CBM would increase the rate of productive recognition of the entire enzyme. In comparison with the minimum complex models formed on the hydrophobic (110) crystal surface, less stable complexes formed on the hydrophilic (100) and (010) crystal surfaces. The lowest-energy minima reached no less than about -40 kcal/mol, and the obviously second lowest minimum could not even be detected from the parallel maps. The eight minimum complex models, four from the (110) maps and four from the A minima of both the (100) and (010) maps, were selected for the molecular dynamics calculations with an explicit solvent system. Solvated Molecular Dynamics Simulations of the Minimum Complex Models. The square crystal model consisting of 30-mer cellulose chains was adopted for construction of the (110) complex models, where the CBM was placed on either of the two minimum positions. As was reported previously for solvated native cellulose crystal models,12,33,34 the model developed a fiber twist resulting in a curved (110) surface. The CBM remained to bind with the (110) surface in the four complex models throughout the 2 ns simulation in spite of the significant fiber twist. To examine the effect of the fiber twist on complex formation, positional constrained solvated MD calculations were carried out for the (110) complex models. A constraint force of 100 kcal/(mol Å2) was imposed on the atomic positions of the two inner chain layers in the course of the MD simulation, and this held the surface layers flat while allowing local motion of the surface atoms. Similarly, the hydrophilic (100) and (010) complex models were constructed by placing the CBM at the minimum positions on the bc and ac crystal surfaces of the 30-mer diagonal model, respectively; only the A-minimum models were examined. The trajectories of the GB/ SA binding free energies (∆Gbind-GB/SA) for these 12 complex models are included in the Supporting Information, and Table 3 lists the average binding free energies and their root-meansquare deviations (rmsd) calculated from the final 1 ns MD trajectories. The stability ranking of the complex models in the unconstrained dynamics state was similar to that of the minimized models listed in Table 2. Among the eight unconstrained models, the two antiparallel (110) and the parallel and antiparallel (010) complex models lowered their binding free energies the most (by >10 kcal/mol) from the starting minimized values. As a result, the antiparallel CBM orientation was obviously preferred among the four (110) complex models. The antiparallel complex models at the A and B minima showed the two lowest average binding energies of -68 and -66 kcal/ mol, respectively. Their parallel counterparts gave values of -59 and -53 kcal/mol. Since the curved (110) surface acts like a

54

J. Phys. Chem. B, Vol. 114, No. 1, 2010

chiral interface for the docking face of the CBM, it could provide greater affinity to antiparallel CBMs over parallel ones. The four minimized complex models involving the hydrophilic interfaces gave comparable binding energies (-42 to -37 kcal/ mol, Table 2). In comparison, the two (010) complex models were more effectively stabilized by the solvated MD calculation and gave average binding energies of -53 and -48 kcal/mol (Table 3) for the parallel and antiparallel complex models, respectively. Both the hydrophilic crystal surfaces have variable and complicated topologies due to rotatable substituents such as hydroxyl and hydroxymethyl groups. As a result of dynamic motion, this may have allowed the (010) surface to fit better with the docking face of the CBM. In the case of the constrained (110) complex models, both the parallel and antiparallel B-minimum models were significantly stabilized on the flat, constrained surface, while the two A-minimum models were destabilized. Molecular visualizations of the MD trajectories for all four constrained complex models consistently displayed significant z-rotations of about 14-16° for the CBM as a result of dynamic motion. The direction of the rotation was found to be clockwise when looking down on the (110) face from the outside. In contrast, the z-rotations of the CBM were diverse in the unconstrained complex models. The CBMs of the parallel and antiparallel A-minimum models rotated to be anticlockwise, and the two B-minimum models exhibited clockwise rotations. Although the determination of z-rotation was less definite due to the curved (110) plane, it appeared to be around 10°. The largest rotation was observed in the CBM of the antiparallel A-minimum model. After docking of the whole Cel7A system, spatial movement of the actual CBM would be constrained by the catalytic module through the linker peptide. Obviously, the constraint operates on the C-terminal of the CBM in an anisotropic fashion. The z-rotation to an undesired direction would cause the CBM to quickly lose affinity to the crystal surface. The docking grid searches for the square model successfully assigned the hydrophobic (110) surface to a binding face of the CBM. However, the actual crystalline microfibrils of native cellulose have mostly been observed with diagonal-type constituent crystal faces. An intact microfibril involves the (110) edges consisting of a single cellulose chain, and they can be expanded as an enzymatic hydrolysis proceeds. To construct a more realistic (110) complex model, a hexagonal model with two equivalent (110) faces (Face 1 and 2) of three cellulose chains was examined as a partially hydrolyzed microfibril. The CBM was then placed on both of the (110) surfaces, and the two equivalent interfacial structures were examined simultaneously. Figure 4 compares the 5 ns trajectories of the GB/SA binding free energy calculated for the four hexagonal complex models. Both of the trajectories obtained from each interface displayed substantially similar behaviors in all of the four complex models, and consequently only those of the one interface (i.e., Face 1) are displayed in Figure 4. Compared with the square complex models, there was remarkable stabilization observed in the trajectory of the hexagonal antiparallel Aminimum model, with binding energies ranging from -100 to -70 kcal/mol. The affinity of the CBM appears to be gradually loosened after the 3 ns simulation. The rest of the hexagonal complex models exhibited energies comparable to their square model counterparts. Our previous MD studies of solvated native cellulose crystal models reported that the extent of fiber twist depended on the crystal model dimensions.33,34 Therefore, the hexagonal model had more moderate fiber twist than the square model, which may have happened to be the surface with a better

Yui et al.

Figure 4. Trajectories of the GB/SA binding free energy calculated for the four complex models, the two minimum, A and B, (110) complex models having the CBM in parallel (P) and antiparallel (AP) orientations. Only the data obtained from the Face 1 interface are shown.

topological fit to the CBM residing at the antiparallel Aminimum position. This also suggested that the affinity of the CBM depends on the size, or number of cellulose chains composing the surface, of the (110) surface with respect to the horizontal direction. Slight clockwise z-rotations were observed only in the two parallel models. MD studies of the hexagonal complex models reasonably characterized the substrate specificity of the CBM in terms of its orientation preference and crystal binding surface. However, the predicted values of the binding free energy were not realistic in light of the fact that Cel7ACBM1 reversibly adsorbs on the cellulose crystal substrate.35 These absolute values, up to 100 kcal/mol in the antiparallel A-minimum model, seemed too large to allow reversible adsorption. Consequently, the MD trajectories and solvation free energy were re-evaluated using the 3D-RISM (or RISM) solvation model. This alternative approach had more sophisticated background theory but was more costly than the GB/SA model. Figure 5 shows the GB/SA and 3D-RISM binding free energies of the parallel and antiparallel A-minimum models at the Face 1 interface. Due to the large computational cost of the 3D-RISM calculations, only 1 ns interval snap shots were extracted from the MD trajectories. Compared with the GB/SA binding free energies, the 3D-RISM models gave more realistic absolute values of about 13-28 kcal/mol. Both types of the binding free energies shared the same set of AMBER

Carbohydrate Binding Module Protein of Cel7A

Figure 5. GB/SA (∆Gbind-GB/SA) and 3D-RISM binding free energies (∆Gbind-RISM) calculated for the parallel (O) and antiparallel (b) A-minimum models at five 1 ns interval snap shots.

Figure 6. Differences in the noncovalent interactions of the potential energy (∆Ebind) and the 3D-RISM solvation free energy (∆∆Gsolv-RISM) calculated for the parallel (O) and antiparallel (b) A-minimum complex models at the five 1 ns interval snap shots.

potential energies at each snap shot as described in eq 5. The 3D-RISM model clearly estimated larger positive solvation free energies than the GB/SA model. Additionally, differences in the 3D-RISM energy between the parallel and antiparallel complex models became indistinct, and the parallel model gives lower energies at three of the five snap shots. Considering that the binding free energy trajectory can vary by up to 20 kcal/ mol as estimated from the trajectories in Figure 4, it should be noted that these selected snap shots do not necessarily represent the whole trajectory. Besides, there is no theoretical reason that a pair of the solvation energies predicted by the two methods should be linearly related to each other. Figure 6 displays the noncovalent interactions of the potential energy and the 3DRISM solvation free energy for interactions, each partitioned from the 3D-RISM binding free energy. The potential energy showed large exothermic values throughout the trajectories, which were mostly compensated by the endothermic solvation free energy, so that the two energy trajectories appeared to be correlated. During the initial 3 ns, the CBM of the antiparallel complex model was more strongly bound to the crystal surface. After about 3 ns, it gradually loosened its noncovalent interaction with the surface, and the solvation free energy began to dominate. The energy graphs of the parallel complex model displayed contrasting trends to the antiparallel complex model. Unfortunately, we were not able to propose a reasonable interpretation for the significant differences observed in the solvation free energy between the GB/SA and 3D-RISM methods. Further analyses based on the systematic calculations would be required for this purpose. However, the 3D-RISM approach did successfully predict binding free energies of a reasonable amount. This was achieved without introducing any adjustable parameter and by largely offsetting the AMBER potential energy for the studied systems. Default use of the GB/ SA method tends to underestimate the solvation effect. Our previous studies using the GB/SA method on cellulose crystal models predicted a realistic lattice energy difference of about 3 kcal/mol between the cellulose IR and Iβ models.34 In that study,

J. Phys. Chem. B, Vol. 114, No. 1, 2010 55 since the differential energy was derived between the same molecular systems of the different allomorphs, it was quite likely that any implicit error was canceled out. In contrast, the present docking study included two molecular components with largely different cavity volumes. Consequently, we wish to remark on the advantages and disadvantages of the 3D-RISM theory and the continuum solvent model as represented by the GB model. The greatest advantage of the continuum solvent model is its low computational cost; in a standard continuum solvent model, the computational load is in the order of O(N2), where N represents the number of atoms belonging to the solute, and is the same order as those for calculating noncovalent interaction energies such as Lennard-Jones and/or coulomb interactions inside the solute. However, there are at least two disadvantages for a continuum solvent model: the complete omission of the solvation structure and the introduction of additional adjustable parameters that are mainly related to the cavity volume in the solvent or the dielectric boundary in which the solute molecules are accommodated. Although the cavity volume or the dielectric boundary is empirically determined on the basis of the van der Waals radius, the above choice is somewhat theoretically arbitrary. Furthermore, especially in the case of the GB model, it has been reported that the choice of the effective Born radius is crucial,36 where it strongly depends on the choice of dielectric boundary. In terms of practical application of the GB model, for example to a polypeptide, it has been highlighted that GB models tend to overestimate the stability of salt-bridging structures.37 This would be due to an insufficient electrostatic screening in the GB model. As for the 3D-RISM theory, the computational load is much heavier, but in spite of this it provides a number of advantages over the continuum solvent model. Since the 3D-RISM theory is a statistical mechanical theory based on 3D distribution functions, we can directly obtain an image of the 3D solvation structure around the solute. The 3D-RISM is the theory for molecular liquids, and it assumes the molecular model of the solvent at the classical level, such as SPC/E,38 TIP3P,20 and so on. In this sense, the 3D-RISM theory can be categorized with modeling with explicit solvent molecules. An outstanding feature of the 3D-RISM theory is that it provides the solvation free energy in a straightforward manner by simply applying eq 4. In other models with explicit solvents, the free energy calculation requires an elaborate procedure involving free energy perturbation, thermodynamic integration, and so on. We note that eq 4 is originally based on the thermodynamic integration, and integration over the “coupling parameter” or the interaction-switching parameter has been completed.24,27,39 As a consequence, eq 4 is an exact expression of the solvation free energy under the assumption of eqs 1 and 2. Another point is that the 3D-RISM theory does not include additional empirical parameters, except for the dielectric constant and the interaction parameters such as Lennard-Jones parameters or partial charges. The DRISM theory requires the dielectric constant when calculating hvv νγ(r), which is part of the input data for the 3D-RISM calculation (see eq 1). The 3D-RISM theory can provide quite accurate results on solvation structures around large solute molecules including proteins and thermodynamics of the solutions.24,40-44 Dynamics Structural Features of the Antiparallel AMinimum Complex Model. Averages and distributions of the representative structure parameters of the antiparallel Aminimum model were obtained from the final 3 ns MD trajectory of the hexagonal complex model with the two CBMs. In contrast to the binding free energy behavior, some of the structural features were found to be quiet diverse between Face 1 and 2.

56

J. Phys. Chem. B, Vol. 114, No. 1, 2010

Figure 7. Distribution profiles of the side group orientations of Tyr5 (χ1 and χ2) of the CBM residing on the Face 1 and 2 interfaces on the (110) surfaces of the antiparallel A-minimum complex model (solid line) and of the isolated CBM (broken line).

Distribution profiles for the side group χ1 and χ2 angles of Tyr5 were compared for the CBMs at the Face 1 and 2 interfaces of the complex model and the isolated CBM (Figure 7). The profiles of the isolated CBM were considerably broader and even separated into two peaks, indicating that the side group of Tyr5 underwent a conformational conversion between the two states. While the χ1 peaks appeared around either the trans or gauche positions in the isolated CBM, the two states of χ2 differed by 180° and were obviously stereochemically equivalent. In the complex models, the rotational motions of the side group of Tyr5 gave single peaks and were restricted at the interfaces. Peak widths indicated that the motions were restricted more at Face 1 than at Face 2. The Tyr31 and Tyr32 peaks appeared essentially as a single profile with similar peak width between the isolated CBM and the complex models (data not included in the figure). Figure 8 compares the combined distribution profiles for the hydroxymethyl group orientation (ω) of the cellulose chains of the complex model and of the isolated cellulose crystal. Each profile was synthesized from the ω angles of 18 glucose residues from the center of the (110) face, namely, those present underneath and at the periphery of the CBM in the complex model. The starting orientation of the hydroxymethyl group was set at the tg position (ω ≈ 180°) according to the crystal structure.11 During the MD simulations, some of the hydroxymethyl groups rotated into the more stable gt (ω ≈ 60°) conformation under aqueous conditions. Interestingly, the hydroxymethyl groups belonging to the Face 1 residues more frequently changed into the gt conformation in both the complex models and the isolated crystal model. A possible interpretation of this nonequivalency would be that the fiber twist led to variation in the local structure between the two faces. This may account for the differences in the side group conformations of Tyr5 between the two interfaces (Figure 7). When comparing the isolated cellulose crystal and the complex models, the tg conformation population was fairly reduced at both the faces of the latter model, probably as a result of interactions with the CBM. Hydroxymethyl groups in the gt or gg conformations prevent the cellulose chains from forming the stable interchain hydrogen bond on the (110) face. The surface chains may be readily swollen to depart from the crystal surface. A longer simulation is required to confirm this speculation. The

Yui et al.

Figure 8. Combined distribution profiles of the hydroxymethyl group orientation (ω) for the 16 glucose residues belonging to the Face 1 and 2 interfaces on the (110) surfaces of the antiparallel A-minimum complex models (solid line) and isolated cellulose crystal (broken line).

use of the default 1-4 scaling parameters may have degraded the accuracy of the distribution profiles of the hydroxymethyl orientations in Figure 8. The profiles of the complex model should be discussed only in comparison with those of the isolated crystal model as a control system. We observed that, while the hydroxymethyl group of the solvated glucose molecule displayed more frequently transitions between the gg and gt conformers in the MD calculation with the default 1-4 scaling parameters than that without the 1-4 scaling, the extent of fiber twist observed for the diagonal crystal model was not seriously affected by introducing the 1-4 scaling (unpublished data). Table 4 lists the hydrogen bonds detected at the interfaces between the CBM and the cellulose crystal model, and the glucose residue position labels from this table are depicted in Figure 9. Only those with an occupation time of >10% are listed in the table. Despite the similarity in the binding free energy between the two interfaces, rather diverse hydrogen bonding schemes were observed at the Face 1 interface. The residues Gln7, Leu11, Ser14, and Gln34 commonly participated in the hydrogen bonding schemes at both the interfaces. Gln7 in particular involves the largest number of hydrogen bonds. Among the five additional amino acid residues listed in the Face 1 entry, stable hydrogen bonds having an occupation time of about 50% were observed for His4, Tyr5, and Asn29. The diversity of the hydrogen bonding scheme at the two interfaces was in reasonable accord with experimental results suggesting an absence of a strong hydrogen bond.45 The Gln7 side chain flanks the sequence of the three aromatic rings and lies in parallel to the hydrophilic groove present between the neighboring cellulose chains, which may play an important role in its hydrogen bond formation. Conclusions The present study successfully characterized the substrate specificity of the CBM to the cellulose crystal surfaces. One can readily infer the binding preference of the CBM to the hydrophobic (110) surface based on its docking face consisting of three tyrosine residues. The MD simulation allowed quantitative assessment of the binding interactions at the interface,

Carbohydrate Binding Module Protein of Cel7A

J. Phys. Chem. B, Vol. 114, No. 1, 2010 57

TABLE 4: Representative Hydrogen Bonding Schemea at the Interface between the CBM and the Cellulose Crystal Surface of the Antiparallel A-Minimum Complex Model CBM residue Gln2 Ser3 His4 Tyr5 Gln7

Ile11 Ser14 Asn29 Gln34 Gln7 Ile11 Ser14 Gln34 a

cellulose crystal surface

average length

occupation time

atom

residue

atom

(nm)

(%)

Nε2-H Oγ3-H O O Oε1 Oε1 Nε2-H N-H Hε2-H O O O Oγ-H Nδ2-H Nδ2-H Nδ2-H Nε2-H

c15 c14 c15 c14 b15 c16 b15 b14 c16 c17 c18 b16 c17 a14 b16 a15 b14

Face 1 O2 O5 O2 O5 H-O2 H-O6 O3 O6 O6 H-O3 H-O6 H-O6 O6 O3 O2 O6 O2

0.304 0.301 0.275 0.278 0.273 0.265 0.308 0.308 0.306 0.291 0.293 0.290 0.287 0.295 0.301 0.310 0.309

24 17 57 52 95 84 62 35 18 60 57 24 15 45 45 15 15

a15 b15 a15 a15 a15 c14

Face 2 O6 O2 H-O6 H-O2 H-O3 O6

0.271 0.278 0.288 0.278 0.304 0.304

79 69 10 68 82 87

Oε1 Oε1 O O N Nε2-H

Those with an occupation time >10%.

of the Cel7A.32 This specificity may partly assist Cel7A in achieving a productive recognition at the catalytic module, which could account for the high efficiency of the enzyme. However, during real enzyme action, it is probable that the CBM is randomly adsorbed on the various positions including at regions other than the minima. In addition, the present 3D-RISM calculations suggested solvent effects contributed considerably to complex formation. To confirm unambiguously the directional specificity, the partition function should be compared by performing the solvated MD calculations at grid points all over the crystal surface. Obviously, this type of MD search is costly and requires use of extended ensemble models. Assuming that the features obtained from the minimum complex models represent those of all the possible configurations of the complex, then two additional features were suggested by the results of this study. These were that complementary recognition of the CBM may be partly developed by the right-handed fiber twist of the cellulose crystal and that the hydrogen bonding scheme may not play an essential role in the substrate specificity. Only Gln7 may be likely to form the specific hydrogen bonds with the cellulose chains. Acknowledgment. This work was partly supported by the Strategic Research Fund 2006 from the University of Miyazaki and the Strategic Research Promotion Fund 2009-2014 from the Ministry of Education, Culture, Sport, Science and Technology, Japan. Supporting Information Available: The GB binding free energy maps for the (110) crystal surface with the (60° z-rotated CBM and those for the (100) and (010) crystal surfaces and the trajectories of the GB/SA binding free energies for the 12 complex models. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes

Figure 9. Labels and their corresponding positions of the glucose residues on the (110) crystal surface.

suggesting greater affinity of the CBM toward the (110) surface than the other two hydrophilic surfaces. We also proposed two novel features facilitating CBM substrate recognition toward the (110) surface: the two distinguished minima in the binding potential energy surface with the unit cell dimensions, and the directional specificity of the CBM at the minima with respect to the fiber axis. The two minima, separated by about 4-4.5 Å along the fiber axis, were related by the 2-fold helical symmetry of the cellulose chain. Recently, Bu et al. reported similar potential energy surfaces for the CBM-cellulose Iβ crystal model complex system by adopting a coarse-grained cellulose model in which a glucose residue was reduced into three coarsegrained beads for approximation.46 In the study, the potential energy surface exhibited alternative appearances of deeper and shallower minima at about every 1 nm along the fiber axis. The predicted directional specificity of the CBM at the minimum positions was consistent with the observed processing direction

(1) Kuhls, K.; Lieckfeldt, E.; Samuels, G. J.; Kovacs, W.; Meyer, W.; Petrini, O.; Gams, W.; Bo¨rner, T.; Kubicek, C. P. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 7755–7760. (2) Teeri, T.; Reinikainen, T.; Ruohonen, L.; Jones, T. A.; Knowles, J. K. C. J. Biotechonol. 1992, 24, 169–176. (3) Cantarel, B. L.; Coutinho, P. M.; Rancurel, C.; Bernard, T.; Lombard, V.; Henrissat, B. Nucleic Acids Res. 2009; D233-238. (4) Boraston, A. B.; Bolam, D. N.; Gilbert, H. J.; Davies, G. J. Biochem. J. 2004, 382, 769–781. (5) Kaulis, P. J.; Clore, G. M.; Nilges, M.; Jones, T. A.; Patterson, G.; Knowles, J.; Gronenborn, A. M. Biochemistry 1989, 28, 7241–7257. (6) Lehtio¨, J.; Sugiyama, J.; Gustavsson, M.; Fransson, L.; Linder, M.; Teeri, T. Proc. Natl. Acad. Sci. 2003, 100 (2), 484–489. (7) Mulakala, C.; Reilly, P. J. Proteins 2005, 60, 598–606. (8) Morris, G. M.; Goodsell, D. S.; Halliday, R. S.; Huey, R.; Hart, W. E.; Belew, B. K.; Olson, A. J. J. Comput. Chem. 1998, 19, 1639–1662. (9) Nimlos, M. R.; Matthews, J. F.; Crowley, M. F.; Walker, R. C.; Chukkapalli, G.; Brady, J. W.; Adney, W. S.; Cleary, J. M.; Zhong, L.; Himmel, M. E. Protein Eng., Des. Sel. 2007, 20, 179–187. (10) Zhong, L.; Matthews, J. F.; Crowley, M. F.; Rignall, T.; Talo´n, C.; Nimlos, M. R.; Brooks III, C. L.; Himmel, M. E.; Brady, J. W. Cellulose 2008, 15, 261–273. (11) Nishiyama, Y.; Sugiyama, J.; Chanzy, H.; Langan, P. J. Am. Chem. Soc. 2003, 125, 14300–14306. (12) 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. (13) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Merz, K. M.; Wang, B.; Pearlman, D. A.; Crowley, M.; Brozell, S.; Tsui, V.; Gohlke, H.; Mongan, J.; Homak, V.; Cui, G.; Beroza, P.; Schafmeister, C.; Caldwell, J. W.; Ross, W. S.; Kollman, P. AMBER 8; University of California: San Francisco, 2004. (14) Wang, J.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049–1074. (15) Kirchner, K. N.; Woods, R. J. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 10541–10545.

58

J. Phys. Chem. B, Vol. 114, No. 1, 2010

(16) Basma, M.; Sundara, S.; Calgan, D.; Venali, T.; Woods, R. J. J. Comput. Chem. 2001, 22, 1125–1137. (17) Kirchner, K. N.; Woods, R. J. J. Phys. Chem. A 2001, 105, 4150– 4155. (18) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. Chem. Phys. Lett. 1995, 246, 122–129. (19) Weiser, J.; Shenkin, P. S.; Still, W. C. J. Comput. Chem. 1999, 20, 217–230. (20) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926–935. (21) Verlet, L. Phys. ReV. 1967, 159, 98–103. (22) Ryckaert, J. P.; Cicotti, G.; Brendsen, H. J. C. J. Compt. Phys. 1977, 23, 327–341. (23) Essmann, U.; Pereta, L.; Berkowitz, M. L.; Darden, T. A.; Lee, H.; Pedersen, L. G. J. Comput. Phys. 1995, 14, 33–38. (24) Hirata, F. Molecular Theory of SolVation; Kluwer: Dordrecht, 2003. (25) Perkyns, J. S.; Pettitt, B. M. Chem. Phys. Lett. 1992, 190, 626– 630. (26) Perkyns, J.; Pettitt, B. M. J. Chem. Phys. 1992, 97, 7656–7666. (27) Kovalenko, A.; Hirata, F. J. Chem. Phys. 1999, 110, 10095–10112. (28) Case, D. A.; Cheatham, T. E., III; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M., Jr.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. Comput. Chem. 2005, 26, 1668–1688. (29) Weiser, J.; Shenkin, P. S.; Still, W. C. J. Comput. Chem. 1999, 20, 217–230. (30) Kovalenko, A.; Ten-no, S.; Hirata, F. J. Comput. Chem. 1999, 20, 928–936. (31) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33–38.

Yui et al. (32) Imai, T.; Boisset, C.; Samejima, M.; Igarashi, K.; Sugiyama, J. FEBS Lett. 1998, 432, 113–116. (33) Yui, T.; Nishimura, S.; Akiba, S.; Hayashi, S. Carbohydr. Res. 2006, 341, 2521–2530. (34) Yui, T.; Hayashi, S. Biomacromolecules 2007, 8, 817–824. (35) Markus, L.; Reeri, T. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 12251– 12255. (36) Onufriev, A.; Case, D. A.; Bashford, D. J. Comput. Chem. 2002, 23, 1297–1304. (37) Zhou, R. Proteins 2003, 53, 148–161. (38) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269–6271. (39) Singer, S. J.; Chandler, D. Mol. Phys. 1985, 55, 621–625. (40) Imai, T.; Kovalenko, A.; Hirata, F. Chem. Phys. Lett. 2004, 395, 1–6. (41) Imai, T.; Hiraoka, R.; Kovalenko, A.; Hirata, F. J. Am. Chem. Soc. 2005, 127, 15334–15335. (42) Yoshida, N.; Phongphanphanee, S.; Maruyama, Y.; Imai, T.; Hirata, F. J. Am. Chem. Soc. 2006, 128, 12042–12043. (43) Phongphanphanee, S.; Yoshida, N.; Hirata, F. J. Am. Chem. Soc. 2008, 130, 1540–1541. (44) Kiyota, Y.; Hiraoka, R.; Yoshida, N.; Maruyama, Y.; Imai, T.; Hirata, F. J. Am. Chem. Soc. 2009, 131, 3852–3853. (45) Boraston, A. B.; Bolam, D. N.; Gilbert, H. J.; Davies, G. J. Biochem. J. 2004, 382, 769–781. (46) Bu, L.; Beckham, G. T.; Crowley, M. F.; Chang, C. H.; Matthews, J. M.; Bomble, Y. J.; Adney, W. S.; Himmel, M. E.; Nimlos, M. R. J. Phys. Chem. B 2009, 113, 10994–11002.

JP908249R