Enzyme Molecular Mechanism as a Starting Point to Design New

DOI: 10.1021/jp202079e. Publication Date (Web): May 4, 2011. Copyright © 2011 American Chemical Society. *(V.M.) E-mail: [email protected]. Phone: (+34)...
0 downloads 4 Views 4MB Size
ARTICLE pubs.acs.org/JPCB

Enzyme Molecular Mechanism as a Starting Point to Design New Inhibitors: A Theoretical Study of O-GlcNAcase Jer^onimo Lameira,† Claudio Nahum Alves,*,† I~naki Tu~non,‡ Sergio Martí,§ and Vicent Moliner*,§ †

Laboratorio de Planejamento e Desenvolvimento de Farmacos, Instituto de Ci^encias Exatas e Naturais, Universidade Federal do Para, CP 11101, 66075-110 Belem, PA, Brazil ‡ Departament de Química Física, Universitat de Valencia, 46100 Burjassot, Valencia, Spain § Departament de Química Física i Analítica, Universitat Jaume I, 12071 Castellon, Spain

bS Supporting Information ABSTRACT: O-Glycoprotein 2-acetamino-2-deoxy-β-D-glucopyranosidase (O-GlcNAcase) hydrolyzes O-linked 2-acetamido-2deoxy-β-D-glucopyranoside (O-GlcNAc) residues from post-translationally modified serine/threonine residues of nucleocytoplasmic protein. The chemical process involves substrate-assisted catalysis, where two aspartate residues have been identified as the two key catalytic residues of O-GlcNAcase. In this report, the first step of the catalytic mechanism used by O-GlcNAcase involving substrateassisted catalysis has been studied using a hybrid quantum mechanical/molecular mechanical (QM/MM) Molecular Dynamics (MD) calculations. The free energy profile shows that the formation of the oxazoline intermediate in the O-GlcNAcase catalytic reaction takes place by means of a stepwise mechanism. The first step would be a cyclization of the acetomide group, which seems to be dependent on the proton transfer from a conserved aspartate, Asp298 in Clostridium perfringens O-GlcNAcase. From this new intermediate, a proton is transferred from the azoline ring to another conserved aspartate (Asp297) thus forming the oxazoline ion and departure of the aglycone. In addition, averaged values of proteinsubstrate interaction energy along the reaction path shows that, in fact, the transition states present the highest binding affinities. A deeper analysis of the binding contribution of the individual residues shows that Asp297, Asp298, and Asp401 are basically responsible of the stabilization of these complexes. These results would explain why O-(2-acetamido2deoxy-D-glucopyranosylidene)amino-N-phenycarbamate (PUGNAc), 1,2-dideoxy-20 -methyl-R-D-glucopyranoso-[2,1-d]-Δ20 -thiazoline (NAG-thiazoline), and GlcNAcstatin derivatives are potent inhibitors of this enzyme, resembling the two transition states of the O-GlcNAcase catalytic reaction path. These results may be useful to rational design compounds with more interesting inhibitory activity.

’ INTRODUCTION Typically about 1% of the genome of any organism encodes glycoside hydrolases (GHs), and collectively, many thousands of genes encoding GHs have been annotated, with structural similarities into over 100 families of these proteins.1 The field of glycobiology was turned inside-out more than 20 years ago by the discovery of O-linked 2-acetamido-2-deoxy-β-Dglucopyranoside (O-GlcNAc) from serine and threonine residues of post-translationally modified nuclear and cytoplasmic proteins.2,3 On several proteins, O-GlcNAc and O-phosphate alternatively occupy the same or adjacent sites, leading to the hypothesis that one function of this saccharide is to transiently block phosphorylation.4 Recently, perturbations in the regulation of O-GlcNAc have been related with type II diabetes5,6 and neurodegenerative disorders such as Parkinson dystonia7 and Alzheimer’s disease,8 which is the reason that O-GlcNAcase was used as target for therapeutic agents.911 O-Glycoprotein 2-acetamino-2-deoxy-β-D-glucopyranosidase (O-GlcNAcase) is a member of the family 84 glycoside hydrolases r 2011 American Chemical Society

(GH 84),12 as defined by the CAZY database, that hydrolyzes O-GlcNAc residues from post-translationally modified serine/ threonine residues of nucleocytoplasmic protein. O-GlcNAcase catalyzed reaction takes place through a two-step mechanism, first forming a transient oxazoline intermediate (P structure in Scheme 1) that is subsequently broken. Oxazoline formation, which is the step studied in this paper, is a substrate-assisted catalytic process, where the 2-acetomide group is activated by an aspartate residue while departure of the aglycone leaving group is assisted by means of a proton transfer from a different protonated aspartate residue (see the hypothetical transition structure, TS, in Scheme 1). These two key residues have been identified as Asp174 and Asp175 in human O-GlcNAcase13,14 and Asp297 and Asp298 in Clostridium perfringens.15 It has been proposed that in this first step of the O-GlcNAcase catalyzed Received: March 4, 2011 Revised: April 14, 2011 Published: May 04, 2011 6764

dx.doi.org/10.1021/jp202079e | J. Phys. Chem. B 2011, 115, 6764–6775

The Journal of Physical Chemistry B

ARTICLE

Scheme 1. Proposed Mechanism of the First Step of the Reaction Catalyzed by GlcNAcasea

a

Key distances employed to explore the molecular mechanism are depicted in Rboat.

Scheme 2. Structures of Known Inhibitors of GlcNAcase

reaction Asp174 directs and polarizes the 2-acetamido group to act as a nucleophile and forming the oxazoline intermediate.13 Asp175 meanwhile acts as a general acid, encouraging departure of the aglycone leaving group.13 In the second step, ringopening, Asp174 would facilitate departure of the 2-acetamido group, whereas Asp175 would act as a general base, promoting the attack of a water molecule to yield the β-hemiacetal product.13 According to this mechanism, oxocarbenium ionlike species with sp2-hybridized anomeric center and delocalized positive charge between the endocyclic ring oxygen and the anomeric center have been proposed as the transition state (TS) for enzymes using substrate-assisted mechanisms.1618 The hypothetical TS for the O-GlcNAcase (TS in Scheme 1) has been used as a model for rational design of inhibitors of this enzyme.19,20 Among these inhibitors, O-(2-acetamido-2deoxyD-glucopyranosylidene)amino-N-phenycarbamate (PUGNAc), 1,2-dideoxy-20 -methyl-R-D-glucopyranoso-[2,1-d]-Δ20 -thiazoline (NAG-thiazoline), and GlcNAcstatin B (see Scheme 2 for detail

of the structures) are the best-characterized inhibitors of human O-GlcNAcase.2022 Recently, the crystal structure of a bacterial Clostridium perfringens homologue (CpGH84H) in complex with PUGNAc has been reported, with significant sequence homology to the human O-GlcNAcase.14 Moreover, similarities between the inhibition constants of human O-GlcNAcase and CpGH84H were observed, making this bacterial enzyme a good model of the human enzyme.14 An analysis of the enzyme-PUGNAc complex carried out by Rao et al.15 suggests that this molecule must inhibit O-GlcNAcase by mimicking the oxocarbenium ion-like TS of the O-GlcNAcase-catalyzed hydrolysis of N-acetylglucosaminide by virtue of its sp2 anomeric C1, whereas NAG-thiazoline would be geometrically similar to the oxazoline intermediate (P in Scheme 1). Nevertheless, Vocadlo and coworkers have shown that NAG-thiazoline is a TS analogue (TSA) for O-GlcNAcase, whereas PUGNAc is a poor TSA or a serendipitous binding inhibitor.21 Finally, GlcNAcstatin B derivative, the most potent human O-GlcNAcase inhibitor reported to date, is a 6765

