MM Study of Mechanisms for Compound I Formation in the

In mechanism II, an initial O−O bond cleavage followed by a concomitant proton and electron .... Journal of Chemical Theory and Computation 2012 8 (...
0 downloads 0 Views 510KB Size
Download PDF
Published on Web 09/14/2006

QM/MM Study of Mechanisms for Compound I Formation in the Catalytic Cycle of Cytochrome P450cam Jingjing Zheng,† Dongqi Wang,† Walter Thiel,*,† and Sason Shaik*,‡ Contribution from the 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 May 17, 2006; E-mail: [email protected]; [email protected]

Abstract: In the catalytic cycle of cytochrome P450cam, after molecular oxygen binds as a ligand to the heme iron atom to yield a ferrous dioxygen complex, there are fast proton transfers that lead to the formation of the active species, Compound I (Cpd I), which are not well understood because they occur so rapidly. In the present work, the conversion of the ferric hydroperoxo complex (Cpd 0) to Cpd I has been investigated by combined quantum-mechanical/molecular-mechanical (QM/MM) calculations. The residues Asp251 and Glu366 are considered as proton sources. In mechanism I, a proton is transported to the distal oxygen atom of the hydroperoxo group via a hydrogen bonding network to form protonated Cpd 0 (prot-Cpd0: FeOOH2), followed by heterolytic O-O bond cleavage that generates Cpd I and water. Although a local minimum is found for prot-Cpd0 in the Glu366 channel, it is very high in energy (more than 20 kcal/mol above Cpd 0) and the barriers for its decay are only 3-4 kcal/mol (both toward Cpd 0 and Cpd I). In mechanism II, an initial O-O bond cleavage followed by a concomitant proton and electron transfer yields Cpd I and water. The rate-limiting step in mechanism II is O-O cleavage with a barrier of about 13-14 kcal/mol. According to the QM/MM calculations, the favored low-energy pathway to Cpd I is provided by mechanism II in the Asp251 channel. Cpd 0 and Cpd I are of similar energies, with a slight preference for Cpd I.

1. Introduction

The cytochromes P450 constitute a superfamily of monooxygenases found in nearly all living species. These enzymes perform a variety of functions, such as detoxification of xenobiotics as well as biosynthesis of sex hormones, musclerelaxing, antiflammatory, and antihypertensive compounds,1 and they catalyze many types of reactions including aliphatic and aromatic hydroxylation, epoxidation, and heteroatom oxidation. In these enzymes the active site is the heme unit that consists of an iron protoporphyrin IX complex with a cysteine as the proximal axial ligand of the heme iron. It is commonly accepted that the iron-oxo species called Compound I (Cpd I) is the active oxidant. This species has not been trapped experimentally in P450cam because of the very fast reaction steps occurring after O2 binding. Reports of an observation of Cpd I by means of rapid scan absorption spectroscopy2 and cryogenic X-ray diffraction3 have been questioned by ENDOR and EPR spec† ‡

Max-Planck-Institut fu¨r Kohlenforschung. The Hebrew University.

(1) (a) Ortiz de Montellano, P. R., Ed. Cytochrome P450: Structure, Mechanisms, and Biochemistry, 2nd ed.; Plenum Press: New York, 1995. (b) Ortiz de Montellano, P. R., Ed. Cytochrome P450: Structure, Mechanisms, and Biochemistry, 3rd ed.; Plenum Press: New York, 2005. (c) Makris, T. M.; Denisov, I.; Schlichting, I.; Sligar, S. G. In ref 1b, pp 149-182. (2) Egawa, T.; Shimada, H.; Ishimura, Y. Biochem. Biophys. Commun. 1994, 201, 1464-1469. (3) Schlichting, I.; Berendzen, J.; Chu, K.; Stock, A. M.; Maves, S. A.; Benson, D. A.; Sweet, R. M.; Ringe, D.; Petsko, G. A.; Sligar, S. G. Science 2000, 287, 1615-1622. 13204

9

J. AM. CHEM. SOC. 2006, 128, 13204-13215

troscopic studies4 which identified the ferric hydroperoxo complex (so-called Compound 0, Cpd 0) as the last observable intermediate in the catalytic cycle and showed that Cpd I of P450cam does not accumulate in a frozen solution at a detection temperature of 200 K (or even lower). However, Cpd I has been observed spectroscopically in other P450 enzymes (CYP1195 and substrate-free CYP1016). On the theoretical side, the protonation steps from the ferricperoxo complex to Cpd 0 and subsequently Cpd I have been studied extensively.7-12 It is generally believed that there are three possible proton sources, namely Glu366, Asp251, and the propionate side chains of the heme.9 The carboxylic acid side chain of Glu366 is linked to the hydroxyl group of Thr252 by a chain of water molecules in the available crystal structures (e.g., 1DZ8).3 Asp251 is connected to the protein surface via a salt bridge to Arg186 so that its carboxylic acid side chain points away from the heme and the active site.3 However, as proposed (4) Davydov, R.; Makris, T. M.; Kofman, V.; Werst, D. E.; Sligar, S. G.; Hoffman, B. M. J. Am. Chem. Soc. 2001, 123, 1403-1415. (5) Kellner, D. G.; Hung, S.-C.; Weiss, K. E.; Sligar, S. G. J. Biol. Chem. 2002, 277, 9641-9644. (6) Spolitak, T.; Dawson, J. H.; Ballou, D. P. J. Biol. Chem. 2005, 280, 2030020309. (7) Harris, D. L.; Loew, G. H. J. Am. Chem. Soc. 1998, 120, 8941-8948. (8) Kamachi, T.; Shiota, Y.; Ohta, T.; Yoshizawa, K. Bull. Chem. Soc. Jpn. 2003, 76, 721-732. (9) Taraphder, S.; Hummer, G. J. Am. Chem. Soc. 2003, 125, 3931-3940. (10) Guallar, V.; Friesner, R. A. J. Am. Chem. Soc. 2004, 126, 8501-8508. (11) Kumar, D.; Hirao, H.; de Visser, S. P.; Zheng, J. J.; Wang, D. Q.; Thiel, W.; Shaik, S. J. Phys. Chem. B 2005, 109, 19946-19951. (12) Shaik, S.; Kumar, D.; de Visser, S. P.; Altun, A.; Thiel, W. Chem. ReV. 2005, 105, 2279-3935. 10.1021/ja063439l CCC: $33.50 © 2006 American Chemical Society

QM/MM Study of Mechanisms for Compound I Formation

ARTICLES

Scheme 1. Two Mechanisms for Cpd I Formationa

a

The intermediate in mechanism II is a hydrogen-bonded complex.

in early experimental work,13,14 the carboxylic acid side chain of Asp251 may become linked via Wat901 to the hydroxyl group of Thr252 by an internal rotation of reasonable amplitude, and the previous9 and present molecular dynamics (MD) simulations show that this is feasible with low energetic cost. By contrast, there is no direct water channel connecting the propionate side chains of the heme with the FeOOH group since the substrate (camphor) blocks the channel. Therefore, we consider only Glu366 and Asp251 as possible proton sources in the present work. In most of the available literature, the mechanism of Cpd I formation from Cpd 0 is thought to involve the initial protonation of the distal oxygen atom followed by heterolytic O-O bond cleavage to generate Cpd I and water (Mechanism I in Scheme 1). This O-O activation mechanism is a general paradigm in the P450 family of enzymes where a pair of functional acid/alcohol groups is highly conserved and believed to cause heterolytic O-O cleavage.1c The conversion of Cpd 0 to Cpd I has been studied theoretically by several groups,7,10,12,15,16 but the results are not clear-cut and partly inconsistent with each other. Recently, we have found protonated Cpd 0 (prot-Cpd0) as an intermediate in DFT and preliminary QM/MM calculations.11 Here we present a more complete QM/MM study of Cpd I formation in P450cam and propose an alternative mechanism, where an initial O-O bond cleavage that is homolytic in nature is instantaneously followed by a proton transfer with a concomitant electron transfer to yield Cpd I and water (Mechanism II in Scheme 1). The issue of heterolytic vs homolytic dissociation of the O-O bond in the FeOOH moiety has been discussed for a number of heme and non-heme systems.17,18 Several DFT model calcula(13) Vidakovic, M.; Sligar, S. G.; Li, H.; Poulos, T. L. Biochemistry 1998, 37, 9211-9219. (14) Gerber, N. C.; Sligar, S. G. J. Biol. Chem. 1994, 269, 4260-4266. (15) Ogliaro, F.; de Visser, S. P.; Cohen, S.; Sharma, P. K.; Shaik, S. J. Am. Chem. Soc. 2002, 124, 2806-2817. (16) Kamachi, T.; Yoshizawa, K. J. Am. Chem. Soc. 2003, 125, 4652-4661. (17) For recent reviews, see: (a) Denisov, I. G.; Makris, T. M.; Sligar, S. G.; Schlichting, I. Chem. ReV. 2005, 105, 2253-2277. (b) Decker, A.; Solomon, E. I. Curr. Opin. Chem. Biol. 2005, 9, 152-163.

tions in the literature report on homolytic O-O bond scission. For example, porphyrin degradation in heme oxygenase was predicted to proceed via O-O bond homolysis,19,20 and it was suggested that, in the absence of Cpd I, the bound radical resulting from O-O homolysis of Cpd 0 could also lead to minor amounts of substrate oxidation (the major reaction would be heme degradation).21 A recent DFT model study proposed that all substrate hydroxylation originates in Cpd 0 via a “somersault” mechanism which involves an initial O-O homolysis of Cpd 0 yielding a hydrogen-bonded “inverted metastable hydroperoxide” [FeO‚‚‚HO] that is supposed to act as an oxidant in the following hydroxylation step.22 In this paper, we address the question of heterolytic vs homolytic O-O cleavage during the conversion of Cpd 0 to Cpd I at the QM/ MM level, in a realistic enzyme environment. Another issue concerns the energy of Cpd I relative to Cpd 0. In early work,7,15 a proton was added to a minimal Cpd 0 model which led to spontaneous conversion to Cpd I and water in a highly exothermic process (by more than 300 kcal/mol). This basic model is clearly not realistic and needs to be refined to account for the situation in the enzyme. Subsequent DFT calculations on a model including the ferric hydroperoxo active site, Asp251, Thr252, and two water molecules favored Cpd I by only a small margin (5.5 kcal/mol).16 In a later QM/MM study of the full enzyme without solvent water, constrained geometry (18) For O-O homolysis in P450 and model systems, see: (a) White, R. E.; Sligar, S. G.; Coon, M. J. J. Biol. Chem. 1980, 255, 11108-11111. For competition between O-O homolysis and heterolysis in heme enzymes, see: (b) Allentoff, A. J.; Bolton, J. L.; Wilks, A.; Thompson, J. A.; Ortiz de Montellano, P. R. J. Am. Chem. Soc. 1992, 114, 9744-9749. (c) Groves, J. T.; Watananbe, Y. J. Am. Chem. Soc. 1988, 110, 8443-8452. (d) Yamaguchi, K.; Watanabe, Y.; Morishima, I. J. Am. Chem. Soc. 1993, 115, 4058-4065. (e) Shimizu, T.; Murakami, Y.; Hatano, M. J. Biol. Chem. 1994, 269, 13296-13304. (f) Stephenson, N. A.; Bell, A. T. Inorg. Chem. 2006, 45, 2758-2766. (19) Sharma, P. K.; Kevorkiants, R.; de Visser, S. P.; Kumar, D.; Shaik, S. Angew. Chem., Int. Ed. 2004, 43, 1129-1132. (20) Kumar, D.; de Visser, S. P.; Shaik, S. J. Am. Chem. Soc. 2005, 127, 82048213. (21) Derat, E.; Kumar, D.; Hirao, H.; Shaik, S. J. Am. Chem. Soc. 2006, 128, 973-989. (22) Bach, R. D.; Dmitrenko, O. J. Am. Chem. Soc. 2006, 128, 1474-1488. J. AM. CHEM. SOC.

9

VOL. 128, NO. 40, 2006 13205

Zheng et al.

ARTICLES

optimizations followed by single-point calculations with a QM region consisting of the hydroperoxo active site, Thr252, Glu366, and a chain of four water molecules predicted Cpd I to be more stable than Cpd 0 by 80 kcal/mol.10 A more recent DFT study employed a larger QM model which included the heme and its ligands, Asp251, Thr252, Glu366, and five crystallographic water molecules: Cpd I and Cpd 0 were found to be very close in energy, in agreement with preliminary QM/MM results.11 It is obvious from this brief overview that the calculated relative energies depend on the chosen model system and the computational methodology. In view of the conflicting results in the published literature,7,10-12,15,16 the present paper also reports a detailed QM/MM investigation of this issue. 2. Methods The initial coordinates for our Cpd 0 systems were taken from the X-ray crystal structure 1DZ83 (ferrous-dioxygen complex). The same solvation and protonation schemes were applied as before.23-25 Since Glu366 and Asp251 are possible proton sources,9,13,14 the corresponding two protonation schemes were adopted in the framework of the standard setup used previously.23-25 One scheme employs protonated Glu366 and deprotonated Asp251, and the other one deprotonated Glu366 and protonated Asp251, to study the proton-transfer process in the Glu366 and Asp251 channels, respectively. Both systems consist of 24 988 atoms and include 5895 TIP3P water molecules.26 These solvated systems were relaxed by performing energy minimizations and molecular dynamic (MD) simulations at the MM level by using the CHARMM22 force field27 as implemented in the CHARMM program.28 The heme unit with the Cys357 and OOH ligands and the outer 8 Å of solvent layer were kept fixed during the classical energy minimizations and MD simulations. Classical MD simulations were also performed to calculate free energy profiles for Asp251 side chain rotation and for the torsion of the salt bridge between Asp251 and Arg186 by using the adaptive umbrella sampling technique29 as implemented in the CHARMM program.28 The chosen QM/MM methodology is analogous to that used in our previous studies.23-25 Here we briefly mention some aspects relevant to the present work. Minimized snapshots from the MD trajectories were taken as the initial structures for QM/MM optimization. In the QM/MM calculations, the QM part was treated by the UB3LYP density functional method, and the MM part was described by the CHARMM22 force field. An electronic embedding scheme30 was adopted in the QM/ MM calculations; i.e. the MM charges were incorporated into the oneelectron Hamiltonian of the QM calculation, and the QM/MM electrostatic interactions were evaluated from the QM electrostatic potential and the MM partial charges. No cutoffs were introduced for the nonbonding MM and QM/MM interactions. Hydrogen link atoms with the charge shift model31,32 were employed to treat the QM/MM boundary. The TURBOMOLE program33 was used for the QM treatment in the QM/MM as well as in the pure QM calculations. The (23) Scho¨neboom, J. C.; Lin, H.; Reuter, N.; Thiel, W.; Cohen, S.; Ogliaro, F.; Shaik, S. J. Am. Chem. Soc. 2002, 124, 8142-8151. (24) Scho¨neboom, J. C.; Thiel, W. J. Phys. Chem. B 2004, 108, 7468-7478. (25) Altun, A.; Thiel, W. J. Phys. Chem. B 2005, 109, 1268-1280. (26) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (27) MacKerell, A. D., Jr., et al. J. Phys. Chem. B 1998, 102, 3586. (28) Brooks, B. R.; Burccoleri, R. E.; Olafson, B. D.; States, D. J.; Karplus, M. J. Comput. Chem. 1983, 4, 187-217. (29) (a) Bartels, C.; Karplus, M. J. Comput. Chem. 1997, 18, 1450-1462. (b) Bartels, C.; Karplus, M. J. Phys. Chem. B 1998, 102, 865-880. (c) Schaefer, M.; Bartels, C.; Karplus, M. J. Mol. Biol. 1998, 284, 835-848. (30) Bakowies, D.; Thiel, W. J. Phys. Chem. 1996, 100, 10580-10594. (31) 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. (32) 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. (33) Ahlrichs, R.; Ba¨r, M.; Ha¨ser, M.; Horn, H.; Ko¨lmel, C. Chem. Phys. Lett. 1989, 162, 165-169. 13206 J. AM. CHEM. SOC.

9

VOL. 128, NO. 40, 2006

Table 1. List of QM Regions QM regionsa

EQ

DQ1

DQ2

DQ3

DQ4

porphyrin-FeOOH SH ligand Cys357ligandb Glu366(CH3COOH) Asp251(CH3COOH) Arg186(CH3CH2NHC(NH2)2) Thr252(CH3CH2OH) Wat901 four crystallographic watersc total number of atoms

x x

x x

x

x x

x

x

x

x x

x x

x x x x

x x x x

62

78

79

95

x x x 71

x

x

a In this table, the first column refers to the reactant, Cpd 0. The representation of the residues (Glu366, Asp251, Arg186, and Thr252) is shown in parentheses. b This proximal ligand includes the whole Cys357, the CO group of Leu356, and the NH-CRH unit of Leu358. c The four crystallographic waters are Wat523, Wat566, Wat687, and Wat902 in 1DZ8.

CHARMM22 force field was run through the DL_POLY34 program to handle the MM part of the systems. The QM/MM calculations were performed with the ChemShell package35 that integrates the TURBOMOLE and DL_POLY programs and also performs geometry optimization with the HDLC optimizer.36 Five QM regions were adopted in the QM/MM calculations (Table 1). The QM region EQ represents the Glu366 channel, while DQ1DQ4 represent the Asp251 channel. All the energy scans were performed with either EQ or DQ1. The QM regions DQ2-DQ4 were only used for optimizations and single-point energy evaluations of species obtained from DQ1 energy scans. Four basis sets were employed in this study. B1 [LACVP(Fe)/6-31G(rest)] was used for geometry optimization, and B1a [LACVP*(Fe)/6-31+G*(S, pyrrole N and OOH2 ligand)/6-31G(rest)], for improving the B1 geometry of prot-Cpd0. The other two basis sets, B2 [LACVP*(Fe)/6-31+G*(rest)] and B3 [LACV3P++**(Fe)/6-311++G**(rest)] only served for single-point calculations. Further computational details are documented in the Supporting Information.

3. Results and Discussion

3.1. Rotation of Asp251 Side Chain. In the P450cam crystal structure 1DZ8,3 Asp251 and Arg186 are engaged in a salt bridge, and there is no direct connection between Asp251 and Wat901. In the D251N mutant, the replacement of a carboxylate oxygen atom in Asp251 by an NH2 group in Asn251 destroys this salt bridge, and the corresponding X-ray structure 2A1N37 shows a 25° rotation of the Asn251 side chain toward Wat901, thus establishing a connection to the active site. This suggests that protonation of Asp251 in the P450cam case (with concomitant breaking of the salt bridge) may lead to an analogous rotation and the formation of a proton-transfer channel from Asp251 via the only crystallographic water molecule in this region (Wat901) to the active site.9,13 For a more quantitative assessment, we have performed classical MD simulations to get free energy profiles for the relevant torsional motions in the systems with protonated and deprotonated Asp251 (see Section 2) which are labeled as 1AspP and 2Asp, respectively. Three dihedral angles were sampled (see Figure 1): χ1 and χ2 describe the orientation of the Asp251 side chain, while χ3 indicates the arrangement in the salt bridge. The (34) Smith, W.; Forester, T. J. Mol. Graph. 1996, 14, 136-141. (35) (a) Sherwood, P., et al. J. Mol. Struct. (THEOCHEM) 2003, 632, 1. (b) ChemShell is a modular QM/MM program developed in the European QUASI project under the coordination of P. Sherwood (see http:// www.cse.clrc.ac.uk/qcg/chemshell). (36) Billeter, S. R.; Turner, A. J.; Thiel, W. Phys. Chem. Chem. Phys. 2000, 2, 2177-2186. (37) Nagano, S.; Poulos, T. L. J. Biol. Chem. 2005, 280, 31659-31663.

QM/MM Study of Mechanisms for Compound I Formation

ARTICLES

Figure 1. Definition of three dihedral angles χ1, χ2, and χ3.

resulting free energy profiles are shown in Figure 2, and the equilibrium dihedral angles are listed in Table 2. The main conclusions from the MD simulations are as follows: (i) The experimental χ1 value from 1DZ8 is reproduced very well in the 2Asp system with deprotonated Asp251, while there is a discrepancy of about 10° in the protonated 1AspP system. (ii) Considering both χ1 and χ2, the X-ray structure 1DZ8 and the 2Asp system have a similar Asp251 side chain orientation, whereas 1AspP is quite different. (iii) The free energy cost for χ1 torsion by 10° away from its optimum value is negligibly small (less than 1 kcal/mol), while the free energies for the χ2 torsion between -100° and -150° remain within 4 kcal/mol in each case (Figure 2). Hence, Asp251 side chain rotation is quite facile within reasonable amplitudes, both for the protonated and deprotonated case. (iv) Disruption of the salt bridge by protonation (2Asp f 1AspP) and by mutation in the D251N mutant (1DZ8 f 2A1N) causes a similar reorientation (χ2) of the side chain in residue 251 toward Wat901. For the purpose of comparison, we have also considered two analogous systems with two additional water molecules placed in the vicinity of Asp251 and Wat901 (where there seems to be some empty space in 1DZ8). Corresponding MD simulations such as those for 1AspP and 2Asp were performed for these two analogous systems and led to similar conclusions. In summary, the MD simulations show that the salt bridge between deprotonated Asp251 and Arg186 observed in the P450cam X-ray structure 1DZ8 is not as rigid as one might expect. Asp251 side chain rotations with considerable amplitude (χ2) are energetically feasible that move the side chain toward Wat901, and upon protonation of Asp251 the system automatically relaxes to such an orientation, establishing a connection of the carboxylic acid side chain via Wat901 to the active site. The QM/MM calculations for the Asp251 channel will therefore start from such conformations with a rotated side chain. 3.2. Mechanism I: Proton Assisted Heterolytic O-O Bond Cleavage. Mechanism I assumes that a proton is transported to the distal oxygen atom of the hydroperoxo group via a hydrogen bonding network to form a protonated Cpd 0 (prot-Cpd0: FeO-

Figure 2. Free energy profiles for Asp251 side chain rotation and for Asp251‚ ‚‚Arg186 salt bridge torsion (see text).

OH2),11 followed by heterolytic O-O bond cleavage that generates Cpd I and water (see Scheme 1). A key species in this mechanism is prot-Cpd0 which has frequently been discussed in the literature and has recently been identified computationally.11 In the present work we have investigated mechanism I in both channels (Glu366 and Asp251) by QM/MM calculations, but prot-Cpd0 was found as an intermediate only in the Glu366 channel (QM region EQ). In the Asp251 channel, we have chosen several different reaction coordinates to convert Cpd 0 to prot-Cpd0 by proton transfer from Asp251 to the distal oxygen atom of the hydroperoxo group, but all energy scans led to continuously increasing energy profiles. By applying J. AM. CHEM. SOC.

9

VOL. 128, NO. 40, 2006 13207

Zheng et al.

ARTICLES Table 2. Dihedral Angles χ1, χ2, and χ3 in the Crystal Structures 1DZ8 and 2A1N (D251N Mutant) and in the Free-Energy Minima of the Two Model Systems (Obtained by Adaptive Umbrella Sampling)

b

χ1 χ2b χ3b

1DZ8

2A1Na

1AspP

2Asp

75.7 -150.8 38.5

75.9 -123.7

65.0 -105.0 5.0

75.0 -135.0 15.0

a 2A1N is the ferrous dioxygen complex of the D251N mutant.37 b See Figure 1: χ1 CG-CB-CA-C, χ2 OD2-CG-CB-CA (ND2-CG-CBCA in 2A1N), χ3 OD1-OD2-HH12-HH22 (OD1-OD2-NH1-NH2 in 1DZ8).

suitable constraints, it was possible to generate artificial “protCpd0-like” structures which, however, always decayed to Cpd 0 during subsequent unconstrained optimizations. Calculations were performed with various QM regions (DQ1 and DQ3) and conformations in the Asp251 channel, but the desired prot-Cpd0 minimum could not be located in any case. Similar problems have also been reported in previous QM/MM calculations.10 Given this situation, we shall not consider this Asp251 channel any more in the remainder of section 3.2 but focus on the Glu366 channel. Cpd 0 has a doublet ground state which is more stable than the quartet state by about 13-14 kcal/mol both in the gas phase15,38 and the enzyme environment. Therefore, we only treat the doublet state. In the Glu366 channel, a proton needs to be transferred from the Glu366 carboxylic group through a chain of three water molecules (Wat523, Wat566, and Wat902) and Thr252, to the FeOOH group of Cpd 0. We performed energy scans for successive local proton transfers that generate hydronium ion species. The corresponding energy profiles were uphill, and the formed hydronium ion species were not stable. By applying judicious constraints, it was possible to construct a sensible path, but unconstrained optimization of species along this path did not yield any local minima of intermediates but led either back to the reactant (Cpd 0) or to the desired product (prot-Cpd0). This would seem to indicate a concerted proton transfer in the Glu366 channel. Contrary to the Asp251 channel, prot-Cpd0 is found to be a genuine minimum in the Glu366 channel. Its geometry has been fully optimized with basis sets B1 and B1a (see Figure 3). When enlarging the basis from B1 to B1a, the O-O and Fe-S bond lengths are shortened by 0.06 and 0.03 Å, respectively. Energetically, prot-Cpd0 lies much higher than Cpd 0 by 33.1 (35.9) kcal/mol with B1 (B1a). The unpaired electron is mainly located on iron (spin density FFe ) 1.120). The water moiety in FeOOH2 has almost zero spin density (-0.005) but a high positive charge (0.454). It is stabilized by a strong hydrogen bond with Thr252 as well as long-range electrostatic interactions with the negatively charged Glu366 and Asp251 residues. Singlepoint QM/MM relative energies with the larger B2 and B3 basis sets are 41.8 and 41.1 kcal/mol, respectively. Single-point QM energy calculations in the gas phase at QM/MM optimized geometries yield relative energies of 22.0/30.4/30.1 kcal/mol with B1/B2/B3. Hence, the enzyme environment destabilizes prot-Cpd0 by about 10 kcal/mol relative to Cpd 0, which is consistent with the previous gas-phase DFT results.11 Although prot-Cpd0 is a genuine minimum structure, it can easily decay to Cpd I via heterolytic O-O bond cleavage with (38) Kamachi, T.; Shiota, Y.; Ohta, T.; Yoshizawa, K. Bull. Chem. Soc. Jpn. 2003, 76, 721-732. 13208 J. AM. CHEM. SOC.

9

VOL. 128, NO. 40, 2006

a barrier of 4.0 kcal/mol. The fully optimized transition state for O-O cleavage is shown in Figure 3 as TSOO. It is structurally quite similar to prot-Cpd0; the main changes concern the breaking O-O bond (lengthened from 1.558 to 1.813 Å) and the strengthening FeO bond (shortened from 1.996 to 1.765 Å). The transition state for the conversion of prot-Cpd0 back to Cpd 0 has not been located precisely, but we know from energy scans (see above) that this process must also be facile, with an estimated barrier of 3-4 kcal/mol. The QM/MM results for a second snapshot with a different MM environment are documented in the Supporting Information (SI). The optimized geometry of Cpd 0 is quite similar in both snapshots (Figure S2), with a chain of three water molecules connecting Glu366 and Thr252. This hydrogen bond network changed, however, during the energy scans and optimizations that were performed to convert Cpd 0 to prot-Cpd0 in the case of the second snapshot (contrary to the first one). The optimized geometry of prot-Cpd0 (Figure S4) now has two hydrogen bonds between the Glu366 carboxylate oxygen atoms and the water molecules (rather than only one), and the shortest connecting chain now involves only two water molecules between Glu366 and Thr252. This reorientation of the hydrogen bond network is apparently stabilizing since the energy of pro-Cpd0 relative to Cpd 0 is now computed to be 23.1 kcal/mol (rather than 33.1 kcal/mol) with basis B1. The network topology remains unchanged in the optimized geometry of Cpd I for the second snapshot (Figure S6). Energy scans starting from prot-Cpd0 show that its conversion to both Cpd 0 and Cpd I is very facile also in the case of the second snapshot (again with barriers of less than 5 kcal/mol). For the sake of completeness, we note that we have studied (less extensively) a third snapshot that is structually analogous to the second one and has a shallow protCpd0 minimum which lies 24.6 kcal/mol above Cpd 0 with basis B1. A comparison between the present QM/MM results for the Glu366 channel and our previous QM results for a large 96atom model system11 is not straightforward. One major difference is that the QM/MM treatment incorporates the structural constraints imposed by the enzyme environment, whereas the published QM results11 refer to fully optimized geometries of the model system. For example, Asp251 remains essentially fixed in a salt bridge with Arg186 in all QM/MM optimized structures for the Glu366 channel, while it moves toward the active site in the QM model system to enhance hydrogen bonding and electrostatic attractions. As a consequence, the hydrogen bonding patterns are not identical in the QM and QM/MM structures, and targeted attempts to reproduce the QM patterns11 at the QM/ MM level were not successful because unconstrained QM/MM optimizations caused reorientation. Compared with the hydrogen bonding networks found in the QM model study,11 the second QM/MM snapshot (see above) shows the closest similarities, but there are still some differences. Focusing on prot-Cpd0 (see Figure S4 in SI and Figure 1b in ref 11), the H2O moiety in the FeOOH2 group has a different orientation, but more importantly it is much further away from the Asp251 carboxylate group in the QM/MM optimized structures, the relevant OO distances being 9.34/8.23 Å in the second QM/MM snapshot and 5.73/ 5.04 Å in the QM model structure.11 Since the net charges are approximately +0.5e on the H2O moiety11 (Figure 3) and -1.0e on the Asp251 carboxylate, this translates into a differential

QM/MM Study of Mechanisms for Compound I Formation

ARTICLES

Figure 3. QM/MM optimized geometries (UB3LYP/CHARMM) of Cpd 0, prot-Cpd0, TSOO, and Cpd I in mechanism I with basis B1 (distances in Å). The numerical values refer to basis B1 (in parentheses: basis B1a). In the case of prot-Cpd0, the net Mulliken charge Q and the spin density F of the H2O moiety in FeOOH2 are also given.

electrostatic stabilization of about 12 kcal/mol of prot-Cpd0 (QM vs QM/MM) as estimated from classical point-charge interactions. The energy of prot-Cpd0 relative to Cpd 0 is 4.5 kcal/ mol in the QM model system (LACVP* basis)11 and 23.1 kcal/ mol in the second QM/MM snapshot (B1 basis, see above). Much of this difference can be attributed to the differential electrostatic stabilization effect discussed above, and the remainder may arise from the difference in the hydrogen bonded networks. It should be noted in this context, for example, that the hydrogen bond between Asp251 and Wat901 is shortened significantly in the QM model system when going from Cpd 0 to prot-Cpd0 (see Figure 1 of ref 11), whereas there is no such effect at the QM/MM level. Hence, the stabilization of protCpd0 is less effective in the QM/MM than in the QM model calculations, mainly because of the structural constraints imposed by the protein environment. In summary, the present QM/MM study yields a local minimum for prot-Cpd0 only in the Glu366 channel, and not in the Asp251 channel. Whether prot-Cpd0 exists in a shallow potential well or not thus seems to depend on the protonation state of the enzyme environment. Even in the Glu366 channel, the energy of prot-Cpd0 relative to Cpd 0 is quite high (more than 20 or 30 kcal/mol for different hydrogen bond networks), and the barriers for decay to Cpd 0 and Cpd I are quite low (less than 5 kcal/mol). In the absence of more extensive sampling, we conclude that prot-Cpd0 is at best a reactive highenergy species. Mechanism I (see Scheme 1) thus does not provide a facile pathway for the conversion of Cpd 0 to Cpd I as long as the formation of prot-Cpd0-like species is involved. We cannot exclude the existence of a concerted protonation-

dissociation pathway from Cpd 0 to Cpd I that would bypass prot-Cpd0-like species, but we have not found any evidence for a low-energy mechanism of this kind in our extensive QM/ MM explorations. 3.3. Mechanism II: O-O Cleavage Followed by Proton Transfer with Concomitant Electron Transfer. This mechanism involves an initial O-O bond cleavage followed by a proton transfer with a concomitant electron transfer to yield Cpd I and water (Scheme 1). We have performed corresponding QM/ MM calculations for both the Glu366 and Asp251 channels. We discuss the latter (more extensive) results first (see Figures 4-6 and Table 3) Two reaction coordinates were adopted for O-O bond cleavage. The straightforward choice is the distance between two oxygen atoms in FeOOH. The alternative is the difference between this OO distance and the distance between the distal hydrogen and proximal oxygen atoms, which guides the system toward a hydrogen bond FeO‚‚‚HO. Both reaction coordinates led to the same intermediate complex (denoted as IC1), in which OH is a hydrogen bond donor to the FedO moiety and a hydrogen bond acceptor from the Thr252 hydroxyl group. The fully optimized structures of the intermediate complex (IC1) and the transition state for its formation (TS1) are shown in Figure 4. Going from Cpd 0 to TS1, the FeO bond shortens (1.894 f 1.680 Å), the OO bond is breaking (1.528 f 2.063 Å), and the hydrogen bond between Thr252 and the distal oxygen atom becomes stronger (1.843 f 1.640 Å) due to the displacement of the OH moiety toward Thr252. In IC1, the FeO bond remains short (1.680 Å), the OO bond is completely broken (2.628 Å), the hydrogen bond HO‚‚‚Thr252 is further strengthJ. AM. CHEM. SOC.

9

VOL. 128, NO. 40, 2006 13209

Zheng et al.

ARTICLES

Figure 4. Optimized geometries of Cpd 0, TS1, IC1, and Cpd I (UB3LYP/B1/CHARMM, DQ1 QM region). Only the QM region is shown. The dihedral angles OD2-CG-CB-CA and CG-CB-CA-C specify the orientation of the Asp251 side chain; the changes in the latter indicate reformation of the salt bridge with Arg186 after proton transfer. For Cpd I, the second set of values given refers to another local minimum (DQ1b in Table S4, see also Figure S17).

ened (1.602 Å), and there is a new strong hydrogen bond FeO‚‚‚HO (1.616 Å). These structural features are remarkably similar to those reported in DFT model studies19,22 of the analogous rearrangement in isolated Por(SH)FeOOH (see Figure S2 in ref 19 and Figure 6 in ref 22). The corresponding transition state has qualitatively the same shape as TS1, with a short FeO bond (1.685, 1.658 Å) and a long OO distance (2.329, 2.186 Å), and the intermediate has, like IC1, a short and essentially linear O‚‚‚HO bond (1.661, 1.662 Å).19,22 There are also obvious structural differences, of course, due to the presence of the protein environment in our QM/MM calculations. The most evident one concerns the four-atom FeO‚‚‚HO moiety in the intermediate which is linear in one of the reported gas-phase structures22 and bent toward Thr252 to form a hydrogen bond in the enzyme (Figure 4). 13210 J. AM. CHEM. SOC.

9

VOL. 128, NO. 40, 2006

The computed QM/MM barrier for the homolytic O-O cleavage (Cpd 0 f IC1) is 14.4 kcal/mol (UB3LYP/B1/ CHARMM, DQ1 QM region). Single-point QM/MM calculations with larger basis sets give slightly higher barriers (B2/ B3: 17.1/15.4 kcal/mol). The QM barriers reported for the gasphase model system Por(SH)FeOOH are significantly higher (21.6 kcal/mol,19 20.3-24.8 kcal/mol22). One obvious reason is that the gas-phase model calculations do not account for the transition state stabilization provided by the strengthened hydrogen bond between the distal oxygen atom and Thr252 (see above). The protein environment thus accelerates the initially homolytic O-O cleavage in Cpd 0. In the intermediate IC1, the FedO moiety carries two unpaired electrons, and the third unpaired electron is mainly located on the OH moiety. Antiferromagnetic and ferromagnetic

QM/MM Study of Mechanisms for Compound I Formation

ARTICLES

Figure 5. Spin density on OH in the intermediate IC1 (doublet state) from QM/MM calculations with different basis sets (B1, B2, and B3) and QM regions (DQ1, DQ3, and Por(SH)FeOOH). See text.

Figure 6. Energy profiles of mechanism II in the Asp251 channel. The solid line connects energies obtained with the DQ1 QM region and the B1 basis set; the given numerical values refer to the B1/B2/B3 basis sets (B1 optimized, B2 and B3 single-point energies). TS2 results estimated from energy scan. The dashed line and the values in italics are based on an energy scan with the DQ3 QM region and the B1 basis set (see text). Table 3. Computed Spin Densities and Mulliken Charge Distributions for Doublet States in the Asp251 Channel with the DQ1 QM Region Fe

OH/H2Oa

Op

SH

Por

Cpd 0 spin 1.042 -0.004 0.057 -0.030 -0.066 charge 0.492 -0.082 -0.364 -0.426 -0.542

Thr252

Asp251

0.000 0.000 0.025 -0.171

TS1

spin 1.352 -0.622 0.422 -0.052 -0.092 -0.008 0.000 charge 0.467 -0.158 -0.354 -0.449 -0.385 -0.012 -0.177

IC1

spin 1.433 -0.834 0.680 -0.055 -0.231 0.007 0.000 charge 0.471 -0.094 -0.524 -0.430 -0.286 -0.026 -0.180

Cpd I spin 1.319 -0.001 0.804 -0.186 -0.936 0.000 0.000 charge 0.464 -0.045 -0.436 -0.367 0.321 -0.052 -0.901 a

H2O refers to Cpd I, and OH, to the other species.

coupling leads to doublet and quartet states, respectively, which are predicted to be almost degenerate in IC1 (11.71 vs 11.70 kcal/mol relative to Cpd 0). Concerning the low-spin/high-spin

energy separation, IC1 thus resembles Cpd I (almost degenerate)23 rather than Cpd 0 (doublet/quartet gap of ca. 13-14 kcal/ mol, see above). Figure 5 shows the spin density on the OH moiety in the doublet state of IC1 (see also Table 3) calculated with three different basis sets (B1, B2, and B3) and QM regions (Por(SH)FeOOH, DQ1, DQ3). It is obvious that the absolute values of the spin density are less than 1 for the larger QM regions DQ1 and DQ3 indicating that OH will not behave as a “perfect” radical in IC1. This is due to the strong hydrogen-bonding interactions of the OH moiety with both the heme (through FeO) and the protein environment (through Thr252). It has already been argued22 that the former will diminish the radical reactivity of OH in the gas-phase model intermediate Por(SH)FeO‚‚‚HO and the latter should have a similar effect in the enzyme, on analogous grounds. For the minimal Por(SH)FeOOH QM J. AM. CHEM. SOC.

9

VOL. 128, NO. 40, 2006 13211

Zheng et al.

ARTICLES

Figure 7. Energy profile of mechanism II in the Glu366 channel. Energies were obtained with the EQ QM region; the given numerical values refer to the B1/B2/B3 basis sets (B1 optimized, B2 and B3 single-point energies).

region, the computed QM/MM spin density of OH is around -1.00e, in agreement with previous QM model results.19 This value is reduced to about -0.80e for the DQ1 QM region (inclusion of Thr252, Wat901, Asp251) and further to around -0.75e for the DQ3 QM region (additional inclusion of Arg186, see Table 1). This trend reflects the enhanced interaction of the OH moiety with the protein environment in the larger QM regions (along the Asp251 channel) which should be accompanied by a corresponding “dilution” of OH radical reactivity. The hydrogen-bonding interaction in the Asp251 channel (DQ1) and the electrostatic attraction between the positively charged Arg186 and the negatively charged active site (DQ3) exert a “pull” effect that facilitates electron transfer from the heme to OH and thus reduces the OH spin density. For the conversion of IC1 to Cpd I, a proton needs to be transported from Asp251 via Wat901 and Thr252 to OH, with a concomitant electron transfer from the heme, to form water and Cpd I. We find this proton transfer to be a concerted process through the hydrogen bond network Asp251-Wat901-Thr252OH. The barrier is only about 2.5 kcal/mol which is an estimate from an energy scan with the DQ1 QM region and the B1 basis set. The corresponding transition state could not be located precisely. A second energy scan was carried out using the larger DQ3 QM region, which includes the Arg186 side chain and is thus neutral. In this case, the reaction is barrierless (see Figure S1 and Table S6 in the SI). The first six steps on the DQ3 path involve a shortening of the hydrogen bond by 0.2 Å and a reduction of the OH spin density from -0.8 to -0.5 at almost constant energy (within 0.4 kcal/mol), and after this initial passage toward a heterolytic situation, the reaction smoothly proceeds downhill to the heterolytic products (note that the heterolytic character develops faster with the larger B3 basis; see Table S7 in the SI). A reduction of the barrier when going from DQ1 to DQ3 is plausible because of the enhanced interactions along the Asp251 channel due to the “pull” effect of Arg186 (see above). The overall energy landscape is shown in Figure 6. Single-point energies with the larger basis sets B2 and B3 are included. According to Figure 6, IC1 occupies a very shallow potential well and can easily decay to Cpd I and 13212 J. AM. CHEM. SOC.

9

VOL. 128, NO. 40, 2006

water, possibly even in a barrierless fashion. This process should be very fast and give the OH moiety little chance for radical reactivity. If the second step in Figure 6 is indeed barrierless, the overall conversion of Cpd 0 to Cpd I can be viewed as one concerted process, i.e., an initial homolytic O-O bond cleavage triggering a concomitant proton and electron transfer. In this case, there is no real “intermediate” and hence no significant side reactions. After the proton transfer in the Asp251 channel is completed, the side chain of Asp251 rotates back to form a salt bridge with Arg186 (see Figure 4). The mechanistic view that emerges from the present QM/ MM calculations is thus qualitatively different from the somersault mechanism for P450 hydroxylation of hydrocarbons that has recently been proposed on the basis of QM model calculations.22 In both cases, the conversion of Cpd 0 starts with an initially homolytic O-O cleavage, although the QM/MM study indicates some hybrid homolytic-heterolytic character already at an early stage of the reaction. In the postulated somersault mechanism, the formed “hydrogen-bonded hydroxyl radical” abstracts hydrogen from the substrate, and subsequent rearrangements then yield the products without any involvement of Cpd I.22 Contrary to the QM model study,22 the present QM/ MM calculations include the environment, and the initially formed species can thus accept a proton through the Asp251 channel, together with an electron from the heme, in an essentially barrierless process. Hence, the protein environment opens a facile path to the heterolytic formation of Cpd I which is favored over other conceivable channels.22 In this scenario, Cpd I retains its role as the decisive oxidant. Up to this point, we have only discussed the Asp251 channel in this section. We now briefly summarize our QM/MM results for mechanism II operating in the Glu366 channel (see Figures 7 and 8). For the initial homolytic O-O bond cleavage, a barrier of 13.0 kcal/mol is found with the EQ QM region and the B1 basis set. Single-point calculations with larger basis sets give similar values (B2/B3: 15.8/14.3 kcal/mol). The OH moiety in the resulting intermediate IC1 can simultaneously get a proton from Thr252 and an electron from the heme to form the a second intermediate (IC2) containing Thr252-O-, water, Cpd I, and a

QM/MM Study of Mechanisms for Compound I Formation

ARTICLES

Figure 8. Optimized geometries of Cpd 0, TS1, IC1, TS2, IC2, TS3, and Cpd I (UB3LYP/B1/CHARMM) in the Glu366 channel for mechanism II. Only the QM region is shown.

water chain. The corresponding barrier (TS2) is only 1.1 kcal/ mol (EQ QM region, B1 basis set) and even lower in singlepoint calculations with larger basis sets (B2/B3: -0.7/0.9 kcal/

mol) implying that the formation of Cpd I from IC1 is essentially barrierless in the Glu366 channel. To complete the reaction, a proton must still be transported from Glu366 to Thr252-O- via J. AM. CHEM. SOC.

9

VOL. 128, NO. 40, 2006 13213

Zheng et al.

ARTICLES

the chain Glu366-Wat523-Wat566-Wat902-Thr252. This proton transfer is concerted with a barrier (TS3) of 6.6 kcal/ mol (B1 basis) which becomes significantly higher in singlepoint calculations with larger basis sets (B2/B3: 12.9/12.5 kcal/ mol). Structurally, the initial transition state TS1 for O-O cleavage and the resulting intermediate IC1 are fairly similar in the Asp251 and Glu366 channels (see Figures 4 and 8). The concerted nature of the proton transfer along the Glu366 water chain is reflected in the evenly distributed hydrogen bond lengths in TS3 (Figure 8). Comparing mechanism II in the Asp251 and Glu366 channels, the first stage is essentially the same, i.e., O-O cleavage with a rate-limiting barrier of 13-14 kcal/mol for basis B1 (slightly higher with larger basis sets). The second stage is however quite different, since the concomitant proton and electron transfer is essentially barrierless in the Asp251 channel, but encounters a sizable barrier in the Glu366 channel. The latter should thus play no significant role. The Asp251 channel for formation of Cpd I also has the conceptual advantage that it is connected to the surface of the protein and can thus deliver protons easily, whereas an additional separate proton-transfer process from the surface would be required in the Glu366 channel to reprotonate Glu366 in its deeply buried hydrophobic environment. The cost of protonating Asp251 may be estimated from the classical MD simulations described in section 3.1, which predict Asp251 side chain rotation to be facile, with a free energy cost of less than 4 kcal/mol to go from the salt-bridge conformation with deprotonated Asp251 to the twisted conformation characteristic of protonated Asp251 (for both 1AspP and 2Asp, see section 3.1). At the optimized QM/MM geometry of the latter conformation of Cpd 0 (see Figure 4), the empirical PROPKA procedure39 yields a pKa value of 6.3 for Asp251 indicating an equilibrium protonation of almost 20% at pH 7. At the same geometry, the PROPKA pKa estimate is 10.5 for Glu366 implying an essentially complete equilibrium protonation at pH 7 for this residue in its hydrophobic environment. Thermodynamically, it is thus more favorable to protonate Cpd 0 at Glu366 than at Asp251. However, there is no proton source close to Glu366, and proton transfer from the surface of the protein to Glu366 via a hydrogen-bonding network will probably require a nonnegligible activation and may well proceed with involvement of Asp251 as a proton shuttle residue. It is thus unlikely that there is any kinetic advantage of protonating Glu366 rather than Asp251 under turnover conditions. Accounting for the cost of protonating Asp251 and Glu366 should therefore not change the preference for the Asp-M-II mechanism (compared with GluM-II), but will increase the effective overall barrier for the formation of Cpd I in the Asp-M-II mechanism (by around 4 kcal/mol). Comparing with other systems, it is known that O-O bond homolysis is not too sensitive to the proximal ligand at iron since no oxidation of the ligand is required.17b It is therefore not surprising that DFT model calculations on Por(L)FeOOH give similar barriers and transition states for L ) imidazole (representing histidine in heme oxygenase)20 and L ) SH (representing cysteine in P450).18,19,22 The nature of the proximal ligand becomes important, however, in the subsequent electron transfer from the heme toward OH that occurs in the proposed mechanism II concomitant with proton transfer. Relative to the histidine ligand in heme oxygenase, the cysteine ligand in P450 (39) Li, H.; Robertson, A. D.; Jensen, J. H. Proteins 2005, 61, 704-721. 13214 J. AM. CHEM. SOC.

