Quantum Mechanics/Molecular Mechanics Studies on the Relative Reactivities of Compound I and II in Cytochrome P450 Enzymes

The cytochromes P450 are drug metabolizing enzymes in the body that typically react with substrates through a monoxygenation reaction. During the catalytic cycle two reduction and protonation steps generate a high-valent iron (IV)-oxo heme cation radical species called Compound I. However, with sufficient reduction equivalents present, the catalytic cycle should be able to continue to the reduced species of Compound I, called Compound II, rather than a reaction of Compound I with substrate. In particular, since electron transfer is usually on faster timescales than atom transfer, we considered this process feasible and decided to investigate the reaction computationally. In this work we present a computational study using density functional theory methods on active site model complexes alongside quantum mechanics/molecular mechanics calculations on full enzyme structures of cytochrome P450 enzymes. Specifically, we focus on the relative reactivity of Compound I and II with a model substrate for O–H bond activation. We show that generally the barrier heights for hydrogen atom abstraction are higher in energy for Compound II than Compound I for O–H bond activation. Nevertheless, for the activation of such bonds, Compound II should still be an active oxidant under enzymatic conditions. As such, our computational modelling predicts that under high-reduction environments the cytochromes P450 can react with substrates via Compound II but the rates will be much slower.


Introduction
The cytochromes P450 are important heme monoxygenases in the body and particularly prevalent in the liver, where they catalyze the biosynthesis of several hormones including estrogen and also take part in the biodegradation and metabolism of toxic compounds [1][2][3][4][5][6][7][8][9]. In general, the P450s bind molecular oxygen on an iron(III)-heme group and transfer one of its oxygen atoms to a substrate through either an aliphatic C-H hydroxylation, aromatic C-H hydroxylation, C=C epoxidation or The consensus catalytic cycle of P450 enzymes is schematically depicted in Figure 2 [20][21][22]. The cycle starts from an iron (III)-water bound heme structure, the resting state, that upon substrate (SubH) binding in the substrate binding pocket releases the water molecule. Then, iron (III) (heme)-  Despite the fact that Cpd I has been detected and characterized for one specific P450 isozyme, there are still on-going discussions on what the actual reactive species in P450 enzymes is and as such, a multiple oxidant hypothesis has been proposed [28]. These studies initially pointed to Cpd 0 as a possible "second oxidant", but a series of computational studies by model complexes [29,30] as well as using quantum mechanics/molecular mechanics (QM/MM) calculations [31] ruled this hypothesis out. Subsequent experimental studies on biomimetic model complexes indeed confirmed the computational prediction and found sluggish reactivity by Cpd 0 as compared to Cpd I [32].
Since then research has moved to alternative oxidants that could act as a "second oxidant" in P450 catalysis. In this respect, van Eldik and co-workers used biomimetic models of Cpd I and Cpd 0 as well as the one-electron reduced form of Cpd I, designated Cpd II [33][34][35]. They studied different oxidative processes (epoxidation, sulfoxidation, hydroxylation, O-H hydrogen atom abstraction, and hydride-transfer reactions) of Cpd I, Cpd II and Cpd 0 toward a variety of substrates. Interestingly, it was found that Cpd II exhibits competent reactivity in O-H hydrogen atom abstraction reactions, whereas it is the most efficient oxidant in hydride-transfer processes [33]. A recent computational study from our groups reasoned that an additional reduction step of Cpd I to give Cpd II would keep the oxidant active in hydrogen atom abstraction reactions, although probably with lesser oxidative power than Cpd I [36]. Furthermore, due to the rising interest in the reactivity of Cpd II, other computational studies that compare the reactivity of Cpd I and Cpd II mimics in other oxidation processes (aldehyde oxidation and alcohol oxidation) [37] or that exclusively study the reactivity pattern of Cpd II [38] have also been reported. In this work, we expand on the studies reported in Ref [36]   Despite the fact that Cpd I has been detected and characterized for one specific P450 isozyme, there are still on-going discussions on what the actual reactive species in P450 enzymes is and as such, a multiple oxidant hypothesis has been proposed [28]. These studies initially pointed to Cpd 0 as a possible "second oxidant", but a series of computational studies by model complexes [29,30] as well as using quantum mechanics/molecular mechanics (QM/MM) calculations [31] ruled this hypothesis out. Subsequent experimental studies on biomimetic model complexes indeed confirmed the computational prediction and found sluggish reactivity by Cpd 0 as compared to Cpd I [32].
Since then research has moved to alternative oxidants that could act as a "second oxidant" in P450 catalysis. In this respect, van Eldik and co-workers used biomimetic models of Cpd I and Cpd 0 as well as the one-electron reduced form of Cpd I, designated Cpd II [33][34][35]. They studied different oxidative processes (epoxidation, sulfoxidation, hydroxylation, O-H hydrogen atom abstraction, and hydride-transfer reactions) of Cpd I, Cpd II and Cpd 0 toward a variety of substrates. Interestingly, it was found that Cpd II exhibits competent reactivity in O-H hydrogen atom abstraction reactions, whereas it is the most efficient oxidant in hydride-transfer processes [33]. A recent computational study from our groups reasoned that an additional reduction step of Cpd I to give Cpd II would keep the oxidant active in hydrogen atom abstraction reactions, although probably with lesser oxidative power than Cpd I [36]. Furthermore, due to the rising interest in the reactivity of Cpd II, other computational studies that compare the reactivity of Cpd I and Cpd II mimics in other oxidation processes (aldehyde oxidation and alcohol oxidation) [37] or that exclusively study the reactivity pattern of Cpd II [38] have also been reported. In this work, we expand on the studies reported in Ref [36] and investigate whether Cpd II could be a viable oxidant of O-H hydrogen atom abstraction reactions. In particular, we combine calculations on synthetic model complexes with QM/MM studies on the full enzyme that take the effect of the protein into consideration. Specifically, we compare the reactivity of Cpd I and Cpd II as oxidants.

