The Metabolic Rainbow: Deep Learning Phase I Metabolism in Five Colors

Current Approaches to Modeling Metabolism.

The existing in silico approaches to modeling phase I metabolism fall into three categories: (1) generic site of metabolism (SOM) prediction, (2) reaction-type specific SOM prediction, and (3) multiple reaction-type SOM prediction. All existing approaches fail to accurately predict bioactivation due to limitations in their ability to unambiguously identify metabolite structures from model predictions.

In generic SOM prediction, models predict which atoms or bonds are modified by metabolism. Sites are labeled as metabolized or not metabolized, but the specific type of metabolic reaction occurring at the site is not specified. SMARTCYP, RSPredictor, SOMP, MetaPrint2D, FAME 2, and XenoSite 1.0 are examples of models that use the generic SOM approach.^24-30 These models are very effective at identifying metabolically labile sites, and they have become widely used in industry as a guide in designing molecules with increased half-lives. However, many sites have the potential to undergo several different types of reactions, yielding distinct metabolites. Consequently, these models do not provide enough information to unambiguously predict the metabolite structures and cannot model complete bioactivation pathways.

In reaction-type specific SOM prediction, models predict a single type of reaction to eliminate issues with ambiguous labels. For example, our group has published quinone formation, epoxidation, N-dealkylation, reduction, and sulfur-oxidation models.^31-34 This approach has four significant advantages. First, the nuances of each reaction type can be carefully considered in isolation. Second, expert knowledge can guide design of input features, which can substantially improve accuracy. Third, curation of data sets of a uniform reaction-type encourages consistent annotation and identifies when data must be supplemented for reasonable performance. Fourth, the metabolic transformation is unambiguous, and the metabolite structure can be predicted. This approach has been invaluable in studying bioactivation pathways, including modeling the bioactivation of structural alerts.^33,34 However, this approach has three key limitations. First, each model is typically trained on disjoint data sets, so their outputs are not scaled to predict the relative likelihood of different reaction-types. Second, information common to different reaction-types is not exploited to improve accuracy. Third, only a small minority of reaction-types have been modeled this way.

In multiple reaction-type SOM prediction, a single model is used to simultaneously predict multiple phase I reactions types. These models attempt to address the limitations of both generic SOM and reaction-type specific SOM models. To our knowledge, only two reaction-type SOM models—SOME³⁵ and He et al.³⁶—have been reported in the literature. In this work, we show that these published models only cover 42–48% phase I metabolism and leave out several reaction types that frequently generate reactive metabolites like epoxidation³⁷ and quinone formation.³⁸ In addition, we argue that these existing models are not probabilistic and are thus difficult to interpret. Nevertheless, the multiple reaction-type approach has potential if three key issues can be addressed: (1) coverage of phase I reactions is limited, (2) existing models are not probabilistic and do not indicate the relative likelihood of different metabolic transformations, and (3) no effective, scalable method for visualizing and interpreting these models is available. As the number of reaction types modeled proliferates, our ability to visualize and coherently organize results is diminished.

Unambiguous Reaction-Type Specific Site of Metabolism Model.

The central innovation put forward in this study is a “five color”, metabolic-rainbow classification of reaction types (Table 1). In this framework, bonds and lone-pairs (sites) are labeled with one of five colors, each corresponding to a set of possible reactions at that site. Given the context of the labeled sites, the colors can be unambiguously decoded into one of the 21 reaction types. We grouped these diverse phase I reaction types into five classes based the structural changes that they introduce to the substrates: stable oxygenation, unstable oxygenation, dehydrogenation, hydrolysis, and reduction. Two design goals guided the classification of reaction types into five colors. First, knowledge of the site and color of a reaction, should almost always be sufficient to determine the expected metabolite. Second, semantically related reaction types should be grouped into the same color. Detailed definitions of each color (or class) can be found in the Data and Methods section.

Table 1.

Five Colors of Phase I Metabolism

			Number
Color	Definition	Reaction Types	Molecule	Site
Red	Stable Oxygenation (SO)	aromatic hydroxylation aliphatic hydroxylation aromatic epoxidation aromatic epoxidation N-oxidation S-oxidation	3098	10280
Orange	Unstable Oxygenation (UO)	N-dealkylation O-dealkylation S-dealkylation C-dealkylation oxidative deamination oxidative dehalogenation	3749	5811
Green	Dehydrogenation (DH)	dehydrogenation alcohol to aldehyde or ketone single- to double-bond double- to triple-bond quinone formation iminium formation	762	2794
Blue	Hydrolysis (HD)	ester hydrolysis amide hydrolysis ether hydrolysis cyanide hydrolysis	3188	3969
Purple	Reduction (RD)	carbonyl reduction nitro reduction sulfo reduction reductive dehalogenation hydrogenation	1213	1590

Open in a new tab