9

VOL. 128, NO. 40, 2006

exerts a larger “push” effect and will thus facilitate such an electron transfer which tends to diminish the reactivity of the bound OH radical and to guide the system toward the “heterolytic” products, Cpd I and water. The “pull” effect in the Asp251 channel (see above) will reinforce this tendency and lead to an essentially concerted conversion which is initially of homolytic and later of heterolytic character. It is conceivable that radical side reactions may be avoided in this manner. 3.4. Cpd I Relative Energy. As discussed in section 1, there is no consensus in the literature concerning the energy of Cpd I relative to Cpd 0. Here we summarize our current QM/MM results obtained from UB3LYP/CHARMM geometry optimizations with the B1 basis and from subsequent single-point calculations with the larger B2 and B3 basis sets. In addition to these QM/MM energies, we report QM-ee and QM energies at the optimized QM/MM geometries obtained from UB3LYP calculations with and without inclusion of the MM charges, respectively. The results for the Glu366 channel are quite uniform (see Tables S1 and S2 in the SI). In this case, the salt bridge between Asp251 and Arg186 is present in Cpd 0 and Cpd I since Glu366 serves as proton source. Cpd I is found to be more stable than Cpd 0 by 2-4 (1-2) kcal/mol when studying mechanism I (II) with basis sets B1-B3; the QM-ee results differ from the QM/MM energies by less than 1 kcal/mol in each case. The results for the Asp251 channel fluctuate more strongly. Here the salt bridge between Asp251 and Arg186 is present only in Cpd I since Asp251 is protonated in Cpd 0, and it is therefore not surprising that a balanced treatment of the corresponding interactions is required to obtain reliable results. Figure 9 shows the computed energies of Cpd I relative to Cpd 0 for various QM regions, basis sets, and conformations (see the SI for numerical results). According to the QM/MM optimizations with the B1 basis and the DQ1-DQ4 QM regions, Cpd I is usually slightly more stable than Cpd 0 (in six of seven cases, by 0.2-4.4 kcal/mol), except for DQ2 where it is less stable by 2.3 kcal/mol. The QM/MM and QM-ee results are consistent with each other (deviation e 1 kcal/mol) only for the DQ3 and DQ4 QM regions which include Asp251 as well as Arg186 and thus treat both partners in the salt bridge on equal footing (deviations around 10 kcal/mol for DQ1 and DQ2). The pure QM results without external MM charges are unbalanced and unreliable in this case: they predict Cpd I to be less stable than Cpd 0 by ca. 10 kcal/ mol for DQ1 and DQ2 (consistent with analogous previous work11) and more stable by ca. 25 kcal/mol for DQ3 and DQ4, the extra stabilization coming from the electrostatic attraction in the salt bridge. Going from the B1 to the larger B2 or B3 basis sets favors Cpd I relative to Cpd 0 by typically 5-8 kcal/ mol (see Tables S3 and S4 in the SI). Our best estimate from the QM/MM calculations for the Asp251 channel is that Cpd I is almost 10 kcal/mol more stable than Cpd 0 (QM region DQ3 with the large B3 basis). In the Glu366 channel, Cpd I is favored by a smaller margin (2-4 kcal/ mol, B3 basis). Our present QM/MM results demonstrate that the relative energy of Cpd 0 and Cpd I depends on the protein environment, in particular on the salt bridge between Asp251 and Arg186. Realistic predictions require a careful choice of the QM region and a balanced treatment of all relevant electrostatic interactions. In an overall assessment, Cpd 0 and Cpd I appear to be of similar energies, with a slight preference for Cpd I in the currently studied systems.

