Quantum and Molecular Mechanical Study of the First Proton Transfer

5126

J. Phys. Chem. B 2008, 112, 5126-5138

Quantum and Molecular Mechanical Study of the First Proton Transfer in the Catalytic Cycle of Cytochrome P450cam and Its Mutant D251N Dongqi Wang,† Jingjing Zheng,† Sason Shaik,‡ and Walter Thiel*,† Max-Planck-Institut fu¨r Kohlenforschung, Kaiser-Wilhelm-Platz 1, D-45470 Mu¨lheim an der Ruhr, Germany, and Department of Organic Chemistry and the Lise Meitner Center for Computational Quantum Chemistry, The Hebrew UniVersity, 91904 Jerusalem, Israel ReceiVed: June 26, 2007; In Final Form: January 30, 2008

In the catalytic cycle of cytochrome P450cam, the hydroperoxo intermediate (Cpd 0) is formed by proton transfer from a reduced oxyheme complex (S5). This process is drastically slowed down when Asp251 is mutated to Asn (D251N). We report quantum mechanical/molecular mechanical (QM/MM) calculations that address this proton delivery in the doublet state through a hydrogen-bond network in the Asp251 channel, both for the wild-type enzyme and the D251N mutant, using four different active-site models. For the wildtype, we find a facile concerted mechanism for proton transfer from protonated Asp251 via Wat901 and Thr252 to the FeOO moiety, with a barrier of about 1 kcal/mol and a high exothermicity of more than 20 kcal/mol. In the D251N mutant with a neutral Asn251 residue, the proton transfer is almost thermoneutral or slightly exothermic in the three models considered. It is still very facile when the Asn251 residue adopts a conformation analogous to Asp251 in the wild-type enzyme, but the barrier increases significantly when the Asn251 side chain flips (as indicated by classical molecular dynamics simulations). This flip disrupts the hydrogen-bond network and hence the proton-transfer pathway, which causes a longer lifetime of S5 in the D251N mutant (consistent with experimental observations). The entry of an additional water molecule into the active site of D251N with flipped Asn251 regenerates the hydrogen-bond network and provides a viable mechanism for proton delivery in the mutant, with a moderate barrier of about 7 kcal/mol.

I. Introduction Proton delivery is of crucial importance in many condensedphase processes. Pioneering work was reported as early as 1806 by de Grotthuss1 who envisioned a cooperative mechanism for passing a proton through a chain of water molecules. The microscopic details of proton transfer in water and in aqueous solution have recently been studied both by computer simulations2 and by femtosecond spectroscopy.3,4 Hydrogen-bond networks composed of water molecules are essential for proton transfer not only in aqueous solution but also in enzymes.5 The protein-mediated long-range proton transfer normally involves two steps: translocation of a proton from bulk solvent into the active site and reorientation of the hydrogen-bond network to facilitate the migration of the proton. It is thus necessary to identify and characterize possible aqueducts in the enzyme. Titratable amino acid residues along the conducting pathways may play a special role because they can adopt different protonation states and because the protein environment can tune their proton affinities. In general, several such pathways may coexist. Evidence on their location can be obtained from X-ray crystal structures which, however, offer only a static picture and a limited resolution (hydrogen atom positions normally not available). A detailed understanding of proton-transfer pathways in enzymes thus calls for the use of other experimental techniques such as kinetic and mutagenesis studies as well as computer simulations. * To whom correspondence should be addressed. E-mail: thiel@ mpi-muelheim.mpg.de. † Max-Planck-Institut fu ¨ r Kohlenforschung. ‡ The Hebrew University.

Cytochrome P450 enzymes constitute a superfamily of monooxygenases with a heme-containing active site.6 They are crucial for a variety of vital processes including biosynthesis, biodegradation of xenobiotics, and drug metabolism. The best characterized member of this family, P450cam,7 catalyzes the hydroxylation of nonactivated camphor to 5-exo-hydroxycamphor in a regio- and stereospecific manner at physiological temperatures. The consensus mechanism for the catalytic cycle of P450cam involves two proton transfers: the ferrous dioxygen complex (S4) is first reduced to the ferric peroxo complex (S5) and then protonated to yield the hydroperoxo compound (Cpd 0) which undergoes a proton-assisted heterolytic cleavage of the O-O bond to produce the putative active oxoferryl species (Cpd I). The formation of Cpd 0 and Cpd I thus requires two subsequent proton transfers. The available crystal structures7 indicate three possible proton delivery pathways through the Asp251,5,8,9 Glu366,7,10 and Arg29911 channels. The latter offers a connection from the active site via the propionate-A side chain of heme and Arg299 to bulk solvent water in the case of the resting state but is blocked after the entry of the substrate and will therefore not to be considered further (see Figure 1). The Glu366 channel has received much attention due to the presence of a hydrogen-bond network between the carboxyl group of Glu366 and the distal oxygen atom (Od) of the O2 moiety (1DZ87) that involves bridging water molecules (Wat523, Wat566, Wat687, Wat902) and the hydroxyl group of Thr252. This network remained stable during molecular dynamics (MD) simulations (on P450eryF)12 and during an annealing procedure with a frozen heme-Cys357 unit.10 According to QM/MM calculations,10 the proton transfer from Thr252 to Od has to

10.1021/jp074958t CCC: $40.75 © 2008 American Chemical Society Published on Web 04/03/2008

Proton Transfer in P450cam

J. Phys. Chem. B, Vol. 112, No. 16, 2008 5127

Figure 1. Active-site environment in cytochrome P450cam from X-ray structure 1DZ8.7 Red dots indicate crystal water molecules.

overcome a barrier of 4.0 kcal/mol and is exothermic by 2.8 kcal/mol, and the subsequent proton transfer from nearby Wat902 to Thr252 is facile with a barrier of 1.0 kcal/mol (remaining steps for proton delivery from Glu366 not studied in ref 10). However, there is also some compelling evidence against the relevance of the Glu366 channel for proton delivery in P450cam. Most notable are mutation experiments13,14 where much of the enzymatic activity is conserved when Glu366 is replaced by methionine (70% for O2 consumption and 85% for coupling, compared with the wild-type). Moreover, Glu366 does not have a direct connection with the bulk solvent, and there is no obvious mechanism for reprotonating Glu366 after proton delivery (estimated pKa value of 7.55 by PropKa15). The Asp251 channel does not suffer from this problem. Asp251 is connected to the protein surface, and hence the bulk solvent, through a salt bridge with Arg186 that is linked to a catalytic triad composed of Arg186-Lys178-Asp18216 or a tetrad of Lys178-Asp182-Arg186-Asp251.17 Having access to bulk water, deprotonated Asp251 can in principle accept a proton from bulk water. Asp251 controls the entry of water molecules into the active site through a “carboxylate switch” that blocks the entry when Asp251 is deprotonated.9 This role of Asp251 is supported by mutagenesis experiments. Substitution of Asp251 (D251) by Asn (N) disrupts the salt link with Arg186, as well as the hydrogen bond between Asp251 and Thr185, which governs the orientation of the carboxyl group of Asp251. In the mutant D251N, Asn251 rotates its amide group out to the bulk solvent (crystal structure 6CP49,18) and forms a hydrogen bond with the carbonyl group of Asp182, which will make the active site more accessible for solvent water molecules and lead to unspecific proton delivery into the protein pocket. The observed9 significant increase of the solvent kinetic isotope effect (SKIE) in D251N supports this notion and suggests that a poor proton source (e.g., water) may be involved in the proton delivery in the mutant. The presence of extra solvent molecules around the Fe-O2 active site leads to unspecific product (uncoupling to peroxide or heterolysis to water). The major concern about the Asp251 channel comes from the implied rupture of the Asp251-Arg186 salt bridge. It is wellknown that salt bridges are important for stabilizing proteins in

