Theoretical Study of Hydroxylation of α- and β-Pinene by a Cytochrome P450 Monooxygenase Model

Previous studies on biocatalytic transformations of pinenes by cytochrome P450 (CYP) enzymes reveal the formation of different oxygenated products from a single substrate due to the multistate reactivity of CYP and the many reactive sites in the pinene scaffold. Up until now, the detailed mechanism of these biocatalytic transformations of pinenes have not been reported. Hereby, we report a systematic theoretical study of the plausible hydrogen abstraction and hydroxylation reactions of α- and β-pinenes by CYP using the density functional theory (DFT) method. All DFT calculations in this study were based on B3LYP/LAN computational methodology using the Gaussian09 software. We used the B3LYP functional with corrections for dispersive forces, BSSE, and anharmonicity to study the mechanism and thermodynamic properties of these reactions using a bare model (without CYP) and a pinene-CYP model. According to the potential energy surface and Boltzmann distribution for radical conformers, the major reaction products of CYP-catalyzed hydrogen abstraction from β-pinene are the doublet trans (53.4%) and doublet cis (46.1%) radical conformer at delta site. The formation of doublet cis/trans hydroxylated products released a total Gibbs free energy of about 48 kcal/mol. As for alpha pinene, the most stable radicals were trans-doublet (86.4%) and cis-doublet (13.6%) at epsilon sites, and their hydroxylation products released a total of ~50 kcal/mol Gibbs free energy. Our results highlight the likely C-H abstraction and oxygen rebounding sites accounting for the multi-state of CYP (doublet, quartet, and sextet spin states) and the formation of different conformers due to the presence of cis/trans allylic hydrogen in α-pinene and β-pinene molecules.


