Visualization and figure preparation (I–V)

Visual analyses were performed in the Bodil modelling environment (Lehtonen et al. 2004). Molecular figures were prepared using MOLSCRIPT (Kraulis 1991), Raster3D (Merritt and Murphy 1994) and VMD (Humphrey et al. 1996). Rocker (Lätti et al. 2016) was used for ROC curves.

5.1 Finding substrate binding mode and site of metabolism (I, II) The prediction of small-molecule binding mode and SOM is a widely considered part of CYP metabolism prediction. Enzyme structure-based methods provide insight into the atomic interactions of ligand-protein complexes. Here, MD simulations and binding free energy calculations using MMGBSA were utilized for binding mode and SOM prediction of seven cases of ligand-CYP enzyme complexes (Table 4). First, on the simpler approach, only the MMGBSA binding energy was used to select the most favoured binding pose. Atomic distances of the ligands to the heme iron were used to predict the SOM. As shortly discussed, the MD/MMGBSA method was also tested on two other protein target cases to evaluate the method on the recognition of crystal structure-like ligand poses and on the correlation of MMGBSA with experimentally determined IC50 values. On the second approach, the binding energy-based binding mode and SOM prediction of CYP substrates were augmented with enhanced stability and atomic distance analysis, and the method was experimented with two CYP substrates.

Four molecular docking poses of each small molecule were subjected to short 2.4 ns MD simulations. Twenty snapshots of each simulation were analyzed as is (II) and/or further minimized (I, II). For each pose, the binding energy was calculated. For the CYP ligands, the atomic distances of the ligands to the catalyzing heme iron were calculated. The stabilities of the poses were analyzed visually (I) or, in addition, using ligand and heme-based RMSDLH and the standard deviation (STD) of the atomic ligand-iron distances (II).

Binding mode prediction of CYP ligands based on MMGBSA values was successful in four out of seven cases (Table 4). The second approach, with the incorporation of the enhanced stability and accessibility analysis, could predict the binding mode and SOM correctly for both applied cases, one of which failed in the binding energy-based prediction (Table 4). Thus, the dynamic stability and accessibility of potential SOMs were crucial to consider in binding mode and SOM prediction in addition to the binding energy. Although the simulations

36

were relatively short, it was found that a single binding mode can be obtained from multiple different docking poses in the MD simulations. Consequently, the correct predictions required also manual analysis. A post-minimization of the simulation snapshots did not enhance the prediction.

TABLE 4 Binding mode and SOM prediction based on MD metrics.

Ligand CYP

isoform Experimentally determined

SOMs or metabolites Predicted SOMs

Top 1 (Top 2) Utilized

1*: Binding mode prediction based on MMGBSA binding energy, SOM prediction based on atomic distances of substrate positions to heme iron. In addition, binding mode stability evaluated visually.

2*: Binding mode and SOM prediction based on MMGBSA binding energy, RMSDLH, atomic distances of substrate positions to heme iron and the atomic distance stability (STD).

5.1.1 Binding mode prediction based on MMGBSA (I, II)

The binding modes of 6MC, 7MC and 7FC in complex with CYP2A6 and 2A5 (I), and coumarin in complex with CYP2A6 (II), were attempted to be predicted based on calculated MMGBSA binding energy from 2.4 ns MD simulations of four initial docking poses of each ligand. The computational prediction was compared to experimental results.

The metabolites of 6MC were determined experimentally only for 7-hydroxylation, and the 7MC metabolites are tentative as no standards were available for the acquired mass spectrometry data. 6MC is converted to fluorescent 7-hydroxy-6-methylcoumarin by both CYP2A6 and CYP2A5 with similar rate to coumarin 7-hydroxylation by CYP2A6. 7MC was found to be a mechanism-based

inactivator for CYP2A6 but not for CYP2A5. CYP2A6 converts 7MC to the tautomer of the molecule and hydroxylates the molecule at positions 3 and/or 4.

CYP2A5 hydroxylates the molecule at 7-methyl and at position(s) 3 and/or 4. 7FC was determined to be a mechanism-based inactivator for both CYP2A6 and 2A5, but possible metabolites were not determined experimentally. (I)

The MMGBSA binding energy-based prediction of a single binding mode of 6MC and 7MC in complex with CYP2A6 and CYP2A5 was in accordance with experimental results (Table 4). Other SOM(s) than the experimentally confirmed ones were suggested by the MD simulations of 6MC and 7MC in complex with CYP2A6. The MMGBSA binding energy of the energetically preferred ligand pose was often within error limits of at least one other pose. In addition, it was common that a ligand could acquire an identical binding pose in the simulations from multiple differing docking poses. The binding mode of 7FC in CYP2A6 and 2A5 could not be predicted as the binding energy values of the poses were very close and the SOM(s) or metabolite(s) of 7FC were not determined experimentally.