the condensed phase. According to a numerical PoissonBoltzmann study,17 there are two strong salt links involving Asp251 in P450cam (∆G: Lys178-Asp251, -7.3 kcal/mol; Arg186-Asp251, -6.0 kcal/mol), while the corresponding interactions are much weaker in the D251N mutant (∆G > -2.4 kcal/mol). Multiple Asp251 side chains sampling using an allatom force field (OPLS-AA) combined with a continuum solvation model10 did not show any significant motion of Asp251 that would be large enough to rupture the salt bridge. By contrast, recent studies on network models5 and MD simulations19 revealed that Asp251 could move its side chain sufficiently to assist in the proton delivery. Previous theoretical studies have mainly addressed the second proton transfer from Cpd 0 to Cpd I using models of increasing complexity.20-24 The latest such density functional theory (DFT) investigation25 employed an active-site model with 96 atoms (including 5 water molecules and key residues from the protein environment). Treating the real system rather than model systems requires the use of combined quantum mechanical/ molecular mechanical (QM/MM) methods, which describe the active site at the QM level and the environment at the MM level to provide realistic results at affordable costs. Two QM/ MM studies have been reported for the second proton transfer in P450cam: the first one covered selected steps in the Glu366 channel using Thr252 as the immediate proton source10 (see above), while the second one investigated both the Glu366 and Asp251 channels and proposed a novel mixed homolyticheterolytic mechanism in the Asp251 channel as the most favorable mechanism.19 Much less theoretical work has been done on the first proton transfer in P450cam that leads from S5 to Cpd 0 (Scheme 1), in spite of the fact that Cpd 0 is the last observable species in the native enzyme and has been invoked as an alternative oxidant either in its own right or in cases where the formation of Cpd I is difficult.21,26-29 Here we report a QM/MM study of the first proton delivery that produces Cpd 0 in the wild-type P450cam enzyme and its mutant D251N. The paper is organized as follows. Section II summarizes the computational methods and protocols employed. Section III presents classical MD and QM/MM results for different proton-transfer pathways. Section

5128 J. Phys. Chem. B, Vol. 112, No. 16, 2008 SCHEME 1: Cytochrome P450cam Catalytic Cycle

IV discusses our results and relates them to experiment. Section V offers conclusions. Throughout the paper, we use the notation shown in Scheme 1. II. Computational Methods Initial coordinates were taken from an available crystal structure7 (PDB code 1DZ8, resolution 1.9 Å). The solvation and protonation procedures were the same as in our previous work.30 The titratable residue Asp297 was protonated so that it can form a hydrogen bond with the nearby A-propionate side chain of the heme. All histidine residues were protonated on Nδ, or N, or on both N atoms, based on visual inspection of their chemical environment. Asp251 was also protonated to serve as proton source in our present study. Glu366 was tested as alternative proton source, but the resulting pathways were not competitive with those in the Asp251 channel and will therefore not be discussed. The solvated system with ∼25 000 atoms was relaxed by energy minimizations and classical MD simulations using the CHARMM22 force field31 as implemented in the CHARMM32 program. The oxyheme unit (O2-418 and heme), Cys357, and the outer 8-Å solvent layer were fixed during these MM calculations. A randomly selected snapshot, hereafter referred to as Model A, was used as the starting structure for the following QM/ MM studies on proton transfer. To investigate the D251N mutant, we manually replaced residue Asp251 by Asn in model

Wang et al. A, and subjected it first to 3600 steps of ABNR (adopted-basis Newton-Raphson) minimization, and then to QM/MM optimization. This resulting system is denoted as Model B. During a classical 1-ns MD simulation of model B, the side chain of Asn251 was found to flip outside. A randomly chosen snapshot with such a flipped side chain was optimized first at the MM and then at the QM/MM level to prepare Model C. After the flip of the Asn251 amide side chain, there is some empty space between Arg186 and Wat901, which is bordered by both hydrophilic (Asn251:OD1) and hydrophobic groups (Thr185:CG2). We manually added another water molecule (WatS) in this region to bridge Wat901 and Arg186, and performed an ABNR minimization to generate Model D1. In the following classical 1-ns MD simulation, WatS remained in a stable topology, forming three hydrogen bonds with Arg186, Asn251, and Wat901. A randomly selected snapshot from this trajectory was subjected to MM minimization to produce Model D2. The solvent/protein environment is rather similar in models C and D1 and is more adapted to the extra water molecule (WatS) in model D2. The applied QM/MM methodology is analogous to that in our previous work,19 and therefore we only briefly describe some aspects relevant to the present work. An electronic embedding scheme33 was adopted in the QM/MM calculations, i.e., the MM charges were included in the one-electron Hamiltonian of the QM part, and the QM/MM electrostatic interactions were evaluated from the QM electrostatic potential and MM partial charges. No cutoffs were introduced for the nonbonding MM and QM/MM interactions. Hydrogen link atoms and the charge shift model34 were employed to treat the QM/MM boundary. The TURBOMOLE35 program was used for the QM calculations while the CHARMM22 force field was run through the DL_POLY36 code for the treatment of the MM part. The QM/ MM calculations were performed with the ChemShell package37 that integrates the TURBOMOLE and DL_POLY programs, and geometry optimizations were carried out using the HDLC optimizer38 implemented in ChemShell. All minima and transition states (TS) reported in this paper have been fully optimized. The QM part was described by unrestricted DFT. We used the (U)B3LYP39 hybrid functional (TURBOMOLE default version, correlation being specified as 0.19 VWN(V) + 0.81 LYP), in combination with three basis sets BS1, BS2, and BS3 that are defined as follows: BS1:LACVP40 small-core ECP basis set (double-ξ quality) for Fe and 6-31G41 for the other atoms (used in pathway calculations, 391-406 basis functions for models A-D); BS2:DZpdf basis42,43 for Fe and 6-31G* for the other atoms (used in single-point calculations, 600-620 basis functions for models A-D). BS3: def2-TZVP basis44 (TURBOMOLE notation) for Fe and 6-311++G(2d,2p)45 for the other atoms (used in single point calculations, 1248-1305 basis functions for models A-D). QM Region. In our previous studies, SH was found to be a good mimic of the Cys357 residue coordinated to Fe at the proximal position.30 This was confirmed in our recent QM/MM work on the oxyheme complexes S4 and S5. Consequently, we adopted in our present QM/MM calculations the following QM region (61/62 atoms for WT/D251N in models A-C, three more atoms in model D, see Figure 2): Iron porphine (without heme side chains), sulfur atom of Cys357, axial O2 moiety, the Cγ2H3sCβHsOγ1H unit of Thr252, Wat901, WatS (for model D), the CβH2sCγ(dOδ1)(sOδ2H) unit of Asp251 (or the CβH2sCγ(dOδ1)(sNδ2H2) unit of Asn251 in the case of the mutant).



