Molecular Mechanical Molecular Dynamics

Apr 29, 2010 - Structures of O-GlcNAcase-PUGNAc complexes obtained after 1.5 ns of AM1/MM MD simulations: (a) wild-type, form A; (b) wild- type, form ...
1 downloads 0 Views 2MB Size
Download PDF
J. Phys. Chem. B 2010, 114, 7029–7036

7029

Quantum Mechanical/Molecular Mechanical Molecular Dynamics Simulation of Wild-Type and Seven Mutants of CpNagJ in Complex with PUGNAc Jeronimo Lameira,†,‡ Cla´udio Nahum Alves,*,† Vicent Moliner,*,‡,§ Sergio Martı´,‡ Raquel Castillo,‡ and In˜aki Tun˜o´n| Laborato´rio de Planejamento e DesenVolVimento de Fa´rmacos, Instituto de Cieˆncias Exatas e Naturais, UniVersidade Federal do Para´, CP 11101, 66075-110, Bele´m, PA, Brazil, Departament de Quı´mica Fı´sica i Analı´tica, UniVersitat Jaume I, 12071 Castello´n, Spain, Institute of Applied Radiation Chemistry, Technical UniVersity of Lodz, 90-924 Lodz, Poland, and Departament de Quı´mica Fı´sica, UniVersitat de Vale`ncia, 46100 Burjassot, Vale`ncia, Spain ReceiVed: December 6, 2009; ReVised Manuscript ReceiVed: March 5, 2010

