Article pubs.acs.org/JPCB
Effect of Protein Environment within Cytochrome P450cam Evaluated Using a Polarizable-Embedding QM/MM Method Nandun M. Thellamurege and Hajime Hirao* Division of Chemistry and Biological Chemistry, School of Physical and Mathematical Sciences, Nanyang Technological University, 21 Nanyang Link, Singapore 637371 S Supporting Information *
ABSTRACT: Metalloenzymes accommodate cofactors and substrates in their active sites, thereby exerting powerful catalytic effects. Understanding the key elements of the mechanism via which such binding is accomplished using a number of atoms in a protein is a fundamental challenge. To address this issue computationally, here we used mechanical-embedding (ME), electronicembedding (EE), and polarizable-embedding (PE) hybrid quantum mechanics and molecular mechanics (QM/MM) methods and performed an energy decomposition analysis (EDA) of the nonbonding protein environmental effect in the “compound I” intermediate state of cytochrome P450cam. The B3LYP and AMBER99/QP302 methods were used to deal with the QM and MM subsystems, respectively, and the nonbonding interaction energy between these subsystems was decomposed into electrostatic, van der Waals, and polarization contributions. The PE-QM/MM calculation was performed using polarizable force fields that were capable of describing induced dipoles within the MM subsystem, which arose in response to the electric field generated by QM electron density, QM nuclei, and MM point charges. The present QM/MM EDA revealed that the electrostatic term constituted the largest stabilizing interaction between the QM and MM subsystems. When proper adjustment was made for the point charges of the MM atoms located at the QM−MM boundary, EE-QM/MM and PE-QM/MM calculations yielded similar QM electron density distributions, indicating that the MM polarization effect does not have a significant influence on the extent of QM polarization in this particular enzyme system.
■
INTRODUCTION Ubiquitously distributed cytochrome P450 enzymes (P450s) have the remarkable ability to catalyze a wide range of oxidation reactions that are involved in physiologically essential processes, such as the metabolism of organic compounds and the biosynthesis of steroid hormones.1−5 Over the past few decades, the reactive species and the metabolic reactions of P450s have garnered considerable interest. It is now widely believed that an oxoiron(IV) porphyrin π-cation radical intermediate, called compound I (Cpd I), is the reactive species that is responsible for the powerful catalytic activity of P450s. Although trapping the Cpd I intermediate in the normal turnover condition may still be a difficult challenge, accumulating experimental evidence supports the role of Cpd I as the reactive species.6−9 Computational studies using density functional theory (DFT) and hybrid quantum mechanics and molecular mechanics (QM/MM) methods have also played invaluable roles in advancing our understanding of the properties of elusive Cpd I.10−22 In addition to providing insights into the intrinsic reactivity of Cpd I, previous theoretical studies, particularly QM/MM studies, demonstrated that the stability and reactivity of Cpd I are significantly influenced by the protein environment. In recent years, we have attempted to analyze the effects of protein environment systematically.23−25 In general, the interior © 2014 American Chemical Society
environment of a protein is highly complicated because of the large number of surrounding atoms from amino acid residues. However, it is our expectation that some kind of “energy decomposition analysis (EDA)”, which was originally developed for smaller molecules,26,27 may provide intuitively interpretable descriptions of protein environments. In one study, one of the authors (H.H.) decomposed the effect of the protein environment in bacterial cytochrome P450cam (P450cam, see Figure 1) into electrostatic, van der Waals (vdW), and polarization terms.23 In so doing, the features of the mechanical-embedding (ME) and electronic-embedding (EE) schemes of the ONIOM (our own n-layered integrated molecular orbital and molecular mechanics) QM/MM method were exploited effectively.23,28,29 In another study, we used DFT-based EDA schemes to analyze the interaction between Cpd I and a water molecule,25 the importance of which had been suggested by QM/MM and DFT studies.18,19,30 Although many studies have applied QM/MM methods to enzymatic systems and their chemical reactions, almost all previous QM/MM calculations used an EE scheme of QM/ MM that takes into account the polarization of the QM Received: December 22, 2013 Revised: January 30, 2014 Published: February 2, 2014 2084
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
Article
polarizabilities for the amino acid residues in the MM subsystem and the K+ ion, whereas the QP302 polarizable force field in QuanPol was used for the crystallographic water molecules.31 The nonbonding interaction energy between the MM subsystem and the QM subsystem, ΔE, consists of electrostatic, vdW, and polarization contributions (eq 1): ΔE = Ees + Evdw + Epol
(1)
The electrostatic term (Ees) arises from the interaction of the QM electrons and nuclei with the partial charges of the MM subsystem. The effect of MM charges can be included in the Kohn−Sham one-electron operator. The DFT energy of the QM system in the gas phase is represented by the equation E0[ρ0 ] = Ts[ρ0 ] + J[ρ0 ] + E XC[ρ0 ] +
∫ v(r)ρ0 (r) dr + Enucl
(2)
where TS, J, and EXC are the kinetic energy of noninteracting electrons, the Coulomb potential energy, and the exchangecorrelation energy, respectively, and v, ρ0, and Enucl are the external potential caused by nuclei, the gas-phase electron density, and the nuclear repulsion energy, respectively. When the QM atoms are immersed in MM protein atoms without allowing the QM electron density to relax, the energy of the QM subsystem may change to E1 = E0[ρ0 (r)] −
∑∫ i
qi ·ρ0 (r) |ri − r|
dr +
∑∑ i
A
qi ·ZA |ri − rA| (3)
Figure 1. P450cam in the Cpd I state from the Protein Data Bank (PDB) (PDB code 1DZ9):7 (a) the entire enzyme; (b) a closeup view of the active site.
where qi and ri are the atomic charge and the position of MM atom i, respectively, r is the position of an electron, and ZA and rA are the charge and the position of nucleus A, respectively. Because the gas-phase QM density was used here, E1 is regarded as the QM energy of ME-QM/MM. The pure electrostatic effect on the QM subsystem by the MM subsystem, Ees, can be evaluated as the difference between E1 and E0
electron density in the presence of fixed atomic charges of classically described protein atoms. This is particularly true for studies of P450s. However, because the protein environment will be polarized by other atoms in reality, it would be interesting to perform QM/MM and EDA calculations for a P450 enzyme while taking into account the polarization effect on the MM side. In this paper, we performed such calculations for P450cam Cpd I (Figure 1) using polarizable-embedding (PE) QM/MM calculations. ME- and EE-QM/MM were also employed in combination to derive insights into the environmental effect within this enzyme.
Ees = E1 − E0
(4)
The vdW interaction between the QM and MM subsystems (EvdW) is expressed in the same way as that presented in ref 23. The MM partial charges would polarize the density from ρ0 to ρ1, and accordingly, the energy would change to E2. E2 can be obtained from iterative EE-QM/MM self-consistent field (SCF) calculations
■
THEORETICAL METHODS The geometry of P450cam in the doublet spin state, which was previously optimized using the ONIOM(DFT:MM) method, was used for the current QM/MM study.23 Units comprising the QM subsystem were also chosen in the same way: namely, protoporphyrin IX, oxoiron(IV), deprotonated Cys357, and camphor. To describe the polarization of the protein surroundings in PE-QM/MM calculations, the QuanPol polarizable force field program was used,31 which was recently implemented by Thellamurege et al. in the GAMESS quantum chemistry package.32,33 The unrestricted B3LYP method and the 6-31G* basis set were used for QM calculations. The remaining amino acid residues of the enzyme, water molecules, and the K+ ion were described by a polarizable force field. The AMBER9934 force field was used in combination with atomic
E2 = E0[ρ1(r)] −
∑∫ i
qi ·ρ1(r) |ri − r|
dr +
∑∑ i
A
qi ·ZA |ri − rA| (5)
The difference between E2 and E1 gives the polarization energy of the QM subsystem (Epol,EE) Epol,EE = E2 − E1
(6)
Several authors have pointed out the importance of using polarizable force fields in QM/MM calculations for the description of biological molecules.35−46 The current study used the induced dipole method implemented in QuanPol to describe the polarization effects of the MM atoms. This method utilizes polarizabilities placed on atomic sites, thereby allowing 2085
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
Article
the induction of dipoles.31,47,48 A dipole induced at a point depends on the polarizability and the electric field felt at that point, i.e.
d i = αi Fi
Scheme 1. Flowchart of Our QM/MM EDA Scheme
(7)
Here, di is the induced dipole vector generated at point i, Fi is the total field felt at point i, and αi is the nine-component polarizability tensor placed at point i ⎛ αi , xx αi , xy αi , xz ⎞ ⎜ ⎟ αi = ⎜ αi , yx αi , yy αi , yz ⎟ ⎜α ⎟ ⎝ i , zx αi , zy αi , zz ⎠
(8)
The polarization energy associated with the induced dipoles can be calculated using eq 9, as the interaction between the induced dipoles and the electric field
1 E id = − FT ·d 2
(9)
An H-link atom was appended at the QM-MM covalent boundary to satisfy the valency of the QM subsystem. As Leu356 and Leu358 are the amino acid residues that are directly bonded to Cys357 in the QM subsystem, the carbonyl C of Leu356 and the amide N of Leu358 were replaced by H atoms. Subsequently, the N−H and C−H bond lengths were shortened to 0.947 and 1.063 Å from the original N−C and C− N bond lengths (1.353 and 1.352 Å, respectively) to prevent instability of energy calculation. The distance scaling adopted here is the same as that used in ref 23. An MM point charge potentially gives rise to an electrostatic interaction that is too large; therefore, it may overpolarize the electron density of the QM subsystem, especially when it is located very close to the QM subsystem.46 One common approach that is used to circumvent this problem is the adjustment of the charges of MM atoms at the QM-MM boundary. In this study, four slightly different methods of charge adjustment (methods 1−4) were examined and compared. In method 1, the original MM charges were used without making any adjustment. In method 2, the charges of MM atoms located up to three bonds away from the QM atom at the boundary were eliminated, and the eliminated charges were then redistributed evenly to the adjacent atoms within the same amino acid residue. For example, as shown in Scheme 2a, the point charges of carbonyl C, carbonyl O, Cα, amide N, Cβ, and Hα of Leu356 were made zero. The total charge eliminated was −0.5047e; thus, a charge of −0.1262e was added to each of the four adjacent atoms, i.e., H, two Hβ, and Cγ. The charge redistribution allowed the total charge of a given residue to be conserved even after it was subjected to the charge adjustment. In method 3, however, this charge redistribution step was omitted. In method 4, only the charges of MM atoms that bind covalently to QM atoms were made zero, whereas other charges were left unchanged. Scheme 2 illustrates the MM partial charges of Leu356 and Leu358 after applying these four methods to the MM charges. Zero polarizabilities were placed on atoms of which the charges were made zero.
where d is the vector of induced dipoles of MM atoms and FT is the transpose of the total field vector at the MM atoms. When the effect of the QM atoms is included in the F vector, mutual polarization of the QM and MM subsystems can be described. As the induced dipoles are dependent on each other, they must be determined in an iterative manner. Typically, in a PE-QM/MM single point energy calculation, this iteration is executed at each SCF step. The induced dipoles polarize the QM electron density, which in turn changes the magnitudes of the induced dipoles of the MM atoms. These two types of iterative calculations are continued until self-consistency is fulfilled. We employed a damping function, as implemented in QuanPol, to scale down the polarizability of an MM atom that was located within a predefined distance of 3.0 Å from any of the QM atoms (see the Supporting Information).31 This polarization energy may also be calculated without including any contributions of QM atoms to F. Using the energy thus calculated (Eid,MM) and Eid, the change in polarization energy, which arises as a result of the interaction between the QM and MM subsystems, Epol,id, is calculated as Epol,id = E id − E id,MM
(10)
The PE-QM/MM calculation yields the QM density ρ2, and the QM energy E3 can be obtained by replacing ρ1 with ρ2 in eq 5. The polarization energy of the QM subsystem for the PE-QM/ MM calculation (Epol,PE) can be obtained as the difference between E3 and E1. Epol,PE = E3 − E1
(11)
The procedure used for these energy decomposition calculations is summarized in Scheme 1. To calculate the QM-MM electrostatic energy, a previous study used the classical Coulomb’s law formula.23 Although that procedure was in line with the treatment of electrostatic interaction energy in the ONIOM-ME method, it required the determination of gas-phase point charges of the QM atoms in a manner that reproduces the electrostatic potential. To skip this process and simplify the calculation of the electrostatic energy, instead we used eqs 2−4 to calculate Ees. Thus, the electrostatic interaction was calculated as the interaction between the point charges in the MM subsystem and the electron density and nuclei in the QM subsystem.
■
RESULTS AND DISCUSSION Table 1 summarizes the electrostatic, polarization, and vdW interaction energies for P450cam Cpd I, which were calculated using the EDA procedure and applying the four QM/MM charge adjustment methods. 2086
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
Article
Scheme 2. MM Charges (e) on Leu356 and Leu358 after Applying the QM/MM-Boundary Charge Redistribution Methods: (a) Method 1, (b) Method 2, (c) Method 3, and (d) Method 4a
a
Adjusted charges are colored blue.
Epol,EE was calculated as −105.5 kcal/mol, which was close to the previous value. Figure 2a illustrates how methods 1 and 2 gave different Mulliken charge distributions in the QM subsystem. One can see clearly in the figure that the charges of Cys357, which is bonded directly to Leu356 and Leu358 in the MM subsystem (Scheme 2), exhibited large differences. The overpolarization also affected the spin distribution in the QM subsystem (Figure 2b). A particularly large difference in spin population was seen for the sulfur atom in Cys357. Modest differences were observed for the atoms of the porphyrin ring, despite the fact that the porphyrin was not very close to Leu356 and Leu358. This result suggests a significant influence of the overpolarization at the sulfur atom on the entire singly occupied a2u-type orbital, which has significant amplitudes on the porphyrin as well.14,15 The electrostatic interaction was diminished by 108−123 kcal/mol, when three different charge adjustment methods (2− 4) were applied (see Table 1). Nevertheless, consistent with a previous conclusion,23 the role played by the electrostatic effect in the nonbonding interaction between the QM and MM subsystems was the largest. A comparison between the results obtained using methods 2 and 3 showed that the electrostatic stabilization and the polarization were slightly larger when the eliminated charges were not redistributed to neighboring atoms. In contrast, method 4, which simply eliminates the charges of only the atoms at the QM-MM boundary, gave electrostatic and polarization stabilization effects that are slightly smaller than those obtained using method 2. Overall, the three methods yielded more or less similar energy values, and the values were not very sensitive to the charge adjustment method employed. Therefore, henceforth, we will mainly use the results obtained using method 2 for further discussions.
Table 1. Decomposed QM/MM Interaction Energy Terms (kcal/mol) for P450cam Cpd I Obtained Using Various Charge-adjustment Methods term
method 1
method 2
method 3
method 4
Ees EvdW Epol,EE Epol,PE Epol,id
−331.3 −105.8 −294.1 −293.4 +62.1
−214.4 −105.8 −105.5 −105.0 +65.8
−222.9 −105.8 −112.7 −112.1 +69.1
−208.7 −105.8 −99.4 −99.0 +59.6
The electrostatic energy was calculated as −331.3 kcal/mol using method 1. This value agrees reasonably well with the value obtained previously (−289.9 kcal/mol).23 The difference between these values can be attributed partly to the different DFT methods and MM force fields used. A more important reason may be that the classical Coulomb’s law formula was used previously for the calculation of Ees; thus, the pairwise interaction energy between MM and QM point charges was scaled down when two atoms were located not more than three bonds away from each other, whereas this sort of scaling was not applied in the present approach. The polarization stabilization of the QM density (Epol,EE: −294.1 kcal/mol) obtained using method 1 was much larger than the stabilization estimated previously (−101.6 kcal/mol). This significant discrepancy arose because the QM subsystem feels MM charges in the immediate proximity (Scheme 2a) in the current EE-QM/MM calculation using method 1, whereas the effect of a few point charges at the QM-MM covalent boundary on the QM density are ignored in the ONIOM-EE calculation.28 To avoid the overpolarization observed in method 1, some type of charge adjustment was deemed appropriate. In fact, the overpolarization was mitigated significantly using method 2: 2087
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
Article
Figure 2. Deviation of the Mulliken charge (a) and Mulliken spin population (b) values of the QM atoms obtained using method 1 from those obtained using method 2. See Scheme S1 for the specific atom numbering.
Since Ees is additive, it can be further decomposed into contributions from individual residues as presented in Figure 3. Such decomposition can be performed by repeating the calculation according to eq 4, while assigning charges of only one residue at a time. The calculated energy contributions were very similar to those reported in ref 23, despite the fact that the DFT methods, force fields, and Ees formulas used were different in these two studies. For example, Arg112, Arg299, and His355, which are located in the close proximity of the carboxylate groups of the protoporphyrin ligand (see Figure 1b), showed large stabilizing electrostatic effects. Tyr96, which forms a
hydrogen bond with the carbonyl group of camphor, had an Ees value of −9.5 kcal/mol, which was comparable with the value calculated in ref 23 (−8.9 kcal/mol). However, large differences were observed for Leu356 (see also Figure S1). This discrepancy arose because, in the current approach, a zeroed charge (e.g., that of the carbonyl carbon of Leu356) in method 2 could never have any electrostatic interaction with the QM subsystem. In contrast, in the previous classical approach, the same carbonyl carbon, for which no MM charge adjustment was made, could have a nonzero pairwise electrostatic interaction with atoms in the QM subsystem. This explains 2088
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
Article
Figure 3. Comparison of electrostatic energy contributions from amino acid residues 10−199 (a) and 200−414 and the potassium ion (#415) (b) obtained in this work with method 2 (blue) and in a previous work (red).23
why the Ees value between the current approach (with method 2, −214.4 kcal/mol) and the previous approach (−289.9 kcal/ mol) differed by ∼76 kcal/mol. The calculated vdW interaction energy was equal to the value obtained previously because the same method was used; thus, a further explanation is not provided here. Table 1 shows that Epol,PE was smaller in magnitude than Epol,EE by 0.5 kcal/mol, which suggests that the polarization of the QM subsystem was slightly suppressed in the PE-QM/MM calculation. However, the difference between these polarization
energies was fairly small; thus, the MM polarization seems to have almost no impact on the QM electron density. Even in the absence of the QM atoms, the MM atoms can have induced dipoles because of the existing electric fields generated by many MM atoms. The stabilization gained by the MM atoms (Eid,MM) amounted to >2900 kcal/mol (Table S2). The polarization energy caused by induced dipoles will change when QM atoms come into the active site. This differential effect of Eid can be quantified by Epol,id (see eq 10); the resulting value, +65.8 kcal/mol (with method 2), was positive. Therefore, it appears that the inclusion of the QM atoms 2089
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
Article
to 69.6 D in the PE-QM/MM calculation. This similarity in dipole moment provides another indication that the extent of polarization of the QM density is nearly equal in the EE-QM/ MM and PE-QM/MM calculations. We reported above that method 1 caused overpolarization. Using method 1, the dipole moment was calculated as 70.3 and 70.5 D with EE-QM/MM and PE-QM/MM, respectively. The values were somewhat larger than those obtained using method 2 (Table 2) because of the strong polarization effect exerted on the QM atoms for method 1.
somewhat diminished the polarization stabilization associated with induced dipoles. The key atomic spin populations obtained from different calculations are compared in Figure 4 (see also Figure S2 for
■
CONCLUSIONS In this work, we performed a QM/MM study of cytochrome P450cam in the Cpd I state to quantify the protein environmental effect of the enzyme. The B3LYP method was used for the QM subsystem, whereas AMBER99 and QP302 were used for the MM subsystem. ME-, EE-, and PE-QM/MM methods were applied in a stepwise manner to decompose the nonbonding QM−MM interaction energy into electrostatic, vdW, and polarization contributions. In the absence of adjustment for the atomic charges of the MM atoms, an overpolarization of the QM density was observed because the MM charges at the QM−MM boundary interacted with the QM atoms too strongly. However, this overpolarization effect was alleviated effectively by adopting charge-adjustment methods. The QM/MM energy decomposition analysis showed that the electrostatic interaction was the major factor that stabilized the binding of the QM subsystem to the protein. The polarization and vdW effects also contributed significantly to the binding process. However, the polarization stabilization gained by the MM atoms from induced dipoles was reduced when the protein accommodated the QM atoms. The QM polarization energies, distributions of charge and spin, and dipole moments all suggest that the electron density of the QM subsystem of this enzyme is polarized to a similar extent in EEQM/MM and PE-QM/MM calculations.
Figure 4. Key Mulliken atomic spin populations obtained with MEQM/MM (gas phase), EE-QM/MM, and PE-QM/MM. Some of the atom labels were taken from the PDB file (PDB code 1DZ9).
full details of the spin and charge distribution). As observed previously in ref 15, there was noticeable and unnatural spin localization at the propionate oxygen atoms in the gas-phase or ME-QM/MM calculation. Consequently, no significant spin population was seen on the porphyrin ring. In contrast, the spin populations on propionates, which interact with arginine residues (see Figure 1b), were shifted back to the porphyrin moiety in the EE-QM/MM calculation.15 As this shifted spin originates from the unpaired electron occupying the a2u-type orbital,11 significant spin populations were observed on the meso carbon and pyrrole nitrogen atoms of the porphyrin. The PE-QM/MM calculation produced a spin distribution that was similar to that of the EE-QM/MM calculation. This result indicates that EE-QM/MM and PE-QM/MM yielded similarly polarized QM electron densities. Table 2 lists the calculated dipole moments of the QM subsystem. Although the QM subsystem prefers a lower polarization in the gas-phase or ME-QM/MM calculation, the polarization was enhanced in the EE-QM/MM calculation because the polar nature of the protein environment was taken into account. As a result, the dipole moment increased from 27.3 to 69.9 D. The dipole moment was only slightly reduced
■
S Supporting Information *
QM energies at each step of Scheme 1, Eid and Eid,MM values, differences in Ees of a few residues between the current and previous work, differences in Mulliken charge and spin distributions between the ME- and EE-/PE-QM/MM methods, numbering of QM atoms, spin distributions obtained with a few more residues included in the QM region, and the damping function employed. This material is available free of charge via the Internet at http://pubs.acs.org.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail: [email protected] (H.H.). Notes
The authors declare no competing financial interest.
Table 2. Comparison of the Calculated Dipole Moments (D) of the QM Subsystem Obtained with Method 2
■
ACKNOWLEDGMENTS This work was supported by a Nanyang Assistant Professorship and an AcRF Tier 1 grant (RG14/12). We thank the High Performance Computing Centre at Nanyang Technological University for computer resources.
dipole moment (D) ME-QM/MM EE-QM/MM PE-QM/MM
ASSOCIATED CONTENT
27.3 69.9 69.6 2090
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
■
Article
Cytochrome c Pperoxidase, and Ascorbate Peroxidase. J. Comput. Chem. 2006, 27, 1352−1362. (21) Porro, C. S.; Sutcliffe, M. J.; de Visser, S. P. Quantum Mechanics/Molecular Mechanics Studies on the Sulfoxidation of Dimethyl Sulfide by Compound I and Compound 0 of Cytochrome P450: Which is the Better Oxidant? J. Phys. Chem. A 2009, 113, 11635−11642. (22) de Visser, S. P.; Ogliaro, F.; Sharma, P. K.; Shaik, S. What Factors Affect the Regioselectivity of Oxidation by Cytochrome P450? A DFT Study of Allylic Hydroxylation and Double Bond Epoxidation in a Model Reaction. J. Am. Chem. Soc. 2002, 124, 11809−11826. (23) Hirao, H. Energy Decomposition Analysis of the Protein Environmental Effect: the Case of Cytochrome P450cam Compound I. Chem. Lett. 2011, 40, 1179−1181. (24) Hirao, H. The Effects of Protein Environment and Dispersion on the Formation of Ferric-Superoxide Species in Myo-Inositol Oxygenase (MIOX): A Combined ONIOM(DFT:MM) and Energy Decomposition Analysis. J. Phys. Chem. B 2011, 115, 11278−11285. (25) Thellamurege, N.; Hirao, H. Water Complexes of Cytochrome P450: Insights from Energy Decomposition Analysis. Molecules 2013, 18, 6782−6791. (26) Kitaura, K.; Morokuma, K. A New Energy Decomposition Scheme for Molecular Interactions within the Hartree-Fock Approximation. Int. J. Quantum Chem. 1976, 18, 1558−1565. (27) von Hopffgarten, M.; Frenking, G. Energy Decomposition Analysis. WIREs Comput. Mol. Sci. 2012, 2, 43−62. (28) Vreven, T.; Byun, K. S.; Komaromi, I.; Dapprich, S.; Montgomry, J. A., Jr.; Morokuma, K.; Frisch, M. J. Combining Quantum Mechanics Methods with Molecular Mechanics Methods in ONIOM. J. Chem. Theory. Comput. 2006, 2, 815−826. (29) Chung, L. W.; Hirao, H.; Li, X.; Morokuma, K. The ONIOM Method: Its Foundation and Application to Metalloenzymes and Photobiology. WIREs Comput. Mol. Sci. 2012, 2, 327−350. (30) Kumar, D.; Altun, A.; Shaik, S.; Thiel, W. Water as Biocatalyst in Cytochrome P450. Faraday Discuss. 2011, 148, 373−383. (31) Thellamurege, N. M.; Si, D.; Cui, F.; Zhu, H.; Lai, R.; Li, H. QuanPol: a Full Spectrum and Seamless QM/MM Program. J. Comput. Chem. 2013, 34, 2816−2833. (32) Gordon, M. S.; Schmidt, M. W. Theory and Applications of Computational Chemistry; Dykstra, C. E., Frenking, K. S., Scuseria, G. E., Eds.; Elsevier: Amsterdam, The Netherlands, 2005. (33) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347−1363. (34) Wang, J.; Cieplak, P.; Kollman, P. A. How Well Does a Restrained Electrostatic Potential (RESP) Model Perform in Calculating Conformational Energies of Organic and Biological Molecules? J. Comput. Chem. 2000, 21, 1049−1074. (35) Warshel, A.; Levitt, M. Theoretical Studies of Enzymic Reactions: Dielectric, Electrostatic and Steric Stabilization of the Carbonium Ion in the Reaction of Lysozyme. J. Mol. Biol. 1976, 103, 227−249. (36) Russell, S. T.; Warshel, A. Calculations of Electrostatic Energies in Proteins: The Energetics of Ionized Groups in Bovine Pancreatic Trypsin Inhibitor. J. Mol. Biol. 1985, 185, 389−404. (37) Warshel, A.; Sharma, P. K.; Kato, M.; Parson, W. W. Modeling Electrostatic Effects in Proteins. Biochim. Biophys. Acta 2006, 1764, 1647−1676. (38) Thompson, M. A.; Schenter, G. K. Excited States of the Bacteriochlorophyll b Dimer of Rhodopseudomonas viridis: A QM/ MM Study of the Photosynthetic Reaction Center That Includes MM Polarization. J. Phys. Chem. 1995, 99, 6374−6386. (39) Illingworth, C. J. R.; Gooding, S. R.; Winn, P. J.; Jones, G. A.; Ferenczy, G. G.; Reynolds, C. A. Classical Polarization in Hybrid QM/ MM Methods. J. Phys. Chem. A 2006, 110, 6487−6497. (40) Zhang, Y.; Lin, H.; Truhlar, D. G. Self-Consistent Polarization of the Boundary in the Redistributed Charge and Dipole Scheme for
REFERENCES
(1) Omura, T.; Sato, R. A New Cytochrome in Liver Microsomes. J. Biol. Chem. 1962, 237, PC1375−PC1376. (2) Denisov, I. G.; Makris, T. M.; Sligar, S. G.; Schlichting, I. Structure and Chemistry of Cytochrome P450. Chem. Rev. 2005, 105, 2253−2278. (3) Ortiz de Montellano, P. R. Cytochrome P450: Structure Mechanism and Biochemistry, 3rd ed.; Ortiz de Montellano, P. R., Ed.; Kluwer Academic/Plenum Publishers: New York, 2005. (4) Sono, M.; Roach, M. P.; Coulter, E. D.; Dawson, J. H. Heme Containing Oxygenases. Chem. Rev. 1996, 96, 2841−2887. (5) Guengerich, F. P. Cytochrome P450 and Chemical Toxicology. Chem. Res. Toxicol. 2008, 21, 70−83. (6) Groves, J. T. High-Valent Iron in Chemical and Biological Oxidations. J. Inorg. Biochem. 2006, 100, 434−447. (7) Schlichting, I.; Berendzen, J.; Chu, K.; Stock, A. M.; Maves, S. A.; Benson, D. E.; Sweet, R. M.; Ringe, D.; Petsko, G. A.; Sligar, S. G. The Catalytic Pathway of Cytochrome P450cam at Atomic Resolution. Science 2000, 287, 1615−1622. (8) Davydov, R.; Makris, T. M.; Kofman, V.; Werst, D. E.; Sligar, S. G.; Hoffman, B. M. Hydroxylation of Camphor by Oxy-cytochrome P450cam: Mechanistic Implications of EPR and ENDOR Studies of Catalytic Intermediates in Native and Mutant Enzymes. J. Am. Chem. Soc. 2001, 123, 1403−1415. (9) Rittle, J.; Green, M. T. Cytochrome P450 Compund I: Capture, Characterization, and C-H Bond Activation Kinetics. Science 2010, 330, 933−937. (10) Meunier, B.; de Visser, S. P.; Shaik, S. Mechanism of Oxidation Reactions Catalyzed by Cytochrome P450 Enzymes. Chem. Rev. 2004, 104, 3947−3980. (11) Shaik, S.; Kumar, D.; de Visser, S. P.; Altun, A.; Thiel, W. Theoretical Perspective on the Structure and Mechanism of the Cytochrome P450 Enzymes. Chem. Rev. 2005, 105, 2279−2328. (12) Shaik, S.; Hirao, H.; Kumar, D. Reactivity Patterns of Cytochrome P450 Enzymes: Multifunctionality of the Active Species, and Two States-Two Oxidants Conundrum. Nat. Prod. Rep. 2007, 24, 533−552. (13) Shaik, S.; Cohen, S.; Wang, Y.; Chen, H.; Kumar, D.; Thiel, W. P450 Enzymes; Their Structure Reactivity and Selectivity − Modeled by QM/MM Calculations. Chem. Rev. 2010, 110, 949−1017. (14) Schöneboom, J. C.; Lin, H.; Reuter, N.; Thiel, W.; Cohen, S.; Ogliaro, F.; Shaik, S. The Elusive Oxidant Species of Cytochrome P450 Enzymes: Characterization by Combined Quantum Mechanical/ Molecular Mechanical (QM/MM) Calculations. J. Am. Chem. Soc. 2002, 124, 8142−8151. (15) Schöneboom, J. C.; Cohen, S.; Lin, H.; Shaik, S.; Thiel, W. Quantum Mechanical/Molecular Mechanical Investigation of the Mechanism of C-H Hydroxylation of Camphor by Cytochrome P450cam: Theory Supports a Two-State Rebound Mechanism. J. Am. Chem. Soc. 2004, 126, 4017−4034. (16) Cho, K.-B.; Hirao, H.; Chen, H.; Carvajal, M. A.; Cohen, S.; Derat, E.; Thiel, W.; Shaik, S. Compound I in the Heme Thiolate Enzymes: A Comparative QM/MM Study. J. Phys. Chem. A 2008, 112, 13128−13138. (17) Schöneboom, J. C.; Neese, F.; Thiel, W. Toward Identification of the Compound I Reactive Intermediate in Cytochrome P450 Chemistry: A QM/MM Study of Its EPR and Mossbaur Parameters. J. Am. Chem. Soc. 2005, 127, 5840−5853. (18) Altun, A.; Guallar, V.; Friesner, R. A.; Shaik, S.; Thiel, W. The Effect of Heme Environment on the Hydrogen Abstraction Reaction of Camphor in P450cam Catalysis: A QM/MM Study. J. Am. Chem. Soc. 2006, 128, 3924−3925. (19) Altun, A.; Shaik, S.; Thiel, W. What is the Active Species of Cytochrome P450 during Camphor Hydroxylation? QM/MM Studies of the Different Electronic States of Compound I and of Reduced and Oxidized Iron-Oxo Intermediates. J. Am. Chem. Soc. 2007, 129, 8978− 8987. (20) Harvey, J. N.; Bathelt, C. M.; Mulholland, A. J. QM/MM Modeling of Compound I Active Species in Cytochrome P450, 2091
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
Article
Combined Quantum-Mechanical and Molecular-Mechanical Calculations. J. Chem. Theory Comput. 2007, 3, 1378−1398. (41) Zhang, Y.; Lin, H. Flexible-Boundary QM/MM Calculations: II. Partial Charge Transfer across the QM/MM Boundary That Passes Through a Covalent Bond. Theor. Chem. Acc. 2010, 126, 315−322. (42) Wanko, M.; Hoffmann, M.; Frähmcke, J.; Frauenheim, T.; Elstner, M. Effect of Polarization on the Opsin Shift in Rhodopsins. 2. Empirical Polarization Models for Proteins. J. Phys. Chem. B 2008, 112, 11468−11478. (43) Frähmcke, J. S.; Wanko, M.; Phatak, P.; Mroginski, M. A.; Elstner, M. The Protonation State of Glu181 in Rhodopsin Revisited: Interpretation of Experimental Data on the Basis of QM/MM Calculations. J. Phys. Chem. B 2010, 114, 11338−11352. (44) Welke, K.; Frähmcke, J. S.; Watanabe, H. C.; Hegemann, P.; Elstner, M. Color Tuning in Binding Pocket Models of the Chlamydomonas-Type Channelrhodopsins. J. Phys. Chem. B 2011, 115, 15119−15128. (45) Wolter, T.; Welke, K.; Phatak, P.; Bondar, A. N.; Elstner, M. Excitation Energies of a Water-Bridged Twisted Retinal Structure in the Bacteriorhodopsin Proton Pump: A Theoretical Intvestigation. Phys. Chem. Chem. Phys. 2013, 15, 12582−12590. (46) Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem., Int. Ed. 2009, 48, 1198−1229. (47) Li, H. Quantum Mechanical/Molecular Mechanical/Continuum Style Solvation Model: Linear Response Theory, Variational Treatment and Nuclear Gradients. J. Chem. Phys. 2009, 131, 184103. (48) Li, H.; Gordon, M. S. Gradients of the Polarization Energy in the Effective Fragment Potential Method. J. Chem. Phys. 2006, 125, 194103−194112.
2092
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
Effect of Protein Environment within Cytochrome P450cam Evaluated Using a Polarizable-Embedding QM/MM Method Nandun M. Thellamurege and Hajime Hirao* Division of Chemistry and Biological Chemistry, School of Physical and Mathematical Sciences, Nanyang Technological University, 21 Nanyang Link, Singapore 637371 S Supporting Information *
ABSTRACT: Metalloenzymes accommodate cofactors and substrates in their active sites, thereby exerting powerful catalytic effects. Understanding the key elements of the mechanism via which such binding is accomplished using a number of atoms in a protein is a fundamental challenge. To address this issue computationally, here we used mechanical-embedding (ME), electronicembedding (EE), and polarizable-embedding (PE) hybrid quantum mechanics and molecular mechanics (QM/MM) methods and performed an energy decomposition analysis (EDA) of the nonbonding protein environmental effect in the “compound I” intermediate state of cytochrome P450cam. The B3LYP and AMBER99/QP302 methods were used to deal with the QM and MM subsystems, respectively, and the nonbonding interaction energy between these subsystems was decomposed into electrostatic, van der Waals, and polarization contributions. The PE-QM/MM calculation was performed using polarizable force fields that were capable of describing induced dipoles within the MM subsystem, which arose in response to the electric field generated by QM electron density, QM nuclei, and MM point charges. The present QM/MM EDA revealed that the electrostatic term constituted the largest stabilizing interaction between the QM and MM subsystems. When proper adjustment was made for the point charges of the MM atoms located at the QM−MM boundary, EE-QM/MM and PE-QM/MM calculations yielded similar QM electron density distributions, indicating that the MM polarization effect does not have a significant influence on the extent of QM polarization in this particular enzyme system.
■
INTRODUCTION Ubiquitously distributed cytochrome P450 enzymes (P450s) have the remarkable ability to catalyze a wide range of oxidation reactions that are involved in physiologically essential processes, such as the metabolism of organic compounds and the biosynthesis of steroid hormones.1−5 Over the past few decades, the reactive species and the metabolic reactions of P450s have garnered considerable interest. It is now widely believed that an oxoiron(IV) porphyrin π-cation radical intermediate, called compound I (Cpd I), is the reactive species that is responsible for the powerful catalytic activity of P450s. Although trapping the Cpd I intermediate in the normal turnover condition may still be a difficult challenge, accumulating experimental evidence supports the role of Cpd I as the reactive species.6−9 Computational studies using density functional theory (DFT) and hybrid quantum mechanics and molecular mechanics (QM/MM) methods have also played invaluable roles in advancing our understanding of the properties of elusive Cpd I.10−22 In addition to providing insights into the intrinsic reactivity of Cpd I, previous theoretical studies, particularly QM/MM studies, demonstrated that the stability and reactivity of Cpd I are significantly influenced by the protein environment. In recent years, we have attempted to analyze the effects of protein environment systematically.23−25 In general, the interior © 2014 American Chemical Society
environment of a protein is highly complicated because of the large number of surrounding atoms from amino acid residues. However, it is our expectation that some kind of “energy decomposition analysis (EDA)”, which was originally developed for smaller molecules,26,27 may provide intuitively interpretable descriptions of protein environments. In one study, one of the authors (H.H.) decomposed the effect of the protein environment in bacterial cytochrome P450cam (P450cam, see Figure 1) into electrostatic, van der Waals (vdW), and polarization terms.23 In so doing, the features of the mechanical-embedding (ME) and electronic-embedding (EE) schemes of the ONIOM (our own n-layered integrated molecular orbital and molecular mechanics) QM/MM method were exploited effectively.23,28,29 In another study, we used DFT-based EDA schemes to analyze the interaction between Cpd I and a water molecule,25 the importance of which had been suggested by QM/MM and DFT studies.18,19,30 Although many studies have applied QM/MM methods to enzymatic systems and their chemical reactions, almost all previous QM/MM calculations used an EE scheme of QM/ MM that takes into account the polarization of the QM Received: December 22, 2013 Revised: January 30, 2014 Published: February 2, 2014 2084
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
Article
polarizabilities for the amino acid residues in the MM subsystem and the K+ ion, whereas the QP302 polarizable force field in QuanPol was used for the crystallographic water molecules.31 The nonbonding interaction energy between the MM subsystem and the QM subsystem, ΔE, consists of electrostatic, vdW, and polarization contributions (eq 1): ΔE = Ees + Evdw + Epol
(1)
The electrostatic term (Ees) arises from the interaction of the QM electrons and nuclei with the partial charges of the MM subsystem. The effect of MM charges can be included in the Kohn−Sham one-electron operator. The DFT energy of the QM system in the gas phase is represented by the equation E0[ρ0 ] = Ts[ρ0 ] + J[ρ0 ] + E XC[ρ0 ] +
∫ v(r)ρ0 (r) dr + Enucl
(2)
where TS, J, and EXC are the kinetic energy of noninteracting electrons, the Coulomb potential energy, and the exchangecorrelation energy, respectively, and v, ρ0, and Enucl are the external potential caused by nuclei, the gas-phase electron density, and the nuclear repulsion energy, respectively. When the QM atoms are immersed in MM protein atoms without allowing the QM electron density to relax, the energy of the QM subsystem may change to E1 = E0[ρ0 (r)] −
∑∫ i
qi ·ρ0 (r) |ri − r|
dr +
∑∑ i
A
qi ·ZA |ri − rA| (3)
Figure 1. P450cam in the Cpd I state from the Protein Data Bank (PDB) (PDB code 1DZ9):7 (a) the entire enzyme; (b) a closeup view of the active site.
where qi and ri are the atomic charge and the position of MM atom i, respectively, r is the position of an electron, and ZA and rA are the charge and the position of nucleus A, respectively. Because the gas-phase QM density was used here, E1 is regarded as the QM energy of ME-QM/MM. The pure electrostatic effect on the QM subsystem by the MM subsystem, Ees, can be evaluated as the difference between E1 and E0
electron density in the presence of fixed atomic charges of classically described protein atoms. This is particularly true for studies of P450s. However, because the protein environment will be polarized by other atoms in reality, it would be interesting to perform QM/MM and EDA calculations for a P450 enzyme while taking into account the polarization effect on the MM side. In this paper, we performed such calculations for P450cam Cpd I (Figure 1) using polarizable-embedding (PE) QM/MM calculations. ME- and EE-QM/MM were also employed in combination to derive insights into the environmental effect within this enzyme.
Ees = E1 − E0
(4)
The vdW interaction between the QM and MM subsystems (EvdW) is expressed in the same way as that presented in ref 23. The MM partial charges would polarize the density from ρ0 to ρ1, and accordingly, the energy would change to E2. E2 can be obtained from iterative EE-QM/MM self-consistent field (SCF) calculations
■
THEORETICAL METHODS The geometry of P450cam in the doublet spin state, which was previously optimized using the ONIOM(DFT:MM) method, was used for the current QM/MM study.23 Units comprising the QM subsystem were also chosen in the same way: namely, protoporphyrin IX, oxoiron(IV), deprotonated Cys357, and camphor. To describe the polarization of the protein surroundings in PE-QM/MM calculations, the QuanPol polarizable force field program was used,31 which was recently implemented by Thellamurege et al. in the GAMESS quantum chemistry package.32,33 The unrestricted B3LYP method and the 6-31G* basis set were used for QM calculations. The remaining amino acid residues of the enzyme, water molecules, and the K+ ion were described by a polarizable force field. The AMBER9934 force field was used in combination with atomic
E2 = E0[ρ1(r)] −
∑∫ i
qi ·ρ1(r) |ri − r|
dr +
∑∑ i
A
qi ·ZA |ri − rA| (5)
The difference between E2 and E1 gives the polarization energy of the QM subsystem (Epol,EE) Epol,EE = E2 − E1
(6)
Several authors have pointed out the importance of using polarizable force fields in QM/MM calculations for the description of biological molecules.35−46 The current study used the induced dipole method implemented in QuanPol to describe the polarization effects of the MM atoms. This method utilizes polarizabilities placed on atomic sites, thereby allowing 2085
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
Article
the induction of dipoles.31,47,48 A dipole induced at a point depends on the polarizability and the electric field felt at that point, i.e.
d i = αi Fi
Scheme 1. Flowchart of Our QM/MM EDA Scheme
(7)
Here, di is the induced dipole vector generated at point i, Fi is the total field felt at point i, and αi is the nine-component polarizability tensor placed at point i ⎛ αi , xx αi , xy αi , xz ⎞ ⎜ ⎟ αi = ⎜ αi , yx αi , yy αi , yz ⎟ ⎜α ⎟ ⎝ i , zx αi , zy αi , zz ⎠
(8)
The polarization energy associated with the induced dipoles can be calculated using eq 9, as the interaction between the induced dipoles and the electric field
1 E id = − FT ·d 2
(9)
An H-link atom was appended at the QM-MM covalent boundary to satisfy the valency of the QM subsystem. As Leu356 and Leu358 are the amino acid residues that are directly bonded to Cys357 in the QM subsystem, the carbonyl C of Leu356 and the amide N of Leu358 were replaced by H atoms. Subsequently, the N−H and C−H bond lengths were shortened to 0.947 and 1.063 Å from the original N−C and C− N bond lengths (1.353 and 1.352 Å, respectively) to prevent instability of energy calculation. The distance scaling adopted here is the same as that used in ref 23. An MM point charge potentially gives rise to an electrostatic interaction that is too large; therefore, it may overpolarize the electron density of the QM subsystem, especially when it is located very close to the QM subsystem.46 One common approach that is used to circumvent this problem is the adjustment of the charges of MM atoms at the QM-MM boundary. In this study, four slightly different methods of charge adjustment (methods 1−4) were examined and compared. In method 1, the original MM charges were used without making any adjustment. In method 2, the charges of MM atoms located up to three bonds away from the QM atom at the boundary were eliminated, and the eliminated charges were then redistributed evenly to the adjacent atoms within the same amino acid residue. For example, as shown in Scheme 2a, the point charges of carbonyl C, carbonyl O, Cα, amide N, Cβ, and Hα of Leu356 were made zero. The total charge eliminated was −0.5047e; thus, a charge of −0.1262e was added to each of the four adjacent atoms, i.e., H, two Hβ, and Cγ. The charge redistribution allowed the total charge of a given residue to be conserved even after it was subjected to the charge adjustment. In method 3, however, this charge redistribution step was omitted. In method 4, only the charges of MM atoms that bind covalently to QM atoms were made zero, whereas other charges were left unchanged. Scheme 2 illustrates the MM partial charges of Leu356 and Leu358 after applying these four methods to the MM charges. Zero polarizabilities were placed on atoms of which the charges were made zero.
where d is the vector of induced dipoles of MM atoms and FT is the transpose of the total field vector at the MM atoms. When the effect of the QM atoms is included in the F vector, mutual polarization of the QM and MM subsystems can be described. As the induced dipoles are dependent on each other, they must be determined in an iterative manner. Typically, in a PE-QM/MM single point energy calculation, this iteration is executed at each SCF step. The induced dipoles polarize the QM electron density, which in turn changes the magnitudes of the induced dipoles of the MM atoms. These two types of iterative calculations are continued until self-consistency is fulfilled. We employed a damping function, as implemented in QuanPol, to scale down the polarizability of an MM atom that was located within a predefined distance of 3.0 Å from any of the QM atoms (see the Supporting Information).31 This polarization energy may also be calculated without including any contributions of QM atoms to F. Using the energy thus calculated (Eid,MM) and Eid, the change in polarization energy, which arises as a result of the interaction between the QM and MM subsystems, Epol,id, is calculated as Epol,id = E id − E id,MM
(10)
The PE-QM/MM calculation yields the QM density ρ2, and the QM energy E3 can be obtained by replacing ρ1 with ρ2 in eq 5. The polarization energy of the QM subsystem for the PE-QM/ MM calculation (Epol,PE) can be obtained as the difference between E3 and E1. Epol,PE = E3 − E1
(11)
The procedure used for these energy decomposition calculations is summarized in Scheme 1. To calculate the QM-MM electrostatic energy, a previous study used the classical Coulomb’s law formula.23 Although that procedure was in line with the treatment of electrostatic interaction energy in the ONIOM-ME method, it required the determination of gas-phase point charges of the QM atoms in a manner that reproduces the electrostatic potential. To skip this process and simplify the calculation of the electrostatic energy, instead we used eqs 2−4 to calculate Ees. Thus, the electrostatic interaction was calculated as the interaction between the point charges in the MM subsystem and the electron density and nuclei in the QM subsystem.
■
RESULTS AND DISCUSSION Table 1 summarizes the electrostatic, polarization, and vdW interaction energies for P450cam Cpd I, which were calculated using the EDA procedure and applying the four QM/MM charge adjustment methods. 2086
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
Article
Scheme 2. MM Charges (e) on Leu356 and Leu358 after Applying the QM/MM-Boundary Charge Redistribution Methods: (a) Method 1, (b) Method 2, (c) Method 3, and (d) Method 4a
a
Adjusted charges are colored blue.
Epol,EE was calculated as −105.5 kcal/mol, which was close to the previous value. Figure 2a illustrates how methods 1 and 2 gave different Mulliken charge distributions in the QM subsystem. One can see clearly in the figure that the charges of Cys357, which is bonded directly to Leu356 and Leu358 in the MM subsystem (Scheme 2), exhibited large differences. The overpolarization also affected the spin distribution in the QM subsystem (Figure 2b). A particularly large difference in spin population was seen for the sulfur atom in Cys357. Modest differences were observed for the atoms of the porphyrin ring, despite the fact that the porphyrin was not very close to Leu356 and Leu358. This result suggests a significant influence of the overpolarization at the sulfur atom on the entire singly occupied a2u-type orbital, which has significant amplitudes on the porphyrin as well.14,15 The electrostatic interaction was diminished by 108−123 kcal/mol, when three different charge adjustment methods (2− 4) were applied (see Table 1). Nevertheless, consistent with a previous conclusion,23 the role played by the electrostatic effect in the nonbonding interaction between the QM and MM subsystems was the largest. A comparison between the results obtained using methods 2 and 3 showed that the electrostatic stabilization and the polarization were slightly larger when the eliminated charges were not redistributed to neighboring atoms. In contrast, method 4, which simply eliminates the charges of only the atoms at the QM-MM boundary, gave electrostatic and polarization stabilization effects that are slightly smaller than those obtained using method 2. Overall, the three methods yielded more or less similar energy values, and the values were not very sensitive to the charge adjustment method employed. Therefore, henceforth, we will mainly use the results obtained using method 2 for further discussions.
Table 1. Decomposed QM/MM Interaction Energy Terms (kcal/mol) for P450cam Cpd I Obtained Using Various Charge-adjustment Methods term
method 1
method 2
method 3
method 4
Ees EvdW Epol,EE Epol,PE Epol,id
−331.3 −105.8 −294.1 −293.4 +62.1
−214.4 −105.8 −105.5 −105.0 +65.8
−222.9 −105.8 −112.7 −112.1 +69.1
−208.7 −105.8 −99.4 −99.0 +59.6
The electrostatic energy was calculated as −331.3 kcal/mol using method 1. This value agrees reasonably well with the value obtained previously (−289.9 kcal/mol).23 The difference between these values can be attributed partly to the different DFT methods and MM force fields used. A more important reason may be that the classical Coulomb’s law formula was used previously for the calculation of Ees; thus, the pairwise interaction energy between MM and QM point charges was scaled down when two atoms were located not more than three bonds away from each other, whereas this sort of scaling was not applied in the present approach. The polarization stabilization of the QM density (Epol,EE: −294.1 kcal/mol) obtained using method 1 was much larger than the stabilization estimated previously (−101.6 kcal/mol). This significant discrepancy arose because the QM subsystem feels MM charges in the immediate proximity (Scheme 2a) in the current EE-QM/MM calculation using method 1, whereas the effect of a few point charges at the QM-MM covalent boundary on the QM density are ignored in the ONIOM-EE calculation.28 To avoid the overpolarization observed in method 1, some type of charge adjustment was deemed appropriate. In fact, the overpolarization was mitigated significantly using method 2: 2087
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
Article
Figure 2. Deviation of the Mulliken charge (a) and Mulliken spin population (b) values of the QM atoms obtained using method 1 from those obtained using method 2. See Scheme S1 for the specific atom numbering.
Since Ees is additive, it can be further decomposed into contributions from individual residues as presented in Figure 3. Such decomposition can be performed by repeating the calculation according to eq 4, while assigning charges of only one residue at a time. The calculated energy contributions were very similar to those reported in ref 23, despite the fact that the DFT methods, force fields, and Ees formulas used were different in these two studies. For example, Arg112, Arg299, and His355, which are located in the close proximity of the carboxylate groups of the protoporphyrin ligand (see Figure 1b), showed large stabilizing electrostatic effects. Tyr96, which forms a
hydrogen bond with the carbonyl group of camphor, had an Ees value of −9.5 kcal/mol, which was comparable with the value calculated in ref 23 (−8.9 kcal/mol). However, large differences were observed for Leu356 (see also Figure S1). This discrepancy arose because, in the current approach, a zeroed charge (e.g., that of the carbonyl carbon of Leu356) in method 2 could never have any electrostatic interaction with the QM subsystem. In contrast, in the previous classical approach, the same carbonyl carbon, for which no MM charge adjustment was made, could have a nonzero pairwise electrostatic interaction with atoms in the QM subsystem. This explains 2088
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
Article
Figure 3. Comparison of electrostatic energy contributions from amino acid residues 10−199 (a) and 200−414 and the potassium ion (#415) (b) obtained in this work with method 2 (blue) and in a previous work (red).23
why the Ees value between the current approach (with method 2, −214.4 kcal/mol) and the previous approach (−289.9 kcal/ mol) differed by ∼76 kcal/mol. The calculated vdW interaction energy was equal to the value obtained previously because the same method was used; thus, a further explanation is not provided here. Table 1 shows that Epol,PE was smaller in magnitude than Epol,EE by 0.5 kcal/mol, which suggests that the polarization of the QM subsystem was slightly suppressed in the PE-QM/MM calculation. However, the difference between these polarization
energies was fairly small; thus, the MM polarization seems to have almost no impact on the QM electron density. Even in the absence of the QM atoms, the MM atoms can have induced dipoles because of the existing electric fields generated by many MM atoms. The stabilization gained by the MM atoms (Eid,MM) amounted to >2900 kcal/mol (Table S2). The polarization energy caused by induced dipoles will change when QM atoms come into the active site. This differential effect of Eid can be quantified by Epol,id (see eq 10); the resulting value, +65.8 kcal/mol (with method 2), was positive. Therefore, it appears that the inclusion of the QM atoms 2089
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
Article
to 69.6 D in the PE-QM/MM calculation. This similarity in dipole moment provides another indication that the extent of polarization of the QM density is nearly equal in the EE-QM/ MM and PE-QM/MM calculations. We reported above that method 1 caused overpolarization. Using method 1, the dipole moment was calculated as 70.3 and 70.5 D with EE-QM/MM and PE-QM/MM, respectively. The values were somewhat larger than those obtained using method 2 (Table 2) because of the strong polarization effect exerted on the QM atoms for method 1.
somewhat diminished the polarization stabilization associated with induced dipoles. The key atomic spin populations obtained from different calculations are compared in Figure 4 (see also Figure S2 for
■
CONCLUSIONS In this work, we performed a QM/MM study of cytochrome P450cam in the Cpd I state to quantify the protein environmental effect of the enzyme. The B3LYP method was used for the QM subsystem, whereas AMBER99 and QP302 were used for the MM subsystem. ME-, EE-, and PE-QM/MM methods were applied in a stepwise manner to decompose the nonbonding QM−MM interaction energy into electrostatic, vdW, and polarization contributions. In the absence of adjustment for the atomic charges of the MM atoms, an overpolarization of the QM density was observed because the MM charges at the QM−MM boundary interacted with the QM atoms too strongly. However, this overpolarization effect was alleviated effectively by adopting charge-adjustment methods. The QM/MM energy decomposition analysis showed that the electrostatic interaction was the major factor that stabilized the binding of the QM subsystem to the protein. The polarization and vdW effects also contributed significantly to the binding process. However, the polarization stabilization gained by the MM atoms from induced dipoles was reduced when the protein accommodated the QM atoms. The QM polarization energies, distributions of charge and spin, and dipole moments all suggest that the electron density of the QM subsystem of this enzyme is polarized to a similar extent in EEQM/MM and PE-QM/MM calculations.
Figure 4. Key Mulliken atomic spin populations obtained with MEQM/MM (gas phase), EE-QM/MM, and PE-QM/MM. Some of the atom labels were taken from the PDB file (PDB code 1DZ9).
full details of the spin and charge distribution). As observed previously in ref 15, there was noticeable and unnatural spin localization at the propionate oxygen atoms in the gas-phase or ME-QM/MM calculation. Consequently, no significant spin population was seen on the porphyrin ring. In contrast, the spin populations on propionates, which interact with arginine residues (see Figure 1b), were shifted back to the porphyrin moiety in the EE-QM/MM calculation.15 As this shifted spin originates from the unpaired electron occupying the a2u-type orbital,11 significant spin populations were observed on the meso carbon and pyrrole nitrogen atoms of the porphyrin. The PE-QM/MM calculation produced a spin distribution that was similar to that of the EE-QM/MM calculation. This result indicates that EE-QM/MM and PE-QM/MM yielded similarly polarized QM electron densities. Table 2 lists the calculated dipole moments of the QM subsystem. Although the QM subsystem prefers a lower polarization in the gas-phase or ME-QM/MM calculation, the polarization was enhanced in the EE-QM/MM calculation because the polar nature of the protein environment was taken into account. As a result, the dipole moment increased from 27.3 to 69.9 D. The dipole moment was only slightly reduced
■
S Supporting Information *
QM energies at each step of Scheme 1, Eid and Eid,MM values, differences in Ees of a few residues between the current and previous work, differences in Mulliken charge and spin distributions between the ME- and EE-/PE-QM/MM methods, numbering of QM atoms, spin distributions obtained with a few more residues included in the QM region, and the damping function employed. This material is available free of charge via the Internet at http://pubs.acs.org.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail: [email protected] (H.H.). Notes
The authors declare no competing financial interest.
Table 2. Comparison of the Calculated Dipole Moments (D) of the QM Subsystem Obtained with Method 2
■
ACKNOWLEDGMENTS This work was supported by a Nanyang Assistant Professorship and an AcRF Tier 1 grant (RG14/12). We thank the High Performance Computing Centre at Nanyang Technological University for computer resources.
dipole moment (D) ME-QM/MM EE-QM/MM PE-QM/MM
ASSOCIATED CONTENT
27.3 69.9 69.6 2090
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
■
Article
Cytochrome c Pperoxidase, and Ascorbate Peroxidase. J. Comput. Chem. 2006, 27, 1352−1362. (21) Porro, C. S.; Sutcliffe, M. J.; de Visser, S. P. Quantum Mechanics/Molecular Mechanics Studies on the Sulfoxidation of Dimethyl Sulfide by Compound I and Compound 0 of Cytochrome P450: Which is the Better Oxidant? J. Phys. Chem. A 2009, 113, 11635−11642. (22) de Visser, S. P.; Ogliaro, F.; Sharma, P. K.; Shaik, S. What Factors Affect the Regioselectivity of Oxidation by Cytochrome P450? A DFT Study of Allylic Hydroxylation and Double Bond Epoxidation in a Model Reaction. J. Am. Chem. Soc. 2002, 124, 11809−11826. (23) Hirao, H. Energy Decomposition Analysis of the Protein Environmental Effect: the Case of Cytochrome P450cam Compound I. Chem. Lett. 2011, 40, 1179−1181. (24) Hirao, H. The Effects of Protein Environment and Dispersion on the Formation of Ferric-Superoxide Species in Myo-Inositol Oxygenase (MIOX): A Combined ONIOM(DFT:MM) and Energy Decomposition Analysis. J. Phys. Chem. B 2011, 115, 11278−11285. (25) Thellamurege, N.; Hirao, H. Water Complexes of Cytochrome P450: Insights from Energy Decomposition Analysis. Molecules 2013, 18, 6782−6791. (26) Kitaura, K.; Morokuma, K. A New Energy Decomposition Scheme for Molecular Interactions within the Hartree-Fock Approximation. Int. J. Quantum Chem. 1976, 18, 1558−1565. (27) von Hopffgarten, M.; Frenking, G. Energy Decomposition Analysis. WIREs Comput. Mol. Sci. 2012, 2, 43−62. (28) Vreven, T.; Byun, K. S.; Komaromi, I.; Dapprich, S.; Montgomry, J. A., Jr.; Morokuma, K.; Frisch, M. J. Combining Quantum Mechanics Methods with Molecular Mechanics Methods in ONIOM. J. Chem. Theory. Comput. 2006, 2, 815−826. (29) Chung, L. W.; Hirao, H.; Li, X.; Morokuma, K. The ONIOM Method: Its Foundation and Application to Metalloenzymes and Photobiology. WIREs Comput. Mol. Sci. 2012, 2, 327−350. (30) Kumar, D.; Altun, A.; Shaik, S.; Thiel, W. Water as Biocatalyst in Cytochrome P450. Faraday Discuss. 2011, 148, 373−383. (31) Thellamurege, N. M.; Si, D.; Cui, F.; Zhu, H.; Lai, R.; Li, H. QuanPol: a Full Spectrum and Seamless QM/MM Program. J. Comput. Chem. 2013, 34, 2816−2833. (32) Gordon, M. S.; Schmidt, M. W. Theory and Applications of Computational Chemistry; Dykstra, C. E., Frenking, K. S., Scuseria, G. E., Eds.; Elsevier: Amsterdam, The Netherlands, 2005. (33) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347−1363. (34) Wang, J.; Cieplak, P.; Kollman, P. A. How Well Does a Restrained Electrostatic Potential (RESP) Model Perform in Calculating Conformational Energies of Organic and Biological Molecules? J. Comput. Chem. 2000, 21, 1049−1074. (35) Warshel, A.; Levitt, M. Theoretical Studies of Enzymic Reactions: Dielectric, Electrostatic and Steric Stabilization of the Carbonium Ion in the Reaction of Lysozyme. J. Mol. Biol. 1976, 103, 227−249. (36) Russell, S. T.; Warshel, A. Calculations of Electrostatic Energies in Proteins: The Energetics of Ionized Groups in Bovine Pancreatic Trypsin Inhibitor. J. Mol. Biol. 1985, 185, 389−404. (37) Warshel, A.; Sharma, P. K.; Kato, M.; Parson, W. W. Modeling Electrostatic Effects in Proteins. Biochim. Biophys. Acta 2006, 1764, 1647−1676. (38) Thompson, M. A.; Schenter, G. K. Excited States of the Bacteriochlorophyll b Dimer of Rhodopseudomonas viridis: A QM/ MM Study of the Photosynthetic Reaction Center That Includes MM Polarization. J. Phys. Chem. 1995, 99, 6374−6386. (39) Illingworth, C. J. R.; Gooding, S. R.; Winn, P. J.; Jones, G. A.; Ferenczy, G. G.; Reynolds, C. A. Classical Polarization in Hybrid QM/ MM Methods. J. Phys. Chem. A 2006, 110, 6487−6497. (40) Zhang, Y.; Lin, H.; Truhlar, D. G. Self-Consistent Polarization of the Boundary in the Redistributed Charge and Dipole Scheme for
REFERENCES
(1) Omura, T.; Sato, R. A New Cytochrome in Liver Microsomes. J. Biol. Chem. 1962, 237, PC1375−PC1376. (2) Denisov, I. G.; Makris, T. M.; Sligar, S. G.; Schlichting, I. Structure and Chemistry of Cytochrome P450. Chem. Rev. 2005, 105, 2253−2278. (3) Ortiz de Montellano, P. R. Cytochrome P450: Structure Mechanism and Biochemistry, 3rd ed.; Ortiz de Montellano, P. R., Ed.; Kluwer Academic/Plenum Publishers: New York, 2005. (4) Sono, M.; Roach, M. P.; Coulter, E. D.; Dawson, J. H. Heme Containing Oxygenases. Chem. Rev. 1996, 96, 2841−2887. (5) Guengerich, F. P. Cytochrome P450 and Chemical Toxicology. Chem. Res. Toxicol. 2008, 21, 70−83. (6) Groves, J. T. High-Valent Iron in Chemical and Biological Oxidations. J. Inorg. Biochem. 2006, 100, 434−447. (7) Schlichting, I.; Berendzen, J.; Chu, K.; Stock, A. M.; Maves, S. A.; Benson, D. E.; Sweet, R. M.; Ringe, D.; Petsko, G. A.; Sligar, S. G. The Catalytic Pathway of Cytochrome P450cam at Atomic Resolution. Science 2000, 287, 1615−1622. (8) Davydov, R.; Makris, T. M.; Kofman, V.; Werst, D. E.; Sligar, S. G.; Hoffman, B. M. Hydroxylation of Camphor by Oxy-cytochrome P450cam: Mechanistic Implications of EPR and ENDOR Studies of Catalytic Intermediates in Native and Mutant Enzymes. J. Am. Chem. Soc. 2001, 123, 1403−1415. (9) Rittle, J.; Green, M. T. Cytochrome P450 Compund I: Capture, Characterization, and C-H Bond Activation Kinetics. Science 2010, 330, 933−937. (10) Meunier, B.; de Visser, S. P.; Shaik, S. Mechanism of Oxidation Reactions Catalyzed by Cytochrome P450 Enzymes. Chem. Rev. 2004, 104, 3947−3980. (11) Shaik, S.; Kumar, D.; de Visser, S. P.; Altun, A.; Thiel, W. Theoretical Perspective on the Structure and Mechanism of the Cytochrome P450 Enzymes. Chem. Rev. 2005, 105, 2279−2328. (12) Shaik, S.; Hirao, H.; Kumar, D. Reactivity Patterns of Cytochrome P450 Enzymes: Multifunctionality of the Active Species, and Two States-Two Oxidants Conundrum. Nat. Prod. Rep. 2007, 24, 533−552. (13) Shaik, S.; Cohen, S.; Wang, Y.; Chen, H.; Kumar, D.; Thiel, W. P450 Enzymes; Their Structure Reactivity and Selectivity − Modeled by QM/MM Calculations. Chem. Rev. 2010, 110, 949−1017. (14) Schöneboom, J. C.; Lin, H.; Reuter, N.; Thiel, W.; Cohen, S.; Ogliaro, F.; Shaik, S. The Elusive Oxidant Species of Cytochrome P450 Enzymes: Characterization by Combined Quantum Mechanical/ Molecular Mechanical (QM/MM) Calculations. J. Am. Chem. Soc. 2002, 124, 8142−8151. (15) Schöneboom, J. C.; Cohen, S.; Lin, H.; Shaik, S.; Thiel, W. Quantum Mechanical/Molecular Mechanical Investigation of the Mechanism of C-H Hydroxylation of Camphor by Cytochrome P450cam: Theory Supports a Two-State Rebound Mechanism. J. Am. Chem. Soc. 2004, 126, 4017−4034. (16) Cho, K.-B.; Hirao, H.; Chen, H.; Carvajal, M. A.; Cohen, S.; Derat, E.; Thiel, W.; Shaik, S. Compound I in the Heme Thiolate Enzymes: A Comparative QM/MM Study. J. Phys. Chem. A 2008, 112, 13128−13138. (17) Schöneboom, J. C.; Neese, F.; Thiel, W. Toward Identification of the Compound I Reactive Intermediate in Cytochrome P450 Chemistry: A QM/MM Study of Its EPR and Mossbaur Parameters. J. Am. Chem. Soc. 2005, 127, 5840−5853. (18) Altun, A.; Guallar, V.; Friesner, R. A.; Shaik, S.; Thiel, W. The Effect of Heme Environment on the Hydrogen Abstraction Reaction of Camphor in P450cam Catalysis: A QM/MM Study. J. Am. Chem. Soc. 2006, 128, 3924−3925. (19) Altun, A.; Shaik, S.; Thiel, W. What is the Active Species of Cytochrome P450 during Camphor Hydroxylation? QM/MM Studies of the Different Electronic States of Compound I and of Reduced and Oxidized Iron-Oxo Intermediates. J. Am. Chem. Soc. 2007, 129, 8978− 8987. (20) Harvey, J. N.; Bathelt, C. M.; Mulholland, A. J. QM/MM Modeling of Compound I Active Species in Cytochrome P450, 2091
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
The Journal of Physical Chemistry B
Article
Combined Quantum-Mechanical and Molecular-Mechanical Calculations. J. Chem. Theory Comput. 2007, 3, 1378−1398. (41) Zhang, Y.; Lin, H. Flexible-Boundary QM/MM Calculations: II. Partial Charge Transfer across the QM/MM Boundary That Passes Through a Covalent Bond. Theor. Chem. Acc. 2010, 126, 315−322. (42) Wanko, M.; Hoffmann, M.; Frähmcke, J.; Frauenheim, T.; Elstner, M. Effect of Polarization on the Opsin Shift in Rhodopsins. 2. Empirical Polarization Models for Proteins. J. Phys. Chem. B 2008, 112, 11468−11478. (43) Frähmcke, J. S.; Wanko, M.; Phatak, P.; Mroginski, M. A.; Elstner, M. The Protonation State of Glu181 in Rhodopsin Revisited: Interpretation of Experimental Data on the Basis of QM/MM Calculations. J. Phys. Chem. B 2010, 114, 11338−11352. (44) Welke, K.; Frähmcke, J. S.; Watanabe, H. C.; Hegemann, P.; Elstner, M. Color Tuning in Binding Pocket Models of the Chlamydomonas-Type Channelrhodopsins. J. Phys. Chem. B 2011, 115, 15119−15128. (45) Wolter, T.; Welke, K.; Phatak, P.; Bondar, A. N.; Elstner, M. Excitation Energies of a Water-Bridged Twisted Retinal Structure in the Bacteriorhodopsin Proton Pump: A Theoretical Intvestigation. Phys. Chem. Chem. Phys. 2013, 15, 12582−12590. (46) Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem., Int. Ed. 2009, 48, 1198−1229. (47) Li, H. Quantum Mechanical/Molecular Mechanical/Continuum Style Solvation Model: Linear Response Theory, Variational Treatment and Nuclear Gradients. J. Chem. Phys. 2009, 131, 184103. (48) Li, H.; Gordon, M. S. Gradients of the Polarization Energy in the Effective Fragment Potential Method. J. Chem. Phys. 2006, 125, 194103−194112.
2092
dx.doi.org/10.1021/jp412538n | J. Phys. Chem. B 2014, 118, 2084−2092
