Molecular modeling and docking of small molecule inhibitors against NEK2

Aberrant expression of NEK2 (NIMA-related kinase 2) is indicated in a wide variety of human cancers. NEK2 is highly correlated to multi drug resistance by activating drug efflux activity. Identification of new small molecule inhibitors targeted against NEK2 therefore, facilitates to increase drug sensitivity of cancer cells, by stabilizing drug influx and minimizes the dose of therapeutic drug. Our work investigates to screen for optimal small molecule inhibitors against NEK2. In this study, we used a computational approach by modeling NEK2 protein using I-TASSER (Iterative Threading ASSEmbly Refinement) software. The modeled structure was subjected to protein preparation wizard; to add hydrogens and to optimize the protonation states of His, Gln and Asn residues. Active site of the modeled protein was identified using SiteMap tool of Schrodinger package. We further carried out docking studies by means of Glide, with various ligands downloaded from EDULISS database. Based on glide score, potential ligands were screened and their interaction with NEK2 was identified. The best hits were further screened for Lipinski’s rule for drug-likeliness, bioactivity scoring and ADME properties. Thus, we report two (didemethylchlorpromazine and 2-[5-fluoro-1Hindol- 3-yl] propan-1-amine) compounds that have successfully satisfied all in silico parameters, necessitating further in vitro and in vivo studies.