Introduction
Cytochrome P450 (CYP) enzymes represent a notable family of biocatalysts involved in the metabolism of both exogenous and endogenous compounds using thiolate-ligated heme scaffolds ( Figure 1) [1,2]. For instance, cytochrome P450 2E1 (CYP2E1) is one of the major human hepatic CYP enzymes with direct relevance in alcoholism [3]. CYPs are responsible for the metabolism of about 90% of administered drugs in men [4]. CYP Biocatalytic transformation of monoterpenes, such as those catalysed by CYPs, is an active area of research. Monoterpenes and their oxygenated products are important raw materials for a variety of value-added products such as pharmaceuticals, flavors, and fragrances [7,15,16]. Bicyclic monoterpenes such as α-pinene and β-pinene are major components of turpentine oil, the distillate of certain pine trees [17]. Important pharmacological effects [18] of α-and β-pinenes were demonstrated in several studies such as their antitumor [17,19], antimicrobial [20], and anticonvulsant [21] activities. In addition, monoterpenes represent around 11% of the volatile organic compounds emitted into the atmosphere, with α-and β-pinenes being among the most abundant terpenes in the troposphere. Monoterpenes thus play considerable roles in atmospheric chemistry where oxidation of Biocatalytic transformation of monoterpenes, such as those catalysed by CYPs, is an active area of research. Monoterpenes and their oxygenated products are important raw materials for a variety of value-added products such as pharmaceuticals, flavors, and fragrances [7,15,16]. Bicyclic monoterpenes such as α-pinene and β-pinene are major components of turpentine oil, the distillate of certain pine trees [17]. Important pharmacological effects [18] of αand β-pinenes were demonstrated in several studies such as their antitumor [17,19], antimicrobial [20], and anticonvulsant [21] activities. In addition, monoterpenes represent around 11% of the volatile organic compounds emitted into the atmosphere, with αand β-pinenes being among the most abundant terpenes in the troposphere. Monoterpenes thus play considerable roles in atmospheric chemistry where oxidation of these products can take place by atmospheric oxidants such as nitrate radicals and ozone molecules [22][23][24].
Several mechanistic studies on biocatalytic transformations by CYP enzymes revealed multi-state reactivity in many cases with different spin states and reaction routes that are close in energy. These results explain the generation of various products from the same initial substrate [2,37,38]. A study on CYP125 enzyme presents an interesting example of these multi-transformations caused by the competition among various possible catalytic routes. CYP125 sequentially oxidizes the terminal methyl of cholest-4-en-3-one to its alcohol, aldehyde, and lastly cholest-4en-3-one-26-acid, following the consensus catalytic cycle of most oxidizing CYPs [39]. Five additional metabolites were also identified in this reaction produced by the deformylation of the aldehyde via peroxohemiacetal formation [39]. Other studies showed that CYP-catalyzed epoxidation reaction generates one product in the doublet state and a mixture of products in the quartet state [4,40].
Thus, regarding the reactivity of CYP, both theoretical and experimental studies show no clear-cut mechanism. In this study, we focus on the hydrogen-abstraction and hydroxylation steps of pinenes catalyzed by CYP. Understanding H-abstraction reactions has a key role in deciphering complex physical organic concepts and reaction systems such as oxidation of biological systems, autoignition, and combustion of hydrocarbons, as well as oxidations in atmospheric chemistry. On the other hand, hydroxylation reactions and rebound mechanisms by CYP are vital steps in drug metabolism [4,[41][42][43][44]. A detailed investigation of the possible mechanisms and interactions of the CYP/pinene system is important for predicting the various possible products observed experimentally. Up to the present, such mechanistic studies have not been reported. Herein, we present a systematic mechanistic study on the different possible hydrogen-abstraction and hydroxylation reactions of CYP on αand β-pinenes in multi-state reactivity (different spin states) on the basis of density functional theory (DFT).

Alpha-and Beta-Pinene Bare Systems
We first started by computing the relative energies of the radicals and the hydroxylated products of the α-pinene and β-pinene without a CYP450 model using the B3LYP/LAN methodology. We arbitrarily named the carbon atoms in these compounds with Greek letters to identify the sites from which hydrogen atoms can be abstracted. The structures in Tables 1 and 2 illustrate the employed nomenclature. The nu (ï) site was employed to differentiate the reaction pathways of the cis hydrogen (same side) and trans hydrogen (opposite site) on the epsilon C. Starting with α-pinene (Table 1), our calculations showed that the most stable radical was formed with the removal of the hydrogen atom from the epsilon site. All other radicals, with the exception of the one formed at the gamma site (~2 kcal/mol), had considerably higher energies (17-29 kcal/mol). The steric strain of the bridgehead alpha radical might have contributed to making the alpha radical considerably less stable (21.7 kcal/mol) with lack of resonance [41,42]. The order of stability of the allyl Table 1. Computed relative energies of the radicals and hydroxylated products of α-pinene. The term "cor" refers to the dispersion, BSSE, and Grimme quasi-harmonic corrections to the calculated absolute Gibbs free energies (in Hartree). Energies in kcal/mol were calculated with reference to the most stable conformation for the hydroxylated radical and product. employed to differentiate the reaction pathways of the cis hydrogen (same side) and trans hydrogen (opposite site) on the epsilon C. Starting with α-pinene (Table 1), our calculations showed that the most stable radical was formed with the removal of the hydrogen atom from the epsilon site. All other radicals, with the exception of the one formed at the gamma site (~2 kcal/mol), had considerably higher energies (17-29 kcal/mol). The steric strain of the bridgehead alpha radical might have contributed to making the alpha radical considerably less stable (21.7 kcal/mol) with lack of resonance [41,42]. The order of stability of the allyl radicals was as follows: epsilon > gamma > alpha. Considering the computed energies of the hydroxylated structures, the products formed at the zeta, delta, and alpha sites were found to be more stable (0-2 kcal/mol) than the other hydroxylated ones. Table 1. Computed relative energies of the radicals and hydroxylated products of α-pinene. The term "cor" refers to the dispersion, BSSE, and Grimme quasi-harmonic corrections to the calculated absolute Gibbs free energies (in Hartree). Energies in kcal/mol were calculated with reference to the most stable conformation for the hydroxylated radical and product.

Radical
Gcor (Eh) With respect to β-pinene (Table 2), our geometry optimizations showed that the most stable radical was formed at the delta allylic site. Similar to the α-pinene, all the other radicals were far above this structure in the potential energy surface, with the gamma radical being the least stable (~31 kcal/mol). Also in line with the α-pinene, the tertiary bridgehead alpha radical was considerably less stable (21.9 kcal/mol) than the other allylic radical at delta carbon. The most stable hydroxylated products were also found to be the tertiary hydroxylated carbons at zeta and alpha sites.

CYP Monooxygenase Model
Modelling biochemical reactions is demanding, particularly with diverse enzyme families such as CYP. Several model systems have been used for CYP monooxygenases. In this study, we employed one model ( Figure 1) of CYP that showed good efficiency in previous DFT calculation studies. It is important to mention that this model is representative of the CYP family but does not account for all CYP classes. The chosen Fe 4+ O 2− (C20N4H12) − (SH) − model [3,4] is composed of an unsubstituted porphyrin ring instead of a porphyrin IX ring and an SH-group instead of SMe or protein side chains as axial Cys437 [45].

R and S Enantiomers Cis/Trans Reaction
We first computed the cis/trans hydrogen transfer to the CYP reaction using both R and S enantiomers of α-pinene ( Figure 2 and Table 3) and β-pinene ( Figure S1 and Table  S1, Supplementary Materials). The results show that there are no significant structural differences and consequently no significant energy differences between the two enantiomers. This result can be explained by the absence of any electronic or steric clashes differences introduced in the hydrogen transfer reaction mechanism by the R or S pinene isomers, as expected.  With respect to β-pinene (Table 2), our geometry optimizations showed that the most stable radical was formed at the delta allylic site. Similar to the α-pinene, all the other radicals were far above this structure in the potential energy surface, with the gamma radical being the least stable (~31 kcal/mol). Also in line with the α-pinene, the tertiary bridgehead alpha radical was considerably less stable (21.9 kcal/mol) than the other allylic radical at delta carbon. The most stable hydroxylated products were also found to be the tertiary hydroxylated carbons at zeta and alpha sites.

CYP Monooxygenase Model
Modelling biochemical reactions is demanding, particularly with diverse enzyme families such as CYP. Several model systems have been used for CYP monooxygenases. In this study, we employed one model ( Figure 1) of CYP that showed good efficiency in previous DFT calculation studies. It is important to mention that this model is representative of the CYP family but does not account for all CYP classes. The chosen Fe 4+ O 2− (C 20 N 4 H 12 ) − (SH) − model [3,4] is composed of an unsubstituted porphyrin ring instead of a porphyrin IX ring and an SH-group instead of SMe or protein side chains as axial Cys 437 [45].

R and S Enantiomers Cis/Trans Reaction
We first computed the cis/trans hydrogen transfer to the CYP reaction using both R and S enantiomers of α-pinene ( Figure 2 and Table 3) and β-pinene ( Figure S1 and Table S1, Supplementary Materials). The results show that there are no significant structural differences and consequently no significant energy differences between the two enantiomers. This result can be explained by the absence of any electronic or steric clashes differences introduced in the hydrogen transfer reaction mechanism by the R or S pinene isomers, as expected. [3,4] is composed of an unsubstituted porphyrin ring instead of a porphyrin IX ring and an SH-group instead of SMe or protein side chains as axial Cys437 [45].

R and S Enantiomers Cis/Trans Reaction
We first computed the cis/trans hydrogen transfer to the CYP reaction using both R and S enantiomers of α-pinene ( Figure 2 and Table 3) and β-pinene ( Figure S1 and Table  S1, Supplementary Materials). The results show that there are no significant structural differences and consequently no significant energy differences between the two enantiomers. This result can be explained by the absence of any electronic or steric clashes differences introduced in the hydrogen transfer reaction mechanism by the R or S pinene isomers, as expected.   Table 3. Gibbs free energy coordinates of the reaction R/S and cis/trans paths for the α-pinene. Data are provided in kcal/mol and relative to the most stable structure formed in reaction coordinate 1, the gamma conformer. RC is reaction coordinate.

CYP-Catalyzed Hydroxylation of β-Pinene
Based on the preliminary results confirming the stability of the allyl radical from the delta site of the β-pinene and the interesting high energy of the allyl radical generated at the alpha site, we focused on studying the catalytic interactions of the CYP model with allyl-type hydrogens only. The higher stability of allyl radicals versus other non-delocalized radicals can be rationalized in terms of the delocalization of the paired electrons of the olefin bond and the non-paired electron [41,42,46].
First, we studied the catalytic C-H abstraction followed by hydroxylation of the delta site of β-pinene by CYP in the doublet, quartet, and sextet spin states since transition metal complexes (and CYP) can react through multistate reactivity patterns with several close-lying spin states [2,37,40,47,48]. Figures 3 and 4 illustrate the computed mechanisms using the B3LYP/LAN methodology. Table 4 presents absolute energies, kinetic, and thermodynamic calculated parameters (see Tables S2-S10, S13 and Figure S2 in the Supplementary Materials for the coordinates of the energy profile and the calculations in detail). Since there are two hydrogen atoms in the delta carbon, the hydrogen abstraction can occur in a cis/trans fashion with respect to the η site of β-pinene. We set up to use one configuration for cis (R-cis) and the other for trans (S-trans) since our preliminary calculations of the R and S structures of the pinene system showed similar electronic energies (vide supra). The sextet spin state was found to be considerably high in energy ( Figure 4B), in line with previous studies in the literature [49], so it was excluded in further calculations. delta site of the β-pinene and the interesting high energy of the allyl radical generated at the alpha site, we focused on studying the catalytic interactions of the CYP model with allyl-type hydrogens only. The higher stability of allyl radicals versus other non-delocalized radicals can be rationalized in terms of the delocalization of the paired electrons of the olefin bond and the non-paired electron [41,42,46].
First, we studied the catalytic C-H abstraction followed by hydroxylation of the delta site of β-pinene by CYP in the doublet, quartet, and sextet spin states since transition metal complexes (and CYP) can react through multistate reactivity patterns with several closelying spin states [2,37,40,47,48]. Figures 3 and 4 illustrate the computed mechanisms using the B3LYP/LAN methodology. Table 4 presents absolute energies, kinetic, and thermodynamic calculated parameters (see Tables S2-S10, S13 and Figure S2 in the Supplementary Materials for the coordinates of the energy profile and the calculations in detail). Since there are two hydrogen atoms in the delta carbon, the hydrogen abstraction can occur in a cis/trans fashion with respect to the η site of β-pinene. We set up to use one configuration for cis (R-cis) and the other for trans (S-trans) since our preliminary calculations of the R and S structures of the pinene system showed similar electronic energies (vide supra). The sextet spin state was found to be considerably high in energy ( Figure 4B), in line with previous studies in the literature [49], so it was excluded in further calculations.    We observed that all cis transition state structures have slightly more elevated ener- We observed that all cis transition state structures have slightly more elevated energies than the trans structures. Analysis of the structures suggests that this could be due to the steric clashes between the oxygen atom in the CYP and the β-pinene scaffold. According to our calculations, the doublet and quartet reaction paths (coordinate I in Figure 4B) have similar energies, while the paths of the sextet spin state were found to be high in energy, as mentioned. For the doublet and quartet reaction paths, the hydrogen transfer reactions are characterized by exothermicity and relatively low electronic barriers. As can be observed in Table 4, the general trend in the application of vibrational corrections to the electronic energies led to a decrease in the enthalpy of activation and in the enthalpy of reaction (except to sextet delta trans and alpha quartet structures, in which the effect is nearly null). For the entropy corrections, a less clear trend was observed. Except for the sextet high-spin states and quartet delta trans conformer, the entropic factor led to increasing the activation energy and reaction Gibbs free energy. It is largely known that the two main DFT sources of error are the missing dispersion forces and the basis set superposition error [50]. Thus, as expected, the effects of corrections added to our DFT/LAN method led to lower barriers and more favorable thermodynamics (again, except for the high-spin thermodynamic parameters) for the hydrogen transfer reaction. We observed that the addition of empirical corrections for the dispersive forces lead to a great decrease in the CYP-catalyzed reaction barriers ( Figure S3). The same trend was observed by Mulholland and colleagues [51]. These lower reaction barriers are reasonable in a context of catalyzed biochemical reactions where barriers are expected to be no greater than 18 kcal/mol [52]. Moreover, we observed that the stability order of the radical products is marginally affected by these corrections, in contrast to the stability of the hydroxylated products. Therefore, the applications of these corrections are essential for more accurate predictions of the selectivity and reaction barriers. Table 4. Thermodynamic and kinetic data computed for the hydrogen capture reaction by CYP in β-pinene. The terms ∆ r E, ∆ r H, and ∆ r G refer to electronic, enthalpic, and Gibbs free energy variation in the reaction, respectively. The terms ∆E ‡ , ∆H ‡ , and ∆G ‡ refer to the electronic, enthalpic, and Gibbs free energy reaction barriers, respectively, and the term "cor" refers to the dispersion, BSSE, and Grimme quasi-harmonic corrections to the calculated absolute Gibbs free energies. Hyd ∆ r G cor refers to the hydroxylation free energies computed as the difference between the hydroxylated product and the reactant in coordinate 1. Data are displayed in kcal/mol.

Site
Spin The calculation of the allylic alpha hydrogen abstraction by CYP showed the unexpected instability of the alpha radical that was observed in the calculations using the bare β-pinene system ( Figure 4A). It appears that alpha-H abstraction is destabilizing the bicyclic ring. At reaction coordinate I, the alpha site CYP + β-pinene complex in the doublet state was found to be the most stable structure (0 kcal/mol) over the potential energy surface. Nevertheless, the energy differences among other optimized structures in coordinate I were very small, excluding the sextet state. We are able to identify several other conformers in the doublet and quartet with slightly more elevated energies, as the β-pinene can be freely rotated/translated around the CYP due to its weak interaction with the enzyme. To the alpha site, our calculations showed that the activation barriers are more elevated at the doublet state than the quartet state. However, the hydrogen abstraction from this site is unlikely due to the unfavorable thermodynamics for the alpha radical formation.
We then investigated the rebound mechanism from the β-pinene radical after the hydrogen abstraction by CYP. The transition states connecting the radical and the hydroxylated product for the low and high spin states could not be found in an exhaustive potential energy surfaces search. In the quartet state, we were successful in locating these transition states structures only at a lower theoretical level (B3LYP/6-31G). Searches in the potential energy surface through scan calculations suggested a small electronics barrier of~3.0 kcal/mol. However, at a higher theoretical level, the potential energy surface turns flat and all attempts to optimize these transition states at B3LYP/LAN level were not fruitful. Therefore, in line with previous studies, we consider that the hydroxylation rebound mechanism is barrierless at our theoretical level, with the hydrogen-abstraction by the CYP enzyme as the rate-limiting step [12].
The bare system alpha-delta trend in the relative stability of hydroxylated product repeats in the CYP + β-pinene system. The doublet alpha hydroxylated product is about 3 kcal/mol more stable than the delta hydroxylated product. We observed that highspin sextet structures ( Figure 5) do not give classical rebounded hydroxylated products. The increase of electronic repulsion around the iron atom in higher spin states repels the ligands above and below the plane occupied by the porphyrin system, evidenced by the increased interaction distance between Fe and O and the slight distortion of the Fe coordination sphere.
flat and all attempts to optimize these transition states at B3LYP/LAN level were not fruitful. Therefore, in line with previous studies, we consider that the hydroxylation rebound mechanism is barrierless at our theoretical level, with the hydrogen-abstraction by the CYP enzyme as the rate-limiting step [12].
The bare system alpha-delta trend in the relative stability of hydroxylated product repeats in the CYP + β-pinene system. The doublet alpha hydroxylated product is about ~3 kcal/mol more stable than the delta hydroxylated product. We observed that high-spin sextet structures ( Figure 5) do not give classical rebounded hydroxylated products. The increase of electronic repulsion around the iron atom in higher spin states repels the ligands above and below the plane occupied by the porphyrin system, evidenced by the increased interaction distance between Fe and O and the slight distortion of the Fe coordination sphere.  With these results, we can reason a dynamic equilibrium scenario for the most probable reaction path in a theoretical mixture. Table 5 shows the computed Boltzmann distribution based on the free energies of the reactant structures in the first reaction coordinate. The alpha site reactant in the doublet potential energy surface is the most stable structure, responding by about half of the conformation states in the distribution. The other half is mainly composed of delta cis/trans doublet and delta cis quartet conformers. Due to the complex electronic nature of the transition metal-containing systems, the reaction can proceed through multiple potential energy surfaces. As the structures in doublet and quartet reaction paths are very similar, nearly as degenerate states or with very close energies, we can assume that spin crossover can occur through these spin states. Indeed, experimental and theoretical evidence suggests facile interconversion between these spin states in iron porphyrin compounds [12,53]. Lastly, the calculated Boltzmann distribution for the product conformers of the Habstraction (coordinate III) predicts approximately a 50:50 distribution around cis/trans conformers in the doublet state (Table 6). Using the Eyring equation, we can estimate the reaction through the delta trans doublet path as three and two times faster than the delta cis quartet and doublet paths, respectively. Despite the initially lower concentration in the reaction coordinate I for the trans doublet conformer, the lower activation barrier suggests greater kinetic favor to the reaction through this path. Thus, it is expected that β-pinene CYP-catalyzed hydrogen abstraction led to the formation of the doublet trans radical conformer as the major reaction product. After the formation of the cis/trans doublet radical, the rebound mechanism can form the highly stable hydroxylated product. The formation of these structures is accompanied by the release of substantial amount of energy. The total Gibbs free energies released by the formation of doublet cis and trans hydroxylated products can be estimated as approximately −48 kcal/mol. Reversal of the hydroxylation reaction demands overcoming barriers as large as 40 kcal/mol. Thus, our results suggest that the hydroxylated product formation acts as a "thermodynamic trap" for the reaction system since the reverse reaction is kinetically hindered.

CYP-Catalyzed Hydroxylation of α-Pinene
Considering the α-pinene molecule, we studied the CYP-catalyzed hydrogen abstraction followed by hydroxylation of the epsilon, alpha, and gamma allylic sites in the doublet, quartet, and sextet spin (for epsilon only) states. Figures 6 and 7 (and Figure S3 and S4) show the calculated reaction coordinates. Table 7 (and Table S11, S12 and S14) presents the calculated absolute energies and the kinetic and thermodynamic parameters. As in the case of the β-pinene delta-site, the epsilon-site has two cis/trans hydrogen atoms (regarding the η site), liable for being abstracted by the CYP. The R and S configurations also had similar electronic energies. The cis transition states of hydrogen abstraction (coordinate II) had slightly higher energies than the trans transition states, except for the doublet state. The sextet state potential energy surface (PES) lies above the doublet and quartet PES. We thus excluded sextet PES from calculations at alpha and gamma sites. The reactions through the epsilon and gamma sites had favorable thermodynamic profiles and low activation barriers, whereas the reaction at the alpha site was found to be thermodynamically unfavorable. As the β-pinene, the application of vibrational corrections to the electronic energy led to a decrease in the enthalpic activation energies (Table 7) for the α-pinene reaction. The effect in the reaction enthalpy is less accentuated, but trends to an increase in the liberated energy. Except for the doublet epsilon cis conformation, the entropic effect led to an increase in the activation energy values regarding the enthalpy of activation. The results of the empirical corrections added in our DFT/LAN methodology over the pure DFT repeat the trend observed in the β-pinene, resulting in better thermodynamic and kinetic profiles for the CYP-catalyzed hydrogen transfer reaction.
tronic energy led to a decrease in the enthalpic activation energies (Table 7) for the αpinene reaction. The effect in the reaction enthalpy is less accentuated, but trends to an increase in the liberated energy. Except for the doublet epsilon cis conformation, the entropic effect led to an increase in the activation energy values regarding the enthalpy of activation. The results of the empirical corrections added in our DFT/LAN methodology over the pure DFT repeat the trend observed in the β-pinene, resulting in better thermodynamic and kinetic profiles for the CYP-catalyzed hydrogen transfer reaction.   Our calculations at B3LYP/LAN level showed that the most stable products were formed in the doublet state by hydrogen abstraction at the epsilon site in agreement with the calculations of the bare system. The trans product is slightly more stable; however, the reaction through the cis path has a lower activation barrier (9.00 vs. 12.2 kcal/mol) and lower Gibbs free energy (−15 vs. −13.8 kcal/mol).
As in the case of β-pinene, the transition state for the rebound mechanism was not identified. So, we consider that this mechanism is barrierless at our theoretical level with the hydrogen abstraction by the CYP enzyme as the rate-limiting step. The most stable hydroxylated rebounded product was formed at the alpha site followed by the epsilon trans product, which was only 0.7 kcal/mol less stable. The same order of the allylic products (alpha > epsilon > gamma) was obtained as in the bare system. However, the difference between alpha and epsilon is only 1 kcal/mol rather than those 5-6 kcal/mol obtained from the bare system. This indicates that the presence of CYP450 and the rebounded structure formation linked to the iron center slightly shifted the equilibrium in relation to the bare model. Table 7. Thermodynamic and kinetic data computed for the hydrogen capture reaction by CYP in α-pinene. The terms ΔrE, ΔrH, and ΔrG refer to electronic, enthalpic, and Gibbs free energy variation in the reaction, respectively. The terms ΔE ‡ , ΔH ‡ , and ΔG ‡ refer to the electronic, enthalpic, and Gibbs free energy reaction barriers, respectively, and the term "cor" refers to the dispersion, BSSE, A B Our calculations at B3LYP/LAN level showed that the most stable products were formed in the doublet state by hydrogen abstraction at the epsilon site in agreement with the calculations of the bare system. The trans product is slightly more stable; however, the reaction through the cis path has a lower activation barrier (9.00 vs. 12.2 kcal/mol) and lower Gibbs free energy (−15 vs. −13.8 kcal/mol).
As in the case of β-pinene, the transition state for the rebound mechanism was not identified. So, we consider that this mechanism is barrierless at our theoretical level with the hydrogen abstraction by the CYP enzyme as the rate-limiting step. The most stable hydroxylated rebounded product was formed at the alpha site followed by the epsilon trans product, which was only 0.7 kcal/mol less stable. The same order of the allylic products (alpha > epsilon > gamma) was obtained as in the bare system. However, the difference between alpha and epsilon is only 1 kcal/mol rather than those 5-6 kcal/mol obtained from the bare system. This indicates that the presence of CYP450 and the rebounded structure formation linked to the iron center slightly shifted the equilibrium in relation to the bare model. Table 7. Thermodynamic and kinetic data computed for the hydrogen capture reaction by CYP in α-pinene. The terms ∆ r E, ∆ r H, and ∆ r G refer to electronic, enthalpic, and Gibbs free energy variation in the reaction, respectively. The terms ∆E ‡ , ∆H ‡ , and ∆G ‡ refer to the electronic, enthalpic, and Gibbs free energy reaction barriers, respectively, and the term "cor" refers to the dispersion, BSSE, and Grimme quasi-harmonic corrections to the calculated absolute Gibbs free energies. Hyd ∆ r G cor refers to the hydroxylation free energies computed as the difference between the hydroxylated product and the reactant in coordinate 1. Data are displayed in kcal/mol.

Site
Spin With this data, we can reason the possible paths that the reaction can proceed in this theoretical scenario. The Boltzmann distribution (Table 8) related to the structures in coordinate I shows the existence of a mixture of conformers in the initial reaction mass. There is no hydrogen transfer reaction through the sextet spin state and the reaction through the alpha site is thermodynamically unfavorable. The reactions through the remaining sites have relatively low barriers and are thermodynamically feasible. However, our calculations (Table 9) showed that the most stable radicals are formed at the epsilon cis and trans path in the doublet state. The Boltzmann distribution indicates that the equilibrium is totally shifted to the formation of these products. From these structures, the hydroxylation reaction of α-pinene can proceed fast with a substantial release of~50 kcal/mol, a little more than the ones computed for β-pinene. To get more insight into the hydrogen transfer to the CYP model reaction, we conducted an electronic analysis of the α-pinene. We analyzed the Mulliken and NPA chargers, as well as the Mulliken and NBO spin population in doublet and quartet states for epsilon cis and trans paths (Table S12). The orbital analysis shows that the lowest orbitals shown are the bonding type orbitals along the Fe-O bond and include the σz 2 for the overlap of 3dz 2 on iron with 2pz on oxygen and the degenerate pair of πxz and πyz molecular orbitals for the bonding interaction between the 3dxz (or 3dyz) on iron with the 2px (or 2py) atomic orbital on oxygen [45,54]. The antibonding combination of this pair of orbital π*xz reflects the antibonding interactions of the 3dxz orbitals on Fe with the 2py orbitals on O along the Fe-O bond [54]. The δx 2 -y 2 orbital is non-bonding and located in the plane of the heme/porphyrin. The two σ* antibonding orbitals are high in energy and virtual-one along the O-Fe-S axis (the z-axis) [2,55], namely σ*z2, is seen in Figure S4.

Materials and Methods
All DFT calculations were conducted in the Gaussian09 software (dftd4 version 3.2.0 from GitHub) [56] with the well-known three-parameter hybrid functional of Becke (B3) [57,58] in combination with the Lee-Yang-Parr (LYP) [59,60] correlation functional. The B3LYP hybrid functional with unrestricted formalism has been largely used in the theoretical studies for modeling active sites and analyzing reaction mechanisms in many systems, particularly CYP [3,61] and Heme [62] models of enzymes, where it was shown to give correct spin-state energetics for these species [2,4,38]. Studies have shown that this functional can provide accurate geometries and reasonable agreement with Coupled Cluster Singles and Doubles Theory (CCSD(T)) calculations of the relative spin states energies for porphyrin [63] complexes and heme-related [62,64] models. For geometry optimizations, vibrational analyses, and Intrinsic Reaction Coordinate (IRC) calculations, we employed the Pople double-zeta polarized 6-31G(d) basis set [65,66] for all atoms, except iron. In this atom, the orbitals are expanded with the LANL2DZ [67][68][69] ECP basis set. The default Gaussian algorithm (Berny algorithm using GEDIIS) was adopted to minimize energy in the minimum and maximum stationary points. In the optimization process, we consider that structural convergence occurs only when the maximum force and its root mean square (RMS) fall below 4.5 × 10 −4 and 3 × 10 −4 hartree/bohr. We fix a threshold of 1.8 × 10 −3 and 4 × 10 −3 bohr for the maximum displacement, and its RMS with the default grid for the numerical integration step in all calculations. Low-lying harmonic vibrational mode's contribution to entropy frequencies obtained from calculations were treated using a free-rotor approximation [70] (quasi-rigid-rotor-harmonic-oscillator) implemented in the GoodVibes software (GoodVibes version 3.1.1 from GitHub) [49]. To address the basis set superposition error (BSSE), we added Grimme geometrical counterpoise gCP [71,72] correction for all computed energies. Dispersion effects may play a relevant role in CYP-like systems [73,74]. Therefore, to assure the correct dispersion computation, we added the D4 Grimme's empirical corrections [75,76] for all computed energies via external software (Gcp version 2.3.1 from GitHub). The Polarizable Continuum Model (PCM) [77,78] was used to simulate the water effect in this molecular system. This computational methodology will be referred to as B3LYP/LAN in the manuscript.

Conclusions
To conclude, we report a systematic mechanistic study of biocatalytic hydrogen abstraction and hydroxylation of αand β-pinenes by multi-state cytochrome P450 enzymes (CYP) based on the density functional theory (DFT) method using the Gaussian09 software and B3LYP/LAN computational methodology. CYPs are an interesting class of biocatalysts as they can generate several products from the same substrate due to the presence of multi-states (different spin states). In our DFT study, we started by studying the hydrogen abstraction and hydroxylation of α-pinene and β-pinene without a CYP450 model using the B3LYP/LAN methodology (bare system) as a reference to the subsequent CYP-catalyzed hydroxylation of β-pinene and CYP-catalyzed hydroxylation of α-pinene accounting for the doublet, quartet, and sextet spin states of CYP. According to the potential energy surface and Boltzmann distribution for radical conformers, the major reaction products of CYPcatalyzed hydrogen abstraction from β-pinene were the doublet trans (53.4%) and doublet cis (46.1%) radical conformers at the delta site, and the formation of doublet cis/trans hydroxylated products released a total Gibbs free energy of~48 kcal/mol. As for α-pinene, the most stable radicals were trans-doublet (86.4%) and cis-doublet (13.6%) at epsilon sites, and their hydroxylation products released a total of~50 kcal/mol Gibbs free energy. The transition state connecting the radical (from hydrogen abstraction) and the hydroxylated product (from oxygen rebound) was not located. The rebound mechanism is considered barrierless at our mentioned theoretical level, in line with previous studies in the literature, and the hydrogen abstraction by CYP is the rate-limiting step.  Data Availability Statement: The raw and processed data required to reproduce the above findings are available in the manuscript and the Supplementary Materials.

Conflicts of Interest:
The authors declare no conflict of interest.