Differentiating Inhibitors of Closely Related Protein Kinases with Single- or Multi-Target Activity via Explainable Machine Learning and Feature Analysis

Protein kinases are major drug targets. Most kinase inhibitors are directed against the adenosine triphosphate (ATP) cofactor binding site, which is largely conserved across the human kinome. Hence, such kinase inhibitors are often thought to be promiscuous. However, experimental evidence and activity data for publicly available kinase inhibitors indicate that this is not generally the case. We have investigated whether inhibitors of closely related human kinases with single- or multi-kinase activity can be differentiated on the basis of chemical structure. Therefore, a test system consisting of two distinct kinase triplets has been devised for which inhibitors with reported triple-kinase activities and corresponding single-kinase activities were assembled. Machine learning models derived on the basis of chemical structure distinguished between these multi- and single-kinase inhibitors with high accuracy. A model-independent explanatory approach was applied to identify structural features determining accurate predictions. For both kinase triplets, the analysis revealed decisive features contained in multi-kinase inhibitors. These features were found to be absent in corresponding single-kinase inhibitors, thus providing a rationale for successful machine learning. Mapping of features determining accurate predictions revealed that they formed coherent and chemically meaningful substructures that were characteristic of multi-kinase inhibitors compared with single-kinase inhibitors.


Introduction
In drug discovery, the ability of small molecules to interact with more than one protein in a well-defined manner provides the basis of polypharmacology; that is, the induction of the desired (or undesired) in vivo effects of drugs through the engagement of multiple targets [1][2][3][4][5]. Such multi-target activities of small molecules are a topic of intense investigation, from different perspectives [1][2][3][4][5][6][7][8][9]. Multi-target compounds (MT-CPDs) might be identified, for example, using profiling assays or proteomics techniques [10][11][12]. However, rationalizing multi-target activities of compounds (also referred to as promiscuity) at the molecular level of detail, distinguishing true activities from assay artifacts, and understanding how MT-CPDs might differ from single-target compounds (ST-CPDs) are far from being trivial tasks [13][14][15][16][17].
In addition to experimental methods, computational data analysis and predictive modeling are also applicable to aid in the analysis of the multi-target activities of small molecules [6][7][8][9]. For example, machine learning (ML) on the basis of chemical structure has been successfully used to systematically distinguish between MT-CPDs and corresponding ST-CPDs from medicinal chemistry or biological screening, thus providing evidence for the presence of structural features that differentiate MT-and ST-CPDs [18,19]. Furthermore, it mass of 1000 Da or more and potential assay interference compounds were removed using public tools and filters [16,37,38].

Target Selection
A search was carried out for kinase triplets including at least two closely related kinases from the same family for which sufficient numbers of MT-and corresponding ST-CPDs for meaningful ML and feature analysis were available (see Results and Discussion). Based on available inhibitors with high-confidence activity data and the applied selection criteria, two kinase triplets were prioritized, as reported in Table 1. In addition to focal adhesion kinase 1 (FAK1), triplet 1 contained two Januskinases, that is, Janus kinase 2 (JAK2) and 3 (JAK3). Triplet 2 comprised three closely related dual specificity tyrosinephosphorylation-regulated kinases (DYRK1A, DYRK1B, and DYRK2). For triplet 1, larger numbers of compounds were available than for triplet 2. A generally limiting factor for triplet assembly was the limited availability of sufficient numbers of MT-CPDs with highconfidence activity data, which were essential for the analysis. Nine of the ST-CPDs for FAK1 (triplet 1) were designated allosteric compounds. All other inhibitors for both triplets, including all MT-CPDs, were ATP-competitive.

Molecular Representation
As a molecular representation, atom environments were selected as preferred topological features [20,39]. The RDKit [37] implementation of the Morgan fingerprint corresponding to the extended connectivity fingerprint [39] was utilized to generate hash values of molecule-specific layered atom environments (up to a bond radius of 2, corresponding to a bond diameter of 4) for each atom in a compound. Obtained feature hashes were assigned to unique positions in the final feature vector to avoid bit collisions, thereby ensuring the interpretability of calculated feature importance values.