We demonstrate that the five color approach overcomes several limitations of prior approaches. First, predictions corresponding to 21 phase I reaction types can be unambiguously displayed on just five molecular diagrams, each corresponding to a different color. Second, because the model is trained on all reaction types and colors at the same time, information can be shared across reaction types, which reduces the total number of trainable weights and could improve the performance of the model on rare reaction types. Third, by using the same molecules to train all reaction types, the model outputs are expected to predict the relative probabilities between different reaction-types with greater accuracy.

Using the five color metabolism framework, we constructed a deep neural network model and trained it to produce probabilistic predictions of both the site and the type of human phase I metabolism for small, drug-like molecules. The model computes raw site- and molecule-level prediction scores, each corresponding to the probability that the site or molecule is metabolized by one of the five colors of metabolic reaction. The raw site-level scores are then aggregated and scaled to make general phase I and 21 reaction-type specific predictions. We trained the model on a diverse data set of 20 736 literature-derived human phase I reactions comprising 21 phase I reaction types ranging from the most common types like hydroxylation and dealkylation to less common yet important reactive-metabolite generating reaction types like epoxidation and quinone formation.

Metabolism Training Data.

We collected a chemically diverse set of molecules and metabolic reactions from the literature-derived Accelrys Metabolite Database (AMD). A total of 20 736 human in vitro and in vivo records were gathered. Each record was composed of one or more chemical transformations. A diverse set of phase I reaction types, including hydroxylation, dealkylation, deamination, dehalogenation, and hydrolysis, among others were present in our training data set (Table 1). Our study accounted for 92.3% of all human phase I metabolism in the AMD data set. The remaining 7.7% phase I reaction types excluded from the current study, in order from most to least common, were tautomerization, isomerization, rearrangement, radical formation, hydration, deacylation, denitrogenation, and decarbonylation. Reactions starting from the same substrate were mapped onto a single training molecule using maximum common substructure matching.³⁹ Our final data set contained 9 674 unique molecules. Under the AMD licensing agreement, we were not allowed to disclose the structures of molecules in the data set. However, to enable rebuilding of our complete database and reproduction of our results, we provide all reaction and molecule AMD registry numbers in the Supporting Information.

Five Colors of Metabolism.

Using the five color classification approach, we identified 21 reaction types, and grouped them into five classes: red, orange, green, blue, and purple (Table 1). The first class is red, for stable oxygenation. In a red reaction, an oxygen is added to the substrate, yielding a single observable metabolite. Red reactions include hydroxylation, aromatic/aliphatic epoxidation, nitrogen oxidation, and sulfur oxidation (Table 2). The second class is orange, for unstable oxygenation. In an orange reaction, addition of an oxygen to the substrate destabilizes one of the molecule’s heavy bond (bond connecting between two non-hydrogen atoms). Eventually this bond breaks and gives rise to two metabolites, one contains the added oxygen and the other does not. Orange reactions include nitrogen-, oxygen-, sulfur-, and carbon-dealkylation and oxidative dehalogenation (Table 3). The third class is green, for dehydrogenation. In a green reaction, hydrogen(s) is extracted to form double bond, triple bond, iminium, or quinone. Green reactions include double/triple bond formation and quinone/imine/methide formation (Table 4). The forth class is blue, for hydrolysis. In a blue reaction, addition of a water molecule causes amide-, ester-, or ether-bond breakage and gives rise to two metabolites. Blue reactions include amide, ester, and ether hydrolysis (Table 5). The fifth class is purple, for reduction. In a purple reaction, an oxygen is removed or hydrogen(s) is(are) added to the substrate. Purple reactions include nitro-, carbonyl-, and sulfo-reduction, reductive dehalogenation, and hydrogenation (Table 6).

Table 2.

Red Reactions: Stable Oxygenations (SOs)^a

Reaction Type	Example	SOM Labeling	Site
Hydroxylation		Bond between a heavy atom and a hydrogen	8475
N-oxidation		Lone-pair on nitrogen atom	665
Epoxidation		Double/aromatic bond between two carbon atoms	636
S-oxidation		Lone-pair on surfur atom	424

Open in a new tab

^a

Examples, labeling strategies, and reaction type counts for each “red”, stable-oxygenation reaction type.

Table 3.

Orange Reactions: Unstable Oxygenations (UOs)^a

Reaction Type	Example	SOM Labeling	Site
N-dealkylation		Bond between nitrogen and carbon	2970
O-dealkylation		Bond between oxygen and carbon	1712
S-dealkylation		Bond between sulfur and carbon	93
C-dealkylation		Bond between two carbon	577
Oxidative Dehalogenation		Bond between halogen and carbon	309

Open in a new tab

^a

Examples, labeling strategies, and reaction type counts for each “orange”, unstable-oxygenation reaction type.

Table 4.

Green Reactions: Dehydrogenations (DHs)^a

Reaction Type	Example	SOM Labeling	Site
Double/triple bond formation		Bond between heavy atom and hydrogen	2741
Quinone/Imine/Methide formation		Bond between heavy atom and hydrogen	90

Open in a new tab

^a

Examples, labeling strategies, and reaction type counts for each “green”, dehydrogenation reaction type.

Table 5.

Blue Reactions: Hydrolyses (HDs)^a

