Computational identification of potential inhibitors targeting cdk1 in colorectal cancer

Introduction: Despite improved treatment options, colorectal cancer (CRC) remains a huge public health concern with a significant impact on affected individuals. Cell cycle dysregulation and overexpression of certain regulators and checkpoint activators are important recurring events in the progression of cancer. Cyclin-dependent kinase 1 (CDK1), a key regulator of the cell cycle component central to the uncontrolled proliferation of malignant cells, has been reportedly implicated in CRC. This study aimed to identify CDK1 inhibitors with potential for clinical drug research in CRC. Methods: Ten thousand (10,000) naturally occurring compounds were evaluated for their inhibitory efficacies against CDK1 through molecular docking studies. The stability of the lead compounds in complex with CDK1 was evaluated using molecular dynamics simulation for one thousand (1,000) nanoseconds. The top-scoring candidates’ ADME characteristics and drug-likeness were profiled using SwissADME. Results: Four hit compounds, namely, spiraeoside, robinetin, 6-hydroxyluteolin, and quercetagetin were identified from molecular docking analysis to possess the least binding scores. Molecular dynamics simulation revealed that robinetin and 6-hydroxyluteolin complexes were stable within the binding pocket of the CDK1 protein. Discussion: The findings from this study provide insight into novel candidates with specific inhibitory CDK1 activities that can be further investigated through animal testing, clinical trials, and drug development research for CRC treatment.


Introduction
Colorectal cancer (CRC) is the third most common cancer next to lung and breast cancers with an estimated diagnosis of 2 million individuals and a mortality of approximately 900,000 annually (Bray et al., 2020;Holowatyj et al., 2020).It is thought to be a disease in high-income countries; however, studies have shown the incidence of CRC to increase at a projected rate of 70% by 2030 for low and middle-income countries (Ferlay et al., 2015) with a 5% lifetime average-risk incidence globally (Weitz et al., 2005).These gory statistics have significant implications for the public health systems both in the country and globally hence the need for a more dynamic scientific approach to the CRC burden.
The onset of CRC follows a pattern of abnormal growth of the epithelial cells of the colon and/or rectum leading to a tumor mass that fills the intestinal lumen, metastasizing into nearby organs of the abdominal area, lymph nodes, and other parts of the body (Adeoti et al., 2016).It usually begins as a small, adenomatous mass of cells known as polyps that form on the inside of the colon, which may become malignant or remain benign (Ikwu et al., 2020).As with many cancers, the cause of CRC is not known (Simon et al., 2017), however, risk factors such as age, lifestyle habits (alcohol and tobacco use), high-fat diet, predisposing medical conditions such as type 2 diabetes, and a family history of CRC have been associated with the disease course.Although traditional chemotherapy options have yielded remarkable results in CRC treatment with possibilities of relapses, limited molecular understanding of the oncogenetic initiation and interplay of biochemical mechanisms/factors have impaired the development of novel pharmacologically therapeutic anti-cancer drug agents.
Cyclin-dependent kinase 1 (CDK1) enzyme has been implicated in the oncogenesis of CRC (Sofi et al., 2022a).It is a serine/threonine kinase known to be the catalytic subunit of the M-phase promoting factor (MPF) complex (Kalous et al., 2020) and belongs to a family of cyclic cell regulatory proteins involved in maintaining cell cycle efficiency (Sofi et al., 2022a;Sofi et al., 2022b).Though involved in many cellular events that ensure cell survival, CDK1 specifically acts by regulating the G 2 /M and G 1 /S phases oscillating in activity between the M transition stages to bring about replication, activating checkpoints, and more recently has been essential to replication fork maintenance (Adhikari et al., 2012;Liao et al., 2017;Matthews et al., 2022).Prior to initiating cell replication, CDK1 is activated by binding to cyclin B 1 to form CDK1-cyclin B complex (Figure 1) (Brown et al., 2015).However, the over expression of CDK1 has been linked to unrestricted cell proliferation.Its over expression has been linked to unrestricted cell proliferation, tumor development, and dysregulation associated with the initiation of many cancers including colorectal cancer (Li et al., 2020;Sofi et al., 2022a;Sofi et al., 2022b).Gene expression has been associated with biomarkers in complex diseases like prostate cancer (Chikwambi et al., 2023), epilepsy (El Abed et al., 2023), malaria/COVID-19 (Nzungize et al., 2022) and SARS-CoV-2 infection (Nyamari et al., 2023).Several studies suggest that overexpression of CDK1, which is the hallmark event in many malignancies is associated with 5 = fluorouracil resistance through altered activation of Wnt/β-catechin signaling, mTOR, and p53 pathways leading to poor prognosis in CRC (Zhu et al., 2019;Zhu et al., 2020a;Zhu et al., 2020b;Sofi et al., 2022a).
Natural products have drawn a lot of attention recently as potential sources of treatments for a variety of illnesses.In recent years, bioactive compounds from plants have been proven to have opportunities in the pharmaceutical industry for therapeutic purposes (Sasidharan et al., 2011).Resveratrol, a polyphenol has been documented for its ability to reduce or inhibit cell proliferation, reduce apoptosis and suppress angiogenesis and metastasis in cancer cells (Greenlee et al., 2022;Dar et al., 2023).Also, flavopiridol, a flavonoid based substance has undergone clinical trials for its anticancer benefits and as a potent Cyclin-Dependent Kinase (CDK) inhibitor (Deep et al., 2018).By specifically targeting CDK1, cell cycle arrest is initiated, slowing down the growth of tumors.Compared to synthetic alternatives, natural products have fewer adverse effects.
Despite recent advances in diagnostics and treatments for CRC as well as the failure of conventional cancer treatments including chemotherapy and radiotherapy to sustain lower relapse rates, there remains a gap in treatment and increasing incidence of CRC cases.With the vast plethora of computer-aided drug design and combinatorial chemistry technologies available, a potentially wide array of therapeutic options from small molecule compounds that can target significant aspects of CDK1 activity exists.Our study investigated potential inhibitors of CDK1 in silico and explored their pharmacokinetic properties using a number of computational techniques.