Results
To understand the reactivity differences of Cpd I and Cpd II, we decided to do a computational study and used two approaches, namely (1) density functional theory (DFT) studies on enzyme active site models and (2) quantum mechanics/molecular mechanics (QM/MM) calculations on the full protein with solvent layer. Figure 3 displays the calculated reaction mechanism and the definition of the individual structures. As previous work [36] showed that Cpd II may have difficulties with radical rebound steps, we decided to use TEMPOH (2,2,6,6-tetramethyl-piperidine-1-ol) as a substrate as it is a one-electron transfer substrate only. The reaction starts with isolated reactants, i.e., Cpd I + TEMPOH (Re CpdI ) or Cpd II + TEMPOH (Re CpdII ), and proceeds via hydrogen atom abstraction via transition state TS HA to form a radical intermediate Int, which is an iron(IV)-hydroxo heme with TEMPO • or iron (III)-hydroxo heme with TEMPO • . In a subsequent step two TEMPO • radicals pair up to form a dimer. on the full enzyme that take the effect of the protein into consideration. Specifically, we compare the reactivity of Cpd I and Cpd II as oxidants.

Results
To understand the reactivity differences of Cpd I and Cpd II, we decided to do a computational study and used two approaches, namely (1) density functional theory (DFT) studies on enzyme active site models and (2) quantum mechanics/molecular mechanics (QM/MM) calculations on the full protein with solvent layer. Figure 3 displays the calculated reaction mechanism and the definition of the individual structures. As previous work [36] showed that Cpd II may have difficulties with radical rebound steps, we decided to use TEMPOH (2,2,6,6-tetramethyl-piperidine-1-ol) as a substrate as it is a one-electron transfer substrate only. The reaction starts with isolated reactants, i.e., Cpd I + TEMPOH (ReCpdI) or Cpd II + TEMPOH (ReCpdII), and proceeds via hydrogen atom abstraction via transition state TSHA to form a radical intermediate Int, which is an iron(IV)-hydroxo heme with TEMPO • or iron (III)-hydroxo heme with TEMPO • . In a subsequent step two TEMPO • radicals pair up to form a dimer. We initially studied the reaction described in Figure 3 with a gas-phase DFT model (Section 2.1) but followed it up with a QM/MM study (Section 2.2) that takes the shape and size of the protein into consideration.

DFT on Model Complexes
We started with a small gas-phase DFT model that contains an iron (IV)-oxo porphyrin cation radical model used previously [39][40][41] that has all substituents on the porphyrin replaced by hydrogen atoms and an axial thiolate ligand. We calculated the mechanism on the lowest lying doublet and quartet spin state surfaces. Figure 4 displays the potential energy profile for hydrogen atom abstraction by 4,2 Cpd I ( 4,2 1) or 4,2 [Fe IV (O)(Por +• )SH] from TEMPOH. As can be seen on both spin states, the hydrogen atom abstraction barrier is negligible and collapses to the iron (IV)-hydroxo product rapidly. The structure shows short TEMPO-H and long H-OFe distances and as such the transition states are early on the potential energy landscape. These distances match with previously reported hydrogen atom abstraction structures by P450 Cpd I well [42], although different substrates We initially studied the reaction described in Figure 3 with a gas-phase DFT model (Section 2.1) but followed it up with a QM/MM study (Section 2.2) that takes the shape and size of the protein into consideration.