Machine Learning
Compounds were classified using a balanced random forest (BRF) model consisting of an ensemble of decision trees [40,41]. For each tree, a unique bootstrap sample of the training set was drawn and subsequently balanced by randomly under-sampling the majority class. This approach allowed us to utilize the majority of training data under the condition of class balance; an important criterion, given the presence of significantly different numbers of MT-and corresponding ST-CPDs. Predicted probabilities for multi-target activity were calculated as the mean probability over individual trees, which estimated the class probability as the fraction of samples of the given class in the final leaf node. Hyperparameters such as number of decision trees ("n_estimators": 25, 50, 100, 200, 400), minimal number of samples for a split ("min_samples_split": 2, 3, 5, 10), and minimum number of samples for a leaf node ("min_samples_leaf": 1, 2, 5, 10) were assessed via internal 10-fold shuffle-split cross validation on the training set. The final MT-vs. ST-CPD classifier was trained with the best-performing hyperparameter combination using the complete training set, representing a random sample of 75% of the compounds. Model performances were assessed using a balanced sample of the remaining 25% of MT-and ST-CPDs as test instances over 10 individual trials. As performance measures, balanced accuracy (BA) [42], F1-score (F1) [43], precision, recall, and Matthews correlation coefficient (MCC) [44] values were calculated. These performance measures are defined as follows: BA = 1/2(TPR + TNR) TP, TN, FP, and FN stand for true positives, true negatives, false positives, and false negatives, respectively.

SHAP Analysis and Feature Extraction
The use of BRF models enabled the accurate calculation of SHAP feature importance values for individual predictions using the TreeExplainer algorithm with tree-path dependent feature perturbations [30]. SHAP theory is provided in the Supplementary Methods, and cumulative SHAP feature contributions yielding a class label probability are illustrated in Supplementary Figure S1.
A feature extraction scheme was devised for correctly predicting instances taking into account that SHAP feature importance values might differ from test instance to test instance [31]. This feature extraction scheme bridges instance predictions and SHAP feature importance across all test instances, as follows: (i) For each correctly predicted MT-CPD, the top-ranked N features with the highest SHAP values were pre-selected and these features were pooled. (ii) The pool of the top-ranked N features was re-ranked by the feature frequency of occurrence in correctly predicted MT-CPDs, and the top M most frequent features were selected.
For the calculations reported herein, N = 5 and M = 10 settings were consistently applied.

Study Design
The newly generated kinase triplet test system with available MT-and ST-CPDs, as described below, enabled us to first address the key question if multi-kinase inhibitors could be systematically distinguished from single-kinase inhibitors by ML on the basis of chemical structure information. Given the frequent assumption that kinase inhibitors tend to be promiscuous, this was not necessarily likely. If accurate classification of MT-vs. ST-CPDs is possible, however, then structural features detectable by ML must exist that differentiate MT-and ST-CPDs and hence determine accurate predictions. If so, the second step of the analysis then aims at identifying these features via an independent explanatory approach (SHAP). Whether or not features determining algorithmic predictions might be chemically relevant and explainable in chemical terms was another open question. Therefore, in the third step, we aimed at rationalizing distinguishing features (provided they were identified) from a chemical perspective. Hence, the analysis was designed to identify the structural features driving the correct prediction of ST-and/or MT-CPDs, which might also be implicated in kinase selectivity or promiscuity, respectively.

