The complete conformational free energy landscape of β-xylose reveals a two-fold catalytic itinerary for β-xylanases† †Electronic supplementary information (ESI) available. See DOI: 10.1039/c4sc02240h Click here for additional data file.

Ab initio conformational free energy landscapes, together with molecular dynamics simulations, enable to predict the catalytic itineraries of β-xylanase enzymes.


Introduction
Carbohydrates are the most abundant biomolecules on Earth. They have a huge diversity of roles, ranging from structural elements in the cell walls of bacteria and plants to cell-cell recognition processes of nonphotosynthetic cells. The vast amount and diversity of carbohydrate-based structures present in nature requires a large group of enzymes responsible for their metabolism. 1 This is the function of glycoside hydrolases (GHs), polysaccharide lyases (PLs) or glycosyltransferases (GTs), which constitute approximately 1-2% of the genome of any organism. 2,3 GHs are highly specic enzymes that catalyze the cleavage of glycosidic bonds in carbohydrates. These enzymes are systematically classied on the basis of amino acid sequence similarity into 133 families (see http://www.cazy.org and http:// www.cazypedia.org), 4 which generally share a common fold and reaction mechanism: retention or inversion of the anomeric conguration. With some exceptions, 5 retaining GHs follow a double displacement mechanism, with the formation of a covalent glycosyl-enzyme intermediate, whereas inverting GHs operate on one unique step. 6 Regardless of the mechanism, the transition state of the reaction has oxocarbenium ion-like character. 6 It is nowadays well established that the carbohydrate substrate undergoes critical conformational changes upon binding to the GH active site. 7,8 As observed in NMR and X-ray experiments, [9][10][11] the sugar ring located at the À1 enzyme subsite distorts away from its 4 C 1 conformation in solution towards a high-energy one (e.g., boat or skew-boat) upon binding to GHs. A fascinating thread of GH research, and one with major impact on the design of enzyme inhibitors, is the conformational analysis of reaction pathways within the diverse enzyme families. 7,8,12 Different enzyme families can harness different reaction pathways and, therefore, transition states. In general, GHs with a common type of mechanism (retention or inversion of the anomeric conguration) that act on similar substrates follow a common catalytic itinerary. 7,8 For instance, retaining b-glucosidases follow a 1 S 3 / [ 4 H 3 ] ǂ / 4 C 1 itinerary, whereas inverting b-glucosidases use 2 S O / [ 2,5 B] ǂ / 5 S 1 . These itineraries can be drawn on the Cremer-Pople 13 puckering sphere (a representation of all possible ring conformations dened by the puckering or polar coordinates Q, f and q; see Fig. 1). 14 However, simpler 2D projections, such as Stoddart or Mercator diagrams, are normally used to draw catalytic itineraries (Fig. 2). 8,15 b-xylanases, responsible for the hydrolysis of glycosidic bonds in b-xylans, are particularly interesting for their wide range of industrial applications, such as bioenergy and biofuel production from plant cell-wall biomass, in bread production, in pulp and paper industry and in the enhancement of animal feedstocks. [16][17][18] An intriguing aspect of b-xylanases is that, unlike other types of GHs, they have been proposed to follow several conformational itineraries (Fig. 2). In the case of retaining b-xylanases, three itineraries have been predicted: one that involves a 4 21 In addition, xylobio-imidazole and xylobiolactam oxime TS-like inhibitors have been found to adopt conformations close to 4 H 3 ( 4 E) in a GH10 xylanase 22 (this is particularly relevant in the light of a recent work on mannosidase inhibitors showing that imidazole derivatives exhibit good TS shape mimic properties). 23 Experiments on a GH39 xylanase also agree with a 1 S 3 / [ 4 H 3 ] ǂ / 4 C 1 itinerary. 24 Support for the second itinerary ( 2 S O / [ 2,5 B] ǂ / 5 S 1 ) mainly originates from reports on covalently bound substrates of family 11 GHs in which the xylose saccharide located at the À1 subsite (2-deoxy-2-uoroxylopyranosyl substrate) adopts a 2,5 B conformation, 25,26 as well as experiments on synthetic xyloside analogues locked on a 2,5 B conformation. 27 Computational modeling studies on B. circulans GH11 xylanase addressed this question, concluding that the conformation of the xylose saccharide in the Michaelis complex is either 2 S O or 2,5 B. 28,29 To add to the confusion, a recent structural study reports an almost undistorted conformation ( 4 C 1 / O E) for a xylohexaose substrate in complex with a GH11 xylanase mutant, predicting an unusual O E / [ O S 2 ] ǂ / B 2,5 itinerary (Fig. 2). 30 This is particularly surprising since the O S 2 conformation does not fulll the stereochemical requirements of an oxocarbenium ion-like transition state (ideally, C5, O5, C1 and C2 atoms should be coplanar). Finally, inverting b-xylanases are expected to follow 2 S O / [ 2,5 B] ǂ / 5 S 1 itinerary, 31,32 although only family 43 has been so far characterized. Therefore, the precise conformations followed by the substrate during catalysis in b-xylanases, especially those of family 11, have not been unambiguously resolved.
Detailed information concerning the energetic, structural and electronic relations of xylose is necessary for the understanding of xylosidase reaction mechanisms, especially to elucidate whether several catalytic itineraries for b-xylanases are possible. It has been previously shown that the conformational free energy landscape (FEL) of isolated simple sugars, obtained by ab initio metadynamics, informs about the conformations being preactivated for catalysis in GHs acting on the given sugar (e.g. bglucose in b-glucosidases, b-mannose for b-mannosidases, etc.) and, therefore, can be used to predict catalytic conformational itineraries. 8 Recent metadynamics studies using force-elds have contributed to popularize the representation of sugar conformational maps, 33,34 in the spirit of the seminal work by Dowd, French and Reilly. 35,36 Pederiva et al. [37][38][39] raised controversy about the inuence of the choice of the collective variables (i.e. Cartesian coordinates, q x and q y , or polar coordinates, q and f; see Fig. 1) in metadynamics simulations of sugar puckering. In addition, several works on six-membered ring conformations have also appeared in recent years (see e.g. ref. 40). Of particular interest is the study of Mayes et al., in which several sixmembered rings are studied using standard optimization techniques within Density Functional Theory (DFT). 41 Unfortunately, the relevant transition between 1 S 3 and 4 C 1 in b-xylose, which is one of the main pathways proposed experimentally for b-xylosidases, was not resolved. Therefore, a complete free energy landscape for b-xylose is still missing.
In this work, we rst demonstrate that the conformational free energy landscape obtained by ab initio metadynamics is invariant with respect of the type of collective variables (polar or Cartesian) used in the calculations. In a second step, we provide the FEL of b-xylose and use it to assess the catalytic itineraries that have been proposed for retaining b-xylosidases. Finally, we perform classical and quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations on two selected E$S complexes to evaluate the effect of enzyme mutation on substrate conformation.

Computational details
Stoddart and Mercator representations of the puckering sphere Any of the possible conformations of a six-atom sugar ring can be unequivocally assigned using the 3 Cremer and Pople puckering coordinates Q, f and q (Fig. 1). 13 The Q coordinate is the sum of the perpendicular distance of each ring atom (j) to Since these are polar coordinates, any ring distortion would fall within the puckering sphere-like volume (Fig. 1). While Q may differ among different conformations, q and f are sufficient to differentiate between all the conformers of Fig. 1. On the poles (q ¼ 0 or p) are located the two chair conformers ( 4 C 1 and 1 C 4 , respectively); on the equatorial region (q ¼ p/2) the 6 boat and 6 skew-boat structures are sequentially placed in steps of f ¼ p/6.
Two representations have been historically used to map the Cremer and Pople 3D plot into a simpler two dimensional plot. Stoddart diagram 42 corresponds to the projection of the polar coordinates onto the equatorial plane and can be dened by the Cartesian coordinates q x and q y (Fig. 2). 43 On another hand, the so-called plate carrée or Mercator representation is an equidistant cylindrical projection that results in a rectangular map with respect to q and f. 35 Which representation to use is essentially a matter of choice. Catalytic GH itineraries are easy to interpret when they are drawn on Stoddart representation, as there are no discontinuities between conformations. Moreover, since GHs normally follow straight itineraries either on the Northern or Southern hemispheres, but not on both (i.e. they do not follow transitions from 4 C 1 to 1 C 4 or vice versa), only one diagram is sufficient to sample the experimentally relevant conformations, which further simplies the problem. On another hand, the Mercator representation represents all conformations on a single diagram, thus mapping the complete Cremer-Pople sphere. This can be most convenient in a more general context, if one is interested in sampling all conformations of the sugar ring.

Ab initio molecular dynamics simulations
All the rst principles simulations at room temperature reported in this work were performed within the Car-Parrinello approach, 44 as implemented in the CPMD 3.15.1 program. 45 The analyzed systems consist of a single cyclohexane or b-xylose unit enclosed in an orthorhombic box of size 11Å Â 11Å Â 11.5Å and 12.5Å Â 13.5Å Â 11.6Å, respectively. The electronic structure was computed within the DFT, using the Perdew, Burke, and Ernzerhoff generalized gradient-corrected approximation (PBE). 46 This functional gave reliable results in previous Car-Parrinello simulations of isolated carbohydrates and GHs. In particular, the error on relative energies due to the DFT functional employed is AE0.6 kcal mol À1 for b-glucose. 43 Kohn-Sham orbitals were expanded in a plane wave (PW) basis set with a kinetic energy cutoff of 70 Ry. Norm-conserving pseudopotentials were employed, generated within the Troullier-Martins scheme. 47 The ctitious mass for the electronic degrees of freedom was set to 1000 au and 850 au for the cyclohexane and xylose systems, respectively, as used in previous studies of isolated sugars, 43,48 bglucose in solution 49 and carbohydrate-active enzymes. 32,50-53

Metadynamics simulations
Metadynamics (MTD) is a molecular dynamics based technique that enables the phase space exploration and the estimation of the FEL. This method is based on the reduction of the phase space dimensionality to a smaller set of so-called "collective variables" (CVs) which enclose the slowest modes that are relevant to the process of interest. The addition of small repulsive potentials on this reduced phase space during the molecular dynamics simulation enhances the exploration of the phase space. When convergence is reached, the inverse of the sum of the added potentials corresponds to the underlying free energy prole on the CVs space. 54 A detailed explanation of the method can be found elsewhere. 55 Metadynamics was used to sample all conformations of the puckering sphere of the 6membered rings analyzed here (cyclohexane and b-xylose) and obtain their conformational free energy landscapes.
The representation of the FEL with respect to the Mercator diagram ("Mercator FEL") was obtained using directly the f and q polar puckering coordinates as collective variables in the metadynamics simulation. On another hand, the "Stoddard FEL" was obtained from the Cartesian coordinates q x and q y (Fig. 2), dened as Q sin(q)sin(f) and Q sin(q)cos(f), respectively. Both types of collective variables (polar and Cartesian) were used for the simulations of cyclohexane, whereas only polar coordinates were used for b-xylose.
The height of the Gaussian terms was set to 0.13 kcal mol À1 (cyclohexane) and 0.18 kcal mol À1 (xylose), which ensures sufficient accuracy for the reconstruction of the FEL. The width of the Gaussian terms was set to 0.1 CV units in both cases, according to the oscillations of the selected collective variables observed in the free dynamics. The extended Car-Parrinello Lagrangian was used to describe the dynamics of the collective variables for the simulation using q x /q y variables, as done in our previous work. 43,48 The mass of the ctitious particle and the force constant of the coupling potential were tested to ensure that the coupled particle follows the value of the associated collective variable in the real system. Values of 5.0 amu for the mass of the ctitious particle and 0.5 au for the force constant fulll these conditions. A direct implementation of the MTD algorithm was used to obtain the Mercator FEL. A new Gaussian potential was added every 200-400 MD steps at the rst stage of the simulation, increasing it up to 1000 MD steps to ensure a proper convergence of the simulation. A total number of 3500/9000 Gaussian functions were added to completely explore the free energy landscape of cyclohexane (Mercator/Stoddart representations) and 4500 for b-xylose (Mercator representation). The convergence of the metadynamics simulations was assessed by checking the invariance of energy differences and the free energy landscape with the progression of the simulation, following ref. 56. The convergence error in the energies was found to be lower than 0.5 kcal mol À1 .

Modeling of E$S complexes
Classical MD and ab initio QM/MM MD simulations were performed to test the effect of enzyme mutations on the Michaelis complex of Trichoderma ressei GH11 xylanase in complex with xylohexaose (PDB entry 4HK8), 30 as well as the Michaelis complex of Streptomyces olivaceoviridis GH10 xylanase in complex with xylopentaose (PDB entry 2D24). 21 The AMBER 57 and CPMD 45 codes were used for the classical MD and ab initio QM/MM simulations, 58 respectively. The set up of the systems and simulation details are reported in the ESI (Sections 4 and 5 †).

Results and discussion
The conformational free energy landscape of cyclohexane Cyclohexane, the most simple six-membered ring for which experimental information of the relative energy among conformations is available, was used to test the dependence of the conformational free energy landscape with respect of the type of collective variables used in the metadynamics simulation (polar variables, q and f, leading to the Mercator representation, or Cartesian variables, q x and q y , leading to the Stoddart representation).
The conformational free energy landscape of cyclohexane, reconstructed from the metadynamics simulation using q and f as collective variables (i.e. Mercator representation), is shown in Fig. 3a. The most stable conformers are 4 C 1 and 1 C 4 , which are equivalent. 59 The skew-boat conformers correspond to local minima, 6.3 kcal mol À1 higher in energy than the chair conformers, whereas boat conformations lay 1 kcal mol À1 higher in energy and correspond to transition states between skew-boats. Interconversion from chair to any skew-boat crosses a half-chair conformer with an energy barrier of 10.5 kcal mol À1 . Half-chair and envelope conformations are energetically equivalent. These values are in good agreement with experimental estimates (DG ǂ chair/skew-boat ¼ 10.4-10.8 kcal mol À1 ; DG chair4skew-boat ¼ 4.7-6.2 kcal mol À1 ), 60-62 as well as previous theoretical predictions (see for instance ref. 40). It is noteworthy that the transition states (TS) between skew-boat and chair conformers occur at q angles of 60 (Northern hemisphere) and 120 (Southern hemisphere), thus they deviate from pure half-chair conformations. The last are located at 50.8 and 129.2 , respectively, and display four consecutive coplanar atoms. 13 A separate analysis by standard TS search methods using various levels of theory (page S2 of the ESI †) leads to similar results. The non-planarity of the transition state, already observed by Hendrickson in 1967, 63 has oen been overlooked in the literature. 64,65 The FEL obtained by using Cartesian coordinates (q x and q y ) as collective variables is shown in Fig. 3b. Similarly to the Mercator FEL, it reects the expected hexagonal symmetry, within statistical errors due to nite sampling of the phase space. 66 Remarkably, it exhibits the same features as the simulation using polar coordinates as collective variables, resulting in very similar free energy differences (DG ǂ chair/skew-boat ¼ 10.5 kcal mol À1 and DG chair4skew-boat ¼ 6.0 kcal mol À1 ). These results demonstrate that metadynamics simulations performed either with polar or Cartesian coordinates give equivalent results.
It is interesting to analyze not only the variation of the q and f coordinates, which determine ring conformation, but also the puckering amplitude Q (Fig. 1), which gives an idea of the degree of "puckering" of the ring. Analysis of the structures obtained in the metadynamics simulation evidences that the distribution of Q values, P(Q), features a bimodal behavior (Fig. 4), with peaks around 0.54Å and 0.72Å for the chair and skew-boat conformations, respectively. Therefore, the puckering volume is not spherical but ellipsoidal, with the Q radii increasing slowly from the poles to the equator.
The bimodal distribution of the puckering amplitude contrast with the results obtained in a previous study using force-elds, 37 which show a unique narrow peak around 0.5Å for the probability distribution. This type of unimodal distribution would lead to anomalous values of the <CCC angles, creating substantial ring strain on the boat and skew-boat conformations. In the above mentioned work, 37 it was also argued that calculations using the Cartesian coordinates q x and q y can affect the range of sampled conformations at different Q. This conclusion was based on the plot of the puckering amplitude Q vs. q r , where q r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi q x 2 þ q y 2 p is the projected puckered vector (see Fig. 1), which shows a strong correlation for q r values larger than 0.5Å. However, as it is shown in Fig. S2, † the two metadynamics simulations performed in this work display a very similar dual correlation between q r and Q, which is consistent with the puckering amplitude displaying a bimodal distribution. Therefore, both Cartesian and polar CVs lead to a similar correlation between Q and q r . 67 The conformational free energy landscape of b-xylose The free energy landscape of b-D-xylose in the Mercator representation is shown in Fig. 5. It contains several local minima. Two of them at the poles, corresponding to the two chair conformers ( 4 C 1 and 1 C 4 , at q ¼ 0 and 180 , respectively), as well as a two local minima in the equator (q ¼ 90 ): a wide minimum centered at 2 S O /B 3,O (f ¼ 150-180 ) and a small minimum centered at 3 S 1 / 3,O B (f ¼ 0-30 ). As previously observed for other sugar molecules, 43,48 not all minima have a direct correspondence with the ideal conformations of the Mercator diagram. For instance, conformations B 1,4 and B 2,5 do not correspond to energy minima but lay in high energy regions of the FEL. On another hand, there are minima in between two canonical conformations (e.g. 2

S O /B 3,O ).
It is clear from Fig. 5 that the 4 C 1 chair is the global minimum. The inverted chair, 1 C 4 , is 4 kcal mol À1 higher in energy, closely followed by the mixed 2 S O /B 3,O conformation, which is 5 kcal mol À1 above 4 C 1 . The remaining minima ( 3,O B/ 3 S 1 , 1 S 5 , 5 S 1 ) are much higher in energy ($7 kcal mol À1 ). A noticeable feature of the FEL is the presence of a wide valley on the equator, covering conformations 1 S 3 -B 3,O -2 S O -2,5 B-5 S 1 , which is connected with the 4 C 1 global minimum by a low energy region at q ¼ 75 .
One might wonder whether the particular shape of the FEL is related with the occurrence of distorted conformations in bxylosidase complexes. To this aim, we obtained the puckering coordinates of the sugar ring at the À1 subsite in all available structures of b-xylosidase complexes (Table 1), either Michaelis complexes (E$S) or glycosyl-enzyme covalent intermediates. Most of these complexes were obtained from enzyme mutants. It turned out that all experimental structures of Michaelis complexes (Fig. 5) and their corresponding TS predicted structures are located in low-energy regions of the free energy   (Table 1) correspond to complexes of GH10 and GH39 b-xylosidases, for which a 4 H 3 transition state is expected. Even though none of the experimental points and the predicted TS lie precisely at the center of a FEL minimum, they fall into a low energy region of the FEL. Structure 7 (see Table 1) is a Michaelis complex of a GH11 xylanase. Structure 8, corresponding to a glutamine mutant of the acid/base residue, is anomalous, as it features a À1 xylose saccharide in a practically undistorted conformation, which is intermediate between 4 C 1 and O E (the origin of this peculiar conformation will be discussed later on). Structures 9-15 (see Table 1 and blue stars near 2,5 B), correspond to covalent intermediates of the hydrolysis reaction of GH11 b-xylosidases, for which a 2,5 B transition state has been proposed. Structure 16 is a Michaelis complex of a GH120 b-xylanase from Thermoanaerobacterium saccharolyticum. Overall, Michaelis complex conformations and predicted transition states fall in low energy regions of the FEL, whereas there are no experimental structures in high energy regions. Therefore, we can conclude that the calculations of the isolated sugar molecules reveal lower free energy regions that agree with the possible conformational routes followed by different enzyme families. This does not mean, of course, that a given particular enzyme environment will not inuence the conformational landscape reported here. In fact, it has been shown for an a-mannosidase 15 and a bglucosidase 51 that the enzyme further restricts the accessible conformations of the sugar residue with respect to the ones observed in the gas phase. However, most of the conformations accessible on-enzyme are already reected on the FEL of the isolated molecule, which can be considered as a ngerprint of the given sugar. As mentioned before, catalytic conformational itineraries have been oen drawn with respect of Stoddart diagram (Fig. 2, le panel), 7,68 in which all conformations are interconnected. Therefore, we transformed the Mercator map of Fig. 5 into a Stoddart representation (Fig. 6a) 69 to identify conformational itineraries. It is easy to see that two of the itineraries that have been experimentally proposed ( 1 S 3 / [ 4 H 3 ] ǂ / 4 C 1 and 2 S O / [ 2,5 B] ǂ / 5 S 1 ) cover the low energy regions of the FEL. In  contrast, the O E / [ O S 2 ] ǂ / B 2,5 itinerary, recently proposed for Trichoderma ressei GH11 xylanase, passes through a high energy region of the FEL. In the light of this and our previous studies on isolated sugars (b-D-glucose, 43 aand b-D-mannose, 15,48 a-Lfucose 70 ) this is unprecedented and it might indicate that trends derived from the computed FEL break down in this case. However, the fact that the O S 2 conformation proposed as TS does not fulll the stereochemical requirements of an oxocarbenium ion-like transition state (C5, O5, C1 and C2 atoms should be coplanar) suggests that the complex obtained with the modied enzyme might correspond to a non-productive complex (or a pre-Michaelis complex). To test the above hypothesis, we performed classical and ab initio QM/MM molecular dynamics simulations on the structure reported by Wan et al., 30 which corresponds to the complex of the E177Q enzyme mutant with xylohexaose (E177 is the catalytic acid/base residue). To assess the effect of the mutation on the substrate conformation, the simulations were performed on both the mutant and the WT enzyme (in the last case, we manually reverted the Glu / Gln mutation of the acid/base residue). The simulations for the enzyme mutant complex fairly reproduced the experimental conformation of the xylan hexasaccharide: the À1 xylose ring adopts an intermediate O E/ 4 C 1 conformation (see Fig. 7a). However, once the mutation was reverted, the conformation of the À1 xylose ring started oscillating between 4 C 1 / O E and 2 S O (Fig. 7b), which is precisely the conformation that has been proposed for GH11 and it is supported by our computed FEL for b-xylose. In addition, QM/MM calculations on the WT enzyme complex show that the 2 S O conformation is energetically more stable than 4 C 1 / O E, whereas the opposite is found for the mutant complex (see ESI † for details). This is probably due to the structure of the modied enzyme not being able to reproduce the hydrogen bond pattern of the WT enzyme. While a neighboring Asn44 interacts with the acid/base residue acting as a hydrogen bond donor in the WT enzyme, the opposite is found in the mutant enzyme (i.e. Asn44 acts as a hydrogen bond acceptor, see Fig. S3 †). As a control calculation, we considered an opposite case in which the complex obtained with a modied enzyme conforms with our computed FEL. We selected a family 10 modied retaining xylanase from Streptomyces olivaceoviridis (PDB code 2D24, Table 1), which was crystallized with a natural substrate, as the previously analyzed complex. The simulations show that there are two stable substrate conformations in the active site ( 1 S 3 /B 3,O and 4 C 1 / 4 E), with the distorted one ( 1 S 3 /B 3,O ) being the most stable (Fig. S4 of the ESI †). Most importantly, both the modied and the wild type enzyme behave similarly with respect to the conformation of the xyloside at the À1 subsite. Therefore, the complex of modied Streptomyces olivaceoviridis GH10 with its natural substrate is a good mimic of the WT Michaelis complex (the "true" Michaelis complex), being in the pathway towards the transition state. In contrast, the complex of Trichoderma ressei GH11 xylanase with its natural substrate corresponds to a non-productive complex, in which the substrate has been enforced to adopt a different conformation than the one of the true Michaelis complex. The 4 C 1 / O E conformation observed in this complex cannot be taken as informative of conformational itinerary.
To analyze how different sugars stabilize different distorted conformations, thus favoring different catalytic itineraries, we compared the computed FEL of b-xylose with the ones previously reported for b-glucose and b-mannose ( Fig. 6b and c,  respectively). In all cases, only the Northern hemisphere of the Cremer-Pople sphere (Fig. 1) is mechanistically relevant, as found in all experimental structures of b-xylosidases, b-glucosidases and b-mannosidases. Therefore, only the Northern hemisphere is shown in Fig. 6. The landscapes of b-xylose and bglucose share common features. In both cases, the lowest energy region is in the bottom part of the diagram. However, whereas b-glucose exhibits a second low energy region in the Northwestern quadrant (around B 2,5 ), b-xylose contains high energy conformers in this region. These differences can be traced back to changes in sugar ring substitutions. Xylose differs from glucose (and mannose) in the absence of the hydroxymethyl (CH 2 OH) substituent. In glucose and mannose, the CH 2 OH forms intramolecular hydrogen bonds in conformations in the vicinity of B 2,5 (Fig. 8), something that is not possible in b-xylose. The stability of the Northwestern quadrant is even more pronounced for b-mannose (Fig. 6c) compared to b-glucose and b-xylose. This is due to the particular orientation of the 2-OH substituent in b-mannose (trans with respect to H-5, preventing steric interaction between the two) compared to b-glucose and b-xylose (cis with respect to H-5). Similarly, the Southeastern quadrant of mannose is higher in energy compared to b-xylose and b-glucose due to unfavorable interactions between the C2 exocyclic group (2-OH) and the hydroxymethyl (CH 2 OH) substituent.
Previous studies of GH conformational maps have revealed that properties other than the free energy alone, 48 such as the partial charges of C1, O1, and O5, the C1-O1 and C1-O5 bond distances and the orientation of the C1-O1 bond (axial or equatorial, U), together with the relative free energy, should be taken into account to identify conformations that are preactivated for catalysis. Following a similar procedure as in ref. 48, we computed all these properties for a number of optimized structures (350) taken from the metadynamics simulation and combined them into a unique index (x) (see ESI, † pages S6-S12). This index should reect the likelihood that a given sugar conformation is present in the Michaelis complex of a b-xylosidase, i.e. the best preactivated conformation for catalysis. The calculations show that conformations 2 S O and 1 S 3 have high x values (Fig. 9), whereas O S 2 , 4 C 1 , 3

Summary and conclusions
In this work, we have computed the conformational free energy landscape for cyclohexane and b-xylose. The structure and relative free energy of the eight energy minima of cyclohexane (two chair structures, 4 C 1 and 1 C 4 , which are the global minima, and six skew-boat structures) and their interconversion energy barriers agree well with the experimental information available. Interestingly, we found out that the structure of the TS between the chair and the skew-boat conformers deviates by 9.2 in the q puckering coordinate with respect to a "pure" half-chair conformer (q of 50.8 and 129.2 ), which exhibits four consecutive coplanar atoms. Our simulations also show that, contrary to the conclusions of a force-eld-based metadynamics investigation, 37 the free energy landscape remains unaltered by the use of polar or Cartesian collective variables. As already reported in previous ab initio metadynamics simulations of pyranoses, 15 but not reproduced in the above force-eld-based study, 37 the puckering amplitude sphere is not spherical but elliptical, with a smaller radius at the poles compared to the equator.
The conformational FEL of b-xylose contains two local minima, corresponding to the two chair conformers ( 4 C 1 and   1 C 4 ), as well as wide minimum in between B 3,O and O S 2 . As previously found for other pyranoses, the puckering coordinates of all the À1 xylose saccharides extracted from experimental structures of b-xylosidase Michaelis complexes are located, without exceptions, in low-energy regions of the FEL. This nding, which is likely to be general in GH complexes, suggests that the factors governing the distortions present in these complexes are dictated by the intrinsic properties of a single saccharide unit (b-glucose in b-glucosidases, 43 b-mannose in bmannosidases, 48 a-mannose for a-mannosidases, 15 a-fucose for a-fucosidases 70 and b-xylose for b-xylosidases) and that enzymesubstrate interactions have evolved to fulll all those criteria for efficient catalysis.
Two itineraries, 1 S 3 / [ 4 H 3 ] ǂ / 4 C 1 and 2 S O / [ 2,5 B] ǂ / 5 S 1 , which pass through low energy regions of the FEL (Fig. 6a), are consistent with our simulations. Conformations 2 S O and 1 S 3 are preactivated for catalysis in terms of free energy, anomeric charge and bond distances (Fig. 9) route passes through a high energy region of the FEL. Additional classical and QM/MM molecular dynamics simulations on the complex of a xylose hexasaccharide with Trichoderma ressei GH11 xylanase show that the 4 C 1 / O E conformation observed in the reported enzyme complex is due to enzyme mutation. Although this cannot be taken as general (similar calculations on the complex of modied Streptomyces olivaceoviridis GH10 with its natural substrate show that it is a good mimic of the true Michaelis complex), it indicates that complexes obtained from modied enzymes do not necessarily represent the structures that take place in the reaction coordinate.
In summary, the simulations of the FEL of b-xylose reveal lower free energy regions that agree with two of the three conformational routes previously proposed for different bxylosidase families. We also show that the preactivation index is a good measure of the general preferences of b-xylosidases to accommodate specic b-xylose conformations at the À1 enzyme subsite. Most likely, the higher exibility of the xylose unit in comparison with other pyranoses, due to the lack of the hydrogen-bond interactions with the CH 2 OH group, can cause the xylose ring to adopt different conformations ( 2 S O and 1 S 3 ) with small changes of the active site structure, as encountered in different xylanase families.