Various drugs that can either inhibit or by-pass drug efflux pumps include annamycin [3], tetrandrine [4], and toremifene [5] as well as utilizing nanocariers [6]. A new trend of poor substrates for efflux pumps was also developed [7]. In spite of this, the outcomes were often poor with agents against multidrug resistance and most of them have failed due to heavy toxicity or drug metabolism related issues, limited significance, non-efficacious in clinic and/ or have not been tried even in clinical settings [2]. Hence, new inhibitors focused on novel therapeutic targets may throw light into more molecular targeted therapies that could specifically and potentially bring down MDR. Figure 1: a) Modeled NEK2 protein, using I-TASSER software; Visualization of modeled protein NEK2 using Maestro, depiction of secondary structures like helix (red), sheet (Cyan) and loop (Gray); b) Active site of NEK2; Ligand binding pocket of NEK2 is shown with predicted hydrophobic regions in yellow mesh, hydrogen bond donor surfaces in blue mesh and acceptor surfaces in red mesh NEK2 (NIMA-related kinase 2) is a MDR protein, that has gained significant importance in recent years, of controlling several MDR pathways of the ABC (ATP-binding cassette) transporter family, such as ABCB1, ABCG2, ABCC1 and was found to increase chromosomal instability, cell proliferation. NEK2 was found to act mainly through efflux pumps to induce drug resistance [8,9]. Regarding the essentiality of NEK2, NEK2 is a mitotic protein kinase (Serine/Threonine) of the NEK family and has a role in regulating centrosome separation. But compared to other mitotic kinases, NEK function is subtle and neither suppression nor silencing dramatically affected cell cycle and play only a supportive role in centrosome separation [10,11]. In normal cells, NEK2 is expressed in traces or mostly undetectable, but over expressed and found localized in the nucleus of cancer cells [12,13]. NEK2 controls key MDR proteins in a variety of cancers such as multiple myeloma, myeloid leukemia, lung adenocarcinoma, mantle cell lymphoma, mesothelioma, head and neck squamous cell carcinoma, bladder carcinoma, glioblastoma, T-cell acute lymphoblastic leukemia, colon carcinoma, hepatocellular carcinoma, melanoma, and ovarian adenocarcinoma [8,14] ,breast cancer [15,16], non-small cell lung cancer [17], diffuse large B cell lymphoma [18,19]. Several in vivo studies have shown that ablation of NEK2 with siRNA has improved sensitivity of tumor to respond to drug treatment [20,21]. Concerning clinical studies, many reports suggests that, abnormal expression of Nek2 at the protein level might be one of the mechanisms of tumorigenesis and indicate that, Nek2 may represent a new potential target for therapeutic intervention [8,17,[21][22][23][24][25][26]. In this study, we used a computational approach using I-TASSER to model the target protein NEK2 and GLIDE tool (Schrodinger) to identify and to validate reliable small molecule inhibitors against NEK2.

Protein modeling
The amino acid sequence of human NEK2 was retrieved from UniProt (P51955). It is a Serine/threonine protein kinase containing 445 amino acids. The FASTA sequence of NEK2 was subjected to BLAST, choosing Protein Data Bank database to identify appropriate template. Full length protein template was not identified. Hence, using Iterative threading assembly refinement (I-TASSER) server, we modeled the protein. I-TASSER server is an integrated platform for automated protein structure and function prediction based on the sequence to structure paradigm. It generates three dimensional atomic models from multiple threading alignments. An estimate of accuracy, of the predictions was provided based on the confidence score of the modeling and the best scored model was identified and utilized for model validation. The protein was modeled by using Protein Preparation Wizard of Schrodinger Suite; the protein structure was prepared by adding hydrogen atoms, optimizing hydrogen bonds and verifying the protonation states of His, Gln and Asn. Energy minimization was carried out using default constraint of 0.3 Å RMSD and OPLS 2005 force field. SiteMap tool was used to identify binding pockets of NEK2. The site with a site score above 1, was selected for grid generation using Glide grid and docking studies were carried out.

Molecular docking
The ligands from different databases (ChemBridge, Maybridge, Pubchem, Sigma Aldrich, Specs) were downloaded from EDinburgh University LIgand Selection System (EDULISS) [27]. Schrodinger ligand preparation product, Ligprep was used to prepare high quality, all atom 3D structures. The ligand preparation included 2D-3D conversions, generating variations, correction, verification and optimization of the ©2016 structures. Receptor grid was generated using Receptor grid generation in the Glide application (Glide, version 5.8, Schrödinger, LLC, and New York-2) of Maestro (Schrödinger, LLC, New York, NY, 2014-2). The receptor grid for NEK2 was generated by specifying the binding (active) site residues, which was identified by SiteMap tool. Once the receptor grid is generated, the ligands are docked to the protein (NEK2) using Glide version 5.8 (Grid based LIgand Docking with Energetics) docking protocol. The ligands were docked using "Extra precision mode" (XP). The docked conformers were evaluated using Glide (G) Score. The G Score is calculated as follows: G Score = a*vdW+b*Coul + Lipo + Hbond + Metal + BuryP + RotB + Site Wherein vdW denotes van der Waals energy, Coul denotes Coulomb energy, Lipo denotes lipophilic contact, HBond indicates hydrogen-bonding, Metal indicates metal-binding, BuryP indicates penalty for buried polar groups, RotB indicates penalty for freezing rotatable bonds, Site denotes polar interactions in the active site and the a=0.065 and b=0.130 are coefficients of vdW and Coul.  Arg129. Arg129 and Arg130 shows Pi-Cation interaction; c). Interaction of NSC132835 with NEK2 is shown; the compound forms hydrogen bonds with Glu126 and Glu292. Arg130 was engaged in Pi-Cation interaction ADME property analysis ADME properties of selected ligands were analyzed using QikProp tool of Schrodinger suite. The tool predicts the physiochemical properties with a detailed analysis of: (i) predicted IC50 value for blockage of HERG K+ channels (ii) Caco-2 cell permeability (iii) brain/blood partition coefficient (iv) apparent MDCK cell permeability (v) prediction of binding to human serum albumin (vi) percentage of human oral absorption.

Modeling of nek2 and docking
In the present study, in silico docking studies were carried out, to screen and identify appropriate ligands that can fit into the most favorable binding mode, against NEK2 modeled protein.
The amino acid sequence of NEK2 was retrieved from Uniprot with accession number P51955. The structure was predicted by I-TASSER. This software implements a hierarchical approach to protein structure prediction. Structural templates are identified using multiple threading approach and full length atomic ©2016 models were constructed by iterative template fragment assembly simulation. The modeled structure was subjected to protein preparation to add missing hydrogen, optimize the protonation states of the His residues and orientation of hydroxyl groups. The structure was also subjected to energy minimization using the OPLS_2005 force field to constrain the heavy atoms. Figure 1a shows the structure of NEK2 after protein preparation.
Active site of NEK2 was obtained (Figure 1b) using SiteMap tool, which provides a fast and effective means of identifying potential binding pockets of proteins. SiteMap identifies the character of binding sites using novel search and assesses each site by calculating various properties like size, volume, amino acid exposure, enclosure, contact, hydrophobicity, hydrophilicity and donor/acceptor ratio. As a result, three binding sites with site score more than 1 were identified. Site Table 1. Docking score for the compound NSC636674 was -9.46, this compound establishes two hydrogen bonds with Glu126, Glu292 and pi-pi stacking with Arg130. The binding mode of interaction is shown in Figure 2a.