Systematic Analysis of Kinase Triplets
The 489 human kinases with available active compounds were systematically organized into triplets, yielding nearly 20 million (19,368,964) unique combinations. Of these possible combinations, 6,132,688 were found to share at least one inhibitor. Figure 1 reports the distribution of the number of MT-CPDs for these triplets, revealing that~75% of all triplets had no more than two MT-CPDs and only 64 triplets had at least 50 MT-CPDs. Hence, confirmed MT-CPDs were generally rare.
step of the analysis then aims at identifying these features via an independent explanatory approach (SHAP). Whether or not features determining algorithmic predictions might be chemically relevant and explainable in chemical terms was another open question. Therefore, in the third step, we aimed at rationalizing distinguishing features (provided they were identified) from a chemical perspective. Hence, the analysis was designed to identify the structural features driving the correct prediction of ST-and/or MT-CPDs, which might also be implicated in kinase selectivity or promiscuity, respectively.

Systematic Analysis of Kinase Triplets
The 489 human kinases with available active compounds were systematically organized into triplets, yielding nearly 20 million (19,368,964) unique combinations. Of these possible combinations, 6,132,688 were found to share at least one inhibitor. Figure 1 reports the distribution of the number of MT-CPDs for these triplets, revealing that ~75% of all triplets had no more than two MT-CPDs and only 64 triplets had at least 50 MT-CPDs. Hence, confirmed MT-CPDs were generally rare. As a minimal amount of negative data, triplets were required to have at least 17 ST-CPDs for each kinase, reducing the number of triplets to 57. For 35 of these triplets, a mean balanced accuracy of BRF models greater than 80% was observed (see below) and SHAP calculations prioritized features determining the predictions. For the 35 triplets, the number of MT-CPDs ranged from 51 to 310. As representative triplets for subsequent analysis, triplet 1 with very high prediction accuracy and a large number of available MT-CPDs, and triplet 2 with lower prediction accuracy and a smaller number of MT-CPDs were chosen. In both instances, features decisive for the predictions were clearly interpretable in chemical terms (which is not necessarily the case in ML). Figure 2 summarizes the performance of our BRF models. For both kinase triplets, MT-and ST-CPDs were distinguished with surprisingly high accuracy, as determined on the basis of different performance measures. The calculations were generally stable, as reflected by the narrow distributions of the results obtained over independent trials. For triplet 2, prediction accuracy was consistently above 80%. However, for triplet 1, the predictions were nearly perfect, with values of all performance measures approaching 1.0 (the trial set-ups and results were thoroughly re-examined, excluding the presence of artifacts for triplet 1). Taken together, these findings provided evidence for the presence of distinguishing structural features and an unexpectedly solid foundation for subsequent feature analysis. As a minimal amount of negative data, triplets were required to have at least 17 ST-CPDs for each kinase, reducing the number of triplets to 57. For 35 of these triplets, a mean balanced accuracy of BRF models greater than 80% was observed (see below) and SHAP calculations prioritized features determining the predictions. For the 35 triplets, the number of MT-CPDs ranged from 51 to 310. As representative triplets for subsequent analysis, triplet 1 with very high prediction accuracy and a large number of available MT-CPDs, and triplet 2 with lower prediction accuracy and a smaller number of MT-CPDs were chosen. In both instances, features decisive for the predictions were clearly interpretable in chemical terms (which is not necessarily the case in ML). Figure 2 summarizes the performance of our BRF models. For both kinase triplets, MT-and ST-CPDs were distinguished with surprisingly high accuracy, as determined on the basis of different performance measures. The calculations were generally stable, as reflected by the narrow distributions of the results obtained over independent trials. For triplet 2, prediction accuracy was consistently above 80%. However, for triplet 1, the predictions were nearly perfect, with values of all performance measures approaching 1.0 (the trial set-ups and results were thoroughly re-examined, excluding the presence of artifacts for triplet 1). Taken together, these findings provided evidence for the presence of distinguishing structural features and an unexpectedly solid foundation for subsequent feature analysis. Biomolecules 2022, 12, x 6 of 12

