Due to relatively high costs and labor required for experimental profiling of the full target space of chemical compounds, various machine learning models have been proposed as cost-effective means to advance this process in terms of predicting the most potent compound-target interactions for subsequent verification. However, most of the model predictions lack direct experimental validation in the laboratory, making their practical benefits for drug discovery or repurposing applications largely unknown. Here, we therefore introduce and carefully test a systematic computational-experimental framework for the prediction and pre-clinical verification of drug-target interactions using a well-established kernel-based regression algorithm as the prediction model. To evaluate its performance, we first predicted unmeasured binding affinities in a large-scale kinase inhibitor profiling study, and then experimentally tested 100 compound-kinase pairs. The relatively high correlation of 0. 77 (p < 0. 0001) between the predicted and measured bioactivities supports the potential of the model for filling the experimental gaps in existing compound-target interaction maps. Further, we subjected the model to a more challenging task of predicting target interactions for such a new candidate drug compound that lacks prior binding profile information. As a specific case study, we used tivozanib, an investigational VEGF receptor inhibitor with currently unknown off-target profile. Among 7 kinases with high predicted affinity, we experimentally validated 4 new off-targets of tivozanib, namely the Src-family kinases FRK and FYN A, the non-receptor tyrosine kinase ABL1, and the serine/threonine kinase SLK. Our sub-sequent experimental validation protocol effectively avoids any possible information leakage between the training and validation data, and therefore enables rigorous model validation for practical applications. These results demonstrate that the kernel-based modeling approach offers practical benefits for probing novel insights into the mode of action of investigational compounds, and for the identification of new target selectivities for drug repurposing applications.
Deregulated kinase activity plays a role in many diseases, hence calling for therapeutic compounds that could effectively inhibit specific members of the protein kinome. Although kinase inhibitors form the largest group of new drugs approved for cancer treatment [1], a majority of them are ATP-competitive, and therefore present a highly promiscuous mechanism of action (MoA), due to the high evolutionary conservation of the kinase ATP-binding pockets [2,3]. The polypharmacological interactions contribute both to therapeutic and toxic responses seen in clinically-approved and investigational kinase inhibitors. Thus, improved knowledge of the complex compound-target binding interactions across the full protein kinome, including both on- and off-target effects, is of high clinical relevance for future drug discovery applications. Recent technological advances in chemoproteomic approaches, such as thermal profiling [4], have enabled efficient determination of kinome-wide compound potency. Several commercial providers are available for preclinical kinase inhibitor testing in vitro, including DiscoverX, Millipore and Reaction Biology. Even though the experimental compound-target interaction mapping is critical to characterizing a compound’s MoA, computational methods provide a complementary and cost-effective approach with the potential to accelerate the exploration of the enormous size of the chemical universe, estimated to consist of approximately 1020 molecules exhibiting good pharmacological properties [5]. The hypothesis is that in silico models could provide fast, large-scale and systematic pre-screening of chemical probes, toward prioritization of the most potent interactions for further in vitro or ex vivo verification in the laboratory [6–10]. In particular, a lot of work has been devoted to compound-based interaction prediction methods, including quantitative structure-activity relationship (QSAR) models, which aim to relate structural properties of the chemical molecules to their bioactivity profiles [11,12]. Another class of machine learning methods, so called target-based methods, focus on evaluating similarities between amino acid sequences or three-dimensional structures of protein targets [13]. In these supervised learning approaches, models are trained using available bioactivity data, together with either compound or protein information, which allows then predicting either new target interactions for a given drug or new drugs targeting a given protein. Furthermore, such methods typically focus on a limited set of molecules of interest. As a more recent class of computational modelling approaches, systems-based frameworks take advantage of the information available on both compounds and targets. For instance, Yamanishi et al. proposed a supervised machine learning approach for categorizing drug-target pairs as interacting or non-interacting based on an integrated model of chemical and genomic molecular profiles [14]. Since this seminal work, a wide variety of systems-based prediction methods have been developed that utilize various molecular descriptors and learning techniques, including random forest, neural networks and kernel learning [15–32]. Even though such models may hold a great potential, their computational predictions are rarely being directly verified in the laboratory and, consequently, their practical benefits for the drug discovery or repurposing applications remain largely unknown. Toward testing the practical potential of systems-based machine learning models, we implemented a computational-experimental framework for prediction and verification of compound-target bioactivity profiles (Fig 1). We focused on a regression problem, where the task is to predict the actual binding affinities, instead of the standard bioactivity classification setting that treats molecular interactions as simple on/off relationships. As a prediction model, we applied a well-established kernel-based regularized least squares learning algorithm (KronRLS [33]), because kernels, in addition to offering a computationally efficient means for increasing the power of linear learning algorithms, are particularly well-suited for capturing and learning complex molecular properties for prediction purposes [34,35]. The specific contributions of the work are the following. First, we evaluated a large number of molecular descriptors in the form of kernels, including our novel, extended target profile-based protein kernel, and a generic string kernel that has not previously been used in the context of compound-protein interaction prediction. Second, we show how these kernels guided us in filling the gaps in a large-scale compound-kinase interaction map [3]. Third, we experimentally tested a subset of 100 predicted binding affinities, achieving a high correlation of 0. 77 between the measured and predicted bioactivities. Finally, we demonstrate the potential of the modelling approach in a more challenging task of predicting target selectivities for such a new candidate compound that has no bioactivity data available for model training. As a specific case study, we used an investigational tyrosine kinase inhibitor tivozanib whose established target profile consists only of 3 on-targets. We experimentally tested and validated 4 out of 7 kinases predicted as tivozanib’s off-targets, providing novel insights into its MoA, and thereby extending the potential therapeutic target space of tivozanib.
A key assumption in the systems-based compound-target interaction prediction algorithm is that similar drug compounds are likely to bind to similar protein targets, and therefore the first challenge lies in the representation and use of molecular similarities in the most predictive way. We encoded here similarities between drugs and similarities between proteins using different types of kernels, constructed based on chemical two- and three-dimensional structures, amino acid sequences, protein structures, and molecular interaction profiles (see Materials and Methods for details). Such systematic construction of chemical and genomic molecular descriptors resulted in 12 drug kernels and 8 protein kernels. To predict compound-protein binding affinities using the regression setup, we applied a regularized least squares (RLS) model for each pair of drug kernel and protein kernel (KronRLS algorithm, see Materials and Methods). For computational evaluation of the predictive performance of various molecular descriptors and optimization of model parameters under separate prediction scenarios, we carried out two systematic nested cross-validation (CV) procedures, i. e. , leave-one-out cross-validation (LOO-CV, S12 Fig) and leave-drug-out cross-validation (LDO-CV, S13 Fig), using 16,265 known binding affinities (pKi values) between 152 kinase inhibitors and 138 protein kinases measured in a large-scale functional bioassay by Metz et al. [3] (S9 Fig). We applied LOO-CV to tune model parameters and to evaluate its predictive performance when filling experimental gaps in large-scale target profiling studies (the Bioactivity Imputation scenario, Fig 2A), and LDO-CV in the inference of target interactions for a new candidate drug compound (the New Drug scenario, Fig 2B). LOO-CV corresponds to the design where scattered missing values are present in otherwise known compound-protein bioactivity matrix, and the aim is to predict the missing entries within the training data. LDO-CV, on the other hand, simulates more challenging inference problem, in which the aim is to predict targets of an investigational drug compound, not encountered in the training data (Materials and Methods, Tables 1 and S1). Next, we trained the KronRLS algorithm with 16,265 bioactivities between 152 kinase inhibitors and 138 kinases measured in the study by Metz et al. [3], together with the best-performing under the Bioactivity Imputation scenario drug interaction profile kernel (KD-GIP) and kinase domain-based generic string protein kernel (KP-GS-domain). We then used the optimized model to predict the remaining 4,711 binding affinities that were missing in this experimentally-measured compound-kinase interaction map (S9 Fig). To assess the model’s practical utility, we experimentally tested a set of 100 predicted binding affinities between 5 drug compounds (cediranib, lapatinib, gefitinib, pazopanib and vx-745) and 20 kinases (ABL1, AXL, BRK, BTK, EGFR, FAK, FYN A, HER2, HER4, IGF1R, InsR, ITK, JAK3, KDR, LCK, LYN B, PYK2, SRC, SYK, TRKA). Among these, new potential interactions, not present in the Metz et al. dataset, were predicted for cediranib, lapatinib and gefitinib (S4 Fig). We note that pazopanib presented very high prediction accuracy, despite of having a sparse binding affinity profile available for the model training. On the other hand, vx-745 had no potent activities, either measured or predicted, against any of the kinases, and therefore it served as a negative control in the validation (S4 Fig). We tested the predicted bioactivities using a cell-free ADP-Glo Kinase Assay (see Materials and Methods for details). We observed a relatively high Pearson correlation of 0. 774 (p < 0. 0001) between the model-predicted (pKi) and experimentally-measured (pIC50) bioactivities among the 100 compound-kinase pairs (Fig 4A). The IC50 readout from our assay, similarly to the inhibition constant Ki in the Metz et al. study, indicates the concentration of the compound needed to inhibit enzymatic activity of a kinase by 50%. Even though IC50 is known to depend on the concentration of the enzyme, inhibitor, and substrate, along with other experimental conditions, whereas Ki is an intrinsic, thermodynamic quantity independent of the substrate [36], recent studies have shown a sufficiently high level of association between pIC50 and pKi readouts, permitting their reliable comparison [37,38]. We also observed a strong technical correlation of 0. 769 (p < 0. 0001) between pIC50 readouts from our profiling assay and pKi values measured in the study by Metz et al. (S5 Fig), which supports the feasibility of our experimental validations. Furthermore, based on the published information [39], the ATP concentration used in our assay (10 μM) is expected to be below, or in some cases equal to, the ATP Km values of the kinases tested, suggesting that the IC50 values should be very close to the respective Ki values. We also compared both the model-predicted and experimentally-measured interaction mapping to the results from another large-scale binding assay by Davis et al [40]. In this study, 72 clinically relevant kinase inhibitors were profiled against 442 kinases, providing, for each compound-kinase pair tested, dissociation constant Kd, indicating the tendency of a larger molecular complex to dissociate reversibly into the component molecules. We again observed a very good agreement (correlation of 0. 796, p < 0. 0001) between the computationally-predicted (pKi) and measured (pKd) binding affinities across the overlapping 95 compound-kinase pairs (S5 Fig). We noted even a higher technical correlation of 0. 916 (p < 0. 0001) between pKd values from Davis et al. study and pIC50 values from our experimental assay (S5 Fig). Notably, a comparison between predicted pKi values and measured pKd readouts across a larger set of 2,662 compound-kinase pairs overlapping between Metz et al. [3] and Davis et al. [40] studies resulted in a lower correlation (0. 642, p < 0. 0001, S6 Fig), compared to that when considering only the pairs included in our experimental assay (0. 796, p < 0. 0001, S5 Fig). A high pKd indicates that a substrate is more likely to be bound to an enzyme, whereas pKi measures the potency of a drug. Even though both high pKi and pKd values are considered as indicators of drug activity, a drug with high pKi does not necessarily result in a high pKd (S6 Fig). As expected, the correlation between the model-predicted and measured by Metz et al. pKi values in the training data (excluding pairs blinded in the model training, marked with an orange cross points in Fig 4A) was somewhat higher than that for the missing compound-kinase pairs (correlation of 0. 802, p < 0. 0001, S5 Fig). Of note, our experimental assay confirmed computationally-predicted high binding affinities between cediranib-KDR, lapatinib-EGFR and pazopanib-KDR, the two first of which were even not measured in the study of Metz et al. Further, although some bioactivities missing in Metz et al. dataset corresponded to compound-kinase pairs already tested in other studies (e. g. lapatinib-EGFR measured in the assay by Davis et al.), these were not used in the training of our model. Taken together, the observed high agreement between the predicted and experimentally-measured bioactivities demonstrates the potential of the kernel-based modeling framework with appropriately chosen kernels for filling the gaps in existing compound-target interaction maps. Finally, we tested whether the optimized model can also predict target interactions for a new chemical probe, which has no profiling data available for the model training. We used here tivozanib as an example of an investigational tyrosine kinase inhibitor, known to be potent towards all three VEGF receptors (FLT1, KDR, FLT4) [41]. Beyond VEGFRs, however, the target profile of tivozanib has otherwise remained poorly characterized, including its potential off-targets. We therefore again used 16,265 binding affinities between 152 compounds and 138 kinases measured in the study by Metz et al. to train the KronRLS model with shortest paths between atoms-based drug kernel (KD-sp) and amino acid sequence-based generic string protein kernel (KP-GS), which were found to perform best under the New Drug scenario. Since the model should always be tuned separately under distinct prediction scenarios, even if the training dataset is the same, the model used here and the one described in the previous section differ in their chosen kernels and optimized value of the regularization parameter. With the optimized model, we predicted the bioactivity of tivozanib against the set of 138 kinases (S3 Table). As the first positive control, the model correctly predicted high potency of tivozanib against its known on-targets FLT1, KDR and FLT4 (Fig 5A and S3 Table). To further assess the quality of the predictions, we used publicly available bioactivity data from the study of Gao et al. who profiled 158 kinase inhibitors, including tivozanib, for their inhibitory activity at 1 μM and 10 μM against 234 kinases [42]. Although the concentrations adopted in this screen were too high for pre-clinical testing of positive interactions, we used these data to evaluate the negative predictions from the model. In total, 64 out of 82 kinases with low predicted affinities (pKi < 6 M) were screened by Gao et al. Among these, 59 kinases (92%) have at least 50% of the activity remaining at the high compound concentration of 1 μM (S3 Table), thus effectively validating the model’s negative predictions (Fig 5B). We next went on and tested experimentally 7 predicted off-target interactions (ABL1, Aurora A, FRK, FYN A, HIPK4, RPS6KB1, SLK). These 7 kinases were selected among the set of 25 kinases with the highest predicted binding affinities by focusing on off-targets unique to tivozanib. Specifically, we compared the predicted target interaction profile of tivozanib to other VEGFR inhibitors found in the ChEMBL database [43]. For instance, RET was not selected, even though it was predicted to have high potency towards tivozanib (pIC50 = 6. 9 M), since it is targeted by 76% of VEGFR inhibitors in ChEMBL (potency of at most 100 nM), whereas FYN A was included in our experimental assay because it is targeted only by 26% of VEGFR inhibitors (S3 Table). We tested the predicted bioactivities using a cell-free ADP-Glo Kinase Assay (see Materials and Methods for details). Among the pre-selected off-target predictions, our experiments confirmed strong binding affinity between tivozanib and 4 out of the 7 tested kinases (57%), namely FRK, ABL1, SLK and FYN A (Fig 5A). Statistical significance of this success rate depends on the underlying distribution of the true target space of tivozanib, which is unknown. However, if one assumes that no more than 18 of 138 considered kinases are actual targets of tivozanib (13%), then the observed overlap is significant (p < 0. 05, hypergeometric distribution). The observed correlation of 0. 668 (p = 0. 035, Fig 5A) between the predicted and measured binding affinities can be considered relatively high, given the rather limited spectrum of kinases tested and the fact that instead of selecting the top predicted off-targets only, we focused on the set of kinases that were unique to tivozanib among 25 kinases with the highest predicted binding affinities against it. Our experimental results provide not only novel insights into the MoA of tivozanib, but also demonstrate how the in silico framework offers a cost-effective tool for prioritizing the most promising target interactions of an investigational compound for further experimental evaluation.
Recently, a lot of effort has been placed on the development of systems-based machine learning models that could aid drug discovery process in terms of providing cost-effective compound-target bioactivity predictions. Their main differences lie in the way how the models construct and treat molecular descriptors, and utilize various learning techniques, including those based on random forest [16,19,21], kernel learning [15,22,23,32], recommender systems [26], matrix factorization [20,25,26], Boltzmann machines [27], deep neural networks [16,17], logistic regression [28], learning to rank [29], and ensemble learning [30,31]. Although such models have been shown to perform well in cross-validation setups, their practical benefits still remain largely unknown due to the lack of systematic verification, using targeted experimental assays, carried out sub-sequent to the prediction phase. In model-guided mapping applications, the validation experiments are performed based on the model predictions. Such experimental validation setup effectively avoids any possible information leakage between the training and validation data, since the validation data does not exist at the time of making the predictions. The computational-experimental approach, implemented in this study, therefore makes it impossible to overfit the model to the training data. Here, we used the approach to evaluate the predictive power of a well-established kernel-based learning technique [33]. We chose this model family since kernel regression approaches have proved good performance in recent computational studies, including prediction of drug-target interactions [32], peptide-protein binding affinities [44], drug sensitivities in cancer cell lines [45], as well as in metabolite identification [46] and QSAR modelling [35]. Taking into account that molecular interactions are not simple on/off relationships, we focused on a binding affinity prediction problem, using the RLS regression model with a Kronecker product kernel (KronRLS). Our systematic evaluation of the predictive performance of various descriptors in the form of kernels revealed that their choice has a critical impact on the prediction accuracy. This is expected because kernel matrix is a central component of the kernel-based learning algorithm as it should capture our prior belief on the relationships between the input objects. In particular, known binding affinities, even if sparse, constitute an important information source not only for model training but also for the kernel matrix construction. Purely structure-based chemical descriptors were not able to fully capture the changes in compounds’ activity caused by minor structural differences. Furthermore, we introduced a novel protein kernel (KP-SW+), based on extended target profile, and showed how it consistently outperformed its commonly-used counterpart, in which the Smith-Waterman amino acid sequence alignment is adapted exclusively to proteins included in the dataset of interest (KP-SW). This was evident particularly under the New Target setup (Fig 2C), where the aim is to predict compounds targeting a new protein not encountered in the training data (correlation of 0. 669 for KP-SW+ vs. 0. 506 for KP-SW; S7 Fig). Under the same setting, we also observed a clear advantage of using, for the first time in the context of drug-protein interaction inference, generic string kernel applied to kinase domains and ATP-binding pockets over full protein sequences (correlations of 0. 651 for KP-GS-domain, 0. 628 for KP-GS-atp, 0. 508 for KP-GS; S7 Fig). The majority of kinase inhibitors, including those considered here, bind to ATP-binding pockets, and short sequences of these pockets are included within the kinase domain sequences, thus capturing also the neighbouring context. However, polypharmacological activities of kinase inhibitors, originating from the conservation of kinase ATP-binding pockets, make the prediction problem highly challenging, and better accuracies are likely obtained with compounds having more distinct target profiles. The methodology introduced here could equally well incorporate other compound and protein classes, such as ion channels or G-protein-coupled receptors (GPCRs), but further work will be required to investigate its practical performance under various scenarios using both computational and experimental validations. For instance, it remains an open question which kernels should be calculated to best represent such extended pharmacological spaces. Our results also demonstrate the importance of a proper evaluation procedure of the in silico models. A rigorous computational CV protocol is critical to ensure realistic performance estimates for the optimized models. In particular, the lack of the nested CV strategy in the model selection may lead to over-optimistic prediction results [32]. It is also important that CV design reflects the practical application use case of the model. Given a query drug-protein pair (dx, px), four different prediction scenarios can be distinguished, depending on whether there exist other compounds with measured bioactivities against px, or proteins with measured bioactivities against dx (Fig 2 and S1 Table). In turn, different types of CV designs need to be implemented in order to tune the model parameters and to evaluate its predictive performance. Here, we focused on the two most common and practical scenarios of the Bioactivity Imputation (Fig 2A) and New Drug (Fig 2B). Additionally, we provided the CV results under the New Target setup (Fig 2C) in S7 Fig. We first adopted LOO-CV, a design that simulates scattered missing values in otherwise known compound-target bioactivity map, to optimize the model and asses its performance in filling the experimental gaps (Bioactivity Imputation). Next, LDO-CV was applied in predicting target interactions for new candidate compounds having no measured bioactivity data for the model training (New Drug). The latter is much more challenging task, which was also demonstrated in our results; the correlation between the measured and predicted binding affinities under the Bioactivity Imputation setup was much higher (0. 829) than under the New Drug scenario (0. 653). However, even in the latter scenario, we still obtained a large number of statistically significant correlation values (S3 Fig). We further observed a high average AUC values under the Bioactivity Imputation setup (0. 945, S2 Fig), but also under the New Drug prediction scenario (0. 853, S2 Fig), which indicates that the model is able to discriminate well the interacting from non-interacting compound-kinase pairs. As expected, the classification accuracy increased with the increasing activity threshold as the true positive set includes a growing number of most likely interactions. In practical applications, the method requires features extracted from both compounds and proteins, such as readily available chemical two-dimensional structures and amino acid sequences, respectively, based on which kernels can be then calculated. However, if one is interested in a bioactivity of an investigational drug against a protein with unknown sequence, a reasonable prediction accuracy can still be achieved if there exist other compounds with measured bioactivities against a query protein (the New Drug setup). The practical solution is to replace the protein kernel matrix with an identity matrix, which implies that each protein will be considered similar to itself only and, effectively, the model will use just known bioactivity data and drug-drug similarities during both training and prediction phase. In particular, we noted only a small drop in Pearson correlation after replacing the protein kernel (KP-GS) with the identity matrix (0. 645 vs. 0. 653). This observation is expected in multitask or transfer learning problems, such as the New Drug or New Target setups, where one of the similarities is essential for generalizing to new instances (drug-drug similarities under the New Drug setting, protein-protein similarities under the New Target scenario). Ultimately, even though proper CV design is crucial to tune the model and assess its performance, subsequent experimental verification in the laboratory is the only way to really demonstrate the practical utility of the model predictions for drug discovery applications. The relatively good agreement between the computationally-predicted and experimentally-measured bioactivities validated the potential of the kernel-based algorithm, not only for filling the experimental gaps in existing drug-target interaction maps, but also in later stages of the drug development process, including prioritizing new target interactions of investigational compounds for further experimental evaluation, hence assisting in understanding of their MoA. Even though in silico inference of target interactions for new candidate drug compounds is a highly challenging task, our results with tivozanib suggest that, given enough-representative and high-quality training data, reliable off-target interaction predictions can be made. In addition to tivozanib, we initially considered also three other investigational kinase inhibitors, namely fedratinib, vx11e and ulixertinib (a compound derived from vx11e). However, we finally selected tivozanib because, unlike for the other compounds, its known on-targets were placed among the strongest predicted target interactions (S3 Table). This indicates that the primary on-target space of tivozanib is well-represented in our training data (S8 Fig). In the future, it is therefore important to profile and build up more diverse training data sets, including more examples of compounds targeting different kinase and other target classes. Tivozanib was originally-developed as a VEGFR inhibitor meant to block angiogenesis by targeting endothelial cells in the tumor vasculature. However, its MoA has not yet been fully elucidated. Based on the model predictions, our experimental assay confirmed the two Src family kinases FRK and FYN A, as well as the non-receptor tyrosine kinase ABL1 and serine/threonine kinase SLK as tivozanib’s off-targets. Our results highlight that tivozanib has an unusual target spectrum beyond the VEGFR family of kinases, and this suggests that the best anti-cancer use of this compound may not be in diseases where other VEGFR inhibitors with different target profiles have proven effective, but rather in ones where the target spectrum of tivozanib is more unique. For example, it can be hypothesized that tivozanib may have powerful activity in Src-family kinase addicted cancers, where it would target both angiogenesis and the cancer cells directly. Tivozanib has been shown to have a better safety profile than other marketed tyrosine kinase inhibitors and, currently, it is undergoing several clinical trials for the treatment of renal cell carcinoma (NCT03136627), refractory advanced renal cell carcinoma (NCT02627963), metastatic and non-resectable soft tissue sarcomas (NCT01782313), advanced liver cancer (NCT01835223), recurrent ovarian, fallopian tube, or primary peritoneal cancer (NCT01853644), and advanced prostate cancer (NCT01885949; July 2017). It will be interesting to see which of these trials will report successful treatment outcomes, and which will be terminated due to insufficient efficacy or toxicity. Although presented here results are promising, there is much room for improvement. For instance, we formulated the predictive model using only one pair of chemical and genomic descriptors at a time. However, even better accuracies could be obtained with a multiple kernel learning framework, which integrates multiple biological and molecular data sources, along with learning their importance for the prediction task [47]. Additional improvement of the predictive performance could be achieved also by creating more sophisticated kernelized molecular descriptors, for instance, by comparing three-dimensional structures of protein binding pockets. Furthermore, we used here as the training data a single yet very comprehensive kinase inhibitor profiling assay containing large number (16,265) of measured compound-kinase binding affinities, spanning different kinase branches (S11 Fig). However, as a future direction, we plan to work on integrating bioactivity values originating from various target profiling experiments and bioactivity end-points into a single model [3,19,40,48]. Recently-initiated community-driven efforts, such as Drug Target Commons (https: //drugtargetcommons. fimm. fi), which aim to collectively extract, manage and curate high-quality compound-target bioactivity data from public databases, literature and other resources, as well as annotate them with a common ontology, will be essential to facilitate the data standardization and computational modelling purposes. Nevertheless, we hope the current work provides a useful starting point and a practical guide on how to computationally prioritize the most promising target interactions for further experimental evaluation.
We used publicly available compound-target interaction map generated by Metz et al. using a large-scale functional bioassay, which measured the concentration of a compound needed to inhibit the reaction catalysed by a kinase enzyme of interest by 50% [3]. The readout corresponds to an inhibition constant Ki, typically expressed in the logarithmic scale as pKi = -log10Ki. Although the universal activity threshold cannot be explicitly defined for each compound-kinase pair, the higher the pKi value, the stronger the binding affinity between the compound and kinase. Among molecules included in the screen, 201 compounds are present in ChEMBL [43], and 169 proteins belong to the group of catalytically active human protein kinases [37]. The study is not complete, and therefore we used here a subset of these data: kinases and compounds for which at least 30% of the binding affinity values were measured, resulting in 152 drug compounds and 138 kinase targets. In total, there are 16,265 binding affinities in this selected interaction map (S9 and S10 Figs). Of note, most of the compounds constitute investigational, not yet FDA-approved, chemical probes. In supervised learning tasks, training data has the form { (xi, yi) }i=1N, where N denotes the number of training examples, xi ∈ X is an input object represented as a vector with the feature values (e. g. a compound represented as a fingerprint vector) and yi ∈ Y is its known associated label value (e. g. a potency of a compound against a certain protein). The aim is to find a prediction function f that models the relationship between xi’s and yi’s, and which can then be used to predict the label values for new instances outside the training space. Classical algorithms search for linear dependencies but often the actual relations underlying the data are highly nonlinear. Kernels offer the advantage of increasing the power of the linear learning machines by providing a computationally efficient way of projecting the input data into a high-dimensional feature space. A linear model in this implicit feature space corresponds to a nonlinear model in the original space. A separation of the statistical learning technique and the data representation is another convenient attribute of kernels. Formally, a kernel is a function k that for all x, z ∈ X satisfies k (x, z) = ⟨ϕ (x), ϕ (z) ⟩, where ϕ is a mapping from the input space X to an inner product high-dimensional feature space F: ϕ: x ∈ X → ϕ (x) ∈ F, and it can be considered as a similarity measure between two objects x and z. It is, however, often possible to avoid the explicit computation of the mapping ϕ, and define the kernel directly in terms of the original input data items by replacing the inner product ⟨∙, ∙⟩ with an appropriately chosen kernel function satisfying certain mathematical properties (so-called kernel trick). Kernels are particularly handy for calculating similarities between structured objects, including molecules. Here, we focused on a regression problem with the objective of predicting real-valued compound-target binding affinities. We used Kronecker regularized least-squares model (KronRLS) [33,49], a special variant of kernel ridge regression (KRR) which combines linear least squares with L2-norm regularization (ridge regression) and the kernel trick [50]. In KRR, given a set of N compound-protein pairs as training inputs xi’s, i = 1, …, N, and associated labels yi’s indicating binding affinities between them, we aim to find the minimizer of the following objective function J: J (f) =∑i=1N (f (xi) −yi) 2+λ‖f‖k2, (1) where f indicates the prediction function, f (xi) is the predicted binding affinity of ith compound-protein pair xi, λ denotes a regularisation parameter controlling the balance between training error and model complexity (λ > 0), and ‖f‖k is a norm of f on the space associated to kernel function k. In Eq (1), the first term corresponds to the training error, and the second, controlled by λ, is the penalty term that is larger for complex models that are more likely to overfit to training data but not generalize well to new instances. According to the representer theorem [51], the prediction function that minimizes J (f) can be expressed in terms of linear combination of the training examples: f (x) =∑i=1Nαik (xi, x) =αTk, (2) where k is a vector with kernel values k (xi, x) between each training point xi and test point x for which the prediction is made. The squared norm of f is therefore written as ‖f‖k2=∑i=1N∑j=1Nαiαjk (xi, xj). (3) A vector α, consisting of parameters αi that define the solution to KRR, is found by solving the following system of linear equations (K+λI) α=y, (4) where I is the N×N identity matrix, and y is the vector consisting of labels yi. K denotes N×N pairwise kernel matrix constructed for all training examples x1, x2, …, xN, and thus containing similarities between all compound-protein pairs. However, the size of K makes the training of the model computationally very heavy, even for moderate number of compounds and proteins. KronRLS is a special variant of KRR, where one assumes each data point xito consist of two separate parts, such as compound and protein, each equipped with its own kernel function, which enables to speed up the model training. Indeed, pairwise kernel K is computed as the Kronecker product of compound kernel KD of size nD×nD and protein kernel KP of size nP×nP (N = nD×nP): K=KD⨂KP. (5) Using Eq (5), the solution to KronRLS can be calculated from: α=vec (UPCUDT), (6) where vec (·) is the vectorization operator that arranges the columns of a matrix into a vector, UD and UP are orthogonal matrices with eigenvectors of drug kernel KD and protein kernel KP, respectively: KD=UDΣDUDT, (7) KP=UPΣPUPT, (8) and C is a matrix for which it holds that: vec (C) = (ΣD⨂ΣP+λI) −1vec (UPTYTUD). (9) Here, ΣD and ΣP denote diagonal matrices containing eigenvalues of KD and KP. Label matrix Y stores binding affinities between nD drug compounds (rows) and nP protein targets (columns). This way, we completely avoid the computation of the large pairwise kernel K, and therefore significantly shorten the training time. After applying the well-known property of the Kronecker product, (A⊗B) vec (D) = vec (BDAT), the prediction for test point x can be calculated as f (x) = (kD⨂kP) α= (kD⨂kP) vec (UPCUDT) =kP (UPCUDT) kDT. (10) Of note, the above shortcuts work only if there are no missing values present in the label matrix Y. Thus, as the pre-processing step, we imputed the missing binding affinities in Y by the weighted row (compound) average. The contribution of each protein was weighted by its similarity (normalized Smith Waterman score) to the protein for which the binding affinity was missing. Such imputed values were discarded when assessing the predictive performance of the model. The implementation of KronRLS is available at https: //github. com/aatapa/RLScore. We tuned the regularization parameter λ of KronRLS algorithm using nested CV (Table 1). We computed several types of drug compound and protein target molecular descriptors in the form of kernel matrices KD and KP, respectively. The summary is presented in Fig 3C. We used nested cross-validation (CV) procedure for model selection, and we assessed the predictive power with Pearson correlation and root mean squared error (RMSE) between original and predicted pKi values. In k-fold CV, the dataset is randomly divided into k subsamples of equal size, and the model is trained based on k-1 of them (training data). Then, the remaining subsample (test data) is used to assess how well the model that has been found generalizes to new instances, i. e. to calculate the predictive performance. The procedure is repeated k times, such that each subsample is used once as the test data, and the average error over the k folds gives the final estimate. Nested CV consists of two loops, outer and inner one. In the outer CV loop, each of the k folds is kept as a test set at a time. The remaining k-1 training folds of the outer CV loop are further divided into training and test set of the inner CV loop. Here, the inner CV was performed during each round of the outer CV, with the aim of selecting the regularization parameter λ of KronRLS algorithm as well as kernel parameters (S4 Table). We performed a grid search in order to select the most suitable combination of all parameters. Then, training folds of the outer CV loop were used to train the model with selected parameters, and the predictive performance was evaluated on the test set. We applied two different CV strategies, i. e. leave-one-out nested cross-validation (LOO-CV) and leave-drug-out nested cross-validation (LDO-CV), summarized in Table 1 and S12 and S13 Figs. We note that it is critical to discard binding affinities of compound-protein pairs belonging to the test fold of both inner and outer CV prior to computing Gaussian interaction profile kernels (KD-GIP, KP-GIP) in order to avoid significant model overfitting. Here, we used interaction profile kernels with LOO-CV, and we removed from the compound-protein interaction matrix Y the whole column (row) containing the test point before computing KD-GIP (KP-GIP) kernel.
Significant efforts have been devoted in recent years to the development of machine learning models to support different stages of drug development process. Given the enormous size of the chemical universe, such models could offer a complementary and cost-effective means to experimental determination of drug-target interactions, toward prioritization of the most potent ones for further verification in the laboratory. In order to demonstrate the benefits of the prediction models in practical application cases, we carefully evaluated the predictive power of a well-established machine learning model in filling the gaps in existing profiling studies and prediction of target interactions for a new drug candidate. As a specific case study, we focused on kinase inhibitors, which form the largest class of new drugs approved for cancer treatment, but are also known to have wide multi-target activities contributing both to their therapeutic and toxic responses. The high agreement observed between the predicted and experimentally-measured drug-target bioactivities under the implemented rigorous validation setup demonstrates the potential of the machine learning approach, not only for filling the gaps in existing drug-target interaction maps, but also toward off-target interaction prediction for investigational drugs, and finding potential new uses for already approved drugs (drug repurposing).
sequencing techniques
medicine and health sciences
protein interactions
applied mathematics
enzymology
simulation and modeling
algorithms
mathematics
forecasting
statistics (mathematics)
artificial intelligence
pharmacology
molecular biology techniques
enzyme inhibitors
research and analysis methods
sequence analysis
computer and information sciences
sequence alignment
bioinformatics
proteins
mathematical and statistical techniques
statistical methods
molecular biology
biochemistry
kinase inhibitors
kernel methods
database and informatics methods
biology and life sciences
physical sciences
drug interactions
protein sequencing
machine learning
Computational-experimental approach to drug-target interaction mapping: A case study on kinase inhibitors