Substrate Distortion in the Michaelis Complex of Bacillus 1,3–1,4- -Glucanase INSIGHT FROM FIRST PRINCIPLES MOLECULAR DYNAMICS SIMULATIONS*

The structure and dynamics of the enzyme-substrate complex of Bacillus 1,3–1,4-glucanase, one of the most active glycoside hydrolases, is investigated by means of Car-Parrinello molecular dynamics simulations (CPMD) combinedwith force fieldmolecular dynamics (QM/MM CPMD). It is found that the substrate sugar ring located at the 1 subsite adopts a distorted S3 skew-boat conformation upon binding to the enzyme. With respect to the undistorted C1 chair conformation, the S3 skew-boat conformation is characterized by: (a) an increase of charge at the anomeric carbon (C1), (b) an increase of the distance between C1 and the leaving group, and (c) a decrease of the intraring O5-C1 distance. Therefore, our results clearly show that the distorted conformation resembles both structurally and electronically the transition state of the reaction in which the substrate acquires oxocarbenium ion character, and the glycosidic bond is partially broken. Together with analysis of the substrate conformational dynamics, it is concluded that the main determinants of substrate distortion have a structural origin. To fit into the binding pocket, it is necessary that the aglycon leaving group is oriented toward the region, and the skew-boat conformation naturally fulfills this premise. Only when the aglycon is removed from the calculation the substrate recovers the all-chair conformation, in agreement with the recent determination of the enzyme product structure. The QM/MM protocol developed here is able to predict the conformational distortion of substrate binding in glycoside hydrolases because it accounts for polarization and charge reorganization at the 1 sugar ring. It thus provides a powerful tool to model E S complexes for which experimental information is not yet available.

