To identify a plausible binding structure of PTX in CYP2C8, we initially performed molecular docking simulations. We utilized MOE 2020 software (Molecular Operating Environment (MOE), 2022) and docked PTX into several crystal structures of CYP2C8 (PDB codes 2NNH, 2NNI, 2VN0, and 2NNJ). Unlike CYP3A4, these crystal structures exhibit minimal structural variations, leading to comparable docking results. Consequently, we proceeded with the structure from 2NNJ. To further refine the CYP2C8–PTX complex structure, we performed molecular dynamics (MD) simulations using AMBER (Case et al., 2020). The MD simulations essentially followed the same protocol employed in our prior study of CYP3A4 (Yue and Hirao, 2023), including clustering to select a representative structure.

To investigate the reaction mechanism, we employed the ONIOM(QM:MM) method, a subtractive QM/MM scheme (Chung et al., 2012; Chung et al., 2015). Gaussian 16 software was used for QM/MM computations (Frisch et al., 2016). The CYP2C8–PTX complex geometry, derived from MD simulations, served as the starting point. To simplify the system, a significant portion of the solvent water molecules outside the enzyme was removed. The QM region encompassed the porphine and iron-oxo units, the C_βH₃S⁻segment of the cysteine axial ligand, and the PTX molecule. DFT methods were applied to describe the QM atoms. Geometry optimization and vibrational frequency calculations were performed using the ONIOM mechanical-embedding scheme and the B3LYP/6-31G(d) QM method. The ONIOM energy obtained with this basis set is referred to as E1. For single-point energy calculations, the B3LYP/def2-TZVP method and the electronic-embedding scheme were utilized (Becke, 1993; Lee et al., 1988; Vosko et al., 1980; Hehre et al., 1972; Weigend and Ahlrichs, 2005). MM calculations within the ONIOM framework employed AMBER and TIP3P parameters (Cornell et al., 1995; Jorgensen et al., 1983). Force field parameters for PTX and the heme group were consistent with those in the MD simulations. Free energies were determined by summing the ONIOM single-point energy with a larger basis set (E2), the DFT-D3BJ dispersion correction to the QM energy (E _disp) (Grimme et al., 2011; Becke and Johnson, 2005), and the free energy correction obtained from ONIOM vibrational frequency analysis (G _corr). The sum of these energy values (G) was subsequently used to construct reaction energy diagrams.

Human P450s generally possess spacious active sites, enabling them to accommodate a wide range of substrates. CYP2C8 (PDB code 2NNJ), with a 746 Å³ active site, readily binds PTX (Figure 1A). However, an α-helix overlying the heme restricts ligand access to this region. This structural feature is also observed in CYP2C9, which shares 78% sequence identity with CYP2C8 (Williams et al., 2003). As a result, molecular docking identified a binding mode (Pose A) avoiding the α-helix region and positioning the metabolic site’s H_β closer to the heme iron compared to H_α (Figure 1B). At the C6 position of PTX in Scheme 2, H_α and H_β are potential sites for hydroxylation. While C–H_α hydroxylation would produce the experimentally observed 6α-hydroxypaclitaxel metabolite (Harris et al., 1994a; Kumar et al., 1994), C–H_β hydroxylation would yield the unobserved 6β-hydroxypaclitaxel. Thus, Pose A is inconsistent with experimental findings. Further analysis revealed that the Pose A conformation of PTX is not very stable in its isolated form, which could lead to an overestimation of its binding affinity in docking simulations. Consequently, a lower-ranked binding mode, involving a substrate with an intrinsically more stable conformation, should be a more probable candidate for C–H_α hydroxylation (Yue and Hirao, submitted). Nevertheless, in this study, we delved deeper into the reactivity of Pose A, which unexpectedly led us to uncover intriguing insights into high-spin reactivity.

Following MD refinement of the Pose A structure, we conducted QM/MM mechanistic studies. P450-catalyzed alkane hydroxylation typically initiates with H-abstraction from a C–H bond (Schöneboom et al., 2004). As illustrated in Figure 2A, the geometry-optimized pre-reaction complex between Cpd I and PTX (1) in the doublet ground state (² 1) exhibits an O–H distance of 2.28 Å between the H_β atom and the Cpd I oxo group. This proximity suggests a favorable conformation for subsequent H-abstraction. Figure 2B highlights key interactions between PTX and surrounding amino-acid residues. PTX is securely anchored within the CYP2C8 active site via non-bonded interactions including hydrogen bonds and π-π interactions. Specifically, Asn99, Ser103, Asn209, and Gly365 form hydrogen bonds with PTX’s polar groups. Moreover, a hydrogen-bond network involving Asn99, Ser100, Gln214, Ser103, Leu208, and Asn209 stabilizes the Cpd I–PTX complex. Phe201 and Phe205 engage in hydrophobic interactions with the phenyl moiety of the benzoyloxy group in PTX. Figure 2C depicts the available space around PTX using a surface model, clearly demonstrating that PTX’s excellent fit within the active site.

