Developing a machine learning model to identify protein–protein interaction hotspots to facilitate drug discovery

Dataset collection and feature aggregation

As a starting point, the dataset and codebase from the SpotOn study (Moreira et al., 2017) were acquired. This study was selected as the starting point for its high effectiveness in identifying potential hotspots that could aid in drug discovery. The SpotOn database already has information regarding amino acid composition, solvent-accessible surface area (SASA) information, position-specific scoring matrices, the number of amino acids at 2.5 and 4.0 Angstrom, the number of nearby hydrophobic residues, the total change in SASA, the number of interfacial residues, pseudo-amino acid composition, and scales-based descriptors of 2D and 3D descriptors from the protr R package (see below) for a total of 881 features.

In order to add more information to this dataset to better aid model prediction, the protr R package (Xiao et al., 2015) was used to add more features related to amino acid composition, dipeptide composition, etc., to the already pre-existing data. Additionally, data related to pair potential, complex/monomer accessible surface area, residue information, amino acid information, etc. were extracted from the HotPoint database (Tuncbag, Keskin & Gursoy, 2010) and then added to the pre-existing dataset. This data was added to add more information regarding the entire protein complex, as evidenced by most of protr’s features, and to add residue specific features such as pair potential that could improve predictive power. The addition of new features in the protr R package and the HotPoint database led to a total of 2,323 features.

Upon further investigation of the SpotOn dataset, we found that chains I of proteins with PDB code 2FTL, 3SG8, and 1CH0 do not exist as specified in the Protein Data Bank. In the SpotOn study, these chains are specified, and features were derived for these chains; however, in this study, as additional features are added and these chains could not be identified, these chains have been removed from our dataset. This leads to a total of 520 protein residues, lower than SpotOn’s 534 protein residues. A total of 398 residues are labeled as non-hotspots, and 122 residues are labeled as hotspots.

In order to derive features on our prediction dataset with the EphB2-ephrinB2 complex (PDB code: 1KGY), we first downloaded the structure from the Protein Data Bank, and ran this structure through the SpotOn’s codebase/pipeline to collect features specific to the SpotOn study. Then, we sequentially added additional features unique to this study, such as from the protr’s R package and features from the HotPoint database. Missing values are assigned the average value of all non-missing values in a feature.

Preprocessing and feature engineering

Similar to the SpotOn study, both the training and testing sets were normalized, and the testing set was normalized using mean and standard deviation of the training set. In addition, before the model was run, data balance had to be accounted for, and oversampling was performed in order to retain the properties of the majority class without sacrificing the information available in this class (More, 2016). SMOTE, or synthetic minority oversampling technique, was performed with k = 5 nearest neighbors (Chawla et al., 2002). To account for multicollinearity, principal component analysis was also performed. This leads to four different combinations: a pipeline without any changes to the training data, a pipeline with only SMOTE applied, a pipeline with only PCA applied, and a pipeline with both SMOTE and PCA applied.