Reaction Type	Example	SOM Labeling	Site
Amide hydrolysis		Amide bond	1118
Ester hydrolysis		Ester bond	2112
Ether hydrolysis		Ether bond	47

Open in a new tab

^a

Examples, labeling strategies, and reaction type counts for each “blue”, hydrolysis reaction type.

Table 6.

Purple Reactions: Reductions (RDs)^a

Reaction Type	Example	SOM Labeling	Site
Nitro reduction		Bond between nitrogen and attached oxygen	313
Carbonyl reduction		Carbonyl bond	869
Sulfo reduction		Bond between sulfur and attached oxygen	46
Reductive dehalogenation		Bond between halogen and attached carbon	2
Hydrogenation		Double, triple bond between pair of hydrogenated atoms	1435

Open in a new tab

^a

Examples, labeling strategies, and reaction type counts for each “purple”, reduction reaction type.

While abundant data is available for phase I reactions like hydroxylation and dealkylations, reactions such as quinone formation and cyanide hydrolysis have only a few examples in our data set. Were individual models constructed for each reaction type, some of these models would be susceptible to ascertainment bias due to data limitations.⁴⁰

Colored Sites of Metabolism.

We implemented automated labeling rules that use the structures of substrate and product and the reaction type to identify the specific bond(s) or lone-pair(s), called site(s) of metabolism (SOM(s)), at which the metabolic transformation takes place in each reaction record. Although labeling bonds as SOMs is obvious, labeling lone-pairs as SOMs is unconventional. We define a lone-pair as a pair of unshared valence electrons that do not participate in a covalent bond. For example, there are two lone pairs on the oxygen atom in a water molecule. In phase I metabolic reactions, a lone pair can donate electrons to form a covalent bond. When this occurs, we consider the lone pair is the site of metabolism. Reaction-type specific SOMs are labeled via a set of rules defined for each reaction type (Table 2-6). The final data set had 10 280, 5811, 2794, 3869, and 1590 sites of stable oxygenation, unstable oxygenation, dehydrogenation, hydrolysis, and reduction, respectively (Table 1).

Our labeling rules could identify SOMs for the majority of AMD reaction records, but they were not perfect. Each reaction record was manually inspected to identify mislabeled SOMs. About 15% of all automatically labeled SOMs were identified as incorrect and manually corrected. As with other manual processes, there were human errors in the label inspection. We approximated such errors by randomly selecting a subset of the manually corrected reaction records for further evaluation. For these selected reaction records, multiple independent annotators repeated the manual inspection. Next, all annotators examined and discussed each reaction record together to reach consensus on the correct labeling. Of the 200 sampled reaction records, ten records had annotations that are deviated from the established correct labeling. As a result, approximately 5% of the manually annotated reaction records contained errors.

Identification of Potential Sites of Metabolism.

Each type of metabolic reaction occurs at some but not all sites. For example, epoxidation only occurs at double and aromatic bonds while hydrolysis takes place at amide, ester, and ether bonds. These sets of potential sites could be used to heuristically estimate specific chemical domains for reaction types—where in the molecule a specific reaction happens. To assess our model’s ability to predict sites of metabolism for individual reaction types within their specific chemical domains, we identified all reaction-type specific potential sites (Table S6) using SMARTS patterns.⁴¹ While useful for identifying potential sites, SMARTS patterns alone cannot specifically identify SOM because only a small fraction of potential sites are metabolized. Of 691, 349 bonds and lone-pairs (pseudobonds) in our data set, 461, 876 (66.81%) sites were identified as potential sites for at least one of the covered phase I reaction types. On average, 99.97% of class-specific labeled SOMs were in the corresponding set of potential sites. Class-specific labeled SOMs that did not fall within the corresponding domain for their class were in fact mislabeled sites.

Descriptors.

To predict the susceptibility to metabolism of each bond/lone pair in the data set, our model used structural, physical, and chemical information encoded in numerical descriptors. Each bond was represented by a vector of 404 descriptors that describe its properties at atom, bond, and molecule levels. All of these descriptors are calculated by in-house software from 2D Open Babel SDF files.⁴² The majority of the descriptors used in this project are atom-level descriptors (e.g., atom identity, charge, hybridization) previously developed for XenoSite metabolism models, as described in previous studies from our group.^26,31,43 For each bond, there were 189 descriptors for each atom on either side of the bond (Table S1), 10 bond descriptors (e.g., bond type, bond length) (Table S2), and 16 molecule descriptors (e.g., molecular weight, topological polar surface area, molar refractivity) (Table S3). The complete list of descriptors used in this study is provided in the Supporting Information Tables S1, S2, and S3.

Multitarget Deep Neural Network Model.

Our phase I model was a multitarget, deep neural network with one input layer, four hidden layers, and two output layers (Figure 2). Multitask neural networks can reuse features learned for one class of reactions to predict another, data-limited class, and can significantly reduce the number of model parameters required.^23,44,45