dx.doi.org/10.1021/jp202079e |J. Phys. Chem. B 2011, 115, 6764–6775

The Journal of Physical Chemistry B novel scaffold obtained by exploiting the structural similarity of Z-PUGNAc and the naturally occurring potent human hexosaminidases A/B (HexA/B) inhibitor nagstatin.23 Recently, we have carried out molecular dynamics (MD) simulations using a hybrid quantum mechanics/molecular mechanics (QM/MM) approach to determine the binding of PUGNAc and NAG-thiazoline, with a bacterial O-GlcNAcase.24 Our calculations allowed identifying the key residues in the proteininhibitor interaction, concluding that both complexes presented very close values, and in any case, the closer the structure was from the hypothetical TS, the better the inhibitor would be.24 Transition state theory (TST) has been used to explain why enzymes can catalyze chemical reactions, and one of its main goals is the successful application of the concept of TSA.25 The ideas of molecular structure and energetics are initially unified by potential energy surfaces (PES). An understanding of a protein PES is expected to give a good insight into the protein chemical reaction.26 In order to carry out these kinds of studies for such large molecular systems, hybrid QM/MM schemes can be used, where a small part of the system is described by quantum mechanics and the rest of the system by a classical force field.25 Free energy profiles can thus be obtained which render values directly compared with experimental data. The free energy profile of the transformation from reactants to products can be obtained as a potential of mean force (PMF) from biased MD simulations using an adequate distinguished reaction coordinate. The maximum of this profile corresponds to the TS of the reaction while the free energy difference between this maximum and the minimum, which corresponds to reactants, will be related with the chemical rate constant. This free energy barrier can be related and compared with experimental rate constants since averaged structures of reactants, TS, and possible intermediates of the chemical reaction are obtained under the protein environment effect by means of a flexible and realistic model.28 In this study, we have used both potential energy profiles and free energy profiles for studying the catalytic mechanism used by O-GlcNAcase involving substrate-assisted catalysis. In addition, a detailed analysis of the interactions of the substrate with the key residues inside the binding pocket in reactants, TS and intermediate states, has been carried out in order to decide which structure can present the most efficient interactions with the protein, thus becoming the most effective model to design efficient inhibitors.

’ METHODS The System. We have recently carried out MD simulations, using a combined QM/MM approach, for studying protein inhibitors interaction energies.24,29 In these calculations, the small part of the system described quantum mechanically was the ligand/substrate species, whereas the protein and solvent environment were represented by MM force fields. This hybrid methodology, by comparison with methods based on just MM, avoids most of the work needed to obtain new force field parameters for each new species. Treating the ligand quantum mechanically and the protein molecular mechanically has the additional advantage of the inclusion of quantum effects such as ligand polarization upon binding.29 Moreover, as the largest part of the system is described classically, enough sampling can be obtained at reasonable computational cost.

ARTICLE

In the work reported here, the initial coordinates for the QM/ MM MD calculations were taken from the trimeric crystal structure of CpGH84H-PUGNAc complex at 2.25 Å resolution, PDB code 2CBJ.14 This crystal structure (from Clostridium perfringens) presents significant sequence homology to the hOGA N-terminus (CpNagJ, 34% sequence identity and 51% sequence similarity with hOGA) and thus it is a suitable model for studies into the mechanisms of recognition/specificity of O-GlcNAcylated peptides.14 Since the standard pKa values of ionizable groups can be shifted by local protein environments,30 an accurate assignment of the protonation states of all these residues at pH = 7 was carried out. Recalculation of the pKa values of the titratable amino acids has been done using the “cluster method” as implemented by Field and co-workers.31 In this method, each titratable residue is perturbed by the electrostatic effect of the protein environment. We have also calculated the pKa values of aminoacids within the empirical propKa program of Jensen et al.32 Results obtained with both methods are qualitatively in agreement, predicting the same protonation state for all residues at the selected pH of the simulations. Thus, according to these results, most residues were found at their standard protonation state, except for the Asp-298 residue that was protonated. In fact, as demonstrated from a previous study where a QM/MM PMF corresponding to the proton transfer from the carboxylate oxygen atom of Asp298 to the oxime nitrogen atom of PUGNAc was traced, the PUGNAc-protein complex was found to be more stable with the proton on the Asp298.33 After adding the hydrogen atoms to the structure, series of optimization algorithms (steepest descent conjugated gradient and L-BFGS-B34) were applied. To avoid a denaturation of the protein structure, all the heavy atoms of the protein and the inhibitor were restrained by means of a Cartesian harmonic umbrella with a force constant of 1000 kJ mol1 Å2. Afterward, the system was fully relaxed, but the peptidic backbone was restrained with a lower constant of 100 kJ mol1 Å2. The optimized protein was placed in a cubic box of preequilibrated waters (80 Å side), using the principal axis of the protein-inhibitor complex as the geometrical center. Any water with an oxygen atom lying in a radius of 2.8 Å from a heavy atom of the protein was deleted. The remaining water molecules were then relaxed using optimization algorithms. Finally, 100 ps of hybrid QM/MM Langevin-Verlet MD (NVT) at 300K were used to equilibrate the solvent. During the MD simulations, the atoms of the substrate and the side chains of Asp297 and Asp298 (present in active site) were selected to be treated by QM, using a semiempirial AM1 Hamiltonian.35 The rest of the system (protein plus water molecules) were described using the OPLS-AA36 and TIP3P37 force fields, respectively, as implemented in the fDYNAMO library.38 To saturate the valence of the QM/MM frontier atoms, link atoms between the CR and Cβ atoms of Asp297 and Asp298 residues were used. The final system contains a total of 23 475 atoms, 50 of them in the QM region. Due to the large amount of degrees of freedom, any residue 20 Å apart from any of the atoms of the initial inhibitor was kept frozen in the remaining calculations. Cut-offs for the nonbonding interactions were applied using a switching scheme, within a range radius from 14.5 to 16 Å. Afterward, the system was equilibrated by means of 1.1 ns of QM/MM MD at temperature of 300 K. The computed rmsd for the protein during the last 500 ps renders a value always below 6766

dx.doi.org/10.1021/jp202079e |J. Phys. Chem. B 2011, 115, 6764–6775

The Journal of Physical Chemistry B

ARTICLE

0.9 Å. Furthermore, the rms of the temperature along the different equilibration steps was always lower than 2.5 K and the variation coefficient of the potential energy during the dynamics simulations was never higher than 0.3%. The Energy Function. In the work reported here, we have made use of a hybrid AM1/OPLS-AA method, based on a semiclassical approach. The potential energy of our scheme is derived from the standard QM/MM formulation ^ o jΨæ EQM=MM ¼ ÆΨjH 0 *   + q    MM þ@ Ψ Ψ þ re, MM 



þ EvdW QM=MM

1

∑∑

The activation free energy can be then expressed as46

ZQM qMM A rQM, MM

þ EMM

ΔG‡ ðξÞ ¼ Wðξ‡ Þ  ½WðξR Þ þ Gξ ðξR Þ ð1Þ

vdW EQM=MM ¼ Evac þ Eelect QM=MM þ EQM=MM þ EMM

ð2Þ