5.1.1.1 Successful primary predictions

Selection of a single 6MC binding pose based on MMGBSA binding energy was in accordance with the experimental result of 6MC 7-hydroxylation by both CYP2A6 and 2A5. In the MD simulations of 6MC in complex with both enzymes, the energetically preferred pose had the position 7 and 6-methyl closest to the heme iron (I Fig. 2 B). The poses were also the most stable ones based on visual analysis. In addition, in the CYP2A5 simulations, the position 7 is closest to the heme iron in two of the most preferred binding poses (I Fig. 2 B). Both poses are stable by visual analysis and the binding energies are within error limits of each other.

The simulations of 6MC in complex with CYP2A6 suggested that other than the experimentally determined 7-hydroxylated metabolite could be possible. The MMGBSA values of two poses are within error limits of the most favored pose.

In one of these two poses, the position 5 stays close to the heme iron and is thus available for oxidation. However, the experimental studies tested 6MC only for 7-hydroxylation and thus a metabolite catalyzed at the position 5 has not been experimentally confirmed.

The top binding pose from the 7MC simulations in complex with CYP2A6 and 2A5 agrees with experimentally determined metabolites. In addition, the comparison of the CYP2A6 and 2A5 simulations indicated a difference in the compound’s binding to these enzymes, which could explain the mechanism-based inactivation of CYP2A6 by 7MC. The pose in which the 7-methyl is the closest position to the heme iron is favored energetically and by stability in both CYP2A6 and 2A5. The same pose was acquired in three out of four docking poses in the CYP2A5 simulations, and in all there is space available for oxygen to bind the heme iron (I Fig. 2 D, F). Accordingly, a metabolite hydroxylated at the 7-methyl was indicated by the experimental results. In contrast to the acquired 7MC-CYP2A5 complexes, the 7-methyl group lies at a very close proximity to the CYP2A6 heme iron and leaves no space for oxygen to bind the iron (I Fig. 2 C, F).

38

This could explain the mechanism-based inactivation of CYP2A6 by 7MC. The reaction mechanism is outside the scope of MD simulations, but a hypothesis was be derived based on the coordinates of the favored binding mode. The low distance between the methyl and the heme iron could lead to formation of a 7-methylene-7H- chromen-2-ol tautomer, which the experimental data suggests for CYP2A6 but not for CYP2A5. The formed tautomer could bind tightly with the iron by the tautomer’s pi-electrons or react irreversibly with the binding site environment in an unidentified way.

The other simulated poses of 7MC in complex with CYP2A6 and 2A5 suggested also other metabolites for the compound. In CYP2A6, other SOMs for 7MC could be the positions 4 and 5, but the experimental results suggested a 3- and/or 4-hydroxylated product. In CYP2A5, the single simulation where the 7-methyl was not closest to the heme, positions 3 and 4 lie close to the heme.

Although the pose was energetically clearly worst, it showed great stability and could thus be considered a secondary binding mode of 7MC in CYP2A5. This is in line with the experimental indication of a 3- and/or 4-hydroxylated product of 7MC by CYP2A5.

5.1.1.2 Wrong or ambiguous predictions

The prediction of coumarin binding mode, based on MMGBSA binding energies, did not agree with experimentally determined 7-hydroxylation by CYP2A6. The molecular docking yielded coumarin poses, where the closest sites to the heme iron are (1) 2-carbonyl, (2) position 7, (3) 2-carbonyl and position 3 and (4) positions 5 and 6 (II Fig. 1 B). The calculated MMGBSA values are within error limits of each other in the poses (II Table 1). The pose 4 is the most preferred, followed by pose 2, 3 and 1. Accordingly, the binding energy estimations and the atomic distances suggest that position 5 would be the most preferred SOM, although the correct pose 2 is also within error limits of the top pose. From four differing docking poses, two general coumarin orientations were acquired in the simulations where either the position 7 or 5 is closest to the heme iron (II Table 1). The mechanism-based inactivation of CYP2A6 and 2A5 by 7FC could not be deduced from the MD simulations. Based on the binding energy, the preferred binding mode of 7FC in complex with CYP2A6 and 2A5 is not clear as the MMGBSA values of the poses are close and steadily within error limits of each other. However, the 7-formyl of 7FC orientates to very close proximity of the CYP2A6 heme iron at the end of the simulation of the pose 4 (I Fig. 2 H). This could indicate that the 7-formyl group plays a role in the mechanism-based inactivation, which was backed up by literature of commonly observed mechanism-based inactivation of CYP enzymes by aldehydes.

5.1.2 Enhanced stability and distance analysis (II)

