Identification and antitumor activity of a novel inhibitor of the NIMA-related kinase NEK6

The NIMA (never in mitosis, gene A)-related kinase-6 (NEK6), which is implicated in cell cycle control and plays significant roles in tumorigenesis, is an attractive target for the development of novel anti-cancer drugs. Here we describe the discovery of a potent ATP site-directed inhibitor of NEK6 identified by virtual screening, adopting both structure- and ligand-based techniques. Using a homology-built model of NEK6 as well as the pharmacophoric features of known NEK6 inhibitors we identified novel binding scaffolds. Twenty-five compounds from the top ranking hits were subjected to in vitro kinase assays. The best compound, i.e. compound 8 ((5Z)-2-hydroxy-4-methyl-6-oxo-5-[(5-phenylfuran-2-yl)methylidene]-5,6-dihydropyridine-3-carbonitrile), was able to inhibit NEK6 with low micromolar IC50 value, also displaying antiproliferative activity against a panel of human cancer cell lines. Our results suggest that the identified inhibitor can be used as lead candidate for the development of novel anti-cancer agents, thus opening the possibility of new therapeutic strategies.

SCIENtIFIC RepoRts | (2018) 8:16047 | DOI: 10.1038/s41598-018-34471-y drug design (CADD) strategies also have been attempted to rationally design novel NEK6 inhibitors, but in silico-identified hit compounds were not tested for activity in cell lines 17,18 . Overall, to date no inhibitor of NEKs has entered clinical trials for the treatment of cancer. The need to further investigate the role of NEK6 and the lack of adequate drug discovery efforts prompted us to identify novel compounds with NEK6 inhibitory activity. Here we report a search strategy for new putative inhibitors of NEK6 based on a sequential approach that involves a structure based virtual screening via docking simulations followed by the application of a pharmacophore-based screening to select the best candidates. Twenty-five compounds were identified and in vitro kinase assays demonstrated that the best compound (8) was able to inhibit NEK6 at low micromolar concentrations. Cellular assays subsequently demonstrated antiproliferative activity for the same compound.