where EMM is the energy of the MM subsystem terms, EvdW QM/MM the van der Waals interaction energy between the QM and MM subsystems and Eelect QM/MM includes both the Coulombic interaction of the QM nuclei (ZQM) and the electrostatic interaction of the polarized electronic wave function (Ψ) with the charges of the protein (qMM). Finally, the interaction energy between the inhibitor and the environment is evaluated as the difference between the QM/MM energy and the energies of the separated, noninteracting, QM and MM subsystems with the same geometry. Considering that the MM part is described using a nonpolarizable potential, the interaction energy contribution of each residue of the protein is given by the following expression:  + *  q    MM Ψ EInt Ψ QM=MM ¼ re, MM 



þ

QM MM þ EvdW ∑∑ rQM QM=MM , MM

Z

q

coordinate variable, ξ, is constrained around particular values.45 The values of the variables sampled during the simulations are then pieced together to construct a distribution function from which the PMF is obtained as a function of the distinguished reaction coordinate (W(ξ)). The PMF is related to the normalized probability of finding the system at a particular value of the chosen coordinate by eq 4 Z ð4Þ WðξÞ ¼ C  kT ln FðrN ÞδðξðrN Þ  ξÞ drN  1

where the superscripts indicate the value of the reaction coordinate at the reactants (R) and TS and Gξ(ξR) is the free energy associated with setting the reaction coordinate to a specific value at the reactant state. Normally this last term makes a small contribution47 and the activation free energy is directly estimated from the PMF change between the maximum of the profile and the reactant’s minimum ΔG‡ ðξÞ  Wðξ‡ Þ  WðξR Þ ¼ ΔW ‡ ðξÞ

ð6Þ

The selection of the reaction coordinate is usually trivial when the mechanism can be driven by a single internal coordinate or a simple combination (as the antisymmetric combination of two interatomic distances). However this is not the case for the reaction subject of study in this paper where at least 6 coordinates are participating. Instead we were compelled to obtain a much more computationally demanding 2D-PMF using two coordinates: ζ1 and ζ2. The 2D-PMF is related to the probability of finding the system at particular values of these two coordinates Z Wðζ1 , ζ2 Þ ¼ C0  kTln FðrN Þδðζ1 ðrN Þ

ð3Þ

This interaction energy can be exactly decomposed in a sum over residues provided that the polarized wave function (Ψ) is employed to evaluate this energy contribution. The global polarization effect can be obtained from the gas phase energy difference between the polarized, Ψ, and unpolarized, Ψ0, wave functions. The Potential Energy Surface (PES). The fDYNAMO library38 was used to explore the different PESs as a function of the distances R1, R2, R3, R4, R5, and R6 (see Scheme 1). All PESs were obtained using two antisymmetric combinations of these distances: the R1-R2 corresponds to the proton transfer from Asp298 to the linking anomeric oxygen, R3-R4 corresponds to the attack of the 2-acetamido carbonyl oxygen on the anomeric center and departure of aglycone leaving group, and R5-R6 corresponds to the proton transfer from linking N-acetyl amide hydrogen to Asp297. The Potential of Mean Force (PMF). In order to obtain the free energy associated to the proton transfers and cyclization, we have traced several PMF.4143 All PMFs have been calculated using the weighted histogram analysis method (WHAM) combined with the umbrella sampling approach44,45 as implemented in fDYNAMO.38 The procedure for the PMF calculation is straightforward and requires series of molecular dynamics simulations in which the distinguished reaction

ð5Þ

 ξ1 Þδðζ2 ðrN Þ  ξ2 Þ drN  2

ð7Þ

To estimate the activation free energy from this quantity we recovered one-dimensional PMF changes tracing a maximum probability reaction path on the 2D-PMF surface and integrating over the perpendicular coordinate. In constructing the 2D-PMF, the distinguished reaction coordinates were the antisymmetric combination of the distances R1-R2 and R3-R4, corresponding to proton transfer from Asp298 to linking anomeric oxygen and coordinates corresponds to attack of the 2-acetamido carbonyl oxygen on the anomeric center, respectively. A total of 39 simulations were performed at different values of R1-R2 (ranging from 2.0 to þ2.0 Å), with an umbrella force constant of 2800 J 3 mol1 3 Å2 applied to this distinguished reaction coordinate. In addition, 30 simulations were performed at different values of R3-R4 (ranging from 1.5 to þ1.5 Å), also with an umbrella force constant of 2800 J 3 mol1 3 Å2 on this antisimmetric combination of distances. Consequently, 1170 simulation windows were needed to obtain the 2D-PMF. The values of the variables sampled during the simulations were then pieced together to construct a full distribution function from which the 2D-PMF was obtained. In each window, 5 ps of relaxation were followed by 10 ps of production with a time step of 0.5 fs due to the nature of the chemical step involving a hydrogen transfer. The Verlet algorithm was used to update the velocities. 6767

dx.doi.org/10.1021/jp202079e |J. Phys. Chem. B 2011, 115, 6764–6775

The Journal of Physical Chemistry B

ARTICLE

Figure 1. Time evolution of total substrateprotein interaction energy and R2 distance during the 1 ns of QM/MM MD simulations on reactants state.

Three monodimensional PMFs were also computed being the distinguished reaction coordinates the antisymmetric combination of the distances R1-R2, R3-R4, or R5-R6 coordinates. A total of 60 simulations were performed at different values of R1-R2 (in a range from 2.0 to þ1.0 Å), R3-R4 (in a range from 2.0 to þ2.0 Å) and R5-R6 (in a range from 1.5 to þ1.5 Å), with an umbrella force constant of 2800 kJ 3 mol1 3 Å2. On each window the protocol of the simulation was as described above for 2D-PMF.

’ RESULTS AND DISCUSSION Substrate Distortion and PES. The substrate distortion during the enzymatic hydrolysis of glycosides has been related in some reports.4853 Recently, the reaction mechanisms of different glycoside hydrolases that follow different conformational itineraries have been related by Vocadlo and Davies.1 On the other hand, Brameld and Goddard III have examined the plausible conformations for hexaNAG substrate bound to the active site of Chitinase by means of classical MD simulation.54 The authors proposed that the hydrolysis mechanism of Chitinase A involves substrate distortion and that the protonation of the linking anomeric oxygen requires a boat conformation for the GlcNac residue at binding subsite 1.54 In this proposed mechanism, the first step of the reaction starts with the substrate in the boat conformation.12,20 In order to analyze the substrate distortion, 1 ns of QM/MM MD was carried out using the chair conformation as starting point. After 400 ps of MD simulation a change conformation from chair to boat was observed but the conformation of the ring was back to chairlike after 550 ps, a conformation that is maintained until 800 ps of the MD simulations, where a change conformation from chair to boat was observed again. This result indicates that both chair and boat conformations are close in energy, and as a consequence, significant population of both conformers can be observed at 300 K. This is reflected in the time-evolution of the interaction energy of the substrate with the enzyme and of some key intermolecular distances (see Figure 1, panels a and b). Analysis of time evolution of distances along the MD shows that R2, the hydrogen bond distance from protonated Asp298 to the oxygen atom of the leaving group, is shorter when the substrate adopts the boat conformation than in the chair conformation (see Figure 1b). This relative movement between Asp298 and the substrate is associated with a change on the charge of the ether oxygen that varies from 0.55 au (chair conformation) to 0.61 au (boat conformation). Shorter hydrogen bond distances leads to higher interaction energies, as deduced by combining Figure 1, panels a and b. These results