DFT on Model Complexes
We started with a small gas-phase DFT model that contains an iron (IV)-oxo porphyrin cation radical model used previously [39][40][41] that has all substituents on the porphyrin replaced by hydrogen atoms and an axial thiolate ligand. We calculated the mechanism on the lowest lying doublet and quartet spin state surfaces. Figure 4 displays the potential energy profile for hydrogen atom abstraction by 4,2 Cpd I ( 4,2 1) or 4,2 [Fe IV (O)(Por +• )SH] from TEMPOH. As can be seen on both spin states, the hydrogen atom abstraction barrier is negligible and collapses to the iron (IV)-hydroxo product rapidly. The structure shows short TEMPO-H and long H-OFe distances and as such the transition states are early on the potential energy landscape. These distances match with previously reported hydrogen atom abstraction structures by P450 Cpd I well [42], although different substrates were investigated.  Subsequently, we calculated the same mechanism but from Cpd II ( 5,3 2) on the triplet and quintet spin states. On the triplet spin state, a small barrier is encountered, whereas it is barrierless on the somewhat higher lying quintet spin state. Therefore, both Cpd I and Cpd II should be able to easily activate TEMPOH at room temperature. The driving force is slightly less exothermic for Cpd II than for Cpd I and follows the order of the hydrogen atom abstraction barriers well. The geometries of the Cpd I and Cpd II hydrogen atom abstraction transition states are dramatically different. As mentioned above for Cpd I the transition states are early with short TEMPO-H distances; however, for the Cpd II transition states these are well larger, in particular for 3 TSHA,CpdII where we find a TEMPO-H distance of 1.329 Å as compared to a value of 1.029 Å for 2 TSHA,CpdI. At the same time the accepting O-H distance is considerably shorter for the Cpd II transition states as compared to the ones for Cpd I.  Subsequently, we calculated the same mechanism but from Cpd II ( 5,3 2) on the triplet and quintet spin states. On the triplet spin state, a small barrier is encountered, whereas it is barrierless on the somewhat higher lying quintet spin state. Therefore, both Cpd I and Cpd II should be able to easily activate TEMPOH at room temperature. The driving force is slightly less exothermic for Cpd II than for Cpd I and follows the order of the hydrogen atom abstraction barriers well. The geometries of the Cpd I and Cpd II hydrogen atom abstraction transition states are dramatically different. As mentioned above for Cpd I the transition states are early with short TEMPO-H distances; however, for the Cpd II transition states these are well larger, in particular for 3 TS HA,CpdII where we find a TEMPO-H distance of 1.329 Å as compared to a value of 1.029 Å for 2 TS HA,CpdI . At the same time the accepting O-H distance is considerably shorter for the Cpd II transition states as compared to the ones for Cpd I.
The most dramatic difference, as seen from Figures 4 and 5 is that the substrate approach on the oxo group is from the side for Cpd I, while it is from the top in Cpd II. The substrate approach previously [43][44][45][46][47] was correlated with the electron transfer process that happens in the transition state. In Cpd I the oxidant has electronic configuration of δ x2-y2 2 π* xz 1 π* yz 1 a 2u 1 , whereby the three unpaired electrons are either ferromagnetically coupled into a quartet spin state or have the heme a 2u electron with down-spin and the π* xz and π* yz electrons with up-spin into an overall doublet spin state. Figure 6 gives the orbital shapes and Cpd I and Cpd II occupation numbers. The lowest energy orbital of the set shown in Figure 6 is the δ x2-y2 orbital, which is non-bonding and in the plane of the porphyrin ring. The degenerate π* xz and π* yz molecular orbitals are higher in energy for the antibonding interactions of metal and oxo group through 3d and 2p atomic orbital components in these planes. Two virtual orbitals complete the set of metal-type orbitals, namely the σ* z2 and the σ* xy . The former is the σ* z2 antibonding orbital for the metal with oxo interaction, whereas the σ* xy orbital represents the antibonding interactions of the metal with the heme nitrogen atoms. The heme also has several high-lying π*-type orbitals and one of them, the a 2u orbital, is close in energy to the π* xz and π* yz orbitals and is singly occupied. Upon reduction of Cpd I to Cpd II, the additional electron fills the a 2u orbital with a second electron in the triplet spin state, whereas in the quintet spin state also a promotion from δ x2-y2 to σ* xy takes place. Consequently, the electronic configuration of 3 Cpd II ( 3 2) is δ x2-y2 2 π* xz 1 π* yz 1 a 2u 2 and 5 Cpd II ( 5 2) is δ x2-y2 1 π* xz 1 π* yz 1 σ* xy 1 a 2u 2 .
Subsequently, we calculated the same mechanism but from Cpd II ( 5,3 2) on the triplet and quintet spin states. On the triplet spin state, a small barrier is encountered, whereas it is barrierless on the somewhat higher lying quintet spin state. Therefore, both Cpd I and Cpd II should be able to easily activate TEMPOH at room temperature. The driving force is slightly less exothermic for Cpd II than for Cpd I and follows the order of the hydrogen atom abstraction barriers well. The geometries of the Cpd I and Cpd II hydrogen atom abstraction transition states are dramatically different. As mentioned above for Cpd I the transition states are early with short TEMPO-H distances; however, for the Cpd II transition states these are well larger, in particular for 3 TSHA,CpdII where we find a TEMPO-H distance of 1.329 Å as compared to a value of 1.029 Å for 2 TSHA,CpdI. At the same time the accepting O-H distance is considerably shorter for the Cpd II transition states as compared to the ones for Cpd I.  The most dramatic difference, as seen from Figures 4 and 5 is that the substrate approach on the oxo group is from the side for Cpd I, while it is from the top in Cpd II. The substrate approach previously [43][44][45][46][47] was correlated with the electron transfer process that happens in the transition state. In Cpd I the oxidant has electronic configuration of δx2-y2 2 π*xz 1 π*yz 1 a2u 1 , whereby the three unpaired electrons are either ferromagnetically coupled into a quartet spin state or have the heme a2u electron with down-spin and the π*xz and π*yz electrons with up-spin into an overall doublet spin state. Figure 6 gives the orbital shapes and Cpd I and Cpd II occupation numbers. The lowest energy orbital of the set shown in Figure 6 is the δx2-y2 orbital, which is non-bonding and in the plane of the porphyrin ring. The degenerate π*xz and π*yz molecular orbitals are higher in energy for the antibonding interactions of metal and oxo group through 3d and 2p atomic orbital components in these planes. Two virtual orbitals complete the set of metal-type orbitals, namely the σ*z2 and the σ*xy. The former is the σ*z2 antibonding orbital for the metal with oxo interaction, whereas the σ*xy orbital represents the antibonding interactions of the metal with the heme nitrogen atoms. The heme also has several high-lying π*-type orbitals and one of them, the a2u orbital, is close in energy to the π*xz and π*yz orbitals and is singly occupied. Upon reduction of Cpd I to Cpd II, the additional electron fills the a2u orbital with a second electron in the triplet spin state, whereas in the quintet spin state also a promotion from δx2-y2 to σ*xy takes place. Consequently, the electronic configuration of 3 Cpd II ( 3 2) is δx2-y2 2 π*xz 1 π*yz 1 a2u 2 and 5 Cpd II ( 5 2) is δx2-y2 1 π*xz 1 π*yz 1 σ*xy 1 a2u 2 . Hydrogen atom abstraction from substrate by Cpd I, typically fills the a2u orbital with a second electron. This usually [48][49][50][51][52] creates a hydrogen atom abstraction structure with a Fe-O-H angle of around 120°, i.e., gives a side-on approach of substrate on oxidant. Indeed, the transition states shown in Figure 4 for the reactivity of Cpd I show this side-on approach in line with the literature. An iron (IV)-hydroxo (porphyrin) complex is then formed in the doublet and quartet spin states with orbital occupation δx2-y2 2 π*xz 1 π*yz 1 a2u 2 ΦSub 1 , whereby ΦSub is the radical on the substrate which is up-spin in the quartet spin state but down-spin in the doublet spin state. Indeed, the group spin densities show that the spin on the porphyrin ring drops from −1/ + 1 in 2 Cpd I/ 4 Cpd I in the reactants to −0.85/0.73 in 2 TSHA,CpdI/ 4 TSHA,CpdI, respectively ( Figure 4). These values drop further to −0.15 and −0.13 for the radical intermediates with concomitant increase of spin density on the substrate moiety.
The reaction of Cpd II with TEMPOH also leads to hydrogen atom abstraction and the donation of another electron into the oxidant set of orbitals. The lowest available virtual orbital for Cpd II is the σ*z2 orbital, which is located along the Fe-O bond. This will lead to an iron (III)hydroxo(porphyrin) with orbital occupation of δx2-y2 2 π*xz 1 π*yz 1 σ*z2 1 a2u 2 ΦSub 1 in the triplet spin state. The substrate radical in the triplet spin state is down-spin, while the metal-based orbitals have an Hydrogen atom abstraction from substrate by Cpd I, typically fills the a 2u orbital with a second electron. This usually [48][49][50][51][52] creates a hydrogen atom abstraction structure with a Fe-O-H angle of around 120 • , i.e., gives a side-on approach of substrate on oxidant. Indeed, the transition states shown in Figure 4 for the reactivity of Cpd I show this side-on approach in line with the literature. An iron (IV)-hydroxo (porphyrin) complex is then formed in the doublet and quartet spin states with orbital occupation δ x2-y2 2 π* xz 1 π* yz 1 a 2u 2 Φ Sub 1 , whereby Φ Sub is the radical on the substrate which is up-spin in the quartet spin state but down-spin in the doublet spin state. Indeed, the group spin densities show that the spin on the porphyrin ring drops from −1/ + 1 in 2 Cpd I/ 4 Cpd I in the reactants to −0.85/0.73 in 2 TS HA,CpdI / 4 TS HA,CpdI , respectively ( Figure 4). These values drop further to −0.15 and −0.13 for the radical intermediates with concomitant increase of spin density on the substrate moiety.
The reaction of Cpd II with TEMPOH also leads to hydrogen atom abstraction and the donation of another electron into the oxidant set of orbitals. The lowest available virtual orbital for Cpd II is the σ* z2 orbital, which is located along the Fe-O bond. This will lead to an iron (III)-hydroxo(porphyrin) with orbital occupation of δ x2-y2 2 π* xz 1 π* yz 1 σ* z2 1 a 2u 2 Φ Sub 1 in the triplet spin state. The substrate radical in the triplet spin state is down-spin, while the metal-based orbitals have an unpaired up-spin electron. The transition state for hydrogen atom abstraction ( Figure 5) indeed gives negative spin density on the substrate group (−0.41) and a large spin density on the FeO group (2.54) and implicates electron transfer of an up-spin electron from the substrate into the σ* z2 orbital along the iron-oxo bond. In the quintet spin state, the σ* xy orbital is also filled with one electron and the hydrogen atom abstraction creates a quintet spin iron (III)-hydroxo (porphyrin) complex with orbital occupation δ x2-y2 1 π* xz 1 π* yz 1 σ* z2 1 σ* xy 1 a 2u 2 Φ Sub 1 . The group spin densities show negative spin density arising on the substrate moiety (−0.28 in 5 TS HA,CpdII ), while the spin on the FeO group increases to 4.22 ( Figure 5). Thus, in both triplet and quintet Cpd II, the electron transfer from substrate to oxidant is into the σ* z2 orbital, which is located along the Fe-O axis. Therefore, the hydrogen atom abstraction is aligned with this axis and we see the substrate approaching from the top. Consequently, the structures of the hydrogen atom transition states depend on the electron transfer processes, as seen previously for iron (IV)-oxo intermediates.