Figure 2. QM region of models A-D (note different orientations of Asn251 and added WatS).

Figure 3. Motion of the Asn251 side chain during the MD simulation. Left, changes of torsion angles OD1-CG-CB-CA and ND2-CG-CB-CA; right, definition of atom labels.

Active Region. We specify a core region consisting of the non-hydrogen atoms of heme, Cys357, Glu366, Asp251, Thr252, Wat523, Wat566, Wat687, Wat901, and Wat902. The optimized active region was defined to include all residues and water molecules within 6 Å of any atom of the core region. This results in ca. 1400 atoms to be optimized, which belong to the irondioxygen-porphyrin complex, the amino acid residues in the active site, camphor, and (crystallographic and solvent) water. Residue notations and numbering schemes in this paper are the same as in the crystal structure (1DZ8). Further computational details are documented in the Supporting Information. III. Results 1. Asn251 Side-Chain Flip. Models B and C of the mutant D251N differ in the orientation of the amide group of Asn251. They exhibit the side-chain orientation found in the X-ray structures of the wild-type enzyme (1DZ87) and the mutant (6CP49), respectively. A more recent X-ray study46 reports a refined X-ray structure of the mutant (2A1N) that contains two different Asn251 conformations in chains A and B (flipped or not). We have chosen two torsion angles between the amide side chain and the backbone of Asn251 (OD1-CG-CB-CA and ND2-CG-CB-CA) as measures of the orientation of the amide side chain and monitored the change of these two torsion angles

during a 1-ns classical MD simulation starting from a model B structure (Figure 3). It is obvious from Figure 3 that the side chain flips after 45 ps, since both torsion angles change abruptly from 30°/-150° to -100°/100° (for OD1-CG-CB-CA/ ND2-CG-CB-CA). In the new conformation, the flipped amide side chain of Asn251 forms strong hydrogen bonds with nearby Thr181:OG1 (1.949 Å from Asn251:HD22), Asp182:OD1 (1.576 Å from Asn251:HD21), Arg186:HH22 (1.707 Å from Asn251:OD1), and Arg186:HH12 (1.878 Å from Asn251:OD1). In consideration of the crucial role of Wat901 in the proton delivery channel, we have checked its mobility during the MD simulation. The flip of Asn251 generates some empty space in the distal pocket so that Wat901 could move around (between Thr252 and Arg186). By contrast, the Thr252:OH group is constrained by the protein environment to remain essentially in its position. We have therefore used the distance between Wat901:OH2 and Thr252:OG1 (see Figure 3) to monitor the motion of Wat901 during the MD simulation (Figure 4). Evidently, Wat901 does not escape from the protein pocket during the 1-ns simulation but stays close to Thr252 for most of the time. The distance between Thr252:OG1 and Wat901: OH2 remains around 3 Å, indicating that the hydrogen bond between Thr252 and Wat901 is conserved. During the limited simulation time, no extra water molecule enters the distal pocket

5130 J. Phys. Chem. B, Vol. 112, No. 16, 2008

Figure 4. Monitoring the mobility of Wat901 (atom labels see Figure 3).

Figure 5. Superimposition of QM/MM optimized structures (QM region). Left, models A (green: wild-type), B (blue: D251N, no flip), and C (red: D251N, flipped); right, models C (red: without WatS) and D1 (cyan: with WatS).

even though there seems to be enough space in this region after the Asn251 flip. This may be due to the partially hydrophobic environment and/or an insufficient sampling time. To explore its influence, we have manually inserted an additional water molecule (WatS) into this region, which can form hydrogen bonds with Wat901, Asn251:OD1, and Arg186 (model D1, see section II). Figure 5 presents overlays of the QM/MM optimized structures for the QM regions of models A-D1. Compared with the Asp251 side chain of the wild-type enzyme (model A), there is a notable rotation of the Asn251 amide side chain of the D251N mutant into the active site for model B, possibly due to the repulsion from the nearby Arg186 residue. The flip of the Asn251 side chain in model C relieves this repulsion so that the side chain moves out of the pocket again. The overlay of the structures for models C and D1 indicates that the addition of the extra water molecule and the subsequent geometry optimization only cause minor structural changes (slight shift of Asn251). 2. Proton-Transfer Pathways. In this section we present the QM/MM results for the proton-transfer pathways obtained for models A-D. In each case, the proton is delivered to the O2 moiety through the hydrogen-bond network in the Asp251/ Asn251 channel. We discuss the QM/MM-optimized structures of the reactant complex (RC, S5), transition state (TS), and product complex (PC, Cpd 0). The corresponding energy profiles

Wang et al.

Figure 6. QM/MM energy energy profiles for proton transfer in model A (wild-type) and models B-D (D251N mutant). Note that model A contains protonated Asp251, while models B-D contain neutral Asn251. The energy profiles are thus directly comparable only among models B-D (D251N mutant).

are collected in Figure 6. Unless noted otherwise all results refer to the doublet state and basis set BS1. A. Wild-Type P450cam: Model A. The QM/MM calculations predict a facile concerted proton transfer with a barrier of only 1.2 kcal/mol and an exothermicity of 22.0 kcal/mol (Figure 6). Such an energy profile is consistent with the experimental finding that the reactant oxyheme complex S5 cannot be detected at or above 77 K in wild-type P450cam.22,47 Upon cooling to liquid He temperature, an electron paramagnetic resonance signal was recorded that has been ascribed to this elusive species.22,48 The optimized geometries of the stationary points (Figure 7) show a well-organized hydrogen-bond network that extends from Asp251 (which is connected to the surface of the protein) via Wat901 and Thr252 to Od, which needs to be protonated. Along this network, the proton from the Asp251 carboxyl group migrates to the bridging Wat901 molecule that donates one of its other protons to Thr252, which in turn delivers the proton from its hydroxyl group to the distal Od atom of the O2 moiety to form the hydroperoxo species. The TS geometry is much closer to the reactant S5 than to the product Cpd 0 (see, for example, the bond lengths for OdOp and Fe-Op in Figure 7, which are almost the same for RC and TS). We thus have an early transition state, in accordance with the Hammond postulate49 for strongly exothermic reactions. In the product, the Fe-Op and Fe-S bonds are shorter (stronger) than in the reactant. Thr252 now acts as a hydrogenbond acceptor that helps to stabilize the newly formed hydroperoxo species, as suggested previously.46 B. Mutant D251N: Model B (No Asn251 Flip). In model B, the hydrogen-bond network between the proton source (Asn251 in D251N) and the O2 unit is conserved (Figure 8). The barrier for the proton transfer is calculated to be 2.2 kcal/mol (i.e., 1.0 kcal/mol higher than in the wild-type enzyme), and the reaction is only mildly exothermic, by 4.8 kcal/mol (Figure 6). Inspection of the optimized geometries indicates that the transition state occurs somewhat later in the mutant than in the wild-type enzyme (cf. Figures 7 and 8). Going from reactant to transition state and product in model B, the Od-Op bond lengthens from 1.395 to 1.444 and 1.511 Å (indicating loss of superoxo character), while the Fe-Op and Fe-S bonds shorten from 2.055 to 1.903 and 1.831 Å and from 2.601 to 2.585 and 2.505 Å, respectively. These geometrical changes are ac-