show that the substrate must adopt a boat conformation for the reaction to proceed, in accordance with results reported in literature.12,20 This conformation will be the starting point to explore the different PESs described below. The first PES (Figure 2a) was obtained using two antisymmetric combinations; R1-R2 and R3-R4. R1-R2 corresponds to the proton transfer from Asp298 to linking anomeric oxygen, while R3-R4 corresponds to the attack of the 2-acetamido carbonyl oxygen on the anomeric center and cleavage of the leaving group (see Scheme 1). First conclusion that can be derived from this PES is that while R structure corresponds to reactive structure in the boat conformation, the product of this step, INT on the PES, corresponds to an oxazoline ion/methanol complex; a structure where the proton of the N-acetyl amide has not been still transferred to the Asp297 (see Scheme 1). It can be also observed that formation of INT takes place in a concerted but very asynchronous way; first the proton transfer from Asp298 to O-GlcNAc and then cyclization can be started via attack of the 2-acetamido carbonyl oxygen on the anomeric center. In the TS structure, TS1, the bond distance for R1 corresponds to 1.90 Å, whereas R2 corresponds to 0.99 Å, indicating a very advanced stage for the proton transfer. On the other hand, the bond distance for R3 in the TS corresponds to 1.91 Å, whereas R4 corresponds to 2.52 Å, indicating that CO bond breaking, R3, and approximation of the 2-acetamido carbonyl oxygen on the anomeric center, R4, are in an earlier stage of the process. The energy barrier at the AM1/MM level was 146.3 kJ mol1 (35.0 kcal/mol). These observations are in agreement with results of the Brameld and Goddard for Chitinase, who, using ab initio gas phase calculations with a reduced model, obtained a spontaneous glycosidic bond cleavage with subsequent formation of an oxazoline ion/methanol complex from a geometry optimization of the protonated linking anomeric oxygen of O-GlcNac in its boat conformation.54 A PES from INT to P was generated using just the R5-R6 antisymmetric combination of distances, that corresponds to the proton transfer to Asp297, as the distinguished reaction coordinate (Figure 2b). The resulting potential energy path render an exothermic chemical reaction with an energy barrier of 50.3 kJ mol1 (12.0 kcal/mol) thus confirming the viability of a stepwise mechanism kinetically determined by the first step. In order to explore a possible concerted mechanism, a PES was constructed using the R3-R4 and the R5-R6 antisymmetric combinations of distances. The proton transfer from Asp298 to the linking anomeric oxygen is not explicitly shown in this PES, but analysis of the PES structures reveals that it takes place spontaneously when these coordinates are used as distinguished reaction coordinates thus being the product of the reaction an 6768

dx.doi.org/10.1021/jp202079e |J. Phys. Chem. B 2011, 115, 6764–6775

The Journal of Physical Chemistry B

ARTICLE

Figure 2. Potential energy surfaces calculated at AM1/MM level obtained by means of (a) R1-R2 and R3-R4 antisymmetric combinations of distances, (b) R5-R6 antisymmetric combination, and (c) R3-R4 and R5-R6 antisymmetric combinations of distances. Values of contour potential energy lines are reported in kJ mol1 and coordinates in Å.

Figure 3. (a) 2D-PMF (in kJ mol1) obtained for proton transfer from Asp-298 to substrate (R-R2, in Å) and attack of the 2-acetamido carbonyl oxygen on the anomeric center (R3-R4, in Å). (b) PMF corresponding to proton transfer from the substrate to Asp-297 (R5-R6, in Å). Values on isoenergetic lines are reported each 10 kJ mol1.

oxazoline intermediate, P. The R5 and R6 distances of the TS located on this PES (Figure 2c) are 1.00 and 2.04 Å, respectively, and R1 and R2 distances, defining the position of the proton to be transferred from Asp298, are 0.99 and 2.05 Å. These results indicate that the proton from N-acetyl amide, as well as the proton from Asp298, are not yet transferred. The values found for R3 and R4 are 2.44 and 1.79 Å, respectively. All together, these results indicate that in this mechanism the proton transfer from linking N-acetyl amide hydrogen to Asp297 occur after the cycle formation in this mechanism. R structure corresponds to reactive structure in the boat conformation. The energy barrier between R and TS in the PES calculated at AM1/MM level was 246.6 kJ.mol1 (58.9 kcal/mol). This result suggests that the reaction could also proceed by a concerted mechanism, although with an energy barrier higher than the rate limiting step of the stepwise mechanism. Potentials of Mean Force. Based on the conclusions derived from the PESs presented above, only the free energy profile of the stepwise mechanism will be studied. The first step will be the proton transfer from Asp298 to linking anomeric oxygen with the concomitant attack of the 2-acetamido carbonyl oxygen on the anomeric center, and the second step would be the proton transfer from linking N-acetyl amide hydrogen to Asp297. In order to study the first step that involves forming and breaking of four bonds, a 2D-PMF using the same coordinates as the ones employed to generate the PES of Figure 2a was computed. The

topology of the resulting 2D-PMF surface, depicted in Figure 3a, results to be quite similar to the corresponding PES of Figure 2a, where the products of this step is the oxazoline ion/methanol intermediate, INT. In order to check this result, free downhill MD simulations taking TS1 structure (protonated GlcNAc) as starting point have been performed showing that the timeevolution of the system leads to spontaneous glycosidic bond cleavage with subsequent formation of an oxazoline ion-methanol complex. We have also performed 500 ps of MD simulation taking INT structure as starting point and the result shows that the oxazoline ion intermediate is a stable structure at AM1/MM level, in agreement with the work of Brameld and Goddard III.54 Structures of R in Figure 3a correspond to reactive structures initially in a chair conformation that change to a boat conformation during the MD simulations at 300 K. These structures can then evolve toward TS1. In TS1, the averaged bond distance for R1 was found to be 1.65 Å, whereas for R2 was 1.02 Å, indicating that the proton is already almost completely transferred. On the other hand, the averaged bond distance for R3 was 2.02 Å, whereas for R4 was 2.27 Å. Other averaged key distances obtained from the 2D-PMF are showed in Table 1. Using as starting point the INT structure (ion oxazoline/ methanol), the second PMF (Figure 3b) was traced using R5-R6 as distinguished coordinate, which corresponds to the proton transfer from linking N-acetyl amide nitrogen atom to one of the oxygen atoms of the carbocylate group of Asp297. The product 6769

dx.doi.org/10.1021/jp202079e |J. Phys. Chem. B 2011, 115, 6764–6775

The Journal of Physical Chemistry B

ARTICLE

Table 1. Averaged Distances of Key Distances for Rchair, Rboat, TS1, INT, TS2, and P Structures (in Å)a

a

distance

Rchair

Rboat

TS1

INT

TS2

P

R1

0.98 (0.03)

0.98 (0.03)

1.65 (0.04)

5.14 (1.15)

4.96 (0.65)

5.45 (0.53)

R2

4.06 (0.56)

3.61 (0.65)

1.02 (0.03)

0.97 (0.03)

0.98 (0.01)

0.98 (0.01)

R3

1.42 (0.03)

1.43 (0.03)

2.02 (0.15)

6.60 (1.46)

4.95 (0.51)

5.30 (0.64)

R4

2.95 (0.13)

2.97 (0.14)

2.27 (0.14)

1.51 (0.04)

1.48 (0.03)

1.46 (0.03)

R5

1.00 (0.03)

1.00 (0.03)

1.00 (0.03)

1.01 (0.03)

1.25 (0.03)

2.09 (0.03)

R6

2.08 (0.13)

2.08 (0.11)

2.08 (0.13)

1.99 (0.13)

1.27 (0.03)

0.98 (0.02)

Standard deviations, computed during the last 10 ps of the MD simualtions, are reported in brackets (in Å) from 2D-PMF.

