Identification of NEK7 inhibitors: structure based virtual screening, molecular docking, density functional theory calculations and molecular dynamics simulations

Abstract NEK7 is a NIMA related-protein kinase that plays a crucial role in spindle assembly and cell division. Dysregulation of NEK7 protein leads to development and progression of different types of malignancies including colon and breast cancers. Therefore, NEK7 could be considered as an attractive target for anti-cancer drug discovery. However, few efforts have been made for the development of selective inhibitors of NIMA-related kinase but still no FDA approved drug is known to selectively inhibit the NEK7 protein. Dacomitinib and Neratinib are two Enamide derivatives that were approved for treatment against non-small cell lung cancer and breast cancer respectively. Drug repurposing is a time and cost-efficient method for re-evaluating the activities of previously authorized medications. Thus, the present research involves the repurposing of two FDA-approved medications via comprehensive in silico approach including Density functional theory (DFTs) studies which were conducted to determine the electronic properties of the Dacomitinib and Neratinib. Afterward, binding orientation of selected drugs inside NEK7 activation loop was evaluated through molecular docking approach. Selected drugs exhibited potential molecular interactions engaging important amino acid residues of active site. The docking score of Dacomitinib and Neratinib was –30.77 and –26.78 kJ/mol, respectively. The top ranked pose obtained from molecular docking was subjected to Molecular Dynamics (MD) Simulations for investigating the stability of protein-ligand complex. The RMSD pattern revealed the stability of protein-ligand complex throughout simulated trajectory. In conclusion, both drugs displayed inhibitory efficacy against NEK7 protein and provide a prospective therapy option for malignant malignancies linked with NEK7. Communicated by Ramaswamy H. Sarma