The enzyme O-glycoprotein 2-acetamino-2-deoxy-β-D-glucopyranosidase (O-GlcNAcase) is responsible for the removal of N-acetylglucosamine moieties from 2-acetamido-2-deoxy-β-D-glucopyranose (O-GlcNAc) residues of serine/threonine residues of modified proteins. We herein present results of hybrid quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations applied to the study of the interactions established between a bacterial Clostridium perfringens homologue (CpNagJ) and PUGNAc, a potent known inhibitor of this enzyme. Electrostatic binding free energy and energy term decomposition have been computed for the wild-type CpNagJ and several mutants: D297N, D298N, Y335F, N390A, N396A, D401A, and W490A. The theoretical results have been compared with recently experimental data based on crystallographic and mutation studies on the same system. Our results reveal that, first, there is a strong interaction between Asp401, Asp298, and Asp297 residues and the PUGNAc inhibitor; and, second, the electrostatic substrate binding free energy is higher in wild-type, N390A and W490A mutants than in D297N, D298N, Y335F, N396A, and D401A ones, in accordance with the experimental results. Finally, both our theoretical predictions and the experimental data are compatible with a substrate-assisted reaction mechanism, involving two conserved aspartate residues. Introduction O-Glycoprotein 2-acetamino-2-deoxy-β-D-glucopyranosidase (O-GlcNAcase)1,2 is a member of the family of glycoside hydrolases-84 (GH84)3 that is able to cleave 2-acetamido-2deoxy-β-D-glucopyranose (O-GlcNAc) from post-translationally modified protein serine and threonine residues.1-3 It has been proposed that this enzyme uses a catalytic mechanism involving substrate-assisted catalysis.1-3 O-GlcNAc is abundant in the brain, particularly on cytoskeletal proteins,4,5 which has supported the hypothesis that disruptions of O-GlcNAc proteins may contribute to certain neurodegenerative disorders such as Parkinson dystonia6 and Alzheimer’s disease.5 O-GlcNAc is also present in many transcription regulatory proteins. Recently, it has been shown that O-(2-acetamido-2-deoxy-D-glucopyranosylidene)amino-N-phenylcarbamate (PUGNAc, see Scheme 1) is an effective inhibitor of O-GlcNAcase in skeletal muscle, dramatically reducing insulin sensitivity in adipocytes.7 To gain insight into which amino acid residues are crucial for substrate recognition and catalysis, the study of the crystal structure of a GH84, from Clostridium perfringens (CpNagJ), and also the crystal structure in complex with the inhibitor PUGNAc were first reported by Rao and co-workers, with * To whom correspondence should be addressed. (V.M.) E-mail: [email protected]. Telephone: (+34) 964-728084. Fax: (+34) 964-728066. (C.N.A.) E-mail: [email protected]. Telephone: (+55) 91-32018235. Fax: (+55) 91-32017633. † Universidade Federal do Para´. ‡ Universitat Jaume I. § Technical University of Lodz. | Universitat de Vale`ncia.

significant sequence homology to the human O-GlcNAcase.8 Rao and co-workers8 have also performed mutagenesis experiments for seven mutants of CpNagJ, and the results showed that the residues Asp401, Asp297, Asp298, Tyr335, Trp490, and Asn396 are important to substrate/inhibitor binding and catalysis, while the mutation of Asn390 showed kinetic parameters indistinguishable from those of the wild-type enzyme. Rao and co-workers suggest that PUGNAc inhibits OGlcNAcase by mimicking the transition state of the OGlcNAcase-catalyzed hydrolysis of N-acetylglucosaminide by virtue of its sp2 anomeric C1, similar to the oxocarbenium ionlike transition state.8 On the other hand, Whitworth and coworkers have suggested that NAG-thiazoline is a transition state analogue for O-GlcNAcase, while PUGNAc is a poor transition state analogue.9 In a previous paper,10 we employed hybrid quantum mechanics/molecular mechanics (QM/MM)11 molecular dynamics (MD) simulations to determine electrostatic binding free energy of CpNagJ complexed with both PUGNAc and NAG-thiazoline. We found that NAG-thiazoline-enzyme complex has a binding free energy, which is close to the value calculated for the PUGNAc-enzyme complex. The small difference observed between these two inhibitors was in agreement with experimental data of inhibitory activity. Inhibitors of O-GlcNAcase activity have been developed trying to mimic the geometrical and/or electronic properties of structures appearing during the reaction process, either transition structures or intermediates. In this sense, the main aim of this study is to investigate structural and dynamical properties of CpNagJ complexed with PUGNAc ligand, which will be useful for designing inhibitors of O-GlcNAcase. We used hybrid QM/

10.1021/jp9115673  2010 American Chemical Society Published on Web 04/29/2010

7030

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

Lameira et al.

SCHEME 1: Schematic Representation of the Two Models Used to Model the CpNagJ-PUGNAc Complex: A Form, Where the Asp298 is Protonated, and B Form, Where the Oxime Nitrogen of PUGNAc is Protonated

MM MD simulations to study the details of the CpNagJ-PUGNAc interaction for wild-type and D297N, D298N, Y335F, N390A, N396A, D401A, and W490A CpNagJ mutants. Mutations represent a fascinating experimental challenge, and theoretical simulations are invited to disclose still unexplored features on enzyme systems. In this regard, our computational protocol allows a detailed analysis of the interactions of PUGNAc inhibitor with the key residues inside the binding pocket, as well as an estimation of the inhibitor-protein electrostatic binding free energy that can be used to explain the difference in inhibitory potency of the wild-type and CpNagJ mutants. Computational Methods As successfully done in previous studies based on MD simulations within a combined QM/MM approach,10,12 the ligand/substrate species may be described by a QM model, while the protein and solvent environments are represented by MM force fields. This hybrid methodology 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.13 Moreover, as the largest part of the system is described classically, enough sampling can be obtained at reasonable computational cost. In the work reported here, the computational model for the QM/MM MD calculations was taken from the trimeric crystal structure of the CpNagJ-PUGNAc complex, with PDB code 2CBJ.8 Since standard pKa values of ionizable groups can be shifted by local protein environments,14 an accurate assignment of the protonation states of all these residues at pH ) 7 was carried out by recalculating the standard pKa values of the titratable amino acids using the “cluster method”,14 as implemented by Field and co-workers.15 According to this method, each titratable residue in the protein is perturbed by the electrostatic effect of the protein environment. As a result, most residues were found at their standard protonation state, except for the Asp298 residue that seems to be protonated at this pH. Interestingly, Asp298 has been proposed to be the catalytic acid, protonating the glycosidic bond in the first step of the OGlcNAcase-catalyzed hydrolysis of N-acetylglucosaminides, while Asp297 would be a general base in a typical acid/base catalysis.16 Therefore, in this report, we have modeled the

PUGNAc-enzyme complex protonating the Asp298 (A form in Scheme 1) or protonating the oxime nitrogen of PUGNAc (B form in Scheme 1). Complex B would resemble a more advanced stage of the reaction process, with the leaving group being protonated. After adding the hydrogen atoms to the full structure, series of optimization algorithms (steepest descent conjugated gradient and L-BFGS-B17) were applied. Then the optimized protein was placed in a cubic box of pre-equilibrated waters (80 Å side), using the center of mass of the 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 the same optimization algorithms. Afterward, 100 ps of hybrid QM/MM Langevin-Verlet MD at 300 K and in a canonical thermodynamic ensemble (NVT) were used to equilibrate de model. For the hybrid QM/MM MD calculations, we have two molecular models of simulation, as explained above. In the first model, A form in Scheme 1, the atoms of the inhibitor and Asp298 were selected to be treated by QM, while in the second model, B form in Scheme 1, only the atoms of the inhibitor were included in the QM region. To saturate the valence of the QM/MM frontier in the first model, we used a link atom,18,19 placed between the CR and Cβ atoms of the Asp298 residue. The number of QM atoms resulted then to be 53 and 45 for models A and B, respectively, while the final system contains a total of 23 486 atoms in both models. In both cases, the semiempirical AM1 Hamiltonian was employed to describe the QM part,20 while the rest of the system was described using the OPLS-AA21 and TIP3P22 force fields for protein and water molecules, respectively, as implemented in the DYNAMO library.23 Hybrid QM/MM MD simulations for wild-type and D297N, D298N, Y335F, N390A, N396A, D401A, and W490A CpNagJ mutants complexed with PUGNAc were performed for both models. Due to the amount of degrees of freedom, any residue 20 Å apart from any of the atoms of the initial inhibitor was selected to be kept frozen (which represents 19 247 atoms) in the remaining calculations in order to make the model computationally feasible. Cutoffs for the nonbonding interactions were applied using a switching scheme, within a radius range from 14.5 to 16 Å. Once the systems were pre-equilibrated, 1500 ps of QM/MM MD were run at a temperature of 300 K. The computed RMSD

QM/MM MD Simulation of Wild-Type and Mutant CpNagJ for the protein during the last 400 ps renders a value always below 0.9 Å, while 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 potential energy of our scheme is derived from the standard QM/MM formulation:

(∑〈 | | 〉 ∑∑ )

ˆ o |ψ〉 + EQM/MM ) 〈ψ|H

ψ

qMM ψ + re,MM

ZQMqMM vdW + EQM/MM + EMM rQM,MM

elect vdW EQM/MM ) Evac + EQM/MM + EQM/MM + EMM

(1)

(2)

vdW where EMM is the energy of the MM subsystem, EQM/MM is the van der Waals interaction energy between the QM and MM subsystems, Evac is the gas phase energy of the polarized QM subsystem, 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). It is important to point out that in order to take advantage of powerful optimizations algorithms (see for example ref 24), especially when higher level QM hamiltonians were used, the electrostatic interaction term can be substituted by a pure coulombic one, in which the charges of the QM atoms are obtained by means of electrostatic potential fit methods based on the electronic density.25 The interaction energy between the inhibitor and the environment, computed by residue, 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 (i) of the protein is given by the following expression:

Int elect vdW EQM/MM,i ) EQM/MM,i + EQM/MM,i ) qMM ZQMqMM Ψ Ψ + + r r e,MM MM∈i QM QM,MM σQM,MM 12 σQM,MM 4εQM,MM rQM,MM rQM,MM MM∈i QM

∑

∑

[〈 | | 〉 ∑ ] ∑ [( ) ( ) ] 6

(3)

Finally, in order to evaluate the electrostatic contribution to the interaction free energy, series of QM/MM MD simulations have been carried out, introducing a λ parameter in the electrostatic QM/MM interaction:

(∑〈 | | 〉 )

ˆ o |ψ〉 + λ EQM/MM(λ) ) 〈ψ|H

∑∑

ψ

qMM ψ + re,MM

ZQMqMM vdW + EQM/MM + EMM rQM,MM

(4)

Then λ smoothly changes between values of 1, corresponding to a full MM charge-wave function interaction, and 0, whereas no electrostatic interaction with the force field is introduced. Note that the van der Waals interaction is always calculated. The calculation of the free energy difference of two consecutive

J. Phys. Chem. B, Vol. 114, No. 20, 2010 7031 windows is performed by means of free energy perturbation (FEP) methods, and then the total free energy change is evaluated as the sum of all the windows covering the full transformation from the initial to the final state:

∆GElect-QM/MM ≈ -

1 β

∑ ln〈e-β[E(λ i

i+1)

- E(λi)]

〉λi

(5)

We have used a total of 100 windows, from λ ) 0 (no electrostatic interaction) to λ ) 1 (full interaction), which turns into a δλ of 0.01. The procedure can be applied both forward (i f i + 1) and backward (i + 1 f i) as long as the comparison provides information about the convergence of the process. In each window, a total of 5 ps of relaxation followed by 10 ps of production QM/MM MD has been performed using the NVT ensemble at the reference temperature of 300 K. A note of caution has to be mentioned at this point, since, for computing the electrostatic contribution to the free energy of interaction, as well as to evaluate the interaction energy between the inhibitor and the environment by residue, QMMM partitioning of the system has been done including in the QM region just the inhibitor in both models A and B (Asp298 was treated by MM in the A model for these calculations). Moreover, it must be pointed out that, since the use of high level Hamiltonians for performing hybrid MD is actually prohibitive, the electrostatic free energy has been evaluated using the semiempirical AM1 Hamiltoninan. In this regard, the B3LYP functional, within the 6-31G(d) basis set, has been employed to perform single-point QM/MM internal energy calculations in order to evaluate the goodness of the AM1/MM potentials to describe the nonbonded interactions between the inhibitors and the protein. In this case, and since electrostatic calculations strongly depend on the cutoff selected to compute the nonbonding electrostatic interactions in charged systems, the whole system within the minimum image convention has been used and the same protocol has been applied for the AM1/MM both in the A and B forms (see the Supporting Information). In order to obtain the free energy change associated to the proton transfer between the oxime nitrogen of PUGNAc hydrogen and the carboxylate oxygen atom Asp298 (i.e., the transformation between A and B forms), we have computed the corresponding potential of mean force (PMF),26,27 using the weighted histogram analysis method (WHAM) combined with the umbrella sampling approach28-30 as implemented in DYNAMO.23 In this PMF, the distinguished reaction coordinates were the antisymmetric combination of the distances HE2-OD2 (R1) and N15-HE2 (R2) (Scheme 1) that describes the proton transfer from Asp298 to oxime nitrogen hydrogen. A total of 60 simulations were performed at different values of R1 and R2 (in a range from -1.5 to +1.5 Å), with an umbrella force constant of 2800 kJ · mol-1 · A-2. In each window, 5 ps of relaxation was followed by 10 ps of production with a time step of 0.5 ps due to the nature of the chemical step involving a hydrogen transfer. The Verlet algorithm was used to update the velocities at 300 K. Results and Discussion As explained in a previous section, hybrid QM/MM MD simulations of 1.5 ns for wild-type and D297N, D298N, Y335F, N390A, N396A, D401A, and W490A CpNagJ mutants complexed with PUGNAc have been performed. Averaged key substrate-protein distances for wild type and mutants, computed over 100 000 configurations from the last 100 ps of the QM/ MM MD simulations, are listed in Tables 1 and 2 for models

7032

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

Lameira et al.

TABLE 1: Key Interatomic Distances (in Å) between A Form of PUGNAc (I) and Some Residues for Wild-Type and D297N, D298N, Y335F, N390A, N396A, D401A, and W490A Mutantsa

WT D297N D298N Y335F N390A N396A D401A W490A a

I-D298

I-D298

I-D297

I-D401

I-D401

I-Y335

I-Y335

I-N396

N15-HE2 3.64 (0.55) 6.14 (0.32)

O18-HE2 2.70 (0.38) 4.75 (0.53)

H29-OD2 2.16 (0.31)

H39-OH 5.01 (0.37) 5.51 (0.24) 4.04 (0.47)

4.05 4.89 6.38 6.17 2.60

H27-OD2 1.86 (0.19) 2.46 (0.26) 1.78 (0.14) 1.77 (0.14) 1.85 (0.17) 1.91 (0.23)

N15-HH 2.18 (0.24) 2.83 (0.32) 2.40 (0.35)

4.20 3.70 4.24 5.79 4.24

H26-OD1 2.23 (0.51) 3.77 (0.32) 3.59 (1.36) 2.00 (0.29) 1.83 (0.16) 4.19 (1.14) 2.40 (0.42)

1.90 (0.20)

2.73 2.76 2.41 2.21

3.23 3.85 3.56 5.11

O3-HD21 2.07 (0.17) 2.08 (0.19) 2.15 (0.22) 2.14 (0.24) 2.06 (0.18)

(0.31) (0.52) (0.48) (0.72) (0.57)

(0.36) (0.50) (0.41) (0.70) (0.40)

2.35 2.06 2.13 2.21 1.93 2.19

(0.62) (0.26) (0.28) (0.32) (0.18) (0.33)

(0.33) (0.40) (0.34) (0.22)

(0.35) (0.50) (0.39) (0.32)

2.11 (0.20) 2.15 (0.19)

Standard deviations (in Å) computed during the last 100 ps are reported in parentheses.

TABLE 2: Key Interatomic Distances (in Å) between B Form of PUGNAc (I) and Some Residues for Wild-Type and D297N, D298N, Y335F, N390A, N396A, D401A, and W490A Mutantsa

WT D297N D298N Y335F N390A N396A D401A W490A a

I-D298

I-D297

I-D401

I-D401

I-Y335

I-Y335

I-N396

HE2-OD2 1.62 (0.12) 1.81 (0.30)

H29-OD2 1.76 (0.13)

H26-OD1 1.88 (0.18) 1.85 (0.17) 1.79 (0.16) 1.87 (0.26) 3.17 (1.14) 1.86 (0.18)

H27-OD2 1.87 (0.21) 1.80 (0.15) 1.75 (0.13) 1.79 (0.14) 2.41 (0.26) 1.81 (0.17)

N15-HH 5.17 (0.18) 3.49 (0.33) 3.62 (0.51)

H39-OH 2.14 (0.25) 4.70 (0.62) 4.14 (0.54)

1.83 (0.15)

1.80 (0.15)

5.34 3.84 3.75 3.45

4.48 4.26 3.49 4.89

O3-HD21 2.13 (0.24) 2.21 (0.23) 2.28 (0.27) 2.16 (0.23) 2.06 (0.18)

1.60 1.62 1.82 1.86 1.75

(0.11) (0.13) (0.23) (0.58) (0.19)

1.96 1.90 1.96 1.94 1.97 2.00

(0.20) (0.18) (0.22) (0.19) (0.25) (0.23)

(0.29) (0.31) (0.51) (0.43)

(0.49) (0.45) (0.66) (0.45)

2.25 (0.27) 2.19 (0.26)

Standard deviations (in Å) computed during the last 100 ps are reported in parentheses.

A and B, respectively. Representative snapshots of the averaged structures of wild-type and the D401A mutant (which is the one presenting the largest structural deviations with respect the wild-type) obtained from the 1.5 ns MD simulations are shown in Figure 1 for both A and B forms. As mentioned in a previous section, contributions of individual residues to the total inhibitor-protein interaction energy have been computed and averaged during the last 100 ps of MD simulations and the results are shown in Figure 2 for both models A and B. In this figure, negative values correspond to stabilizing effects. Electrostatic contributions to the interaction free energy obtained from the electrostatic FEP are reported in Table 3 for all studied systems. As mentioned in the Computational Methods section, B3LYP/MM calculations have been carried out in order to evaluate the goodness of the AM1/MM potentials to describe the nonbonded interactions between the inhibitors and the protein. Comparative analysis of the total interaction energies for A and B forms and the contributions of individual residues computed at both levels of theory are displayed in the Supporting Information. It is important to point out that the pattern of interactions that can be deduced from Figure 2 is qualitatively equivalent to the corresponding calculations at higher level DFT/ MM, thus confirming the reliability of the AM1 semiempirical Hamiltonian to describe this kind of interactions for stable molecules in their ground state. Wild-Type CpNagj-PUGNAc Complex. As can be observed in Figure 1 and Tables 1 and 2, the proton of oxime nitrogen of PUGNAc interacts through a hydrogen bond with OD2 oxygen atom of Asp298 in B form (See Scheme 1 for numbering). Meanwhile, in A form, this hydrogen bond interaction is not maintained and the unprotonated oxime nitrogen of PUGNAc interacts with the Tyr335 residue. In this model A, Asp298 presents a weaker hydrogen bond interaction with the oxygen atom of the N-phenyl carbamoyl group of PUGNAc, O18 (2.70 Å), instead. Interestingly, Cleveland and co-workers showed how the N-acetylglucosamino-1,5-lactone oxime (LOGNAc), which lacks the N-phenyl carbamoyl group of PUG-

NAc, binds 30-fold more poorly than PUGNAc,4 providing an experimental indication of the importance of this interaction. The strong interaction with Asp298 in B form (see Figures 1b and 2b) can be one of the major reasons of the stabilization of the PUGNAc-enzyme complex as can be observed by comparing with A form, where this interaction does not exist (see Figures 1a and 2a), in accordance with the dramatic diminution of the electrostatic interaction free energies from B form to A form (see Table 3), as will be explained afterward. These results are also in accordance with the observation that the protonated nitrogen engaged in an adventitious electrostatic interaction with the general acid/base catalytic residue,10 and gives importance to the N-phenyl carbamoyl group in PUGNAc. Correspondingly, the same argument can be used to explain the electrostatic repulsion observed between positive Lys218 and the PUGNAc inhibitor in B form, which is not in A form (see Figure 2a and b). X-ray crystallographic data8 are, in general, not far from both structures A and B. However, the interaction between Asp298 and PUGNAc is better described in the B model, while the interaction between Tyr335 and PUGNAc is better described in the A one, as can be seen in Figure 3a and b where overlay of the PUGNAc-CpNagJ complex from X-ray crystallographic data,8 and A and B model, respectively, is depicted. As it can be seen, the position of residues Asp298 and Tyr335 deviates with respect to X-ray crystallographic data in models A and B, respectively.8 In fact, the distance between OD2 of Asp298 and oxime nitrogen of PUGNAc is 2.88, 3.33, and 2.66 Å for X-ray crystallographic structure, A model, and B model, respectively. Meanwhile, the distance between OH of Tyr335 and oxime nitrogen of PUGNAc is 2.68, 3.06, and 4.54 Å for X-ray crystallographic structure, A model, and B model, respectively. As observed in Figure 2a and b, other important residues presenting a favorable interaction with the ligand are Asp297, Asn396, and Asp401. Asp297 presents hydrogen bond interactions with nitrogen of amide in A and B forms, with this interaction being responsible for stabilizing the conformation

QM/MM MD Simulation of Wild-Type and Mutant CpNagJ

J. Phys. Chem. B, Vol. 114, No. 20, 2010 7033

Figure 1. Structures of O-GlcNAcase-PUGNAc complexes obtained after 1.5 ns of AM1/MM MD simulations: (a) wild-type, form A; (b) wildtype, form B; (c) D401A, form A; and (d) D401A, form B.

of the acetamido group. Regarding Asp401, it presents two hydrogen bond interactions with O4 and O6 hydroxyl oxygen atoms of the ligand in both A and B forms. Interestingly, these results would be in agreement with experimental measurements of Rao et al., who observed that mutation of Asp401 has significant effects on both the ability of the enzyme to bind the inhibitor and catalysis.8 Lys218 forms a hydrogen bond with the O3 sugar hydroxyl oxygen atom in A form, in agreement with X-ray structure,8 which is not found in B form. Asn396 forms a hydrogen bond with the oxygen atom of acetamido group in both forms, which is reflected in the favorable interaction observed in Figure 2a and b. Tyr335, Trp394, and Trp490 seem to be stabilizing the cavity created by the three aforementioned residues. Finally, due to the change of substrate charge from model A to B, Arg224 (a residue which is more than 6.46 Å away from the inhibitor) presents an electrostatically unfavorable interaction with the inhibitor in B form but not in A form. A similar argument could be applied to Lys184′ behavior: a residue belonging to a different protein chain that is interacting with the carbonyl oxygen atom of PUGNAc in A form but is more than 5.84 Å from the inhibitor in B form. This change in the relative position of Lys184′ is reflected in the inhibitor stabilization, as observed in Figure 2. Mutant CpNagj-PUGNAc Complex. As mentioned in the Introduction section, in order to get a deeper insight into the role of the different residues of the CpNagj active site, seven

different mutations have been set up and the same computational protocol as the one performed for the wild-type enzymePUGNAc complex has been followed. After the QM/MM MD simulations, we can observe that in all seven CpNagJ mutants the inhibitor remains bound in the enzyme active site, which is in accordance with data reported by Rao and co-workers.8 The analysis of the structural results reveals that the largest difference between mutants and the wild-type enzyme was observed for D401A mutant. The absence of the interaction with Asp401 provokes some changes in the interactions established with other residues (see Figures 1c,d and 2c,d). The lack of the hydrogen bond interaction between Asp401 and O4 and O6 hydroxyl oxygen atoms of the ligands decreases the favorable interaction between this residue and the inhibitor, as deduced when comparing Figure 2a with 2c, or 2b with 2d. The averaged distance between Asp298 and PUGNAc increases dramatically from the wild-type to D401A in A form, while it remains almost invariant in the B form case. On the other hand, the interactions between the inhibitor and Asp297 are not perturbed after D401A mutations in both forms. Therefore, the loss of ligand-protein interaction in D401A CpNagJ, mainly through Asp401 and Asp298 residues, may be responsible for the significant consequences of the mutation on the ability of the enzyme to bind the inhibitor and on the catalytic efficiency.8 Mutation of Asp297 (D297N) provokes a weaker interaction with Asp401, Asp298, Tyr335, and Asn396 residues in A form,

7034

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

Lameira et al. TABLE 3: Binding Free Energies (∆GElect - QM/MM) Obtained from Electrostatic FEP Calculations at the AM1/MM Level (Values in kJ/mol) for A and B Formsa ∆GElect -QM/MM wild-type D297N D298N Y335F N390A N396A D401A W490A

A

B

kcat/Km (µM-1 · s-1)

IC50 (nM)b

-142.1 -113.6 -130.3 -125.6 -150.3 -120.7 -111.8 -149.5

-367.4 -315.9 -291.7 -362.2 -373.6 -388.2 -252.7 -348.1

3.6 0.001 0.0003 0.0001 3.6 0.0002 0.0003 0.1

8.4 ND ND ND 10.4 ND ND 24.0

a Experimental ligand-protein kinetic data (kcat/Km in µM-1 · s-1) and biological activities (IC50 in nM) reported for comparison purposes (obtained from ref 8). b ND ) not determined.

the inhibitor and Asp297, Asp298, Tyr335, and Asn396 residues (see Tables 1 and 2) causes a less favorable protein-ligand interaction which can be the origin of the lower ability of the mutants to bind the inhibitor and to catalyze the reaction (see the Supporting Information). Finally, the structures of N390A and W490A mutants do not differ dramatically from the structure of wild-type enzyme, and consequently, the kinetic parameters are almost indistinguishable for N390A and have minor effects for W490A.8 Electrostatic Free Energy Calculations. The ratio between the catalytic rate constant (kcat) and the Michaelis constant (Km) can be related to the binding free energy of the transition state. According to Scheme 2, it is possible to write the following:

-RT ln

Figure 2. Contributions of individual amino acid residues to inhibitor interaction energy (in kJ/mol) averaged over the last 100 ps of AM1/ MM MD simulations of optimized structures obtained after 1.5 ns of AM1/MM MD simulations for the different O-GlcNAcase-PUGNAc structures: (a) wild-type, form A; (b) wild-type, form B; (c) D401A, form A; and (d) D401A, form B.

as shown by the increase of the interatomic distances (see Table 1) and decreased interaction energies by residues (see the Supporting Information). The most dramatic effect of this mutation on B form is a weaker interaction with Asp298 and Tyr335. Similar trends were observed in D298N, Y335F, and N396A mutants: the lack of interactions established between

kcat ‡ TS ) ∆Guncat + ∆Gbind Km

(6)

This thermodynamic cycle was already applied by Warshel and co-workers31 to rationalize the binding of a series of phenolate ions to the active site of kesteroid isomerase (KSI) by comparison of experimental data obtained by Kraut and co-workers.32 For a series of mutant proteins catalyzing the same reaction, the variation in the ratio kcat/Km should be related to differences in the binding free energy of the TS. If the substrate used (PUGNAc in our case) is a good mimic of the real TS, the previous relationship should also hold for the binding free energy of this transition state analogue, TSA. Table 3 shows the electrostatic contribution to the proteinsubstrate interactions, free energies obtained from the FEP calculations at AM1/MM level, for A and B forms, together with experimental kinetic and biological activities values, for all systems under study. Obviously, the free energy estimated in our FEP calculations is just one of the electrostatic contributions to the total binding free energy. Other terms such as the electrostatic contributions to the desolvation of the substrate and the protein can be important. Nevertheless, the first one is constant for the series of experimental data analyzed here, but the desolvation of the protein can make important differences among mutants with different polarity or charge in the active site. The first conclusion that can be derived from the calculations is that the protein-ligand interactions obtained at AM1/MM and DFT/MM levels render the same tendencies (see the Supporting Information for a graphical comparative analysis). The protein-ligand electrostatic binding free energies obtained at AM1/MM level are listed in Table 3 for the wild-type and all mutants. As observed, absolute values obtained in model B are dramatically larger than the ones obtained in model A. These values suggest that the interaction between the negatively charged Asp297 and Asp298 with the positively charged inhibitor

QM/MM MD Simulation of Wild-Type and Mutant CpNagJ

J. Phys. Chem. B, Vol. 114, No. 20, 2010 7035

Figure 3. (a) Overlay of the PUGNAc-CpNagJ complex from X-ray crystallographic data8 (gray) with (a) PUGNAc-CpNagJ complex from model A (in green) and (b) PUGNAc-CpNagJ complex from model B (in green).

SCHEME 2: Thermodynamic Cycle, Where Raq, Renz, TSaq, and TSenz, Correspond to Reactants in Solution, Reactants in Enzymatic Site, Transition State in Solution, and Transition State in Enzymatic Site, Respectively

could be responsible for the larger values obtained in model B, a result that does not necessarily reflect a more stable system since the complete analysis requires the simulation of the proton-transfer step between Asp298 and the substrate, as discussed below. Anyway, when comparing the values for the B form, the results for the Y335F, N390A, N396A, and W490A mutants are close to the result of the wild-type. This trend is not observed for A form, where just the N390A and W490A keep the binding free energy contribution close to the wild-type value. These results for the B form are in agreement with the previous geometrical and substrate-protein interaction analysis that reveal small effects observed after these mutations. In fact, by comparison with the experimental values, the mutations N390A and W490A do not change the ability of the enzyme to bind the inhibitor and catalysis. Differences in A form are observed for the D297N, D298N, Y335F, N396A, and D401A mutants of CpNagJ, which have significantly lower values of the electrostatic free energies. Interestingly, it has been observed that mutation of Asp297, Asp298, Tyr335, Asn396, and Asp401 residues has significant effects on both the ability of the enzyme to bind the inhibitor and catalysis.8 The loss of ligand-protein interaction, mainly through Asp401, Asp297, Asp298, Asn396, and Tyr335 residues in PUGNAc complex, may be responsible for the lost of the activity of this inhibitor. These results suggest that design of an effective inhibitor should resemble the TS of the enzyme catalyzed reaction and theoretical calculations can be used to guide its synthesis. By comparison between the theoretical predictions on both forms and the experimentally measured kinetic data, it seems that the free energy calculations are much more reliable for the A form than for the B form. A direct relationship between the binding free energies and kcat/Km can be observed for the A form,

Figure 4. Potential of mean force obtained for proton transfer from Asp298 to PUGNAc. Relative energies in kJ/mol and reaction coordinate (R2-R1) in Å.

which could be used as a guide in the design of new inhibitors. The correlation coefficients obtained for the linear fit of ln(kcat/ Km) versus the free energy contributions are 0.68 and 0.17 for the A and B forms, respectively (see the Supporting Information for linear regression plots). The larger difference between the correlation coefficients found for the B form than for the A form can be interpreted as an argument in favor of the latter as the most representative of the protonation state of the inhibitor in the active site. Protonation State of PUGNAc. It is interesting to note that using a comparison between the molecular electrostatic potentials of PUGNAc and Thizoline inhibitors we concluded in a previous work10 that the form having a protonated oxime nitrogen could be a more reliable model of PUGNAc in the active site of the wild-type enzyme, while the opposite conclusion is reached here using a more complete set of simulations including different mutants. On the basis of the results and observations described above, we decided to trace the PMF corresponding to the proton transfer from the carboxylate oxygen atom of Asp298 to the oxime nitrogen atom of PUGNAc (this is, from A to B). As shown in the PMF profile depicted in Figure 4, the PUGNAc-protein complex is more stable in the A form than in the B form by about 94 kJ/mol. Thus, we conclude that a protonated Asp298 and a deprotonated PUGNAc seem to be a better model for this inhibitor in the active site of CpNagJ. The free energy barrier for the A to B conversion is only 14 kJ/mol.

7036

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

Conclusions The work reported here provides an insight into the activity of PUGNAc inhibitor against human O-GlcNAcase. Hybrid AM1/ MM MD calculations have been carried out for wild-type and D297N, D298N, Y335F, N390A, N396A, D401A, and W490A CpNagJ mutants complexed with PUGNAc using two different forms for the PUGNAc-enzyme complex: protonating the Asp298 (A form) or protonating the oxime nitrogen of PUGNAc (B form). The calculations allow determination of the electrostatic contribution to the protein-ligand interaction free energies. B3LYP/MM calculations have been carried out in order to evaluate the goodness of the AM1/MM potentials to describe the nonbonded interactions between the inhibitors and the protein showing that the higher level calculations render results qualitatively equivalent to the ones obtained by means of AM1/MM. The analysis of individual interactions between the inhibitor and the amino acids of the enzyme active site, averaged along the AM1/MM MD simulation, reveals how the influence of Asp401, Asp297, Asp298, Asn396, and Tyr335 seems to be crucial, with the interactions established between the inhibitor and these residues being especially important. Moreover, it has been found in form A that Asp298 residue interacts with the carbonyl oxygen atom of the N-phenyl carbamoyl group of PUGNAc, which would explain the experimentally observed poor binding of LOGNAc (that lacks this group) by comparison with PUGNAc. Asp298 residue interacts with the protonated oxime nitrogen of PUGNAc in B form, which seems to be crucial to keep a structure that allows other important interactions with the key residues mentioned above but presents also unfavorable interactions with other residues such as Lys218, Arg224, or Lys184′. These latter interactions were not observed in the A form. On the other hand, the orientation of the inhibitor in CpNagJ mutants are close to the one adopted in the wild-type, except for D401A that showed more conformational change for inhibitor in CpNagJ. The structures found for the seven mutants were similar to that of the wild-type but the absence of the interaction with the residues Asp401, Asp297, Asp298, Asn396, and Tyr335 provoked an important loss in the interactions. This may be responsible for the high resistance observed experimentally to the PUGNAc. Finally, by comparison between the theoretical predictions on both forms and the experimentally measured kinetic data, kcat/KM, it seems that the calculated binding free energies are much more reliable when assuming that Asp298 is protonated, form A, than when oxime nitrogen of PUGNAC is supposed to be protonated, form B. The prevalence of the A form has been also confirmed tracing the free energy profile corresponding to the proton transfer between the oxime nitrogen atom of PUGNAc and the carbocylate group of Asp298. Acknowledgment. We thank the Spanish Ministry Ministerio de Ciencia e InoVacio´n for Project CTQ2009-14541, Universitat Jaume I - BANCAIXA Foundation for Projects P1 · 1B2005-13, P1 · 1B2005-15, and P1 · 1B2005-27, Generalitat Valenciana for Prometeo/2009/053 project, and SEUI for financial support of a Hispano-Brasilen˜o collaboration project (PHB2005-0091-PC). The authors also acknowledge the Servei d’Informatica, Universitat Jaume I for generous allotment of computer time. J.L. would like to thank 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 Educacio´n for mobility financial support.

Lameira et al. Supporting Information Available: Structures and contributions of individual amino acid residues to the interaction energy (in kJ/mol) after 1.5 ns of AM1/MM MD simulations for the mutants and linear regression plots for the lineal fit of ln(kcat/Km) versus the electrostatic free energy contributions for both A and B forms. Contributions of individual amino acid residues to the interaction energy computed at DFT/MM level for the different O-GlcNAcase-PUGNAc structures. Graphical comparativeanalysisofAM1/MMversusDFT/MMprotein-ligand total energy interactions for both A and B forms. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Torres, C. R.; Hart, G. W. J. Biol. Chem. 1984, 259, 3308–3317. (2) Wells, L.; Vosseller, K.; Hart, G. W. Science 2001, 291, 2376– 2378. (3) Hanover, J. A. FASEB J. 2001, 15, 1865–1876. (4) Dong, D. L.; Xu, Z. S.; Hart, G. W.; Cleveland, D. W. J. Biol. Chem. 1996, 271, 20845–20852. (5) Arnold, C. S.; Johnson, G. V.; Cole, R. N.; Dong, D.; Lee, M.; Hart, G. W. J. Biol. Chem. 1996, 271, 28741–28744. (6) 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. U.S.A. 2000, 97, 5735–5739. (7) Park, S. Y.; Ryu, J. W.; Lee, W. Exp. Mol. Med. 2005, 37, 220– 229. (8) Rao, F. V.; Dorfmueller, H. C.; Villa, F.; Allwood, M.; Eggleston, I.; van Aalten, D. M. EMBO J. 2006, 25, 1569–1578. (9) 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. (10) Lameira, J.; Alves, C. N.; Kanaan, N.; Moliner, V.; Tun˜o´n, I.; Martı´, S. J. Phys. Chem. B 2008, 112, 14260–14266. (11) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227–249. (12) Alves, C. N.; Martı´, S.; Castillo, R.; Andre´s, J.; Moliner, V.; Tun˜o´n, I.; Silla, E. Chem.sEur. J. 2007, 13, 7715–7724. (13) Gao, J.; Xia, X. Science 1992, 258, 631–635. (14) Antosiewicz, J.; McCammon, J. A.; Gilson, M. K. J. Mol. Biol. 1994, 238, 415–436. (15) Field, M.; David, L.; Rinaldo, D. Personal Communication. (16) C¸etinbas, N.; Macauley, M. S.; Stubbs, K. A.; Drapala, R.; Vocadlo, D. J. Biochemistry 2006, 45, 3835–3844. (17) Byrd, R. H.; Lu, P. H.; Nocedal, J.; Zhu, C. Y. SIAM J. Sci. Comput. 1995, 16, 1190–1208. (18) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1986, 7, 718–730. (19) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11, 700–733. (20) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902–3909. (21) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225–11236. (22) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926–935. (23) (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. (24) (a) Zhang, Y.; Liu, H.; Yang, W. J. Chem. Phys. 2000, 112, 3483– 3492. (b) Martı´, S.; Moliner, V.; Tun˜o´n, I. J. Chem. Theory Comput. 2005, 1, 1008–1016. (25) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361– 373. (26) Kirkwood, J. G. J. Chem. Phys. 1935, 3, 300–313. (27) McQuarrie, D. A. Statistical Mechanics; Harper & Row: New York, 1976; p 266. (28) Roux, B. Comput. Phys. Commun. 1995, 91, 275–282. (29) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. J. Comput. Chem. 1992, 13, 1011–1021. (30) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187–199. (31) Warshel, A.; Sharma, P. K.; Chu, Z. T.; Åqvist, J. Biochemistry 2007, 46, 1466–1476. (32) Kraut, D. A.; Sigala, P. A.; Pybus, B.; Liu, C. W.; Ringe, D.; Petsko, G. A.; Herschlag, D. PLoS Biol. 2006, 4, 0501-0519.

JP9115673