The top output layer computed molecule prediction vectors, each contained five scores corresponding to the five classes of phase I reactions. Each of these molecule prediction scores correlated to the probability that a given molecule was observed to be metabolized via the corresponding reaction class. Similarly, the bottom output layer computed bond prediction vectors, each containing five scores. These scores correlated to the probabilities that a given bond/lone-pair was an actual site of metabolism for each of the five reaction classes.

The network was trained in two stages. First, the bond-level network was trained to produce accurate bond-level scores. For an input molecule, atom-level, bond-level, and molecule-level descriptors were calculated for all bonds and lone-pairs. These descriptors were inputs to a neural network with two hidden layers (each contains 25 nodes, ReLU activation) that generated vectors of five scores corresponding to the five classes of phase I reactions. Each element in the vector represented the probability that a bond/lone-pair was a phase I reaction-class SOM. The weights of the model were calibrated during training by performing gradient descent on the cross-entropy error (with L2 penalty coefficient of 0.3) of the difference between the predicted and actual response values of each bond using TensorFlow and SciPy LBFGS Optimizer.^46-48 Second, the molecule-level network of two hidden layers (each contains five nodes, ReLu activation) was trained. In this training step, each row of the data matrix represented a molecule, and each column was descriptor. Descriptors included the top five bond-level scores for each reaction class (25 in total) and all molecule-level descriptors. Similar to the first stage, we trained the weights of the networks using gradient descent on the cross-entropy error (with L2 penalty coefficient of 0.3) so that molecules observed being metabolized by certain class of reaction received higher score for that class than other molecules.

The accuracy of the model was evaluated using a 10-fold cross-validation procedure. The data set was divided into ten groups with an equal number of substrates in each group. Predictions were generated for each tenth of the data by a model trained using the remaining nine tenths of the data. The resulting predictions covered the complete data set and were used to compute the accuracy metrics in this study.

Visualizing Metabolism in Colors.

There is no established way of displaying SOM data from many reaction types in a coherent or efficient way. One motive in defining reaction classes as “colors” is to introduce a mnemonic metaphor as a foundation for efficient visualizations of metabolism. Reaction color becomes a visual cue for reaction classes and reaction types.

Specifying the five reaction colors more precisely, each was mapped to colorblind safe hex codes.^49,50 Generally, rainbow coloring is often deemed inappropriate for visualization.^51-53 In this case, however, it is appropriate because there are no natural ordering to the reaction classes to which the color is keyed.

In displaying metabolism data and predictions the 2D depiction of the molecule is always used, and it is annotated with both (1) white circles with a black shadow to denote each atom and bond that are metabolized and (2) background shading with color and size that denotes the strength of predictions at each atom and bond. Usually, the 2D depiction of the molecule is black. Here, we use color in the 2D depiction to indicate reaction color. The white circles and rainbow shading on a certain colored 2D depiction are the observed SOMs and predictions of that specific metabolism color.

For the most part, SOMs are visualized exactly as used for modeling described in the prior section, but in some cases, adjustments are made. For a SOM labeled at the bond between a hydrogen and a heavy atom, the circle is displayed at the corresponding heavy atom. Similarly, a lone-pair SOM is visualized at the atom bearing the lone-pair. Likewise, model predictions for a bond between a hydrogen and a heavy atom or for a lone-pair is visualized at the corresponding heavy atom. These modifications allow for displaying all known SOMs and model predictions on a standard 2D molecular diagram, commonly used by organic chemists, which only contains heavy atoms.

Three Levels of Visualization.

We define three levels of visualizations at which to display observed SOMs and model predictions. The three levels of visualization work together (Figure 3). A large number of molecules can be displayed with just the highest visualization level, and lower level visualizations can clarify details when they are requested. Also, multiple levels can be displayed at the same time to display metabolism data at different resolutions (Figure 3).

At the highest level, we show predictions as a single molecule, so this may be most useful for screening large numbers of molecules in the first pass (Figure 3A). Predictions of each of the five reaction types are aggregated into a single score for each SOM, as outlined in the following section (Prediction Aggregation and Scaling). This visualization uses a black molecule depiction with white circles. In the future, we might add color to the circles to indicate reaction class, perhaps including multicolored circles when multiple reaction colors are observed at a single site.

At the middle level, each row has five depictions, one for each reaction color (Figure 3B). This visualization makes immediately apparent the reaction type of most predictions and SOMs. Here, each of the five molecules depicts a reaction color, and observed sites of metabolism are annotated with white circles. The model predictions for each reaction color are displayed as shading behind the molecule of the corresponding color. The color of the depicted molecule is a critical visual cue. With practice, one can distinguish reaction-types based on structural context and the reaction color.

At the lowest level, all reaction types are displayed in separate depictions (Figure 3C). The visualization still colors the molecule depiction with the class color, but there are multiple molecules with the same color. Reported sites of metabolism, again, are marked with white circles. The reaction type name is included too, to clarify which depiction corresponds to which reaction type. This is the most expansive visualization, which is most suitable for drilling down into the details of molecules being examined most rigorously.