found for this profile is, indeed, the oxazoline/methanol complex, P in Scheme 1. The calculated activation free energies, ΔG‡, for TS1 and TS2 are 195 and 40 kJ mol1, respectively (46.6 and 9.6 kcal mol1). On the base of these results, it is reasonable to suggest that the first step, the proton transfer and glycosidic bond cleavage occurring simultaneously, would be the rate limiting step of the process. A note of caution must be introduced at this point, since these barriers appear to be overestimated keeping in mind that they correspond to an enzyme catalytic process. Obviously, the low level of the AM1 semiempirical Hamiltonian, required to allow performing long hybrid QM/MM MD simulations, must be the main source of error and absolute energetic values must be considered at relative terms. Thus, corrections at B3LYP(6-31G(d))/MM level have been carried out by adding the single-point gas phase energy differences between the relative energies obtained at a high level (B3LYP) and a low level (AM1) of theory, using the structures from the AM1/MM potential energy surface. These corrections are then incorporated into the simulations using a two-dimensional cubic spline function.55 The new free energy surface render a barrier for the first step of 155 kJ mol1 (37.0 kcal mol1). This value, while still being too high for the kind of process under study, is significantly lower than the one obtained at AM1/MM level. The remaining error could be due to the use of AM1 geometries in our high level calculations. Interaction Energy. A problem of central importance in computational biology is the quantitative determination of absolute binding affinities in diverse and complex systems. Predicting the binding free energy of ligands to macromolecules can have great practical values in identifying novel molecules that can bind to target receptors and act as therapeutic drugs.56 In order to get a deep insight into de proteinsubstrate interactions, 500 ps of MD were run for all substrateprotein complexes appearing along the free energy profile (Rchair, Rboat, TS1, INT, TS2, and P), constraining the corresponding distinguished reaction coordinate in the two transition states. A note of caution has to be mentioned at this point since the methanol formed in the first step of reaction is eliminated from the system after formation of INT. Thus, methanol is not included in QM/ MM MD simulations to compute proteinsubstrate interactions for INTenzyme, TS2enzyme, and Penzyme. Averaged substrateprotein interaction energies for all these complexes at AM1/MM level are listed in Table 2. It must also be pointed out that Asp298 and Asp297 were included in the QM region in all of these calculations to obtain comparable results of substrateprotein total interaction energies for all the stationary states, whatever the degree of proton transfer to or from the enzyme. As a consequence, the reported values in Table 2 are not substrateprotein interaction energies, strictly speaking, but the

Table 2. Averaged AM1/MM ProteinSubstrate Interaction Energies for Rchair, Rboat, TS1, INT, TS2, and P Structures Obtained from the PMFs (in kJ mol1)a structures

EAM1/MM (kJ mol1)

Rchair

