Comparative Theoretical Studies of the Phosphomonoester Hydrolysis

Jun 15, 2010 - We present here the first ONIOM (our own n-layered integrated molecular orbital + molecular mechanics method) studies of a purple acid ...
4 downloads 4 Views 4MB Size
Download PDF
7110

J. Phys. Chem. A 2010, 114, 7110–7116

Comparative Theoretical Studies of the Phosphomonoester Hydrolysis Mechanism by Purple Acid Phosphatases M. Retegan, A. Milet, and H. Jamet* DCM, Equipe Chimie The´orique, UMR CNRS 5250, ICMG, FR CNRS, UniVersite´ J. Fourier, BP. 53, 38041 Grenoble Cedex, France ReceiVed: January 18, 2010; ReVised Manuscript ReceiVed: April 21, 2010

We present here the first ONIOM (our own n-layered integrated molecular orbital + molecular mechanics method) studies of a purple acid phosphatase enzyme. Our study focused on the structures of the red kidney bean PAP (kbPAP) complexed with phosphate and with phenyl phosphate and on the mechanism of the phenyl phosphate hydrolysis by the enzyme. Density functional theory (DFT) calculations were also performed using models of different sizes for comparison purpose. Results show that the inclusion of three histidine residues, His202, His295, and His296, with their protein surrounding, is crucial to properly describe the coordination of the substrates. They induce a conformation with the substrate closer to the nucleophilic µ-hydroxyde bridge. In the mechanistic study, a transition state is stabilized by a strong hydrogen bond between His202 and the leaving group of the substrate. Consequently, a smaller value for the activation energy barrier is obtained from DFT calculations including this histidine to the same calculations without this histidine. Using the ONIOM method, this activation energy barrier is even more reduced. So the mechanism, which considers the hydroxo group bridging the two metal ions as nucleophile, becomes really convincing, contrary to the results obtained with a small model at the DFT level. Introduction Purple acid phosphatases (PAPs) constitute a new class of nonheme metalloenzymes, which catalyze the hydrolysis of a wide range of phosphate esters and anhydrides like ATP under acid conditions (Figure 1a). These enzymes are distinguished from other acid phosphatases by their characteristic purple color, at around 560 nm, due to the presence of a tyrosine residue ligated to the ferric ion. They have been isolated from a variety of sources including mammals, plants, and microorganisms. Many X-ray structures are available as, for example, the resolution of the PAPs from kidney bean (kbPAP),1 sweet potato,2 porcine uterine fluids (uteroferrin, Uf),3 or again rat bone (recRPAP).4 All of the structures exhibit the same dinuclear framework Fe3+-M2+ (with M ) Fe, Zn, or Mn) and seven invariant metal ligands in their active sites. The equivalence of their active site structures has also been demonstrated by metal ion replacement studies. They have shown that the Zn(II)derivatives of PAPs from pig, rat, and human are kinetically indistinguishable from native di-iron forms,5-10 while the diiron derivative of PAP from red kidney bean has properties similar to the native Fe(III)-Zn(II) form.11,12 The Fe(III) site in PAPs is coordinated to the oxygen atom of a deprotonated tyrosine (Tyr167), the nitrogen atom of a histidine (His325), and the oxygen atoms of two aspartic acid residues (Asp135, Asp164). Similarly, the Zn(II) site is coordinated to the oxygen atom of the bridging aspartic acid (Asp164), to the nitrogen atom of two histidines residues (His286, His323), and to an asparagine oxygen atom (Asn201) (Figure 1b). The presence of a hydroxo group bridging the two metal ions was observed in the electron density maps of PAP from pig and rat and in a more recent crystal structure of PAP complexed with sulfate from red kidney bean.4,13,14 ENDOR15 and EXAFS16 studies of * Corresponding author. Telephone: (+33) 4-76-51-43-74. Fax: (+33) 4-76-51-49-46. E-mail: [email protected].

Figure 1. (a) The phosphomonoester hydrolysis catalyzed by the purple acid phosphatase. (b) Schematic representation of the active site of red kidney bean purple acid phosphatase.