Results and Discussion
NEK6 homology modeling. NEK6 and the homologous NEK7 (313 and 302 residues, respectively) are the shortest members of the NIMA family, lacking the regulatory domain and consisting only of a catalytic domain with a very short N-terminal extension to the catalytic domain (NTE, residues 20-33) whose contribute to NEK7 activity was demonstrated 19 . Interestingly, although most of the functions described for NEK6 and NEK7 are very similar, the majority of NEK6 and NEK7 substrates identified to date are specific for one or other kinase 20 .
To obtain a reliable structural model of NEK6, whose crystal structure has not yet been resolved, we applied two different homology modeling approaches: SWISS-MODEL and MODELLER. Initial screening for possible templates was performed using a PSI-BLAST 21 analysis of the amino acid sequence of NEK6 against the PDB resolved structures. The currently available crystal structure of human NEK7, which shows the highest sequence identity (82%), was identified in both approaches as suitable template and the highest resolution structure was selected (PDB: 2WQN) 19 (Fig. 1A). To select the models we relied on energy and stereochemical geometry. The overall stereochemical quality of the models was assessed by PROCHECK. The model structure obtained from SWISS-MODEL featured the majority of residues (98.0%) in the most favored region of the Ramachandran plot (85.9 and 12.1% in the core and allowed region, respectively), with only few residues localized in the generously allowed region (1.6%), and disallowed region of the diagram (0.4%). The corresponding percentages for the MODELLER structure were 98.4%, 1.2% and 0.4%. In order to verify whether the interaction energy of each residue with the remainder of the protein was negative (typically positive values correspond to tricky parts of a model) the protein structures were submitted to the ProSA web server available at https://prosa.services.came.sbg.ac.at/prosa. php. The low ProSA z-score values obtained (−7.37 and −7.33 for SWISS-MODEL and MODELLER models, respectively) are within the range of scores typically found for native proteins of similar size, thus confirming the good quality of our model. The reliability of the modeled fold was also checked with VERIFY-3D, which evaluates the compatibility of a given residue in a certain three-dimensional environment. A score below zero for a given residue means that the conformation adopted by that residue in the model is not compatible with its surrounding environment. The VERIFY-3D scores indicated that the 94.2% (SWISS-MODEL) and 86.6% (MODELLER) of the residues have a 3D-1D averaged score higher than 0.2 (Supplementary Information Figures S1 and S2). To measure the average distance between the backbone atoms of the generated models and template, the Cα-based superimposition RMSD values were calculated using Discovery Studio. Thus superposition of the generated models on the NEK7 template showed good overlap with RMSD values of 0.5 Å (SWISS-MODEL) and 0.6 Å (MODELLER), respectively. On the whole, no major differences were observed between the two generated model structures in these validation analyses. The protein, in inactive conformation 19 , revealed a bilobal fold consisting of a smaller N-terminal and a larger C-terminal lobe connected by a "hinge" (Leu124-Ala127). The N-lobe comprises a five-stranded β sheet and the "αC-helix" which shows an outward rotation characteristics of the inactive state, whereas the C-lobe is mostly α-helical (Fig. 1A,B). The C-terminal domain contains a flexible activation loop (T-loop), 20 amino acids (192-209) marked by a conserved Asp-Leu-Gly ("DLG") motif at the start. Notably in the SWISS-MODEL structure, region G192-T202 adopts a α-helical conformation as observed in NIMArelated, EGRF and Src/Hck kinases [22][23][24] . This short helix, which has been already predicted for NEK6 25 , is not generated using MODELLER. For this reason the SWISS-MODEL-based homology-modeled structure of NEK6 was selected to be used in the subsequent docking simulations for virtual screening of inhibitors.
Virtual screening. To validate the reliability of the docking approach adopted in this work, AutoDock's ability to predict binding poses and discriminate active from inactive compounds were explored. Validation of the accuracy of AutoDock was performed docking the crystallographic ADP cofactor, the unique available native co-crystallized ligand of NEK7, to the 2WQN binding site. A root mean square deviation (RMSD) of the of heavy atoms between the top-ranked pose and the crystallographic one of 1.5 Å was calculated, indicative of a successful scoring function 26 . A data set consisting of 11 NEK6 inhibitors from ChEMBL database (Supplementary Information Table S1) and 612 decoys generated with the Enhanced Directory of Useful Decoys resource (DUD-E) 27 was used to evaluate sensitivity of the docking score using statistical parameters. The receiver operating characteristic curve (ROC), generated to reveal the overall screening performance, showed that the area under the curve (AUC) was 0.79 which represents a good predictive capacity (Supplementary Information Figure S3). We also measured the enrichment of inhibitors among the top-ranking scores of the simulated database through the enrichment factor (EF). We found an EF value of 3.2% for the first 20% of the database and 63.6% of active compounds in the top 20% of results. AutoDock 4.2 was then used for structure based virtual screening. At the time of this study Asinex and Maybridge libraries, as implemented in the ZINC12 database, consisted of 522100 compounds. To focus the virtual screening on compounds that could be promising for further development, we selected a subset of drug-like molecules using an online tool of ZINC12. This resulted in a subset of 6121 compounds. The general workflow of the multistep virtual screening approach implemented in this work is presented in Fig. 2 Figures S4 and S5), structure based virtual screening targeting the ATP-binding pocket and the adiacent allosteric pocket 28 , was carried out. The rank of each compound was determined by the free energy of binding of the lowest energy conformation. The focused library consisting of the top 1000 AutoDock solutions was then screened against the four high ranking pharmacophore hypotheses generated with LigandScout 4.1 which is within the best performing algoritms in compound library enrichment 29 . No three dimensional structure has so far been reported for NEK6, therefore ligand-based pharmacophore models were developed for virtual screening.
The eleven compounds with experimentally tested inhibitory activity to NEK6 available in ChEMBL database (Supplementary Information Table S1) were selected as a training set for pharmacophore model generation. Ten pharmacophore models (Model1-10) were generated by LigandScout, and the results are shown in Table 1. The features selected for this pharmacophore generation were hydrogen bond acceptor (HBA), hydrogen bond donor (HBD), hydrophobic (H) and exclusion volumes (XVOLs) into which a molecule is not allowed to protrude to avoid steric clashes. The default "pharmacophore fit and atom overlap" scoring function (values ranging from 0.79 to 0.74) was used to rank the generated 10 pharmacophore models. The hypotheses were validated using a data set of 12 known inhibitors and 650 decoys to determine how well they were capable to distinguish active from inactive compounds (Supplementary Information Table S2). The first four pharmacophore models, achieving the highest "pharmacophore fit and atom overlap" score (Table 1), were able to cover 91.7% of the active compounds and were selected for virtual screening of our library (Fig. 3A,B). For each pharmacophore model, the fifteen compounds with the highest fit values were selected. The hits were visually inspected for similarity to test a set of structurally diverse compounds. Finally, 25 virtual hits were selected, purchased and subjected to biological testing (Supplementary Information Table S3).
Identification of NEK6 inhibitors. The twenty-five compounds resulting from virtual screening were firstly tested for the ability to inhibit the enzymatic activity of NEK6 at a single concentration of 30 µM, using LANCE technology. Quercetin, able to inhibit a broad panel of kinase, was used as positive control of inhibition 12,30 . This initial screening allowed the identification of two compounds determining >70% NEK6 inhibition at 30 µM, i.e. compound 8 (6-hydroxy-4-methyl-2-oxo-5-[(5-phenyl (2-furyl)) methylene] pyridine-3-carbonitrile), and 21    Figure S6). Thus, while substantially confirming IC 50 value for 8, this latter technical approach gave a value significantly higher for compound 21. Indeed, as it is widely accepted 31,32 , factors related to differences in methodologies of detection and experimental conditions (e.g. buffer composition, pH, enzyme quality, incubation temperatures, ATP concentrations) might account for the differences observed in IC 50  were suggestive of a low intrinsic activity and/or solubility of the molecule: due to this reason, we focused on 8 as a potential lead compound for further investigation. The possible binding mode of compound 8 in the ATP-binding site of NEK6 as well as nonbonded interactions across the complex interface are shown in Fig. 6. Compound 8 may use its phenyl ring and furyl ring to form hydrophobic Pi-alkyl interactions with Ala72 (strand β3), Ala125 (hinge region) and Val59 (strand β2) of the ATP binding site (Fig. 6B). Notably, the dihydropyridinic ring of compound 8, besides interacting with Lys174 of the catalytic loop, forms a short C α -H···O hydrogen bond with Gly192 of the conserved DLG motif (Fig. 6B). This is particulary relevant when considering that C-H···O hydrogen bonds play an important role in the ligand-binding process [33][34][35] .   Kinase selectivity. The structure of NEK6 exhibits a significant identity to NEK7 and thus developing NEK6-selective inhibitors is arduous because of the highly similar (85% identical) catalytic domains 3 . Therefore, we decided to evaluate the selectivity of compound 8 against a panel of NEK kinases at 10 µM, using the technique of Off-chip Mobility shift assay. This concentration was chosen according to the guidelines proposed by Knight and Shokat, since the probability of off-target effects increases dramatically above 10 µM 36 . Furthermore, ATP concentration was set to approximate the K m,ATP for each kinase, so that the resulting selectivity profile directly reflects the intrinsic affinities of the inhibitor 37 . Notably, results obtained showed that compound 8 is very selective against NEK1 and NEK6, with comparable inhibition of 47% and 42.5%, respectively while no significant activity was observed against NEK2, NEK7, and NEK9 (≤0, ≤0 and 3.1% respectively). Considering the high similarity between NEK6 and NEK7, these results may be related to the presence of a glutamic acid rather than lysine in position 200 of NEK6 activation loop (Fig. 1). A major consequence of this substitution may be the redistribution of surface charges at the flexible activation loop resulting in an altered electrostatic potential and ability of NEK6 kinase to interact with other molecules.