Glycoside hydrolases (GHs) 2 are the enzymes responsible for the hydrolysis of glycosidic bonds and play important biological functions such as glycan processing in glycoproteins, remodeling the cell walls, and polysaccharide modification and degradation. The reaction mechanism, a classical textbook example of enzymatic reaction, has attracted much interest because genetically inherited disorders of glycoside hydrolysis often occur and because inhibitors of these enzymes can act as new therapeutic agents for the treatment of viral infections (1,2). Despite the large number of GHs known, classified into more than 90 families (1), the catalytic mechanism is similar. They typically operate by means of acid/base catalysis with retention or inversion of the anomeric configuration, although a different mechanism has recently been proposed for the GH family 4 (3). The acid/base reaction is assisted by two essential residues: a proton donor and a nucleophile or general base residue (4). Inverting enzymes operate by a single nucleophilic substitution, whereas retaining glycosidases follow a double displacement mechanism via formation and hydrolysis of a covalent glycosyl-enzyme intermediate. Both steps involve oxocarbenium ion-like transition states (Fig. 1). A current issue in the understanding of GH mechanisms is the conformational itinerary that the substrate follows during the reaction (5,6), in which substrate distortion is induced upon binding to the enzyme to reach a transition state with sp 2 geometry at the anomeric carbon.
Bacterial 1,3-1,4-␤-glucanases are highly active retaining endoglycosidases (7) belonging to family 16 retaining glycoside hydrolases (7,8). These enzymes act on linear ␤-glucans containing ␤-1,3 and ␤-1,4 linkages such as cereal ␤-glucans and lichenan, with a strict cleavage specificity for ␤-1,4 glycosidic bonds on 3-O-substituted glucosyl residues. Two Glu residues act as nucleophile and general acid/base catalyst, respectively (9). The Michaelis complex of neither 1,3-1,4-␤-glucanases nor another member of the GH 16 enzymes has yet to be characterized. However, during the last decade a growing number of crystallographic studies on retaining ␤-glycoside hydrolases have shown that the substrate binds to the enzyme in a distorted conformation (5). In particular, the saccharide unit binding at subsite 3 Ϫ1 is found to adopt a boat ( 1,4 B) or skew-boat ( 1 S 3 or 1 S 5 ) type conformation instead of the relaxed 4 C 1 chair conformation (Fig. 2a). Ring-distorted conformations have been observed in the complex of endoglucanase I from Fusarium oxysporum with a non-hydrolyzable inhibitor (10), as well as cellulase Cel5A from Bacillus agaradhaerens (11), chitobiase from Serratia marcescens (12), chitinases (13,14), ␤-galactosidase (15), and ␤-mannanase. Similar distorted structures have been characterized for inverting glycosidases and carbohydrate-bound biological receptors (16 -19, 21-23). Recent crystallographic studies on Thermobifida fusca endoglucanase Cel6A show a significant decrease in activity upon mutation of the closest Tyr residue, a feature, which has been attributed to a conformation change of the substrate in subsite Ϫ1 from skew-boat to chair (24,25). It is also worth noting that a 4 C 1 -2 S 0 equilibrium has been recently detected in the binding of glycosylaminoglycans to polypep-tides by NMR techniques (26). This kind of saccharide ring distortion has favorable mechanistic consequences in glycoside hydrolysis. As shown in Fig. 2b, the 1 S 3 distortion places the glycosidic oxygen near the acid/base residue. It also reduces the steric interaction between the hydrogen at the anomeric carbon and the nucleophile, and places the aglycon (i.e. the leaving group) in a pseudo-axial position that facilitates nucleophile attack on the anomeric carbon (27). These distortions in the Michaelis complex are therefore in the pathway to reach the transition state of the reaction.
The fact that all ␤-glycoside hydrolase enzyme-substrate complexes so far characterized by x-ray crystallography bear a distorted substrate suggests that substrate distortion is a general feature of ␤-glycoside hydrolases. However, in most cases these distortions are encountered in inhibitor-enzyme complexes (i.e. modified forms of the substrate) or in complexes with inactive enzyme mutants, which also raises the question whether the inhibitor or the mutant could affect the substrate conformation (28).
Theoretical calculations can be very helpful to solve these issues, as they can be performed directly on the "native" substrate-bound enzyme. Classical molecular dynamics simulations demonstrated that the boat conformation at subsite Ϫ1 is critical in the mechanism of family 18 chitinases (29). Recent studies confirmed that the Ϫ1 sugar moiety of the substrate in cellulase Cel6A from Trichoderma reesei, an inverting glycosidase, adopts a skew-boat conformation ( 2 S 0 ) (30). Similarly, modeling studies of ␤-galactosidase and xylanases provided evidence of substrate distortion (15,31). All these studies rely on parametrized expressions (force fields) to describe the interaction among atoms, and thus the interplay of electronic/structural factors on the substrate conformation cannot be analyzed.
To overcome these limitations, we have taken a step forward in accuracy and predictive power by using first principles methods, thus taking into account electronic effects and charge rearrangements in the active site. In the framework of our structure/function studies of bacterial 1,3-1,4-␤-glucanases (7,(32)(33)(34), we investigate here the conformation of the substrate in the Michaelis complex of Bacillus 1,3-1,4-␤-glucanase by means of hybrid quantum mechanics/molecular mechanics (QM/MM) simulations (35). In this approach, the atoms of the QM region evolve in time under the effect of the quantum mechanical forces, computed using density functional theory (DFT), and the electronic cloud adapts instantaneously to the chemical environment, whereas the forces on the MM region are ruled by a force field. The effect of the enzyme on the properties of the substrate is analyzed, providing insight on the factors leading to the substrate conformation we predict for this enzyme. The substrate chosen for the analysis is the 4-methylumbelliferyl tetrasaccharide shown in Fig. 3, which is a good substrate extensively used in enzyme kinetics (32)(33)(34). To the best of our knowledge, the conformational itinerary of the substrate in GHs has not yet been investigated by the first principles methods.

MATERIALS AND METHODS
Enzyme-Substrate Initial Structure-The only structures available for 1,3-1,4-␤-glucanase, a family GH16 enzyme, are that of the native enzyme (36,37) and that of the covalent enzyme-inhibitor complex with an epoxybutyl saccharide (36). In addition, the structure of the enzyme product complex has been recently characterized (PDB accession code 1U0A) (69). This structure was used to build the initial structure of the Michaelis complex. For this purpose, the missing methylumbelliferyl aglycon (MU in Fig. 3) was inserted manually to generate the enzymesubstrate complex. This structure was submitted to a preliminary molecular dynamics simulation (33) using the Cornell et al. force field  (38), as implemented in the HYPERCHEM package, 4 with the restraint that only the atoms 6 Å from the anomeric carbon were allowed to move. During this procedure the substrate maintained the original chair conformation ( 4 C 1 ) in the four sugar rings. To force a structure similar to that of an oxocarbenium ion, the charge on the anomeric carbon was increased by 0.6 electrons. After this replacement, the sugar ring adopted a distorted 1 S 3 conformation.
For the sake of simplification, hereafter we will use the notations 4 C 1substrate and 1 S 3 -substrate to refer to the substrate isomer in which the sugar ring of the Ϫ1 subsite adopts either 4 C 1 or 1 S 3 conformation, respectively. These two structures were used for the subsequent calculations.
First Principles Molecular Dynamics Simulations-First principles molecular dynamics simulations were performed to analyze the dynamics of the isolated substrate and that of the substrate in the presence of the catalytic residues (Glu 109 and Glu 105 ). The calculations were performed using the Car-Parrinello method (CP) (40,41), which is based on DFT. Both DFT and the CP method have already been employed with success in a number of investigations of biological processes (see for instance Refs. [42][43][44][45][46][47][48] including carbohydrate structure and dynamics (49 -52) and carbohydrate reactivity (see for instance Refs. [53][54][55][56][57][58]. The generalized gradient-corrected approximation of DFT, following the prescription of Perdew, Burke, and Ernzerhoff (59), was used. The choice of this functional is based on its reliability in the description of hydrogen bonds (60). We employed ab initio pseudopotentials, generated within the Troullier-Martins scheme (61). The Kohn-Sham orbitals (62) are expanded in a plane wave basis set with the kinetic energy cutoff of 70 Ry. Structural optimizations were performed by means of molecular dynamics with annealing of the atomic velocities, using a time step of 0.12 fs, and the fictitious mass of the electrons was set at 1200 a.u. With this setup the total energy and the fictitious kinetic energy of the electrons were conserved within 1.01 ⅐ 10 Ϫ6 a.u./ps⅐atom and 3.6 ⅐ 10 Ϫ5 a.u./ps⅐atom, respectively. The Nosé-Hoover thermostat for the nuclear degrees of freedom was used to maintain the temperature as constant as possible (63). The systems were enclosed in supercells of size 17.5 ϫ 12.5 ϫ 12.5 Å 3 (isolated substrate) and 15.8 ϫ 13.2 ϫ 17.0 Å 3 (substrate ϩ catalytic residues). The calculations were performed with the CPMD program, 5 and structure analysis was performed with VMD (65). Atomic charges were computed from the electrostatic potential (ESP). Interaction energies between the substrate and each of the catalytic residues were obtained by subtracting the energy of the complex from the energy of the isolated fragments in their corresponding optimized structure. The energy of the 4 H 3 transition conformation (the transition state for a conversion between the 1 S 3 and 4 C 1 conformers) was obtained by fixing the dihedral angle defined by the C2, C1, O5, and C5 ring atoms (hereafter referred as angle, shown in Fig.  3a) to zero degrees and optimizing all other degrees of freedom. Simi- larly, the 1 S 3 conformation (which is not a local minimum for the isolated substrate) was optimized by fixing the C2-O5-C5-C4 dihedral angle.
Hybrid QM/MM Molecular Dynamics Simulations-Hybrid QM/MM simulations on the complete protein were performed for each of the two substrate conformations described above ( 4 C 1 -substrate and 1 S 3 -substrate).
Before starting the QM/MM simulations, a classical molecular dynamics simulation was performed to equilibrate the protein and allow the substrate to accommodate in the binding cavity. The following parameters were used in the classical simulations. The Cornell et al. force field (38), as implemented in the AMBER 7.0 program (66), and the GLYCAM parameter set (67) were used for the protein residues and for the substrate, respectively. The MU aglycon was parameterized using the antechamber module. The atomic charges of the substrate were obtained from a first principles (Car-Parrinello) calculation of the isolated substrate. All His residues (located in the protein surface) were taken as protonated (i.e. positive charge), and all Asp and Glu residues were taken as deprotonated (i.e. negative charge) except Glu 109 (the acid-base residue) and Asp 107 , which is hydrogen-bonded to the nucleophile Glu 105 . Six chlorine atoms were added to achieve neutrality of the protein structure. The system was enveloped in a 52 Å ϫ 40 Å ϫ 66 Å box of equilibrated TIP3P water molecules and was equilibrated in several steps. First, all water molecules were relaxed with a gradient minimizer and then equilibrated for 20 ps at 150 K (protein constrained). Next, the whole system was minimized and subsequently equilibrated for 20 ps at 300 K. During equilibration, the system was coupled to a heat bath to achieve the desired temperature of 300 K. The simulation was continued for 20 ps at constant pressure, allowing the cell volume to evolve until equilibration. Analysis of the trajectories was carried out by using standard tools of AMBER. Several scenarios were tested for the initial protonation state of the acid/base (Glu 109 ) and the Asp 107 residues. The configuration that better maintained the interaction of the catalytic residues with the substrate upon a short test MD run was chosen for the production runs. In this configuration, the OH group of Glu 109 interacts with O gly and the OH group of Asp 107 interacts with Glu 105 .
Two separate classical simulations were performed, one with the substrate in the chair conformation ( 4 C 1 -substrate) and another one with the substrate in the distorted skew-boat conformation ( 1 S 3 -substrate). While the chair conformer was found to be stable, the skew-boat conformer evolved toward the chair one during the optimization process, unless a larger atomic charge for the anomeric carbon is used. For this reason, the simulation of the skew-boat conformer was performed by using a different charge in C1 (see above).
Once the system was equilibrated and the relative position of the substrate and the enzyme did not change (i.e. the root mean-square deviation variations were stabilized), QM/MM simulations were initiated. The method developed by Laio et al. (68) was used. This method combines the first principles molecular dynamics method of Car and Parrinello (CPMD) (40) with a force field molecular dynamics methodology (i.e QM/MM CPMD). In this approach, the system is partitioned into a QM fragment and an MM fragment. The dynamics of the atoms on the QM fragment depends on the electronic density, (r), computed with DFT, whereas that of the atoms on the MM fragment is ruled by an empirical force field. The QM-MM interface is modeled by the use of a monovalent pseudopotential that saturates the QM region (64). The electrostatic interactions between the QM and MM regions are handled via a fully Hamiltonian coupling scheme (68) where the short range electrostatic interactions between the QM and the MM regions were explicitly taken into account for all atoms. An appropriately modified Coulomb potential was used to ensure that no unphysical escape of the electronic density from the QM to the MM region occurs. The electrostatic interactions with the more distant MM atoms are treated via a multipole expansion. Bonded and van der Waals interactions between the QM and the MM region are treated with the standard AMBER force field. Long range electrostatic interactions between MM atoms have been described with P3M implementation (66). The mesh used for P3M was 64 ϫ 64 ϫ 64. A time step of 5 a.u. and a fictitious electron mass of 850 a.u. were used. The variation of the total energy and the fictitious kinetic energy of the electrons was less than 1.19 ⅐ 10 Ϫ6 a.u./ps⅐atom and 1.07 ⅐ 10 Ϫ6 a.u./ps⅐atom, respectively. Constant temperature was achieved by coupling the systems to a Nosé-Hoover thermostat (63) of 3500 cm Ϫ1 frequency. As an average, one QM/MM (i.e. 0.12 fs) step took 23.4 s using 16 processors of an IBM-SP4 computer.
Two types of QM-MM partitions were tested. In the first one, the QM fragment included the Ϫ1 sugar subsite and the MU aglycon (42 atoms). A second model including the carboxylic groups of Glu 109 , Glu 105 (the catalytic residues), and Asp 107 (the residue that is hydrogenbonded to Glu 105 ) in the QM region was also tested (62 atoms). However, no significant differences in terms of structure and dynamical behavior were found with respect of the use of a small QM region. For instance, the hydrogen bond distances between the catalytic residues and the saccharide (Glu 109 . . . O5 and Glu 105 . . . H-O 2 ) vary by only 1%. Therefore, only the latter results (small QM region) will be discussed here. Further details on the QM/MM methodology used here can be found in Ref. 68.

RESULTS
The Isolated Substrate-Before modeling the enzyme-substrate complex, it is useful to analyze the dynamics of the Ϫ1 sugar ring for the isolated substrate. Therefore, we built two models consisting of the Ϫ1 sugar ring in either the 1 S 3 (distorted) or 4 C 1 (undistorted) conformation, plus the MU aglycon. Only the 4 C 1 conformation was found to be a local minimum; whereas the 1 S 3 conformer evolved toward a 1,4 B boat conformer during the structure optimization.
Molecular dynamics simulations were performed starting from each of the two substrate conformations ( 1 S 3 and 4 C 1 ). As a way to sketch the flexibility of the Ϫ1 pyranose ring of the substrate, we monitored the dihedral angle (defined under "Materials and Methods") during the CPMD simulation. The value of this dihedral is positive (45°) for the skew-boat conformation and negative (Ϫ45°) for the chair conformation. Fig. 4 shows the ring conformations sampled during the CPMD simulation. As observed in Fig. 5a, the 1 S 3 -substrate is unstable and evolves toward the 4 C 1 (chair) conformation in 5-6 ps, sampling different conformations in its itinerary ( 1 S 3 , 1,4 B and 1 S 5 ). The 1 S 3 3 4 C 1 transition takes place via the 4 H 3 (half-chair) conformation ( ϭ 0). In contrast, a simulation starting from the 4 C 1 conformation (not represented) shows that the pyranose ring keeps the conformation during the dynamics.
It is worth noting that, for an isolated substrate, the 1 S 3 3 4 C 1 transition involves a large motion of the aglycon. As shown in Fig. 6a, the axial to equatorial change of the C1-O gly bond moves the MU aglycon from the ␤ region to the ␣ region with respect to the average sugar plane. Such transition would be impeded in the protein unless rotation around the C 1 -O gly bond takes place simultaneously (Fig. 6b). In this case, the aglycon would remain in the ␤ region during the transition. It will be seen later on that this is indeed the most likely path for a possible 1 S 3 3 4 C 1 transition in the protein.
Effect of the Catalytic Residues-To investigate the effect of the nearby protein environment on the conformational dynamics of the substrate, we extended our model including the closest residues to the Ϫ1 sugar ring. These are the Glu 109 and Glu 105 catalytic residues. Glu 109 interacts with the glycosidic oxygen and Glu 105 interacts with the 2-OH group of the saccharide. These two residues were modeled with propionic acid/propionate molecules, respectively. Their terminal methyl groups were kept fixed during the calculation to simulate the anchoring of these groups to the protein backbone.
The first difference we observed with respect to the calculations of the isolated substrate is that the 1 S 3 conformation is a local minimum (i.e. it is stable under optimization). However, it evolves toward the chair conformer when temperature effects are taken into account. As shown in Fig. 5b, the pyranose ring adopts several conformations until it undergoes a transition toward the most stable 4 C 1 chair conformer via the 4 H 3 transition conformation. The 1 S 3 3 4 C 1 transition takes place sooner (2.7 ps) than in the simulation of the isolated substrate (5.3 ps). Analysis of the boat to chair transformation in both cases provides an explanation for this fact. It appears that the 4 H 3 conformation is stabilized by a hydrogen bond between the pyranose oxygen and the acid/base residue (Fig. 2c). Thus, its relative stabilization is expected to lower the energy barrier for the skew to chair transition. In fact, the relative energy between the 1 S 3 and 4 H 3 conformers decreases by 6 kcal/mol (from 11.2 to 4.7 kcal/mol) when the catalytic residues are present (Fig. 7a), whereas the energy difference between the 1 S 3 and 4 C 1 conformations (Ϸ3 kcal/mol) remains unchanged. Therefore, even though the catalytic residues do not change the main dynamic features of the system (the sugar ring evolves toward the most stable chair conformation in a few ps) they affect the type of ring conformations accessible at a given temperature. In other words, they change quantitatively the shape of the potential energy surface with respect to ring distortions.
Another effect of the catalytic residues is the polarization of the electronic charge distribution of the substrate and, in particular, the charge on the anomeric carbon, q(C1). As shown in Table 1, the (positive) charge on C1 is larger for the skew-boat (0.28) than for the chair conformation (0.15). Moreover, it becomes larger in the presence of the two catalytic residues (0.32). This increase of charge in q(C 1 ) reminds the electronic structure of the transition state of the hydrolysis reaction, in which a oxocarbenium ion is formed (second structure in Fig. 1) (4). These results also explain why the substrate distorts toward a skew-boat conformation by increasing the charge on the anomeric carbon in the force field simulation ("Materials and Methods"). Such change just adapts the force field to the distorted conformation, at the cost of worsening the description of the chair conformation. Because the charge on C1 changes significantly with ring conformation, it would be difficult to describe 1 S 3 / 4 C 1 substrate conformational changes using standard force fields.
Analysis of the MD trajectory shows that the nucleophile residue (Glu 105 ) is permanently involved in a hydrogen bond interaction with the 2-OH group of the sugar ring (H . . . O ϭ 1.57 Ϯ 0.12 Å during the simulation starting from the 1 S 3 -substrate). This is expected in view of the large strength of the 2-OH . . . Ϫ OOC-Glu 105 hydrogen bond (33 kcal/mol). Instead, the acid/base residue (Glu 109 ), which interacts  weakly with the glycosidic oxygen (O gly . . . HOOC-Glu 109 ϭ 3 kcal/mol) appears to be quite mobile. The carboxyl group undergoes frequent oscillations around the C-C bonds and the OH group remains close to the glycosidic oxygen only at the beginning of the simulation for 3 ps. It is caused by this flexibility that the skew-boat to chair transition takes place, because it requires a large motion of the MU aglycon upon changing from an axial to equatorial position with respect to the sugar ring (Fig. 6a). The acid/base residue accommodates the aglycon in this motion. It is worth noting that such movement of the acid/base residue could not take place in the protein environment, because of steric interactions with the neighboring residues. As a way to mimic the steric restraint of the protein, we repeated the simulation fixing the Glu 109 residue. The transition from skew-boat to a chair conformation was suppressed, and the 1 S 3 substrate conformation was maintained during the whole simulation. Therefore, the mobility of the Glu 109 residue appears to affect significantly the energy barrier for the skew-boat to chair transition in this model. In the protein environment, other factors could play a role because a putative 1 S 3 3 4 C 1 transition requires a simultaneous rotation around the C 1 -O gly bond (Fig. 6b).
In summary, the calculations for a small model including the substrate and the catalytic residues show that these residues have an effect in increasing the charge of the anomeric carbon, and this effect is more pronounced for the 1 S 3 form than for the 4 C 1 form. The catalytic residues affect the energy barrier of the 1 S 3 3 4 C 1 transition and the relative energy among different conformations. However, this is not enough to stabilize a distorted conformation of the substrate and, therefore, we now consider the possible role of the protein environment.
Effect of the Protein Environment-Whereas gas phase calculations in small models are instructive, a more realistic approach must be based on calculations taking into account the complete protein environment. Thus, we performed additional calculations using the mixed first principles/classical MD methodology within the Car-Parrinello approach (QM/MM CPMD). The initial structure of the enzyme-substrate complex was taken from a classical MD simulation (see "Materials and Methods"). Two separate calculations were performed, starting with either the chair or the skew-boat conformer. The systems were optimized and then used as starting points for molecular dynamics simulation at 300 K for a total time of 15 ps ( 1 S 3 -substrate) and 13 ps ( 4 C 1 -substrate).
The simulation starting from the 1 S 3 -substrate shows that the substrate does not evolve toward the chair conformation (Fig. 5c, upper  line). Instead, it alternates between the skew-boat 1 S 3 and the boat 1,4 B conformations at a 1 S 3 / 1,4 B ratio of 2:1. Therefore, unlike what we found in small models including only a few protein residues, the distorted conformation is stable in the presence of the complete protein environment. The simulation starting from the 4 C 1 -substrate also results in a stable structure (Fig. 5c, lower line), which did not distort toward a skew-boat or boat conformation, but the 4 C 1 -substrate is 11 kcal/mol higher in energy than the 1 S 3 -substrate (Fig. 7b). Therefore, our results show that the 1 S 3 -substrate is the most favored form, but both chair and skew-boat conformers correspond to energy minima. Clearly, the effect of the protein is that of stabilizing more the skew-boat than the chair conformation, as the distorted conformation is not an energy minimum for the isolated substrate.
As mentioned before, a putative conformational change from 1 S 3 into 4 C 1 is expected to have a large energetic cost, as it requires a synchronized change of the and ⍀ dihedral angles (Fig. 6b). In this way the MU aglycon remains in the same region of the space during the skew-boat to chair transition. We indeed observed this transition in the classical MD simulation ("Materials and Methods"). As a way to estimate the energy barrier for a hypothetical 1 S 3 3 4 C 1 conversion, we optimized the structure of the Michaelis complex taking a snap-shot from the classical simulation in which the substrate is in the 4 H 3 transition state conformation and optimizing the structure fixing the angle. Because of the large steric hindrance associated with the 4 H 3 conformation, the energy barrier for the conversion from chair to skew-boat was found to be very large (80 kcal/mol). This value should be considered as an upper limit because the substrate conformation that we have optimized might not mach exactly the true transition state for the 1 S 3 3 4 C 1 conversion at the QM/MM level. Nevertheless, a sizable barrier is to be expected, since conversion between the two forms was not observed during the QM/MM simulations. Fig. 8 shows the optimized structures of the 4 C 1 -substrate and 1 S 3substrate complexes. Table 1 lists their most relevant structural parameters, in comparison with the previous models considered (the isolated substrate and with the catalytic residues). For the same type of substrate conformer, the internal structure is not significantly affected by the interaction with the protein. For instance, considering the chair conformation, the C1-O gly distance is similar in the isolated substrate than in the presence of the catalytic residues or the full protein. Only the H-bond distances show some variation, reflecting the strength of these interactions. For instance, the O 2 -H distance lengthens in the presence of the catalytic residues because of the strong 2-OH . . . Ϫ OOC-Glu 105 interaction. In contrast, the O6-H, which is involved in a weak hydrogen  Complex of 1,3-1,4 JANUARY 20, 2006 • VOLUME 281 • NUMBER 3 bond interaction with the acid/base residue, practically does not change. The most interesting changes occur upon substrate distortion from the 4 C 1 to the 1 S 3 conformation. It is apparent that the C1-O gly distance increases upon distortion (ϩ0.06 Å for the isolated substrate, ϩ0.07 Å when the catalytic residues are present, and ϩ0.06 Å in the protein). At the same time, the distance between the anomeric carbon and the sugar oxygen (C1-O5) decreases (Ϫ0.05Å for the isolated substrate, Ϫ0.04 Å when the catalytic residues are present and Ϫ0.03 Å in the protein). These changes are reminiscent of the ones occurring during the enzymatic reaction: the C1-O5 acquires partial double bond character and the C1-O gly distance lengthens in the oxocarbenium ionlike transition state (57). Therefore, our results once more show that the distortion brings the substrate to a structure closer to the transition state of the reaction. Additional calculations removing the MU aglycon show that the Ϫ1 sugar ring distortion is not stable any more and evolves toward the 4 C 1 conformation in Ϸ6 ps.

DISCUSSION
The substrate conformation in the Michaelis complex of GHs influences the conformational itinerary of the substrate during the catalytic reaction (5, 6); thus, it is a topic of ongoing interest in glycobiology. The substrate sugar ring of the Ϫ1 subset of the Michaelis complex of several GH enzymes (GHs) is known to adopt a distorted conformation that facilitates the hydrolysis reaction. However, the determinants of the substrate conformation, as well as its electronic/structural implications, are far from being understood.
In this work we have investigated the Michaelis complex of Bacillus 1,3-1,4-␤-glucanase, whose structure is not yet known, by means of hybrid first principles/classical molecular dynamics approaches. We have quantified for the first time the electronic/structural changes associated with a distorted substrate conformation and provided insight into the factors governing the conformation of the substrate.
A distorted (skew-boat) conformation is found to be unstable for the isolated substrate. This conformation is 3 kcal/mol higher in energy than the undistorted 4 C 1 conformer and does not correspond to a local minimum. In consequence, during the MD simulation the 1 S 3 conformation quickly evolves toward the chair conformation in a few picoseconds ( 1 S 3 3 4 C 1 transition in (Fig. 5a). The effect of the catalytic residues (Glu 109 and Glu 105 ) does not change this picture qualitatively, although the skew-boat conformation becomes a local minimum. Clearly, the effect of the catalytic residues is not enough to stabilize a distorted conformation of the substrate.
In contrast, the 1 S 3 skew-boat conformation is the most stable form for the substrate in the protein. The 4 C 1 -substrate corresponds to a higher energy configuration (11 kcal/mol), and both substrate isomers are separated by a large barrier (80 kcal/mol). Therefore, the protein changes qualitatively the relative energy among the 1 S 3 and 4 C 1 conformers with respect to an isolated substrate.
Analysis of the atomic charges on the substrate atoms show that the distortion increases the charge of the anomeric carbon and this effect is enhanced by the catalytic residues. In addition, the distortion of the substrate elongates the glycosidic C1-O gly bond and shortens the intraring C1-O5 bond. Again, these changes are in the direction of the transition state of the reaction, in which the glycosidic bond is partially broken, and the C1-O gly bond acquires partial double bond character (second structure in Fig. 1). Therefore, the substrate preorganization in a 1 S 3 conformation prepares it for the enzymatic reaction both from a structural point of view (the distorted structure is more similar to the transition state of the reaction than the undistorted one) as well as from an electronic point of view (positive charge is developed in the anomeric carbon). At this stage, the nucleophile is still far from the anomeric carbon (Glu 109 -O . . . C1 ϭ 3.50 Å, see Fig. 8b) but the substrate is already prepared for the enzymatic reaction.
The higher stability of the 1 S 3 -substrate can be rationalized in terms of enzyme-substrate interactions and the orientation of the aglycon. Any substrate conformation that fits into the binding cavity (Figs. 9 and 10) needs to have the aglycon in the ␤ region (Fig. 7b). The enzyme is thus engineered to select only those conformations presenting the aglycon in the ␤ region, as it is the case for the 1 S 3 skew-boat conformation. Other conformations that could in principle fulfill this condition are 1,4 B a All polar residues located within 5 Å distance from the anomeric carbon are considered. and 1 S 5 (Fig. 4). The former is indeed found during the simulation of the Michaelis complex (with a residence time that is half-times that of the 1 S 3 conformer). The latter is very unlikely because the strong interaction with the nucleophile keeps the 2-OH oriented toward it. In consequence, the 1 S 5 conformation is only observed for the isolated substrate.
It should be pointed out that the condition that the leaving group is in the ␤ region (Fig. 7b) is a steric effect independent of the nature of the leaving group. Therefore, we predict that a distorted conformation will be always found for the Michaelis complex of Bacillus 1,3-1,4-␤glucanase.
Independently of the initial conformer ( 1 S 3 -substrate or 4 C 1 -substrate), we found that the substrate penetrates into the binding cavity in the course of the initial classical MD simulation (see "Materials and Methods"). At the end of the simulation, the anomeric carbon is displaced by Ϸ1 Å with respect to the initial position. This indicates that in the enzyme-product complex (from which our E⅐S complex was built) the product is located slightly backwards (in the direction of the Ϫ2 subset) with respect to its position in the E⅐S complex. Indeed, it is consistent with the idea that the substrate displaces backwards once the strain, because of the terminal aglycon, is released after hydrolysis.
Analysis of the enzyme-substrate intermolecular interactions does not show significant differences between the 1 S 3 and 4 C 1 substrate conformers. The most salient features are the significant decrease of the C2OH . . . (Ϫ) OOC-Glu 105 distance (Ϫ0.26 Å) and the distance of the MU carbonyl group with two neighboring Asn residues (Fig. 9). Because of the different orientation of the 2-OH group, the C2OH . . . (Ϫ) -OOC-Glu 105 distance is 0.06 Å shorter in the 1 S 3 form than in the 4 C 1 form (Table 2) and this contributes to a better binding for the 1 S 3substrate. The stacking interaction of MU with Trp 192 is more efficient for the 1 S 3 -substrate than for the 4 C 1 -substrate. In addition, two hydrogen bond interactions of aglycon-C ϭ O . . . HN(Asn) type are found for the 1 S 3 -substrate and only one for the 4 C 1 -substrate. All these factors make the 1 S 3 conformer to be more tightly held in the protein than the 4 C 1 conformer. Another factor that raises the relative energy of the 4 C 1 -substrate is probably the fact that the aglycon is in the ␤ region ( Fig.  7b) instead of its preferred ␣ orientation (Fig. 7a) in the gas phase.
A completely opposite situation would be the absence of the leaving group. Calculations removing the MU aglycon show that the Ϫ1 sugar ring adopts the 4 C 1 conformation, and the whole substrate displaces backwards in the binding cavity. This again reinforces that the shape and interactions in the ϩ1 subsite are responsible for the substrate distortion. In fact, additional calculations in which the side chains of residues in the ␣ region of ϩ1 subsite are removed (Tyr 123 and Asn 129 replaced by Gly) are not enough to affect the substrate conformation. On the other hand, our results are consistent with the fact that in the x-ray structure of the enzyme-product complex (i.e. the complex of the enzyme with a 1,3-1,4-␤-glucan tetrasaccharide) all saccharide units adopt the 4

C 1 conformation.
An open question of the enzymatic mechanism of GHs is how the distorted conformation is achieved, i.e. whether the substrate distorts once it is accommodated in the active site, or it is already distorted before reaching the active site. This is a fundamental question that goes back to the general discussion on the induced fit effect of enzymes (39). Our study can provide some insight into this problem, for the particular case of Bacillus 1,3-1,4-␤-glucanase. The calculations show that the first step of the reaction is more favored for the 1 S 3 substrate than for the FIGURE 8. a, optimized structure for the 4 C 1 -substrate, obtained from the QM/MM calculation. b, optimized structure for the 1 S 3 -substrate, obtained from the QM/MM calculation. 1 C 4 -substrate. Assuming that only the 1 S 3 -substrate is reactive, the large barrier obtained for the 4 C 1 3 1 S 3 conversion (80 kcal/mol) indicates that a conformational change toward the reactive skew-boat conformer is not possible in the enzyme. Therefore, we suggest that the most effective reaction path corresponds to the situation in which the substrate is being distorted during the binding process.
In summary, our calculations predict that the substrate of Bacillus 1,3-1,4-␤-glucanase adopts a distorted 1 S 3 conformation upon binding to the enzyme. This distortion is caused by the shape and interactions in ϩ1 subsite and prepares the substrate for the enzymatic reaction both from a structural point of view as well as from an electronic point of view. Finally, we have shown that mixed QM/MM simulations are a useful tool to predict the conformation of the substrate in E⅐S complexes of GHs and to analyze its mechanistic implications.  . Stereoview of the accessible volume in the binding cavity for the most stable conformer of the ؊1 sugar ring ( 1 S 3 ). Areas corresponding to polar and neutral residues are represented in green and white, respectively, whereas the substrate is represented in ball and sticks. The solvent water molecules have been omitted for clarity.