the Fe(III)-Fe(II) and Fe(III)-Zn(II) forms of PAP from pig have confirmed the presence of the bridging hydroxide ligand in the enzyme. Despite the availability of structural information, the catalytic mechanism of PAPs remains ambiguous. Experiments with bovine spleen PAP (bsPAP) and a substrate containing a chiral phosphate group showed that the hydrolysis occurs via a nucleophile substitution mechanism with full inversion of configuration at the phosphorus center,17 ruling out the mechanism with a phosphoenzyme intermediate that has been proposed earlier.18 This supports a mechanism in which the substrate is directly attacked by water or by a hydroxide group. To date, the nature of this nucleophile has eluded identification. Studies assign it to the hydroxo group bridging the two metal ions despite the weak nucleophilicity of this group (see Figure

10.1021/jp100478f  2010 American Chemical Society Published on Web 06/15/2010

QM/MM Studies of Purple Acid Phosphatase 1).14,15,19-21 To verify the validity of this hypothesis, theoretical models can be used. Theoretical chemistry can give molecular-level insights into enzymatic reactivity through determinations of intermediates and transition states’ structures and energies. To describe metalloenzymes, a quantum method is required, typically by hybrid density functional theory (DFT). Because of the large computational cost of these methods, a lot of studies have been limited to small models of the active site. Even though these approaches have had considerable success,22 it is not always obvious to define an adequate size model. In addition, the effect of the protein surrounding is not taken into account. The use of quantum mechanics/molecular mechanics (QM/MM) methods has made the inclusion of this effect possible.23 In particular, the ONIOM24 (our own n-layered integrated molecular orbital + molecular/mechanics method) method has been previously employed with great success in studies of different enzymes.25-27 Briefly, in the ONIOM approach, the system is divided into several layers, each one treated at a different level of theory. In this work, we present the first ONIOM studies of PAP. In the crystal structure of (kbPAP) complexed with phosphate,1 three histidines (His202, His295, and His296) are located near to the dimetallic center of kbPAP. These residues that are in positions where they can interact with the free phosphate can be correctly included in ONIOM models. The first part of this Article is devoted to describe the crystal structure of (kbPAP) complexed with phosphate.1 ONIOM calculations on the whole enzyme were performed and compared to DFT calculations on models of different sizes. Subsequently, the mechanism involving the bridging hydroxo group with a phenylphosphate as substrate is analyzed. Large deviations are obtained between the different calculations, showing the necessity of using large models and not only models containing the atoms, which contribute mostly to the reactivity, to describe these structures. Methods Setup of the System and Molecular Mechanics Calculations. All calculations in this work are based on the crystal structure of kbPAP complexed with phosphate (PDB: 4KBP).28 Because the four subunits presented in the crystallographic data are equivalent, only subunit D was included. Hydrogen atoms were placed according to the predicted pKa of the amino acids at the crystallization pH by employing the PROPKA 2.0 software.29 In the case of the three histidines that are close to the active site, His295 was taken protonated, while for the histidines His296 and His202, a detailed discussion is presented in the Results and Discussion because they are hydrogen bonded to the substrate. The system was solvated in a cubic water box of about 93 × 93 × 93 Å3 and was neutralized by randomly adding 13 chlorine ions. All molecular mechanics calculations have been performed using the Gromacs package30 together with the AMBER-94 force field31 for the protein and a TIP3P model32 for water. The missing charges were derived by the ab initio electrostatic surface potential (ESP) calculation (HF/6-31G*), using the RESP methodology.33 The van der Waals parameters for zinc atom are taken from Hoops et al.,34 and for iron atom from Oda et al.,35 but their values have no effect on binding energies. All simulations were done with the amino acids bound to the metal ions, the substrate (phosphate or phenyl phosphate for the mechanistic study) and three histidine residues close to the active site (His202, His295, and His296) fixed to their crystallographic positions. The system was first energy minimized for