QM/MM Study of Mechanisms for Compound I Formation

ARTICLES

Figure 9. Energy of Cpd I relative to Cpd 0 for various QM regions (DQ1-DQ4) and conformations in the Asp251 channel. Labels a, b, and c denote different conformations in a given QM region. QM-ee and QM-gas refer to QM energies computed with and without MM point charges, respectively. B1 data unless noted otherwise; B2 and B3 data are marked by an outer square.

4. Conclusions

Two mechanisms for the conversion of Cpd 0 to Cpd I in P450cam have been investigated at the QM/MM level. The required proton comes either from Glu366 or from Asp251. Classical free energy simulations for the torsion of the Asp251 side chain confirm that the initial conformations in the Asp251 channel can be reached with small energetic cost (ca. 4 kcal/ mol). Mechanism I assumes an initial protonation of the distal oxygen atom in Cpd 0. A local minimum is found for protCpd0 only in the Glu366 channel, but it is very high in energy (more than 20 kcal/mol above Cpd 0) and the barriers for its decay are very low (only 3-4 kcal/mol both toward Cpd 0 and Cpd I). A corresponding prot-Cpd0 minimum could not be found in the Asp251 channel, and all of the pathways considered for mechanism I required large activation. In mechanism II, an initial O-O bond cleavage followed by a concomitant proton and electron transfer yields Cpd I and water. The rate-limiting step in both channels is O-O bond cleavage with a barrier of about 13-14 kcal/mol. This generates a hydrogen-bonded species containing a bound OH moiety with diminished radical reactivity. In the Asp251 channel, the “push” effect of a thiolate ligand and the “pull” effect of Arg186 combine to facilitate a concomitant electron and proton transfer that yields Cpd I and water in an essentially barrierless reaction. The corresponding conversion encounters a sizable barrier in the Glu366 channel. Hence, according to the current QM/MM calculations, mechanism II in the Asp251 channel provides the preferred pathway for the conversion of Cpd 0 to Cpd I. Electronically, this reaction exhibits a hybrid character since it is homolytic initially and becomes heterolytic soon after the first transition state. Accounting for the cost of protonation is not expected to change the preference for the Asp251 channel but will increase the effective overall barrier for this pathway. It is however possible that in the P450cam mutants which involve Asp251, e.g., D251N, the Glu366 channel would operate. Indeed, the D251N mutant is known to retain its activity toward camphor but to react at a significantly slower rate compared with the wild-type enzyme.1c