1186.9 ( 71.7

Rboat TS1

1260.6 ( 78.2 1622.6 ( 27.4

INT

1520.3 ( 35.9

TS2

1605.2 ( 36.5

P

1215.1 ( 56.3

a

Values reported for INT, TS2 and P have been computed treating the formed methanol molecule in the MM region.

analysis of the trend along the reaction progress can be used to obtain an idea of the most adequate structures to design new inhibitors. Analysis of the data reported in Table 2 reveals that the structures presenting the highest total substrateprotein interaction energy, in absolute value, are TS1 and TS2. This conclusion confirms these two structures as the candidates to be used as a mold to design new inhibitors, which is reasonable keeping in mind that most of the inhibitors are transition state analogs (TSA). The origin of the substrateenzyme interaction can be analyzed decomposing the total interaction energy in a sum of contributions due to each of the enzymatic residues. Averaged proteinsubstrate interaction energies by residue, computed including now both aspartate residues (Asp297 and Asp298) in the MM region, are displayed in Figure 4. In this figure, positive values correspond to unfavorable interactions, whereas a negative values means that the interaction is stabilizing the substrate in that particular state. The pattern of interactions of Rchair enzyme complex is very close to the one obtained for the Rboatenzyme complex. Both show the interaction with Asp401 as the most important one for both chair and boat conformation. This result is in agreement with the observations of Rao et al. demonstrating that Asp401 has significant effect on the ability of the enzyme to bind the substrate/inhibitor (5 10-fold increase in Km).14 Other important residues interacting with reactants are Tyr-335 and Asn-396. These two residues would then be well positioned to stabilize any partial negative charge that develops at the glycosidic oxygen in the TS, consistent with the role of these residues assigned for the catalysis.14 However, these interactions can be also important for substrate fixation in both of its conformations. The major difference between the Rboatenzyme and Rchairenzyme complex is in interaction with the Asp298 residue. In fact, as observed in Figure 4 Asp 298 presents favorable interaction in Rboatenzyme complex, where this interaction is not appreciable 6770

dx.doi.org/10.1021/jp202079e |J. Phys. Chem. B 2011, 115, 6764–6775

The Journal of Physical Chemistry B

ARTICLE

Figure 4. Contributions of individual amino acids to substrateprotein interaction energy (in kJ mol1) computed for Rchair, Rboat, TS1, INT, TS2, and P.

in Rchairenzyme complex. This explains the higher interaction energies obtained when substrate adopts a boat conformation (see Figure 1a). The pattern of interactions of Penzyme complex reveals very weak interactions, in accordance with the total value of the substrateprotein interaction energy as reported in Table 2. INTenzyme complex presents a pattern of interactions very close to the one observed for TS2enzyme but with weaker interaction with Asn396 and unfavorable interaction with Lys’184. Moreover, the unfavorable interaction with Lys218, observed in both INT and TS2, is weaker in the later. Thus, our analysis will be focused on the transition states TS1 and TS2. A detailed analysis of the contribution of individual residues to the total TS1enzyme and TS2enzyme interaction energy shows that while the strongest interaction of TS1 is established with Asp298 residue, while in the case of TS2, it is with Asp297 residue. The other significant difference between these structures is the interaction with Lys218 residue, which is favorable for TS1 but unfavorable for TS2. Interaction with Lys184’ stabilizes both transition states, and Asn396 presents favorable interactions not only in these two complexes but also, as observed previously, in both reactants conformers. The fact that Asp297, Asp298, together with Asp401, stabilize all species created along the reaction profile is in agreement with mutagenesis and structural studies reported for O-GlcNacase from bacterial Clostridium perfringens by Rao et al. showing that Asp297, Asp298 and Asp401 are the key catalytic residues of O-GlcNAcase.14 These findings are also in agreement with our previous QM/ MM studies of interactions of two inhibitors (PUGNAc and NAG-thiazoline) with CpGH84H protein.24 Nevertheless, in the present report, analysis of contribution of individual residues to the total proteinsubstrate interaction energy has been carried

out along the reaction profile showing that Asp298 has important stabilizing interaction with TS1 complex, Asp297 with TS2 and oxazoline ion/methanol complex, and Asp401 presents favorable interaction with the substrate in all states along the reaction. Inhibitors are generally designed to resemble the structure of a molecular species occurring during the chemical reaction and, in particular, transition states structures, if accepting the seminal hypothesis of Pauling that enzymes catalyze reactions by preferentially binding (and stabilizing) the transition state.57 In our case, NAG-thiazoline inhibitor structure is close to TS2enzyme structure, whereas PUGNAc and GlcNAcstatin B inhibitors seem to be closer to TS1enzyme structure. The high resemble between the TS2 structure and NAG inhibitor, and between TS1 and PUGNAG and GlcNAcstatin B can be seen in Figure 5, where superposition of the structures of the two transition states with their corresponding inhibitors are displayed. As mentioned in the Introduction, PUGNAc, NAG-thiazoline and GlcNAcstatin B are well characterized inhibitors of human O-GlcNAcase that present the highest binding values (Ki = 46, Ki = 70, and Ki = 0.4 nM, respectively).22 Table 2 shows that total interaction energies for TS1 and TS2 are also close and, in fact, they are the structures presenting the highest stabilizing values. All these observations suggest that the three studied inhibitors might be considered as TSA, which would explain their potency, by comparison with other ground state analogues based inhibitors (i.e., methyl 2-acetamido-2-deoxy-β-D-glucopyranoside, Ki = 11 μM)21 or simply serendipitous binders (i.e., iminosugar inhibitors as the nojirimycin class)58 As mentioned by Vocadlo and co-workers,21 NAG-thiazoline molecule has an obvious geometrical resemblance to the intermediate of the reaction (P in scheme 1), but we herein are demonstrating that it tightly binds the active site of the protein by its resemblance to the TS2, that is, in turn, close to the oxazoline intermediate. 6771

dx.doi.org/10.1021/jp202079e |J. Phys. Chem. B 2011, 115, 6764–6775

The Journal of Physical Chemistry B

ARTICLE

Figure 5. Detail of the active site of structures obtained in (a) TS1 and (b) TS2. (c) Superposition of TS1 (green) and PUGNAc (silver) obtained from crystal structure of CpGH84H-PUGNAc complex at, PDB code 2CBJ. (d) Superposition of TS1 (green) and GlcNAcstatin (silver) obtained from cpOGA in complex with GlcNAcstatin inhibitor, code PDB 2WB5. (e) Superposition of TS2 (green) and N-butyl-thiazoline (silver), obtained from structure of the B. thetaiotaomicrometer GH84 O-GlcNAcase, in complex with N-butyl-thiazoline.

Figure 6. Map electrostatic potential (MEP) surfaces derived from B3LYP/6-31G* calculations for the Rboat, TS1, INT, TS2, PUGNAc, GlcNAcstatin (in its protonated state), and NAG-thiazoline polarized by the charges of the protein environment in the corresponding states. The increase of negative charges goes from the blue (positive) to red (negative).

Molecular Electrostatic Potential (MEP). The molecular electrostatic potential (MEP) is defined as the interaction energy between the charge distribution of a molecule and a unit positive

charge. The MEP is highly informative of the nuclear and electronic charge distribution of a given molecule. The MEP has been applied to a range of fields, such as the study of 6772

dx.doi.org/10.1021/jp202079e |J. Phys. Chem. B 2011, 115, 6764–6775

The Journal of Physical Chemistry B biological interactions, topographical analysis of the electronic structure of molecules, and definition of molecular reactivity patterns.5962 The potential applications of the MEP as a tool for interpretations and prediction of chemical reactivity, as well as to probe the similarity of transition states and inhibitors, has long been recognized.63 Recently, electrostatic potentials have been calculated for PUGNAc and NAG-thiazoline.24 Vocadlo and coworkers carried out a similar study exploring MEP of PUGNAc and NAG-thiazoline and three species found along the gas phase reaction path of cyclization step involving O-GlcNAcase and OGlcNac: ground state, transition state model and oxazoline.21 Three dimensional MEP surfaces for Rboat, TS1, INT, TS2, P, PUGNAc, NAG-thiazoline, and GlcNAcstatin B are depicted in Figure 6. Our results provide a more realistic picture since the reaction profile is more complete and the electronic density distribution of the substrate along the reaction path is polarized, and computed, under the effect of the protein environment. The MEP surfaces were derived from B3LYP/6-31G (d) calculations using Gaussian 0364 and visualized in Gaussview 3.07.65 These surfaces correspond to an the isodensity value of 0.002 au. It is interesting to highlight the noncovalent interactions occurring at the molecular surface. In the figure, the most nucleophilic regions (negative electronic potential) are in red, and the most electrophilic regions (positive electrostatic potential) appear in blue. TS1 displays large positive (blue) regions matching the positions of Asp-298, and TS2 displays large negative (red) regions matching the positions of Asp-297. The negative electrostatic potential for TS1 can be found around the oxygen atoms, mainly on carbonyl oxygen. This oxygen is ready to attack on the anomeric center to form a covalent bicyclic oxazoline. The positive electrostatic potential for TS1 and INT can be found around the hydrogen atoms, mainly on atom H linked on nitrogen atom, this hydrogen is ready to be transferred to Asp297 that can afford electrostatic stabilization of the TS1 and INT (oxazoline ion) and acts as a general base which accepting the proton from INT. The pKa of the nitrogen atom changes during the cyclization step from around 15 to about 5.5,21 a result that is in agreement with the differences of the MEP between Rboat, TS1, INT, TS2, and P. Analysis of the MEPs displayed in Figure 6 reveals that PUGNAc is closer to TS1 than TS2. The negative electronic potential for PUGNAc found on carbonyl oxygen is also found on carbonyl oxygen of TS1 and Rboat. In the MEP of both PUGNAc and TS1, a positive potential is found in the region of hydrogen linked on nitrogen atom of N-acetamido group. In accordance with these observations, we can say that the MEP of PUGNAc is more like MEP of TS1 than TS2 or INT. On the other hand, Glycoside hydrolases use alternative conformational pathways involving the half-chair 4H3 or boat 2,5B/B2,5 transition state conformation.1 In this study, the TS1 obtained adopts 4H3 conformation in agreement with literature.1 Other important characteristic of TS1 is the double bond character of the C1O5 bond. The value of this distance (1.32 Å) indicates an sp2hybridized center installed at C1 of the pyranose ring in agreement with previous works.1,13 A double bond character in the C1O5 is also observed in PUGNAc that adopts 4E conformational.15 In accordance with these results PUGNAc mimic both the geometrical and electrostatic features of TS1. Nevertheless, although PUGNAc resembles TS1, better inhibitors could be proposed. TS1 display large positive (blue) regions matching the positions of Asp-298, On the other hand, PUGNAc display large negative (red) in this region. In previous report, we

ARTICLE

have proposed that the oxime nitrogen atom of PUGNAc is not able to accept a proton from Asp-298.31 PUGNAc would be a better inhibitor if its oxime nitrogen atom was protonated. In this regard, analysis of GlcNActatin B, significantly more potent inhibitor than PUGNAc and NAG-thiazoline inhibitors of human O-GlcNAcase, revealed a tight interaction with catalytic acid (Asp-298) and the (presumably protonated) imidazole.22 Heightman and Vasella had suggested that a lateral protonation of the exocyclic nitrogen atom of the imidazole ring should account for the excellent inhibiting of these compounds.66 The protonation of the exocyclic nitrogen atom of the imidazole ring leads to positive potential in this region matching positions of Asp-298. As observed in Figure 6, the MEP of protonated GlcNActatin B is clearly closer to the MEP of TS1 than PUGNAc, which would explain its higher inhibition activity. The MEP for the last of the studied inhibitors, NAG-thiazoline, reveals that it is close to TS2 and P. NAG-thiazoline adopts 4 C1 conformation likely TS2 and P. It is important to point out that MEP of TS2 was computed without H of amide group since the R5 (distance between the hydrogen and nitrogen atoms of the N-acetyl amide group, see scheme 1) in TS2 corresponds to 1.25 Å. Thus we can say that NAG-thiazoline might mimic both electrostatic and conformational features of TS2 and P (oxazoline intermediate) indicating that it is likely to be a mimic of species that lies on the reaction coordinate. However, following the discussion of Heightman and Vasella66 on the activity dependency on the protonation state of the exocyclic nitrogen atom of the imidazole ring, a protonation state of NAG-thiazoline should be studied in the future, since NAGthiazolinium ion might be more likely to oxazolinium ion (INT). It is important to point out that this new structure appeared as a result of an exhaustive hybrid QM/MM exploration of the first step of the full catalytic mechanism used by O-GlcNAcase; a step that has appeared to be not concerted but stepwise (the oxazoline intermediate described in the literature correspond to P in our complete reaction mechanism).

’ CONCLUSIONS We have applied QM/MM MD techniques to study the first step of the catalytic mechanism used by O-GlcNAcase to hydrolyze O-GlcNAc. The free energy profile shows that the formation of the oxazoline intermediate in the O-GlcNAcase catalytic reaction takes place by means of a stepwise mechanism. The first step would be a cyclization of the acetomide group, which seems to be dependent on the proton transfer from a conserved aspartate, Asp298 in Clostridium perfringens O-GlcNAcase. From this new intermediate, a proton is transferred from the azoline ring to another conserved aspartate (Asp297) thus forming the oxazoline ion and departure of the aglycone (a methanol molecule in our model). The obtained results are consistent with experimental and previous theoretical results, where the cyclization proceeds via attack of the 2-acetamido carbonyl oxygen on the anomeric center to form a covalent bicyclic oxazoline intermediate. Nevertheless, the presence of a new high-energy intermediate and a second transition state in the catalyzed pathway could better explain some of the experimental observations. O-GlcNAc is distorted to a boat conformation before a proton from Asp298 was transferred to the linking anomeric oxygen. Analysis of contribution of individual residues to the total proteinsubstrate 6773

dx.doi.org/10.1021/jp202079e |J. Phys. Chem. B 2011, 115, 6764–6775

The Journal of Physical Chemistry B interaction energy shows that Asp297 or Asp298 and Asp401 have important interaction with substrate in all the states along the reaction but, in particular, in the transition states (TS1 and TS2) and the oxazolyne intermediate (INT). The global interaction energy between protein and ligand shows that TS1 and TS2 present the highest binding affinity, the latter being close to the values observed for the intermediate. These observations suggest that the most efficient inhibitors for OGlcNAcase catalyzed glycoside hydrolysis, NAG-thiazoline, PUGNAc, and the GlcNAcstatin derivatives, might be considered as TSA, in contrast with previous predictions. Our results would explain the similar potency of NAG-thiazoline and PUGNAc inhibitors, by comparison with other ground state analogues based inhibitors or simply adventitious binders to the enzyme. NAG-thiazoline molecule has an obvious geometrical resemblance to P (the oxazoline intermediate described in the literature), but we herein are demonstrating that it tightly binds the active site of the protein by its resemblance to the TS2, that is, in turn, close to the oxazoline intermediate. GlcNAcstatin B, which is the most potent human O-GlcNAcase inhibitor reported so far, presents the closest geometrical and electronic structure to the TS1, in agreement with our main statement that in order to design efficient inhibitors of a certain enzyme, these must resemble the TS of the catalyzed reaction as much as possible. Thus, the present theoretical study not only is important to get an insight into the reaction, but the knowledge of the molecular mechanism can be used to get structural and electronic information of the species appearing along an enzymatic reaction path, which is of great value for a systematic protocol for synthesis of inhibitors.

’ ASSOCIATED CONTENT

bS

Supporting Information. Complete ref 64. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*(V.M.) E-mail: [email protected]. Phone: (þ34) 964-728084. Fax: (þ34) 964-728066. (C.N.A.) E-mail: [email protected]. Phone: (þ55) 91-32018235. Fax: (þ55) 91-32017633.

’ ACKNOWLEDGMENT We thank the Spanish Ministry Ministerio de Ciencia e Inovacion for Project CTQ2009-14541-C2, Universitat Jaume I - BANCAIXA Foundation for Projects P1 3 1B2008-36, P1 3 1B2008-37, and P1 3 1B2008-38, Generalitat Valenciana for Prometeo/2009/053 project and SEUI for financial support of a Hispano-Brasile~ no collaboration project (PHB2005-0091-PC). The authors also acknowledge the Servei dInformatica, Universitat Jaume I for generous allotment of computer time. J.L. thanks CAPES for their financial support and the warm hospitality during the research stay at Departament de Química Física i Analítica, Universitat Jaume I. V.M. thanks the Spanish Ministry Ministerio de Educacion for traveling financial support, Project PR2009-0539 ’ REFERENCES (1) Vocadlo, D. J.; Davies, G. J. Curr. Opin. Chem. Biol. 2008, 12, 539–555.

ARTICLE

(2) Torres, C. R.; Hart, G. W. J. Biol. Chem. 1984, 259, 3308–3317. (3) Dong, D. L; Hart, G. W. J. Biol. Chem. 1994, 269, 19321–19330. (4) Wells, L.; Vosseller, K.; Hart, G. W. Science 2001, 291, 2376–2378. (5) Akimoto, Y.; Hart, G. W.; Hirano, H.; Kawakami, H. Med. Mol. Morphol. 2005, 38, 84–91. (6) Vosseller, K.; Wells, L.; Lane, M. D.; Hart, G. W. Proc. Natl. Acad. Sci. 2002, 99, 5313–5318. (7) Shafi, R.; Iyer, S. P. N. L.; Ellies, G.; O’Donnell, N.; Marek, K. W.; Chui, D.; Hart, G. W.; Marth, J. D. Proc. Natl. Acad. Sci. 2000, 97, 5735–5739. (8) Arnold, C. S.; Johnson, G. V.; Cole, R. N.; Dong, D.; Lee, M.; Hart, G. W. J. Biol. Chem. 1996, 271, 28741–28744. (9) Butters, T. D.; Dwek, R. A.; Platt, F. M. Adv. Exp. Med. Biol. 2003, 535, 219–226. (10) von Itzstein, M; Wu, W-Y; Kok, G B; Pegg, M S; Dyason, J C; Jin, B; Van Phan, T; Smythe, M L; White, H F; Oliver, S W; Colman, P M; Varghese, J N; Ryan, D M; Woods, J M; Bethell, R C; Hotham, V J; Cameron, J M; Penn, C R. Nature 1993, 363, 418–423. (11) Mendel, D. B.; Tai, C. Y.; Escarpe, P, A.; Li, W.; Sidwell, R. W.; Huffman, J. H.; Sweet, C.; Jakeman, K. J.; Merson, J.; Lacy, S. A.; Lew, W.; Williams, M. A.; Zhang, L.; Chen, M. S.; Bischofberger, N.; Kim, C. U. Antimicrob. Agents Chemother. 1998, 42, 640–646. (12) Coutinho P. M., Henrissat B. Carbohydrate-Active Enzymes server at URL: http://afmb.cnrs-mrs.fr/CAZY, 1999. (13) C-etinbas, N.; Macauley, M. S.; Stubbs, K. A.; Drapala, R.; Vocadlo, D. J. Biochemistry 2006, 45, 3835–3844. (14) Toleman, C.; Paterson, A. J.; Kudlow, J. E. Biochim. Biophys. Acta 2006, 1760, 829–839. (15) Rao, F. V.; Dorfmueller, H. C.; Villa, F.; Allwood, M.; Eggleston, I.; van Aalten, D. M. EMBO J. 2006, 25, 1569–1578. (16) Sinnott, M. L. Chem. Rev. 1990, 90, 1171–1202. (17) Davies, G. J.; Sinnott, M. L.; Withers, S. G. Comprehensive Biological Catalysis; Sinnot, M. L., Ed.; Academic Press: San Diego, CA, 1998; Vol 1, pp 119209. (18) Sinnott, M. L.; Souchard, I. J. Biochem. J. 1973, 133, 89–98. (19) Greig, I. R.; Williams, I. H. Chem. Commun. 2007, 3747–3749. (20) Macauley, M. S.; Whitworth, G. E.; Debowski, A. W.; Chin, D.; Vocadlo, D. J. J. Biol. Chem. 2005, 280, 25313–25322. (21) Whitworth, G. E.; Macauley, M. S.; Stubbs, K. A.; Dennis, R. J.; Taylor, E. J.; Davies, G. J.; Greig, I. R.; Vocadlo, D. J. J. Am. Chem. Soc. 2007, 129, 635–644. (22) Dorfmueller, H. C.; Borodkin, V. S.; Schimpl, M.; van Aalten, D. M. F. Biochem. J. 2009, 420, 221–227. (23) Dorfmueller, H. C.; Borodkin, V. S.; Schimpl, M.; Shepherd, S. M.; Shpiro, N. A.; van Aalten, D. M. F. J. Am. Chem. Soc. 2006, 128, 16484–16485. (24) Lameira, J.; Alves, C. N.; Moliner, V.; Martí, S.; Kanaan, N.; Tu~ non, I. J. Phys. Chem. B 2008, 112, 14260–14266. (25) Wolfenden, R. Nature 1969, 223, 704–705. (26) Ma, B.; Kumar, S.; Tsai, C.; Hu, Z.; Nussinov, R. J. Theor. Biol. 2000, 203, 383–397. (27) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227–249. (28) Roca, M.; Martí, S.; Andres, J.; Moliner, V.; Silla, E.; Tu~ non, I.; Bertran, J. Chem. Soc. Rev. 2004, 33, 98–107. (29) Alves, C. N.; Martí, S.; Castillo, R.; Andres, J.; Moliner, V.; Tu~ non, I.; Silla, E. Chem.—Eur. J. 2007, 13, 7715–7724. (30) Antosiewicz, J.; McCammon, J. A.; Gilson, M. K. J. Mol. Biol. 1994, 238, 415–436. (31) Field, M. J.; Amara, P.; David, L.; Rinaldo, D. Personal Communication (32) Hui, L.; Robertson, A. D.; Jensen, J. H. Proteins 2005, 61, 704–721. (33) Lameira, J.; Alves, C. N.; Moliner, V.; Martí, S.; Castillo, R.; Tu~ non, I. J. Phys. Chem. B 2010, 114, 7029–7036. (34) Byrd, R. H.; Lu, P.; Nocedal, J.; Zhu, C. J. Sci. Comp. 1995, 16, 1190–1208. (35) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F. J. Am. Chem. Soc. 1985, 107, 3902–3909. 6774

dx.doi.org/10.1021/jp202079e |J. Phys. Chem. B 2011, 115, 6764–6775

The Journal of Physical Chemistry B

ARTICLE

(36) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225–11236. (37) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926–935. (38) (a) Field, M. J. A practical Introduction to the Simulation of Molecular Systems; Cambridge University Press: Cambridge, U.K., 1999. (b) Field, M. J.; Albe, M.; Bret, C.; Proust-de Martin, F.; Thomas, A. J. Comput. Chem. 2000, 21, 1088–1100. (39) (a) Zhang, Y.; Liu, H.; Yang, W. J. Chem. Phys. 2000, 112, 3483–3492. (b) Martí, S.; Moliner, V.; Tu~non, I. J. Chem. Theory Comput. 2005, 1, 1008–1016. (40) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361–373. (41) Kirkwood, J. G. J. Chem. Phys. 1935, 3, 300–313. (42) McQuarrie, D. A. Statistical Mechanics; Harper & Row: New York, 1976; p 266. (43) Roux, B. Comput. Phys. Commun. 1995, 91, 275–282. (44) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187–199. (45) Roux, B. Comput. Phys. Commun. 1995, 91, 275–282. (46) Schenter, G. K.; Garrett, B. C.; Truhlar, D. G. J. Chem. Phys. 2003, 119, 5828–5833. (47) Roca, M.; Moliner, V.; Ruiz-Pernía, J. J.; Silla, E.; Tu~ non, I. J. Phys. Chem. A 2006, 110, 503–509. (48) Sulzenbacher, G.; Driguez, H.; Henrissat, B.; Sch€ulein, M; Davies, G. J. Biochemistry 1996, 35, 15280–15287. (49) Tews, I.; Perrakis, A.; Oppenheim, A.; Dauter, Z.; Wilson, K. S.; Vorgias, C. E. Nat. Struct. Biol. 1996, 3, 638–648. (50) Davies, G. J.; Mackenzie, L.; Varrot, A.; Dauter, M.; Brzozowski, A. M.; Sch€ulein, M.; Withers, S. G. Hydrolase Biochem. 1998, 37, 11707–11713. (51) Money, V. A.; Smith, N. L.; Scaffidi, A.; Stick, R. V.; Gilbert, H. J.; Davies, G. J. Angew. Chem. Int. Ed. 2006, 45, 5136–5140. (52) Tailford, L. E.; Offen, W. A.; Smith, N. L.; Dumon, C.; Morland, C.; Gratien, J.; Heck, M. P.; Stick, R. V.; Bleriot, Y.; Vasella, A.; Gilbert, H. J.; Davies, G. J. Nat. Chem. Biol. 2008, 4, 306–312. (53) Numao, S.; Kuntz, D. A.; Withers, S. G.; Rose, D. R. J. Biol. Chem. 2003, 278, 48074–48083. (54) Brameld, K. A.; Goddard, W. A., III J. Am. Chem. Soc. 1998, 120, 3571–3580. (55) Ruiz-Pernía, J. J.; Silla, E.; Tu~non, I.; Martí, S. J. Phys. Chem. B 2006, 110, 17663–17670. (56) Wlodawer, A. Annu. Rev. Med. 2002, 53, 595–614. (57) Pauling, L. Nature 1948, 161, 707–709. (58) (a) Stutz, A. E. Iminosugars as Glycosidase Inhibitors, Nojirimycin and Beyond; John Wiley and Sons: Weinheim, Germany, 1999. (b) Berland, C. R.; Sigurskjold, B. W.; Stoffer, B.; Frandsen, T. P.; Svensson, B. Biochemistry 1995, 34, 10153–10161. (59) Gavezotti, A. J. Am. Chem. Soc. 1983, 105, 5220–5224. (60) Murray, J. S.; Brinck., T.; Lane, P.; Paulsen, K.; Politzer, P. THEOCHEM 1994, 307, 55–64. (61) van Duyne, G. D.; Standaert, R. F.; Karplus, P. A.; Schreiber, S. L.; Clardy, J. Science 1991, 252, 839–842. (62) Cubero, E.; Luque, F. J.; Orozco, M. Proc. Natl. Acad. Sci. 1998, 95, 5976–5980. (63) Bagdassarian, C. K.; Schramm, V. L.; Schwartz, S. D. J. Am. Chem. Soc. 1996, 118, 8825–8836. (64) Frisch, M. J.; Gaussian 03, revision C.02; Gaussian, Inc.; Wallingford, CT, 2004. (65) Dennington, R.; Keith, T.; Millam, J.; Eppinnett, K.; Hovell, W. L.; Gilliland, R. GaussView, V3.0; Semichem Inc.: Shawnee Mission, KS, 2003. (66) Heightman, T. D.; Vasella, A. T. Angew. Chem., Int. Ed. 1999, 38, 750–770.

6775

dx.doi.org/10.1021/jp202079e |J. Phys. Chem. B 2011, 115, 6764–6775