Antiproliferative activity of compound 8. A panel of cancer cell lines representative of human tumors,
including breast (MDA-MB-231, MCF-7), ovary (PEO1, COV318), lung (NCI-H1975, NCI-H1299) and colon (HCT-15, SW948), were firstly analyzed for NEK6 and NEK1 expression. Western blot analysis showed that both NEK1 and NEK6 protein were present in the selected experimental models, although at a different extent (Fig. 7). Therafter, the ability of compound 8 in inhibiting cell growth was tested in this panel of tumor cell lines. Results obtained demonstrated that 8 was able to inhibit the growth of MDA-MB-231, PEO1, NCI-H1299 and HCT-15 with IC 50 value below 100 µM ( Table 2). Activity did not correlate with the expression level of the target proteins (i.e. NEK6 and NEK1), this suggesting that other factors may contribute in driving cytotoxicity. Indeed, it is possible that the different susceptibilities of human cancer cell lines to compound 8 might also reflect the genomic diversity of human cancer (Supplementary Information Table S4), in keeping with the notion that cell lines faithfully recapitulate oncogenic alterations identified in tumors and that many of these associate with drug sensitivity/resistance 38 .

Combination study of compound 8 with anticancer drugs on ovarian cancer cells PEO1. Platinum-
and taxane-based regimens represent the standard of care for adjuvant treatment in epithelial ovarian cancer 39 . Therefore, the possibility to use the tested compound in combination with standard cytotoxic drug was investigated in PEO1, cell line representative of high grade serous ovarian cancer harbouring a germline mutation in BRCA2, a gene encoding a protein essential for DNA repair by homologous recombination (HR) 40,41 (Supporting Information,  Table S4). Cisplatin IC 50 value obtained in PEO1 in our experimental conditions was 7.9 ± 0.65 µM (mean ± SEM), in line with literature data 26,42 . Compound 8 showed synergism with cisplatin, resulting in a significant reduction of cisplatin IC 50 from 7.9 ± 0.65 to 0.1 ± 0.01 µM, with combination 8 (44 µM) + cisplatin (10 µM) exhibiting the greatest synergistic effect (Table 3). From a mechanistic point of view, our results are coherent with the evidence suggesting that, besides NEK1, also NEK6 may have a role in DNA damage response, with the identification of several novel NEK6 interactors involved in this process 43 . Given this background, the synergistic effect on PEO1 between 8 and cisplatin could derive by genomic instability resulting from different processes, including (a) PEO1 inability to repair DNA damage by homologous repair and nucleotide excision repair 40,44 , (b) inhibition of the DNA damage repair activity exerted by NEK1 and, possibly, by NEK6 and, finally (c) by the impairment of NEK6 and mitotic functions that may lead to further abnormal centrosome functions, spindle defects and chromosome mis-segregations. With regard to taxane, the suggested role of NEK kinases in microtubule dynamics regulation supported our investigation to evaluate the effect of a combination between compound 8 and paclitaxel on PEO1 cells. Indeed, previous findings have shown that NEK6 and NEK7 may directly regulate microtubule organization 45 , an aspect that makes them attractive targets for drugs that would be complementary to microtubule targeting agents 2 . The   42,46 . Results obtained showed a slight synergistic effect at 1 nM paclitaxel, shifting IC 50 paclitaxel from 7.0 ± 0.001 to 0.64 ± 0.057 nM (Table 4). However, at higher paclitaxel concentrations antagonistic effects were observed. Although this is not an unusual finding, since synergism and antagonism can be different at different dose levels 47 , it certainly requirer further investigation to properly identify the best combination drug schedule. Results from the combination studies are in line with our previous findings demonstrating that stable NEK6 silencing in ovarian cancer cells produced a modest but significant increase of cisplatin and paclitaxel sensitivity, while a stable overexpression slightly increased drug resistance 11 .

