J. Phys. Chem. B 2004, 108, 10083-10088
10083
QM/MM Study of the Product-Enzyme Complex in P450cam Catalysis Hai Lin,† Jan C. Scho1 neboom,† Shimrit Cohen,‡ 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: February 12, 2004; In Final Form: April 29, 2004
The enzyme-product complex in P450cam (CYP101) has been studied by combined quantum mechanical/ molecular mechanical (QM/MM) calculations. The central iron(III) porphyrin complex and part of the catalytic product (5-exo-hydroxycamphor) are treated with density functional theory, while the protein/solvent environment is represented by the CHARMM force field. The computations indicate a doublet minimum at an Fe-O distance of ca. 2.2 Å, and a flat, barrierless potential for the dissociation of the Fe-O bond. Comparisons with analogous calculations on the isolated QM system in the gas phase show that inclusion of the protein/solvent environment lowers the activation energy for bond dissociation in the doublet state because of interactions within the binding pocket and accounts for a significant stabilization of the quartet and sextet states. The theoretical results allow for a tentative interpretation of recent ENDOR data (Davydov, R.; Makris, T. M.; Kofman, V.; Werst, D. E.; Sligar, S. G.; Hoffman, B. M. J. Am. Chem. Soc. 2001, 123, 1403).
1. Introduction
SCHEME 1
Cytochrome P450 enzymes are ubiquitous hemo-proteins that are involved in the metabolism of very diverse compounds, such as steroids, fatty acids, and alkaloids. These monooxygenases play an important role in many bioregulatory functions, for example, biosynthesis of hormones or detoxification of xenobiotics.1 Most experimental data are available for the bacterial enzyme P450cam (CYP101) from Pseudomonas putida, including X-ray structures for various intermediates in the catalytic cycle. This enzyme catalyzes the highly stereoselective hydroxylation of camphor, yielding 5-exo-hydroxycamphor (see Scheme 1). X-ray,2 UV-vis,2 and EPR3,4 data for the product complex 1 show that 5-exo-hydroxycamphor is coordinated to the heme iron atom and indicate that a mixture of low-spin (LS) doublet (2A) and high-spin (HS) sextet (6A) species is present, with the LS species being predominant. Scheme 2 shows the corresponding orbital occupations together with that of the alternative quartet (4A) state. According to detailed low-temperature electron nuclear double resonance (ENDOR) studies at 200 K,4 product formation occurs through three distinct hydroxycamphor-bound conformations. The spectroscopic results suggest that the immediate product after hydroxylation corresponds to a nonequilibrium state, involving a more or less normal Fe-O distance (≈2 Å), and that the heme pocket then relaxes in two detectable steps to accommodate the anomalously long Fe-O bond of 2.67 Å, found in the crystal structure.2 The long Fe-O bond in the equilibrium state is thought to be a compromise between favorable Fe-O interactions on one hand and unfavorable desolvation of the polar hydroxy group in the nonpolar protein pocket on the other hand. Recent QM/MM calculations at the B3LYP/CHARMM level considered the 2A and 4A states of 1.5 This study pointed to a * To whom correspondence should be addressed. E-mail:
[email protected]; phone: +49-208-306-2150; fax: +49-208-306-2996. † Max-Planck-Institut fu ¨ r Kohlenforschung. ‡ The Hebrew University.
SCHEME 2
quartet ground state, showing a single minimum with a broken Fe-O bond, that is, r(Fe-O) ) 2.84 Å. The doublet state exhibited two nearly isoenergetic minima approximately 5 kcal/ mol higher in energy, one with a long (2.82 Å) and one with a short (2.26 Å) Fe-O distance. However, because of the high conformational complexity of the system, the presence of another lower-lying doublet minimum could not strictly be ruled
10.1021/jp0493632 CCC: $27.50 © 2004 American Chemical Society Published on Web 06/16/2004
10084 J. Phys. Chem. B, Vol. 108, No. 28, 2004 out. In the present paper, we report potential energy surface (PES) scans of the internal coordinate r(Fe-O) at the QM/MM level, carried out for the 2A, 4A, and 6A states, respectively. This allows us to identify the global minima and the crossing points of the spin surfaces. Apart from checking for the existence of previously unidentified minima of 1, the present study of the dissociation of the Fe-O bond is of interest in the context of product release in P450cam catalysis. The breaking of the Fe-O bond constitutes the first step, before the product eventually is expelled from the binding pocket. Comparative calculations on the isolated QM region in the gas phase were carried out for complex 1 to assess the influence of the protein/ solvent environment on this process. Our results show that the energy cost associated with dissociation of the Fe-O bond is significantly reduced by the protein environment, a feature which is crucial for catalytic efficiency. 2. Computational Details A. Setup of the System. Initial coordinates were derived from QM/MM (B3LYP/CHARMM) optimized structures of the 2A state of Compound I (snapshot 40) reported in ref 6. The entire system consists of 24394 atoms, including 16956 atoms of solvent water. The actual starting geometry for the PES scans was obtained by replacing the camphor atoms and the oxo ligand in the QM/MM optimized structure of Compound I by the 5-exo-hydroxycamphor atoms from the X-ray structure2 (PDB code 1NOO, Brookhaven Protein Database) after a superimposition of these two structures on the basis of a least-squares fit of the heme group. The resulting geometry of the complex is close to the X-ray conformation2 and shows the same prominent protein-substrate interactions. Most notably, the camphor carbonyl forms a H-bond with the hydroxyl of Tyr96, and camphor methyl groups have hydrophobic interactions with residues Phe87, Val295, Ile395, and Val396. B. QM/MM Calculations. The general QM/MM methodology adopted in the present study is documented in refs 6 and 7. In the following, we only describe the aspects relevant to the present work. QM Methods. The QM part was treated with unrestricted Kohn-Sham DFT using the B3LYP density functional.8 As has been pointed out recently for FeII transition-metal complexes,9,10 the B3 exchange functional overestimates the stability of HS states relative to LS configurations because of the admixture of exact (Hartree-Fock) exchange, whereas pure GGA functionals are strongly biased toward LS situations. It has therefore been recommended9 to reduce the amount of 20% exact exchange in the B3 expression to 15%. Instead of choosing this pragmatic approach, we have carried out additional calculations with the BLYP functional (no exact exchange). This may serve to assess the intrinsic accuracy of the DFT approach, and to estimate “upper bounds” for the calculated LS-HS splittings. The QM region comprised iron-porphyrin (without side chains of the heme), sulfur of Cys357, C4H, C5H(OH), and C6H2 of 5-exo-hydroxycamphor. In the gas-phase calculations, this corresponds to [Fe(SH)(porph)] + (CH3)2CHOH (51 QM atoms). The basis set consists of a small-core effective core potential and the associated LACVP basis of a double ζ quality for iron11 and a 6-31G basis for all other atoms.12 This combination of QM region and basis set was shown to reproduce the electronic and geometric features of a more extended QM region (including full camphor and cysteine357) and basis generally well.5,6 QM/MM Methods. An electronic embedding scheme13 was applied, that is, the fixed MM charges were included into the
Lin et al. one-electron Hamiltonian of the QM calculation and the QM/ MM electrostatic interactions were evaluated from the QM electrostatic potential and the MM atomic charges. No cutoffs were introduced for the nonbonding MM and QM/MM interactions. To treat the QM/MM boundary, we used hydrogen link atoms with the charge shift model.14,15 Codes. For the QM treatment in the QM/MM as well as in the pure QM calculations, we employed the TURBOMOLE program.16 All QM/MM calculations were performed with the ChemShell package.7,17 The CHARMM22 force field18,19 run through the DL-POLY program20 was used for the treatment of the MM part of the system. Gas-Phase Calculations. To assess the influence of the protein environment on the QM subsystem, full geometry optimizations of the isolated species in the gas phase were performed. Equilibrium geometries were obtained by unconstrained minimizations of the isolated QM system. Potential energy surfaces for r(Fe-O) were generated by a stepwise scan along this internal coordinate. Geometry Optimizations. The QM/MM optimizations included the following set of residues close to the active site: Pro86, Phe87, Tyr96, Pro100, Thr101, Leu244, Gly248, Asp251, Thr252, Val295, Asp297, Arg299, Phe350, His355, Leu356, Cys357, Leu358, Gly359, Gln360, Ile395, Val396, heme, 5-exohydroxycamphor, Wat61, Wat253, and Wat325. This corresponds to a total of 442 optimized atoms. Geometry optimizations were performed with the HDLC optimizer,21 which is part of ChemShell. 3. Results A. Equilibrium Species. Unconstrained geometry optimizations from the points of lowest energy along the PES scans (see below) lead to the global minima of the three spin state surfaces. Figure 1 summarizes relative energies, selected structural features, and spin densities of the three species. For the full system in the enzyme environment, B3LYP/CHARMM optimizations predict the quartet minimum to have the lowest energy, being more stable than the doublet (sextet) by 4.0 (4.1) kcal/mol. The doublet geometry features a relatively short Fe-O bond of 2.249 Å, compared with experimental datum of 2.67 Å from X-ray data.2 The calculated distance for the quartet species is 2.686 Å, while the bond is essentially broken in the sextet minimum (2.832 Å). The elongation results from occupation of the d(z2) orbital in the two latter states, which makes the dominant contribution to the antibonding σ*(Fe-O) orbital. The calculated spin densities are indicative of the electronic configuration with one (doublet), three (quartet), or five (sextet) unpaired electrons occupying the d-orbitals at the metal; because of the mixing with ligand orbitals, the actual atomic spin density at iron is, however, reduced to 2.69 e in the quartet and 4.02 e in the sextet state. Turning to B3LYP calculations on the isolated QM systems in the gas phase, the order of spin states is reversed: Consistent with previous gas-phase calculations,22 the doublet minimum is calculated to lie 3.2 (5.2) kcal/mol below the quartet (sextet). Comparison of the gas-phase structures with the QM/MM results for the full system shows that the bonds to the axial ligands (Fe-O and Fe-S) are slightly shorter in all of the three minima (Fe-O bond lengths 2A: 2.226 Å, 4A: 2.620 Å, 6A: 2.712 Å). A typical feature upon going to the gas phase is an increase of positive spin density on the sulfur, indicating that the protein environment, by polarization and the presence of H-bond donors in the proximal pocket, effects a stabilization of electron density (decrease of unpaired spin density) on the sulfur.6 The change
QM/MM Study on Product Release from P450cam
J. Phys. Chem. B, Vol. 108, No. 28, 2004 10085
Figure 1. Equilibrium structures obtained by free minimization in the doublet, quartet, and sextet states in the enzyme environment (B3LYP/ CHARMM). Values in parentheses refer to B3LYP optimizations in the gas phase. (a) Selected structural data: distances in Å, angles (oblique numbers) in degrees. (b) Mulliken spin densities in e, B3LYP/CHARMM (B3LYP) relative energies in kcal/mol.
TABLE 1: B3LYP/CHARMM PES Scan of the Coordinate r(Fe-O) in the 2A, 4A, and 6A States in the Enzyme Environmentf EQM/MM r(Fe-O) 2.084 2.155a 2.235 2.249b 2.385 2.445c 2.535 2.685 2.686d 2.832e 2.835 2.984 3.134
2
A
4
A
0.42 0.12 0.06 0.00 0.48 0.72 1.57 1.90
1.90 0.04 -1.41
2.01 2.05 2.45
-3.89 -3.83 -3.58
-3.11 -3.74 -3.97 -3.97
TABLE 2: BLYP/CHARMM PES Scan of the Coordinate r(Fe-O) in the 2A, 4A, and 6A States in the Enzyme Environmenta
EQM 6
A
2
A
4
A
-2.24 -3.90 -4.37
0.98 0.70 0.84 0.54
-0.88 -0.66 -0.49 0.00 0.77 1.24 2.60 3.26
0.08 0.39 0.31 0.22
3.90 4.17 4.44
-4.98 -4.79 -4.30
5.45 2.38
-5.66 -5.48 -5.31 -5.32
EQM/MM 6
A
r(Fe-O)
3.96
1.935 2.084 2.235 2.385 2.535 2.685 2.835 2.984 3.134
1.97 1.11 1.06 1.68 1.86 1.54 1.81 1.97 2.24
a
Close to the crossing point of 2A and 4A QM/MM energy surfaces. b Minimum from an unconstrained optimization in the 2A state. c Close to the crossing point of 2A and 6A QM/MM energy surfaces. d Minimum from an unconstrained optimization in the 4A state. e Minimum from an unconstrained optimization in the 6A state. f Distance in Å; QM/ MM energy E(QM/MM) and QM contributions E(QM) in kcal/mol.
of ground-state multiplicity from doublet to quartet upon inclusion of the protein environment at the B3LYP/MM level is consistent with the findings of our previous QM/MM work.5 B. Potential Energy Surface Scans. Tables 1 and 2 summarize the relative QM/MM energies EQM/MM, and the corresponding QM contributions EQM, obtained from the scans in the enzyme environment with the B3LYP and BLYP functional, respectively. The potential energy profiles for the gas-phase scans (B3LYP and BLYP) are given in Table 3. In the B3LYP calculations, the reference point is the energy of the equilibrium structure in the 2A state. BLYP energies were obtained by reoptimizing starting from the B3LYP geometries; the reference point for relative BLYP energies is the point of lowest energy from the 2A scan. Plots of the energy profiles are shown in Figure 2. On the basis of the B3LYP energy profiles, the Fe-O distances at the crossing points (see Table 1) were estimated and constrained geometry optimizations at these distances were carried out to confirm the near-degeneracy of the spin states.
2
A
4.42 0.91 0.00 0.09 0.35 0.92 1.07 0.98 1.23
4
A
15.62 9.26 5.17 3.04 2.25 1.62 1.08 0.89 0.98
EQM 6
A
27.19 21.59 17.86 15.87 14.91 14.60 14.04 13.65 13.50
2
A
2.32 0.45 0.00 0.38 1.71 3.15 2.95 3.13 3.74
4
A
11.37 6.37 2.96 1.40 1.71 1.47 1.25 1.33 1.19
6
A
25.68 21.03 17.71 16.41 15.80 16.14 15.86 15.76 15.47
a Distance in Å; QM/MM energy E(QM/MM) and QM contributions E(QM) in kcal/mol.
TABLE 3: QM PES Scan of the Coordinate r(Fe-O) in the 2A, 4A, and 6A States in the Gas Phasea B3LYP r(Fe-O) 1.935 2.085 2.235 2.385 2.535 2.685 2.835 2.984 3.134 a
2
A
1.47 0.00 1.61 2.41 3.56 4.45 5.32 5.82 6.33
4
A
15.13 8.83 6.00 5.90 5.61 5.27 5.34 5.41 5.48
BLYP 6
A
15.22 9.46 8.54 7.55 7.57 7.22 7.12 7.07 7.00
2
A
2.16 0.00 0.71 1.19 1.91 2.78 3.37 3.90 4.36
4
A
22.92 14.80 12.84 11.04 10.10 9.68 9.43 9.29 9.15
6
A
33.31 25.34 23.64 22.13 21.31 20.83 20.50 20.22 19.98
Distance in Å; QM energy EQM in kcal/mol.
According to B3LYP/CHARMM PES scan in the enzyme environment, the quartet is the ground state for r(Fe-O) > 2.16 Å and is clearly preferred over the doublet (by 4-6 kcal/mol) for r(Fe-O) > 2.4 Å. The crossing of doublet and quartet surfaces occurs at a distance of r(Fe-O) of 2.16 Å, while the crossing of the 2A and 6A surfaces occurs at r(Fe-O) ) 2.45 Å. The potential energy minimum of the 2A state (see Table 1) is close to the crossing point of doublet and quartet surfaces, both in energetic (∆E < ≈ 0.2 kcal/mol) and geometric (∆r(FeO) ≈ 0.06 Å) terms. This might increase the probability of intersystem crossing between the two spin states, by placing the intersection of the spin surfaces into well-accessible regions
10086 J. Phys. Chem. B, Vol. 108, No. 28, 2004
Lin et al.
Figure 2. PES scans of the coordinate r(Fe-O) in the enzyme environment (QM/MM) and in the gas phase (QM). Energy in kcal/mol, bond length r in Å. 2A state: solid line, squares; 4A state: broken line, stars; 6A state: dashed line, triangles.
of the PES. As expected, QM/MM calculations employing the BLYP functional favor the doublet state, which becomes the ground state for r(Fe-O) < 2.8 Å. No crossing of doublet and sextet surfaces is observed in the range of considered Fe-O distances. Comparing the lowest BLYP/CHARMM energies of the three spin states obtained from the respective scans (see Table 2) indicates a global doublet minimum, with the lowest quartet and sextet species lying 0.9 and 13.5 kcal/mol higher in energy. With both functionals, B3LYP and BLYP, the QM/MM energy of the doublet state increases by only 1-2 kcal/mol when going from the minimum to r(Fe-O) ) 3.134 Å. Hence, the dissociation of the hydroxylated substrate on the doublet surface is facile. Analyzing the QM contributions to the QM/MM energy (Tables 1 and 2) shows that these terms alone account for a dissociation energy of ca. 4 kcal/mol, which we mainly attribute to the loss of Fe-O binding energy. However, the MM terms that reflect the interactions of the substrate with the protein pocket partly compensate for this destabilizing effect. In both gas-phase PES scans (B3LYP and BLYP), the quartet and sextet states are calculated at significantly higher energy as compared to the respective QM/MM points. Hence, B3LYP predicts a doublet (quartet) ground state for r(Fe-O) smaller (larger) than 2.9 Å, while no crossing at all is observed when the BLYP functional is employed. The sextet state in both cases has the highest energy over the whole range of distances r(FeO) considered (1.9-3.1 Å). This clearly shows that the inclusion of the protein environment stabilizes the quartet and sextet states relative to the doublet (see also Figure 2). For example, at r(FeO) ) 2.235 Å, the quartet (sextet) is lowered relative to the gas-phase calculations by 5.9 (4.6) kcal/mol and by 7.0 (5.1) kcal/mol in the B3LYP and BLYP calculations, respectively. Another aspect of the gas-phase study is the significantly larger dissociation energy on the doublet surface, which amounts
to 5.8 (3.9) kcal/mol in the B3LYP (BLYP) case compared to 2.5 (1.2) kcal/mol in the B3LYP/CHARMM (BLYP/CHARMM) calculations. 4. Discussion The present calculations complement the preliminary conclusions regarding the 5-exo-hydroxycamphor-P450cam complex drawn in our recent B3LYP/CHARMM study,5 which considered only local minima from optimizations in the 2A and 4A states. Two nearly isoenergetic doublet minima were found in that study,5 at r(Fe-O) ) 2.26 and 2.82 Å, respectively, while a single minimum (ca. 5 kcal/mol lower in energy) was located on the quartet surface at r(Fe-O) ) 2.84 Å. The PES scans reported herein show that the doublet species at short and long Fe-O distances are interlinked by a flat barrierless potential, permitting facile dissociation of the Fe-O bond. Furthermore, on the basis of the scan carried out at the B3LYP/CHARMM level, we can rule out the presence of another doublet minimum with a significantly lower energy than the quartet species at this level of theory. When comparing our earlier5 and present results in more detail, the protein environment was obtained from reaction path calculations with successive optimizations previously5 and derived from an X-ray structure currently (see section 2A). The local minima computed in these two studies will thus differ slightly because of minor differences in the protein environment. Given these circumstances, there is good agreement between the optimized Fe-O distances in the first doublet minimum (2.25 Å vs 2.26 Å). The larger deviation for the quartet minimum (2.84 Å vs 2.69 Å) reflects the flatness of the potential in this region (see Figure 2). The doublet potential is also extremely flat in this range, and it has indeed been possible previously to find a second local doublet minimum at
QM/MM Study on Product Release from P450cam r(Fe-O) ) 2.82 Å by an optimization starting from the quartet structure5 (which was not attempted presently). The computed QM/MM energy profiles for Fe-O dissociation in the protein environment show that the doublet state is the ground state at short bond distances. The relative stability of the spin state minima and the location of the crossing points of the spin surfaces depend on the amount of exact exchange employed in the density functional: The global minimum is predicted to be of the 4A-type with B3LYP, while BLYP gives it as 2A. Linear interpolation between the B3LYP/CHARMM and BLYP/CHARMM results indicates that the 4A minimum will remain below 2A when reducing the amount of exact exchange from 20% to 15% as recommended,9 but only by a small margin. Hence, definite conclusions regarding the multiplicity of the ground state of 1, and the crossing point of the doublet and higher spin surfaces, are difficult to reach at the DFT level. This issue requires further studies on the basis of ab initio methods, which can offer a more balanced description of exchange and correlation effects. It should be stressed, however, that these variations do not affect the overall mechanistic picture of camphor hydroxylation (in terms of two-state reactivity), since the product release is not the rate-limiting step of the reaction.5 In the mechanistic picture that emerges from the present and preceding studies,5 the initial product alcohol complex is formed in the 2A state after a nearly barrierless oxygen-rebound step5,22 and is characterized by a relatively short Fe-O bond (≈2.2 Å). This conformation is, however, not a thermodynamic sink: Dissociation of the Fe-O bond on the 2A surface requires very little activation. The loss of Fe-O binding energy is (partly) compensated by a stabilization of approximately the same magnitude because of interactions with the protein environment (apparent in the calculations from the MM terms), which could be a result of more favorable interactions of 5-exo-hydroxycamphor with the pocket. Candidates for these energy contributions are dispersive interactions of the camphor methyl groups with residues Phe87, Val295, Ile395, and Val396 and the hydrogen bond of the camphor carbonyl group with Tyr96. In fact, the distance r(O-O) between the camphor carbonyl oxygen and the hydroxo oxygen of Tyr96 decreases as r(Fe-O) is increased (data not shown), and single-point B3LYP calculations of a hydroxycamphor-phenol model complex (phenol representing Tyr96) indicate H-bond stabilization energies of 6.1 (6.6) kcal/mol at the optimized B3LYP/CHARMM geometries of the doublet (quartet) state. Eventually, the system can undergo intersystem crossing to the higher spin state, which represents the ground state of pentacoordinate FeIII systems. It is an intriguing fact that the quartet and sextet minimum structures are stabilized relative to the doublet in the enzyme environment (the quartet even becoming the ground state in the B3LYP/CHARMM calculations). This feature could facilitate intersystem crossing by shifting the intersections of the spin surfaces to well-accessible regions of the doublet PES. The underlying factors that cause this stabilization have been studied in the closely related ferric aquo heme complex, that is, the resting form of the enzyme. These results are presented and discussed in detail elsewhere.23 In short, this study showed that the polarizing effect of the protein/solvent environment causes a stabilization of electronic charge on the sulfur atom, which effectively lowers the energy of the Fe-S σ* orbital. Since this orbital has to accommodate an unpaired electron in the quartet and sextet states but not in the doublet state, the former are stabilized significantlysin the
J. Phys. Chem. B, Vol. 108, No. 28, 2004 10087 ferric aqua complex this accounts for a shift of up to 10 kcal/ mol with respect to the gas phase. Our results allow for a tentative interpretation of lowtemperature ENDOR results on 1,4 which identified two intermediate (nonequilibrium) species during product formation, presumably with shorter Fe-O distances than the X-ray value2 of 2.67 Å. The QM/MM calculations show the presence of at least one such minimum on the 2A surface with an Fe-O bond length around 2.2 Å, which probably corresponds to one of these intermediate species. Because of the flat potential along the dissociation of the Fe-O bond, it is conceivable that additional local minima exist on the way to the equilibrium species, resulting from subtle variations of the binding pocket conformation. 5. Conclusion The present paper presents QM/MM calculations on cytochrome P450cam complexed with its catalytic product 5-exohydroxycamphor. B3LYP/CHARMM and BLYP/CHARMM potential energy scans along the internal coordinate r(Fe-O) in the range from 1.94 to 3.13 Å were performed in the doublet, quartet, and sextet spin state. Free geometry optimizations from the points of lowest energy in the scan located the global minima of the respective spin surfaces. Analogous calculations in the gas phase served to assess the influence of the protein matrix on the QM region. The results show that the protein environment (i) facilitates the cleavage of the Fe-O bond by favorable interactions of hydroxycamphor with the binding pocket and (ii) lowers the quartet and sextet states relative to the doublet, which has consequences for the intersection points of the spin surfaces. The calculations are consistent with ENDOR results, which suggest that the immediate product of the catalytic conversion involves a species with a short Fe-O bond and that the active site subsequently relaxes to accommodate a species with a longer Fe-O bond. Acknowledgment. This research is sponsored by the Binational German Israeli Foundation (GIF). J.C.S. thanks the Fonds der Chemischen Industrie for a Kekule´ scholarship. References and Notes (1) Cytochrome P450: Structure, Mechanisms and Biochemistry, 2nd ed.; Ortiz de Montellano, P. R., Ed.; Plenum Press: New York, 1995; Vol. 2. (2) Li, H.; Narasimhulu, S.; Havran, L. M.; Winkler, J. D.; Poulos, T. L. J. Am. Chem. Soc. 1995, 117, 6297. (3) Lipscomb, J. D. Biochemistry 1980, 19, 3590. (4) Davydov, R.; Makris, T. M.; Kofman, V.; Werst, D. E.; Sligar, S. G.; Hoffman, B. M. J. Am. Chem. Soc. 2001, 123, 1403. (5) Scho¨neboom, J. C.; Cohen, S.; Lin, H.; Shaik, S.; Thiel, W. J. Am. Chem. Soc. 2004, 126, 4017. (6) Scho¨neboom, J. C.; Lin, H.; Reuter, N.; Thiel, W.; Cohen, S.; Ogliaro, F.; Shaik, S. J. Am. Chem. Soc. 2002, 124, 8142. (7) 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.; Sjovoll, M.; Fahmi, A.; Scha¨fer, A.; Lennartz, C. J. Mol. Struct. (THEOCHEM) 2003, 632, 1. (8) (a) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (b) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 85. (c) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (9) Reiher, M.; Salomon, O.; Hess, B. A. Theor. Chem. Acc. 2001, 107, 48. (10) Salomon, O.; Reiher, M.; Hess, B. A. J. Chem. Phys. 2002, 117, 4729. (11) Hay, J. P.; Wadt, W. R. J. Chem. Phys. 1985, 82, 299. (12) Ditchfield, R.; Hehre, W. J.; Pople, J. A. J. Chem. Phys. 1971, 54, 724.
10088 J. Phys. Chem. B, Vol. 108, No. 28, 2004 (13) Bakowies, D.; Thiel, W. J. Phys. Chem. 1996, 100, 10580. (14) Antes, I.; Thiel, W. Hybrid Quantum Mechanical and Molecular Mechanical Methods; Gao, J., Ed.; ACS Symposium Series 712; American Chemical Society; Washington, DC, 1998; pp 50-65. (15) de Vries, A. H.; Sherwood, P.; Collins, S. J.; Rigby, A. M.; Rigutto, M.; Kramer, G. J. J. Phys. Chem. B 1999, 103, 6133. (16) (a) Ahlrichs, R.; Ba¨r, M.; Ha¨ser, M.; Horn, H.; Ko¨lmel, C. Chem. Phys. Lett. 1989, 162, 165. (b) Ahlrichs, R.; Ba¨r, M.; Baron, H.-P.; Bauernschmitt, R.; Bo¨cker, S.; Ehrig, M.; Eichkorn, K.; Elliot, S.; Furche, F.; Ha¨ser, M.; Horn, H.; Ha¨ttig, C.; Huber, C.; Huniar, U.; Kattanneck, M.; Ko¨hn, A.; Ko¨lmel, C.; Kollwitz, M.; May, K.; Ochsenfeld, C.; O ¨ hm, H.; Scha¨fer, A.; Schneider, U.; Treutler, O.; v. Arnim, M.; Weigend, F.; Weis, P.; Weiss, H. TURBOMOLE 5.5, University of Karlsruhe: Karlsruhe, Germany, 2002. (17) 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.
Lin et al. (18) CHARMM22 force field: MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; 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. (19) For details about the generation and geometric features of the snapshot structures, nonstandard force field parameters, and QM/MM geometry optimizations, see ref 6 and Supporting Information thereof. (20) Smith, W.; Forester, T. J. Mol. Graph. 1996, 14, 136. (21) Billeter, S. R.; Turner, A. J.; Thiel, W. Phys. Chem. Chem. Phys. 2000, 2, 2177. (22) Ogliaro, F.; Harris, N.; Cohen, S.; Filatov, M.; de Visser, S. P.; Shaik, S. J. Am. Chem. Soc. 2000, 122, 8977. (23) Scho¨neboom, J. C.; Thiel, W. J. Phys. Chem. B 2004, 108, 7468.