Introduction
Cancer is an uncontrolled division of abnormal cells (Seyfried & Huysentruyt, 2013). These cells proliferate and destroy surrounding cells and organs, resulting in metastasis. A 2008 survey revealed a global cancer prevalence of around 28.8 million people (Bray et al., 2013). Statistical data from World Health Organization (WHO) affirmed that cancer is the second leading element of mortality in the world causing million deaths annually (Shibuya et al., 2002). Every one out of six deaths are due to cancer worldwide (Prager et al., 2018). Protein kinases are majorly involved in normal cell regulation and cell signaling pathways (Duncan et al., 2010). Human genome contain 497 protein kinases, 284 protein kinases have been found to exist experimentally in apo or inhibitor form (Modi & Dunbrack, 2022). The mutation in protein kinases disrupts its genetic makeup and can cause cell dysregulation which results in abnormal cells growth which is the main checkpoint of cancer (Bononi et al., 2011). Other factors that can cause cancer include: age, gender, lifestyle, dietary changes or exposure to toxins (Fenech & Bonassi, 2011). The important families of kinases include the NIMA family of kinases, Polo-like kinases (PLKs), Aurora, and cyclindependent kinases (CDKs) (Fry et al., 2017). These protein kinases are controlling elements in cell cycle regulation (Fry et al., 2012) and control important checkpoints in cell division. Dysregulation of protein kinases lead to metastasis. Therefore, protein kinase inhibitors plays a significant role in cancer treatment (Grant, 2009). The NEK family of kinases plays a critical role in financing the proliferative pathway of metabolic cells (Malumbres, 2011). NEK family comprised of 11 members (Quarmby & Mahjoub, 2005) out of which NEK6, NEK7 play an decisive character in spindle assembly, cell division and proliferation (O'regan et al., 2007), and NEK7 shares 85% sequence identity with NEK6 protein (De Donato et al., 2018) both these proteins plays important functions in cell division but NEK7 role is unique which cannot be replaced by NEK6 (Saloura et al., 2015;de Souza et al., 2014. Established studies have shown that NEK7 is engaged in a spindle assembly, placement of centrosomes, and cell division. NEK7 is abundantly found in a variety of tissues, including the heart, liver, lung, brain, muscle, testis, leukocyte, and spleen, according to a high-throughput transcriptome research (Kim et al., 2007) where it regulates mitotic progression (Fry et al., 2012). Phosphorylation of a protein substrate elevates its activity or interaction with other molecules, resulting in a variety of physiological responses (Hardie et al., 2003). NEK7 phosphorylates the Kinesin 5 protein (Eg5) at SER 1033 prior to nuclear envelop collapse, causing centrosome separation (O'Regan & Fry, 2009). Interfering with NEK7 can result in unregulated phosphorylation, which increases mitotic cells and finally leads to cell death (Salem et al., 2010). Protein kinases are extensively explored targets during the last 20 years in attempt to develop newer anti-neoplastic medicines (Canduri et al., 2007). There are 53 FDA-approved medications (Gagic et al., 2019) with proven anti-cancer potential in the United States, and over 200 leads are being developed into newer anti-cancer therapies (De Falco & De Luca, 2010;Mullard, 2020), but just a handful FDA-approved inhibitors are effective against NIMA related kinases. Only Dabrafenib has been found to inhibit BRAF-mutant melanoma and NRAS-mutant melanoma cell lines with high NEK9 expression. Dabrafenib revealed a wide range of actions against the NEK9 with an IC50 value in nanomolar range (1-9 nM) (Phadke et al., 2018). However, no FDA-approved drugs exhibit NEK7 inhibitory action, implying that NEK7 inhibitors are scarce .
These findings motivated us to revisiting the potential of FDA approved drugs via a comprehensive in silico investigation. These drugs belongs to Enamide derivatives ( Figure 1) that have established activity against breast cancer and noncell small lung cancer. Repurposing technique could be efficient approach in exploring the new indications for already approved drugs. So we have utilized the detailed in silico approach to identify selective NEK7 inhibitors which can be used as a powerful tool for cancer therapy (Bissantz et al., 2000). Density functional theory investigations were used to determine the reactive characteristics and stability of Dacomitinib and Neratinib (DFTs). The optimization and frequency calculation of both drugs were studied using DFTs. Moreover, HOMO/LUMO analysis, molecular reactivity descriptors, and electrostatic potential (ESP) of both drugs were also evaluated. The binding mode of both drugs was explored by using Molecular Operating Environment (MOE) 2015 (Chemical Computing Group, 2008). The top ranked pose was selected on the basis of docking score and binding affinity (predicted inhibitory constant k i ) and subjected to Molecular dynamic simulations (MD). MD simulation studies is an effective approach to establish the substantiality of complex and evaluation of potential molecular interactions between protein-ligand complexes. In order to determine the best candidate for a therapeutic intervention in breast cancer and other associated tumors, this is the first comprehensive in silico technique to guide future molecular investigation of these derivatives.

Methods and methodology
2.1. Computational studies 2.1.1. Density functional theory calculations The electron density of a molecule is predictive of its electronic energy in the ground state (Hohenberg & Kohn, 1964). The electron density characterizes the quantity of electrons, location of nuclei and nuclear charge in a molecule (Calais, 2009;Parr & Yang Density, 1989). Density functional theory approaches aim to design functionals that can connect the electron density of a molecule to its energy (Ziegler, 1991). The electron density of a molecule can be assumed through set of orbitals utilizing some exchange-correlation like B3LYP (Bartolotti & Flurchick, 1996). DFTs methods have high accuracy (St-Amant, 1996) in making assumptions on electron density of molecule and on these basis, DFT techniques are a dependable and effective tool for accurately estimating the electronic characteristics of compounds (Baerends & Gritsenko, 1997;Marom et al., 2008). In the current study, the optimization and frequency calculation studies of selected drugs were executed on the Gaussian 09W software. Visualization of HOMO/LUMO orbitals and optimized structures were carried out using Gauss view 6 program. The 2D structures of selected drugs were retrieved from PubChem and converted to sybyl Mol2 format which is readable by Gaussian software. The geometries of selected drugs were optimized using DFT/B3LYP level of theory and Karlsruhe-Type Basis Sets (SVP) was utilized for accurate assumptions on the electronic properties of both drugs. The Split valance polarized (SVP) set is a [3s2p] contraction of a (7s4p) set of primitive functions reliable for the geometry optimization of similar structures. Moreover, Karlsruhe-Type Basis Sets (SVP) are designed to provide correct assumptions for the elements up to Kr (Sch€ afer et al., 1992;Weigend & Ahlrichs, 2005). SVP basis set is computationally quite reliable, efficient and available for large fractions of elements in periodic table due to which it is widely used for Hatree-fock/DFT and electron correlation methods (Raji et al., 2020;Sch€ afer et al., 1994). In addition, the DFT/B3LYP approach was used to generate frontier molecular orbitals (FMOs), which are required for calculating the LUMO/HOMO energy gap and the reactivity descriptors. After completion of computational calculations, the output log file was visualized using Gauss view 6 to determine the dipole moment, optimization energy, frequency, and polarizability (Dennington et al., 2019).