Before the model was trained, the dataset was first subjected to feature engineering. Three existing features that were selected for further exploration are the number of intermolecular contacts within 4.0 Angstroms (#Dist-4.0), the number of hydrophobic contacts (#Hydrophobic), and the pair potential of a specific residue (Pair Potential). We hypothesized that an increase of hydrophobic contacts would cause a decrease in hydrophobic pair potential due to the attractive interaction because of the hydrophobic effect (Israelachvili & Pashley, 1982). As a result, we multiplied both variables and multiplied by −1 to amplify the effects of this association and accounting for the inverse correlation. In addition, we hypothesized that the number of intermolecular contacts will increase the pair potential as this may lead to many body potentials, which are mostly repulsive at short distances (Byggmästar, Granberg & Nordlund, 2018). To model this association, #Dist-4.0 and #Hydrophobic are multiplied to amplify the effects as well. These two new engineered variables were named #Dist-4.0 * Pair Potential and -#Hydrophobic * Pair Potential. This lead to a total of 2,323 features on the training and testing datasets, as well as our dataset containing residue information on the crystal structure of the EphB2-ephrinB2 complex (PDB code: 1KGY).

Machine learning model selection

Five different machine learning models were selected in order to evaluate and develop a model: Logistic regression (LR), XGBoost (XGB), a balanced random forest classifier (RF), K Nearest Neighbors (KNN), multilayer perceptron neural network (MLP), and a Gaussian Naïve Bayes (GNB). This data was then split into a training:testing set ratio of 80:20. 10-fold cross validation was performed 10 times on the training set to prevent overfitting. GridSearch was performed in order to identify the best combination of hyperparameters/parameters that could yield the best results. The following hyperparameters/parameters were tested: LR, with C equal to 0.01, 0.1, 1, 10, 50, 100, 500, 1,000, and 5,000; RF with balanced class weight, with the number of estimators equal to 50, 100, 150, 250, 350, 500, and maximum depth of 5, 7, 9, 11; XGB, with a learning rate of 0.001, 0.01, 0.1, the number of estimators as 50, 100, 150, 200, and maximum depth of 4, 5, 6; KNN, with n neighbors of 1, 3, 5, 10, 15, 20; a multilayer perceptron model of hidden_layer_sizes (10, 10, 10), (50, 1), (10, 10), (10, 1), (5, 5), (5, 5, 5) and alpha of 0.0001, 0.0002, 0.0005, 0.001; and GNB with variance smoothing of 1e−8, 1e−7, 1e−6, 1e−5, 1e−4, 1e−3, and 1e−2. The metric used to identify the best model from these sets of parameters on the validation set is AUROC (area under receiver operating characteristic), and this was chosen to compare our models more accurately with the AUROC and receiver-operating characteristic curves provided in SpotOn. Four different run conditions on the four different pipelines was also run and the results are compared. The model(s) with the best run conditions on the highest scoring pre-processing dataset will be used to build an ensemble model, similar to the SpotOn study. If the ensemble model has a higher predictive capability than any individual model, the ensemble model will then be used to predict hotspots on the EphB2-ephrinB2 complex, as this complex has been overexpressed in many cancer cells, most notably in prostate, gastric, colorectal and melanoma cancers (Pasquale, 2010). PyMol was utilized to visualize the hotspots predicted on the EphB2-ephrinB2 complex.

Small molecule selection

A cluster of hotspots was identified and LigandScout (Wolber & Langer, 2005) was used to create an apo-site pharmacophore. Virtual screening was then performed on this pharmacophore to identify possible new drug indications. To perform the drug screening, an approved Drugbank (Wishart et al., 2008) database that has a library of all molecules that have molecular weight from 150 to 500 daltons was used. These small molecules were then ranked by the LigandScout software to identify molecules that most strongly conform to the pharmacophore based on the chemical and structural properties of that molecule. The drug-disease associations were then verified with scientific literature to assess the validity and efficacy of the model, and then we identified new drug-disease associations that have not been previously identified by cross-referencing existing scientific literature.

Metric calculation

In context, sensitivity is the ability for the model to identify the hotspots and the specificity/recall is the ability for the model to identify the non-hotspots, and both of these statistics are defined as: $$\mathrm{R}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}=\mathrm{S}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}={\displaystyle \frac{\mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e}\mathrm{P}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}}{\mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e}\mathrm{P}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}+\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e}\mathrm{N}\mathrm{e}\mathrm{g}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}}}$$ $$\mathrm{S}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{f}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{t}\mathrm{y}={\displaystyle \frac{\mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e}\mathrm{N}\mathrm{e}\mathrm{g}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}}{\mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e}\mathrm{N}\mathrm{e}\mathrm{g}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}+\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e}\mathrm{P}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}}}$$

Precision is defined as: $$\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}={\displaystyle \frac{\mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e}\mathrm{P}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}}{\mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e}\mathrm{P}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}+\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e}\mathrm{P}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}}}$$

F1, MCC, Kappa, and Precision-Recall are all metrics that are robust in dealing with data imbalance. They are defined as: $$\mathrm{F}1=2\ast {\displaystyle \frac{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\ast \mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}}{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}+\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}}}$$ $$\mathrm{M}\mathrm{C}\mathrm{C}={\displaystyle \frac{\mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e}\mathrm{P}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\ast \mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e}\mathrm{N}\mathrm{e}\mathrm{g}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}-\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e}\mathrm{P}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\ast \mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e}\mathrm{N}\mathrm{e}\mathrm{g}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}}{\sqrt{\left(\mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e}\mathrm{P}\mathrm{o}\mathrm{s}+\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e}\mathrm{P}\mathrm{o}\mathrm{s}\right)\left(\mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e}\mathrm{P}\mathrm{o}\mathrm{s}+\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e}\mathrm{N}\mathrm{e}\mathrm{g}\right)\left(\mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e}\mathrm{N}\mathrm{e}\mathrm{g}+\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e}\mathrm{P}\mathrm{o}\mathrm{s}\right)\left(\mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e}\mathrm{N}\mathrm{e}\mathrm{g}+\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e}\mathrm{N}\mathrm{e}\mathrm{g}\right)}}}$$

$Kappa=\left(p_o-p_e\right)/\left(1-p_e\right)$ where p_o is the probability of agreementassigned to any sample, and p_e is the expected/hypothetical probability of chance agreement.