Though the predictions in the model can be visualized, it does take effort to learn how to navigate the figure. In the future, we plan to integrate this model as a building block into composite algorithms designed to answer specific questions in drug development. For example, we hope to build algorithms that can predict bioactivation pathways. Ultimately, we expect this visualization may be more useful as a diagnostic for tools that make use of predictions from the model.

Prediction Aggregation and Scaling.

To improve visualization, the model raw predictions ${\mathit{p}}_i^j$ were, at times, aggregated and scaled. The aggregated and scaled predictions are visualized as the background shading behind each molecule using a rainbow- and size-based color map (Figure 3).

The highest visualization described in the prior section contains aggregated predictions (Figure 3A). In the aggregation step, each site i is assigned with a single score p_i by computing the probabilistic-OR over its five raw color-prediction scores ${\mathit{p}}_i^j$:

5

The first row of the middle visualization described in the prior section contains scaled prediction ${\rho }_{\mathit{i}}^{\mathit{j}}$ (Figure 3B). In this scaling step, each raw color-prediction site score ${\mathit{p}}_i^j$ of site i is scaled by the color j predictions of all n sites in the same molecule:

The second row of the middle visualization (Figure 3B) and the lowest visualization (Figure 3C) described in the prior section contains scaled predictions ${\varrho }_i^j$. Because some types of reactions are rare, their raw predictions ${\mathit{p}}_i^j$ from the model are low and their corresponding color on the visualization scale are difficult for human eyes to differentiate. To improve the visibility of probabilistic differences between sites of the same reaction type, the raw prediction ${\mathit{p}}_i^j$ of site i is scaled by predictions on all N potential sites specific to that reaction type across the entire data set using logit and sigmoid (σ) functions:

The scaled predictions ${\varrho }_i^j$ for a reaction type are displayed only at its potential sites. These predictions are visualized using the red–blue color scheme, where red indicates probabilities higher than average and blue indicates probabilities lower than average. In addition, we split epoxidation and dehydrogenation reaction types into aromatic and aliphatic subtypes for visualization because these subtypes have drastically different background probabilities.

Reaction Coverage Comparison.

We compared our five-color metabolism model (Rainbow XenoSite) to several previously published phase I SOM prediction models: SMARTCyp, RSPredictor, SOMP, MetaPrint2D, FAME 2, Site of Metabolism Estimator (SOME), He et al., and a selection of single reaction type models (Table 7).^24-36 Side by side summary of the models’ methodology and data (Table 7) revealed several key observations encompassing interpretability, chemical space coverage, and reaction type coverage.

First, Rainbow XenoSite and single reaction type models are the only models that produce well-scaled, probabilistic outputs (Figure S1). Compared to ranking and binary outputs, probabilistic predictions are more useful for triaging purposes in drug development because medicinal chemists rely on estimates of confidence in a predicted outcome to direct their efforts and resources in the most efficient ways.

Second, Rainbow XenoSite covers the largest proportion of known human phase I metabolic reactions of all phase I models. We calculated the percentage contribution of each reaction type to the totality of phase I reactions in the AMD database. We estimated each model’s coverage of phase I metabolism by summing up the percentages corresponding to each reaction type included in that model. While previously published models cover up to 74.7% of phase I reactions, Rainbow XenoSite covers 92.3% of AMD phase I reactions (Table 7). Specifically, SMARTCyp, Xenosite 1.0, RSPredictor, SOMP, MetaPrint2D, and FAME 2 are five published generic SOM phase I prediction models.^24-30 SMARTCyp, XenoSite 1.0/RSPredictor, SOMP, MetaPrint2D, and FAME 2 account for 45.5%, 48%, 46.1%, 74.7%, and 48% of all human phase I reactions from the AMD database, respectively (Table 7). Similarly, Site of Metabolism Estimator (SOME)³⁵ and He et al.³⁶ are two published multireaction phase I prediction models. SOME and He et al. account for 42.6% and 37.6% of all human phase I reactions from the AMD database, respectively (Table 7). The large difference in reaction coverage and output format makes a direct performance comparison between Rainbow XenoSite and previously published models difficult.

Third, Rainbow XenoSite is a more comprehensive model of reaction types and includes many important bioactivation reactions missed by other models. Previously, we published several models that focused on specific reactions like epoxidation, quinone species formation, N-dealkylation, S-oxidation, and N-reduction.^31-34 These reactions play important roles in the bioactivation of multiple structural alerts leading to the formation of deleterious reactive metabolites.^33,34 Given a molecular structure, each of the published single-reaction models outputs probabilistic prediction scores corresponding to the likelihood that each atom or bond in the molecule is a site of metabolism for that specific reaction. Rainbow XenoSite can make predictions for all of these important bioactivation reactions. In contrast, the SMARTCyp, RSPredictor, SOMP, MetaPrint2D, FAME 2, SOME, and He et al. models do not comprehensively cover these important bioactivation reactions (Table 7).

Ambiguity of Generic SOMs.