J. Phys. Chem. A, Vol. 114, No. 26, 2010 7111 1000 steps using a steepest descent algorithm, and then a 100 ps NVT molecular dynamic simulation was performed with a time step of 2 fs at 300 K. All covalent bonds involving H-atoms were constrained with the LINCS approach.36 A 14 Å cutoff was used for the Lennard-Jones interactions. Coulomb interactions within 10 Å were computed each step, and, beyond this cutoff, the particle mesh Ewald method37 was used with a reciprocal grid spacing of 1.2 Å. The last structure of the molecular dynamic simulation was taken as the initial geometry for the subsequent ONIOM optimizations. Solvation water molecules and counterions introduced for this simulation have been removed in the ONIOM calculations. Consequently, depending on the substrate, the system used for the ONIOM calculations presented in this Article contains 6824 or 6836 atoms. A possibly more accurate approach would have been to run multiple dynamics simulations and systematic ONIOM optimizations on the last molecular dynamics structure. Nevertheless, because our QM region is rather large (6824 and 6836 atoms depending on the substrate), we deemed that the computational cost to perform the ONIOM optimization on different simulations would be prohibitively high. QM and ONIOM Hybrid Calculations. All QM and ONIOM calculations were carried out with the Gaussian 03 program.38 The ONIOM calculations have been performed with the use of a two-layer ONIOM (QM/MM) scheme, in which hydrogen link atoms treat the interface between the QM and MM regions. Electrostatic interactions between the QM and MM layers were calculated by an electronic embedding scheme.39 The MM charges were incorporated into the QM Hamiltonian. Thus, the use of ONIOM method with this scheme takes into account the polarization of the QM part by the MM part. The QM part was described with the B3LYP functional,40 which remains one of the most popular functionals for studying mechanisms of enzymes containing transition metals.41 The atoms were represented with the LANL2DZ basis set,42 whose core parts in the case of iron, zinc, and phosphorus atoms were described by an effective core potential (ECP).43 A set of polarization functions was added for the oxygen (R ) 0.961), nitrogen (R ) 0.736), and phosphorus (R ) 0.364) atoms.44 According to the EPR studies,45 high spin configuration was chosen for the iron, and thus the spin-unrestrictive formalism was used. For the MM region, we used the AMBER-94 force field.31 Because of the large number of degrees of freedom, the position of all amino acids that were more than 12 Å from the phosphorus atom of the substrate was fixed during the optimization. A similar approach was already employed and proved to give reliable results.46 A schematic representation of our ONIOM calculations and the partition of QM and MM regions are depicted in Figure 2. The QM models consist of 113 or 133 atoms (for phosphate or phenyl phosphate as substrate), and all QM-MM borders were placed on a CR-Cβ bond. For the DFT calculations, the same functional and basis set as in the description of the QM layer were used. Two models were built. The first model includes only the amino acids of the first coordination sphere of iron and zinc atoms, while the second also encompasses two of the three histidines close to the active site in the crystal structure of (kbPAP) complexed with phosphate, His202 and His296. The Nε-atom of these histidines interacts through hydrogen bonds with the bound phosphate. For the amino acids, we used the following models for the amino acids: tyrosine was modeled as O-Ph, histidine as imidazole, while the aspartic acid and asparagine as OCOCH3 group and OC(CH3)NH2 groups, respectively. These two models called, respectively, 1DS and 1DL are represented in Figure 3.

7112

J. Phys. Chem. A, Vol. 114, No. 26, 2010

Figure 2. (a) Schematic representation of our ONIOM calculations with a CPK representation of the QM part. (b) Partition of QM and MM regions; the QM-MM borders were placed on CR-Cβ bonds. The atoms of the QM part are given with a ball and stick representation; the MM atoms close to the frontier are given with the licorice representation.

The histidine residues were cut at the CR-Cβ bond, and the free valence was saturated with hydrogen atoms. The optimized structures obtained in the ONIOM calculations were used as initial geometry. In the preliminary calculations on the smaller model (DS) complex with either phosphate or phenyl phosphate, we noticed after optimization a hydrogen bond between atoms H1 and Ob (Figure 3a). This is linked to the fact that this model did not include the residue His296, which interacts with the substrate. To avoid this problem, the dihedral angle H1-O1-P-Ob (see Figure 3a) is frozen to its value obtained in the ONIOM optimizations. The nature of all structures optimized by DFT calculations was characterized by frequency calculations at the same level of theory. Results and Discussion The different labels used in this work follow the subsequent pattern: a number (1 or 2) to differentiate the two substrates studied (phosphate or phenyl phosphate), a letter (D or O) for

Retegan et al.

Figure 3. Schematic representation of the DFT model used in this work: (a) 1DS model, (b) 1DL model.

the level of calculations used (DFT or ONIOM resp), a letter (p or d) for the protonation state of His296 in ONIOM models, and letter (S or L for small and large) for the size of the system in DFT models. For the mechanistic study, an additional letter was used to differentiate the reactant (r) and the transition state (ts). Modelization of the kbPAP Complexed with Phosphate. In the X-ray structure of kbPAP complexed with phosphate,1 we can notice three histidine residues, His202, His295, and His296, close to the active site. They can interact with the free phosphate (Figure 4). To relax the protein and the solvent, a short classical molecular dynamic simulation was carried out prior to the ONIOM optimizations. For the final 20 ps of the simulation, the system is stable, and we calculated an average root-meansquare deviation (rmsd) of 0.8 Å with the X-ray structure as a reference. At 5 Å of iron, zinc, and phosphorus atoms, we found the amino acids bound to the two metal ions and the three histidines mentioned before. The QM part of ONIOM calculations included all of these residues. Different protonation states are possible for these histidines. This choice is crucial for His202 and His296 because the Nε-atom of these histidines interacts through hydrogen bonds with the bound phosphate. Hist202 was taken protonated to be consistent with the mechanism study. Indeed, in the phosphomonoester hydrolysis mechanism, His202

