A pragmatic pharmacophore informatics strategy to discover new potent inhibitors against pim-3

Pim-3 (proviral integration site moloney murine leukemia virus-3) is an oncogene which encodes proteins belonging to serine/threonine kinase family and PIM subfamily. It is generally over-expressed in epithelial and hematological tumors. It is known to involve in numerous cellular functions such as cell growth, differentiation, survival, tumorigenesis, and apoptosis. It also plays a crucial role in regulation of signal transduction cascades. Therefore, it emerged as a hopeful therapeutic target for cancer treatment. In the current study, indole derivatives having potent inhibitory activity against Pim-3 were taken and pharmacophore-based virtual screening was carried out. A 5-point pharmacophore hypothesis with one hydrogen bond acceptor, one hydrogen bond donor, and three aromatic rings, i.e., ADRRR, was developed with acceptable R2and Q2 values of 0.913 and 0.748, respectively. It was employed as a query and screening was conducted against Asinex and Otava lead library databases to screen out potent drug like candidates. The obtained compounds were subjected to SP, XP docking using 3D model of pim-3 which was constructed through comparative homology modelling, and finally binding free energies were calculated for top hits. The docking and binding free energy studies revealed that six hit molecules showed higher binding energy in comparison to the best active molecule. Finally, MD simulations of the top hit with the highest binding energy was carried out which indicated that the obtained hit N1 formed a stable complex with pim-3. We believe that these combined protocols will be helpful and cooperative to discover and design more potent pim-3 inhibitors in the near future.


Introduction
Pim-3 a serine/threonine kinase which belongs to Ca +2 / calmodulin-dependent protein kinase group is a member of provirus integrating site moloney murine leukemia virus (Pim) family. Initially pim-3 was recognized as a novel gene that is induced by depolarization (KID)-1 or forskolin in PC12 cells (a rat pheochromocytoma cell line) [1]. Later, it was renamed as pim-3 owing to its high similarity with that of Pim family of kinases. In the human genome, pim-3 gene is located on chromosome 22q13 encoding a protein with 326 amino acids. Its molecular weight is approximately 35 kD [2,3]. Pim-3 is mainly unregulated in endodermderived tumor tissues, and it is known to inhibit apoptosis and increase cell proliferation in solid and some hematological cancers [4][5][6]. Pim-3 phosphorylates the Ser-112 site of BAD [5,6] in particular and thereby inhibits apoptosis and promotes cancer progression. Additionally by increasing the phosphorylation of STAT3 [7], pim-3 promotes migration and invasion of melanoma cells. Moreover, pim-3 was known to be upregulated in hepatocellular carcinoma (HCC) tumor tissues, and the knockdown of pim-3 in hepatoma cells enhances apoptosis and attenuates cell proliferation [3].
Previous studies have shown that pim kinases including pim-3 are aberrantly expressed in several types of malignancies and helps in inducing tumorigenicity [3,5,6,8,9]. Over-expression of pim-3 mRNA was found in nasopharyngeal carcinoma cell lines and in a group of human Ewing family tumor cell lines. In addition, over-expression of pim was also found in the premalignant and malignant lesions in the stomach, colon, and liver compared with the normal tissues [3,5,6,9]. Furthermore, Pim-3 promotes EWS/ FLI-mediated NIH 3T3 tumorigenesis and hepatocellular carcinoma evident from mice studies [3,10]. According to recent studies, it is also aberrantly expressed in PDAC cells. Pim-3 serves as positive regulator of STAT3 signalling in pancreatic cancer cells [11,12] and is reported to be regulated by transcription factors such as ETS-1.
Computer-aided drug discovery (CADD) techniques [13,14] have paved a new way in the modern era for developing and designing new drugs. CADD methodologies proved to be cost-effective and time-saving in drug discovery processes [15]. Molecular modelling studies such as virtual screening, 3D QSAR, and pharmacophore modelling emerged as promising methods in exploring new leads with drug-like properties. Pharmacophore modelling has been successfully employed in discovering new ALK, JAK2, and D2 inhibitors which are useful in the treatment of cancers such as non-small cell lung cancer (NSCLC) [16], cardiovascular and myeloproliferation disorders [17], and Parkinson's disease, Huntington's chorea, schizophrenia, and tardive dyskinesia [18], respectively. Virtual screening technique gained popularity among scientific community as it decreased the number of hits to be experimentally studied. Furthermore, pharmacophore-based virtual screening abetted the development of specific inhibitors against EGFR and HER2 [19,20]. Therefore, in the present work, we experimented pharmacophore-based virtual screening followed by structure base virtual screening protocols to discover new potent inhibitors against pim-3.