Through detailed QM/MM calculations on 6β-hydroxypaclitaxel formation pathways, we identified four distinct electron-shift patterns for the H-abstraction step. One doublet-state pathway (Path A) involves an electron shift from the substrate orbital (ϕ) to the a_2u-type orbital of Cpd I, yielding intermediate ² 3a with a negative spin density (ρ) value (∼–1.0) on the PTX moiety (Scheme 1B; Supplementary Table S2). Another doublet-state pathway (Path B) features an electron shift from ϕ to the π* orbital of Cpd I, forming intermediate ² 3b with a positive ρ value (∼1.0) on PTX (Supplementary Table S2). In the quartet spin state, an electron migrates from the substrate to the a_2u-type orbital of Cpd I, producing intermediate ⁴ 3 with a positive ρ value (∼1.0) on PTX (Supplementary Table S2). Finally, the sextet spin state involves electron migration from ϕ to the σ*_z 2 orbital of Cpd I, resulting in intermediate ⁶ 3 with a negative ρ value (∼–1.0) on PTX (Supplementary Table S2).

Figure 3 presents the free energy profile for PTX hydroxylation. As described, two electron-shift patterns can occur during the doublet-state H-abstraction step, resulting in two TSs (² 2a* and ² 2b*) for Paths A and B, respectively. TSs were also obtained in the quartet (⁴ 2*) and sextet (⁶ 2*) states. Comparing relative energies reveals ² 2b* as slightly more stable (25.6 kcal/mol) than ² 2a* (26.7 kcal/mol). Without corrections, the uncorrected E1 values produce barriers exceeding 30 kcal/mol for both doublet states (Supplementary Table S1), significantly surpassing typically observed values of around 20 kcal/mol (Schöneboom et al., 2004). These results suggest that the substrate experiences significant strain around the TS in the protein environment. Intermediates ² 3a and ² 3b exhibit reversed stability (5.3 and 9.5 kal/mol, respectively). The stability of the quartet-state TS (⁴ 2*, 25.4 kcal/mol) is comparable to that of ² 2b*, with subsequent intermediate ⁴ 3 at 5.6 kcal/mol.

Surprisingly, the sextet-state TS (⁶ 2*, 19.4 kcal/mol) is significantly more stable than the others, with energy differences of at least 6 kcal/mol. This exceptional stability arises from a remarkably low barrier of 11.3 kcal/mol on the sextet-state free energy surface, contrasting sharply with the higher barriers observed for the other spin states. While the sextet state typically benefits from transition-state stabilization due to exchange enhancement (Hirao et al., 2005), the substantially lower energy of the sextet TS in the current system is an unusual observation. Normally, the sextet state’s stability at the TS is comparable to or slightly inferior to that of lower spin states. In contrast, the CYP3A4-catalyzed aromatic hydroxylation reaction exhibited a slight energetic preference for the sextet TS, with a small energy gap of 0.3 kcal/mol (Yue and Hirao, 2023). The significantly greater stability of the sextet TS in the current CYP2C8-catalyzed reaction is expected to increase the likelihood of this reaction channel being accessed.

Figure 4 displays the optimized geometries of key TSs and intermediates. The larger O–H distance observed in ⁶ 2* (1.25 Å) than those in the other corresponding species (1.16–1.24 Å) indicates an earlier TS in the sextet pathway, consistent with the lower energy barrier observed in the sextet-state energy profile (Figure 3). A closer examination of Fe–O–H angles across different spin states reveals a larger angle for ⁶ 2* (135.1°) compared to the others (around 127.0°). This larger Fe–O–H angle in the sextet state is characteristic of the high-spin H-abstraction mechanism, involving an electron shift from the substrate’s ϕ orbital to the axial σ*_z 2 orbital of iron(IV)-oxo. Optimal orbital overlap is achieved when the substrate approaches the iron (IV)-oxo unit from above, leading to a larger Fe–O–H angle (Hirao et al., 2005).

We further investigated the rebound step leading to cycloalkanol product formation (Figures 3, 4) (Groves, 1985). Starting from intermediates ² 3a and ³ 3b in the doublet spin state, we obtained product ² 5 (−37.0 kcal/mol) via TSs ² 4a* (16.7 kcal/mol) and ² 4b* (18.8 kcal/mol), respectively. While Path A exhibits a lower overall energy, its rebound barrier of 11.3 kcal/mol is unusually high for a P450 reaction in the doublet state compared to other reported values, likely due to steric constraints within the enzyme’s active site. The quartet TS ⁴ 4* (17.9 kcal/mol) lies energetically between ² 4a* and ² 4b*, leading to product ⁴ 5 with a relative free energy of −36.7 kcal/mol, which is less stable than the doublet product (² 5, −37.0 kcal/mol). From the most stable intermediate ⁶ 3 (4.6 kcal/mol), we obtained the most stable product ⁶ 5 (−39.8 kcal/mol) via TS ⁶ 4* (13.2 kcal/mol). Surprisingly, the sextet state exhibits the lowest rebound barrier in the sextet state (8.6 kcal/mol), contrary to the typical trend of higher rebound barriers in the sextet state.