Conclusions
Here, we employed in silico screening techniques to search for novel NEK6 inhibitors. The best identified hit, i.e. ((5Z)-2-hydroxy-4-methyl-6-oxo-5-[(5-phenylfuran-2-yl)methylidene]-5,6-dihydropyridine-3-carbonitrile) is able to inhibit NEK6 at micromolar order of magnitude. Interestingly, 8 shows selective inhibition for NEK6 and NEK1, while not acting on NEK7; the identification of a compound able to inhibit NEK6 with respect to NEK7 is an important achievement, considering that NEK6 and NEK7 have the highest degree of sequence identity and differ only for one residue in the active site 48 . Compound 8 shows antiproliferative activity against a panel of human cancer cell lines, and displays a synergistic effect with cisplatin and paclitaxel in a BRCA2 mutated ovarian cancer cell line, this supporting a possible use for personalized therapy. On the whole, the newly discovered inhibitor deserves consideration for further development by SAR studies to optimize its activity.

Methods
Molecular modeling. We implemented standard computational strategies to construct the tridimensional model of human NEK6, since its experimental structure is not yet available. The sequence was obtained from the UniProt database (UniProt accession number Q9HC98) and the modeling was performed using two different homology modeling approaches: the fully automated structure modeling method implemented in the SWISS-MODEL workspace (http://swissmodel.expasy.org//SWISS-MODEL.html) 49 and the program MODELLER 50 , as implemented in DiscoveryStudio 4.5 (Dassault Systèmes BIOVIA), which generates a protein model through the satisfaction of spatial restraints. The crystallographic structure of human NEK7 (PDB code 2WQN) 19 was identified by PSI-BLAST 21 as tertiary structural template showing 82% sequence identity. The generated models were evaluated using the programs PROCHECK, VERIFY3D and ProSA-Web 51-53 .
Structure based virtual screening. The first virtual screening step was carried out using docking simulations to generate a library containing molecules having appropriate structural complementarity with the search zone defined around the ATP binding site in NEK6. The atomic coordinates of NEK6 obtained from the homology modeling were used as the receptor model in the virtual screening with docking simulations using AutoDock 4.2.6 54 . To verify the predictive ability of AutoDock scoring function for this target, a data set consisting of 11 NEK6 inhibitors from ChEMBL database (http://www.ebi.ac.uk/chembl/) (Supplementary Information  Table S1) and 612 decoys generated with the Enhanced Directory of Useful Decoys resource (DUD-E) 27 was used. Enrichment calculation and receiver operator characteristic (ROC) curve analysis were performed using the Enrichment Calculator tool of Schrödinger suite (Schrödinger LLC). Compounds from suppliers Asinex and Maybridge were downloaded from the ZINC12 database 55 selecting the drug-like subset 56 . The selected descriptors were the molecular weight (150 < MW < 500), the hydrophobicity (−4 < xlogP < 5), the net charge (−5 < NC < 5), the number of rotable bonds (RB < 8), the number of H-bond donors (HBD < 10), the number of H-bond donors acceptors (HBA < 10), the polar surface area (Å 2 ) (PSA < 150), the polar desolvation (kcal/mol) (−400 < PD < 1), the apolar desolvation (kcal/mol) (−100 < AD < 40). The filtering procedure yielded 6121 compounds. Both the NEK6 receptor and the small molecules were prepared for docking using the AutoDockTools (ADT) software package 57   for each atom type in the ligands (A, C, HD, N, NA, OA, SA, Cl, I, F, S, Br) plus one for electrostatic interactions, were computed covering the putative ATP binding site using 60 × 60 × 60 grid points with a spacing of 0.375 Å. A Lamarckian genetic algorithm was chosen, and default parameters were used except "Number of GA runs", "Population size", and "Maximum number of evaluations", which were set to 10, 50, and 2.5 × 10 5 . Results were clustered according to all-atom RMSD (root mean square deviation) with a tolerance of 2 Å and representative binding mode was defined as the lowest-energy complex of the cluster with the largest population.
Pharmacophore model generation and ligand based virtual screening. LigandScout 4.1 59 , via a ligand-based strategy, was used for 3D pharmacophore generation, refinement and virtual screening. This approach searches for a common feature pattern that is shared in an active ligand set and considers the conformational flexibility of the ligands. The training set consisted of the eleven ChEMBL compounds with inhibitory activities to NEK6. Twelve kinase inhibitors from the DrugKiNET database (www.drugkinet.ca) were taken as a test set for pharmacophore model validation using a decoy set of 650 inactive compounds assembled from DUD-E database 27 (Supplementary Information Table S2). Excluded volumes representing the sterically occupied regions by the receptor, were taken into account to increase the selectivity of the model. All LigandScout parameters were kept default and the models were ranked based on the "pharmacophore fit and atom overlap scoring function". The four best-fitting models exhibiting good score values were selected and applied to a 3D virtual focused library of 1000 compounds represented by the top 1000 AutoDock solutions.
NEK6 kinase assay. Compounds  Cell cultures. The following human carcinoma cell lines were used for the study: MDA-MB-231 and MCF-7 (breast cancer; ATCC, Sesto San Giovanni, MI), PEO1 and COV318 (ovarian cancer; ECACC, Salisbury, UK and ATCC), HCT-15 and SW948 (colon cancer; Public Health England, Salisbury, UK and ATCC), NCI-H1975 and NCI-H1299, (lung cancer; ATCC). Cells were routinely tested free of mycoplasma (MycoAlert mycoplasma detection kit, LONZA, Rockland, ME, USA) and validated by STR (Short Tandem Repeat) DNA profiling (BMR Genomics srl, PD). Cells were grown following supplier indication, in a fully humidified atmosphere of 5% CO 2 /95% air. All the reagents were purchased from Sigma-Aldrich (St. Louis, MO), if not otherwise specified.
Western blot analysis. Western blot analysis was carried out as previously described, following cell incubation in lysis buffer containing 20 mM Tris-HCl pH 7.4, 5 mM EDTA, 150 mM sodium chloride, 1% glycerol e 1% Triton X-100, in the presence of a cocktail of both protease and phosphatase inhibitors 11