Dataset
A dataset consisting of 57 indole derivatives reported as pim-3 inhibitors were taken for the present study to develop a 3D pharmacophore. The pim-3 inhibitory activity of these derivatives was determined by fluorescence polarization assay method; the activity of these molecules ranged from 5.295 to 8.187 [21][22][23]. The structures of these compounds were sketched in Maestro build panel and subsequently prepared using LigPrep module [24,25]. OPLS-2005 force field was used to generate different conformers, and low-energy conformer of each ligand was taken for further study. The in vitro inhibitory activity of these compounds reported in the literature were converted to pIC 50 using Eq. 1.
(1) pIC 50 = −log 10 IC 50 Pharmacophore generation PHASE (pharmacophore alignment and scoring engine) in Schrodinger [26] was used for pharmacophore generation which actually uses scoring techniques and fine-grained conformational sampling to find a common pharmacophore hypothesis (CPH) [27]. The prepared compounds along with the pIC 50 values were imported and divided into actives and inactives quondam to generate pharmacophore. The ligands with pIC 50 greater than 7.730 were treated as actives and pIC 50 lower than 6.325 as inactive and remaining as the intermediates with moderate biological activity. To determine common pharmacophore hypothesis, PHASE from active ligands employs a tree-based partitioning technique [28]. PHASE is equipped with a built-in set of six pharmacophoric features (aromatic ring (R), positively ionizable (P), negatively ionizable (N), hydrophobic group (H), hydrogen bond acceptor (A), and hydrogen bond donor (D) [27] based on which pharmacophore hypotheses were developed and further utilized for 3D QSAR generation.

3D QSAR model generation
3D QSAR model was achieved through PHASE module of Schrodinger. A chemo metric technique named partial least square analysis (PLS) was used for achieving 3D QSAR model which quantifies the relation between the structure of molecules and experimental biological activities [29]. For generating a 3D QSAR model, initially the dataset was divided into training and test set randomly; the training set was used for deriving 3D QSAR model as it portrayed structural features and biological activities. The test set molecules were used for validating the QSAR model [30]. The pharmacophoric features of the training set compounds were placed into a typical grid of cube of 1 Å spacing [27], each cubes assigned with 0 (or) 1 bit to account for distinct pharmacophore features in the training set. Hence, a single cube will be engrossed with more than one pharmacophoric site expanding one or more volume bits. Thus, depending upon the occupation of various sites in the cube, a single molecule will be expressed with a string of binary values. The biological activity values will be treated as dependent variables, whereas binary values will be considered as independent variables for developing a 3D QSAR model employing PLS regression [27]. Therefore, a series of regression models will be generated by the PHASE [31]. The obtained 3D QSAR model was then verified using test set of molecules using various statistical parameters such as regression coefficient (R 2 ), cross-validated correlation coefficient (Q 2 ), variance ratio (F), standard deviation (SD), Pearson-R, and RMSE (root mean square error) [27,31].

Screening of 3D databases using pharmacophore model
The developed PHASE model was used as a query, and screening was conducted against Asinex and Otava lead library databases to screen out the novel lead molecules that match with the top pharmacophore hypothesis and gauge the predicted biological activities using the generated 3D QSAR model.

Homology modelling and Dock based virtual screening
The crystal structure of pim-3 is not available to date in the protein data bank; hence, the protein was modelled using the sequence information derived from Uniprot database [32] (Accession code: Q86V86). BLAST (Basic Local Alignment Search Tool) [33] was employed for finding the best template from protein data bank (PDB). Based on E value, sequence identity, and secondary structure similarities, a template was selected, and pairwise alignment was done with Clustal-X [34] in order to find out the similarities, identities, conserved regions, and the differences between the target and the template. The human pim-1 (PDB id: 1XWS) showed a sequence similarity of 96% with pim-3. Modeller 9.13 [35] was used for generating a new three-dimensional structure, and various homology models were built out of which the least energy model was taken for further studies. The generated model was analyzed and validated with PROCHECK [36], Verify 3D [37], and ProSA [38]. The modelled protein is equivalent to the low-resolution X-ray crystal structure since the sequence similarity between the template and target was greater than 80% [39].
GLIDE module [40] of the Schrödinger package was employed for carrying out the docking studies. It is regularly used for molecular docking studies which is a grid-based ligand docking approach with notable amount of success. It can predict the protein-ligand-binding modes with decent accuracy and has significantly contributed in identifying numerous potent drug-like candidates against various targets [41,42]. Initially a grid was created centering on the ligand within a cubic box using receptor grid generation panel available in the GLIDE. The default van der Waals scaling was set to 0.9 [43] for non-polar atoms with a partial charge cut-off of 0.25. Later, the molecules which were screened from Asinex and Otava lead library databases were subjected to docking using standard precision protocol. The top 20% of the ligands which showed high dock scores were taken and passed for extra precision (XP) docking. Finally, binding energies were computed for the top 2% of the ligands retrieved from XP docking. Furthermore, the screened hits were pursued for their ADMET profile using Qikprop module of Schrodinger suite [40] and pkCSM online server [44].

Binding free energy calculations
Molecular docking is not always fruitful (or) ultimate method for lead identification. Furthermore, the ranking of molecules based on vigorous binding energy calculations proved to be more propitious than sheer XP dock score [45]. Therefore, molecular mechanics-generalized born surface area (MM/GBSA) calculations were carried out on the top hits obtained from XP docking using prime module available in the Schrodinger. The equation for calculating binding free energy is given below: ΔE MM = Difference in the energy between the protein ligand complex and sum of energies of the protein with ligand and without ligand. ΔG solv = Difference in the GBSA solvation energy between the protein ligand complex and sum of the solvation energies for the ligand and unliganded protein. ΔG SA = Difference in the surface energy between the protein ligand complex and sum of the surface energies for the ligand and uncomplexed protein.

MD simulations
MD simulations were performed on the hit obtained from screening process (N1) using Desmond 3.8 [46]. The stability of docked complex and the effect of conformational changes on the binding interactions were studied. SPC (simple point change) water molecules [47] were used for soaking the complex in orthorhombic box of 15 Å × 15 Å × 15 Å dimensions. The net charges of the system were balanced by adding oppositely charged ions, and the salt concentration was kept 0.15 mol L −1 during the simulation process. The OPLS 2005 force field [48] was used where the system was minimized two times, initially with restraints on solute and later without any restraints. The limited memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) algorithm and steepest descent integrator were employed during minimization until the value concentrated to 5 kcal mol −1 Å −1 . A small relaxation protocol of 12 and 24 ps was applied on the system before proceeding to the extensive MD simulations. The temperature and pressure were maintained at 310 K and 1.01325 bars, respectively, using Nose-Hoover thermostat and Martyana-Tobaisklein approach. The long-range electrostatic interactions were computed using particle-mesh Ewald method [49], and cut-off distance was curtailed to 9 Å for calculating non-bonded interactions. The SHAKE algorithm [50] was exercised on all bonds enmeshed with ΔG bind = G complex − G ligand(unbond) + G receptor(unbond) = ΔE MM + ΔG solv + ΔG SA H-bonds. Finally 50-ns MD simulations at NPT conditions were carried out and unblemished classical trajectories were recorded. The simulation interaction diagrams and root mean square deviation (RMSD) plots were generated from the trajectories in order to inspect the consistency and stability of various protein-ligand interactions.