QM/MM Studies of Purple Acid Phosphatase

J. Phys. Chem. A, Vol. 114, No. 26, 2010 7113

Figure 4. Cartoon representation of the crystal structure of kbPAP complex with phosphate (subunit D). The amino acids bound to the two metal ions and the three histidine residues found close to the substrate are represented with sticks.

TABLE 1: Selected Active Site Distances (Å) for All of Our Theoretical Models of the kbPAP Complexed with Phosphate metal ion

ligand atom

1DS

1DL

1Od

1Op

Zn(II)

Asp164, Oδ1 Asn201, Oδ1 His286, Nε His323, Nδ phosphate, O µ-OH, O Asp135, Oδ1 Asp164, Oδ1 Tyr167, O His325, Nε phosphate, O µ-OH, O µ-OH, O

2.23 2.16 2.23 2.19 2.21 2.06 2.01 2.10 1.90 2.23 2.15 2.02 3.11

2.21 2.15 2.21 2.18 2.27 2.07 1.96 2.14 1.91 2.22 2.14 2.01 3.06

2.21 2.16 2.19 2.29 2.32 2.00 1.90 2.30 1.92 2.20 2.35 1.94 2.70

2.14 2.16 2.18 2.20 2.56 2.08 1.88 2.28 1.91 2.26 2.16 1.86 2.76

Fe(III)

P(V)

interacts with the leaving group of the phosphomonoester (H2PO3(OR)). The only possibility to have a hydrogen bond between this substrate and the Nε-atom of His202 implies protonation of the nitrogen atom of this histidine. For His296, it is not the case. Experimentally the X-ray structure of the kbPAP complexed with phosphate was determined at pH ) 4.5. Thus, phosphate as well as the phosphate monoester in the mechanistic study are described using their monoanionic form, which is their predominant form in solution at this pH. The hydroxo group of the phosphate can interact with the Nε-atom of His296 by hydrogen bond even if this nitrogen atom is not protonated. Because of this ambiguity, we consider two ONIOM models, the model called 1Op for a protonated His296 and the model called 1Od for a nonprotonated His296. Concerning studies at the DFT level, we used for comparison purpose the two different models described before in the Methods, QM and ONIOM hybrid calculations (see 1DS and 1DL models in Figure 3). Table 1 gathers the relevant optimized distances for the different models. The relative low resolution of the X-ray structure (2.70 Å) does not allow a direct comparison between the optimized structures and the experimental data. However, comparisons between the different models are still meaningful and can give us important information. For ONIOM models, we can compare the 1Od and 1Op models. In the 1Od model, the phosphate group is bound in a nearly symmetrical way, with a Fe-O(phosphate) distance equal to 2.35 Å and a Zn-O(phosphate) distance equal to 2.32 Å. This is not the case for the 1Op model for which we found a phosphate group bound in a rather monodentate coordination, with a long Zn-O(phosphate) distance equal to 2.56 Å. Because in several phosphate-PAP structures, it has been shown

Figure 5. Superposition of some atoms in the optimized geometries of 1DL and the QM part of 1Od. Atoms of 1DL are in blue, while those of the QM part of 1Od are in magenta.

experimentally that the phosphate ion bridges to two metals in a rather symmetrical way,3,26,28,47 only the 1Od model was considered in the following. For the DFT calculations, the two models give close results. The metal-ligand distances are almost similar (see Table 1, 1DL and 1DS). Root mean square deviations (rmsd) were calculated between the two models, taking into account all of the heavy atoms of the small model. A value equal to 0.017 Å has been found. If we compare the DFT models and the ONIOM 1Od model, a longer P-O(µ-OH) distance is obtained for the DFT models. Its value is equal to 3.11 and 3.06 Å in 1DS and 1DL against 2.70 Å in the 1Od model. On the contrary, the distance between the Fe atom and the oxygen atom of the phosphate equals 2.15 and 2.14 Å in 1DM and 1DL against 2.35 Å in 1Od. Similarly, the distance between the Zn atom and the same oxygen atom of the phosphate is shorter at the DFT level: the value equals 2.21 and 2.27 Å in 1DM and 1DL against 2.32 Å in 1Od. All of these values suggest a more tightly bounded phosphate in the DFT models, but the orientation of the phosphate is modified in the ONIOM model, for which we find a smaller P-O(µOH) distance. Figure 5 illustrates these results with the superposition of some atoms of the optimized geometries of 1DL and of the QM part of 1Od. Orientation of the phosphate is linked to the position of His202 and His296, which interact through a hydrogen bond with the phosphate. These positions are noticeably different in the two models, showing the effect of the protein surrounding in ONIOM calculations. Finally, concerning the electronic structure, in DFT as well as in ONIOM calculations, the Fe(III) tends to accept all of the spin density with a value of 4.09, 4.10, and 4.11 R electrons for 1DS, 1Od, and 1DL models, respectively, in agreement with the absence of magnetic coupling between the two metal centers. Values of the Mulliken atomic charge densities are equal to 0.75 and 0.95 e for zinc and iron atoms in the 1DL model, while for the 1OD model we find 0.83 and 0.97 e. This is considerably different from their formal atomic charges, reflecting a strong donation from the ligands coordinated to the metal ions. Mechanistic Studies. Similarly to the study of the coordination of the substrate, both ONIOM and DFT calculations were performed. The phenyl phosphate was modeled in its monoanionic form (PO3H(OPh))-, coordinated to both metal atoms. To consider a nucleophilic attack of the µ-hydroxo group bridging the two metal ions, the leaving group of the substrate, the phenolate group, has to be placed in the axial position of the hydroxo group. The relevant distances for the optimized structures of the three models of reactants called, respectively, 2DSr, 2DLr, and 2Odr are presented in Table 2. After the geometry optimization, the tilted orientation of the substrate was kept, and the P-O(µ-OH) distances are in good agreement with