Drug likeliness and bioactivity score
Molecular properties of the selected compounds are read using Molinspiration software to satisfy lipinski's rule of five, which is essential for rational drug design. All the compounds showed, no violation of all the five rules; not more than 5 hydrogen bond donors, not more than 10 hydrogen bond acceptors, molecular weight of compounds less than 500, partition coefficient (log P) less than 5, rotatable bonds less than 10, topological polar surface area (TPSA) of not greater than 140 (Table 2). Bioactivity scoring indicated, two compounds; didemethylchlorpromazine and 2-[5-fluoro-1H-indol-3-yl] propan-1-amine, to be a significant ion channel modulator with a score of 0.41 and 0.21 respectively ( Table 3). The obtained scores are in the range of significant score value of above 0.0. For each compound, the bioactivity contribution will be calculated for each substructure of fragment; the bioactivity for the entire molecule will then be calculated, as a sum of activity of contributions of all the fragments in a molecule. This provides a molecule activity score (a number, typically between -3 and 3). It has been recommended by Molinspiration that, molecules with the highest activity score have the highest probability to be active.

ADME prediction of ligands
QikProp tool predicted significant ADME properties such as permeability through MDCK cells (QPlogMDCK), log IC50 value for blockage of K+ channels (QPlogHERG), gut-blood barrier (QPPCaco) and predicted brain/blood partition (Qplog BB) ( Table 4) . Percent of Human Oral absorption is based on number of metabolites, binding to human serum albumin (QPLogKhsa), predicted aqueous solubility (QP logS) and cell permeability.
Physiochemical properties and ADME parameters of these ligands, confirms that they can be considered as drug candidates for further studies.