Representation Features Determining Predictions
On the basis of the BFR results, SHAP analysis was carried out for each correctly predicted MT-and ST-CPD, prioritized features were extracted (see Materials and Methods), and their contributions to accurate predictions were quantified. The analysis was carried out for all representation (fingerprint) features that were present in test compounds as well as for features that were absent, thus comprehensively searching for features determining correct predictions. Figure 3 shows the results of SHAP feature importance analysis. For both triplets 1 and 2, a clear and consistent picture emerged from the analysis. Accurate predictions of MT-CPDs were determined by the features that were present in these compounds. These features made large positive contributions, whereas features absent in MT-CPDs made only small positive or negative contributions to the predictions (which essentially canceled out). By contrast, correct predictions of ST-CPDs were largely determined by features that were absent in these compounds (but present in MT-CPDs). In this case, present features only made small supporting contributions (i.e., negative in the case of the ST-CPD class) or opposing (positive) contributions. These observations paralleled our previous findings for MT-and ST-CPDs with activity against pairs of unrelated targets [31]. Thus, for kinase inhibitors, ML successfully distinguished between MT-and ST-CPDs on the basis of structural features that were unique to MT-CPDs.

Representation Features Determining Predictions
On the basis of the BFR results, SHAP analysis was carried out for each correctly predicted MT-and ST-CPD, prioritized features were extracted (see Materials and Methods), and their contributions to accurate predictions were quantified. The analysis was carried out for all representation (fingerprint) features that were present in test compounds as well as for features that were absent, thus comprehensively searching for features determining correct predictions. Figure 3 shows the results of SHAP feature importance analysis. For both triplets 1 and 2, a clear and consistent picture emerged from the analysis. Accurate predictions of MT-CPDs were determined by the features that were present in these compounds. These features made large positive contributions, whereas features absent in MT-CPDs made only small positive or negative contributions to the predictions (which essentially canceled out). By contrast, correct predictions of ST-CPDs were largely determined by features that were absent in these compounds (but present in MT-CPDs). In this case, present features only made small supporting contributions (i.e., negative in the case of the ST-CPD class) or opposing (positive) contributions. These observations paralleled our previous findings for MT-and ST-CPDs with activity against pairs of unrelated targets [31]. Thus, for kinase inhibitors, ML successfully distinguished between MT-and ST-CPDs on the basis of structural features that were unique to MT-CPDs. Biomolecules 2022, 12, x 7 of 12

Feature Mapping and Rationalization
After extracting and mapping individual features determining the accurate prediction of MT-CPDs, we annotated atoms of MT-CPDs with SHAP values from all features present in the respective compounds (including extracted features) for further analysis. Using layered atom environments (consisting of atom sets) as representation features made it possible to unambiguously map these contributions on a per-atom basis. For MT-CPDs from both triplets 1 and 2, highlighted regions were not evenly distributed over the compound structure, but delineated a coherent substructure, as shown in Figure 4 and Figure 5, respectively. Hence, the mapping of structural features determining the accurate prediction of MT-CPDs identified a well-defined structural motif. In both cases, this structural motif was predominantly formed by extracted (prioritized) features. In ST-CPDs, similar structural features influencing the predictions were not detected, as discussed above. For triplet 1, the delineated substructure was a [1,2,4]triazolo[1,5-a]pyridine, as depicted in Figure 4.