7114

J. Phys. Chem. A, Vol. 114, No. 26, 2010

Retegan et al.

TABLE 2: Selected Active Site Distances (Å) for the Reactant and Transition State of the Mechanism Computed with the DFT and ONIOM Models metal ion Zn(II)

Fe(III)

Zn P P

ligand atom δ1

Asp164, O Asn201, Oδ1 His286, Nε His323, Nδ phosphate, O µ-OH, O Asp135, Oδ1 Asp164, Oδ1 Tyr167, O His325, Nε phosphate, O µ-OH, O Fe µ-OH, O O(OPh)

2DSr 2DSts 2DLr 2DLts 2Odr 2Odts 2.23 2.16 2.18 2.24 2.19 2.07 2.01 2.09 1.89 2.24 2.19 2.02 3.29 3.07 1.66

2.22 2.09 2.10 2.24 2.16 2.47 1.96 2.17 1.90 2.23 2.09 2.32 3.72 1.80 2.20

2.28 2.14 2.18 2.19 2.29 2.05 1.98 2.11 1.91 2.22 2.17 2.02 3.30 3.00 1.69

2.25 2.10 2.13 2.20 2.16 2.25 1.96 2.09 1.87 2.22 2.15 2.22 3.56 1.92 2.02

2.17 2.15 2.13 2.27 2.45 2.00 1.89 2.29 1.92 2.21 2.41 1.95 3.32 2.78 1.72

2.14 2.18 2.15 2.23 2.18 2.11 1.89 2.26 1.89 2.21 2.26 2.04 3.41 2.20 2.30

those obtained for the 1DS, 1DL, and 1Od models (3.07 Å in 2DSr, 3.00 Å in 2DLr, and 2.78 Å in 2Odr). Again, we find an orientation of the substrate closer to the nucleophilic µ-hydroxide in the ONIOM model than in the small chemically intuitive DFT models. For the DFT mechanistic studies, we have fully optimized the transition state structure of the reaction. As expected in the case of the transition state structure, the most notable differences as compared to the reactant are observed around the atoms directly implicated in the reaction. The P-O(µ-OH) distance decreases from 3.07 Å in 2DSr to 1.80 Å in 2DSts and from 3.00 Å in 2DLr to 1.92 Å in 2DLts. The P-O(OPh) increases from 1.66 Å in 2DSr to 2.20 Å in 2DSts and from 1.69 Å in 2DLr to 2.02 Å in 2DLts. A direct result of this modification is that the distances of this µ-OH bridge with the metal ions and the Fe-Zn distances are increased. It is particularly clear for the small model where we note an increase of the Fe-Zn distances from 3.29 Å in 2DSr to 3.72 Å in 2DSts. For the activation energy barrier of the reaction, we obtained a high value of 156.8 kJ/mol with the 2DS model, whereas a noticeably smaller value of 91.1 kJ/mol is obtained when the large model, 2DL, is used. Because the level of theory used is the same, the activation energy barrier decrease is due to the model used, and it highly emphasizes the key role of His202 and His296. More specifically, these results can be explained by the hydrogen bond formed between the substrate and His202. In the transition state, the distance between the protonated His202 and the oxygen atom of the phenolate is highly reduced. This value reaches 1.22 Å in 2DLts against 2.03 Å in 2DLr. The transition state structure gains in stabilization thanks to this hydrogen bond. The analysis of the frequency calculations shows that the imaginary frequency at -541 cm-1, which characterizes the transition state, corresponds to a normal mode coupling the attack of the phosphorus atom by the µ-hydroxo group and the leaving of the phenolate group, but also to the hydrogen transfer from His202 to the leaving group. This later aspect highlights the role played by His202 in the reaction. Finally, the hydrogen bond formed in the 2DL reactive complex between this His202 and the substrate increases the acidity of the leaving group. In agreement with the analysis developed by Warshel et al.,48 we should observe a change in the P-O distance with respect to the nucleophile as the pkA of the system is reduced. Effectively, the P-O(µOH) distance increases from 1.80 Å in 2DSts to 1.92 Å in 2DLts. Results similar to the 2DL model should be found with the ONIOM approach. In this case, a PES scan was performed to