We compared the rainbow-SOM labeling-scheme to the generic-SOM labeling-scheme in term of unambiguous mapping from the labels to the reported products. Given a labeled site, we can identify the reaction types that can occur at that site by assessing whether that site matches the reaction types’ potential sites defined by SMARTS strings (see Table S6). For sites that can undergo more than one type of reactions, the ambiguous sites, we quantified how often knowledge of the site of metabolism alone is sufficient or insufficient to infer the structure of the metabolite.

Most sites marked using generic-SOM labeling-scheme are ambiguous. As an exemplar of the generic SOM labels, we used scheme previously adopted by XenoSite (XenoSite 1.0).²⁶ XenoSite 1.0 used a generic-SOM labeling-scheme where atoms of a compound were labeled as SOM based on the structural change from substrate to product. In addition to XenoSite 1.0, published models like SMARTCyp, SOMP, MetaPrint2D, and FAME 2 also generically predict SOM and use a similar labeling-scheme.^24,27-30 A single generic SOM can map to multiple possible reactions of different types, leaving ambiguity as to which metabolite is generated by metabolism of the site. Of 20 472 SOMs identified using the generic atom labeling-scheme (Table S5), 18 544 can undergo more than one reaction type, yielding 90.58% ambiguous SOMs. We identified 20 major groupings of ambiguous SOMs under this labeling-scheme (Figure 4). Of 23 872 SOMs identified using the generic bond labeling-scheme where bonds between two heavy atoms are labeled as sites, 17 061 can undergo more than one reaction type, which modestly reduces the percentage of ambiguous SOMs to 71.47%.

Compared to the generic-SOM labeling-scheme, the rainbow-SOM labeling-scheme significantly reduces the ambiguity in mapping from the labels to reported products. In contrast to the generic-SOM approach, the rainbow-SOM labels specify both sites and reaction types by identifying modified bonds and the reaction classes of reactions observed at those bonds (Table S5). As a result, only 0.3% of rainbow labeled sites are ambiguous. Despite these improvements, the rainbow labels are ambiguous for a small subset of C-dealkylation reactions wherein oxygen can attach to either of the two carbons in the bond. The rainbow-SOM labeling-scheme nearly eliminates all ambiguity in mapping from the labels to reported products, which will enable metabolite structure prediction from effective rainbow-SOM models of metabolism.

Reaction-Class Specific Site-Level Accuracy.

Considering each color separately, Rainbow XenoSite accurately predicts SOMs within metabolized molecules. We used “top-two” and “average site AUC” to evaluate model accuracy. First, we measured top-two performance on a color by color basis, where a molecule was considered to be correctly predicted if any of its reported SOMs were ranked in the top two sites sorted by that class-specific prediction score. The top-two metric is commonly used to evaluate CYP site of metabolism models and is useful for assessing how the model will perform when comparing potential SOMs within individual drug candidates. Stable oxygenation, unstable oxygenation, dehydrogenation, reduction, and hydrolysis top-two accuracies were, respectively, 72.2%, 80.4%, 74.2%, 94.9%, and 94.7% (Figure 5). Across all five classes, the model had an average cross-validated top-two accuracy of 88.2%. Second, we measured average site AUC, a metric based on the area under the ROC curve (AUC). AUC is a standard metric employed in machine learning to quantify how well a model is able to distinguish between positives and negative cases of a specific chemical class. To compute average site AUC, individual AUCs are computed for each molecule in the data set, color by color. These AUCs are then averaged across the whole data set to produce a single summary statistic for each reaction class. Average site AUC accuracies for stable oxygenation, unstable oxygenation, dehydrogenation, reduction, and hydrolysis were 91.1%, 97.6%, 94.7%, 99.5%, and 99.7%, respectively. Across all five classes, the average AUC was 97.5% (Figure 5).

Considering all colors together, Rainbow XenoSite accurately predicts class-specific SOMs within metabolized molecules. In contrast to the previous accuracy metrics, here all colors are considered together, and the model’s ability to correctly predict the color and site at the same time is measured. In this evaluation, each color–site pair is considered an individual entity. All pairs are ranked by the scores assigned by the model and compared with the known metabolism to verify if the color-specific SOMs are separated from the rest. With all reaction color considered together, the model has top-two and average site AUC of 78.2% and 97.3%, respectively.

Rainbow XenoSite can accurately distinguish metabolized sites from all other potential sites across the entire data set. To evaluate this claim, we introduce the “potential site AUC” metric. In this metric, sites that are unable to undergo metabolism via a certain reaction color, i.e. not matching any SMARTS patterns that defining the corresponding potential site set, were removed from the calculation of the potential site AUC for that reaction color because they can be easily predicted as not-metabolized by both human experts and machine learning models. We then calculated AUCs for the remaining potential sites for each metabolism color.⁵⁵ Stable oxygenation, unstable oxygenation, dehydrogenation, reduction, and hydrolysis potential site AUC accuracies were 85.4%, 94.3%, 86.8%, 98.2%, and 93.9%, respectively. Across all five targets, the model had an average potential site AUC of 92.3% (Figure 5).