Discussion:
MDR to functionally targeted therapies against cancer is of current concern. Various strategies have been tried recently, to effectively bring down MDR but with limited success. Hence, new inhibitors focused on novel therapeutic targets of pathways involved in MDR are a compelling requirement, which could specifically and potentially bring down MDR. Several groups have reported to negating, important molecular targets down the line of MDR pathway. In one report, through homology modeling and molecular docking, several chemotherapeutic inhibitors were screened against MDR1 or ABCB1 and it was found that the drug paclitaxel can be used for repurposing against MDR in non-small cell lung cancer [31]. Pharmacogenomics and molecular docking studies have shown that apigenin can inhibit ABC transporters such as Pglycoprotein and BRCP [32]. Nilotinib is shown to potently sensitize specific anticancer agents by blocking the functions of ABCB1/P-glycoprotein, ABCG2/BCRP and ABCC10/MRP7 transporters involved in MDR [33]. NEK2 protein is well known to induce MDR [8]. Although conflicting, NEK2 is not ©2016 essential and plays only a supportive role in centrosome separation [10,11]. NEK2 represents an interesting multifaceted therapeutic target [11,34], that can regulate a multitude of MDR pathways such as ABCB1, ABCG2, ABCC1 [8]. Indeed, inhibitors are designed against NEK2 through structure based design such as aminopyridine derivatives [11] and inhibitors against NEK2 successfully affected proteasome activity of cancer cells [10]. Therefore, it is evident that small molecule inhibitor designing and pharmacological modulation of NEK2 may affect several mechanisms important for tumor growth such as drug resistance, survival, progression, and metastasis [34 -36]. In this study, we employed a computational approach to screen small molecule inhibitors against NEK2. Binding affinities of the ligands with NEK2 are analyzed, evaluated based on the docking scores and the number of intermolecular hydrogen bonding interactions of the resulting receptor-ligand complexes. The top score compounds showed, minimum of three hydrogen bonding interaction with amino acids, present in the binding pocket. Compound NSC636674 showed Pi-Pi stacking. Pi stacking is prevalent in protein crystal structures, and contributes to the interactions between small-molecules and proteins. The shape and electronic properties of aromatic rings responsible for large polarizabilities and a considerable quadrupole moment resulted in preferred interaction geometries [37]. The geometry of these stacking suggests that, electrostatic interaction, play a role in the attraction between aromatic chain of a compound or amino acid with the amino acid. Further studies suggested that the stacking interaction, contributes to the binding energy of receptor-ligand [38]. Compounds showed no violation of all the five rules; not more than 5 hydrogen bond donors, not more than 10 hydrogen bond acceptors, molecular weight of compounds less than 500, partition coefficient (log P) less than 5, rotatable bonds less than 10, topological polar surface area (TPSA) of not greater than 140. AE-641/00584006 0.14 -0.06 0.14 -0.36 0.20 Data obtained by calculation of sum of activity of contributions of all the sub-structure fragments in a molecule. Molecules with the highest activity score have the highest probability to be active. *indicates significant scoring for Ion channel modulator.
Other compounds established cation -Pi bonding with the active site amino acids. Cation-pi bond interaction is a very important contributor of protein architecture and stability [39 -42]. The pi system is provided by the aromatic ring of the small molecule and the protein cation (Arg 130). Energetically, the cation-pi interaction is comparable to or stronger than a hydrogen bond [43]. Studies suggested that the cation-pi interaction is a powerful force that aids in the recognition between proteins and ligands and is a valuable predictor of drug-receptor interactions [44].
In the current study, all the compounds that showed good binding affinity, also exhibited drug like characteristics based on Lipinski's rule of 5 that determines if the compound, has certain pharmacological or biological activity to make it an orally active drug in humans [45]. The molecular weights of all the compounds are below 500 Daltons, with less than 5 hydrogen bond donors and 10 hydrogen bond acceptors. In addition, analysis of pharmacokinetic properties such as the partition coefficient and Water solubility (QPlogS) of the evaluated compounds are within the range. The top 10 compounds also showed good cell permeability (QPlogKhsa), bioavailability (QPCaco, QPMDCK) and high serum protein binding capacity (QPlogBB). All the pharmacokinetic parameters are within the acceptable range defined for human use, which collectively indicated that the screened compounds could be taken forward, for further analysis.

©2016
All the selected 10 ligands showed good docking scores reflecting drug-binding affinities with NEK2. All the selected ligands showed favorable molecular properties by satisfying Lipinski's rule of 5 and ADME profile. Two compounds out of the 10 screened compounds, showed significant bioactivity scores as ion channel modulators; didemethylchlorpromazine and 2-[5-fluoro-1H-indol-3-yl] propan-1-amine. Interestingly, chlorpromazine derivatives are known to be active against MDR [46,47]. We speculate that the activity of didemethylchlorpromazine against MDR could be due to inhibition by binding to NEK2 from our current in silico study.

Conclusion:
NEK2 is associated with several human cancers with known MDR data. Hence, it is important to develop improved inhibitors for NEK2 negating MDR. In this study, we described modeling of NEK2 and docking studies to identify inhibitors against NEK2. Two compounds namely di-demethylchlorpromazine and 2-[5-fluoro-1H-indol-3-yl] propan-1-amine showed good binding affinity and better ADME properties. The results obtained from this study may be worthwhile, to carry out further in vitro and in vivo studies to design novel and potential inhibitors against NEK2. (a) Predicted IC50 value for blockage of HERG K+ channels (Acceptable range limit -above 5.0), (b) Predicted Caco-2 cell permeability in nm/sec (<25 poor, >500 great), (c) Predicted brain/blood partition coefficient (Acceptable range -3.0 to 1.0), (d) Predicted apparent MDCK cell permeability in nm/sec (<25 poor, >500 great), (e) Prediction of aqueous solubility in mol/L (Acceptable range -6.5 to 0.5), (f) Prediction of binding to human serum albumin (Acceptable range -1.5 to 1.2), (g) Percentage of human oral absorption (<25% is poor and >80% is high).

Declaration of interest statement:
The authors report no declarations of interest.