estimate the activation energy barrier. The chosen reaction coordinate is the difference between the two distances R1 and R2, where R1 is the P-O(phenolate) distance and R2 is the P-O(µ-OH) distance. These two distances were fixed during the optimization. The result of this scan is illustrated in Figure 6. The structure of the transition state of the reaction can be estimated by taking the highest energy structure from this scan. It is the optimized geometry for the reaction coordinate ξ equal to 0.1 that corresponds to R1 equal to 2.3 Å and R2 to 2.2 Å. During the rest of the discussion, we will call this structure 2Odts. We find an ONIOM energy difference of 21.1 kJ/mol between 2Odr and 2Odts models. This gives a good approximation of the activation energy barrier of the reaction. Thus, we can easily conclude that the mechanism is feasible and at least must be considered despite the weak nucleophilicity of the hydroxo group bridging the two metallic atoms. Comparison with DFT calculations can be done by considering only the QM energies of the high layer of our ONIOM calculations. In this case, we obtained a value of 42 kJ/mol (see Figure 6, E QM layer). This value is smaller than the value obtained for the DFT model including His202 and His296 (2DL model). The length of the hydrogen bonds formed between His202 (or His296) and the substrate is similar in the reactive complexes of 2DL and ONIOM models. The same results are obtained for the transition state structures of the two models. Consequently, these interactions cannot explain the smaller value found with the ONIOM model as compared to the DFT model. To check if the low barrier obtained with the ONIOM calculations is not an artifact, DFT energy calculations were performed on a model that corresponds exactly to the QM part of the ONIOM model. As before, the free valence was saturated with a hydrogen atom. We obtain an activation energy barrier value of 52.8 kJ/mol, in agreement with the value obtained with the ONIOM approach. Values are, of course, not strictly equal due to the inclusion of electrostatic interactions between the QM part and the protein in ONIOM calculations. It is worth noting that the electrostatic influence of the MM part on the QM wave function produces a decrease of the activation energy of 10.8 kJ/mol, so about 20% of its value, which is not negligible. The QM part of ONIOM model contains the three histidines located near the active site, His202, His295, and His296. So to compare this calculation with calculations with 2DL model, His295 was removed on the QM part and new DFT energy calculations were performed. We obtain now an activation energy barrier value of 83.4 kJ/mol, close to the value of 91.1 kJ/mol found with 2DL model. Comparisons of the two calculations show that the inclusion of the histidine 295 prompts a higher stabilization of the transition state structure as compared to the reactive, leading to the decrease of the activation energy. Figure 7 shows the optimized structures of the QM part of 2ODr and 2Odts. In the transition state, the orientation of the leaving phenolate group is in a good position to favor an electrostatic interaction between the positive charge of the protonated His295 and its aromatic cycle. This interaction is similar to the wellknown cation-π interaction49 and explains the stabilization effect of this histidine on the transition state structure. The shift of the leaving phenolate group in the structure of the transition state structure as compared to the structure of the reactive complex is due to the H transfer between His202 and the leaving group already mentioned previously. To prove this effect, new DFT energies calculations were done in removing the cationic charge on His295, that is, considering His295 not protonated.

QM/MM Studies of Purple Acid Phosphatase

J. Phys. Chem. A, Vol. 114, No. 26, 2010 7115

Figure 6. Results for the ONIOM PES scan for the relative value of the ONIOM energy. The reaction coordinate ξ is defined as ξ ) R1 - R2, where R1 is the P-O(phenolate) distance and R2 is the P-O(µ-OH) distance.

Figure 7. Superposition of the optimized structures of the QM part of 2ODr and 2Odts. Atoms of 2ODr are in blue, while those of 2Odts are in magenta.