Molecular docking studies
Molecular docking is an advance tool for determination of static interactions between ligand and targeted protein. Both FDA approved drugs have promising activity against various cancers (Gagic et al., 2019). In order to determine the binding orientation of these drugs with in active pocket of NEK7 protein, molecular docking approach was undertaken. The Molecular docking studies were performed using Molecular Operating Environment 2015.10 (MOE). Protein Data Bank (PDB ID ¼ 2WQN) was accessed to get the crystal structure of the target protein (Microbiology Department, Faculty of Basic Sciences, Islamic Azad University, Shahr-e-Qods Branch, Tehran, Iran, 2020). Crystallographic structure of NEK7 was co-crystalized with Adenosine Diphosphate (ADP) which can act as a structure guided approach for identification of new inhibitors for NEK7. The unwanted molecules such as water and unnecessary het atoms were removed from the structure using sequence editor utility of the MOE. Protein structure was 3D protonated at pH 7 and 300 K temperature. The protein structure was then minimized to steepest gradient in order to remove any steric clashes. Furthermore, polar hydrogen atoms and gastier charges were added using MMFF94Â forcefield. After the structure preparation, the site finder utility of the MOE was used for identification of active site amino acid residues of the protein. The 2D structures of FDA approved drugs were generated in ChemDraw 12.0 using IUPAC naming obtained from PubChem database (https:// pubchem.ncbi.nlm.nih.gov/) (Gagic et al., 2019). Initially, both drugs were imported to database and saved to MDB format. Afterwards, these drugs were docked into active pocket of protein using two scoring functions. The triangular matcher scoring function was set to London dG whereas placement and refinement parameters were set to GBVI-WSA dG. The induce-fit receptor model was used for docking purpose. After completion of molecular docking, we selected top poses on the basis of Gibbs free binding energy, binding affinity (k i ) and RMSD of the protein-ligand complex. Best poses were further analyzed and visualized using the MOE visualization tool (Prieto-Mart ınez et al., 2018).

Molecular dynamics simulation
The static depiction of the molecular interactions between a protein-ligand complex was generated by the molecular docking method. Molecular dynamic simulations were used to get a dynamic picture of molecular interactions. The Desmond (a Schr€ odinger LLC package) (Desmond Molecular Dynamics System, 2021), was utilized for 100 ns of molecular dynamic simulation. Molecular docking provides the binding orientation of ligand within the active pocket of a corresponding protein (Ferreira et al., 2015), MD simulations, on the other hand, estimate the average movement of atoms relative to some reference. Furthermore, MD simulations give insight into the protein's stability under experimental settings (Hildebrand et al., 2019;Rasheed et al., 2021).
The protein-ligand complex was processed using Maestro or the protein preparation wizard. The system was constructed using the system builder tool of the Desmond suite. The solvation of the complex was carried out using a Monte-Carlo equilibrated water model. Initially, the system was immersed into TIP3P solvation model with a 10.0 Å extension in XYZ direction. To neutralize the system, NaCl ions at a physiological concentration of 0.15 M were introduced as a counter ions. The parametrization of files were carried out using OPLS 2005 (optimized potential for liquid simulation) (Shivakumar et al., 2010) force field (Shivakumar et al., 2010). The 2 ps coupling constant was used in Martyna-Tuckerman-Klein chain coupling technique for pressure control (Martyna et al., 1994). Whereas, for temperature control, Noose-Hoover chain coupling technique was utilized (Hoover, 1985). To eliminate intramolecular steric interactions, the energy minimization process was carried out for 20,000 steps. The system was initially equilibrated by NVT ensemble for 1 ns, followed by a further 1 ns of NPT ensemble at 300 K temperature and 1 bar pressure. Entire MD simulation of 100 ns was executed under periodic boundary conditions. The Particle Mesh Ewald (PME) method was employed (Luty et al., 1994) to account for comprehensive electrostatic interactions observed during MD simulation (Zhang et al., 2020). The algorithm verlet-leapfrog was employed for numerical integration. The time step utilized for minimization was 1 fs, while the time step used for production run was 2 fs (Humphreys et al., 1994). Thermal MM-GBSA.py script (Jacobson et al., 2004;Jacobson et al., 2002) was utilized to compute the ligand strain and putative binding free energy for docked conformations over a period of 100 ns (Bowers, 2006).

SeeSAR analysis
SeeSAR analysis was used to assess the persuasive rationale for binding affinities and hyde energies of FDA approved drugs with NEK7 . The binding affinity and hyde energies of selected compounds with the targeted protein were illustrated using SeeSAR analysis. In SeeSAR software 12.0.1 (Arya, 2020), the top-ranked conformation acquired from molecular docking was subjected to pose generation. Each drug's complex was uploaded to SeeSAR docking mode, which generated a number of poses depending on binding affinity with corresponding protein. The structural features of drugs that contribute favorably to binding affinity were depicted as green coronas, while structural features that contribute unfavorably were depicted as red coronas. The structural characteristics of the drugs that have no bearing on binding affinity were left colorless.

Density functional theory studies (DFTs)
The electronic properties of both FDA approved drugs were determined by DFT calculations including molecular geometry optimizations, HOMO/LUMO analysis, electrostatic potential (ESP), frontier molecular orbitals energies, energetic parameters and global reactivity descriptors. DFT studies were performed using B3LYP method and SVP basis set (Bauernschmitt et al., 1997).

Optimized geometries
In this study, the structural geometries of the FDA-approved drugs Dacomitinib and Neratinib were thoroughly optimized in gas and solvent (water) phase using the DFT/B3LYP approach and the SVP basis set. After optimizing the structural geometry, no imaginary frequencies were discovered, indicating that the present structural geometries are authentic local minima. Dacomitinib's optimization energy was determined to be -1912.731 hartree, indicating that it has reached equilibrium more rapidly. In contrast, Neratinib showed high tendency toward polarizability and dipole moment which indicate more distortion of its electronic cloud as result of external electric field. Optimized structures of both drugs are presented in Figure 2.
In both phases, Neratinib exhibited high values for polarizability and dipole moment, indicating its strong polarity and chemical reactivity. Table 1 provides details on optimization and polarizability.

Frontier molecular orbital analysis (FMOs)
The highest occupied molecular orbitals (HOMO) and lowest unoccupied molecular orbitals (LUMO) analysis provide significant insight into chemical nature of the compounds. In FMOs, electrons are confined to molecular orbitals which may act as electron acceptor or electron donor. The HOMO/ LUMO analysis was performed using SVP basis set in gas and solvent phase. The FMOs of both drugs are depicted in Figure 3.
FMOs have high significance as they determine the reactivity of the compound. The outermost electrons are directly involve in interaction between ligand and target protein so calculation of FMOs become inevitable. The HOMO represents the ability of a compound to donate electrons while LUMO is the ability to accept or withdraw electrons. The small energy gap DE (E LUMO -E HOMO ) between HOMO and LUMO represents efficient charge transfer and making molecule more polarizable. Whereas, molecules having large DE gap (E LUMO -E HOMO ) are least reactive and non-polarizable (Hossen et al., 2021).
The calculated values for HOMO and LUMO and DE gaps (E LUMO -E HOMO ) are tabulated in Table 2.
As shown in Table 2 and Figure 3, the energy band gap between (E LUMO -E HOMO ) for compound Dacomitinib and Neratinib were found to be 0.152 eV and 0.149 eV in gas phase whereas it was 0.178 eV and 0.146 eV in solvent phase respectively. It was observed that small energy gap was corresponding to high reactivity of the compound. Furthermore, chemical hardness and softness is measure of stability of molecular structure (Thanikaivelan et al., 2000). Compound having the large (E LUMO -E HOMO ) energy gap corresponds to chemical hardness of molecule whereas compound having small energy gap are usually more reactive. Dacomitinib showed softness value of 7.40 in solvent phase whereas Neratinib showed chemical softness value of 6.25 which demonstrate that Dacomitinib is more reactive drug than Neratinib. The global reactivity descriptors for both drugs are given in Table 3. Both drugs' ionization energy and electron affinity were expressed using Koopman's theorem (Ehresmann et al., 2004).
The following parameters were examined using their respective formulas: softness: Determination of electron density is inevitable as valance electrons are involved in chemical bonding with other interacting specie. In present study we have calculated the electron density and molecular electrostatic potential (MEP) of both drugs. The electronic density of HOMO orbital in the gas phase for Dacomitinib was focused on the morpholine and piperdinyl rings, while the electron density of LUMO was concentrated on the carbo nitrile and benzocarbazole moiety. The electron density of HOMO and LUMO for Dacomitinib remained nearly same in solvent phase but it was more localized toward carbonitrile ring as compared to gas phase.

Molecular docking
In silico molecular modelling studies were conducted using MOE 2015.10 in order to evaluate and investigate the stability of protein-ligand complex. FDA approved drugs i.e., Dacomitinib and Neratinib were docked into activation loop of NEK7 (PDB ID: 2WQN). Molecular docking studies showed binding affinity and interaction model of drug/protein candidate. The binding affinities of best pose was determined through predicted inhibitory constant (k i ) and docking scores obtained through molecular docking. Activation loop of NEK7 protein is consist of intact HRD and DFG motifs which must remain intact to retain the active conformation of a protein. Both approved drugs were assessed for their binding potential with amino acid residues of active pocket. Discovery Studio Visualizer 17.2 was used for the visualization of the optimal pose (Studio, 2008). The FDA approved drug showed excellent docking scores and produced stable protein-ligand conformation. The detailed interactions produced by both drugs are summarized in Table 4.
The docked conformation of Dacomitinib with in active pocket of NEK7 formed stable protein-ligand complex with least binding energy of -30.77 kJ/mol and predicted inhibitory constant (k i ) 20.37 mM. Predicted inhibitory constant (k i ) is the measure of binding affinity of protein-ligand complex. The amino acid residues which were involved in molecular interactions with Dacomitinib were as follows; ARG42, GLY41, ASP115, ILE40, ILE113, ALA61, LEU111, GLU112, ALA114, VAL48, ILE95, ASP179, LYS63 and PHE168. Dacomitinib showed important molecular interactions including hydrogen bonding, carbon hydrogen bonding, alkyl, p-alkyl and van der Waals interactions. All of these interactions contributed to the stabilization of the proteinligand complex. The 1-ethyl piperidine ring of Dacomitinib  was discovered to form carbon hydrogen bonds with ARG42 and GLY41, with bond lengths of 2.94 and 3.0 angstroms, respectively. The formamide ring was making p-sigma interaction with ILE40, p-alkyl interaction with ALA61 and VAL48, it is also producing an unfavorable acceptor-acceptor interaction with ASP179. Moreover, the ring of dimethylformadine forms a hydrogen bond with ASP179. The length of the hydrogen bond, 1.71 angstrom, represents the strength of the hydrogen bond. The second hydrogen bond was found between the nitrogen atom of the compound and ALA114 of the corresponding protein.
The second hydrogen bond had a bond length of 2.46 angstrom. Fluorobenzene ring is making p-cation interaction with LYS63 and alkyl interaction with VAL48. Furthermore, among non-bonding interactions, van der Waals interactions were observed with following residues; ASP115, LEU113, GLU112, ILE95, PHE168. Figure 4 depicts the probable 2D and 3D binding mechanism of Dacomitinib.  The bonding and non-bonding molecular interactions of drug Neratinib involved important amino acid residues of active pocket of targeted protein. The amino acid residues which were engaged in molecular interactions with Neratinib were as follows; LEU180, LYS63, ASP179, VAL48, ILE40, ALA61, LEU113, GLU112, ALA114, ILE95, LEU111, ARG121, ASP115, GLY117, PHE168, ASP118, ALA165, LYS163 and ASN166. Neratinib showed important molecular interactions which contributed to binding affinity and binding score of protein-ligand complex. Most important interactions were observed between piperidinyl, pyrazol and 2,4 dichloro-4-floro phenyl ring of Crizotinib. Trimethylamine ring was participating in major stabilizing interaction i.e., conventional hydrogen bond with LYS63 and ASP179 respectively of activation loop. The bond length of conventional hydrogen bonding interaction was 4.29 and 2.11 angstroms respectively (Table 4). 5-Methyl-3-methylene-N-vinyl-2,3-dihydropyridin-4-amine was involved in making p-alkyl interaction with ALA65. Moreover, chlorobenzene ring was involved in making p-anion interaction with ASP179, p-alkyl interaction with ALA161 and VAL48. Pyridine ring was forming p-cation interaction with ASP179, p-alkyl interaction with ALA114, ALA61, ILE95 and LEU111. These interactions were majorly contributing toward increasing the binding energy.
Meanwhile, ethoxyethane ring was involved in carbon-hydrogen bonding with ASP115 amino acid having bond length of 3.63 angstrom. Neratinib's binding energy of À26.78 kJ/mol was a result of significant molecular interactions. van der Waals interactions are significant hydrophobic interactions reported with LEU180, ILE40, LEU113, GLU112, ARG121, GLY117, PHE168, and ASN166 amino acid residues. Figure 5 depicts the probable 2D and 3D binding mechanism of Neratinib.

Molecular dynamic simulation
The process of docking is rather quick and inaccurate. The docking defects and flexibility of proteins may impede protein-ligand complex formation. However, while computationally costly and time-consuming, molecular dynamic simulations give a reliable and precise representation of protein displacement. Considering these considerations, simulations of molecular dynamics were conducted using the Desmond software package (Azam et al., 2018;Hospital et al., 2015). The root mean square deviation (RMSD) patterns give important insight into the average change in atomic displacement relative to reference frame. The RMSD trajectory also provides insight into the protein's structural  structure. It is calculated for each trajectory frame. It is necessary to monitor the RMSD of a protein in order to get insight into its structure. Once the protein-ligand complex has been aligned on a reference frame, the RMSD for ligand's heavy atoms can be calculated. If measured values significantly surpass the RMSD of corresponding protein, it is likely that the ligand has dispersed from its initial binding site. In addition, molecular dynamic trajectory analysis is used to calculate the root mean square fluctuation of a targeted protein.

Interpretation of RMSD pattern for NEK7 and protein-ligand complexes
In order to investigate the influence of bound drug on the conformational stability of NEK7 protein, the RMSD patterns for C alpha atoms of NEK7 protein were calculated. The evolution of RMSD values for NEK7 as a function of time is shown in Figure 6. The 2wqn-Dacomitinib combination finds equilibrium after10 nanoseconds of simulation, despite the fact that side chain residues fluctuated within the permitted range of 1-3 Å, which is deemed inconsequential (Choudhary et al., 2020). The NEK7-Dacomitinib complex exhibited minor variations after 40 ns, which stabilized after 60 ns of simulation time and stayed in equilibrium for remaining duration of the simulation. From 40 to 60 ns, the RMSD pattern fluctuated owing to a reduction in the number of connections with amino acid residues. However, after 60 ns, the magnitude of exposures with amino acid residues rose, notably with ASP179, and the RMSD pattern stabilized.  This demonstrates the presence of stable molecular interactions between Dacomitinib and protein residues. After equilibration, the RMSD values of NEK7 remain within 2 Å. After 35 ns, the Protein RMSD exhibited a little fluctuation up to 2.1 Å , then declined again after 50 ns, but stayed below 2 Å and maintained equilibrium for the remainder of the trajectory. These results indicate that the ligands remained tightly attached to the receptor throughout the duration of the simulation. In addition, lower RMSD patterns imply that binding site residues undergo fewer structural rearrangements and conformational modifications (Katari et al., 2016).
Similarly MD simulations for NEK7-Neratinib complex was performed for 100 ns under physiological experimental conditions. RMSD values of NEK7 protein was within permissible range of 1.2 to 2.7 Å (Choudhary et al., 2020), whereas RMSD value of NEK7-Neratinib complex was within 1.5 to 10.5 Å. Average RMSD of NEK7-Neratinib complex was 7.5 Å. It was observed Neratinib bound complex showed slightly higher fluctuations in the beginning and takes up value of almost 10 angstrom. After 20 ns, Neratinib bounded complex attained equilibrium and maintained almost constant value of 7.5 Å for rest of molecular dynamics. The average dynamic oscillation in RMS for complex was ± 0.5 Å. The oscillation in RMS was appeared due to fewer hydrophobic and hydrophilic interactions in the beginning but as dynamics proceeded, Neratinib produced stronger molecular interactions with amino acid residues including ALA61, ILE40, ASP179, PHE168 and LYS63. These findings are indicating weaker molecular interactions in the beginning of dynamics for NEK7-Neratinib complex. These MD simulations studies are complimentary to molecular docking technique and lending testimony on the docking scores obtained from molecular docking as NEK7-Dacomitinib showed relatively good docking score and exhibited greater binding affinity (Figure 7).
Utilizing the root mean square fluctuation to discover local changes in the protein chain is advantageous (RMSF). Figure 8 showing the peaks on the RMSF graph which reflect the areas of the protein that undergo the most variations during simulation. However, the average RMSF for the NEK7 backbone was 1.7 Å, suggesting less structural reorganizations. The N and C-terminal of proteins are more susceptible to modification than any other region. As shown by the RMSF graph, the RMSF value fluctuated significantly between 180 and 220 amino acids throughout the production run, and these residues correspond to the C terminal lobe. The protein's structure, including its alpha and beta helices and strands, is often more rigid and less changeable than its unstructured portion. The highest-peaking residues, according to MD trajectories, are located in loop regions or the N and C-termini (Figure 8). Low RMSF values for binding site residues indicate a stable ligand-protein interaction.
The contact profiles of NEK7-Dacomitinib and NEK7-Neratinib complexes were calculated based on simulated     trajectories, as seen in Figures 9 and 10, respectively. The Dacomitinib interacted with ILE40 and ASP179 through hydrogen bonding and water bridges. Whereas, H-bonding included the amino acid residues ILE40, ALA114, and ASP179. Over 80% of simulation time was spent with ALA114 forming hydrogen bonds. There were hydrophobic interactions among VAL48, LYS63, LEU113, ILE95, ALA61, ARG121, and PHE168 residues. These molecular interactions contributed to protein-ligand complex stability.
The Neratinib produced large number of water bridges contacts with amino acid residues of active pocket. Neratinib interacted with LYS39, LYS63, ALA114, ALA165, ASN166 and ASP179 through water bridges. Neratinib produced fewer Hbond contacts with low interaction fraction. Only LYS163 and LYS39 were involved in H-bonding for lower interaction fractions. It was evident that why RMSD pattern of NEK7-Neratinib complex showed higher fluctuations as compared to NEK7-Dacomitinib complex. In addition, hydrophobic interactions were observed with VAL48, ALA61, LEU111, and PHE168. Hydrophobic interactions with ALA61, and PHE168 was existed for more than 60% percent of simulation time.
The contact map of both complexes were also retrieved form simulated trajectory. Figure 11(A,B) is representing the contact map for NEK7-Dacomitinib and NEK7-Neratinib complex. The   Another important parameter of MD simulation studies is radius of gyration (Rg). Rg indicates the compactness of the protein (Khalid et al., 2022). MD simulated trajectory was used to calculate the radius of gyration of both complexes. Radius of gyration of protein in both cases remained stable through 100 ns trajectory, slight peaks were observed in the trajectory which got stable after few nano seconds of simulated time. Radius of gyration of both complexes is shown in Figure 12(A,B). Important parameters including ligand properties were also computed from 100 ns simulated trajectory. Important parameters and ligand properties as a function of simulated time are shown in Figure 13.

MM-GBSA energy calculations
Molecular docking is a reliable method for assessing the configuration of protein-ligand complex interaction. However, the accurate identification of the binding affinities of docked ligands is still absent. To identify the precise binding energies of docked conformations, MMGBSA energy calculations were done, which is an efficient and trustworthy approach for calculating binding free energies. MMGBSA technique calculates free energy by considering all hydrophobic, hydrophilic, and electrostatic interactions (Vijayakumar et al., 2014). In addition, MMGBSA findings are repeatable since they are insensitive to alterations in solvent model, ligand and protein structural change, and protonation state uncertainty. MMGBSA is applicable for lengthy MD simulations. Compared to the docking score determined from molecular docking, the MMGBSA energy ratings were more negative and indicative of greater binding affinities. Following equation was used to calculate binding free energy (Lyne et al., 2006); Schrodinger's Thermal_mmgbsa script was used to calculate the MMGBSA energies of the complexes. In Table 5, the MMGBSA energies are tabulated.

SeeSAR analysis
The visual depiction of binding affinities of docked conformations were illustrated through SeeSAR analysis. SeeSAR analysis were performed using SeeSAR version 12.0.1. Coronas represented in green and red colors were used to depict the structural features of compounds contributing toward binding affinity. Structural characteristics that improve binding affinity are represented as green colored coronas, while those features that reduce binding affinity are represented as red colored coronas. Structural component of the compounds having no contribution toward binding affinity remained colorless. It was observed that ethyl piperidine ring was contributing favorably toward binding affinity with Hyde energy of -3.7 kJ/mol. Similarly halogen substituted benzene ring was also contributing toward binding affinity of compound. In terms of Neratinib binding affinity, Chlorine substituted benzene ring showed highest positive contribution toward the binding affinity with Hyde energy of -5.7 kJ/ mol. SeeSAR analysis of FDA approved drugs are shown in Figure 14.

Conclusion
In this study, we have investigated the FDA-approved drugs against the NEK7 protein via in silico investigation. DFT studies were conducted to optimize FDA-approved drugs i.e., Dacomitinib and Neratinib, and in-depth electronic and reactivity characteristics were calculated. Both drugs had the potential to interact with NEK7. MOE 2015.10 was used to perform molecular docking studies to determine the fitting of both drugs into the active pocket of NEK7. The results of molecular docking demonstrated a stable conformation with high binding energies and affinity (k i ). Both drugs exhibited strong molecular interactions and had high docking scores. The findings of molecular docking approach was also supported by molecular dynamic simulations. The RMSD trajectory analysis for the protein and protein-ligand complex revealed less fluctuation and a more stable conformation over time. As a result, the current study's findings provide sufficient evidence for the repurposing potential of these two FDA-approved drugs for the alleviation of NEK7-related cancer malignancies. Furthermore, the outcomes of the research imply that these drugs have a selective inhibitory capability against NEK7, and further in-vitro testing of their action against this protein is recommended.

Disclosure statement
There is no conflict of interest.