Our QM/MM mechanistic study on CYP2C8-catalyzed PTX hydroxylation unexpectedly revealed a dominant role for the high-spin sextet state from an early stage of the reaction. To understand the origins of this pronounced sextet state stability, we conducted further theoretical analyses. The ONIOM-based G value can be expressed as follows (Equation 1): G = E QM + E disp + E MM + E pol + G corr (1) where E _QM is the gas-phase energy of the QM atoms at the QM/MM geometry, and E _MM and E _pol are defined by Equations 2, 3, respectively: E MM = E MM , real   –   E MM , model (2) E pol = E 2 EE   –   E 2 ME (3)

E2(EE) and E2(ME) are ONIOM-EE and ONIOM-ME energies obtained from single-point energy calculations with the def2-TZVP basis set (Hirao, 2011a; Hirao, 2011b). The relative free energy (ΔG) of the H-abstraction TS with respect to ² 1 is given by Equation 4: Δ G = Δ E QM + Δ E disp + Δ E MM + Δ E pol + Δ G corr (4)

While ∆E _QM can be obtained from DFT calculations on the QM atoms, the remaining four terms also contribute to the overall ΔG value. Figure 5 displays the latter four terms for the H-abstraction TSs in different spin states. The consistently negative ∆E _disp values indicate stronger dispersion stabilization in TSs compared to ² 1, with the largest stabilization observed for the sextet state. Therefore, although modest, dispersion contributes to the sextet state’s stability. The relatively large ∆E _disp for the sextet TS is attributed to appropriate positioning of the QM atoms, which enhances interatomic dispersion stabilization. The MM energy term (∆E _MM) is also consistently negative, with the largest stabilization for ⁶ 2*. Therefore, the MM effect partly contributes to the pronounced stability of ⁶ 2*. While the polarization effect (∆E _pol) does not enhance the sextet TS’s stability, the ∆G _corr term stabilizes all TSs, with the greatest stabilization observed for the sextet state. This likely results from the sextet TS’s flexible and entropically favorable structure. Among the four terms examined, the free energy correction term (∆G _corr) provides the most significant stabilization for ⁶ 2*.

Despite the insights gained from the above analysis, the substantial relative stability of ⁶ 2* remained only partially explained. Therefore, we further investigated the ΔE _QM values across different spin states. In addition to the QM/MM approach, DFT calculations were performed using a P450 Cpd I model and cyclohexane to represent the reaction. Figure 6 compares the ΔE _QM values for different spin states obtained from both methods. Interestingly, significant discrepancies were observed between the QM/MM and DFT results. QM/MM calculations consistently yielded higher ΔE _QM by 9.1–11.9 kcal/mol compared to DFT. In other words, the TS geometries are highly strained in the protein environment. However, the smallest destabilization was observed for the sextet state (9.1 kcal/mol). Thus, the protein environment destabilizes TSs more in the other spin states, contributing to the enhanced relative stability of ⁶ 2*. The smaller destabilization in the sextet state is likely due to the protein-imposed geometric constraints that orient the substrate more favorably for the sextet-state TS geometry.

Based on the present theoretical analysis and previous studies, several possible strategies for promoting high-spin reactivity of iron(IV)-oxo complexes can be summarized (Scheme 3). The most fundamental approach involves stabilizing the high-spin pre-reaction complex, which lowers the energy of the entire high-spin energy surface (Scheme 3A). This can be achieved through the use of nonheme ligands that reduce ligand field strength. This study highlights additional factors beyond this approach (Scheme 3B-E), with entropy and QM effects being particularly influential (Scheme 3D,E). We propose that these effects are key to enhancing high-spin reactivity. While pronounced entropy benefits for high-spin TSs have been observed in previous studies (Yue and Hirao, 2023), further enhancing the high-spin reactivity requires increasing the quantum mechanical stability of the high-spin TS relative to others. The surrounding environment plays a crucial role in influencing this effect, as the QM effect is largely related to differing degrees of geometric strain of the substrate. If the environment imposes steric constraints on the substrate, enforcing a linear approach toward the iron(IV)-oxo unit from above, a lower degree of destabilization in the high-spin TS compared to others can be expected, as the high-spin state can undergo an efficient electron shift in an exchange enhanced fashion within this configuration. Therefore, even in nonenzymatic environments, ligand design should consider such steric effects.