We obtain an activation energy barrier of 86.2 kJ/mol close to the values found with models without this histidine. This result highlights the role played by this last histidine residue on the mechanism. To conclude, the ONIOM approach allows us to include correctly the three histidines (His202, His296, and His295) in our calculations. The situation is more delicate in the case of the DFT level. Indeed, the optimization, at the DFT level, of the QM part of ONIOM model starting with the optimized structures obtained in ONIOM calculations leads to such a significant displacement of the histidines that the resulting structure is an unsuitable model of the protein.50 So, if a large model has to be considered without the surrounding environment essential to maintain the histidines to a proper position, constraints must be imposed to prevent these histidines from shifting positions.

phosphate have been performed. If the metal-ligands distances around the active site are well reproduced in all calculations for the amino acids of the first coordination sphere of the iron and zinc cations, it is not the case for the orientation of the substrate. Indeed, inclusion of the three histidine residues His202, His295, and His296 with their protein surrounding as in ONIOM calculations is essential to properly describe the coordination of the substrate. They clearly influence its orientation, approaching it to the nucleophilic µ-hydroxyde bridge. Next, the mechanism of the phenyl phosphate hydrolysis by the enzyme was studied. We consider the mechanism involving the hydroxo group bridging the two metal ions as the nucleophile. DFT calculations with different sizes of models highlight the role played by His202. It forms a strong hydrogen bond with the leaving group of the substrate and stabilizes the transition state of the reaction. ONIOM calculations with geometry scans support this result and also show the stabilizing effect of His295. Thus, by gradually improving the description of the protein, we reach an estimation of the activation energy around 21.1 kJ/mol with the ONIOM model to compare to an estimation of 156 kJ/mol with the smallest model at the DFT level. Thus, this Article shows that it is essential to consider a large model, which at least includes the three histidines (His202, His295, and His296) close to the active site and not only the first coordination sphere, for any further theoretical study of the purple acid phosphatases. Acknowledgment. We thank the CECIC center and CIMENT for providing computer facilities. Supporting Information Available: Cartesian coordinates of all structures given in Tables 1 and 2. For ONIOM calculations, only the coordinates of the QM part are given. This material is available free of charge via the Internet at http:// pubs.acs.org.

Conclusion

References and Notes

In summary, comparative DFT and ONIOM calculations on the structure of (kbPAP) complexed with phosphate and phenyl

(1) Stra¨ter, N.; Klabunde, T.; Tucker, P.; Witzel, H.; Krebs, B. Science 1995, 268, 1489.

7116

J. Phys. Chem. A, Vol. 114, No. 26, 2010

