Article pubs.acs.org/JPCB
Hydrogen Abstraction of Camphor Catalyzed by Cytochrome P450cam: A QM/MM Study Published as part of The Journal of Physical Chemistry virtual special issue “Mark S. Gordon Festschrift”. Rui Lai and Hui Li* Department of Chemistry and Nebraska Center for Materials and Nanoscience, University of Nebraska-Lincoln, Lincoln, Nebraska 68588-0304, United States S Supporting Information *
ABSTRACT: A combined quantum mechanics and molecular mechanics (QM/MM, QM = UB3LYP-D3, MM = AMBER) method is used to study the hydrogen abstraction reaction in P450cam catalyzed hydroxylation of camphor in the quartet state. Compared to QM/MM calculations in the literature, this study uses larger basis sets for the most important atoms at the active site and QM/MM Hessian harmonic frequency calculations to determine the standard Gibbs free energy of activation and kinetic isotope effect. The QM/MM covalent boundary is treated with a capping hydrogen atom method, which is simple and robust. An energy barrier of 21.3 kcal/mol and a standard free energy of activation of 16.8 kcal/mol are obtained for this hydrogen abstraction reaction. These values are similar to those reported in the literature, suggesting that when a general protocol is followed, QM/MM results are reproducible. It is found that using a sufficiently large basis set is important to minimize basis set errors.
I. INTRODUCTION Cytochrome P450 enzymes (P450s) are involved in a number of vital biological processes. In humans, for example, they are involved in the drug metabolism, detoxification of xenobiotics, and carcinogenesis.1−3 It is well-known that the hydroxylation of nonactive C−H to form C−O−H catalyzed by P450 enzymes has high stereo- and regio- selectivity at physiological conditions.2 As the first crystallized cytochrome P450 enzyme, P450cam from bacterium Pseudomonas putida can catalyze the 5exo-hydroxylation of camphor with a high efficiency and selectivity (Figure 1A).4 With the availability of high resolution X-ray crystal structures, P450cam has been used extensively in experimental and theoretical studies to understand the general catalytic mechanism of P450 enzymes.5,6 The structures and functions of P450 enzymes are highly similar in different organisms, so the results from P450cam have direct implications for other P450 enzymes, such as human P450 enzyme. Experimental mechanistic studies of P450s have been limited by the short lifetime of the intermediates produced during the catalytic processes.1,5,6 Quantum chemical electronic structure calculation at various approximation levels have provided critical information toward the understanding of their mechanism. A generally accepted mechanism for the hydroxylation reaction is the hydrogen abstraction−oxygen rebound mechanism,2,7,8 in which the hydrogen abstraction step is characterized by a hydroxylation transition state (Figure 1B) formed from an oxo-ferric active compound, named as compound I (Cpd I, Figure 1C). As the key active oxidant in © XXXX American Chemical Society
most of the reactions catalyzed by P450s, Cpd I was originally proposed from quantum chemical calculations, and then was confirmed by experiments.5,9,10 Among various computational quantum chemical methods, combined quantum mechanics, and molecular mechanics (QM/MM) methods11 are probably the best choice to investigate the P450 catalytic mechanisms because QM/MM methods can retain the enzyme environment and include solvent effects. In these QM/MM methods, the large enzyme/ solvent system is divided into two regions, a relatively small QM region representing the active site and a relatively large MM region representing the protein matrix and the solvent. Many QM/MM calculations have been performed to determine the free energy of activation for the C−H hydroxylation by Cpd I of P450cam, and the results vary from 8.2 to 17.7 kcal/mol for different substrates.9,12−22 For the same substrate, the calculated results may also vary in a relatively large range and differ significantly from experimental values. The experimental upper limit of the standard Gibbs free energy of activation is 15.0 kcal/mol for hydroxycamphor formation, as derived by using transition state theory (at 298.15 K) from the rate constant of 66 s−1 reported by Purdy et al.23 for wild type P450cam. Received: September 30, 2016 Revised: November 10, 2016 Published: November 11, 2016 A
DOI: 10.1021/acs.jpcb.6b09923 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
MM Hessian calculation has been used to determine the zeropoint energy and thermal free energy changes in the literature.
II. COMPUTATIONAL METHODS All calculations were performed with the quantum chemistry polarizable force field program (QuanPol),28 which is integrated into the General Atomic and Molecular Electronic Structure System (GAMESS)29,30 package. The general procedure and QM/MM setup are shown in Figure 2.
Figure 1. (A) Oxidation of camphor by P450cam producing 5-exohydrocamphor. (B) Hydrogen abstraction step. (C) Structure of P450 compound I (Cpd I). Figure 2. QM/MM calculation of P450cam catalyzed hydroxylation of camphor. After MM MD equilibration of 48 268 MM atoms, 1869 atoms (106 QM atoms and 1763 MM atoms) were optimized for the reactant state (RS), transition state (TS), and hydroxo intermediate (HYD) in the presence of the rest 46 399 MM atoms. “aDZ” represents the aug-cc-pVDZ basis set.
Clearly, different QM/MM methods may yield different results.6,15,21,22,24 In general, the differences in these theoretical calculations may be due to the use of different starting structure, different QM methods, basis sets, different numbers of QM atoms in the QM/MM calculations, and different QMMM interactions in the QM/MM methods. Possible improvements on the accuracy of QM calculations have been discussed.6 For example, Lonsdale et al.18,25 proposed that the discrepancy between theory and experiments may lie heavily in the inadequate treatment of dispersion interaction in some density functional theory (DFT) methods, such as B3LYP;26,27 Jerome et al.24 suggested that other intrinsic errors in DFT methods, such as lacking of nondynamic electronic correlation, may also affect the results in P450 calculations. In this paper, we report the results of a QM/MM calculation for the standard Gibbs free energy of activation in hydrogen abstraction of camphor catalyzed by Cpd I of P450cam. The calculations were performed with the quantum chemistry polarizable force field program (QuanPol),28 which is a full spectrum QM/MM program that can work with various standard QM methods such as density functional theory methods, perturbation theory methods, multiconfiguration selfconsistent field methods, and time-dependent density functional theory methods. QM/MM molecular dynamics simulation, geometry optimization (including transition state search), and Hessian frequency calculations can be performed with QuanPol. QuanPol requires minimum user interferences for QM/MM setups, especially QM/MM covalent boundary treatment. Compared to published QM/MM calculations of P450cam, this work used much larger basis sets in the QM calculation (especially for Fe, SCys357 and the oxidative O atoms), a simple and robust scheme to treat the QM/MM covalent boundary, and QM/MM Hessian calculations to determine the Gibbs free energy corrections and deuterium kinetic isotope effect. To the best of our knowledge, no QM/
The X-ray crystal structure 1DZ9.PDB10 for P450cam from Pseudomonas putida was obtained from the Protein Data Bank.31 All water molecules in 1DZ9.PDB were kept except for WAT2206, which is believed to be the water molecule formed by O2. It was shown that including a water molecule in the QM region can lower the calculated energy barrier.15 However, in this study we only consider the case when no water molecule is at the active site. Hydrogen atoms were added to the X-ray structure by using the WHAT IF web interface.32 We manually edited the file so Asp297 was neutral (protonated) and His355 was positively charged. The protein was described with the AMBER ff12SB protein force field.33−37 The water molecules in the crystal structure file were described with a three-point flexible and nonpolarizable water model named QP301.28 Heme, Fe, O, and camphor were described by using a simplified universal force field in QuanPol (keyword LOUT=1). This force field uses Lennard-Jones (LJ) parameters very similar to those in the AMBER force field. The atomic charges were determined by using accumulated bond polarization charges, and then were modified for metal ions, metal ligand atoms, and ionized functional groups such as ammonium and carboxylate. The AMBER ff12SB atomic charge for free thiolate SCys357 was −0.8844 e, which was manually reduced to −0.3844 e together with an equal increase of the Fe ion charge (the final Fe charge thus is +3.3000 e). These partial atomic charges are consistent with the following formal charges (or oxidation states): O is −2, Fe is +5, S is −1, and each heme N is −0.5 e. In addition, we note that the atomic charges assigned by using accumulated bond polarization charges followed by manual editing are B
DOI: 10.1021/acs.jpcb.6b09923 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
Figure 3. QM region with 106 atoms in the QM/MM system (Fe, orange red; S, yellow; O, red; C, tan; H, white).
comparable to those in the AMBER force field. With −1 e charge for each of the two carboxylate groups on heme, the total charge of the heme/camphor is −2 e. These two negative charges are stabilized by three nearby positively charged residues: Arg112, Arg299, and His355. Arg112 is at the protein surface so its positive charge is partially screened by the bulk solvent. As mentioned above, the nearby Asp297 is neutral. The total charge of the protein and substrate is −14 e. The AMBER ff12SB force field for the protein and the simplified universal force field file for the substrates were prepared separately and then combined (QuanPol keyword ICOMBIN=1). After the combination, an additional bond stretching potential [force constant k = 300 (kcal/mol)/Å2 for E = k(r − r0)2 and r0 = 2.271 Å] for Fe−SCys357 and an additional bond angle bending potential [force constant k = 50 (kcal/mol)/rad2 for E = k(θ − θ0)2 and θ0 = 111.5°] for Fe−SCys357−CCys357 were added manually. The protein and substrate were then solvated in an 80 Å × 70 Å × 90 Å periodic boundary water box filled with 13 600 water molecules, 30 Na+ ions, and 16 Cl− ions. The total number of atoms was 48 268, with a zero net charge. All water molecules (the crystalline waters in the PDB file and the 13 600 added waters) were described with the three-point flexible and nonpolarizable water model QP301. The system was equilibrated for one million steps (1 ns) by using force field molecular dynamics (MD) simulation. This time length is sufficient for equilibrating the system because no statistic sampling is required. The periodic boundary condition was used in the MD simulation with a shifting function (QuanPol keyword ISHIFT=4 with a cutoff distance of 12.0 Å) for the charge−charge interactions, and a switching function (QuanPol keyword ISWITCH=1 with the switching range from 10.0 to 12.0 Å) for the LJ interactions. The constant particle number, pressure, and temperature (NPT) ensemble was used in the MD simulation with Berendsen barostat and thermostat38 (Pbath = 1.00 bar and Tbath = 298.15 K). The MD simulations was run by using the Beeman integration algorithm.39 After one million MD steps, the volume stabilized at 78.1 Å × 68.3 Å × 87.8 Å, and the overall geometry of the protein and substrate was similar to that in the original X-ray structure. The root-mean-square deviation of all backbone atom
positions in the last snapshot of MD simulation from the original X-ray structure is 0.665 Å. On the basis of the geometry at the last step of the force field MD simulation, QM/MM style spin-unrestricted open-shell B3LYP (with Grimme’s empirical dispersion correction version 3,40 denoted as UB3LYP-D3 in this paper) geometry optimization and transition state search were performed for 1869 atoms (106 QM atoms and 1763 MM atoms) in the presence of the rest 46 399 MM atoms. The product of the hydrogen abstraction reaction is a hydroxo intermediate (HYD) and camphor radical. This HYD state was optimized by using the same QM/MM settings as for the RS and TS. The 1869 atoms were selected by drawing a 16.0 Å radius sphere around the 5-exo-H atom of camphor as in the last step of the MM MD simulation, so some solvent water molecules were also selected. The same 1869 atoms were optimized for the reactant state (RS), transition state (TS), and product state (HYD) so the optimized structures and energies were comparable. In the QM/MM calculation, the QM region had 106 atoms: the heme, Fe, O, camphor, and the side chain of Cys357 as SCH3 (one of the three H atoms is the capping H atom that replaces the Cα of Cys357), as shown in Figure 3. All of these 106 QM atoms were within the 1869 atoms for optimization. The periodic boundary condition was used in the QM/MM system with the box size fixed at the value of the last step of the MM MD simulation. The same shifting and switching functions were used for MM charge−charge and LJ interactions as in the force field MD simulation. The QM-MM interaction used a special switching function41 in the range from 22.0 to 32.0 Å. The default capping H atom method implemented in QuanPol was used to treat the covalent bonds between QM and MM regions (here only the Cys357 residue of P450cam).28 In the QM/MM geometry optimization the augcc-pVDZ42,43 basis set (denoted as aDZ) was used for 12 most important QM atoms: the Fe atom and its ligand atoms (S, N, N, N, N, O), and the two reacting atoms (5-exo-H and C5) of camphor and three closest atoms (5-endo-H, C4, C6 of camphor). The aug-cc-pVDZ basis set for Fe is a 7s6p4d2f basis set (polarizable triple-ζ quality with additional diffuse s, p, d, f functions). The 6-31G* basis set44 was used for all the other 94 C
DOI: 10.1021/acs.jpcb.6b09923 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B QM atoms, which are H, C, and O atoms. The total number of Cartesian Gaussian basis functions was 1121. P450 Cpd I has two low lying electronic states with spin quantum number S being 1/2 (doublet state) or 3/2 (quartet state). We performed QM/MM style UB3LYP-D3 test calculations for both doublet and quartet states. For the doublet state, the expectation value of the square of the spin ⟨S2⟩ should be 3/4 (or 0.75). However, the computed ⟨S2⟩ from QM/MM style UB3LYP-D3 calculation is about 1.80, showing a large spin contamination from higher spin states, mainly the quartet state. This large spin contamination is not sensitive to the geometry of the QM region and is similar in the reactant state, transition state, and product state. The spin restricted open shell ROB3LYP method should not be used for P450 Cpd I because it forces all paired α and β electrons to occupy the same spatial orbitals. Therefore, neither the UB3LYP nor ROB3LYP method should be used to study P450 Cpd I in the doublet state. Instead, multiconfiguration methods should be used. For the quartet state, both ROB3LYP and UB3LYP can be used. However, we found that ROB3LYP self-consistent-field calculation is very difficult to converge when a large basis set is used. Therefore, we used UB3LYP method in this study, which is easy to converge, and allows for spin polarization. For the quartet state, the ⟨S2⟩ should be 15/4 (or 3.75). The QM/MM style UB3LYP-D3 value is about 3.80, showing a small spin contamination. In addition, we found that the QM/MM style UB3LYP-D3 calculated energies for the doublet state (with large spin contamination) are within a few kcal/mol to those for the quartet state (with small spin contamination). Therefore, the UB3LYP calculation for the quartet state can provide meaningful energetics for understanding the abstraction reaction. In QuanPol, QM/MM Hessian calculations (after geometry optimization) were performed only for QM atoms in the presence of MM atoms, with QM-MM interactions being the same as in the precedent QM/MM geometry optimization calculation and TS search. All QM atoms can be displaced to set up the force constance matrix (FCM), but in QuanPol we also implemented a partial Hessian method45 in that a smaller number of QM atoms can be selected to construct part of the FCM. Based on the UB3LYP-D3/AMBER optimized geometries of the RS and TS, partial Hessian calculations were performed for 12 most important QM atoms to estimate the zero-point energy and thermochemical properties. These 12 atoms were the same 12 atoms that used the aug-cc-pVDZ basis set, as described in prior paragrahs. To get the second derivative of the energy, the double displacement method was used with a step size of 0.01 bohr. A total of 73 QM/MM single-point energy and gradient calculations were required in each of the partial Hessian calculations for the RS and TS. The partial Hessian calculation confirmed that the reactant state (RS) had no imaginary frequency and transition state (TS) had only one imaginary frequency, which is 2050.414i cm−1. Because the Hessian calculation is costly, it was not performed for the HYD state.
Table 1. UB3LYP-D3/AMBER Optimized Interatomic Distances (Å) at the P450cam Active Site HYD TS RS
Fe−SCys357
Fe−O
O···H5‑exo
H5‑exo···C5
2.296 2.342 2.553
1.818 1.755 1.629
0.974 1.208 2.338
2.104 1.334 1.091
1.629 to 1.755 and to 1.818 Å, indicating a weakening in the bond strength due to the addition of the H to the O. To compensate for this change, the Fe−SCys357 bond length is shortened from 2.553 to 2.342 and to 2.296 Å, meaning a stronger bond is formed. These bond length changes are accompanied by the redistribution of the α spin electrons (see discussion in the next subsection). It is interesting to note that in the TS the O−H5‑exo distance is 1.208 Å, longer than the 0.974 Å in the HYD by 0.234 Å, and that the H5‑exo−C5 distance is 1.334 Å, longer than the 1.091 Å in RS by 0.243 Å. The TS structure is like an “average” of those for the RS and HYD. The UB3LYP-D3/AMBER optimized RS, TS, and HYD relative energies are 0.0, 21.3, and 9.7 kcal/mol, respectively. Therefore, the calculated energy barrier for the hydrogen abstraction is 21.3 kcal/mol. The free energy of activation can be estimated by including zero-point energy and thermal energies from the Hessian vibration analysis of the 12 atoms, which indicates that the total correction is −4.5 kcal/mol from RS to TS. The zero-point energy change is −4.7 kcal/mol, the major contributor of the correction. With the total correction of −4.5 kcal/mol, the free energy of activation ΔGa is 16.8 kcal/ mol. This QM/MM result, however, does not include the hydrogen tunneling effect, which, as suggested by very large deuterium isotope effects, can often increase the rate by a few times at 298.15 K and decrease the apparent free energy of activation. Zurek et al.21 estimated that the magnitude of tunneling correction could be as large as 2 kcal/mol. By using multiconfiguration molecular mechanics reaction dynamics calculation on testosterone 6β-hydroxylation by P450, Zhang et al. reported that the tunneling coefficients are around 3, indicating a considerable tunneling effect is involved in the reaction.46 However, such a quantitative analysis has not been reported for camphor in the literature. The transition state theory free energy of activation ΔGa for the hydrogen abstraction of camphor catalyzed by P450cam can be estimated by using the experimental kinetic rate constant. Purdy et al.23 measured the hydroxycamphor formation rate constant kcat to be 66 s−1 for P450cam at T = 298 K, from which a ΔGa can be derived as 15.0 kcal/mol. However, this rate constant kcat is the reflection of an electron transfer reaction, which is the rate-determining step of the multiple-step hydroxylation reaction. Therefore, the ΔGa for the hydrogen abstraction step (not rate-limiting) should be lower than 15.0 kcal/mol. Zurek et al.21 guessed that the hydrogen abstraction rate constant could probably be at least 1000 s−1, meaning a ΔGa of ∼13 kcal/mol. Our result 16.8 kcal/mol is 1.8 kcal/mol higher than the experimental upper limit 15.0 kcal/mol derived from a kcat of 66 s−1. Considering the hydrogen tunneling effect, our result could be lower, and closer to experiment. In general, the intrinsic errors in relative electronic energies computed with UB3LYP are typically on the order of a few kcal/mol. To further improve the computational result, more accurate multiconfiguration methods and/or perturbation treatment are required. In this
III. RESULTS AND DISCUSSION III.A. Free Energy of Activation. The coordinates of the 106 QM atoms in the UB3LYP-D3/AMBER optimized RS, TS, and HYD are presented in the Supporting Information. The relevant interatomic distances in the QM/MM optimized RS, TS, and HYD geometries are shown in Table 1. It is clear that from RS to TS and to HYD, the Fe−O bond is elongated from D
DOI: 10.1021/acs.jpcb.6b09923 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B Table 2. Löwdin Spin Population (e) on the Active Site Atomsa HYD TS RS a
cap-H
SCys357
Fe
O
H5‑exo
C5
Cpd I
camphor
0.020 0.031 0.057
0.086 0.350 0.564
1.813 1.415 1.334
0.188 0.551 0.718
0.046 0.081 0.003
0.635 0.307 0.003
2.028 2.473 2.990
0.972 0.527 0.010
A positive value means net α electron; a negative value means net β electron.
study, the starting geometry for the QM/MM geometry optimization is the last snapshot from MM MD simulation so the conformational influence is not examined. The fact that the MD equilibrated structure is in excellent agreement with the original X-ray structure (Computational Methods) suggests that the P450 enzyme (especially the active site) is relatively rigid, so it is expected that conformational influence would not be significant. In addition, 1869 atoms (106 QM and 1763 MM) within 16 Å to the active site were fully optimized, which should have further minimized conformational effects. The computed −4.7 kcal/mol zero-point energy correction is mainly due to the loss of the camphor C−H stretching mode ongoing from RS to TS, so the deuterium kinetic isotope effect would be significant. We performed UB3LYP-D3/AMBER partial Hessian calculation with deuterium mass (2.0141 amu) to estimate the deuterium kinetic isotope effect and found that deuterium substitution of the 5-exo hydrogen led to a kinetic isotope effect of 7.3. This value is obtained by using the harmonic frequencies from the UB3LYP-D3/AMBER calculation and may have errors from various sources such as lacking of anharmonicity. In general, if hydrogen and deuterium tunneling effects are considered, the deuterium kinetic isotope effects can be much larger. Indeed, experimentally measured deuterium kinetic isotope effects in P450 catalyzed hydroxylation reaction of alkanes vary from 7.5 to 19 ± 10.47−49 III.B. Spin Population. In the UB3LYP-D3 calculation a spin quantum number S = 3/2 is used, so there are three unpaired α electrons in the 106-atom QM region. The Löwdin spin populations at the most relevant atoms are listed in Table 2. Like all other electron population analysis methods, there is some arbitrarity in Löwdin spin population analysis. Despite the arbitrarity, it is clear that in the RS the three unpaired α electrons are mainly localized at the axial S−Fe−O moiety. Virtually all three unpaired α electrons are located at the heme (including Fe and O) and Cys357; camphor has a near zero (0.010 e) α spin population. The three α electrons are mainly localized at SCys357 (0.564 e), Fe (1.334 e), and O (0.718 e) atoms, as the sum of their Löwdin spin population is 2.616 e. In addition, the total α spin population on the Fe−O moiety is 2.052 e, the rest is mainly at the porphyrin ring (0.276 e) and Cys357 (0.662 e), with a total of 0.938 e. These results are virtually identical to those reported by other authors.6,50,51 However, we found that the three highest α spin orbitals do not correspond to the net α spin density population. As shown in Figure 4, the highest α spin orbital (HOMOα) is on the porphyrin ring, with no density at the S−Fe−O moiety (so there is a hole). The second highest α spin orbital (HOMOα− 1) is also on the porphyrin ring, but with some density at the SCys357 atom (no density at Fe and O). The third highest α spin orbital (HOMOα−2) is on one side of the porphyrin ring, also with some density at the SCys357 atom. Neither one of these three α spin orbitals has density at the Fe−O atoms. This is caused by the use of spin-unrestricted open shell UB3LYP method, which does not preserve the correct spin symmetry. Our single-point energy calculations using spin-restricted open
Figure 4. Three highest α spin orbitals and energies in the reactant state (RS, Cpd I with camphor) calculated with UB3LYP-D3[aug-ccpVDZ/6-31G*]/AMBER.
shell ROB3LYP method for the RS indeed led to orbitals virtually identical to those reported by other authors.6,50,51 Ongoing from RS to TS there is a shift of the α spin population from the SCys357 and O atoms to the C5 atom of camphor, whereas the α spin population on the Fe atom slightly increases from 1.334 e to 1.415 e. The spin population on the C5 atom is increased from 0.003 e in RS to 0.307 e in TS, whereas the spin populations on the SCys357 and O atoms are decreased from 0.564 e and 0.718 e in RS to 0.350 e and 0.551 e in TS, respectively (Table 2). The total spin population on the Fe−O moiety is 1.966 e. There is a noticeable α spin population, 0.081 e, at the 5-exo-H atom in the TS. The total α spin population on the camphor molecule (including 5-exo-H) in the TS is 0.527 e, with the rest 2.473 e on Cpd I (porphyrin, Fe, O, 2.069 e) and Cys357 (0.404 e). Passing the TS, the hydroxo intermediate HYD is formed with a large α spin population on the Fe atom (1.813 e) but small α spin populations on the ligand SCys357 (0.086 e) and O (0.188 e). The total spin population on the Fe−O moiety is 2.001 e. The C5 atom of camphor also has a large α spin population of 0.635 e. As a result, the total α spin population on the camphor molecule (including 5-exo-H) in the HYD is 0.972 e, with the E
DOI: 10.1021/acs.jpcb.6b09923 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
and no effect on any MM atoms. As a result, the Cys357 Cα− Cβ bond lengths in the QM/MM optimized RS, TS, and HYD structures are 1.540, 1.538, and 1.536 Å, respectively, all very close to the typical C−C single bond length 1.54 Å. The spin populations at the capping H atom (which is bonded to the Cβ atom of Cys357) in RS, TS, and HYD are 0.057 e, 0.031 e and 0.020 e, respectively. One of the normal H atoms on the Cβ atom of Cys357 has spin populations varying from 0.022 to 0.014 and to 0.001 e, from RS to TS and to HYD, suggesting that the capping H atom behaves slightly different from the normal H atoms. This is not surprising because the capping H atom is ∼1.54 Å away from the Cβ atom, whereas the normal H atoms are ∼1.1 Å away. The applied electronic repulsion potential on the capping H atom to achieve the desired ∼1.54 Å bond length would also have effects on the spin density distribution. Nevertheless, given the short distance from the capping H atom to the active site Fe and S atoms, the spin population at the capping H atom is small and should have insignificant influences on the overall spin density distribution in the system. III.D. Comparison to Other Calculations. In Tables 3 and 4 we compare our results with recently published QM/MM
rest 2.028 e on Cpd I (porphyrin, Fe, O, 1.918 e) and Cys357 (0.109 e). The total spin population on the Fe−O moiety remains at ∼2.0 e in the RS, TS, and HYD, although the Fe−O distance varies. The α spin population on camphor changes from 0.010 e in RS to 0.527 e in TS and to 0.972 e in HYD, suggesting that the simple UB3LYP-D3 method can give a fair description of the spin transfer from Cys357 to camphor. The half α electron (0.527 e) on camphor in the TS suggests that the calculated TS is not only the “average” of the molecular structures of the RS and HYD but also an “average” of their electronic structures. Earlier QM/MM calculations performed by Guallar et al.9,13 resulted in a quite low energy barrier of 11.7 kcal/mol for camphor hydroxylation by P450cam, and a shift of spin density to the propionate side chain. However, such a spin shift was actually caused by artifacts in the QM/MM calculation.6 In our QM/MM calculations, there are virtually zero net α spin populations at the propionate side chain atoms in the RS, TS, and HYD (data not shown). III.C. Capping H Atom Method. Both the link atom method and frozen orbital method have been used in QM/MM calculations of P450cam to treat covalent boundaries. In the frozen orbital method some localized orbitals were put in the boundary atoms and were absent from the DFT self-consistent field (SCF) iteration. In the link atom scheme (with charge shifting), additional degrees of freedom are introduced and sophisticated pseudopotentials are required for the QM/MM boundary atoms. Both of these methods can potentially introduce strong interactions that can alter the final QM/ MM results. For example, when frozen orbitals are used for boundary atoms, the QM self-consistent field calculation cannot be extended to the frozen orbitals, so the forces on the boundary atoms cannot be readily evaluated. As a result, the boundary atoms cannot move freely in the QM/MM calculation. In the link atom scheme, very often the charges of the boundary MM atoms must be shifted to avoid overpolarization of the QM wave function at the boundary. In the current work, we used the capping H atom method to treat the QM/MM boundary so no additional degrees of freedom are introduced and no pseudopotentials are required. When the capping H atom method is used, all QM and MM atoms can move freely in the MD simulation and geometry optimization, as well in Hessian frequency calculations. Although a robust QM/MM boundary treatment may not necessarily improve the accuracy of the QM/MM results, which is mainly determined by the QM level of theory, it does help eliminate potential errors and simplify the overall QM/MM preparation. In the QM/MM calculation, we used a capping H atom (denoted as CX in the coordinates listed in Supporting Information) to replace the Cα atom of Cys357 of P450cam. This capping H atom has a standard 6-31G* basis set (i.e., two s type Gaussian functions). This would lead to a C−H bond length of ∼1.1 Å. To maintain the correct Cα−Cβ bond length, an effective Gaussian repulsion potential was placed at the capping H atom in the form of a·exp(−b·R2), with a = 3.0 e/ bohr, b = 3.0 bohr−2, and R being the distance between the capping H atom (i.e., the center of the repulsion potential) and the electronic coordinate when one-electron integrals were computed for all basis functions in the QM calculation. This short-range repulsion potential acts on mainly the basis functions of the capping H atom itself and the Cβ atom of Cys357. It should have very small effect on other QM atoms,
Table 3. Comparison of Interatomic Distances (Å) in the Calculated Transition State this work Guallar et al.13 Schöneboom et al.12 Zurek et al.21 Lai et al.19,20
Fe−SCys357
Fe−O
O···H5‑exo
H5‑exo···C5
2.34 2.41 2.34 2.58 2.48
1.76 1.75 1.76 1.79 1.79
1.21 1.20 1.19 1.22 1.18
1.33 1.38 1.37 1.40 1.36
results for P450cam catalyzed camphor hydroxylation reaction. The TS geometry from this study is similar to those reported in the literature (Table 3),12,13,19−21 especially that reported by Schöneboom et al.12 The computed activation energy in this study is also similar to the literature values (Table 4), and again, especially that reported by Schöneboom et al.12 The similarity of these data suggest that when a general procotol is followed, QM/MM calculations performed with different programs and groups are reproducible. Our result is slightly higher than most of the other values (Table 4). One possible reason is the basis set effect. Unlike multiconfiguration and perturbation theory methods, which explicitly compute electronic dynamic correlation energies, it is usually believed that DFT methods do not necessarily require very large basis set. However, when a moderate basis set and a larger basis set give significantly different results, cautions must be exercised to ensure that the results are reasonably well converged with respect to basis set. Most of the published QM/ MM calculations used relatively small effective core potential basis sets for Fe in P450cam, and the 6-31G family basis sets for H, C, N, O, and S atoms (Table 4). Altun et al.15,22 reported several reaction activation energies for the hydrogen abstraction reaction of camphor catalyzed by P450cam. In one case the QM region had 120 atoms, including camphor, heme, Cys357, the CO group of Leu356, and the NHCαH unit of Leu358 (Wat903 excluded and Asp297 protonated). This case is similar to our current QM region with 106 atoms that include all camphor atoms, all heme atoms, and five atoms in Cys357 (Wat2206 excluded and Asp297 protonated). They reported hydrogen abstraction energy barriers of 18.1, 21.0, and 19.1 F
DOI: 10.1021/acs.jpcb.6b09923 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B Table 4. Comparison of QM/MM Computed Activation Energies for P450cam QM atoms experiment23 this work Guallar et al.13 Schöneboom et al.12
106