QM/MM Studies
As shown in the previous section, the hydrogen atom abstraction transition states by Cpd I and Cpd II models are very different. While in the former case the substrate approaches from the side, in the latter case a top-approach is found. Obviously in enzymatic systems not always an ideal approach is possible due to the size and shape of the substrate binding pocket and often it is seen that enzymes catalyze stereoselective and regiospecific reaction mechanisms [53][54][55][56]. Furthermore, sometimes small model complexes do not reproduce the correct charge and spin distributions of the real P450 active site due to the absence of, for example, hydrogen bonding interactions to the axial thiolate group [57][58][59][60][61], and consequently full enzymatic approaches are more realistic. Indeed, it was shown previously that hydrogen bonding interactions toward thiolate can affect its charge and pull radical density away from the heme. Therefore, we decided to explore the hydrogen atom abstraction pathway by P450 Cpd I and Cpd II using QM/MM methods. Details of the methodology and set-up have been reviewed and explained in detail elsewhere [62,63]. In general, we take a pdb file from the literature, here 4JWU [64] was used, and modify it to create reactant complexes of either Cpd I + TEMPOH or Cpd II + TEMPOH. From the molecular dynamics simulation several snapshots were tested as starting structures for the QM/MM calculations. The choice of QM region is carefully considered and based on the bonds that are broken, hydrogen bonding interactions and salt bridges. A detailed explanation on how to set up a QM/MM calculation and what to consider as QM region is given in Ref [62]. The QM region used in the QM/MM calculations is displayed in Figure 7 and contained the heme without side chains, methylthiolate for the cysteinate axial ligand and the full TEMPOH substrate. In addition, the aromatic side-chain of Tyr 96 , the side chain of Val 247 and a methanol group representing the Thr 252 side chain were included in the QM region. In all places where the border between the QM and MM regions contains a chemical bond, we inserted link atoms (hydrogens).  [65][66][67]. Furthermore, the experimental crystal structure coordinates and spectroscopic analysis also found Fe-O distances of Cpd I and Cpd II of heme enzymes in close proximity of our computational distances for related heme proteins [68].
A detailed explanation on how to set up a QM/MM calculation and what to consider as QM region is given in Ref [62]. The QM region used in the QM/MM calculations is displayed in Figure 7 and contained the heme without side chains, methylthiolate for the cysteinate axial ligand and the full TEMPOH substrate. In addition, the aromatic side-chain of Tyr96, the side chain of Val247 and a methanol group representing the Thr252 side chain were included in the QM region. In all places where the border between the QM and MM regions contains a chemical bond, we inserted link atoms (hydrogens).   [65][66][67]. Furthermore, the experimental crystal structure coordinates and spectroscopic analysis also found Fe-O distances of Cpd I and Cpd II of heme enzymes in close proximity of our computational distances for related heme proteins [68]. Subsequently, we explored the hydrogen atom abstraction barrier from TEMPOH by Cpd I using QM/MM for snapshots Sn250 and Sn400, see Figure 9. Initially, we performed a geometry scan for the approach of TEMPOH to Cpd I by freezing the FeO-H distance at each step and reoptimizing the rest of the coordinates. The results in Figure 9 show that cleaving the O-H bond of TEMPOH by Cpd I behaves in a different way than that found by DFT on model compounds. In particular, the cleavage has a barrier of about 10 kcal mol −1 that builds up during TEMPOH's approach to Cpd I and which is due to the rearrangement of the protein environment. During the approach, between the FeO·HTEMPOH distance of 2.00 to 4.00 Å, an electron transfer from TEMPOH to the porphyrin group takes place that creates the radical TEMPOH +• cation radical and an [Fe IV = O (Por) (Cys)] − compound, i.e., Cpd II. Thus, a long-range electron transfer by TEMPOH to Cpd I gives Cpd II and oxidized substrate (TEMPOH +• ). The electron transfer, which is expected to have a very low barrier, eases the approach of TEMPOH and reduces a hypothetical higher barrier. At the FeO-HTEMPOH distance of 2.00 Å, the TEMPOH's O-H distance is still 1.03 Å and, from this distance, the scan energy starts decreasing Subsequently, we explored the hydrogen atom abstraction barrier from TEMPOH by Cpd I using QM/MM for snapshots Sn 250 and Sn 400 , see Figure 9. Initially, we performed a geometry scan for the approach of TEMPOH to Cpd I by freezing the FeO-H distance at each step and reoptimizing the rest of the coordinates. The results in Figure 9 show that cleaving the O-H bond of TEMPOH by Cpd I behaves in a different way than that found by DFT on model compounds. In particular, the cleavage has a barrier of about 10 kcal mol −1 that builds up during TEMPOH's approach to Cpd I and which is due to the rearrangement of the protein environment. During the approach, between the FeO·H TEMPOH distance of 2.00 to 4.00 Å, an electron transfer from TEMPOH to the porphyrin group takes place that creates the radical TEMPOH +• cation radical and an [Fe IV = O (Por) (Cys)] − compound, i.e., Cpd II. Thus, a long-range electron transfer by TEMPOH to Cpd I gives Cpd II and oxidized substrate (TEMPOH +• ). The electron transfer, which is expected to have a very low barrier, eases the approach of TEMPOH and reduces a hypothetical higher barrier. At the FeO-H TEMPOH distance of 2.00 Å, the TEMPOH's O-H distance is still 1.03 Å and, from this distance, the scan energy starts decreasing until products are reached. Thus, TEMPOH's O-H bond cleavage takes place in a second barrierless stage where a proton is abstracted from TEMPOH +• , creating the TEMPO • radical. Consequently, the global process that occurs in the protein is better described as a proton-coupled-electron transfer (PCET), unlike the hydrogen atom transfer (HAT) process that takes place in model complexes. . QM/MM calculated geometry scans for substrate approach on the iron (IV)-oxo species of 2 Cpd I as calculated for snapshots Sn250 (salmon) and Sn400 (blue). Each point represents a full geometry optimization with fixed FeO···HO distance between oxo and TEMPOH proton. In the rounded boxes, Mulliken spin densities (in atomic units) of the TEMPOH substrate (green) and the porphyrin ring (plum) are represented, as well as the Ooxo-HTEMPOH distance (Å) for selected structures (black).
From the maximum points in the geometry scans of Figure 9, we attempted to locate transition states. However, the energetic fluctuation during the scan does not represent a barrier of a chemical reaction through bond breaking and/or bond forming, but is related to the structural reorganization in the substrate binding-pocket. In particular, during approach of the substrate on the oxidant, hydrogen bonding interactions are weakened with one group and strengthened with other groups. In addition, water molecules migrate in the substrate binding-pocket and also undergo hydrogen bonding interactions. As such caution should be taken when interpreting these geometry scans. Specifically, the maximum energy points are at a relatively long O-H distance of around 2 Å for Sn400 and 2.8 Å for Sn250. The O-H distances in the scan maxima structures are well longer than normally found in hydrogen atom abstraction transition states that are in the order of 1.2-1.4 Å [42,[47][48][49][69][70][71][72][73]. Therefore, no proper hydrogen atom abstraction transition states could be characterized. Since the scan energy profile at short distances gives little evidence of a proper hydrogen atom abstraction barrier, it implies that also in the protein the hydrogen atom abstraction barrier from TEMPOH by Cpd I is negligible. Consequently, the small model complex captures the potential energy profile of hydrogen atom abstraction quite well and predicts negligible hydrogen atom abstraction barriers probably because the substrate fits into the substrate binding pocket well and can approach the oxidant under the ideal angle.
In summary, the hydrogen atom abstraction by Cpd I of TEMPOH encounters a barrier due to electron transfer and protein movement and, of course, varies with the snapshot chosen. During this process the O-H distance in TEMPOH remains short (see black data in the rounded boxes in Figure 9). The actual hydrogen atom abstraction to oxidant, however, is virtually barrierless and no transition From the maximum points in the geometry scans of Figure 9, we attempted to locate transition states. However, the energetic fluctuation during the scan does not represent a barrier of a chemical reaction through bond breaking and/or bond forming, but is related to the structural reorganization in the substrate binding-pocket. In particular, during approach of the substrate on the oxidant, hydrogen bonding interactions are weakened with one group and strengthened with other groups. In addition, water molecules migrate in the substrate binding-pocket and also undergo hydrogen bonding interactions. As such caution should be taken when interpreting these geometry scans. Specifically, the maximum energy points are at a relatively long O-H distance of around 2 Å for Sn 400 and 2.8 Å for Sn 250 . The O-H distances in the scan maxima structures are well longer than normally found in hydrogen atom abstraction transition states that are in the order of 1.2-1.4 Å [42,[47][48][49][69][70][71][72][73]. Therefore, no proper hydrogen atom abstraction transition states could be characterized. Since the scan energy profile at short distances gives little evidence of a proper hydrogen atom abstraction barrier, it implies that also in the protein the hydrogen atom abstraction barrier from TEMPOH by Cpd I is negligible. Consequently, the small model complex captures the potential energy profile of hydrogen atom abstraction quite well and predicts negligible hydrogen atom abstraction barriers probably because the substrate fits into the substrate binding pocket well and can approach the oxidant under the ideal angle.
In summary, the hydrogen atom abstraction by Cpd I of TEMPOH encounters a barrier due to electron transfer and protein movement and, of course, varies with the snapshot chosen. During this process the O-H distance in TEMPOH remains short (see black data in the rounded boxes in Figure 9). The actual hydrogen atom abstraction to oxidant, however, is virtually barrierless and no transition states could be characterized. This process is the same regardless of the snapshot we investigated.
Next, we investigated the hydrogen atom abstraction from TEMPOH by 3 Cpd II and the results are given in Figure 10. For Cpd II a proper hydrogen atom abstraction transition state could be located with an energy barrier of 22.3 kcal mol −1 with respect to a reactant complex. Again, the QM/MM barrier is higher in energy than the gas-phase barrier due to the protein environment that prevents an ideal approach of substrate to oxidant and thereby raises the barriers with respect to the gas-phase model complexes. Consequently, a tight binding pocket will have a barrier-raising effect but on the other hand may trigger a regioselective or stereospecific reaction path. The enzyme, therefore, balances a regioselective reaction mechanism at a cost of a slower reaction rate. For Cpd II the approach of substrate on the oxidant also gives a geometry scan that fluctuates due to the forming and breaking of hydrogen bonding interactions. However, upon close distance a high barrier appears for the bond breaking and forming processes during the reaction. other hand may trigger a regioselective or stereospecific reaction path. The enzyme, therefore, balances a regioselective reaction mechanism at a cost of a slower reaction rate. For Cpd II the approach of substrate on the oxidant also gives a geometry scan that fluctuates due to the forming and breaking of hydrogen bonding interactions. However, upon close distance a high barrier appears for the bond breaking and forming processes during the reaction. The QM/MM optimized geometry for hydrogen atom abstraction from TEMPOH by Cpd II is shown in Figure 10. The structure is reactant-like with short TEMPO-H distance of 1.14 Å, whereas in the DFT model a value of 1.33 Å was found. Similarly, the FeO-H distance is relatively long in QM/MM (1.23 Å) and much shorter in the gas-phase model (1.10 Å). These differences in optimized geometry are due to several hydrogen bonding interactions surrounding the iron(IV)-oxo and substrate groups and, in particular, three crystal water molecules are highlighted in Figure 10 that are involved in hydrogen bonding interactions to the oxo as well as amide groups and thereby affect the hydrogen transfer mechanism. Obviously, hydrogen bonding interactions from crystal water molecules in the substrate binding pocket have strong effects on the position of the substrate with respect to the oxidant as well as a charge-donation that influences transition state structures and ultimately barrier heights. In previous studies the effect of hydrogen bonding interactions to the oxo group were shown to raise hydrogen atom abstraction barriers significantly [74] in line with the difference between DFT and QM/MM as seen here.

Discussion
In this work a series of DFT and QM/MM results are presented on the reactivity of cytochrome P450 Cpd I and Cpd II intermediates with a model substrate, namely TEMPOH. It is seen that the lowest reaction barriers are obtained with Cpd I, although those found for Cpd II are also accessible The QM/MM optimized geometry for hydrogen atom abstraction from TEMPOH by Cpd II is shown in Figure 10. The structure is reactant-like with short TEMPO-H distance of 1.14 Å, whereas in the DFT model a value of 1.33 Å was found. Similarly, the FeO-H distance is relatively long in QM/MM (1.23 Å) and much shorter in the gas-phase model (1.10 Å). These differences in optimized geometry are due to several hydrogen bonding interactions surrounding the iron(IV)-oxo and substrate groups and, in particular, three crystal water molecules are highlighted in Figure 10 that are involved in hydrogen bonding interactions to the oxo as well as amide groups and thereby affect the hydrogen transfer mechanism. Obviously, hydrogen bonding interactions from crystal water molecules in the substrate binding pocket have strong effects on the position of the substrate with respect to the oxidant as well as a charge-donation that influences transition state structures and ultimately barrier heights. In previous studies the effect of hydrogen bonding interactions to the oxo group were shown to raise hydrogen atom abstraction barriers significantly [74] in line with the difference between DFT and QM/MM as seen here.

Discussion
In this work a series of DFT and QM/MM results are presented on the reactivity of cytochrome P450 Cpd I and Cpd II intermediates with a model substrate, namely TEMPOH. It is seen that the lowest reaction barriers are obtained with Cpd I, although those found for Cpd II are also accessible at room temperature depending on the method and model. As such both Cpd I and Cpd II should be able to activate substrates with weak or medium strength C-H bonds efficiently. In the following, we will try to rationalize the reactivity differences and particularly focus on the thermodynamics of the reactions.
In principle, the hydrogen atom abstraction reaction from a substrate (TEMPOH) by Cpd I or [Fe IV (O)(heme +• )Cys] gives an iron(IV)-hydroxo species and a radical rest group (TEMPO • ) as described with Equation (1). Thermodynamically, the energy for Equation (1) (∆H Eq1 ) is the difference in energy of the bond dissociation energy of the substrate O-H bond that is broken (BDE OH, TE ), Equation (2), and the bond dissociation energy of the O-H bond of the iron (IV)-hydroxo complex that is formed (BDE OH,CpdI ), Equation (3) [42,48,49,[75][76][77][78][79][80][81][82]. As such the reaction energy for hydrogen atom abstraction from substrate by Cpd I (Equation (4) To this end we calculated the BDE OH values of TEMPOH, the iron (IV)-hydroxo (2H + ) and iron (III)-hydroxo complexes, i.e., BDE OH,TE , BDE OH,CpdI and BDE OH,CpdII , respectively. TEMPOH has a very weak O-H bond strength of 56.4 kcal mol −1 . Typical values for aliphatic C-H bond strengths calculated at the same level of theory are 81.9 kcal mol −1 for abstraction of a hydrogen atom from the α-position of ethylbenzene and 93.3 kcal mol −1 for cyclohexane [36]. Clearly, the O-H bond of TEMPOH is weak and it should not cost oxidants much energy to abstract its hydrogen atom. Indeed, a negligible barrier is found for Cpd I and a small barrier for Cpd II is found in agreement with the strength of the O-H bond of the substrate. Based on the BDE OH values [36], we predict hydrogen atom abstraction driving forces for Cpd I of −35.3 kcal mol −1 and for Cpd II of −31.7 kcal mol −1 . These values are in good quantitative agreement with the energy differences between reactants and intermediates seen in Figures 4 and 5 above. The small energetic differences result from dispersion and intermolecular interaction energies. Figure 11 summarizes the thermochemistry of possible hydrogen atom abstraction and electron transfer processes for both Cpd I and Cpd II. Pathways from left to right represent hydrogen atom abstraction from TEMPOH that give the iron (IV/III)-hydroxo complexes. Also given on the diagonal axis in Figure 11 are the electron affinities (EA) of Cpd I and protonated Cpd II (2H + ). Note that the electron affinity of Cpd II is 50.1 kcal mol −1 using the same computational methods.
Thus, the BDE OH of Cpd I is larger than that of Cpd II and therefore, Cpd I will react with lower hydrogen atom abstraction barriers with substrates. However, the energy difference is not dramatic, i.e., just a few kcal mol −1 in energy, so that both oxidants should be able to activate substrates as easily. Consequently, if excess reduction appears in P450 isozymes and Cpd II is formed, the enzyme will still be active in substrate activation processes. interaction energies. Figure 11 summarizes the thermochemistry of possible hydrogen atom abstraction and electron transfer processes for both Cpd I and Cpd II. Pathways from left to right represent hydrogen atom abstraction from TEMPOH that give the iron (IV/III)-hydroxo complexes. Also given on the diagonal axis in Figure 11 are the electron affinities (EA) of Cpd I and protonated Cpd II (2H + ). Note that the electron affinity of Cpd II is 50.1 kcal mol −1 using the same computational methods. Thus, the BDEOH of Cpd I is larger than that of Cpd II and therefore, Cpd I will react with lower hydrogen atom abstraction barriers with substrates. However, the energy difference is not dramatic, i.e., just a few kcal mol −1 in energy, so that both oxidants should be able to activate substrates as easily.

DFT Model Complexes
DFT model complexes were calculated in the Gaussian-09 program package [104] and investigated with density functional theory methods. Geometry optimizations, frequency calculations, geometry scans and intrinsic reaction coordinate profiles were performed at the unrestricted B3LYP level of theory [105,106] and utilized an LACVP basis set with core potential [107] on iron and 6-31G on the rest of the atoms; basis set BS1. Energies were improved at the single point level of theory with an LACV3P+ basis set on iron and 6-311+G* on the rest of the atoms; basis set BS2. Previously, we showed that UB3LYP/BS2 optimized geometries and potential energy landscapes were within a few tenths of a kcal mol −1 from those obtained at UB3LYP/BS2//UB3LYP/BS1 [42,78,100] Moreover, spin state orderings were reproduced well with this method [108]. In addition, the calculations included the polarized continuum model with a dielectric constant of ε = 5.697, which is a typical value for the active site within proteins. Note that earlier test calculations with a range of dielectric constants only gave minor changes to spin state orderings and relative energies and the major difference found was upon changing from the gas-phase to a solvent model [109].
The chemical model was based on a typical P450 active site structure and included an iron embedded in protoporphyrin IX (Por) without substituents on the periphery. The axial cysteinate ligand of iron was abbreviated to thiolate as that was shown to be a better mimic than methylthiolate [110,111].

QM/MM Set Up
The QM/MM model was built from the 4JWU protein databank (pdb) file [64], which is a bacterial (Pseudomonas putida) P450 isozyme crystal structure that contains the homodimer of the Cytochrome P450 domain with its redox partner putidaredoxin [112]. We removed the redox partner and selected a complete single strand (both are identical). We inserted TEMPOH substrate manually into the enzyme structure on the distal site of the heme, and, in addition, replaced the iron(III)-water group by an iron(IV)-oxo species with an Fe-O starting bond distance of 1.65 Å as is typically found in heme Cpd I structures in calculations [113,114]. Subsequently, hydrogen atoms were added to the structure using the online tool pdb2pqr at pH = 7 [115]. All acidic and basic protein groups were visually inspected for having the correct protonation state and in our QM/MM models all carboxylate groups (Asp/Glu) were found to be deprotonated and all Arg and Lys amino acid side chains were protonated. Note that this particular P450 isozyme contains no histidine residues in the protein. In the final stage of the structure preparation, we performed an iterative solvation procedure (left graphic in Figure 12) and added TIP3P water molecules to obtain a final structure containing 32,447 atoms including 8657 water molecules. The structure was neutralized with ten Mg 2+ and four Cl − ions on the surface of the protein. In the next stages of the QM/MM set up we applied a stepwise heating protocol to 300 K using the Charmm force field [116] followed by full equilibration whereby the protein backbone was kept fixed. Thereafter, a full molecular dynamics (MD) simulation was run for 500 ps, see Figure 12 (right), from which we selected three low-energy snapshots after 250, 400 and 500 ps for the actual QM/MM calculations: Sn 250 , Sn 400 and Sn 500 . (right), from which we selected three low-energy snapshots after 250, 400 and 500 ps for the actual QM/MM calculations: Sn250, Sn400 and Sn500. The model was split into a QM and MM region, whereby the QM part was described by density functional theory, whereas the MM part was calculated with the Charmm forcefield [116]. The QM region was calculated at the UB3LYP/BS3 level of theory [105,106] with BS3 representing a def-SVP basis set on all atoms [117]. The QM and MM regions are interfaced by ChemShell [118] through DL-Poly [119], whereby the charges of the MM region were incorporated into the QM Hamiltonian through electrostatic embedding [120]. The boundary of the QM and MM region was described through the Link-atom approach [121]. Geometries were optimized in QM/MM without constraints, while the MM region had a flexible component of 5 Å around the heme and substrate moieties, while all other atoms in the MM region were fixed in the original snapshot orientation. Cpd I and Cpd II were calculated in the doublet and triplet spin states, respectively.
The geometry scan for the hydrogen atom transfer from TEMPOH to the oxidant was performed through stepwise geometry optimizations at QM/MM level of theory with fixed O-H bond distances, while all other degrees of freedom were minimized. The maxima of these scans were subjected to a full transition state search.

Conclusions
In this work a computational study is presented in the reaction of P450 Cpd I and Cpd II with TEMPOH using a combined density functional theory and QM/MM approach. These studies The model was split into a QM and MM region, whereby the QM part was described by density functional theory, whereas the MM part was calculated with the Charmm forcefield [116]. The QM region was calculated at the UB3LYP/BS3 level of theory [105,106] with BS3 representing a def-SVP basis set on all atoms [117]. The QM and MM regions are interfaced by ChemShell [118] through DL-Poly [119], whereby the charges of the MM region were incorporated into the QM Hamiltonian through electrostatic embedding [120]. The boundary of the QM and MM region was described through the Link-atom approach [121]. Geometries were optimized in QM/MM without constraints, while the MM region had a flexible component of 5 Å around the heme and substrate moieties, while all other atoms in the MM region were fixed in the original snapshot orientation. Cpd I and Cpd II were calculated in the doublet and triplet spin states, respectively.
The geometry scan for the hydrogen atom transfer from TEMPOH to the oxidant was performed through stepwise geometry optimizations at QM/MM level of theory with fixed O-H bond distances, while all other degrees of freedom were minimized. The maxima of these scans were subjected to a full transition state search.

Conclusions
In this work a computational study is presented in the reaction of P450 Cpd I and Cpd II with TEMPOH using a combined density functional theory and QM/MM approach. These studies investigate O-H hydrogen atom abstractions of the two potential oxidants and discuss their feasibility in an enzymatic setting. The results show that both Cpd I and Cpd II can react with substrates with weak O-H bonds. Interestingly, Cpd I in the protein reacts via a proton-coupled electron transfer process rather than hydrogen atom abstraction. Thus, the protein enables a long-range electron transfer from substrate to Cpd I followed by a proton transfer. These studies show distinct differences between gas-phase DFT models and protein structures. Overall, both Cpd I and Cpd II are found to be suitable oxidants for TEMPOH activation.