Rainbow XenoSite can also accurately identify which of the five colors of metabolism are most likely for a given site (Figure 6). For each metabolized site, we constructed computed AUCs from the scores assigned to each of the five metabolism colors. The “average reaction-specificity AUC” reflects the model’s ability to specify a reaction type for a metabolized site. The average reaction-specificity AUC was 97.4%.

Reaction-Type Specific Site-Level Accuracy.

In addition to making accurate predictions among the five colors of metabolism, Rainbow XenoSite can accurately separate metabolized sites from nonmetabolized sites for each of the 21 phase I reaction types included. To evaluate this claim, we computed metrics on both reaction-type specific potential sites only (PS) and across all training sites (global) (Table 8). Across 21 reaction types, our Rainbow XenoSite had an average cross-validated potential site AUC and global AUC of 85.7% and 99.0%, respectively.

Ideally, we would like to know how well our model performs in predicting reaction-types in comparison to existing models: SOME³⁵ and He et al.³⁶ using the same set of molecules. Unfortunately, there was no secure online web server available for the user to submit test molecules for either existing model. As a result, we can only report accuracies across all training sites from SOME³⁵ and He et al.³⁶ models on the same reaction types (Table 8). This comparison is inexact, since the specific molecules used to calculate performance metrics differ among the models.

Molecule-Level Model Accuracy.

Rainbow XenoSite accurately predicts whether a molecule undergoes specific classes of phase I reaction. In our training data set, each substrate was metabolized through one or more reaction classes but rarely through all five types. As a result, not all training molecules are a substrate of a certain reaction class. For each reaction class, we used the “molecule AUC” metric to evaluate how well the model separated class-specific substrates from nonsubstrates. Rainbow XenoSite has molecule AUCs of 78.3%, 83.9%, 77.3%, 90.4%, and 92.7% for stable oxygenation, unstable oxygenation, dehydrogenation, reduction, and hydrolysis, respectively (Figure 7).

Probability Scaling of Model Output.

The model’s outputs can be interpreted as probabilities. When we binned class-specific sites by the phase I prediction score, the proportion of class-specific SOMs in each bin closely correlated with the bin’s score (Figure S1). Likewise, when we binned molecules by the phase I molecule scores, the proportion of class-specific metabolized molecules in each bin also correlates with the bin score (Figure S2). Quantitatively, the Pearson regression coefficients of site and molecule levels were 0.996 and 0.947, respectively.

Reactive Metabolite Forming Reactions.

Predicting the formation electrophilically reactive metabolites is an important step in avoiding idiosyncratic adverse drug reactions.^56-59 Experimental methods for detecting reactive metabolites, most frequently in vitro incubation with glutathione and cyanide, are well-established.^60-62 However, these methods require significant time and resources, especially when screening thousands of compounds during drug development. Our previous works and others have shown the instrumental role of computational methods in directing these experimental efforts.^23,31,32,34 We published several models that focused on specific reactions like epoxidation,³¹ quinone species formation,³² S-oxidation,³⁴ nitroaromatic reduction,³⁴ and N-dealkylation.³³ These reactions play important roles the bioactivation of multiple structural alerts leading to the formation of deleterious reactive metabolites.³⁴ Given a molecular structure, each of these models outputs probabilistic prediction scores corresponding to the likelihood that each atom or bond in the molecule is a site of metabolism for that specific reaction. The following five sections demonstrate that Rainbow XenoSite can predict reactive metabolite formation reactions with high accuracy, comparable to previous models of individual reaction types.

Epoxidation.

Epoxides, the second most common type of reactive metabolites, are often formed via P450-mediated oxidation of double bonds or aromatic rings.^37,63 Epoxides are detoxified by epoxide hydrolases and the glutathione system.^64-67 When these pathways are overwhelmed, epoxides can covalently bind to proteins and cause toxic side effects and immune reactions.

Rainbow XenoSite can predict the formation of this key reactive metabolite and is comparable in accuracy to our previous epoxidation model.³¹ The epoxidation model was trained on a data set containing 524 molecules, 389 of which contained sites of epoxidation. Our model was able to predict this data set with an accuracy of 92.7% potential site AUC (Figure 8A), compared to the potential site AUC of 94.3% of the previous epoxidation model (two-sided p = 0.097). Reliability plot of unscaled stable oxygenation on potential epoxidation sites demonstrate that the scores can be interpreted as probabilities (Figure S3A).

One-Step Quinone Formation.

Quinone is the most common class of reactive metabolites, accounting for 40% of known reactive metabolites.⁶⁸ Quinone species are frequently highly electrophilically reactive.³⁸ Quinones are formed through one or two metabolic steps.³² In one-step quinone formation, the parent compound undergoes dehydrogenation to form the quinone metabolite. In two-step quinone formation, the parent compound undergoes a hydroxylation step following by a dehydrogenation step to form the quinone metabolite.

Rainbow XenoSite predicted the formation of this important reactive metabolite. Our model predicted a data set of 39 molecules that form quinone through one-step (dehydrogenation) with an accuracy of 90.1% AUC across all quinone formation potential sites (Figure 8B). This performance is comparable the potential site AUC of 92.5% of previous quinone formation model (two-sided p = 0.36).³² Reliability plot of unscaled dehydrogenation scores on potential quinone sites demonstrate that the scores are well-scaled probabilities (Figure S3B).