Figure 7. QM/MM-optimized geometries of reactant (RC), transition state (TS), and product complex (PC) for proton transfer in model A (wildtype enzyme).

Figure 8. QM/MM-optimized geometries of reactant (RC), transition state (TS), and product complex (PC) for proton transfer in model B (mutant D251N, no flip).

companied by an electron flow from Cys357 via Fe to the O2 unit (see analysis in section III.2.F). Given the high basicity of amides, it is surprising at first sight that the proton transfer from the Asn251 amide side chain to the heme is computed to be a rather facile process. Previous work reports pKa values of 3.6-4.0 for HOO-Fe in HRP and CcP,50,51 while the moieties in the distal pocket, Thr252, Wat901, and Asn251, should be intrinsically more basic (pKa values:52 ethanol 16, water 15.74, acetamide 17) and thus reluctant to donate a proton. In the enzyme, however, there is a positively charged polar residue, Arg186, close to Asn251 which may be expected to exert a repulsive force on a leaving proton and hence to facilitate the proton transfer from Asn251 to the heme through this electrostatic repulsion. Figure 9 shows a more detailed view of this region in model B. When going from the reactant to the product, the contact between Asn251 and Arg186 obviously becomes tighter (smaller distances, stronger interaction), concomitant with a lengthening of the amide CO bond by 0.022 Å and a shortening of the amide CN bond by 0.014 Å (implying a more pronounced delocalization in the amide). The computed Mulliken charges become more negative by 0.08e at the amide oxygen atom and more positive by 0.04e at the amide nitrogen atom. These data indicate that the positively charged Arg186 residue is capable of stabilizing an enolate-type electron distribution in Asn251, which

Figure 9. QM/MM-optimized geometries of reactant (RC) and product complex (PC) in model B: detailed view of the region around Asn251 and Arg186.

will make it easier to remove a proton from Asn251. This influence of Arg186 can be represented qualitatively by the resonance structures for the reactant shown in Scheme 2. C. Model D251N: Model C (Flipped). In the flipped conformer, the amide side chain of Asn251 no longer points toward the heme. Having lost the connection with the hydrogen-

5132 J. Phys. Chem. B, Vol. 112, No. 16, 2008 SCHEME 2: Possible Resonance Structures for the Reactant in Model B: Asn251-Arg186 Pair

bond network in the distal pocket, Asn251 cannot serve as proton source any more, and consequently the influence of Arg186 will be diminished as well (see above). In model C, the proton transfer occurs between Thr252 and the FeOO moiety (see Figure 10). It has to overcome a barrier of 6.9 kcal/mol and is endothermic by 5.4 kcal/mol. The reactant (S5) is thus thermodynamically favored over the product (Cpd 0), which can easily rearrange back to S5. The optimized geometries indicate a late transition state (Figure 10). Compared with S5, the Od-Op bond is elongated by 0.08 Å, while the Fe-Op bond is shortened by 0.23 Å, and both are already of similar length as in Cpd 0. The Fe-S bond distance assumes an intermediate value in the transition state. In model C, the formed Cpd 0 contains a deprotonated Thr252 residue. This anion is stabilized by two strong hydrogen bonds with Wat901 and FeOOH (Figure 10) but remains a high-energy species so that the proton transfer is endothermic (see above). We have tried to neutralize the Thr252 anion by abstracting one proton from Wat901, but this proton always moved back to the emerging hydroxide anion upon full QM/MM optimization, indicating that Wat901 is an even weaker Bro¨nsted acid than Thr252 in the given environment. In these attempts, we have also included the Arg186 side chain into the QM region, in the hope of using Arg186 as proton source, and by judicious positioning of the Wat901 hydroxide anion, we have indeed found a proton-transfer path that results in a species with neutral Thr252, Wat901, and Arg186 (see Supporting Information). However, this species lies 11.7 kcal/mol above the reactant (S5) and is thus less stable than the initially formed Cpd 0 with deprotonated Thr252. Hence, in model C with flipped Asn251, the best way to stabilize this Cpd 0 species would seem to be the rearrangement back to the reduced oxyheme complex S5 (via a barrier of only 1.5 kcal/mol). Given this situation, it is conceivable that the proton transfer from the model C conformer follows a different path. The amide side chain of Asn251 could flip back to generate the less favored model B conformer, which could then undergo a facile and exothermic proton transfer (see above, Figures 6 and 8) to form Cpd 0 with deprotonated Asn251. To explore this possibility, we have studied the conversion between conformers B and C by classical MD simulations using the CHARMM force field and adaptive umbrella sampling (see Supporting Information for details). With the chosen sampling parameters, the statistical uncertainty of the resulting free energies is of the order of 1 kcal/mol. The computed free energy profile as a function of the torsion angle OD1-CG-CB-CA (Figure 3) is plotted in Figure 11. We find two minima with average torsion angles of 45° and -95° which correspond to models B and C (QM/MM optimized torsion angles of 71° and -78°, respectively). Conformer B lies in a shallow well (7.3 kcal/mol above conformer C) and can easily rearrange to conformer C (free energy barrier of 0.7 kcal/mol) as observed in our initial