Since the MMGBSA binding energies of different CYP substrate poses were in many of the above cases within error limits of each other, it was reasonable to take a closer look into other variables of the simulated poses. In order to

investigate some of the options to analyze the energetics and stability of substrate poses, the already discussed MD simulations of four initial docking poses of coumarin and 6MC in complex in CYP2A6 were examined. Both coumarin and 6MC are 7-hydroxylated by CYP2A6 with similar rate (I). Other metabolites of 6MC were not determined (I) and the literature does not suggest other metabolites for coumarin.

A set of simulation-derived metrics was suggested for predicting the preferred binding mode and SOM of a CYP substrate. The binding energy comparison of the simulated substrate poses was incorporated with a simple calculation of the binding pose stability and the atomic distances the substrate to the heme iron. Namely, for each simulation, the RMSD of the substrate and the CYP2A6 heme were calculated together (RMSDLH) in order to consider the stability of the substrate pose in relation to the first frame of the simulation. On one hand, the exclusion of the CYP amino acid residues from the RMSD analysis and superimposition removes the effect of the protein macromovement to the substrate superimposition. This allows a more delicate analysis of the substrate stability. On the other hand, when the heme is included, the RMSDLH indicates of the stability and orientation of the simulated pose in relation to the binding site. Atomic distances and stabilities, the latter indicated by the STD of each distance, of the substrate positions to the CYP heme iron were calculated for each simulation. The closer substrate positions to the heme iron are more likely to be oxidized. The closest substrate position should also be well stabilized at the supposed reaction coordinates, which is indicated by the STD of the mean distance. The distance should be low enough for the position to be available for oxidation, but leave space for the reactive oxygen.

In both coumarin and 6MC simulations, the combination of MMGBSA and the above metrics of ligand stability and atomic distances resulted in correct binding mode and SOM prediction. However, especially in the case of coumarin, this required a manual analysis of the ligand trajectories in the individual poses.

Similar or close to identical binding poses were acquired from several simulations although the simulations of both compounds were initiated from four notably differentiating docking poses.

5.1.2.1 Coumarin

Combining the information from MMGBSA binding energy, RMSDLH and the atomic distances together, the coumarin position 7 is preferred at the proximity of the CYP2A6 heme iron, thus in agreement with experimentally determined 7-hydroxylation. The prediction relied partly on manual analysis of the simulations as the orientation of coumarin changed during the individual trajectories. The pose is energetically the second-preferred, but it is the most stable and the position 7 is closer to the heme iron than position 5 in the other binding poses.

The preferred pose is also closest to the crystal structure of coumarin in complex with CYP2A6. The second-preferred binding pose, where the coumarin position 5 is closest to the heme iron, could also be considered as a naturally occurring binding mode as the calculated descriptors are within error limits of the top

40

binding pose. However, coumarin metabolites catalyzed at the position 5 by CYP2A6 are not found in the literature.

5.1.2.2 6-methylcoumarin

The prediction of the 6MC binding mode and SOM agree with the experimentally determined 7-hydroxylation by CYP2A6. Combining the information together, the preferred binding mode of 6MC orients the 6-methyl and position 7 closest to the heme iron. The MMGBSA binding energy is the best for this pose and the positions 6-methyl and 7 are closer to the heme iron and more stable than the position 5 in pose 4, which is second-ranked by MMGBSA and top-ranked by RMSDLH. In the preferred pose 1, the 6-methyl and position 7 are accessible for oxidation. In this pose, difference between the preference of oxidation of the position 6-methyl and 7 cannot be made based on the applied metrics. Another binding mode where the position 5 (pose 4) or 4 (pose 3) is oriented towards the heme is also possible based on the MD simulations. Although the MMGBSA, position 5 distance to the heme iron, and the distance STD are slightly inferior to those of pose 1, the values are very close. However, only data for the 7-hydroxylation of 6MC by CYP2A6 is available and thus the predicted second-preferred binding mode and SOM are not verified experimentally.

5.1.3 MMGBSA in other applications (II)

In VS and molecular docking, it is important that the method recognizes a true binding mode of ligands and that the scores correlate with experimental results.

The ability of molecular docking and MD/MMGBSA to accomplish these objectives was investigated with twelve androgen receptor (AR) crystal structure ligands and 152 phosphodiesterase 4B (PDE4B) inhibitors with experimentally determined IC50 values. With AR, the MD/MMGBSA method did not generally improve the recognition of the crystal structure-like ligand pose compared to docking. In the case of PDE4B, the correlation of experimental IC50 values with MMGBSA binding energies was significantly higher than with docking scores.

5.1.4 The effect of post-minimization (II)