Unfortunately, Rainbow XenoSite does not predict two-step quinone formation reactions. Only a small subset of quinones are formed via one-step of metabolism. This is an intrinsic limitation of this approach, but it might be addressed by an algorithm that chain predictions of the Rainbow XenoSite into multiple steps. This possibility is out of scope of this study and will be considered future work.

N-Dealkylation.

Alkylated amines are a class of compounds with important biological and pharmaceutical functions.^69,70 N-dealkylation, the replacement of an alkyl group by hydrogen, is the main metabolism pathway for many alkylated amines.^71-73 This removed alkyl group is often transformed into a reactive aldehyde. Unfortunately, metabolic studies often focus on the amine metabolite and do not report the aldehyde metabolite because aldehydes are presumed to be quickly detoxified to carboxylic acids and alcohols by multiple enzyme systems. However, some reactive aldehydes escape these detoxification pathways, and form adducts that can give rise to adverse drug events, especially in patients where the detoxification pathways are inhibited by other drugs.^74,75

Rainbow XenoSite accurately predicts sites of N-dealkylation on the data set of 883 molecules used to train our previous N-dealkylation model.³³ Rainbow XenoSite predicted N-dealkylation with an accuracy of 92.2% global AUC (Figure 8C). This accuracy is comparable to the global AUC accuracy of 94.7% of our previous N-dealkylation model (two-sided p = 0.08). Reliability plot of unscaled unstable oxygenation on potential N-dealkylation sites demonstrate that the scores can be interpreted as probabilities (Figure S3C).

Nitroaromatics Reduction.

Nitroaromatics are common in the environment and are important pharmacophores in many drugs.^76-79 Unfortunately, nitroaromatic groups can metabolically transform into electrophilic, highly unstable intermediates such as nitroso and nitro anion radicals through reduction.^34,80-82

Rainbow XenoSite accurately predicts nitroaromatic reduction. Previously, we published the prediction from an early version of Rainbow XenoSite on nitroaromatic-containing drugs. The early version model was trained on a data set containing 3061 molecules, 98 of which contain nitroaromatics.³⁴ In comparison, the training data set of the current Rainbow XenoSite contains 314 nitroaromatic molecules. Rainbow XenoSite predicted nitro-reduction of the 98 molecules with an AUC accuracy of 90.1% (Figure 8D), compared to 93.0% obtained by the early version of the model (two-sided p = 0.16). Reliability plot of unscaled reduction on nitroaromatics demonstrate that the scores can be interpreted as probabilities (Figure S3D).

Thiophene S-Oxidation.

Thiophenes are five-membered sulfur-containing aromatic rings that are important pharmacophores in many drugs.^83-89 Thiophenes can undergo P450-mediated S-oxidation to form a reactive sulfur oxide, which can subsequently bind to proteins causing toxic side effects. For example, suprofen is a thiophene-containing analgesic medication that failed in late-stage clinical trials because of the formation of a short-lived reactive metabolite that caused irreversible inhibition of phase I enzymes.⁹⁰ However, many nontoxic drugs also contain thiophenes that are not bioactivated, such as eprosartan, a closely related analogue to suprofen.⁹¹

Rainbow XenoSite can predict thiophene S-oxidation. Previously, we published the predictions from an early version of Rainbow XenoSite on thiophene-containing compounds.³⁴ The training data set of the earlier model contained 50 thiophene containing compounds,³⁴ compared to 105 in the training data for the current study. Our current Rainbow XenoSite predicted thiophene S-oxidation with a global AUC accuracy of 86.0% (Figure 8D), compared to the previous model’s performance of 88.0% (two-sided p = 0.398). The reliability plot of unscaled stable oxygenation on thiophenes demonstrate that the scores can be interpreted as probabilities (Figure S3E).

Model Limitations.

Rainbow XenoSite has several limitations. First, our model does not cover all phase I reactions. Although rare, excluded reactions include tautomerization, isomerization, and rearrangement, which can have important implications on the efficacy and safety profile of drugs. In the future, improved models may be able to label and group these reaction types with a similar approach. Second, while the current SOM labeling scheme significantly reduced ambiguity, it does not provide sufficient information to predict the exact metabolites of a small number of reactions, about 0.3% of the total. In future studies, we could add an ordering to the labeling of bonds, enabling us to distinguish the exact site of oxygen addition in C-dealkylation reactions. Such an order-dependent SOM label would remove the ambiguity in these rare cases. Third, the model slightly underperformed previously developed models of specific reactive metabolite generating reactions, though the drop in performance was not statistically significant. It is possible that improved descriptors or training protocols could change this general trend. Fourth, Rainbow XenoSite only considers the metabolic aspect of IADRs. However, IADRs are complex phenomenon with other important aspects, including dosage, route of administration, competing bioactive and detoxification pathways, comorbidities, coadministered medications, and genetic variants.^92-94