All chemical compounds with activity against PIM-1 were collected from the literature and the ChEMBL database (Gaulton et al., 2012). Inorganic and duplicate compounds were removed from the list. Generally, compounds with IC₅₀ ≤ 10 μM will likely be “active,” predicting a large number of active molecules. However, such a high fraction of active compounds cannot be expected from any experimental platform. Therefore, in order to make the most efficient use of costly experimental validation, the optimal model should identify compounds with affinity higher than 10 μM. The higher the value, the higher the drug dose needed to achieve the required potency and, thus, the higher the chance of “off-target” activity. To address this issue, we chose to set the decision boundary at IC₅₀ ≤ 1 μM for active molecules. Molecular descriptors were calculated using the PaDEL software (Yap, 2011). A two-tier selection procedure was applied to select the best descriptors. First, we randomly selected one descriptor from a pair showing >0.85 correlation. Second, descriptors were reduced using the Boruta method (Kursa et al., 2010). We used four different machine learning methods, namely, Support Vector Machine (SVM) (Mitchell, 1997), random forest (Breiman, 2001), Extreme Gradient Boosting (XGBoost) (Chen and Guestrin, 2016), and kappa nearest neighbor (kNN) (Voulgaris and Magoulas, 2008), to build the classification models. All the classification experiments and calculations were conducted using the R.3.0.2 environment (http://www.R-project.org/) and Python (http://www.python.org/) platform. The compounds used in training and test sets are given in Supplementary Tables S1 and S2, respectively.

A receiver operating characteristic (ROC) plot and area under the curve (AUC) were used to assess the performance of the model (Hanley and McNeil, 1983; Park et al., 2004). In Table 1, the terms precision (Eq. 1), recall (Eq. 2), accuracy (Eq. 3), and F1 score (Eq. 4) are defined along with their relationships to the statistical performance calculations used to assess the quality of the model. P r e c i s i o n = T r u e   p o s i t i v e T r u e   p o s i t i v e + F a l s e   N e g a t i v e , (1) R e c a l l = T r u e   p o s i t i v e T r u e   p o s i t i v e + F a l s e   N e g a t i v e , (2) A c c u r a c y = T P + T N T P + T N + F P + F N , (3) F 1 = 1 . P r e c i s i o n   X   R e c a l l P r e c i s i o n + R e c a l l . (4)

In order to highlight the region of the chemical space that contains the chemicals for which the model is expected to make accurate predictions, a well-validated predictive model needs to have a defined applicability domain (AD) (Rakhimbekova et al., 2020). Any predictive model must verify its constraints in terms of its structural domain and response space. As a result, determining a model’s AD and evaluating the accuracy of its predictions are both challenging tasks. These QSAR models typically use the training set to cover a certain chemical space. The model’s predictions are accurate if any query compound falls within this definition of AD. If not, the prediction might not conform to the model’s presumptions. Principal component analysis (PCA) (Sushko et al., 2010) has been employed in our work to define the AD of the compounds used in this study.

To test the robustness of the proposed models, y-randomization was applied. This technique involves randomly mixing up the values of the target variable in the training set (Rücker et al., 2007; Lipiński and Szurmak, 2017). The same parameters used in the initial model are then applied to a new prediction generated with the scrambled data. Every estimate of the model’s accuracy was recorded. In total, 50% of the compounds in the training set were resampled and used in a 500-run y-randomization test.

The Tanimoto coefficient (Tc) (Eq. 5) was computed using MACCS-166 fingerprints to quantify chemical similarity. The active and inactive chemicals in the training set were compared against false and true positive compounds in systematic pairwise similarity computations. T c = C A + B − C . (5)

Molecular substructures related to PIM activity were analyzed using the distribution of MACSS fingerprints in active and inactive compounds (Eq. 6). F r e q u e n c y = ∑ i N F P 1 | 0 N X 100 . (6)

Additionally, the probability scores of the developed classification models were examined. In general, a molecule is defined as inactive if its probability score is lower than 0.5, while a compound with a probability score of 0.5 is considered active (Ponzoni et al., 2019). The more this score approaches 1, the more confident we are in our prediction. Here, we examined the probability score distributions for TP (true positive), TN (true negative), FP (false positive), and FN (false negative) results.

The developed models were used to screen the hits against PIM-1. The NCI library and Maybridge databases were used for virtual screening. The National Cancer Institute maintains a repository of compounds that have been evaluated as potential anticancer agents. These compounds represent unique structural diversity based on synthetic and natural products. The Maybridge library consists of a highly diverse set of over 53,000 lead-like compounds. Maybridge Hit-to-Lead was designed for medicinal chemistry, allowing SAR development and hit-to-lead optimization. The following filters were used to select the hits: Filter 1: compounds predicted to be active by all the validated models; Filter 2: compounds having a probability score; and Filter 3: compounds falling within the chemical space of the training set. These compounds were further processed for molecular docking, followed by molecular dynamics simulations. Finally, compounds with the best affinity and conformance within the active site were selected and analyzed.

Molecular docking was implemented to identify the best physical confirmation of inhibitor binding within the active site of PIM-1 kinase. The PIM kinase enzyme structure was taken from the Protein Data Bank (PDB ID: 5KZI). All of the docking simulations for this work were performed using AutoDock Vina (Trott and OlsonAutoDock, 2009) with a 1 spacing, default exhaustiveness, and full ligand flexibility. The grid resolution was internally set to 1Å. We set the number of binding modes to 10 and exhaustiveness to 8. A cubical grid of size 60 × 60 × 60 size with 0.375 Å spacing was used around the active sites of the protein. To acquire the structure in the PDBQT format, polar hydrogen atoms were added using AutoDock Tools 92.

Selected best compounds were further subjected to molecular dynamics (MD) simulations using Groningen Machine for Chemical Simulations (GROMACS v5.1.5) (Pronk et al., 2013). The parameters and coordinate files for PIM-1 kinase and selected potential hit compounds were generated using the CHARMM27 forcefield in GROMACS and PRODRG, respectively. The TIP3P water model was used for each simulation system, which was neutralized by the addition of Na⁺ ions in a dodecahedron periodic box. Energy minimization was performed for 50,000 nstep using the steepest descent algorithm to avoid steric clashes. Equilibration of each system was performed in two stages: the first phase was carried out with a constant number of particles, volume, and temperature (NVT) ensemble for 500 ps at 300 K, using the V-rescale thermostat (Bussi et al., 2007); and in the second phase, the pressure of each system was equilibrated for 500 ps at a constant number of particles, pressure, and temperature (NPT) at 1 bar using a Parrinello–Rahman barostat (Parrinello and Rahman, 1981). Each equilibrated system was simulated for 30 ns under periodic boundary conditions to avoid edge effects. Electrostatic interactions were handled by the particle mesh Ewald (PME) method, while the heavy-atom bonds were restrained using the LINCS algorithm.

In total, 54 descriptors from the set of 240 were eventually selected using the Boruta method (Supplementary Table S3). All these descriptors belonged to 12 different classes. The descriptors include autocorrelation, information content, atom-type electrotopological state, Burden modified eigenvalues, molecular distance edge, carbon type, and molecular linear free energy relation. The models were trained using four machine learning methods (SVM, random forest, XGBoost, and kNN). Evaluation metrics for the developed models are given in Table 1, including accuracy, recall, precision, F1 Score (a measure of a model’s accuracy, which takes into account both precision and recall), and Area Under the Curve (AUC) values. SVM, random forest, and XGBoost performed than kNN according to these metrics in combination with the selected descriptor set. Among the three, random forest achieved the best accuracy, at 0.87 for the test set (with selected descriptors), as compared to SVM (0.78) and XGBoost (0.84). In addition, these models also had significant AUC values (Figure 1).

An applicability domain (AD) analysis was performed to check the reliability of the generated classification models. Figure 2 shows a scatter plot of the PC1 and PC2 coordinates derived from the set of selected PIM-1 compound descriptors. The training and test compounds share similar PC1 and PC2 coordinates, suggesting that predictions were within the applicability domain (AD) of both the training and test sets. To check the robustness of the developed models, y-randomization tests were performed (Rücker et al., 2007). Y-randomization test accuracies were found to be lower, and none of the random trials achieved higher scores than our main models (Figure 3). The average accuracy across all randomly generated models were found to be less than 0.58. This confirms that the selected models are robust and reliable and were not generated by chance correlations. A pairwise comparison of the compounds in each cluster was found to reflect reasonable Tanimoto coefficient similarities between them.

Probability scores of the selected models, reflecting the probability of belonging to each class, were also analyzed. It is known that a compound with a probability score of ≥0.5 is classified as active, whereas a molecule with a probability below <0.5 is classified as inactive. As this score approaches 1, the higher the value, the higher the model’s confidence in the prediction is (Minerali et al., 2020; Esposito et al., 2021). In our study, we analyzed the distribution of probability scores among TN (true negative), FP (false positive), TP (true positive), and FN (false negative) results. For the SVM model, compounds with a probability score of more than 0.80 (an average value) were more likely to be active, whereas compounds with a probability score of 0.36 were more likely to be inactive. In the case of the random forest model, a compound with a probability score of more than 0.87 was more likely to be active, whereas a compound with a probability score of 0.24 was more likely to be inactive. Random forest achieved values of 0.95 and 0.11 for active and inactive compounds, respectively, indicating greater success in predicting compound activity with the desired probability score (Supplementary Figure S1). False positive compounds were predicted with probability scores of 0.63, 0.65, and 0.69 for the random forest, XGBoost, and SVM models, respectively. In contrast, false negative compounds were found to have probability scores of 0.31, 0.42, and 0.14 for the random forest, SVM, and XGBoost models, respectively. Each predictive model’s effectiveness in the early recognition of hits was visually evaluated using a cumulative gain plot (Table 2). The cumulative gain curve is an evaluation curve that evaluates the model’s performance and contrasts the outcomes with a random selection. It displays the percentage of targets identified when taking into account a particular portion of the population that has the highest likelihood of being a target based on the model. The comparison showed that the XGBoost and random forest methods performed better than SVM and kNN in terms of early recognition of hits (Figure 4).

Molecular substructures related to the PIM-1 activity of the compounds can be identified by analyzing the bits in the MACCS fingerprints. We analyzed the MACCS fingerprints showing a reasonable difference between active and inactive compounds (Supplementary Table S4). The occurrence of MACCS fingerprints differed significantly between active and inactive compounds in the training dataset, suggesting that the substructures represented by these features may be closely related to PIM-1 activity. Descriptions and the number of occurrences of these substructures are listed in Supplementary Table S4. It was found that MACCS38, MACCS52, MACCS92, MACCS98, MACCS107, MACCSFP142, etc. are prevalent in active molecules. This is consistent with previous studies, which shows that compounds with such functional groups have therapeutic potential against PIM kinase (Tsuganezawa et al., 2012; El-Hawary et al., 2018; Park et al., 2021).

The NCI and Maybridge databases were used to screen the potential hits from validated models. Commonly predicted active compounds with high probability scores were selected and further filtered out within the applicability domain (AD) of the training set. These compounds were further subjected to molecular docking simulation (Table 2). Finally, four compounds (CHEMBL303779, CHEMBL690270, CHEMBL748285, and N-[(1-ethylbenzimidazol-2-yl)methyl]-3-(4-methoxyphenyl)-1H-pyrazole-4-carboxamide (EBM-MPC)) were observed to have reasonable binding affinity and stable interaction with the catalytic residues in the active site (Table.3 and Figure 5). A literature survey revealed that Leu44, Lys67, Glu121, and Asp186 are crucial for the interaction of inhibitors (Tsuganezawa et al., 2012; El-Hawary et al., 2018; Park et al., 2021). It can be observed in Figure 4 that CHEMBL690270, CHEMBL303779, and EBM-MPC form hydrogen bond interactions with Lys67 and hydrophobic interactions with Asp186 (Figure 6). In contrast, CHEMBL748285 forms hydrogen bonds with Asp186 (Figure 6). The quinazoline ring of compounds was involved in multiple p–alkyl interactions. In addition, a number of hydrophobic contacts, particularly residues Leu44, Gly47, Phe49, Ile104, and Leu120, stabilize interaction with hits. PIM inhibitors fall into two broad categories: ATP mimetics, which form hydrogen bonds with the glutamate residue that serves as the hinge (Glu121), and non-ATP mimetics, which bind far from the hinge or interact with the hinge through hydrophobic interactions with a number of residues in the specific hydrophobic pocket that serves as the hinge environment (El-Hawary et al., 2018; Park et al., 2021). The Tanimoto coefficient (Tc) similarity score of these selected hits was found to be ≤ 0.5 with high-activity compounds (Figure 5B).

By analyzing 100-ns MD trajectories, the structural changes to PIM-1 upon inhibitor binding were studied. We examined the RMSD of the protein backbone and the RMSF of the protein’s alpha-carbon atoms. As shown in Figure 6, all the systems exhibited stability throughout the 100-ns simulation. The average RMSD value for all four systems was observed to be below 0.31 nm, which indicated that simulated complexes displayed RMSD values below the threshold. The average RMSD values further showed that the CHEMBL690270 PIM-1 complex displayed less deviation (0.26 nm), whereas CHEMBL303779 and CHEMBL748285 demonstrated similar average values of 0.34 nm (Figure 7A). RMSF is a significant value, used to characterize each residue’s fluctuation rate upon ligand binding. It was observed that the inhibitor binding residues (Leu44, Phe49, Lys67, Glu121, and Asp186) did not fluctuate significantly (Figure 7B).