Wang et al. unconstrained MD simulations. In the reverse direction, the back flip of the Asn251 side chain requires an activation of 8.0 kcal/ mol (statistical uncertainty of 0.7 kcal/mol). If this back flip is followed by a proton transfer in the resulting model B conformer (QM/MM barrier of 2.2 kcal/mol, see above and Figure 6), the overall effective barrier for proton transfer in the model C conformer will amount to ca. 9.5 kcal/mol. In the absence of an easy way to reprotonate Thr252 (see above), this stepwise mechanism may actually be most favorable for model C. In this case, Asn251 acts as the proton source (which can be replenished from bulk solvent); the overall effective barrier of 9.5 kcal/mol is not excessive, and the proton transfer is slightly endothermic by 2.5 kcal/mol. It should be noted that the most recent X-ray structure (2A1N) of the D251N mutant contains conformer C in chain A and conformer B in chain B, which suggests similar stability of these conformers; however, 2A1N differs from our model structure by the absence of Wat901 and the conformation of helix I, and is thus not directly comparable. D. Mutant D251N: Model D (Flipped with WatS Added). Model D differs from model C by the presence of an extra water molecule (WatS) in the space being vacated by the flip of the Asn251 residue (see section II). The introduction of WatS extends the hydrogen-bond network of model C toward the Arg186 residue, which is expected to facilitate the proton transfer. Model D1 has been constructed from model C with minimum adjustment of the protein/solvent environment (no MD runs, only optimization). QM/MM calculations on the proton-transfer pathways yield Wat901 as the proton source. In a concerted process, Thr252 donates its hydroxyl proton to the distal Od atom and abstracts another proton from Wat901 (Figure 12). The formed hydroxide anion (Wat901), with a Mulliken charge of -0.56e, is stabilized by two hydrogen-bond donors on both sides (Thr252 and WatS). The barrier to the proton transfer is calculated to be 4.6 kcal/mol (i.e., ca. 2 kcal/mol lower than the value obtained for model C), and the whole process is almost thermoneutral with the product being slightly favored (by ca. 1 kcal/mol). The presence of the extra water molecule in the protein pocket thus seems to promote the proton transfer both kinetically and thermodynamically. Model D2 improves on model D1 by an initial MD equilibration that allows the system to relax its protein/solvent environment (see section II). This entails some motion of Arg186 to enlarge the space available to WatS, which forms three hydrogen bonds with Wat901, Asn251, and Arg186. QM/MM pathway calculations with this equilibrated model D2 are sensitive to the definition of the reaction coordinate. When using the same coordinate as in model D1, i.e., the difference d(OThr252-HThr252) - d(Od-HThr252) between the distances in the breaking and forming OH bonds, we observe facile dissociation of the FeOp bond with a barrier of ca. 1 kcal/mol (estimated from the energy profile) and formation of an intermediate, with a superoxide anion situated between Thr252 and Fe (see Figure S7 in Supporting Information, Od-HThr252 ) 1.469 Å, Fe-Op ) 3.376 Å). When choosing a reaction coordinate that disfavors Fe-Op dissociation, i.e., d(OThr252-HThr252) - d(Od-HThr252) - d(Fe-Op), we find proton transfer (see Figure S6). The fully optimized structure of the transition state (Figure 13) is similar to that found in model D1 (Figure 12). The same applies to the product complex where the formed hydroxide anion (Wat901) is again stabilized by two strong hydrogen bonds to Thr252 and WatS. The computed barrier (6.7 kcal/mol) and exothermicity (-1.0 kcal/mol) for proton transfer are also similar to the values obtained in model D1 (see Figure 6). The two models



Figure 10. QM/MM-optimized geometries of reactant (RC), transition state (TS), and product complex (PC) for proton transfer in model C (mutant D251N, flipped).

Figure 11. Free-energy profile for the conversion of conformers B and C from classical MD simulations using CHARMM and umbrella sampling.

thus yield completely analogous results for proton transfer, but they differ in that the relaxation of the protein/solvent environment in model D2 facilitates a competing reaction, i.e., uncoupling with Fe-Op dissociation. This is already indicated in the optimized structure of the reactant where the Fe-Op bond is significantly longer in model D2 (2.183 Å) than in model D1 (2.057 Å) or Model C (2.096 Å). The proton transfers in models D1 and D2 generate a hydroxide anion (Wat901) that needs to be reprotonated from the bulk solvent. This may be mediated by the Arg186 residue whose pKa value is estimated by PropKa to be around 8 for RC and PC in model D2. To test this hypothesis, we have extended model D2 by including the Arg186 side chain into the QM region and studied the entire proton-transfer process between the heme and Arg186 in the resulting model D2′ (Figure 14). Not surprisingly, there is very little change in the geometries of RC, TS, and PC (see Figures 13 and 14), and the energy profile remains essentially the same in model D2′ (barrier of 6.6 kcal/mol, exothermicity of -2.3 kcal/mol). The subsequent reprotonation of the hydroxide anion by Arg186 is extremely facile (barrier of 0.8 kcal/mol) and almost thermoneutral (overall exothermicity of -2.8 kcal/mol). Since Arg186 is in contact with bulk solvent, model D2′ thus offers a viable proton-transfer path in the D251N mutant.

E. Basis Set Effects. Table 1 summarizes the relative energies obtained from geometry optimizations using basis BS1 and compares them with the values obtained from single point calculations with the larger basis sets BS2 and BS3 (see section II). It is obvious that the computed barriers do not change much upon basis set extension: the three values are within 2 kcal/ mol in each case and normally even within 1 kcal/mol. The calculated reaction energies change significantly in two cases (i.e., for models A and B when going from BS1 to BS2) but much less in all other cases, and the BS2 and BS3 values are generally consistent with each other: upon basis set extension the proton transfer becomes more exothermic for models A, B, and D but remains endothermic for model C. Our qualitative conclusions from this section are thus supported by the singlepoint QM/MM results with the larger basis sets. F. Spin Densities and Partial Charges. The changes in the electronic structure during proton transfer can be monitored through the computed spin densities and Mulliken charges on key atoms and groups. We consider Fe, Op, Od, SH, and porphyrin for this purpose. Since the observed trends are similar in all four models, we only present plots of the results for model A (wild-type) and show the others in Supporting Information. Table 2 collects the numerical results from models A-D for the spin densities, while those for the charges are again given in Supporting Information. It is obvious from Figure 15 that the major changes in the spin densities occur in the FeOO moiety. Initially, in S5, the unpaired electron resides in the antibonding π* orbital of the O2 unit. In the TS region, one electron leaves an iron d orbital to fill the orbitals of the O2 unit, which reorganize to accommodate the new bond to the incoming proton. As a result, Fe is oxidized during the reaction, and the O2 unit changes from superoxide to peroxo character. Inspection of Table 2 shows that the spin density on Fe in the product decreases when going from model A to model C (from 0.930 to 0.844), while the spin density on Op increases (from 0.137 to 0.201). This indicates a somewhat more pronounced spin delocalization in the mutant (model C) and hence less closed-shell character of the hydroperoxo moiety, which seems consistent with its lower stability. According to Figure 16, the proton-transfer causes an increase of the positive charge on Fe (by 0.06e) and of the negative charge on the O2 unit (by 0.10e), but the changes are smaller


Wang et al.

Figure 12. QM/MM-optimized geometries of reactant (RC), transition state (TS), and product complex (PC) for proton transfer in model D1 (mutant D251N).

Figure 13. QM/MM-optimized geometries of reactant (RC), transition state (TS) and product complex (PC) for proton transfer in model D2.

than might have been anticipated from the qualitative orbital arguments given above. This is due to significant donation of electron density from the other ligands (i.e., SH and porphyrin) to Fe, which provides substantial stabilization. The computed total charges on Thr252 merit some comment (see Supporting Information for the numerical data). In model A (wild-type), they remain essentially unchanged during the proton transfer, while they become much more negative in model C (D251N mutant) because of the formation of the Thr252 anion. However, even in the product of model C, the total charge on Thr252 does not exceed -0.68e which reflects significant charge delocalization due to strong interactions with the surroundings (see section III.2.C). In summary, the spin density and charge analysis supports the simple qualitative notion that Fe is oxidized during the proton transfer, while the O2 moiety is reduced (going from superoxo to peroxo character), with the proviso that electron donation from the other ligands moderates these changes. In this sense the proton transfer can be described as a one-electron process (Scheme 3).

