Structure-based modeling reveals molecular basis for CYP153A6’s novel activity toward toluene derivatives
Yao Wei, Silvia Donzella, Sara Foiadelli, Francesco Molinari, Uliano Guerrini, Ivano Eberini
TL;DR
This study explores how a bacterial enzyme selectively hydroxylates toluene derivatives and identifies structural and molecular factors influencing its activity.
Contribution
The study reveals the molecular basis for CYP153A6's activity toward toluene derivatives using structure-based modeling and computational methods.
Findings
CYP153A6 efficiently converts toluene derivatives with apolar or slightly polar substituents.
Hydroxyl or nitro groups on the aromatic ring reduce binding affinity and hinder catalytic activity.
Computational simulations confirm the hydrophobic nature of the enzyme's active site.
Abstract
CYP153A6 from Mycobacterium sp. strain HXN-1500 is a monooxygenase that catalyzes the selective hydroxylation of terminal methyl groups in various alkanes. This study explored CYP153A6 activity towards a set of toluene derivatives, along with the underlying molecular recognition. Initial in vivo evaluation of CYP153A6 activity, conducted using both whole cells and cell-free extracts, showed efficient conversion of toluene derivatives with apolar or slightly polar substituents, while no detectable activity was observed for derivatives bearing more polar groups. A homology model of CYP153A6 3D structure was built and validated, revealing key structural features and molecular tunnels. An ensemble docking strategy was used to identify the most effective docking setup. Molecular dynamics simulations and binding free energy calculations further confirmed the hydrophobic nature of the active…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 10
Figure 11
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 12
Figure 13
Figure 14- —https://doi.org/10.13039/100018694HORIZON EUROPE Marie Sklodowska-Curie Actions
- —https://doi.org/10.13039/501100021856Ministero dell’Università e della Ricerca
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsMetal-Catalyzed Oxygenation Mechanisms · Pharmacogenetics and Drug Metabolism · Microbial bioremediation and biosurfactants
Introduction
Selective hydroxylation of methyl groups is a highly valuable transformation in synthetic chemistry with particular attention to green chemistry. This process enables the functionalization of the inert C-H bonds, which is critical for enhancing chemical solubility, and generating novel molecular frameworks^1,2^. Environmentally, hydroxylation of hydrophobic hydrocarbons facilitates detoxification and biodegradation of industrial pollutants, particularly toluene derivatives frequently encountered in fuels, solvents, and chemical manufacturing^3^.
The selective hydroxylation of benzylic positions in toluene derivatives represents a critical transformation in synthetic and pharmaceutical chemistry, offering a sustainable alternative to traditional chemical oxidation methods^4^. Enzymes, particularly cytochrome P450 (CYP450) monooxygenases and flavin-dependent oxygenases, have emerged as powerful biocatalysts for these reactions due to their regio- and stereo-selectivity under mild conditions. A notable example is that CYP102A1 converts ethylbenzene to (S)−1-phenylethanol with high enantioselectivity (95% ee)^5^. Similarly, 4-hydroxyphenylpyruvate dioxygenase (HPPD) catalyzes the hydroxylation of 4-hydroxyphenylpyruvate to homogentisate via benzylic oxidation^6^. Non-heme iron oxygenases can also hydroxylate benzylic carbons in aromatic compounds such as ethylbenzene, yielding (S)−1-phenylethanol^7^. Flavin-dependent monooxygenases, including styrene monooxygenase, oxidize styrene derivatives at the benzylic position to form epoxides oralcohols^8^. Together, these approaches underscore a broader shift toward environmentally benign and functionally selective hydroxylation methods – particularly relevant for the modification of pharmaceuticals and natural products.
CYP153A monooxygenases belong to a distinct subfamily of cytochrome P450 enzymes (CYPs) and are primarily involved in the oxidation of aliphatic and aromatic hydrocarbons, with a strong preference for terminal hydroxylation of n-alkanes and similar substrates^9,10^. Among them, CYP153A6 from Mycobacterium sp. HXN-1500 is a soluble monooxygenase that exhibits selective activity toward terminal methyl groups in a variety of linear alkanes and alicyclic substrates^11,12^. Notably, it selectively hydroxylates various monoterpene derivatives with high regioselectivity at the allylic methyl group attached to the aliphatic ring. This catalytic feature enables the efficient conversion of limonene into perillyl alcohol, a biotransformation of significant industrial relevance^13^. In an earlier study, CYP153A6 was shown to exhibit activity toward p-cymene, although the resulting product was not characterized; this observation suggests a potentially broader substrate scope, including the ability to target benzylic methyl groups. The selective oxidation of benzylic positions remains a significant challenge in both chemical and enzymatic catalysis. Benzylic hydroxylation offers a valuable route to benzyl alcohols, which are important intermediates in the pharmaceutical, fragrance, and fine chemical industries^14,15^. Given its inherent regioselectivity and terminal methyl preference, CYP153A6 emerges as a potential promising candidate for the selective hydroxylation of methyl substituents on the aromatic rings.
In this study, we investigate the catalytic potential of CYP153A6 for benzylic hydroxylation in series of para-, meta-, and ortho-substituted derivatives. Preliminary experimental work was carried out using whole-cell biocatalysis and crude cell-free extracts of recombinant E. coli expressing CYP153A6. To guide the future protein engineering efforts, we adopted a set of computational modeling strategies. This included homology modeling, molecular dynamics (MD) simulations, substrate tunnel calculations, ensemble docking, protein-ligand interaction fingerprints analysis, MM-GBSA energy analysis, and QM/MM calculations. The aim of this study is to elucidate the CYP153A6 molecular recognition mechanism and provide insights into the toluene derivatives binding modes and the driven potential. This will lay the groundwork for expanding the catalytic ability of CYP153A6 toward high-value aromatic compounds through selective toluene derivative hydroxylation.
Results and discussion
Activity on toluene derivatives
Preliminary experiments were performed using whole cells (WC) of recombinant CYP153A6-E. coli and crude cell free extract (CFE). We first evaluated the hydroxylation of p-cymene (1a) using a nominal substrate concentration of 6 mM; reactions were carried out at a formal substrate concentration of 6 mM, exceeding the aqueous solubility of the compound (Fig. 1).
Fig. 1. Biotransformation of p-cymene using WC and CFE (concentration of CYP153A6 1.8 µM) of E. coli at pH 7.0, 6 mM substrate, 28 °C. Reaction with CFE was carried out in the presence of 6 mM NADH.
p-Cymene underwent highly selective hydroxylation to yield p-isopropylbenzyl alcohol as the only detectable product, with conversions of 85% using whole cells and 57% using the cell-free extract (CFE) after 8 h. Kinetic parameters were measured using CFE containing substrate concentrations in the range of 0.5–10.0 mM, showing that CFEs were used in further experiments where the activity of CYP153A6 was investigated using different p-substituted toluene derivatives (1b–1h, Table 1). Because these substrates are highly hydrophobic, the kinetic constants are apparent values calculated from nominal substrate concentrations, potentially leading to an underestimation of the intrinsic enzymatic affinity. Hydroxylation occurred exclusively at the benzylic position, leading to the formation of the corresponding benzyl alcohols. The most reactive substrate was p-methylanisole (1d), which showed the highest specific activity (1.75 U mg⁻¹) and, being the most water-soluble among the active substrates, exhibited the lowest apparent K_m_. Although the activity on toluene (1b) was sluggish, hydroxylation of p-xylene (1c) proceeded at a good rate; notably, only one of the two benzylic methyl groups of 1c was hydroxylated, even after prolonged incubation, highlighting the enzyme’s remarkable selectivity. Good conversion was also observed for p-chlorotoluene (1f), while p-cresol, p-methyl-benzylalcohol, and p-nitrotoluene showed no detectable activity. Because all these substrates are highly hydrophobic, the kinetic constants represent apparent values calculated from nominal substrate concentrations, and the elevated apparent K_m_ values are consistent with the limited aqueous solubility of most of the substrates. These experimental results suggest that CYP153A6 activity was strongly influenced by substituent polarity, favoring apolar or slightly polar toluene derivatives over more polar analogues.
Table 1. Specific activity and kinetic parameters of CYP153A6 (1.8 µM) towards p-substituted toluene derivatives. Specific activity was measured using CFE at pH 7.0, 6 mM substrate, 28 °C.^a^ n.d. = not determined.
CYP153A6 was also employed for the oxidation of a series of meta- and ortho-substituted toluene derivatives (1i–1o, Table 2); the substrates were chosen to have meaningful indication about the influence of the substituent position on the overall activity.
Table 2. Specific activity and kinetic parameters of CYP153A6 (1.8 µM) towards p-substituted toluene derivatives. Specific activity was measured using CFE at pH 7.0, 6 mM substrate, 28 °C.
Substrates contain a hydroxyl group, either as a phenol (1i and 1j) or a primary alcohol (1m and 1n), exhibited minimal or undetectable activity, with a trend like that observed for the p-substituted toluene analogues. The activity observed across the series of xylene derivatives suggests that steric hindrance plays a critical role: o-xylene, due to its steric bulk, is unreactive, whereas p- and m-xylene exhibit comparable levels of reactivity.
Homology modeling of CYP153A6 3D structure
CYP153A6 structure has not yet been experimentally solved. Therefore, its 3D structure (Fig. 2 (a)) was modeled basing on the experimentally solved structure of CYP153A_P.sP_ (PDB ID 7AO7) by X-ray diffraction^16^. This homologous template enzyme shares similar functionality with CYP153A6, and sequence alignment revealed a global identity of 78.6% and a similarity of 86.5%.
Upon superimposing the CYP153A6 and CYP153A_P.sP_ structures, the root mean square deviation (RMSD) of the Cα atoms were calculated to be smaller than 2 Å, indicating a good structural overlap. The Ramachandran plot (Fig. 2 (b)) further supports the high quality of the CYP153A6 homology model, with most residues occupying favorable conformational regions and only a few outliers. Notably, the plot closely resembles that of the experimentally solved template structure, in which a similar number of outliers is also observed. These results suggest that the model is of high structural quality^17^.
Fig. 2**(a)** CYP153A6 homology model: Cofactor heme is shown in gray sticks, 11 α-helices are labeled with A-K, 5 β-sheets are labeled with β1–5. (b) Ramachandran plot of the CYP153A6 homology model. (Figures were produced with the Schrödinger Maestro Suite (release 2024-2)^18^, and protein secondary structures were labeled using Diagram program: https://app.diagrams.net/).
To assess the structural stability of the CYP153A6 homology model, classical molecular dynamics (MD) simulations were performed. The RMSD profile of the protein (Fig. 3), calculated on the Cα atoms of the amino acids, stabilized around ~ 3Å, indicating that the backbone structure reached equilibrium during the simulation. This suggests that the homology model maintains a reasonable overall fold, providing confidence in the generated 3D structure for subsequent analyses. No ligands were included in this simulation, and the primary goal was to verify the structural reliability of the protein model.
Fig. 3Cα RMSD values of CYP153A6 from three replicate MD simulations.
Molecular tunnels contributes to enzyme activity
The active site of CYP450s is deeply buried within the protein core^19^. These enzymes possess molecular tunnels that connect the active site to the surrounding solvent environment, allowing substrate access and product release. These tunnels play a crucial role in determining substrate specificity^20^. Typically, tunnels involved in the uptake of hydrophobic, membrane-associated substrates are oriented toward the membrane, while those associated with the release of hydrophilic molecules – such as oxidized products or water – extend into the solvent^21^.
Tunnel analysis of CYP153A6 revealed three distinct molecular tunnels (Fig. 4): two are likely substrate access tunnels (designated as tunnels 2b and 2c), and one is predicted to work as a water tunnel (designated as tunnel W). Substrate tunnel 2b has a throughput score of 0.48 and a bottleneck radius of 1.15 Å. It involves residues GLY86, GLY87, ILE88, THR89, ALA96, ALA97, SER98, LEU99, MET101, ILE103, and ALA104, with ALA104 identified as the bottleneck residue. Substrate tunnel 2c and water tunnel W share the same set of lining residues: THR89, MET91, ASN94, ALA97, THR186, THR187, ALA188, SER405, PHE407, and VAL408, with MET91 acting as bottleneck residue in both tunnels. The throughput score values for tunnels 2c and W are 0.34 and 0.31, respectively, and both a bottleneck radius of 1.00 Å. To in-silico validate tunnel W as a water pathway, water molecules were generated using the three-dimensional reference interaction site model (3D-RISM). This method provides time-averaged distribution of hydrogen and oxygen densities in water, as well as free-energy maps that help assess solvent stability and solvation contributions to binding free-energy^22^. All generated water molecules were located within the predicted tunnel (Fig. 5), supporting its functional consequence.
Fig. 4. Detected molecular tunnels in CYP153A6. The paths of tunnels 2b (green), 2c (orange), W (cyan) are depicted. Enzyme cofactor heme is shown in magentas sticks, FG loop is shown in blue, BC loop is shown in red (Protein structure was displayed by PyMOL^23^, molecular tunnels were generated using CAVER 3.03 plugin software^24^, protein secondary structures were labeled using Diagram program: https://app.diagrams.net/).
Fig. 5. In silico solvent analysis of the CYP153A6 water tunnel molecular surface using 3D-RISM. Left: hydration free energy ranking of the predicted water molecules. Right: distribution of 3D-RISM-predicted water molecules within the identified water tunnel. The cofactor heme is shown in magentas sticks, water molecules are shown in the gray (H) – red (O) ball-stick representation, the water tunnel is shown as gray molecular surface. (Protein structures and the in-silico solvents were generated using MOE 2022.02^25^.
All three tunnels are structurally framed by the enzyme’s FG and BC loops. Consistently with their role in regulating tunnel access, classical MD simulations revealed that these loops are highly flexible, as indicated by the root mean square fluctuation (RMSF) plot (Fig. 6).
Fig. 6RMSF plot of Cα atoms from the CYP153A6 protein during MD simulations.
CYP153A6 ensemble Docking study
The binding site of CYP153A6 comprises the heme-binding region and an adjacent substrate-binding pocket, a structural hallmark of CYP450 enzymes. This pocket is deeply embedded within the enzyme’s core and features an iron atom coordinated within a porphyrin ring^26^. In its resting ferric (Fe^3+^) oxidation state, the cofactor heme iron is critical for the enzyme’s catalytic function. Several loops surrounding the active site exhibit conformational flexibility, allowing them to adapt to various ligand structures during the binding process^27^. In addition, multiple molecular access/egress tunnels – partially involving the substrate recognition sites (SRSs) – have been proposed for this enzyme^28^. Therefore, we employed an ensemble docking strategy^29^ to further investigate the molecular recognition mechanism of CYP153A6. This docking workflow can be found in Figure S1.
The generation of protein structural ensemble was based on classical MD simulations. Ensemble docking was carried out under four distinct cofactor heme iron charge states and presence or absence of explicit water molecules in the binding site. Water molecules were predicted using 3D-RISM. The four docking scenarios were: (1) ferric iron (Fe^3+^) without explicit water molecules in the binding site; (2) ferric iron (Fe^3+^) with explicit water molecules in the binding site; (3) ferrous iron (Fe^2+^) without explicit water molecules in the binding site; (4) ferrous iron (Fe^2+^) with explicit water molecules in the binding site.
Since the CYP153A6 model was generated via homology modeling, an in silico 3D-RISM solvent analysis was performed to predict solvent distribution. Based on the analysis results, water molecules werethen added to the enzyme binding pocket for the following ensemble docking study. The analysis revealed significant hydration free energy contributions around surface regions within 10 Å of the binding pocket. In contrast, the active site – located between the heme cofactor and the ligand – exhibited markedly lower hydration free energy density, suggesting that the CYP135A6 active site possesses typical hydrophobic properties (Fig. 7).
Fig. 7. In silico solvent analysis of CYP153A6 surface (within a 10 Å cut off from the binding pocket). Hydration free energy density is shown as a solid blue surface. The enzyme cofactor heme is displayed as magentas sticks, and the ligand is shown as green sticks. (Protein structures and the in-silico solvent surface were generated using MOE 2022.02^25^.
A total of 12 already investigated compounds (Table 3)^13^, comprising both experimentally verified substrates and non-substrates of CYP153A6, were used for evaluating the molecular docking workflow and conditions. Each compound was docked in all selected receptor snapshots under all four conditions. The docking results revealed diverse binding poses. In favorable cases, the reactive methyl group of the ligand was oriented toward the cofactor heme iron at a catalytically relevant distance (approx. 6 Å)^30^. In contrast, other poses positioned the methyl group oriented to the cofactor heme iron out of the 6 Å, or in the opposite direction, suggesting non-reactive binding poses.
Table 3. Chemicals for evaluating CYP153A6 Docking conditions.
The docking performance was assessed by calculating the proportion of reactive poses – those with the reactive group oriented toward the heme iron – relative to the total number of poses. To account for conformational variability of the enzyme, molecular docking was performed on three independent MD simulation replicates, and the results are summarized as box plots in Fig. 8. This provided insight for us to select the most suitable docking condition for further CYP153A6 molecular recognition study. Docking poses were ranked using the GlideScore, which is an empirical scoring function that approximates complex-ligand binding free energy. The binding free energy values of the reactive docking poses are around − 6.00 kcal/mol and are more negative than the non-reactive docking poses. The docking score distribution box plots can be found in (Figure S2).
Fig. 8. Box plots showing the distribution of reactive docking poses for published CYP153A6 substrates and non-substrates. Docking conditions: (1) Fe^3+^ without water (blue boxes); (2) Fe^3+^ with water (green boxes); (3) Fe^2+^ without water (red boxes); (4) Fe^2+^ with water (yellow boxes).
Among the four tested docking conditions studied, the setup using the cofactor heme iron in the ferric state (Fe³⁺) without explicit water molecules in the binding site exhibited the best overall performance. This condition most effectively discriminated between CYP153A6 substrates and non-substrates in our test set, as evidenced by narrower box plots, smaller error bars, and fewer outliers for both the ratio of reactive docking poses and the corresponding docking scores. Consequently, this docking setup was selected for further molecular recognition study of CYP153A6.
Toluene derivatives: molecular characteristics of the CYP153A6 binding site and Docking pose analysis
To further investigate the molecular recognition profile of CYP153A6, 14 toluene derivatives (Table S1) were studied. Our validated molecular docking protocol – ferric heme iron (Fe^3+^) without explicit water molecules in the binding site – was utilized to predict the binding modes of these toluene derivatives with CYP153A6.
In substrate-bound models of CYP153A6 with toluene derivatives, some residues were found within 4 Å of the ligand (Fig. 9). These include: ILE88, THR89, MET101, ILE103, THR 186, VAL251, LEU252, VAL255, GLY256, THR260, LEU303, MET306, and PHE407. Among these, hydrophobic residues predominate, such as ILE (2x), MET (2x), VAL (2x), LEU (2x), and PHE (1x), while polar residues such as THR (2x) are less frequent. This composition suggests that hydrophobic interactions are a major driving force in ligand binding within the CYP153A6 active site.
Fig. 9CYP153A6 toluene derivative substrates binding mode. (All the docking poses were produced with the Schrödinger Maestro Suite (release 2024-2)^19^
The binding modes of these toluene derivative substrates, obtained from the top-scoring docking conformations, are illustrated in Fig. 10. To further characterize the interactions between ligands and active site residues, we employed the Schrödinger Protein-Ligand Interaction Fingerprint tool, which automatically generates visual summaries of significant molecular contacts. The resulting fingerprints revealed a variety of backbone and sidechain interactions, highlighting key residues involved in toluene derivative ligand recognition. Specifically, backbone interactions were primarily mediated by ALA97, ILE101, ILE103, LEU252, VAL255, GLY256, LEU303, MET306, and PHE407. While sidechain and hydrophobic interactions involved ILE88, ALA97, MET101, ILE103, THR186, VAL251, LEU252, VAL255, LEU303, MET306, and PHE407. Notably, LEU252 and VAL255 engaged in both backbone and sidechain interactions with almost all the studied ligands, underscoring their dual role in substrate stabilization. The residues involved in the CYP153A6 protein-ligand interaction are shown in Fig. 10, and the detailed analysis of the CYP153A6 protein-ligand interaction (within 3 Å) fingerprints are shown in Table S2. Table S2 shows the toluene derivative structures, and the ligand binding information, as well as the docking scores.
Fig. 10. Analysis of protein-ligand interaction fingerprints of CYP153A6 with toluene derivatives.
Molecular dynamics simulations to verify the molecular recognition of CYP153A6 for toluene derivatives
To further evaluate the binding behavior of toluene derivatives to CYP153A6, we have conducted three 200 ns replica classical MD simulations for the top scoring ligand binding poses derived from the molecular docking. Each simulation generated 1001 frames, and ligand binding free energy was calculated every 10th frame using MM-GBSA methods. The ligand binding efficacy was primarily assessed based on three criteria: (1) the distance of the ligand’s reactive group to the enzyme cofactor heme iron (within a 6.00 Å cut-off), (2) the number of frames meeting this proximity threshold, (3) the reweighted MM-GBSA (Re- ΔG_GBSA_) calculated using only the frames within the reactive range, (4) ligand RMSD. The detailed results are illustrated in Table 4. Additional factors such as the interactions from lipophilic, electrostatic, and ligand packing, were also considered to assess the ligand binding stability and quality, these data can be found in the Table S3. The protein-ligand interaction of CYP153A6 with the toluene derivative during MD simulations is illustrated in Figure S3.
Table 4. Ligand RMSDs, ligand reactive Methyl group to the enzyme cofactor Heme iron distance and binding free energies during MD simulations for the complexes of CYP153A6 and toluene derivatives (values are presented as mean ± standard deviation).Simulation 1Simulation 2Simulation 3Simulation 1Simulation 2Simulation 3Comparatively stronger binders p -cymene
p -xylene RMSD (Å)2.62 ± 0.392.07 ± 0.482.29 ± 0.562.88 ± 1.572.67 ± 0.993.15 ± 1.28Distance (Å)3.85 ± 0.464.58 ± 0.624.11 ± 0.624.30 ± 1.354.82 ± 1.404.11 ± 0.69ΔG_GBSA_ (kcal/mol)−14.14 ± 4.46−27.97 ± 6.54−14.83 ± 5.18−13.03 ± 4.80−15.40 ± 6.14−12.78 ± 5.18Frames1000993994908884983Re-ΔG_GBSA_ (kcal/mol)−14.14 ± 4.46−27.89 ± 6.53−14.83 ± 5.18−12.35 ± 3.77−14.50 ± 5.41−12.72 ± 5.18 Moderate binders
toluene
m -xylene RMSD (Å)3.16 ± 1.094.11 ± 1.243.20 ± 1.073.48 ± 1.502.98 ± 1.184.04 ± 2.20Distance (Å)5.56 ± 2.207.43 ± 2.465.33 ± 1.925.56 ± 1.724.27 ± 0.625.88 ± 2.39ΔG_GBSA_ (kcal/mol)−7.47 ± 3.59−12.79 ± 6.67−9.21 ± 4.05−17.45 ± 5.08−12.47 ± 4.13−17.66 ± 7.28Frames582360645699989582Re-ΔG_GBSA_ (kcal/mol)−6.76 ± 2.51−7.86 ± 5.36−8.63 ± 4.10−15.87 ± 4.67−12.51 ± 4.19−8.41 ± 4.42 p -chlorotoluene
p -methylanisole RMSD (Å)5.04 ± 1.434.20 ± 1.202.84 ± 0.964.25 ± 1.803.30 ± 0.764.12 ± 1.69Distance (Å)9.31 ± 3.088.56 ± 3.204.81 ± 1.817.26 ± 3.284.12 ± 0.596.65 ± 2.87ΔG_GBSA_ (kcal/mol)−22.32 ± 6.00−17.13 ± 8.41−14.22 ± 4.28−13.80 ± 5.24−15.36 ± 4.05−21.85 ± 6.12Frames231284831475991561Re-ΔG_GBSA_ (kcal/mol)−19.62 ± 5.16−19.21 ± 6.74−14.98 ± 3.42−11.82 ± 5.64−15.36 ± 4.05−18.69 ± 5.41 p -nitrotoluene
p -methyl-benzylalcohol RMSD (Å)5.01 ± 2.262.24 ± 0.412.88 ± 1.193.57 ± 1.765.70 ± 1.065.76 ± 0.76Distance (Å)8.97 ± 4.263.93 ± 0.477.15 ± 1.986.25 ± 2.7810.25 ± 2.039.62 ± 0.99ΔG_GBSA_ (kcal/mol)−22.64 ± 7.58−14.73 ± 3.26−23.17 ± 7.69−19.05 ± 9.37−28.05 ± 6.28−36.95 ± 3.74Frames3919972715904920Re-ΔG_GBSA_ (kcal/mol)−15.87 ± 5.39−14.73 ± 3.29−20.63 ± 5.36−12.77 ± 5.45−20.84 ± 7.56−24.98 ± 9.54 Weak binders
o -methylphenol
o -tolymethanol RMSD (Å)2.98 ± 0.564.34 ± 1.863.98 ± 1.431.52 ± 0.321.74 ± 0.651.89 ± 0.39Distance (Å)4.44 ± 0.927.86 ± 3.287.04 ± 1.984.82 ± 0.265.00 ± 0.705.09 ± 0.50ΔG_GBSA_ (kcal/mol)−8.77 ± 4.48−19.82 ± 5.86−14.23 ± 5.30−3.47 ± 6.90−8.30 ± 8.02−7.50 ± 6.24Frames945451441999935929Re-ΔG_GBSA_ (kcal/mol)−8.61 ± 4.50−16.96 ± 3.00−11.71 ± 5.11−3.47 ± 6.90−7.78 ± 7.97−7.23 ± 6.18 o -xylene
m-methyl-benzylalcohol RMSD (Å)4.14 ± 1.344.86 ± 2.405.57 ± 1.286.91 ± 1.924.01 ± 1.474.34 ± 1.11Distance (Å)6.62 ± 2.187.30 ± 2.957.99 ± 1.9211.02 ± 3.137.30 ± 2.007.46 ± 1.65ΔG_GBSA_ (kcal/mol)−16.70 ± 5.71−18.48 ± 6.05−21.95 ± 6.61−28.15 ± 9.15−18.81 ± 10.54−5.57 ± 8.05Frames473422160149316183Re-ΔG_GBSA_ (kcal/mol)−13.79 ± 4.45−14.23 ± 5.00−12.16 ± 6.40−16.13 ± 7.10−20.78 ± 5.09−13.92 ± 7.07 p -cresol
m -cresol RMSD (Å)4.13 ± 1.814.92 ± 0.883.96 ± 1.222.10 ± 0.545.00 ± 0.664.51 ± 1.30Distance (Å)7.30 ± 3.176.06 ± 2.184.73 ± 2.219.98 ± 0.719.29 ± 4.268.11 ± 2.39ΔG_GBSA_ (kcal/mol)−20.55 ± 11.92−13.85 ± 6.47−10.58 ± 6.62−26.24 ± 3.59−19.10 ± 9.69−19.76 ± 8.84Frames3985878930371207Re-ΔG_GBSA_ (kcal/mol)−11.07 ± 8.17−11.76 ± 7.34−8.909 ± 4.390−10.02 ± 5.06−11.65 ± 6.94
The average binding free energy (ΔG_GBSA_) values across simulations revealed significant variation among the ligands, indicating distinct binding profiles. Among the tested toluene derivatives, p-cymene and p-xylene showed consistently strong and stable binding. p-cymene maintained short reactive methyl group distances to heme iron (3.85–4.58 Å) with nearly all frames in this reactive range, and strong Re-ΔG_GBSA_ values (up to −27.89 kcal/mol), indicating a relatively favorable binding mode compared with other ligands. Similarly, p-xylene showed general short reactive distances (4.11–4.82 Å), with over 884 frames. The Re-ΔG_GBSA_ values ranging from − 12.35 to −14.50 kcal/mol.
Moderate binders included toluene, m-xylene, p-methylanisole, p-chlorotoluene, and p-nitrotoluene. Two out of three toluene simulations showed the reactive distance around 5.33–5.56 Å, with more than 582 frames within a 6.00 Å reactive distance. However, in another simulation, the reactive distance was 7.43 ± 2.46 Å, and all Re-ΔG_GBSA_ values ranged from − 13.22 to −2.50 kcal/mol, reflecting relatively weaker binding. m-xylene exhibited consistent reactive contact with the heme iron and moderate binding energy, but with slightly elevated RMSD. p-methylanisole showed a wide variation in reactive distances, ranging from (4.12–7.26 Å), with one simulation with 991 frames within the 6.00 Å reactive distance, while others having 475 and 561 qualifying frames, respectively. In the simulation with good proximity, the binding energy was strong (Re-ΔG_GBSA_ = −15.36 ± 4.05 kcal/mol). p-cholorotoluene displayed mixed results: in one simulation, the reactive distance was favorable with strong energetic binding (−14.98 ± 3.42 kcal/mol), but in others, the average reactive distances exceeded 6.00 Å, significantly reducing their hydroxylation potential. Despite p-nitrotoluene exhibited similar binding energies and reactive distances to p-cholorotoluene, however it is not catalyzed by CYP153A6 in the wet lab experiments. This discrepancy suggests that while binding is similar, the nitro group exerts an unfavorable electronic effect, likely reducing the reactivity of the benzylic C-H bond by withdrawing electron density. Nitro substituents are well-known to decrease the susceptibility of para-position to oxidative reaction due to their strong electron-withdrawing properties^31^.
In contrast, m-cresol, o-xylene, m-methylbenzyl alcohol, o-methylphenol, and o-tolymethanol showed limited binding within the reactive distance. Despite occasionally favorable ΔG_GBSA_ values, they failed to consistently maintain the 6.00 Å reactive distance, with frames often below 391 or reactive distances greater than 8.00 Å, which is too far for effective hydroxylation. m-cresol showed zero frame within range in one simulation and very high average distances (> 9.00 Å) in others. o-xylene and p-methylbenzyl alcohol exhibited poor proximity (reactive distances were consistently greater than 6 Å across all MD simulations) and showed inconsistent energy profiles (the ΔG_GBSA_ values varied substantially among simulations). Similarly, o-methylphenol exhibited variable binding behavior, with one simulation maintaining favorable reactive distance (4.44 ± 0.92 Å and 945 frames within cutoff), while the others exceeded the reactive range (> 6.50 Å) with fewer than 500 qualifying frames. Its binding energies varied moderately (Re-ΔG_GBSA_ from − 8.61 to −16.96 kcal/mol), suggesting sporadic but insufficiently stable interactions for effective hydroxylation. Although o-tolymethanol exhibited relatively closed reactive distances (4.82–5.09 Å) and low RMSD values (< 2.00 Å) across simulations, its Re-ΔG_GBSA_ varied from − 3.47 to −7.78 kcal/mol, indicating an overall unfavorable protein-ligand interaction.
p-cresol had reaction distance below the cutoff only in the replicate characterized by the least favorable Re-ΔG_GBSA_, while other replicates showed larger distances and still not much favorable Re-ΔG_GBSA_. p-methylbenzyl alcohol was a largely ineffective binder, particularly unstable with minimal frames with the reactive distance (20 and 49 frames) in two of the three replicas. Although it presented strong ΔG_GBSA_ values in some frames, the lack of reactive proximity undermines it likelihood of effective hydroxylation.
Overall, for effective protein-ligand interaction with CYP153A6 enzyme, the reactive group of ligands must remain close and stable, positioning to the heme iron over time. Given their strong binding free energies, p-cymene and p-xylene behave as promising candidates for future protein design strategies based on these criteria, while ligands such as m-cresol, o-xylene, m-methylbenzyl alcohol, o-methylphenol, p-cresol and p-methylbenzyl alcohol are less likely to be effective CYP153A6 substrates, due to poor spatial engagement and/or unstable binding.
QM/MM calculations to explore the effect of p-nitro substitution on CYP153A6 catalytical activity
To further investigate how properties of the p-nitro substituent affect CYP153A6 activity, we performed QM/MM calculations. Substrate oxidation in CYP450 enzymes occurs when the enzyme is in the Compound I (Cpd I) intermediate state. The reaction of the iron(IV)-oxo species with the molecule requires an optimal geometry to achieve efficient orbital overlap and effective biotransformation. The reactivity of Cpd I primarily determines theselectivity of CYP450s oxidative reactions^32^.
In this study, we employed QM/MM calculations to evaluate the effect of a p-nitro group on CYP153A6 reactivity by comparing p-nitrotoluene and p-chlorotoluene. As shown in Table 4, MD simulations suggest that both chemicals display similar binding behavior, with p-nitrotoluene even showing slightly better results. However, experimentally, CYP153A6 exhibits high catalytic reactivity toward p-chlorotoluene but almost no reactivity toward p-nitrotoluene.
Referring to previous studies^32^, the quartet and doublet spin states of Cpd I exhibit comparable energetics during hydroxylation and epoxidation steps. Therefore, we only conducted detailed QM/MM calculations on quartet Cpd I here. For both chemicals, the Cpd I complexes were optimized at the QM/MM level. After optimization of p-nitrotoluene-Cpd I complex, the Fe-O bond length is 1.63 Å, the shortest O-H_benz distance is 2.95 Å, the corresponding Fe-O-C angle is 140,9° and O-H-C angle is 94.5°. In the p-chlorotoluene-Cpd I complex has an Fe-O bond length is 1.62 Å, the shortest O-H_benz distance is 2,97 Å, the corresponding Fe-O-C angle is 137,1°, and O-H-C angle is 112,5°. According to Don and colleagues^33^ in another QM/MM calculation for CYP2D6 hydroxylation reaction, the O-H-C angle for reactants with the Cpd I spans from approx. 106° to 144°, and this suggests p-chlorotoluene is closer to a higher-quality near-reactive conformation. (Fig. 11)
Fig. 11. View on the Topology of CYP153A6 Compound I with p-nitrotoluene (left) and p-chlorotoluene (right). Distances are in Å and angles in degrees. (Figures (Figures were produced with the Schrödinger Maestro Suite (release 2025-1)^34^.
Mulliken population analysis was performed during the QM/MM calculations to evaluate the hypothesis that CYP153A6 exhibits poor reactivity toward p-nitrotoluene. The atomic spin density of the heme iron in Cpd I is comparable in the Cpd I – p-nitrotoluene (1.10451) and Cpd I – p-chlorotoluene (1.12181) complexes. Similarly, the spin density of the Cpd I oxo oxygen shows no significant difference between the two systems (nitro, 0.83696 vs. chloro, 0.81862). In contrast, pronounced differences are observed for the substrates. The total spin density of p-nitrotoluene is very small (0.00226), whereas p-chlorotoluene exhibits a substantially higher value (0.35628). This trend is consistent across individual moieties: the aromatic ring (nitro, 0.00063 vs. chloro, 0.26337), methyl group (nitro, 0.00161 vs. chloro, 0.01611), and substituent group (nitro, 0.00002; chloro, 0.07680). The near-zero spin density of p-nitrotoluene indicates a minimal delocalization of the spin density of the porfirin ring to the ligand itself, indicating a much smaller propensity of the site to oxidize the ligand. Differently, p-chlorotoluene exhibits higher spin density, indicating the incoming oxidation of the ligand by the CYP153A6 enzyme.
Atomic numbering for both ligands, cofactor heme iron, and Cpd I oxo oxygen is shown in Figure S4, and the corresponding spin density values are summarized in Table S4 and Table S5. QM/MM calculation convergence plots for geometry optimizations are provided in Figure S5.
Materials and methods
Expression constructs
The synthetic gene sequence encoding for CYP153A6 and its two redox partners, ferredoxin and ferredoxin reductase, has been deposited in NCBI database with accession number OM622424. The pET100-CYP53A6 plasmid was expressed in Escherichia coli BL21 (DE3) using heat shock transformation and the recombinant cells were plated on LB (Luria–Bertani: 10 g L^−1^ tryptone, 5 g L^−1^ yeast extract, 10 g L^−1^ NaCl) agar with 100 µg mL^−1^ ampicillin.
Enzyme preparation
CYP153A6 was produced as recombinant protein as described before. A single colony was selected and inoculated in 20 mL of LB medium supplemented with 100 µg mL^− 1^ ampicillin and incubated at 37 °C at 120 rpm in 100 ml Erlenmeyer flasks for 16 h. The seed cultures were inoculated at an initial cell density of 0.1 OD_600_ in 1 L baffled flasks containing 200 mL of SB medium (Super Broth: 32 g L^− 1^ tryptone, 20 g L^− 1^ yeast extract and 5 g L^− 1^ NaCl) with 100 µg mL^− 1^ ampicillin. The cultures were incubated at 37 °C, 150 rpm until an OD_600_ of 0.6–0.8 was reached, then brought to 4 °C for 15 min and induced with isopropyl-β-d-thiogalactopyranoside (IPTG) 0.5 mM for 4 h at 28 °C and 150 rpm. Cell pellets were recovered by centrifugation at 5000 rpm for 20 min at 4 °C, washed once and resuspended in 100 mM potassium phosphate buffer pH 7. Cell free extracts (CFE) were prepared by sonication for 6 cycles 1 min each at 4 °C. The insoluble fraction was removed by centrifugation at 15,000 rpm for 1 h at 4 °C. The concentration of CYP153A6 was measured by CO-binding affinity using an extinction coefficient of 91.9 mM^− 1^ cm^− 1^ at 450 nm (Choi et al. 2012). Cell growth was measured both as optical density (OD_600nm_, Eppendorf BioSpectrometer) and as cell dry weight (mg mL^− 1^) of washed cells coming from a known volume of culture.
Biotransformations, analyses and kinetics
All chemicals were from Sigma-Aldrich and were used without further purification. Biotransformations with WC were carried out in shaken flasks (100 mL) with 50 mg_dry cells_ mL^− 1^ in 10 mL phosphate buffer (100 mM) containing glucose (50 mM) at pH 7.0. Biotransformations with CFE (CYP153A6 concentration in the cell extract was ~ 1.8 µM) were carried out in 10 ml flasks with tight sealing caps, in a final volume of 1 ml of 100 mM potassium phosphate buffer pH 7, using 10 mg of crude extract and 6 mM NADH. The substrates were prepared as 600 mM stock solutions in ethanol and added at a final concentration of 6 mM. After incubation at 28 °C and 160 rpm, the reactions were stopped by extraction with an equal volume of EtOAc two times; the organic solvent was dried over anhydrous Na_2_SO_4_ and directly injected for GC-FID analysis. GC analyses were performed using a Shimatzu GC-2030 gas chromatographer equipped with a flame ionization detector (250 °C). Chromatographic conditions were as follows: column: Polyethylene Glycol (PEG) Mega-Wax^®^ (30 m × 0.25 mm); injection volume: 1 µL (split (1/20, 250 °C); injection solvent: EtOAc; carrier: H_2_ (0.6 mL/min). Analyses were performed with the following program: (i) gradient from 50 °C to 230 °C (20 °C/min), (ii) isocratic at 230 °C for 3 min. Data were processed with Shimadzu’s LabSolutions software. Product concentrations were determined using conversion factors derived from standard curves of authentic samples. The initial rates used for kinetic analysis were measured at eight substrate concentrations ranging from 0.1 to 10 mM. Reaction rates were plotted against substrate concentration, and the kinetic parameters V_max_ and K_m_ were obtained by non-linear regression of the Michaelis–Menten equation. Errors represent the 95% confidence intervals derived from the variance of the regression fit.
Homology modeling of CYP153A6 protein structure
The 3D structure of CYP153A6 was constructed using the Homology Modeling tool within the Multiple Sequence Viewer/Editor in the Schrödinger Maestro Suite (release 2024-2)^18^. The template was selected based on the highest sequence identity to CYP153A6, the greatest percentage of BLAST positives, and an e-value close to zero. To guide accurate modeling of the holo form, the ligand bound in the template’s active site was retained. After model building, the structure was refined using the Protein Preparation Wizard. Especially, the bond orders were assigned according to the Cambridge Crystallographic Database (CCD) rules for both the protein and the ligand. For the heme cofactor, the four bonds to the central Fe atom were set to zero-order bonds. In addition, a formal charge of − 1 was assigned to the two nitrogen atoms in the porphyrin ring that are bonded through only two single bonds. Model quality was evaluated using a Ramachandranplot. Then hydrogen atoms were added according to the primary sequence and refined based on the predicted ionization states at the experimental pH (7.4). The final homology model was used in the subsequent molecular docking and MD simulations to study the enzyme-substrate recognition.
Molecular tunnel analysis
CYP153A6 molecular tunnels were identified using CAVER 3.02 software^24^, analyzed with CAVER Analyst 2.0^35^, and visualized using PyMOL^23^. The heme iron atom was set as the starting point for tunnel calculations. Detection parameters included a probe radius of 0.9 Å, shell radius of 3.0 Å, shell depth of 4.0 Å, and a clustering threshold of 3.5 Å. Tunnel analysis was performed on frames from a MD simulation on the CAVER software, which considered all protein residues, the heme cofactor, and the bound ligand. Simulations were run with YASARA parameters^36^ (density: 0.997 g/cm^3^, pH: 7.4, temperature: 298.15 K, duration: 10,000 ps, frame interval: 10 ps), generating 1000 frames. The top five tunnel clusters were extracted and analyzed for bottleneck radii and the throughput score.
In Silico solvent analysis
To validate the predicted water tunnel in CYP153A6 and generate water molecules for subsequent ensemble docking studies, an in silico solvent analysis was performed on the CYP153A6 homology model. The model was first processed using the QuickPrep module in Molecular Operating Environment (MOE 2022.02)^25^ for protein structure preparation, energy minimization, and protonation. Subsequently, the ‘Solvent Analysis’ tool in MOE was used to predict water molecule positions within the enzyme structure. Both protein preparation and solvent analysis were conducted using the AMBER10:EHT force field. A cutoff distance of 10 Å from the retained ligand in the homology model was applied to verify the predicted water tunnel. While a separate 10 Å cutoff from the binding pocket was used to investigate how water molecules influence the enzyme’s binding site and to generate water molecules for subsequent ensemble docking studies.
Ensemble Docking study for selecting the best Docking method
An ensemble docking strategy (Figure S1) was utilized to investigate CYP153A6 enzyme ligand recognition under four conditions:
- Cofactor heme iron as Fe^3+^ (ferric), without water molecules in the binding site;
- Cofactor heme iron as Fe^3+^ (ferric), with water molecules in the binding site;
- Cofactor heme iron as Fe^2+^ (ferrous), without water molecules in the binding site;
- Cofactor heme iron as Fe^2+^ (ferrous), with water molecules in the binding site.
Protein structural ensembles were generated from 500 ns classical MD simulations, performed in three replicas using Desmond simulation engine in Maestro with the OPLS4 force field. Systems were built with the System Setup panel using an orthorhombic box (SPC solvent model, 15 Å buffer box size), and neutralized with Na^+^ ions. Trajectory analyses included RMSD and RMSF analyses. Clustering of the MD trajectory (every 10th frame) was performed with Schrödinger’s trajectory_cluster.py (2024-2). The top ten cluster medoids were selected as representative structures for ensemble docking. Protein receptor grid generation was performed using the Grid Generation panel. The Van der Waals radius scaling factor was set to 1.0, with a partial charge cutoff of 0.25. The grid box was centered on the centroid of workspace ligand (the ligand of the holo-CYP153A6 homology model) and defined with a size of 20 Å.
Ligands were prepared with LigPrep panel using OPLS4 force field, generating all relevant ionization states at pH 7.4 ± 2.0 and retaining specified chirality. Docking was performed using Glide in SP (standard precision) mode with flexible ligand sampling. Epik state penalties were included in the scoring, and the post-docking minimization was applied for all the generated docking poses. All the docking poses were performed post-docking minimization. Finally, the docking poses were ranked by docking score.
Molecular Docking and interaction fingerprints analysis of toluene derivatives
The best docking condition (cofactor heme as Fe^3+^ without water molecules in the binding site) identified in Sect. 3.7 was used to study the substrate recognition of toluene derivatives by CYP153A6. Docking was performed on the best target receptor conformation; protein preparation was the same as described in Sect. 3.4. The docking setup followed the same procedures outlined in Sect. 3.7.
Protein-ligand interaction fingerprints were generated using Maestro’s Interaction Fingerprints panel. The “Docked poses” option was selected, and all interaction types were included in the analysis. Residues interacting with the ligand were account to create a heatmap for visualization, the heatmap image was made by Python with the Seaborn and Matplotlib libraries.
Molecular dynamics studies on toluene derivatives
MD simulations of CYP153A6 complexes with various toluene derivatives were conducted for 200 ns in three replicas, following the same simulation protocol described in Sect. 3.1. Ligand RMSDs were analyzed using Schrödinger’s trj_rmsd.py script (version 2024-2), the distances between the ligand reactive methyl group to the enzyme cofactor heme iron during the simulations were monitored using Schrödinger trajectory_asl_monitor.py script (version 2024-2), ligand energy profiles were obtained using the Schrödinger thermal_mmgbsa.py script (version 2024-2) on every 10th simulation frame.
The mathematical formula of MM-GBSA calculating a ligand binding free energy (ΔG_bind_) to a protein are as^37^:
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {G_{bind}}\; = {\text{ }}\Delta H{\text{ }} - {\text{ }}T\Delta S{\text{ }} \approx {\text{ }}\Delta {E_{MM}} + {\text{ }}\Delta {G_{sol}} - {\text{ }}T\Delta S$$\end{document}Where, ∆E_MM_, ∆G_sol_ and -T∆S are the changes in gas-phase molecular mechanics energy, solvation free energy, and conformational entropy upon binding, respectively. The ∆E_MM_ terms consists of:
- ∆E_internal_: bond, angle, and dihedral contributions.
- ∆E_electrostatic_: electrostatic interactions.
- ∆E_vdW_: van der Waals interactions.
The solvation free energy, ∆G_sol_, is the sum of:
- The electrostatic solvation energy (polar contribution), ∆G_GB_.
- The non-electrostatic solvation component (nonpolar contribution), ∆G_SA_.
QM/MM calculations setup
QM/MM calculations were performed using the Schrödinger Suite (2025-1)^34^. The QM region was defined as the ligand, the heme, the oxygen atom and the side chain of Cys367 residue that binds to the heme (with a hydrogen cap added between the Cα and Cβ atoms). The QM region was treated with Jaguar at the DFT level using the B3LYP-D3 functional, while the MM region was modeled with the OPLS4 force field, including the electrostatic field generated by MM atomic charges. For both geometry optimizations and single-point energy calculations, the LACVP+ basis set was applied to the heme iron atom, and the 6-31G** basis set was used for all other atoms.
Conclusions
This study combined experimental biocatalysis and multiscale computational modeling to elucidate the molecular recognition and catalytic behavior of CYP153A6 toward toluene derivatives.
Whole-cell and cell-free assays confirmed that CYP153A6 selectively hydroxylates the benzylic methyl group of several substrates, with p-methylanisole displaying the highest turnover frequency (52.3 min⁻¹). Substrates bearing hydrophobic para-substituents, such as p-cymene, p-xylene, and p-chlorotoluene, were effectively converted, whereas derivatives with strongly polar or electron-withdrawing groups, such as p-nitrotoluene and m-cresol, showed negligible activity.
Computational analyses provided amechanistic explanation for this selectivity. MD simulations and binding free-energy calculations revealed a predominantly hydrophobic active site, stabilized by residues including ALA97, ILE101, ILE103, LEU252, VAL255, LEU303, MET306, and PHE407. Hydrophobic substituents align productively near the heme center, while polar groups disrupt favorable orientations. QM/MM calculations supported the different reactivity observed between p-chlorotoluene and p-nitrotoluene.
Together, these results highlight that hydrophobic complementarity and electronic compatibility with Compound I are the primary determinants of CYP153A6 substrate specificity. This mechanistic insight establishes a rational basis for protein engineering to broaden substrate scope and enhance selective benzylic hydroxylation.
Supplementary Information
Below is the link to the electronic supplementary material.
Supplementary Material 1
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Simonian, N. G., Brutiu, B. R., Kaiser, D., Maulide, N. & Cationic Iodine (III)-Mediated and directed diastereoselective oxidation of inert C–H bonds in Cyclic hydrocarbons. Angew. Chem. Int. Ed.64, e 202421872 (2025).10.1002/anie.202421872 PMC 1193353339895565 · doi ↗ · pubmed ↗
- 2Chovancova, E. et al. CAVER 3.0: a tool for the analysis of transport pathways in dynamic protein structures. Plos Comput Biol.8 (10), e 1002708 (2012)10.1371/journal.pcbi.1002708 PMC 347566923093919 · doi ↗ · pubmed ↗
- 3Land, H. & Humble, M. S. YASARA: A tool to obtain structural guidance in biocatalytic investigations. Methods Mol Biol.1685, 43–67 (2017).10.1007/978-1-4939-7366-8_429086303 · doi ↗ · pubmed ↗