Ligand preparation
A library of ten thousand bioactive plant compounds was compiled and downloaded from the PubChem database (https:// pubchem.ncbi.nlm.nih.gov).Pubchem is a publicly available small molecule database (Li et al., 2010;Gozalbes and Pineda-Lucena, 2011;Cheng et al., 2014) used in drug discovery and other applications.It contains information on over 250 million unique chemical molecules along with their descriptors.Before carrying out docking, the two-dimensional (2D) structure data file (SDF) format of the compounds downloaded from the pubchem repository were imported into Maestro (Schrodinger suite) for transformation in the LigPrep module of the software (Epik, 2017a).Ionization states were generated at a pH range of 7 ± 2.0 and stereoisomers were calculated at default (Sastry et al., 2013).

Protein preparation
The three-dimensional (3D) human crystal X-ray diffraction structure of Cyclin-dependent kinase 1 (http://www.rcsb.org/pdb/home; PDB ID: 4Y72) (Brown et al., 2015), was downloaded as a pdb file and imported into Maestro suite.For appropriate structural preparation of the virtual protein using the Protein Preparation wizard of the Schrodinger software, hydrogen, bond orders, and missing side chains were added appropriately while water molecules beyond 5 Å from het groups were deleted before optimization at pH 7.0 + 2.0.The protein structure was subsequently minimized using the optimized potentials for liquid simulations (OPLS) 2005 force field to a root-mean-square-deviation (RMSD) of 0. 30 Å (Epik, 2017b).

Molecular docking
The compounds were subjected to molecular docking within the active site of CDK1.Using the receptor grid generation feature in the Schrodinger suite, a grid box where each ligand will be docked was generated.The binding pocket where the co-crystallized ligand is located was used to construct the defining box as well as the residue locations within 4 Å radius of the ligand.Molecular docking was performed via the standard precision (SP) and extra precision (XP) of the Glide tool in the software (Friesner et al., 2006;Tripathi et al., 2013).Following docking, the top compounds with optimal binding scores were selected for further studies.

Drug likeness and in silico ADME predictions
The drug-likeness of the lead compounds was evaluated by subjecting them to computational pharmacological and pharmacokinetic analyses on the SwissADME web server (http://www.swissadme.ch)(Daina et al., 2017;Guan et al., 2019), using the generated Canonical SMILES (simplified molecular-input entry system) transformations of the structures (Cheng et al., 2022).Lipinski's rule of five was used to assess the compounds based on the following parameters of molecular weight, number of hydrogen bond donors and acceptors and octanol-water partition (Lipinski et al., 1997;Veber et al., 2002;Yang et al., 2019).Pro-Tox II online server (https://tox-new.charite.de/protox_II)was used to predict the toxicity of the compounds (Banerjee et al., 2018).

Induced fit docking
The top four ligands with the least binding scores were selected and subjected to the mixed molecular docking and dynamic protocol used by the induced-fit docking (IFD) module of the Maestro algorithm.An implicit solvent model with the OPLS 2005 force field was used to apply the standard IFD protocol to the chosen (centroid) amino acid side chains (8-10, 18-20, 31-33, 64, 80-89, 135, 145-146).The IFD procedure included a ring conformational sampling with a 2.5 kcal/mol energy barrier and a non-planar conformation penalty on amide bonds.Each ligand was given a maximum of 20 permitted postures, with the scaling for both the receptor and the ligand fixed at 0.5.The Prime Refinement technique was used to further refine residues that were within 5 Å of the docked ligand.The protein-ligand complexes with the improved structure were ranked using prime energy.For a final round of Glide docking and scoring, the receptor structures that were less than 30 kcal/mol from the minimal energy structure were submitted.In the next second docking stage, all ligands were redocked into all revised low-energy receptor structures using Glide XP's default settings.

Molecular dynamics simulation
The docked complexes were further subjected to molecular dynamics (MD) simulations using the Desmond module of the Schrodinger package on a GPU-enabled Linux computer.To simulate physiological conditions, the model was created in a TIP3P (transferable intermolecular potential-three point), orthorhombic water buffered solvation box of 10 Å dimensions and 0.15 M salt.The simulation model was then minimized and neutralized by adding sodium (Na + ) and chloride (Cl − ) ions, placed at least 20 Å from the ligands.The protein-ligand complexes were simulated for one thousand nanoseconds in an OPLS-2005 force field.At the end of the simulation period, the root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), Ligand RMSD, solvent accessible surface area (SASA), and intramolecular hydrogen bonds (intraHB) were computed and subsequently analyzed (Chow et al., 2008;Bergdorf et al., 2014;Release, 2014).

Molecular mechanics/generalized born and surface area solvation (MMGBSA) calculations
Molecular mechanics/generalized Born solvent area (MM/ GBSA) calculations (Kollman et al., 2000;Massova and Kollman, 2000;Tsui and Case, 2000) were performed to estimate the relative binding energies of the selected ligands to CDK1.This method estimates the amount of free energy involved in interactions of the protein-ligand complex.The binding energy for 2000 independent snapshots of each ligandprotein complex was generated and average binding energies calculated.ΔG of each ligand was calculated following the given equation: Where ΔG RL represents the free energy of protein-ligand complex; ΔG R is free energy of receptor, and ΔG L the free energy of ligand.
The free binding energy establishes the binding affinity of ligands to a given protein.A more negative ΔG value indicates better binding of the ligand to the protein (Gilson et al., 1997).The free energy of each state was calculated, using FF14SB force field, by the given equations. . .
( 2 ) Contribution of each residue to the total binding free energy was also obtained using MM/GBSA method in Amber 18.

Analysis of the binding modes of co-lig and hit compounds in CDK1 active site
We performed molecular docking study of the phytocompounds against CDK1 target protein.We selected the top four hits with the most favourable binding scores for subsequent studies.To validate the docking procedure, the co-crystallized ligand was re-docked into the binding pocket of CDK1.The RMSD between docked and native pose (Figure 2) is given as 0.28 Å. Figure 3 depicts the result from molecular docking studies of the four top-scoring compounds: Spiraeoside (−12.47 kcal/mol), Robinetin (−12.22 kcal/mol), 6-hydroxy luteolin (−12.07 kcal/ mol), Quercetagetin (−12.06 kcal/mol), and the co-crystallized ligand, LZ9 (−13.99 kcal/mol).The 2D protein-ligand interaction diagrams (Figure 4) show the hydrogen and hydrophobic bond network and interactions shared between the compounds and CDK1.In this study, it was observed that all 4 hit compounds formed hydrogen bonds and hydrophobic interactions with the Ile10 residue of CDK1.Hydrogen bond interactions also existed between the compounds and CDK1 via Leu83, Ser84, and/or Gln 132 residues.Hydrophobic interactions were seen between all the compounds and Ile10, Val18, Ala31, Val64, Phe80, Phe82, Leu83, Met85, Leu135 and Ala145 residues of CDK1.LZ9 shared hydrophobic bonding with Tyr15, in addition to the aforementioned residues.Hydrophobic bonds observed from molecular docking were preserved during induced fit docking.Bond interactions for molecular docking, induced fit docking and molecular dynamics simulation are summarized in Table 1.The complexes formed from induced-fit docking are displayed in Figure 5.

In silico ADME predictions
Table 2 gives details of ADME and toxicity predictions.Physicochemical descriptors, ADME (Absorption, Distribution, Metabolism, and Excretion) parameters, and pharmacokinetic properties of the compounds, except spiraeoside, were predicted to be within the permissible parameter ranges.
From toxicity prediction, spiraeoside was predicted to be inactive for all parameters.LZ9 was predicted to be hepatotoxic and carcinogenic, but not mutagenic.Robinetin, 6-hydroxyluteolin, and quercetagetin were predicted to be inactive for hepatotoxicity, while results for carcinogenicity and mutagenicity were predicted to be active.Their probability scores are shown in brackets.

RMSD analysis of CDK1 and hit compounds
The protein backbone RMSD (Ca) & ligand heavy atoms RMSD over the course of simulation time were plotted.The average values of each system were also calculated.The RMSD depicts the dynamic behavior and stability of the complexes over the simulation period and estimates the average change in atom displacement for a given frame relative to a reference frame.All protein frames were first aligned on the reference frame backbone before determining the RMSD based on the atom selection.RMSD fluctuations give insight to the stability of a ligand within the binding pocket of a protein.A gradual increase in RMSD values was observed in the first 200 ns of the apo protein system, after which it stabilised and had an average value of 2.48 Å (Figure 6A).The CDK1-LZ9 system fluctuated for the first 400 ns, after which it stabilised and had an average of 2.36 Å. CDK1-robinetin stystem stabilised after 200 ns with an average value of 2.40 Å, while the CDK1-quercetagetin system fluctuated more, however, had an average value of 2.24 Å.The most fluctuations were observed in the CDK1-6-hydroxyluteolin system which also had the highest mean value of 3.17 Å.From the ligand RMSD plot (Figure 6B), the ligands displayed quite stable trends over the simulation period, except for quercetagetin, which fluctuated significantly after 400 ns, even though it had an average RMSD value of 4.48 Å. LZ9 had an average RMSD value of 1.27 Å, robinetin, 5.10 Å and 6-hydroxyluteolin, 7.00 Å.

RMSF analysis of CDK1 and hit compounds
In general, secondary structure components of a protein like alpha helices and beta strands are more rigid than the unstructured part and fluctuate less than loop regions during molecular dynamics simulation runs.Analysis of the RMSF plots (Figure 6C), show that all the systems displayed similar patterns.Notable distinctions include quercetagetin with the highest peak value at residue 100, robinetin also displayed highest peak at residue 130, whereas LZ9 had highest peak value at residue 240.LZ9 had the least peak at residue 165 while robinetin had the least value between residues 230 and 240 where other system had highest peaks.A comparative 3D representation of the crucial residue regions of the systems which exhibited higher fluctuations has equally been presented.

Radius of gyration, solvent accessible surface area and intra-hydrogen bond analysis
The radius of gyration for all systems is shown in Figure 6D.When comparing the plots of the apo form to the CDK1 bound states, there is no apparent alteration to the global protein compactness in response to ligand binding.The solventaccessible surface area describes the surface area of a molecule accessible to water.In addition to examining the impacts on the molecule, it expresses the degree to which a simulated protein interacts with solvent molecules.Evaluating the SASA value is an important way to estimate conformational changes caused by complex interactions.This is essential for characterizing proteinligand interactions and understanding structural changes during simulations.Solvent Accessible Surface Area (SASA) plots (Figure 7A) show that LZ9 and quercetagetin have high SASA values, while robinetin and 6-hydroxyluteolin have lower SASA values.Intramolecular hydrogen bonds show the number of internal hydrogen bonds (HB) within a ligand molecule and give insight into the behavior of ligands in the active site of a protein.They can also stabilize conformations that are conducive to protein binding, thus increasing the affinity for the receptor.Intra-hydrogen bond (intraHB) plots (Figures 7B-D) reveal that LZ9 had intraHB interactions almost throughout the simulation time, 6hydroxyluteolin and quercetagetin had more intraHB interactions at the beginning of the simulation, while robinetin maintained a value of zero, that is zero intraHB interactions throughout the simulation time.

Protein-ligand contacts
Figure 8 shows a schematic of the detailed ligand atom interactions with the protein residues.The interactions in this trajectory that occur for more than 30% of the simulation period are displayed.LZ9 and robinetin make pi-pi contacts with Tyr15 and Phe80 for more than 50% of the time.All the ligands made contact with Leu83 via hydrogen bonds; LZ9 and robinetin form hydrogen bonds with Glu81, while 6-hydroxy luteolin6-hydroxy luteolin and quercetagetin form hydrogen bonds with Glu51.Robinetin, 6hydroxyluteolin, and quercetagetin form salt bridges with Lys33.All the ligands show relatively high solvent exposure.Protein-ligand contacts as illustrated by bar charts (Figure 9) during the simulation, show that hydrogen bonds, hydrophobic interactions, and water bridges were predominant in LZ9.The other ligands had ionic interactions in addition to the aforementioned.All compounds exhibited strong hydrogen bonding on Leu83 and hydrophobic interactions on Leu135.They also showed a water bridge on residue Asp146.Only the lead compounds showed ionic bonds at residue Lys33.

Potential energy
Figure 10 shows the potential energy plot of the systems in the apo state and in complex with the ligands.When bound to LZ9, robinetin and quercetagetin, we observe energies in the same range as that of the apo state; −115,500 to −117,000 kcal/mol.On the other hand, when bound to 6-hydroxyluteolin, a lower energy state of −144,000 kcal/mol is observed.

Free binding energy
Robinetin and 6-hydroxyluteolin exhibited good ΔG values of −30.40 kcal/mol and −30.63 kcal/mol, respectively (Table 3).The co-lig, LZ9, had a value of −34.96 kcal/mol.Quercetagetin had a binding energy of −18.52 kcal/mol.Energy contributions of active site residues responsible for stabilizing the CDK1-ligand complexes are shown in Figure 11.Worthy of note is Glu81 as the residue showing the highest energy contribution in the CDK1-LZ9 complex, while Asp86 is the residue showing highest binding energy contribution to the CDK1-robinetin and CDK1-6hydroxyluteolin complexes.The CDK1-quercetagetin system has Asp146 as the active site residue contributing the highest energy.Visual inspection of the figure shows that Van der Waals forces contribute more to the interaction between CDKI-LZ9 and CDK1-6-hydroxyluteolin complexes.Robinetin has more electrostatic energy contributing to the stability of the complex it forms with CDK1.Quercetagetin forms weak interaction with CDK1 as a result of the overall low Van der Waals forces and electrostatic energy.

Discussion
The prevalence of colorectal cancer is fast becoming a public health concern (Holowatyj et al., 2020).Development of new pharmacologically therapeutic anti-cancer drug agents has been hampered by the limited molecular understanding of the oncogenetic initiation and the interaction of biochemical mechanisms/factors, although traditional chemotherapy options have produced impressive results in the treatment of CRC with the possibility of relapses.
Natural products are a valuable source of novel therapeutic leads given that they contain a wide variety of chemical components (Cragg and Newman, 2013;Newman and Cragg, 2016).They have been essential in the development of new drugs, particularly for the treatment of cancerous disorders (Atanasov et al., 2015;Harvey et al., 2015).The compounds identified from this study were found to be in the class of flavonoids.There is growing evidence that flavonoids may be useful in the prevention and treatment of cancer and they are also reported to have a variety of anti-cancer effects, including inducing apoptosis in cancer cells, inhibiting cell proliferation, and may enhance chemotherapy (Abotaleb et al., 2018).Flavonoids have also been reported to play a potential role in colorectal cancer therapies (Amado et al., 2014).
Given the importance of cyclin-dependent kinases (CDK) in cancer, targeting CDKs is a viable strategy for combating the cancer menace (Sofi et al., 2022a).Because of the significant role it plays in one of the phases of cell cycle control, which is a crucial one for the growth of tumor cells, specifically targeting CDK1 became a strategic focus for optimized treatment.The drug discovery pipeline greatly benefits from the use of computational screening techniques (Choudhury et al., 2022).In recent years, computer-aided drug development methods have been widely used to find lead compounds with therapeutic value.These technologies include molecular docking, in silico pharmacokinetic profiling, and molecular dynamics simulations (MDS) (Wan, 2013;Salmaso and Moro, 2018).MDS are anticipated to become more significant as new pharmacological therapies are developed (Durrant and   McCammon , 2011).Several studies now use these computational methods to quickly identify potential inhibitors of key proteins in disease progression as treatment options (Choudhary et al., 2020;Sofi et al., 2022a;Rathod et al., 2022;Castrosanto et al., 2023) others have identified experimentally validated leads (Kalani et al., 2013;Yilmaz Goler et al., 2022).Here, we used these computational techniques to screen phyto-compounds and identify possible inhibitors of CDK1.Molecular docking is a computational method used to predict the preferred binding orientation of a ligand to a receptor.The software algorithm used predicts possible binding orientations and provides scores ranked according to the predicted binding affinity between the receptor-ligand complex (Huang and Zou, 2010;Salmaso and Moro, 2018).To better understand the binding dynamics of the four compounds selected in this study to CDK1, they were subjected to molecular dynamics simulations (MDS).Prior to running simulation studies, we carried out induced fit docking (IFD) (Achilonu et al., 2020).Induced fit docking is a more realistic model of how proteins and ligands interact in nature as it considers the flexibility of the protein and ligand, allowing them to  move and rotate during the docking process, resulting in more accurate predictions of the binding of ligands.For traditional rigid docking, the proteins and ligands are treated as rigid bodies, and the docking algorithm attempts to find the best orientation of the ligand in the protein's binding site.However, in reality, proteins are flexible molecules, which can undergo conformational changes upon ligand binding.Both minor backbone relaxations and substantial sidechain conformational changes in the receptor structure are accounted for by the induced-fit docking (IFD) technique (Sherman et al., 2006;Xu and Lill, 2013).All systems ran through the complete one thousand nanoseconds duration during MDS, except spiraeoside, which ran for only 700 ns?It could be that due to conformational changes in the protein, spiraeoside left the binding site of the protein during the simulation run.For this reason, the post-simulation analysis could not be completed, consequently, results for spiraeoside were excluded.
It is important to pay attention to the protein's RMSD value to properly infer structural changes occurring in a protein which could affect ligand binding.Ca-RMSD demonstrated a stable structural conformation across the course of the simulation, showing that it was balanced.The low RMSD values of LZ9, robinetin and quercetagetin indicate low fluctuations and suggests that the system maintained stability during simulation and the protein did not undergo significant structural changes as a result of binding of the ligands.The higher RMSD value of 6hydroxyluteolin system might suggest that the binding of this ligand might have perturbed the structure of the CDK1 protein.
Ligand RMSD provides insight into how stably the bound ligand sits within the protein binding pocket.The Ligand RMSD results show that LZ9, robinetin, and 6-hydroxyluteolin did not move much, therefore forming a stable complex with CDK1, unlike quercetagetin which had a lot of fluctuations.The fluctuations observed with quercetagetin after 400 ns tell us that it is not tightly bound to the protein and has possibly diffused away from the binding pocket, so it forms an unstable protein-ligand complex and is less likely to perform its intended function.
RMSF measures the structural flexibility of a protein and characterizes local perturbations along a protein chain (Cob-Calan et al., 2019).It also provides information on how each amino acid residue is changing relative to its reference position during simulation (Sarkar and Sen, 2022).Here, we analysed the RMSF of the CDK1 protein in the apo and bound states.Peaks represent regions of the protein where the simulation-induced fluctuations were greatest (Tripathi et al., 2022).Illustrative 3D images of the peaks have been presented to clearly elucidate the silent fluctuation dynamics depicted by the represented peaks.The 3D images reveal a comparative conformational variation exhibited by the various system due to the perturbations created by the MD simulation.The 3D images have been color-coded to tally with their respective corresponding RMSF trajectories.The spikes observed from RMSF results may have been caused by increased flexibility within the protein's loop regions.It may be inferred from the RMSF result that CDK1 has a great deal of flexibility to accommodate the ligand at the binding pocket.
The radius of gyration measures how compact a protein structure is and is equivalent to its principal moment of inertia, describing the common separation of each scattered element from the molecule's mass center (Lobanov et al., 2008).During simulation, a more rigid structure is represented by lower Rg values (Jiang et al., 2019).Since not much difference was observed between the Apo and bound states, we can say that the binding of the ligands did not interfere with the compact structure of CDK1.In this study, the Rg, RMSD, and RMSF were correlated, and this relationship shed light on how complex compressibility affected the variations in complex RMSD and the residual fluctuations seen in RMSF.No noticeable changes are observed in the bound and unbound states of CDK1, suggesting that the protein structure retained its compact conformation.Our results indicate that the surface area of LZ9 and quercetagetin have higher exposure to water than robinetin and 6-hydroxyluteolin.This result is in tandem with that observed from the ligand-atom interaction diagram.
Hydrogen and hydrophobic bonds play key roles in proteinligand interactions.Crucial hydrophobic and/or hydrogen bonds, including Ile10, Phe80, Leu83, and Ser84 between CDK1 and the lead compounds that were first seen after molecular docking were maintained during MDS.Phe80 has previously been described as a hydrophobic checkpoint enabling water molecules to gain access to the catalytic portion of the ATP catalytic domain for an additional hydrogen bond (McGrath et al., 2005).Robinetin is seen to form a pi-pi bond with Phe80 during simulation run, which may further account for its strong inhibitory properties.6-hydroxyluteolin also forms a bond with Thr14 for 88% of the simulation time.Inhibition of Thr14 by small molecule compounds has been proposed as a promising option for the regulation of the G2/M transition of the cell cycle (Schmidt et al., 2017).Inhibition of CDK1 due to the phosphorylation of the kinase on Tyr15 residue by flavonoids (Casagrande and Darbon, 2001) as well as induction of G 2 /M phase arrest and apoptosis by the flavonoid tamarixetin on human leukemia cells (Nicolini et al., 2014) have previously been reported.
The hydrophobic and hydrogen bonds formed between the compounds and CDK1 may play a vital role in stabilizing the complexes for the compounds to exert their inhibitory efficacies (Zhang et al., 2011).Ionic bonds observed on only the lead compounds might also play a significant role in their inhibition of CDK1.Interestingly, the lead compounds are all of the flavonoid class; another study by ((Navarro-Retamal and Caballero, 2016)) also gave insight into flavonoids as CDK1 inhibitors.Natural flavonoids such as Luteolin and Daidzein have also been found to inhibit CDK1 (Ravishankar et al., 2013).
Higher free energy values in fact correlate with more preferable binding affinities of a ligand to a protein.Robinetin, which had good binding energy to CDK1, as revealed by the ΔG values, might have been a result of higher electrostatic and van der Waals energy contributions.Interactions between robinetin and Asp86 residue which was observed from molecular docking and induced fit docking results was maintained after the MM/GBSA run.We can say this interaction might be a key factor to stablizing the CDK1-robinetin complex.Energetic analysis reveals that van der Waals interaction and non-polar contributions are favorable in the formation of complexes and amine group of the ligand, which plays a crucial role for binding (Tripathi et al., 2012;He et al., 2023).Potential energy (PE) and RMSD are both vital indicators of protein stability.In previous work done by Yuan et al. (2012) to ascertain if the structures were stabilized, MD simulation was examined using PE and RMSD from the original model structure (Coutsias et al., 2004) The potential energy plot for 6hydroxyluteolin is a reflection of the trajectory pattern exhibited by the Cα-RMSD plot of 6-hydroxyluteolin system between 230 ns and 420 ns, where a perturbation on the 6-hydroxyluteolin system seems to have created conformational alteration within the system, making it assume a new/different trajectory.The new structural conformational deviation of 6-hydroxyluteolin relative to other systems could presumably have produced the obtained PE.
The physicochemical properties of a compound are important parameters in the pipeline of drug development as they directly affect the absorption, distribution, metabolism, excretion, and toxicity profile (ADMET) of a drug candidate (Manallack et al., 2013).The drug-like and ADMET predictions are used in identifying properties of a ligand suitable for formulation into drugs.Even though all the compounds are adequately hydrophilic and lipophilic, robinetin and 6-hydroxyluteolin showed high intestinal absorption.This implies that the compounds are capable of achieving pharmaceutical and therapeutic bioavailability at the intestinal linings and target cells in the colon.This is further corroborated by their bioavailability scores.A study on natural compounds as inhibitors of CDK1 by Saikat (Saikat, 2021) also reported similar ADME properties.
Lipinski's Rule of five is used to assess and predict the druglikeness of compounds based on oral absorption, membrane permeability and bioavailability (Lipinski et al., 1997).Although, parenteral administration may be applicable in most cases of cancer treatment, oral administration is a valuable treatment option for cancer patients as it is more convenient and reduces the frequency of their clinic visits or venipuncture.In recent times, oral anti-cancer drugs such as capecitabine, topotecan and vinorelbin are successfully used in the treatment of a variety of tumours like colon, breast, ovarian or lung cancer (Colomer et al., 2010).With the advent of targeted deliveries, most of the drug molecules approved as oral medications for cancer treatment including hormone therapies such as, dasatinib, sorafenib, imatinib, lapatinib and abiraterone have resulted in good clinical outcomes (Gralow et al., 2008;Rehman and Rosenberg, 2012).
Despite the fact that natural products frequently violate these rules due to their complexity, they still serve as reliable sources for novel drug discovery and development (Doak et al., 2014).Some anti-cancer agents such as vincritine, irinotecan, paclitaxel and etoposide, widely used for decades in cancer treatments originate from natural compounds (Huang et al., 2021).
The ProTox-II server is a web-based tool that predicts the potential toxicity of a given small molecule compound on various organs and systems in the body.The toxicity is predicted using a machine learning model that has been trained on a sizable dataset of toxic and non-toxic molecules.The toxicity prediction and probability score are considered important parameters for predicting a compound's toxicity (Han et al., 2019).The probability score can be used to evaluate the prediction's level of confidence.A forecast with a high probability score is more certain to occur, whereas one with a low probability value may not.Additionally, if the organ toxicity and/or toxicity end points of a compound are classified as active, it suggests that the compound may have harmful effects if it is ingested, inhaled, or otherwise exposed to the body.Conversely, "inactive" could mean that the compound has not been predicted to have a toxic effect.Early identification of toxicity in the drug development steps is critical in reducing late-stage attrition rates (Roberts et al., 2014;Segall and Barber, 2014).The predicted toxicity class exhibited by the compounds indicates that they can be administered within the accepted safety doses.Although the test ligands are predicted to be non-hepatotoxic, lead optimization may be required to modify the other toxicity endpoints.

Conclusion and next steps
In this study, computational techniques including molecular docking, induced fit docking, molecular dynamics simulation and MM/GBSA were used to study the potential inhibitory activities of natural compounds against CDK1.The findings provided insight into molecular interactions of significant residues involved in inhibitory CDK1 activities through molecular docking and molecular dynamics simulation results.The current investigation of small molecule inhibitors in the activity of the CDK1 explored viable options for the roles of Robinetin and 6-hydroxyluteolin in the pathogenesis of CRC.This research, however, provides a basis for further examination of potential phytochemicals for medical intervention of several cancers.This study aimed to identify potential small molecule inhibitors of CDK1 hence the use of molecular modeling approaches was deployed to highlight compounds of strong binding energies, good spatial amino acid residue arrangement, and better pharmacokinetic properties among others.The MD simulations showed that the hit compounds stably interacted with the CDK1 enzyme and maintained firm positions within the binding pocket of the enzyme over the simulation timeframe.The findings of this study should undergo experimental validation to confirm the interaction between the top ligands identified in our research and the target protein.
Additionally, the ADME/Tox properties that were predicted using computational models need to be verified through experimental testing.curation, Formal Analysis, Methodology, Software, writing-review and editing.IA: Data curation, Formal Analysis, Methodology, Project administration, Software, Writing-review and editing.OA: Conceptualization, Data curation, Formal Analysis, Methodology, Project administration, Software, Supervision, Writing-review and editing.

Funding
The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

FIGURE 2
FIGURE 2Validation of docking showing the co-crystallized ligand (in grey) and the re-docked pose (in red) with an RMSD score of 0.28 Å.The low RMSD value obtained demonstrates the efficacy of the docking protocol used.

FIGURE 4
FIGURE 42D interaction diagrams of top scoring compounds in complex with CDK1 from molecular docking studies.Maestro 2D interaction diagram (Schrödinger) was used to generate the image.

FIGURE 6
FIGURE 6 Molecular dynamics simulation trajectories of CDK1 in presence of ligands and absence of ligand (Apo) as a function of 1,000 ns simulation time.(A) Root Mean Square Deviation (RMSD) of Cα atoms backbone.Insert, # shows the mean Cα RMSD values.(B) Root Mean Square Deviation of the ligands with respect to CDK1 protein.The plot displays the RMSD of the ligands after measuring the RMSD of the heavy atoms in the ligand and aligning the protein-ligand complex on the reference protein backbone.If the reported values are noticeably greater than the protein's RMSD, the ligand has probably diffused away from its initial binding site.Insert, * shows the mean RMSD values.(C): Root Mean Square Fluctuation (RMSF) per residue over time.(D): Radius of gyration for the Apo state and CDK1-ligand complexes.

FIGURE 8 A
FIGURE 8A schematic of detailed ligand atom interactions of LZ9 (A), Robinetin (B), 6-hydroxyluteolin (C) and Quercetagetin (D) with the protein residues of CDK1.Interactions that occur more than 30.0% of the simulation time in the selected trajectory (0 through 1000 ns) are shown.

FIGURE 9
FIGURE 9 Protein-Ligand Contacts.(A) LZ9.(B) Robinetin.(C) 6-hydroxyluteolin.(D) Quercetagetin.Contact interactions can be categorized by type and summarized, as shown in the plot.Protein-ligand contacts are categorized into four types: Hydrogen bonds, Hydrophobic, Ionic and Water Bridges.The stacked bar charts are normalized over the course of the trajectory: for example, a value of 0.7 suggests that 70% of the simulation time the specific interaction is maintained.Values over 1.0 are possible as some protein residue may make multiple contacts of same subtype with the ligand.

FIGURE 10
FIGURE 10Potential energy plot of 1,000 ns simulation time.The plots show the energies of the systems in apo and bound states.

TABLE 1
Bond interactions between ligands and CDK1 during molecular docking, induced fit docking and molecular dynamics simulation.
MD, molecular docking; IFD, induced fit docking; MDS, molecular dynamics simulation.FIGURE 52D interaction diagrams of poses from Induced fit docking.

TABLE 2
Predicted physicochemical descriptors, ADME parameters, pharmacokinetic properties, druglike nature, medicinal chemistry friendliness and toxicity of compounds.