Precision-Recall $=\sum _n\left(R_n-R_{n-1}\right)P_n$ where P_n and R_n are precision and recall, respectively, at the n^th threshold.

All of these calculations are calculated using the Scikit-learn package in Python.

Phase 1

The average test metrics of each of the six algorithms tested on the four different pre-processing pipelines are shown in Table 1. As the preprocessing pipeline where SMOTE and PCA are applied has the highest AUROC, the top algorithms from this pipeline are used to create an ensemble model.

In Table 2 are the best individual algorithms tested in the SMOTE and PCA pipeline. The best set of hyperparameters were selected using GridSearch as follows: the logistic regression with C = 1,000, the random forest classifier with maximum depth of 5 trees and the total number of estimators at 50 trees, an XGBoost classifier with learning rate 0.1, maximum depth of 5, and 50 estimators, K-nearest neighbors with 5 neighbors, a multi-layer perceptron classifier with alpha as 0.0002 and three layers of 10 neurons each, and a Gaussian Naïve Bayes of variable smoothing of 1e−5. The best performing algorithm in this pipeline is the logistic regression, as it has the highest AUROC.

The results of the Logistic Regression from the SMOTE and PCA pipeline were compared with SpotOn’s highest performing algorithm from the upsampling pre-processing procedure, as shown in Table 3. In order to more aptly analyze the predictive capabilities of our top performing algorithm, we adjusted the threshold to achieve a ≥0.88 specificity, as that is the specificity of the highest performing individual algorithm in SpotOn, and comparing other metrics based on this threshold change. This reduced the AUROC from 0.842, as identified in Table 2, to 0.840, as identified in Table 3, but was still greater than any individual model in SpotOn. Our algorithms perform better to that of SpotOn’s ScaledUp processing step in all six metrics as demonstrated in Table 3.

The top ranking algorithms in the SMOTE and PCA pipeline are used to develop an ensemble classifier to achieve better performance compared to any single algorithm. Different ensemble algorithms are tested: stacking, where a meta-classifier is used to combine the predictive power multiple base classifiers, and voting, a simple ensemble method where each of the six algorithms tested votes on a specific data point, and a simple majority vote is used to predict the classification of that data point as shown in Table 4. In this case, the meta-classifier used during stacking is a Logistic Regression classifier where C = 5. Each individual model is used as a base model separately with the meta-classifier, and all models are combined with the meta-classifier. All ensemble models are run on the SMOTE and PCA pipeline. In the voting ensemble, hard voting was implemented, and all six algorithms are subjected to majority voting. Here, the best performing classifier was the stacking classifier where all models are combined with the meta-classifier. However, the AUROC of this ensemble method is still lower than that of the top individual model, the logistic regression in the SMOTE and PCA pipeline.

A comparison of the accuracy and performance of the model developed herein, shown in bold, compared with SpotOn as shown in Table 5. In our model, the logistic regression was our top performing algorithm, and was thus used to develop to predict hotspots with high accuracies. The SpotOn study (Moreira et al., 2017) was used in order to identify the testing accuracies of the SpotOn study and those of the other studies as well. The other studies that are compared to are SpotOn, SBHD213, Robetta23, KFC2-A24, KFC2-B, and CPORT25 (Kim, Chivian & Baker, 2004; Martins et al., 2014; De Vries & Bonvin, 2011; Zhu & Mitchell, 2011).

Phase 2

The top logistic regression algorithm from the SMOTE and PCA pipeline was utilized to predict the hotspots from the EphB2-ephrinB2 complex. This resulted in consecutive residues 1122–1126 predicted as viable hotspots. Drug screening performed on these hotspots results in nine potential small molecules that could bind to this hotspot. The pharmacophore fit scores are provided from our drug screening analysis, where the highest ranked drug also has the highest pharmacophore fit score as shown in Table 6. Extensive literature review was conducted to identify which drugs could aid in treating conditions associated with the EphB2-ephrinB2 complex.

Phase 1