Results and discussion
Description of the generated pharmacophore The pharmacophore model was generated using a set of ten active (pIC 50 > 7.730) ligands in the dataset with the help of PHASE module available in Schrodinger suite. These active set of ligands may contain structural features which are essential for binding to the active site of pim-3 [51]. The structural and biological activity of pim-3 inhibitors taken for our analysis is depicted in Table 1. The variant list was determined by assigning maximum and minimum number of site to 5 after creating pharmacophore sites for all pim-3 inhibitors. The pharmacophore hypotheses (CPH's) identified in this process were scored and ranked according to the survival active and inactive scores by a scoring method in PHASE module. The inactive ligands though not used in CPH production are used to eliminate the hypothesis that cannot extricate between actives and inactives [51]. Hence, to achieve this purpose, the obtained CPH's will be mapped onto the inactive ligands and scored. The hypothesis with good difference in survival active and survival inactive scores ADRRR was taken for further analysis as it clearly discriminates between active and inactive ligands. Among the generated models, ADRRR was identified as the best hypothesis by scoring algorithm. The CPH with one H-bond acceptor (A2), one H-bond donor (D3), and three aromatic ring features (R9, R10, and R11) is shown in Fig. 1. The distances and angles between various sites of the model are depicted in Figs. 2 and 3, and tabulated in Supplementary Tables Tab S1 and Tab S2, respectively. Furthermore, a 3D QSAR model was derived by pharmacophore-based alignment of the ligands.

Analysis of 3D QSAR model
The PHASE module of Schrodinger suite was used to derive a pharmacophore-based 3D QSAR model using the best CPH, i.e., ADRRR. The derived 3D QSAR model will help in identifying all the features that are crucial for an inhibitor for its gross inhibitory activity against target pim-3 [51]. The dataset was randomly divided into 38 training and 19 test sets (2:1) prior to the model generation. Subsequently with a grid spacing of 1 Å, the QSAR model was developed using PLS analysis. The t-set or eliminate variable was fixed to < 2.0 during PLS regression analysis in order to eliminate those variables whose regression coefficients are very sensitive to tiny changes in the training set [52]. Three PLS factors in total were taken as further increase in the number did not improve the model statistics but over-fitted [53]. The dominant statistics was observed with third PLS factor; therefore, it was chosen to create the QSAR model. To ensure the robustness of the developed model, it was validated both internally and externally against training and test sets using R 2 and Q 2 , respectively [54,55]. The R 2 and Q 2 were found to be 0.913 and 0.748. The R 2 and Q 2 values greater than 0.6 and 0.5 indicate the model is sound and has high ability to predict inhibitory activities in training and test set [54,56]. The model exhibited good statistical significance with Pearson R value of 0.880, variance ratio; F value of 101.2, low standard deviation; SD of 0.263, and root mean square error (RMSE) value of 0.286 [55]. The results are tabulated in Table 2. The experimental and predicted values of the entire dataset ligands are shown in Table 1. The 3D QSAR model was taken further for in silico screening in order to identify novel compounds with good inhibitory activity against pim-3.

3D QSAR visualization of best active compound
The relationship between structure and activity relationship (SAR) and biological activity can be best understood with the help of three-dimensional envision of best hypothesis ADRRR and selected indole ligands in aspect of developed QSAR. The 3D QSAR envision of best active compound (compound 39) with biological activity of 6.5 nM (pIC 50 = 8.187) is shown Fig. 4. In Fig. 4a, b, and c, blue and pink cubes represent hydrogen bond acceptor favored and disfavored regions, respectively; yellow and orange cubes represent hydrogen bond donor favored and disfavored regions, while violet and green cubes represent hydrophobic favored and disfavored regions, respectively.
The atom-based 3D QSAR model while predicting the biological activity takes into the account steric clashes apart from pharmacophore features, whereas pharmacophorebased 3D QSAR predicts the biological activity depending on pharmacophore sites and their areas. From Fig. 4a, it can be seen that the NH of the indole ring, amino group of pyrimidine ring, NH attached to pyrazine ring, and dimethyl amine groups were superposed onto the H-bond donor favored pink contours. One of the nitrogens of pyrazine ring is projecting into the blue H-bond acceptor favored contours (Fig. 4a). It matched with the H-bond acceptor feature A2 in the hypothesis. The NH of indole ring superposed onto the H-bond donor feature (yellow contour) is depicted in Fig. 4b. It explicitly matched with the pharmacophoric feature D3. The dimethyl amino group was onto the H-bond donor favored yellow contour (Fig. 4b). The nitrogens of pyrazine and pyrimidine rings were placed towards hydrogen bond disfavored orange contours (Fig. 4b). In Fig. 4c,   Fig. 3 Inter-pharmacophore angle measurements of the model ADRRR it can be observed that the hydrophilic favored amine group and nitrogens of the pyrimidine ring, one of the nitrogens of pyrazine, and the dimethyl amine were superposed onto the green hydrophilic favored contours (Fig. 4c). The presence

Pharmacophore model-based virtual screening
Asinex and Otava lead library databases prepared earlier using LigPrep was used to screen in order to identify potential pim-3 inhibitors using ADRRR common pharmacophore hypothesis as the template. In this process of searching, the conformers generated from the database were screened to match the CPH based on the site distance. Five out of five pharmacophore sites without partial matches were matched with the hypothesis in the present study. The search process resulted in the 1507 molecules out of total 186,473 molecules from both Asinex and Otava lead library databases. The obtained hits will definitely act as promising pim-3 inhibitors since the model was developed from the active set of compounds with high inhibitory activity against pim-3. Therefore, screened out ligands were further subjected to virtual screening process.

Validation of 3D model
The modelled protein was validated using Ramachandran plot (PROCHECK), ERRAT, Verify 3D, and ProSA (Protein Fig. 5 Ramachandran plot of the modelled pim-3 protein obtained from PROCHECK server Structure Analysis Server). Initially the detailed stereo chemical quality of each amino acid residue was evaluated by Ramachandran plot employing PROCHECK program. The statistics of the modelled protein is shown in Fig. 5. The modelled protein has 295 amino residues out of which 240 (96%) residues were in the most favored regions, 9 (3.6%) residues were in additionally allowed regions, and one (0.4%) residue was in generously allowed regions. Interestingly none of the amino acid residues was in disallowed regions excluding proline and glycine. The results from PROCHECK were promising; hence, the generated model is stereo chemically sterling. The information pertaining to non-bonded interactions between unlike atoms can be known by ERRAT. It is "overall quality factor"; the higher the score, the higher the quality. The acceptable range is > 50 for a high-quality model [57]. The overall quality factor, i.e., ERRAT score for generated model was 73.801 ( Fig. 6) indicating dependability and good stability of the protein [58]. The compatibility of an atomic model (3D) with its sequence (1D) can be analyzed by Verify 3D. The assigned 3D-1D score should never be below 0 and should be above 0.2. In the present study, the value was above 0.2, most part (Fig. 7) indicating good compatibility of the model with its own sequence of amino acids. ProSA was used to compute the interaction energy per residue. In the present case, the interaction energy of each amino acid residue was determined with respect to remaining protein in order to judge the energy criteria. The overall quality of the model can be assessed with ProSA Z-score. The quality of the protein using the Z-score is actually estimated by comparing the experimentally determined proteins with similar number of amino acid residues that were deposited in the protein databank by X-ray crystallographic method (light blue) and NMR spectroscopy (dark blue) (Fig. 8a). The Z-score of the modelled protein was − 7.1 (Fig. 8b). This negative value indicates that the overall quality of the protein is good. All the above investigations suggest that a high-quality model of pim-3 ( Fig. 9) was obtained which can be further used for docking and MD simulations.

Virtual screening using Glide
The 1507 hits filtered out from the Asinex and Otava lead library databases were docked into the active site of pim-3 in two stages. Initially all the hits restituted from the pharmacophore based screening were taken for docking using SP (standard precision) protocol. Subsequently, the selected top scoring 297 (20%) molecules obtained from SP docking were subjected to more rigorous XP (extra precision) docking. Finally binding energies were calculated and the top 6 (2%) molecules having binding energy values greater than the best active compound in the dataset were retrieved.

Binding free energy calculations
The binding free energies of top 6 hits were computed and tabulated in the Table 3. The binding energies of the hits ranged from − 94.698 to − 96.128 kcal/mol. The binding  energy value of the best active compound was − 94.419 kcal/ mol. All the screened ligands exhibited good binding energies greater than the best active compound. Therefore, it is clear that these molecules can definitely have greater probability of becoming potent inhibitors against pim-3.

ADMET profile
The in silico ADMET properties of the screened hits (N1 to N6) retrieved from Asinex and Otava lead library databases were calculated using QikProp module of the Schrodinger interface and pkCSM online server. Percent human oral absorption greater than 80% indicates high absorption profile, while less than 25% indicates poor absorption. Intestinal absorbance less than 30% indicates poor absorption profile [59]. QlogBB represents predicted brain/blood partition coefficient, and the range or recommended value of this property is − 3.0 to 1.2. Similarly QPPCaco and QPlogPo/w represent predicted apparent caco-2 cell permeability in nm/s and predicted octanol/water partition coefficient, respectively. Caco-2 cells are a model for the gut-blood barrier. The QPPCaco value greater than 500 is considered to be great and less than 25 is considered poor. The recommended range of QPlogPo/w is − 2.0 to 6. A categorical yes for CYP2D6 or CYP3A4 indicates possible cytochrome P450 inhibitor. Drug clearance primarily occurs as a combination of renal clearance and hepatic clearance. It is measured by the proportionality constant CLtot and is very important for determining dosing rates to achieve steady-state concentrations. The predicted total clearance log(CLtot) of the screened hits is given in log(ml/ min/kg) ( Table 4). Organic Cation Transporter 2 (OCT2) is a renal uptake transporter that plays a crucial role in disposition and renal clearance of drugs and endogenous compounds. Assessing a candidate's potential whether it can be transported by OCT2 provides useful information regarding not only its clearance but potential contraindications. Inhibition of the potassium channels encoded by hERG I/II are the principle causes for the development of acquire long QT syndrome leading to fatal ventricular arrhythmia. A categorical yes for hERG I/II indicates possible hERG inhibitor for the compound in question. It is essential to consider the toxic potency of a potential compound. The lethal dosage (LD50) values are a standard measurement of acute toxicity and used to estimate the relative toxicity of different compounds. The LD50 is the amount of a compound all at once that causes the death of 50% of a group of test animals. The LD50 (in mol/kg) values of all the screened hits are tabulated in Table 4. The lethal concentration values (LC50) represents the concentration of a molecule necessary to cause the death of 50% of the flat head minnows. LC50 values less than 0.5 mM, i.e., log LC50 < − 0.3 are regarded as high acute toxicity [59]. The screened hits exhibited reasonably good and acceptable ADMET profile (Table 4). However, the properties such as total clearance and renal OCT2 substrate, QPPCaco, hERG I/II inhibitor, and minnow toxicity for the compounds N2, N3, N4, and N6, respectively, were slightly deviating from the recommended values. Interestingly, compounds N1 and N5 though showed comparable values, compound N5 relatively exhibited lower percent human oral absorption and QPPCaco. Therefore, among the screened hits, compound N1 can be considered rationally the best possible drug candidate.

Interaction studies of screened hits
The interaction patterns of the screened hits were analyzed from Ligand Interaction Diagram (LID) available in the Schrodinger suite. The interactions are shown in Figs. 10 and 11. The purple lines indicate hydrogen bond interactions, green lines indicate π-π stacking interactions, red lines indicate π-cationic interactions, and blue-pink line indicates salt bridges. Initially, the pim-3 best active complex was examined with the help of LID. The result is illustrated in

Molecular dynamics simulations
The molecular interactions which are actually accountable for the stability are explored with the help of molecular dynamics simulation. To investigate the stability and molecular interactions of the screened hit (N1) obtained from Asinex and Otava lead library databases, MD simulations were performed for 50 ns.

RMSD (Root mean square deviation) calculations
The RMSD of screened hit (N1) complexed with pim-3 is shown in Fig. 12  time. After 20 ns, the pim-3 complexed with N1 was completely stabilized and showed RMSD value ranging from 3.39 Å to 3.98 Å which is evident from the width of the RMSD plot (Fig. 12). The stability of the pim-3 complexed with N1 could be due to consistent hydrogen bond interactions (i.e., more than 50% of simulation time) of compound N1 (94%, 56%, and 59%) with active site amino acids Glu 124, Glu 174, and Phe 51 (water bridged hydrogen bond), respectively. Therefore, it is evident from RMSD values that the screened hit N1 formed a stable complex with pim-3.

Protein-ligand interaction analysis
Protein-ligand interactions is one of the major factors that is analyzed during the course of simulation. They are useful in understanding the dynamic changes and binding modes during MD simulations. There are four substantial protein-ligand contacts or interactions in MD simulations, namely, hydrogen bonding, ionic, hydrophobic and water bridges. The binding site of pim-3 consisted the following amino acid residues Leu  The screened hit showed hydrogen bonding interactions with Glu 124 and Glu 174. The NH sandwiched between two carbonyl groups of purine ring showed hydrogen bonding with Glu 124 with a maximum occupancy of 94% (Figs. 13 and 14). The maximum is the occupancy, the greater is the contact time between the ligand and the amino acid residue of the protein. The nitrogen of the imidazole ring displayed hydrogen bonding with Glu 174 with significant occupancy, i.e., 29% of the simulation time in the trajectory (Figs. 13 and 14). The hydrophobic interactions were observed with Leu 46, Gly 47, Val 54, Ala 67, Leu 94, Val 106, Ile 107, Leu 123, Arg 125, Leu 177, showed hydrophobic interaction with a maximum occupancy of 34% and 20%, respectively. The amino acid residues Ile 188, Ile 107, Leu 177, and Leu 94 displayed hydrophobic interactions with occupancy of 18%, 13%, 11%, and 11%, respectively, while Val 54, Val 106, Leu 46, and Phe 190 showed hydrophobic interactions of minimum occupancy of 7%, 7%, 4%, and 3%, respectively. Interestingly Phe 51 displayed hydrophobic interaction more than 100% of occupancy due to the fact that it formed both π-π stacking (43%) and π-cationic interactions (73%). On the other hand, water bridges were observed with Phe 51, Glu 174, Pro 126, Leu 46, Arg 125, and Asp 189 with an occupancy of 59%, 26%, 12%, 6%, 6%, and 5%, respectively. Surprisingly the ionic interactions were lost and not observed in the protein complex N1. Therefore, it is evident from interaction diagrams that the interactions with Phe 51, Glu 124, and Glu 174 are essential and might play a critical role in binding of inhibitors to pim-3 and enhance the inhibitory activity.

Protein RMSF
Root mean square fluctuations (RMSF) represent the fluctuations of each amino acid residue of the protein over the period of simulation time. These fluctuations also describe the interaction between the protein and the ligand. The highest fluctuating amino acid residue in the protein pim-3 complexed with screened hit N1 was Gly 83, RMSF value 1.92 Å. The RMSF fluctuation for the remaining protein including the interaction surface was found to be in the acceptable range of 0.42 to 1.67 Å (Fig. 15). The terminal amino acid residues such as Pro 301, Ser 303, and Leu 306 which are far away from the active site region and flexible in nature showed RMSF values of 4.44 Å, 4.01 Å, and 4.11 Å, respectively.

Ligand RMSF
The ligand RMSF plot depicting the fluctuations in the ligand per atom with reference to the protein pim-3 is shown in Fig. 16. The RMSF plot helps in probing the interactions between the ligand fragments and the amino acid residues of the protein, and also crucial role that these fluctuations play in the process of binding. The greater the fluctuations in the ligand fragments, the lesser is its interactions with the target protein. In Fig. 16, X-and Y-axes represent atom index and RMSF values in Å, respectively. The nitrogen of the imidazole ring in the compound N1 was found to fluctuate the most, and the RMSF value of compound N1 was found to be in the acceptable range of 0.55 to 0.93 Å. The carbon atoms (atom 13 and 17) of the imidazole ring and the oxygen atom (atom 14) of carbonyl group sandwiched between two nitrogens of purine ring displayed significant fluctuations with RMSF values of 0.85 Å, 0.80 Å, and 0.84 Å, respectively.

Radius of gyration (rGyr) and solvent accessible surface area (SASA)
The radius of gyration (rGyr) is regarded as one of the most important indicator of protein's folding state in different conditions. It is a measure of extendedness of a ligand and is equivalent to its principal moment of inertia. In the present study, rGyr of screened hit N1 complexed with pim-3 was calculated to gain an insight into the compactness of the protein during the course of simulation. The rGyr showed significant fluctuations during the initial phase of MD simulations (0-20 ns); it got stabilized after 20 ns and fluctuated between 4.01 and 4.35 Å with an average value of 4.17 Å (Fig. 17).
The solvent accessible surface area (SASA) of the screened hit N1 complexed with pim-3 was also measured during the course of simulation in order to explore its exposure to the solvent molecules and access its conformational stability. The SASA value of N1 complexed with pim-3 was found to be in the range of 38.55 to 113.20 Å 2 (Fig. 17); average value was 67.75 Å 2 . The values of SASA although fluctuated substantially during the initial phase of simulation (0-20 ns)   Fig. 17) and remained within the acceptable error once the favorable contacts were established between the protein and the ligand. Thus, the results obtained from rGyr and SASA confirmed the formation of stable pim-3-N1 protein complex.

Conclusions
The aim of the present study is used to identify new potent and promising hits against pim-3 by employing diverse computational methods. To achieve the purpose, pharmacophore-based virtual screening, homology modelling, molecular docking, binding free energy analysis, and molecular dynamics simulation techniques were used. Initially a pharmacophore model (ADRRR) was derived using reported pim-3 inhibitors which were used to screen the subsets of molecules from Asinex and Otava lead library databases. Later molecular docking and binding free energies were performed on the retrieved hits which ultimately yielded 6 hits with binding free energies greater than that of best active molecule. Actually the crystal structure of pim-3 was not available hence homology modelling was employed to model the protein. Finally MD simulations were carried out on the top hit with the highest binding energy which showed stable binding with that of the target pim-3. All in all, the study suggests that integrated 3D QSAR, homology modelling, molecular docking, binding free energy analysis, and MD simulation techniques can be helpful in identifying new pim-3 inhibitors. The inference drawn in the current study can provide some insights to the research community to discover new potent inhibitors against pim-3 with greater biological activity.