IV. Discussion Experimentally, the reduced oxyheme complex (S5) has been detected in the wild-type enzyme only at liquid He temperatures of 4-7 K.22,48 This is consistent with the very low barriers for proton transfer that are calculated for model A (Table 1). Similarly low barriers are obtained for the D251N mutant in model B where the Asn251 residue has a conformation analogous to that of the Asp251 residue in wild-type P450cam. However, the available X-ray structures and our present MD simulations indicate that the amide side chain of Asn251 flips outside. In the resulting more stable conformer (model C), the proton transfer is computed to require significantly more activation than in the case of the wild-type enzyme, which remains true when an extra water molecule is added to the active site (model D). This is again in agreement with the available experimental evidence: In the D251N mutant, S5 can be trapped under conditions53 where this is not possible in the wild-type enzyme, and it can be characterized spectroscopically at a temperature of 77 K.22,48 The comparison between the QM/ MM results for models A-D thus confirms the notion9 that this



Figure 14. QM/MM-optimized geometries of stationary points (RC, TS, PC, TS′, and PC′) for the proton transfers RC f PC f PC′ in model D2′.

TABLE 1: Relative Energies (kcal/mol) for Models A-D from QM/MM Optimizations with Basis Set BS1 (Single-Point Results with Basis Sets BS2 and BS3 Given in Parentheses)a model

RC

TS

PC

A B C D1 D2 D2

0.0 0.0 0.0 0.0 0.0 0.0

1.2 (2.9, 2.5) 2.2 (2.9, 2.2) 6.9 (7.1, 6.4) 4.6 (6.4, 5.7) 6.7 (7.2, 7.1) 6.6 (6.5, 7.5)

-21.1 (-30.7, -33.6) -4.8 (-11.6, -13.8) 5.4 (4.2, 3.2) -1.8 (-2.4, -5.2) -1.0 (-2.0, -5.5) -2.3 (-5.4, -6.3)

a Note that model A contains protonated Asp251, while models B-D contain neutral Asn251. The relative energies are thus directly comparable only among models B-D (D251N mutant).

difference between the wild-type enzyme and the mutant is caused by different conformations of the Asp251/Asn251 side chain. Experimental solvent kinetic isotope effects (SKIE) provide further mechanistic information. Measurements of the turnover rates in various protium-deuterium mixtures gave SKIE values

TABLE 2: Spin Densities on Selected Atoms and Groups at the Stationary Points model A B C D1 D2

RC TS PC RC TS PC RC TS PC RC TS PC RC TS PC

Fe

Op

Od

S-H

porphyrin

0.070 0.148 0.930 0.051 0.392 0.896 0.036 0.610 0.844 0.043 0.537 0.950 0.010 0.536 0.969

0.528 0.544 0.137 0.528 0.452 0.163 0.533 0.351 0.201 0.515 0.376 0.121 0.543 0.393 0.115

0.415 0.321 0.002 0.429 0.195 0.002 0.435 0.093 0.013 0.455 0.152 -0.002 0.447 0.119 -0.002

0.007 0.008 -0.008 0.007 0.001 -0.007 0.010 -0.004 -0.006 0.006 -0.006 -0.010 0.010 -0.005 -0.012

-0.012 -0.011 -0.062 -0.006 -0.035 -0.055 -0.006 -0.048 -0.053 -0.013 -0.056 -0.060 -0.000 -0.040 -0.061

(H/D) of 1.8 for wild-type P450cam and of 10 for the D251N mutant.9 Since these experiments were done under steady-state conditions, it was not possible to unambiguously identify the


Figure 15. Changes in the spin densities for model A (wild-type).

Figure 16. Changes in the Mulliken charges for model A (wild-type).

SCHEME 3: Orbital Diagram Showing the Electron Transfer from Fe to the O2 Moiety