Feature Mapping and Rationalization
After extracting and mapping individual features determining the accurate prediction of MT-CPDs, we annotated atoms of MT-CPDs with SHAP values from all features present in the respective compounds (including extracted features) for further analysis. Using layered atom environments (consisting of atom sets) as representation features made it possible to unambiguously map these contributions on a per-atom basis. For MT-CPDs from both triplets 1 and 2, highlighted regions were not evenly distributed over the compound structure, but delineated a coherent substructure, as shown in Figures 4 and 5, respectively. Hence, the mapping of structural features determining the accurate prediction of MT-CPDs identified a well-defined structural motif. In both cases, this structural motif was predominantly formed by extracted (prioritized) features. In ST-CPDs, similar structural features influencing the predictions were not detected, as discussed above. For triplet 1, the delineated substructure was a [1,2,4]triazolo[1,5-a]pyridine, as depicted in Figure 4.  Thus, the substructure was characteristic of triplet 1 MT-CPDs. This finding also explained the decisive role of extracted features defining this substructure for accurate predictions. All 210 MT-CPDs containing this substructure originated from a single patent source establishing triple-kinase activity [45]. It follows that the 21 ST-CPDs containing the [1,2,4]triazolo [1,5-a]pyridine might also have triple-kinase activity. On the other hand, since only 21 of 2433 ST-CPDs contained this substructure, it is very likely to represent a chemical signature of MT-CPDs.
For triplet 2, the substructure in MT-CPDs identified by feature mapping was imidazo(1,2-b)pyridazine, as depicted in Figure 5. In this case, features present in MT-CPDs were more widely distributed across the compound structure, while the imidazo(1,2-    Thus, the substructure was characteristic of triplet 1 MT-CPDs. This finding also explained the decisive role of extracted features defining this substructure for accurate predictions. All 210 MT-CPDs containing this substructure originated from a single patent source establishing triple-kinase activity [45]. It follows that the 21 ST-CPDs containing the [1,2,4]triazolo[1,5-a]pyridine might also have triple-kinase activity. On the other hand, since only 21 of 2433 ST-CPDs contained this substructure, it is very likely to represent a chemical signature of MT-CPDs. For triplet 2, the substructure in MT-CPDs identified by feature mapping was imidazo(1,2-b)pyridazine, as depicted in Figure 5. In this case, features present in MT-CPDs were more widely distributed across the compound structure, while the imidazo(1,2- For triplet 2, the substructure in MT-CPDs identified by feature mapping was imidazo(1,2b)pyridazine, as depicted in Figure 5. In this case, features present in MT-CPDs were more widely distributed across the compound structure, while the imidazo(1,2-b)pyridazine was highlighted by mapped SHAP values. This was consistent with the observation that this substructure was contained in 42 of 74 MT-CPDs (which also originated from a single study establishing triple-kinase activity [46]). For the remaining 32 MT-CPDs, extracted features did not delineate another well-defined structural motif. However, the imidazo(1,2-b)pyridazine substructure was also characteristic of the large subset of 42 MT-CPDs because it did not occur in any DYRK1B or DYRK2 ST-CPD, and only in 1 of 342 DYRK1A ST-CPDs.
For both characteristic substructures, we also found an X-ray structure of a complex formed between an inhibitor containing the substructure and a kinase from triple 1 and triple 2, respectively. Figure 6 shows kinase-inhibitor interaction diagrams computed from these X-ray structures. The [1,2,4]triazolo [1,5-a]pyridine moiety in the inhibitor in Figure 6a was involved in multiple interactions with JAK2, forming the center of interactions for this inhibitor in the ATP binding site of the kinase. The imidazo(1,2-b)pyridazine moiety of the inhibitor in Figure 6b also formed interactions with DYRKA1 (but was not an interaction hot spot). Both characteristic substructures were located in the same region of the ATP site (essentially mimicking the adenosine moiety in ATP) and also interacted with residues conserved in the ATP site of other kinases, consistent with the multi-kinase activity of inhibitors containing these substructures.
Biomolecules 2022, 12, x 9 of 12 b)pyridazine was highlighted by mapped SHAP values. This was consistent with the observation that this substructure was contained in 42 of 74 MT-CPDs (which also originated from a single study establishing triple-kinase activity [46]). For the remaining 32 MT-CPDs, extracted features did not delineate another well-defined structural motif. However, the imidazo(1,2-b)pyridazine substructure was also characteristic of the large subset of 42 MT-CPDs because it did not occur in any DYRK1B or DYRK2 ST-CPD, and only in 1 of 342 DYRK1A ST-CPDs. For both characteristic substructures, we also found an X-ray structure of a complex formed between an inhibitor containing the substructure and a kinase from triple 1 and triple 2, respectively. Figure 6 shows kinase-inhibitor interaction diagrams computed from these X-ray structures. The [1,2,4]triazolo [1,5-a]pyridine moiety in the inhibitor in Figure 6a was involved in multiple interactions with JAK2, forming the center of interactions for this inhibitor in the ATP binding site of the kinase. The imidazo(1,2-b)pyridazine moiety of the inhibitor in Figure 6b also formed interactions with DYRKA1 (but was not an interaction hot spot). Both characteristic substructures were located in the same region of the ATP site (essentially mimicking the adenosine moiety in ATP) and also interacted with residues conserved in the ATP site of other kinases, consistent with the multi-kinase activity of inhibitors containing these substructures.

Conclusions
In this work, we attempted to differentiate inhibitors with triple-kinase or corresponding single-kinase activity by ML on the basis of chemical structure, identify features determining accurate predictions, and interpret key features in chemical terms. Exploring the molecular origins of the varying promiscuity of ATP site-directed kinase inhibitors

Conclusions
In this work, we attempted to differentiate inhibitors with triple-kinase or corresponding single-kinase activity by ML on the basis of chemical structure, identify features determining accurate predictions, and interpret key features in chemical terms. Exploring the molecular origins of the varying promiscuity of ATP site-directed kinase inhibitors continues to be a topic of intense investigation in medicinal chemistry. For our analysis, we generated a test system consisting of kinase triplets for which sufficient numbers of MT-and ST-CPDs with high-confidence activity data were available to enable meaningful ML. Moreover, we specifically aimed to investigate closely related kinases from the same family most likely to have similar compound-binding characteristics. MT-and ST-CPDs from kinase triplet 1 and 2 were differentiated with high accuracy using ML models, providing evidence for the presence of distinguishing structural features. SHAP analysis then identified features determining the predictions. An important finding shows that accurate predictions resulted from features that were present in MT-but absent in ST-CPDs. These features were found to be chemically sensible, forming coherent substructures that were characteristic of MT-CPDs. ML prediction accuracy was nearly perfect for kinase triplet 1 and in this case, the characteristic substructure was present in 94% of all MT-and absent in 99% of all ST-CPDs, thus reflecting the high consistency of ML results and feature analysis. Taken together, the findings reported herein have methodological implications as well as implications for kinase inhibitor research. From a methodological point of view, our results clearly support the utility of explainable ML to rationalize predictions from a chemical or biological perspective and reveal structural information important for drug discovery and design, as exemplified by the identification of substructures characteristic of MT-CPDs. However, despite the consistency of the obtained results, they are difficult to generalize for kinase inhibitor research. ML and feature analysis, at least in part, depend on the composition of the investigated data sets and care must be taken not to over-interpret the findings. For example, while we can conclude with certainty from our analysis that the [1,2,4]triazolo[1,5-a]pyridine and imidazo(1,2-b)pyridazine moieties identified herein are signatures of multi-kinase activity, we cannot conclude that the designated ST-CPDs assembled on the basis of currently available activity data are kinase-selective. Here, varying test frequencies among kinase inhibitors might come into play that are for the most part unknown for compounds collected from different sources and thus cannot be considered in the computational analysis. On the other hand, it makes perfect sense that designated ST-CPDs do not share characteristic structural features determining their prediction. The presence of such features in ST-CPDs would principally not be consistent with kinase selectivity, while the absence of shared features supports differences between these compounds. Clearly, such considerations are important for putting the results into perspective. However, explainable ML as presented herein yields, at the very least, experimentally testable hypotheses for distinguishing between inhibitors with single-and multi-kinase activity and for exploring further structural features implicated in promiscuity vs. selectivity.