(2) Schenk, G.; Gahan, L. R.; Carrington, L. E.; Mitic, N.; Valizadeh, M.; Hamilton, S. E.; de Jersey, J.; Guddat, L. W. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 273. (3) Guddat, L. W.; Mc Alpine, A. S.; Hume, D.; Hamilton, S. E.; de Jersey, J.; Martin, J. L. Structure 1999, 7, 757. (4) Lindqvist, Y.; Johansson, E.; Kaija, H.; Vihko, P.; Schneider, G. J. Mol. Biol. 1999, 291, 201. (5) Beck, J. L.; Keough, D. T.; Jersey, J. D.; Zerner, B. Biochim. Biophys. Acta 1984, 791, 357. (6) Funhoff, E. G.; Ljusberg, J.; Wang, Y.; Andersson, G.; Averill, B. A. Biochemistry 2001, 40, 11614. (7) Funhoff, E. G.; Klaassen, C. H. W.; Samyn, B.; Beeumen, J. V.; Averill, B. A. ChemBioChem 2001, 2, 355. (8) Merkx, M.; Averill, B. A. Biochemistry 1998, 37, 11223. (9) Pinkse, M. W. H.; Merkx, M.; Averill, B. A. Biochemistry 1999, 38, 9926. (10) Funhoff, E. G.; Bollen, M.; Averill, B. A. J. Inorg. Biochem. 2005, 99, 521. (11) Beck, J. L.; McArthur, M. J.; Jersey, J. D.; Zerner, B. Inorg. Chim. Acta 1988, 153, 39. (12) Beck, J. L.; Jersey, J. D.; Zerner, B.; Hendrich, M.; Debrunner, P. J. Am. Chem. Soc. 1988, 110, 3317. (13) Uppenberg, J.; Lindqvist, F.; Svensson, C.; Ek-Rylander, B.; Andersson, G. J. Mol. Biol. 1999, 290, 211. (14) Schenk, G.; Elliott, T. W.; Leung, E.; Carrington, L. E.; Mitic, N.; Gahan, L. R.; Guddat, L. W.; Gorman, P. M.; Kim, S.; Guo, M. BMC Struct. Biol. 2008, 8, 6. (15) Smoukov, S. K.; Quaroni, L.; Wang, X.; Doan, P. E.; Hoffman, B. M.; Que, L. J. Am. Chem. Soc. 2002, 124, 2595. (16) Wang, X.; Que, L. J. Biochemistry 1998, 37, 7813. (17) Mueller, E.; Crowder, M.; Averill, B.; Knowles, J. J. Am. Chem. Soc. 1993, 115, 2974. (18) Vincent, J. B.; Crowder, M. W.; Averill, B. A. J. Biol. Chem. 1991, 266, 17737. (19) Kimura, E. Curr. Opin. Chem. Biol. 2000, 4, 207. (20) Wang, X. D.; Ho, R. Y. N.; Whiting, A. K.; Que, L. J. Am. Chem. Soc. 2000, 122, 8103. (21) Dikiy, A.; Funhoff, E. G.; Averill, B. A.; Ciurli, S. J. Am. Chem. Soc. 2002, 124, 13974. (22) Bassan, A.; Blomberg, M. R. A.; Borowski, T.; Siegbahn, P. E. M. J. Inorg. Biochem. 2006, 100, 727–743. (23) Warhsel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227–249. (24) Vreven, T.; Byun, K. S.; Koma´romi, I.; Dapprich, S.; Montgomery, J. A.; Morokuma, K.; Frisch, M. J. J. Chem. Theory Comput. 2006, 2, 815. (25) Hoffman, M.; Khavrutskii, I. V.; Musaev, D. G.; Morokuma, K. Int. J. Quantum Chem. 2004, 99, 972. (26) Li, X.; Chung, L. W.; Paneth, P.; Morokuma, K. J. Am. Chem. Soc. 2009, 131, 5115. (27) Godfrey, E.; Porro, C. S.; de Visser, S. P. J. Phys. Chem. A 2008, 112, 2464. (28) Klabunde, T.; Stra¨ter, N.; Frohlich, R.; Witzel, H.; Krebs, B. J. Mol. Biol. 1996, 259, 737–748. (29) Li, H.; Robertson, A. D.; Jensen, J. H. Proteins 2005, 61, 704. (30) Spoel, D. V. D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701.

Retegan et al. (31) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (32) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (33) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269. (34) Hoops, S.; Anderson, K.; Merz, K. J. Am. Chem. Soc. 1991, 113, 8262. (35) Oda, A.; Yamaotsu, N.; Hirono, S. J. Comput. Chem. 2005, 26, 818. (36) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. J. Comput. Chem. 1997, 18, 1463. (37) Darden, T.; York, T.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089. (38) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A.; Stratmann, R. E., Jr.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G.; Ayala, A. P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Baboul, A. G.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Andres, J. L.; Gonzalez, C.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 03, revision A.1; Gaussian, Inc.: Pittsburgh, PA, 2003. (39) Bakowies, D.; Thiel, W. J. Phys. Chem. 1996, 100, 10580. (40) (a) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (b) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (41) Siegbahn, P. E. M. J. Biol. Inorg. Chem. 2006, 11, 695. (42) Dunning, T. H., Jr.; Hay, P. J. In Method of Electronic Structure Theory; Schaefer, H. F., Ed.; Plenum Press: New York, 1977; pp 1-28. (43) Hay, P. J.; Wadt, W. R. J. Chem. Phys. 1985, 82, 299. (44) First calculations on DFT models were done using the standart lanl2dz basis for all atoms. With the addition of polarization functions on oxygen, nitrogen, and phosphorus atoms, which should give improved results, the geometry of the models changes dramatically. The distance between the phosphorus atom and the oxygen atom of the µ-hydroxo bridge rises considerably to a value of 3.11 Å against 2.29 Å without polarization functions for the small model 1DS. Thus, for all other calculations, a set of polarization functions was added. (45) Davis, J. C.; Averill, B. A. Proc. Natl. Acad. Sci. U.S.A. 1982, 79, 4623. (46) Torrent, M.; Veven, T. G.; Musaev, D.; Morokuma, K. J. Am. Chem. Soc. 2002, 124, 193. (47) True, A. E.; Scarrow, R. C.; Randall, C. R.; Holz, R. C.; Que, L., Jr. J. Am. Chem. Soc. 1993, 115, 4246. (48) Rosta, E.; Kamerlin, S. C. L.; Warshel, A. Biochemistry 2008, 47, 3725. (49) Dougherty, D. A. Science 1996, 271, 163. (50) One test gives an important displacement of His296. The Hbond between this same histidine and the substrate does not exist anymore and is replaced by π-stacking interactions between His295 and His296.

JP100478F