key step responsible for the SKIE.9 Two more recent studies on the wild-type enzyme reported SKIE values (H/D) of 1.21 ( 0.08 at 298 K for the early steps of dioxygen activation (up to S5)54 and 1.8 at 190 K for the second proton transfer (Cpd 0 f Cpd I).55 These SKIE data indicate that, following reduction of S4, the second proton-transfer step is rate limiting for camphor hydroxylation in wild-type P450cam.55 The situation is different in the D251N mutant where the dramatic increase in the measured SKIE value suggests that the first proton-transfer step becomes the rate-limiting process in the catalytic cycle.9 According to proton inventory analysis, the numbers of protons involved in the rate-limiting step appears to be far larger in the mutant than in wild-type P450cam.9 On the other hand, in cryoreduction experiments with the D251N mutant, the conver-

Wang et al. sion of S5 to Cpd 0 was found to occur in D2O at a lower temperature than in H2O, contrary to expectation for a ratelimiting step involving proton transfer.22 This latter observation may, however, be less relevant for kinetics because the annealing temperature in cryoreduction experiments should mainly reflect the ease of diffusion. It is not straightforward to relate our computational results to the experimental evidence outlined above, as long as the experimental data do not refer to well-defined elementary reaction steps but to more complex processes (e.g., steady-state data). In all our models, the reactant is the reduced oxyheme complex (S5), with one (A-C) or two (D) active-site water molecules, and the product corresponds to Cpd 0 in different environments, with the transferred proton coming from Asp251 (A), Asn251 (B), Thr252 (C, direct mechanism), Asn251 (C, stepwise mechanism), Wat901 (D1 and D2), or Arg186 (D2′). We have not investigated how these proton sources are reprotonated from bulk water. We consider it possible that the observed SKIE will be affected by this reprotonation process, which should be facile in the case of Asp251, Asn251, and Arg186, since these residues are in contact with bulk water. In the case of the D251N mutant, such reprotonation would have to be invoked for SKIE in model C where the most probable scenario involves a back-flip of the Asn251 side chain to yield the model B conformer, which then undergoes a facile proton transfer via Wat901, without direct participation of solvent water. If an additional water molecule enters the cavity from the bulk (WatS, model D) a hydrogen-bond network is formed in the D251N mutant that extends from the distal Od atom to the surface of the protein via the chain composed of Thr252, Wat901, WatS, and Arg186. A two-step proton transfer along this chain is possible (model D2′), with facile reprotonation of the initially generated hydroxide anion (Wat901). In this case, the participation of a bulk water molecule (WatS) is expected to cause some SKIE. This scenario is consistent with different SKIE values for the first proton transfer (S5 f Cpd 0) in the wild-type enzyme and the D251N mutant because it is only in the latter case that the side chain of residue 251 flips and creates some space that can accommodate an extra water molecule. One should note, however, that uncoupling might be a competing side reaction under these circumstances, according to our calculations. The latest X-ray structure (2A1N) of the ferrous dioxygen complex (S4) of the D251N mutant lacks the two water molecules that form a hydrogen-bond network with dioxygen in the corresponding wild-type structure.46 This has been attributed to the orientation of the I helix, which does not open as much as in the wild-type structure so that water cannot move easily into the heme pocket of the D251N mutant.46 Our QM/ MM results show that the first proton transfer (S5 f Cpd 0) requires water in the active site. Its lack in the precursor complex (S4) of the D251N mutant46 indicates the need for water uptake during the reaction, implying some additional activation and a lower propensity for proton transfer in the mutant. A comparison of the X-ray structures of P450cam ferrous dioxygen complexes (S4) shows that there are no dramatic overall differences between the wild-type enzyme (1DZ8) and the D251N mutant (2A1N) (see Figure S12 of Supporting Information). We have based our QM/MM modeling both for the wild-type and the mutant on the 1DZ8 structure that already contains a water molecule in the active site (unlike 2A1N). Starting the modeling of proton transfer from the 2A1N structure

Proton Transfer in P450cam would have required the addition of active-site water in the QM/ MM setup phase (see above), which was considered less favorable. We end with a note of caution. As pointed out above, we have not studied the complete transformation S5 f Cpd 0 in models A-D because we have not considered the reprotonation of the proton sources in these models from bulk water. The computed energy profiles thus capture only part of the overall process. In spite of this limitation, we believe that they offer valuable qualitative insight, particularly with regard to the mechanism of proton transfer in the D251N mutant (since models B-D employ similar setups and are directly comparable). V. Conclusions We have investigated the first proton transfer in the catalytic cycle of cytochrome P450cam that leads to the formation of Cpd 0. Our QM/MM calculations address proton delivery in the doublet state through a hydrogen-bond network in the Asp251 channel, both for the wild-type enzyme and the D251N mutant. For wild-type P450cam, we find a concerted mechanism for proton delivery from protonated Asp251 to the FeOO moiety via Wat901 and Thr252, with a low barrier of about 1 kcal/mol and a high exothermicity of more than 20 kcal/mol (model A). In the case of the D251N mutant with neutral Asn251, there is also a facile proton-transfer route to Cpd 0 as long as there is no flip of the Asn251 residue away from the FeOO moiety (model B). Classical MD simulations suggest that such a flip should indeed occur, and according to the QM/MM calculations, the resulting conformer has no easy path to Cpd 0 so that the reduced oxyheme complex S5 may be relatively long-lived in this case (model C). In spite of the fact that no additional water molecule enters the active site during a 1-ns classical MD run, it is conceivable that this may eventually happen (model D). In this situation, the QM/MM calculations identify a viable path for proton transfer from Arg186 to FeOO via WatS, Wat901, and Thr252 (model D2′). The present QM/MM results emphasize the importance of the hydrogen-bond network in the Asp251 channel for efficient proton delivery. The barriers for concerted proton transfer are low as long as the network is intact (models A and B) even if protonated Asp251 is replaced by Asn251 as proton source. In this latter case, the nearby Arg186 residue activates Asn251 for proton transfer through strong electrostatic interactions (model B). Disruption of the network by a flip of the Asn251 side chain causes an increase in the initial proton-transfer barrier to 7 kcal/mol, makes this step endothermic by 5 kcal/mol and prevents a facile proton delivery from the bulk to the formed Thr252 anion (model C). In this case, a stepwise process with back flip and subsequent proton transfer appears to be the only viable option to complete the reaction, via an effective barrier of almost 10 kcal/mol (model C). Restoring the network by the entry of another water molecule into the actiVe site leads to the most faVorable proton deliVery in the D251N mutant, with a moderate barrier of ca. 7 kcal/mol and an exothermicity of -3 kcal/mol (model D2′). The first proton transfer in the P450cam cycle thus appears to be a delicately balanced process that is controlled and tuned by the protein/solvent environment (Asp251 vs Asn251, Arg186, WatS). The QM/MM results are consistent with the longer lifetime of the reduced oxyheme complex S5 in the mutant D251N compared with the wild-type enzyme. If the preferred conformer of the mutant indeed has a flipped Asn251 side chain as

J. Phys. Chem. B, Vol. 112, No. 16, 2008 5137 indicated by the classical MD simulations, the proton transfer in the Asn251 channel requires either a back flip of the Asn251 side chain or the participation of an extra water molecule in the active site, and even then the barriers are higher than in wild-type P450cam. The lifetime of S5 is thus effectively increased in the mutant by blocking a facile proton transfer toward Cpd 0. Acknowledgment. The research at HU was supported by a grant from BMBF (DIP-F.7.1) and the Israeli Science Foundation. We thank S. G. Sligar and T. L. Poulos for helpful comments. Supporting Information Available: Computational details on system setup and QM/MM methodology, energy profiles from pathway calculations, detailed results for spin densities and partial charges, overlay of X-ray structures. This information is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) de Grotthuss, C. J. T. Anal. Chim. 1806, 58, 54-73. (2) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. Nature 1999, 397, 601-604. (3) Rini, M.; Magnes, B.-Z.; Pines, E.; Nibbering, E. T. J. Science 2003, 301, 349-352. (4) Omta, A. W.; Kropman, M. F.; Woutersen, S.; Bakker, H. J. Science 2003, 301, 347-349. (5) Taraphder, S.; Hummer, G. J. Am. Chem. Soc. 2003, 125, 39313940. (6) For reviews see: (a) Denisov, I. G.; Makris, T. M.; Sligar, S. G.; Schlichting, I. Chem. ReV. 2005, 105, 2253-2277. (b) Shaik, S.; Kumar, D.; de Visser, S. P.; Altun, A.; Thiel, W. Chem. ReV. 2005, 105, 22792328. (c) Sono, M.; Roach, M. P.; Coulter, E. D.; Dawson, J. H. Chem. ReV. 1996, 96, 2841-2887. (d) Meunier, B.; de Visser, S. P.; Shaik, S. Chem. ReV. 2004, 104, 3947-3980. (e) Omura, T. Biochem. Biophys. Res. Commun. 1999, 266, 690-698. (f) Loew, G. H.; Harris, D. L. Chem. ReV. 2000, 100, 407-419. (g) Poulos, T. L.; Cupp-Vickery, J.; Li, H. Structural Studies on Prokaryiotiuc Cytochromes P450. In Cytochrome P450: Structure, Mechanism and Biochemistry, 2nd ed.; Ortiz de Montellano, P. R., Ed.; Plenum Press: New York, 1995. (h) Poulos, T. L. Biochem. Biophys. Res. Commun. 2005, 338, 337-345. (i) Poulos, T. L. Philos. Trans. R. Soc. London, Ser. A 2005, 363, 793-806. (7) Schlichting, I.; Berendzen, J.; Chu, K.; Stock, A. M.; Maves, S. A.; Benson, D. E.; Sweet, R. M.; Ringe, D.; Petsko, G. A.; Sligar, S. G. Science 2000, 287, 1615-1622. (8) Kamachi, T.; Yoshizawa, K. J. Am. Chem. Soc. 2003, 125, 46524661. (9) Vidakovic, M.; Sligar, S. G.; Li, H.; Poulos, T. L. Biochemistry 1998, 37, 9211-9219. (10) Guallar, V.; Friesner, R. A. J. Am. Chem. Soc. 2004, 126, 85018508. (11) Oprea, T. I.; Hummer, G.; Garcia, A. E. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 2133-2138. (12) Harris, D. L.; Loew, G. H. J. Am. Chem. Soc. 1996, 116, 63776387. (13) Shimada, H.; Makino, R.; Unno, M.; Horiuchi, T.; Ishimura, Y. In Cytochrome P450 8th International Conference; Libbey, J., Ed.; Eurotext: Paris, 1994; p 299. (14) Sligar, S. G.; Makris, T. M. Private communication, February 2005. (15) Li, H.; Robertson, A. D.; Jensen, J. H. Proteins 2005, 61, 704721. (16) Gerber, N. C.; Sligar, S. G. J. Am. Chem. Soc. 1992, 114, 87428743. (17) Lounnas, V.; Wade, R. C. Biochemistry 1997, 36, 5402-5417. (18) Poulos, T. L.; Finzel, B. C.; Howard, A. J. J. Mol. Biol. 1987, 195, 687-700. (19) Zheng, J.; Wang, D.; Thiel, W.; Shaik, S. J. Am. Chem. Soc. 2006, 128, 13204-13215. (20) Harris, D. L.; Loew, G. H. J. Am. Chem. Soc. 1998, 120, 89418948. (21) Ogliaro, F.; de Visser, S. P.; Cohen, S.; Sharma, P. K.; Shaik, S. J. Am. Chem. Soc. 2002, 124, 2806-2817. (22) Davydov, R.; Makris, T. M.; Kofman, V.; Werst, D. E.; Sligar, S. G.; Hoffman, B. M. J. Am. Chem. Soc. 2001, 123, 1403-1415. (23) Harris, D. L. J. Inorg. Biochem. 2002, 91, 568-585. (24) Kamachi, T.; Yoshizawa, K. J. Am. Chem. Soc. 2003, 125, 46524661.