The effect of minimization of the MD simulation production run snapshots to the MD/MMGBSA approach was investigated in the cases of CYP2A6 in complex with coumarin and 6MC, AR and PDE4B. In the case of CYP2A6 substrates coumarin and 6MC, it was found that the post-minimization procedure does not improve the binding mode and SOM prediction in these two cases. The ranking of the simulated binding poses was generally unchanged after the minimization.

The MMGBSA and RMSDLH values of the different poses are within error limits (STD) of each other similarly as they are in the non-minimized MD simulations.

However, post-minimization brings coumarin closer to the crystal structure and the distance of position 7 is closer to the heme iron, a similar distance as for a typical distance seen in CYP2A6 crystal structures. In the case of AR, the post-minimization slightly improved the crystal structure-like ligand pose

recognition. The post-minimization did not have a significant effect on the correlation of MMGBSA and IC50 values of PDE4B inhibitors.

5.2 Profluorescent tool molecules – one binding mode of interest (III, IV)

Profluorescent substrates can be used to determine catalytic activity and inhibition of CYP enzymes. There are many such tool molecules available, but they often have poor selectivity as they are catalyzed by multiple CYP enzymes.

Coumarin-based molecules are abundant for this purpose as nonfluorescent coumarin derivatives can be metabolized to a fluorescent 7-hydroxylated metabolite in a typical CYP reaction (III Fig. 1). The coumarin core offers a wide array of options for chemical functionalization and thus potential for varying the selectivity profile of the 7-hydroxylation reaction by CYP enzymes.

Previously published and new 3-phenylcoumarin derivatives (III Fig. 2, IV Fig. 1) were computationally modelled (Fig. 5) and experimentally tested for 7-hydroxylation by 13 human xenobiotic-metabolizing CYP enzymes. The CYP enzymes are hepatic except for CYP1A1 and 1B1. The CYP1 family members was identified as the most suitable target proteins to accommodate the 3-phenylcoumarin core and subsequently catalyze the 7-hydroxylation reaction.

Experimental results confirmed that 19 out of 21 tested 3-phenylcoumarin derivatives were 7-hydroxylated with measurable activity by either or by both CYP1A1 and 1A2. Other isoforms to frequently catalyse the reaction were found to be CYP1B1, 2A6, 2D6 and 2C19 (III Fig. 6, Table 1, IV Fig. 2). The interactions of the 3-phenylcoumarin derivatives with selected CYP isoforms were studied with docking and 24.0 ns MD simulations in the binding pose that was hypothesized to be the most prominent to facilitate the 7-hydroxylation reaction.

The computational modelling suggested key ligand-enzyme interactions and properties of the CYP enzymes that affect the preference of the 7-hydroxylation reaction. Some conclusions may also be made of the general properties and key differences of the CYP1 family based on the computational modelling.

FIGURE 5 Molecular modelling workflow to study the 7-hydroxylation selectivity of 3-phenylcoumarin derivatives by CYP enzymes.

42

5.2.1 Identification of target isoforms

The comparison of the CYP 3D structures (III Table S1) revealed that the CYP1 family and CYP2A6 would be most suitable to accommodate the 3-phenylcoumarin core to their respective binding sites for 7-hydroxylation. Both the CYP1 family and CYP2A6 have binding sites that are well suited for hydrophobic and planar ligands. However, the CYP2A6 binding site is smaller than those of the CYP1 family, and is likely not able to accommodate a wide array of 3-phenylcoumarin derivatives. The 3-phenylcoumarin core represents alpha-naphthoflavone that is co-crystallized with the CYP1 family enzymes. Both small molecules are planar, have three to four closely attached ring structures and a single carbonyl group. The CYP2A6 is co-crystallized with coumarin which is obviously highly similar but slightly smaller than 3-phenylcoumarin. In conclusion, the binding sites of the CYP1 family enzymes and CYP2A6 have excellent topology to bind 3-phenylcoumarins.

5.2.2 The CYP1A subfamily and CYP2A6

In the CYP1A subfamily and CYP2A6, the binding poses from docking suggested that 3-phenylcoumarin 2-carbonyl could form a hydrogen bond (H-bond) with Ser122, Thr124 and Asn297 of the 1A1, 1A2 and 2A6 isoforms, respectively (Fig.

6 A, III Fig. 4 A, B and D and S, IV). In this pose, the position 7 can be oriented towards the catalyzing heme. The H-bond from the 2-carbonyl could be sufficient to stabilize the 3-phenylcoumarins to this binding mode. Accordingly, other

6 A, III Fig. 4 A, B and D and S, IV). In this pose, the position 7 can be oriented towards the catalyzing heme. The H-bond from the 2-carbonyl could be sufficient to stabilize the 3-phenylcoumarins to this binding mode. Accordingly, other