As the best performing classifier, the logistic regression from the SMOTE and PCA pipeline is subjected to PCA, it is difficult to analyze exactly which features are the most valuable. As a result, in order to identify the most relevant features, sensitivity analysis using Python library Pytolemaic (https://pypi.org/project/pytolemaic/) was performed on the best performing logistic regression classifier from the SMOTE-only pipeline to understand the significance of adding new features to the existing dataset as provided by the SpotOn study. Three out of the top ten features (Relative Complex ASA, Complex ASA, and the engineered features Dist-4.0 * Pair Potential) as identified by Fig. 1 were added in this study exclusively, and highlights the improvement in predictive capabilities of the addition of these features.

PyMol (Schrödinger, 2015) was used to derive and highlight residues 1122–1126 on chain E of the EphB2-ephrinB2 complex. Predicted druggable hotspot residues are shown as more visible surface markers (in green), and the other residues are shown in pink or light red in Fig. 2. Residues 1122–1126 were selected for drug screening as consecutive residues may be used as initial fragments in drug screening (Modell, Blosser & Arora, 2016).

Phase 2

An apo-site grid, as shown in Fig. 3, was developed and implemented on hotspot residues 1122–1126 as identified via the machine learning model on the EphB2-ephrinB2 complex, where a visual representation of the hotspots on the complex are shown in Fig. 2. This grid was developed by first calculating the pockets of hotspot residues 1122–1126 on LigandScout (Wolber & Langer, 2005).

Figure 4 shows the 26-feature pharmacophore developed using an apo-site grid derived using hotspot residues 1122, 1123, 1124, 1125, and 1126 identified in Fig. 3 via the machine learning model. A pharmacophore identifies the key parts of the molecular features that define the function and shape of a specific ligand, and includes features such as H-bond acceptors and donors, hydrophobic and aromatic rings, etc. This pharmacophore is then used to identify drugs that fit its features. The scoring of this screening procedure follows a pharmacophore-fit scoring function as provided in LigandScout. A maximum number of two features are omitted from this multi-feature pharmacophore to identify small molecule hits, and the best matching conformation is selected.

Cimetidine, currently an acid reflux medication, was identified via virtual screening to potentially bind to the EphB2-ephrinB2 complex associated with cancer cells as shown in Fig. 5. The right image is cimetidine in relation to the 26-feature pharmacophore developed as shown in Fig. 4. A pharmacophore-fit score of 43.86 was achieved during drug screening. Further literature review identified cimetidine as a potential repositioning target for many different types of cancers, including melanoma, gastric, and colorectal cancers (Pantziarka et al., 2014).

Idarubicin, a chemotherapy medication that’s currently used to treat breast cancer, was identified via virtual screening to potentially bind to the EphB2-ephrinB2 complex, where the expression of the complex is associated with cancer cells as shown in Fig. 6. The pharmacophore fit score of this small molecule is 45.46. This drug was also found to treat cancers linked to the EphB2-ephrinB2 complex such as melanoma and leukemia (Martoni et al., 1986; Jabbour et al., 2017). The right image is idarubicin in relation to the pharmacophore developed as shown in Fig. 4.

Pralatrexate, a T-cell lymphoma medication, was identified via virtual screening to potentially bind to the EphB2-ephrinB2 complex, where the expression of the complex is associated with cancer cells as shown in Fig. 7. This small molecule has a pharmacophore fit score of 47.41, and literature review suggests that this drug could potentially treat breast cancer and prostate cancer (Yu, Zhao & Gao, 2018; Serova et al., 2011). Recent research has shown that pralatrexate can treat esophagogastric cancer, which is associated with the various gastrointestinal cancers of the EphB2-ephrinB2 complex (Malhotra et al., 2020). The right image is pralatrexate in relation to the pharmacophore developed as shown in Fig. 4.

Nadolol, a beta blocker, was identified via virtual screening to potentially bind to the EphB2-ephrinB2 complex, where the expression of the complex is associated with cancer cells as shown in Fig. 8. This small molecule has a pharmacophore fit score of 45.97, and literature review suggests that beta blockers could potentially treat a variety of cancers, including breast cancer and pancreatic cancer (Ishida et al., 2016). A close relative of this drug, propranolol, can induce apoptosis in liver cancer cells (Wang et al., 2018). This research suggests nadolol’s potential role in mitigating the effects of other cancers as well. The right image is nadolol in relation to the pharmacophore developed as shown in Fig. 4.

Virtual drug screening identified nine drugs (pralatrexate, chlortetracycline, nadolol, imipenem, idarubicin, valganciclovir, conivaptan, cimetidine, and barnidipine) that bind to the pharmacophore shown in Fig. 4. Further analysis via literature review identified four drug candidates to potentially treat various types of cancers: cimetidine, idarubicin, pralatrexate, and nadolol. Figure 5 shows the possibility for cimetidine, an antacid, to bind with the EphB2-ephrinB2 complex, and scientific literature identified the possibility for this drug to potentially treat melanoma, gastric, and colorectal cancers (Pantziarka et al., 2014). Figure 6 identifies the possibility for idarubicin, a chemotherapy drug used to treat leukemia, to bind with the EphB2-ephrinB2 complex, and literature review identified the possibility for this drug to potentially treat melanoma and leukemia (Martoni et al., 1986; Jabbour et al., 2017). Figure 7 demonstrates the possibility for pralatrexate, a T-cell lymphoma medication to bind to the EphB2-ephrinB2 complex.