5138 J. Phys. Chem. B, Vol. 112, No. 16, 2008 (25) Kumar, D.; Hirao, H.; de Visser, S. P.; Zheng, J.; Wang, D.; Thiel, W.; Shaik S. J. Phys. Chem. B 2005, 109, 19946-19951. (26) Vaz, A. D. N.; Mcginnity, D. F.; Coon, M. J. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 3555-3560. (27) Jin, S.; Makris, T. M.; Bryson, T. A.; Sligar, S. G.; Dawson, J. H. J. Am. Chem. Soc. 2003, 125, 3406-3407. (28) Jin, S.; Bryson, T. A.; Dawson, J. H. J. Biol. Inorg. Chem. 2004, 9, 644-653. (29) Bach, R. D.; Dmitrenko, O. J. Am. Chem. Soc. 2006, 128, 14741488. (30) Scho¨neboom, J. C.; Lin, H.; Reuter, N.; Thiel, W.; Cohen, S.; Ogliaro, F.; Shaik, S. J. Am. Chem. Soc. 2002, 124, 8142-8151. (31) MacKerell, A. D., Jr.; Bashford, D.; Bellott, R. L.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E., III; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586. (32) Brooks, B. R.; Burccoleri, R. E.; Olafson, B. D.; States, D. J.; Karplus, M. J. Comput. Chem. 1983, 4, 187-217. (33) Bakowies, D.; Thiel, W. J. Phys. Chem. 1996, 100, 10580-10594. (34) (a) de Vries, A. H.; Sherwood, P.; Collins, S. J.; Rigby, A. M.; Rigutto, M.; Kramer, G. J. J. Phys. Chem. B 1999, 103, 6133-6141. (b) Sherwood, P.; de Vries, A. H.; Collins, S. J.; Greatbanks, S. P.; Burton, N. A.; Vincent, M. A.; Hillier, I. H. Faraday Discuss. 1997, 106, 79-92. (35) Ahlrichs, R.; Ba¨r, M.; Ha¨ser, M.; Horn, H.; Ko¨lmel, C. Chem. Phys. Lett. 1989, 162, 165-169. (36) Smith, W.; Forester, T. J. Mol. Graph. 1996, 14, 136-141. (37) Sherwood, P.; de Vries, A. H.; Guest, M. F.; Schreckenbach, G.; Catlow, C. R. A.; French, S. A.; Sokol, A. A.; Bromley, S. T.; Thiel, W.; Turner, A. J.; Billeter, S.; Terstegen, F.; Thiel, S.; Kendrick, J.; Rogers, S. C.; Casci, J.; Watson, M.; King, F.; Karlsen, E.; Sjøvoll, M.; Fahmi, A.; Scha¨fer, A.; Lennartz, C. THEOCHEM 2003, 632, 1-28. (38) Billeter, S. R.; Turner, A. J.; Thiel, W. Phys. Chem. Chem. Phys. 2000, 2, 2177-2186. (39) Becke, A. D. J. Chem. Phys. 1993, 98, 5648-5652.

Wang et al. (40) Hay, P. J.; Wadt, W. R. J. Chem. Phys. 1985, 82, 299-310. (41) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257-2261. (42) Scha¨fer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 25712577. (43) Double-ζ basis set, augmented with two p, one d, and one f function (DZpdf, exponents: 0.134915(p), 0.041843(p), 0.1244(d), and 1.339(f)) with the contraction scheme (14s11p6d1f/8s7p4d1f). (44) Scha¨fer, A.; Huber, C.; Ahlrichs, R. J. Chem. Phys. 1994, 100, 5829-5835. Triple-ζ quality basis set, augmented with one p, one d, and one f function (exponents: 0.134915(p), 0.091(d), and 1.598(f)). (45) (a) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650-654. (b) Blaudeau, J.-P.; McGrath, M. P.; Curtiss, L. A.; Radom, L. J. Chem. Phys. 1997, 107, 5016-5021. (46) Nagano, S.; Poulos, T. L. J. Biol. Chem. 2005, 280, 3165931663. (47) Davydov, R.; Satterlee, J. D.; Fujii, H.; Sauer-Masarwa, A.; Busch, D. H.; Hoffman, B. M. J. Am. Chem. Soc. 2003, 125, 16340-16346. (48) Davydov, R.; Perera, R.; Jin, S.; Yang, T.-C.; Bryson, T. A.; Sono, M.; Dawson, J. H.; Hoffman, B. M. J. Am. Chem. Soc. 2005, 127, 14031413. (49) Hammond, G. S. J. Am. Chem. Soc. 1955, 77, 334-338. (50) Erman, J. E.; Vitello, L. B.; Miller, M. A.; Shaw, A.; Brown, K. A.; Kraut, J. Biochemistry 1993, 32, 9798. (51) Rodrı´gues-Lo´pez, J. N.; Lowe, D. J.; Hernańdez-Ruiz, J.; Hiner, A. N. P.; Garcıá-Cańovas, F.; Thorneley, R. N. F. J. Am. Chem. Soc. 2001, 123, 11838-11847. (52) Ionization Constants of Organic Acids in Solution, IUPAC Chemical Data Series No. 23; Serjeant, E. P., Dempsey, B., Eds.; Pergamon Press: Oxford, UK, 1979. (53) Benson, D. E.; Suslick, K. S.; Sligar, S. G Biochemistry 1997, 36, 5104-5107. (54) Purdy, M. M.; Koo, L. S.; Ortiz de Montellano, P. R.; Klinman, J. P. Biochemistry 2006, 45, 15793-15806. (55) Kim, S. H.; Yang, T.-C.; Perera, R.; Jin, S.; Bryson, T. A.; Sono, M.; Davydov, R.; Dawson, J. H.; Hoffman, B. M. Dalton Trans. 2005, 3464-3469.

Quantum and Molecular Mechanical Study of the First Proton Transfer

Recommend Documents