The relative energy of Cpd 0 and Cpd I depends on the protein environment and is tuned by the electrostatic interactions around the active site. In the systems studied presently, Cpd I tends to be slightly more stable according to the QM/MM results (by a few kcal/mol). We end with a cautionary remark: While the QM/MM calculations are more robust than QM model calculations since they account for the structural constraints and the influence of the enzyme environment, they are still sensitive to the chosen setup (protonation states, number of water molecules, etc.). Moreover, the local geometry optimizations performed presently sample only a limited number of possible conformations and hydrogen bonding networks. Given these limitations, there is a clear need for further, more comprehensive studies which should also address two key questions raised by the present work: (a) If the traditional protonation mechanism (I) for the formation of Cpd I does not operate, can the proposed hybrid homolyticheterolytic mechanism (II) account by itself for the experimental fact17,18 that homolytic and heterolytic pathways coexist and can be distinguished through the observation of different products when using probe substrates such as cumene hydroperoxide? And (b) Does the proposed mechanism (II) provide a consistent overall picture of P450 reactivity, especially with regard to the energy profile for Cpd I formation and subsequent camphor hydroxylation by Cpd I? Since at present we cannot answer these questions without some speculative assumptions, we must defer them to future studies. Acknowledgment. The research at HU was supported by a grant from BMBF (DIP-G.7.1). Supporting Information Available: Computational details. Relative energies, energy scans, spin densities, charge distributions, and optimized structures with selected geometrical parameters (8 tables and 23 figures). Complete refs 27 and 35. This material is available free of charge via the Internet at http://pubs.acs.org. JA063439L J. AM. CHEM. SOC.

9

VOL. 128, NO. 40, 2006 13215