Exploring the Activity of Fungal Phenalenone Derivatives as Potential CK2 Inhibitors Using Computational Methods

Cancer represents one of the most prevalent causes of global death. CK2 (casein kinase 2) activation boosted cancer proliferation and progression. Therefore, CK2 inhibition can have a crucial role in prohibiting cancer progression and enhancing apoptosis. Fungi have gained vast interest as a wealthy pool of anticancer metabolites that could particularly target various cancer progression-linked signaling pathways. Phenalenones are a unique class of secondary metabolites that possess diverse bioactivities. In the current work, the CK2 inhibitory capacity of 33 fungal phenalenones was explored using computational studies. After evaluating the usefulness of the compounds as enzyme inhibitors by ADMET prediction, the compounds were prepared for molecular docking in the CK2-α1 crystal structure (PDB: 7BU4). Molecular dynamic simulation was performed on the top two scoring compounds to evaluate their binding affinity and protein stability through a simulated physiological environment. Compound 19 had a superior binding affinity to the co-crystallized ligand (Y49). The improved affinity can be attributed to the fact that the aliphatic chain makes additional contact with Asp120 in a pocket distant from the active site.


Introduction
Cancer is a complicated illness that is featured by uncontrolled cell proliferation and the development of tumors that may remarkably extend to the entire organ or other organs systemically [1]. It represents one of the most global deadly diseases, particularly in western countries. It was estimated that in 2020,~9.9 million people have died due to cancer [2]. Its current treatment strategies include γ-radiation, chemotherapy, suicide gene therapy, or immunotherapy, which are fundamentally mediated by promoting apoptosis [3]. Chemotherapy is the most efficacious method for metastatic tumors treatment. However, the cancer cell's multidrug resistance, high cost, and untoward effects of these drugs represent the main holdbacks to the success of chemotherapeutic treatment [4]. Therefore, searching for sources of safe and effective treatment has become a crucial research area.
Casein kinase 2 (CK2) is one of the first discovered Ser/Thr kinases [5][6][7] that is involved in many cell processes from gene expression and protein synthesis to cell growth, proliferation, and differentiation [8]. The well-studied tetrameric form of this kinase is Computer-aided drug design/discovery (CADD) are useful tools for screening and identifying drug-like molecules in silico and thereby reducing the number of compounds to be tested experimentally. Several software and programs are used to filter and generate a group of compounds based on specified criteria, predict their physicochemical properties, predict suitable targets, and evaluate the binding affinity of the compounds to the predicted targets. One of the drug design and discovery approaches is structure-based drug design (SBDD). This approach relies on the knowledge of the 3D structure of the targets of interest, and it includes two common methods: molecular docking and molecular dynamic simulation. Molecular docking evaluates how tight the compound binds the target, as determined by the predicted binding affinity, while molecular dynamic (MD) simulations assess the behavior of the ligand-protein complex in terms of binding interactions and 3D conformation in aqueous conditions to mimic the physiological environment [43]. Several CADD methods and tools are used for investigating the phenalenone derivatives. Table 1. List of phenalenone derivatives and their fungal source.

Compound Name
Fungal Source Ref.

. AI (Artificial Intelligence)-Based Target Prediction for Phenalenone Derivatives
Choosing the appropriate target to investigate the inhibitory potential of these phenalenone derivatives relied on performing ligand-based-in silico target prediction [38]. The prediction webserver, SuperPred was the tool of choice to perform ATC (anatomicaltherapeutic chemical) code and predict the potential targets for the investigated phenalenones [57]. After analyzing the results of the predicted target proteins (for example, Cathepsin D, mineralocorticoid receptor, and thyroid hormone receptor-α), the kinase CK2α was selected for the study due to having a high percent of probability and model accuracy ( Table 2). After selecting the target and proper crystal structure, the listed phenalenones were docked in the protein crystal structure, after which the docking method was validated by redocking the co-crystallized ligand. Prediction of ADMET properties of the listed metabolites in silico and MD (molecular-dynamic) simulation for the two metabolites with the highest docking scores were followed.

Ligands and Protein Preparation and Molecular Docking
Target identification filtered out > 100 phenalenone derivatives. These derivatives were prepared for docking, where their energy-minimized 3D structures were generated, and all possible ionization and tautomeric states were created.
For docking, the human CK2α1 crystal structure (PDB ID: 7BU4) was chosen for the study due to the structural similarity of the co-crystallized ligand (Y49) to the selected phenalenones. The Y49 (4-(6-aminocarbonyl-8-oxidanylidene-9-phenyl-7H-purin-2yl)benzoic acid) is made of three aromatic rings: a purine ring with two phenyl moieties attached to it, and polar groups present at the rings, with polar groups at different positions. The selected phenalenones have a nucleus made of three-fused rings substituted with multiple OH and carbonyl groups. Based on their structures, both the phenalenones and Y49 were expected to have a similar 3D conformation in the binding pocket. The PDB file of the 7BU4 crystal structure was downloaded from the protein data bank (PDB) [58], which was then prepared and minimized using Schrodinger's Protein Preparation wizard [59][60][61].
The docking process started by generating a grid box around the binding site of the co-crystallized ligand to locate the pocket where the docking of the compounds occurs. A receptor-Grid-Generation tool in Maestro [62] was utilized for that purpose.
Re-docking of the co-crystallized ligand, 4-(6-aminocarbonyl-8-oxidanylidene-9-phenyl-7H-purin-2-yl)benzoic acid (PDB: Y49), was performed to evaluate the docking method. The re-docked reference had an identical pose ( Figure 3C) and binding interactions to the co-crystallized structure ( Figure 3A,B). For both, the backbone of Val116 is H-bonded with the oxygen of the amide moiety, and the purinone nitrogen interacts through the water bridge with Val116, Asn118, as well as with the amide oxygen. On the opposite side of the molecule, the carboxylate makes an ionic interaction with the side chain of the adjacent Lys68 ( Figure 3A,B). The RMSD (root-mean-square deviation) of the re-docked ligand was minimal, with a value of 0.0744 Å, indicating the docking method is valid ( Figure 3C). The molecular surface display in Figure 3D shows the re-docked reference Y49 occupying the binding pocket of the crystal structure.
After docking validation, docking the 3D structures of the > 100 phenalenones that proceeded from the target prediction using the extra-precision (XP) mode was followed. The docking produced derivatives that are ranked based on their score and approximated the free energy of binding; the more negative the value, the stronger the binding. Different docking scores were generated, including the gscore (best for ranking different compounds), emodel (best for ranking conformers), and XP gscore. Glide uses emodel scoring to select the best poses of the docked compounds; then, it ranks the best poses based on the given gscores. XP gscore ranks the poses generated by the Glide XP mode. In general, Glide uses gscore to sort and rank the docked compounds. The 33 derivatives listed in Table 3 are the ones with gscores close to or better than the gscore of the reference Y49 (−9.049 kcal/mol), with the top two compounds, 19 and 31, scoring −12.878 and −12.521 kcal/mol, respectively.
Compound 19, in addition to interacting directly with Val116 and Lys68 in the protein's binding pocket like the reference ligand, had a long chain that extended along the surface of the protein allowing the terminal (R,6E,10E,14E)-2,6,10,14-tetramethylhexadeca-6,10,14triene-2,3-diol group to bind a distant binding pocket ( Figure 4). Besides, compound 31 seems to have similar interactions with the protein; however, the tetrahydropyran at the end of the aliphatic chain did not occupy the distant pocket like 19 and remained exposed to the solvent ( Figure 5). Figure 6 showed compounds 19 and 31 simultaneously superimposed on the reference Y49 inside the binding pocket in the molecular surface display. After docking validation, docking the 3D structures of the > 100 phenalenones that proceeded from the target prediction using the extra-precision (XP) mode was followed. The docking produced derivatives that are ranked based on their score and approximated the free energy of binding; the more negative the value, the stronger the binding. Different docking scores were generated, including the gscore (best for ranking different compounds), emodel (best for ranking conformers), and XP gscore. Glide uses emodel scoring to select the best poses of the docked compounds; then, it ranks the best poses based on the given gscores. XP gscore ranks the poses generated by the Glide XP mode. In general, Glide uses gscore to sort and rank the docked compounds. The 33 derivatives listed in Table 3 are the ones with gscores close to or better than the gscore of the reference Y49 (−9.049 kcal/mol), with the top two compounds, 19 and 31, scoring −12.878 and −12.521 kcal/mol, respectively.

In Silico ADMET Properties of Selected Ligands
The drug-likeness and ADMET properties of the processed compounds were predicted using Maestro's QikProp Schrodinger module in terms of absorption, distribution, metabolism, excretion, and toxicity, among others [63]. The module can quickly and reliably predict many physicochemical properties and other descriptors, such as the number of possible metabolites and number of reactive functional groups, in order to detect and filter compounds that can be problematic during the late stages of drug discovery and development. Therefore, it can eliminate unnecessary tests and experiments that will ultimately fail in clinical trials [64]. The ADMET prediction evaluates the usefulness of the examined compounds by describing and determining their drug-likeness, physicochemical properties, and expected toxicity profiles. Several descriptors were predicted for these derivatives, and most of the predicted values of ADMET descriptors fell within the recommended range. The predicted ADMET properties and descriptors are presented in Table 4.

In Silico ADMET Properties of Selected Ligands
The drug-likeness and ADMET properties of the processed compounds were predicted using Maestro's QikProp Schrodinger module in terms of absorption, distribution, metabolism, excretion, and toxicity, among others [63]. The module can quickly and reliably predict many physicochemical properties and other descriptors, such as the number of possible metabolites and number of reactive functional groups, in order to detect and filter compounds that can be problematic during the late stages of drug discovery and development. Therefore, it can eliminate unnecessary tests and experiments that will ultimately fail in clinical trials [64]. The ADMET prediction evaluates the usefulness of the examined compounds by describing and determining their drug-likeness, physicochemical properties, and expected toxicity profiles. Several descriptors were predicted for these derivatives, and most of the predicted values of ADMET descriptors fell within the recommended range. The predicted ADMET properties and descriptors are presented in Table 4.

Molecular Dynamic (MD) Simulation
The MD simulations are performed using Desmond software [65,66] to simulate the aqueous physiological environment and assess the changes in protein conformation and binding affinity during the simulation time compared to the original affinity and confirmation of the crystal structure [67]. Only the two top-scoring compounds from the docking study, i.e., compounds 19 and 31 along with reference Y49, were analyzed by MD. The root mean square deviation (RMSD) is a calculated value that compares the poses of investigated compounds to that of the co-crystalized ligand [43]. RMSD plots of the selected compounds complexed with the CK2α1 measure the average change in the positions of the atoms of the protein and ligand inside the binding pocket at the end of the simulation period (100 ns) compared to their starting positions before the simulation at 0 ns. The RMSD plot of Y49 showed that both the protein and the reference Y49 were stable, and the observed fluctuations were insignificant since they were within the acceptable range of 1-3 Å ( Figure 7A). For compound 19, the RMSD of the protein and 19 laid over each other, indicating increased binding affinity of 19 to the protein and stability of the CK2α1-19 complex. Additionally, the fluctuation seen for both over the 100 ns was within the range as well ( Figure 8A). A similar RMSD pattern was observed for 31 and CK2α1 complex, despite the sudden, non-significant fluctuation of 31 at around 90 ns, which is potentially a result of the compound adjusting its pose in the pocket ( Figure 9A). When calculating the RMSD for the compounds, it is not uncommon to observe fluctuation in the plot for some time at the beginning of the simulation, as observed in Figure 7A, Figure 8A, and Figure 9A within the first 20 ns of the run. This expected fluctuation happens as the compound keeps adjusting its conformation inside the pocket to assume a pose that has the least free energy.
The secondary structure of the CK2α1 protein (PDB ID: 7BU4) was also evaluated throughout the simulation while it was complexed with each ligand. Figure 7B represented the protein evaluation while it was complexed with ligand Y49. The top plot showed the distribution of the SSE (α-helices and β-sheets) with the protein represented by the residue index. The middle plot summarized the SSE composition for each trajectory frame throughout the simulation, while the bottom plot monitored each residue and its SSE assignment over the simulation time. Both plots indicated that the overall %SSE of the protein was maintained, and each SSE was stable over the course of the simulation. Comparable results were obtained when 19 ( Figure 8B) and 31 ( Figure 9B) were complexed with the protein.
The MD study also evaluated the binding interactions of a protein-ligand complex. For ligand Y49, the interactions between the Y49 and protein are presented in Figure 10A; the interaction types are color-coded. The stacked bar chart is normalized over the course of the trajectory: for example, a value of 0.8 suggested that the specific interaction was maintained during 80% of the simulation time. Values over 1.0 are possible and indicate that some protein residue may make multiple interactions of the same subtype with the ligand. As indicated in Figure 10A, Val116 made direct H-bonding as well as through water bridges with Y49 and had a normalized value of~1.9. The value >1 represented the combined value of >1 type of interaction, and it indicated that these interactions were maintained for~190% of the simulation time. The other key interactions were with Glu114, Asn118, Lys68, and Asp175, having values of~0.9,~0.7,~1.1, and~1.5, respectively. Figure 10B showed only the interactions between Y49 and the protein that occurred ≥ 30% of the simulation time. Figure 10C displayed the total specific interactions between ligand Y49 and the protein (top plot), whilst the bottom panel demonstrated the protein residues that interacted with the ligand at each time point. As mentioned earlier, Val116 made > 1 interaction with the ligand, which was represented by the dark orange color in the plot throughout the trajectory. Other residues: Lys86, Glu114, and Asp175, also made specific interactions with the ligand. range as well ( Figure 8A). A similar RMSD pattern was observed for 31 and CK2α1 complex, despite the sudden, non-significant fluctuation of 31 at around 90 ns, which is potentially a result of the compound adjusting its pose in the pocket ( Figure 9A). When calculating the RMSD for the compounds, it is not uncommon to observe fluctuation in the plot for some time at the beginning of the simulation, as observed in Figures 7A, 8A, and 9A within the first 20 ns of the run. This expected fluctuation happens as the compound keeps adjusting its conformation inside the pocket to assume a pose that has the least free energy.       Figure 11 shows the amino acid residues of the protein binding pocket that interacted with compound 19. The fused ring system of 19 made similar interactions with the pocket residues as Y49, where Val116 and Lys86 interacted through the H-bond with the ring system ( Figure 11B). As previously seen in the molecular surface display (Figure 4A), the extended aliphatic chain occupied a distant pocket and created new interaction points between the two OH groups at the end of the chain and Asp120, where they occurred > 85% during the simulation ( Figure 11B). This interaction was not present with Y49 ( Figure 3B,D). It might be safe to assume that the enhanced binding affinity and stability of the complex were due to this new interaction with Asp120, which can also be inferred from the RMSD plot ( Figure 8A).
Compound 31 also created new interacting points with amino acids in the main binding pocket: the chiral OH of the fused ring system made a strong H-bonding with His160 and Asn161 through bridging water molecules with a value of~1.18 and 0.6, respectively. An enhanced interaction with Arg47 was observed as well (~0.98) ( Figure 12B). The tetrahydropyran ring at the end of the aliphatic chain of 31 also extended along the protein surface but did not occupy the distant pocket, as did compound 19 ( Figure 5A). The reason for that might be the fact that the chain in compound 31 is 6-carbon shorter than that of 19, so the group could not reach the second pocket. Another explanation could be the large size of the substituted tetrahydropyran hindered the group from occupying the pocket. Additionally, there is a high probability that the fluctuation observed in the RMSD of 31 towards the end of the simulation time ( Figure 9A) might be a result of the inability of this group to bind to the second pocket. The L-RMSF (ligand-root-mean-square fluctuation) represents how the atoms of the ligand interact with the protein and the changes in the positions of the ligand atoms. As seen in the L-RMSF plot for compound 31 (Figure 13), the positions of atoms 29-45 were dramatically changed because of the free rotation around the aliphatic chain, which in turn decreased the interaction between this part of the molecule with the protein and was reflected by > 3 Å fluctuation in the RMSF plot. The time-depended representation of the CK2α1-31 interactions showed that residues Arg47, His160, and Asp175 were the ones making specific interactions with the ligand, as indicated by the darker color in the plot ( Figure 12C). The secondary structure of the CK2α1 protein (PDB ID: 7BU4) was also evaluated throughout the simulation while it was complexed with each ligand. Figure 7B represented the protein evaluation while it was complexed with ligand Y49. The top plot showed the distribution of the SSE (α-helices and β-sheets) with the protein represented by the residue index. The middle plot summarized the SSE composition for each trajectory frame throughout the simulation, while the bottom plot monitored each residue and its SSE assignment over the simulation time. Both plots indicated that the overall %SSE of the protein was maintained, and each SSE was stable over the course of the simulation. Comparable results were obtained when 19 ( Figure 8B) and 31 ( Figure 9B) were complexed with the protein.
The MD study also evaluated the binding interactions of a protein-ligand complex. For ligand Y49, the interactions between the Y49 and protein are presented in Figure 10A; the interaction types are color-coded. The stacked bar chart is normalized over the course of the trajectory: for example, a value of 0.8 suggested that the specific interaction was maintained during 80% of the simulation time. Values over 1.0 are possible and indicate that some protein residue may make multiple interactions of the same subtype with the ligand. As indicated in Figure 10A, Val116 made direct H-bonding as well as through water bridges with Y49 and had a normalized value of ~1.9. The value >1 represented the combined value of >1 type of interaction, and it indicated that these interactions were maintained for ~190% of the simulation time. The other key interactions were with Glu114, Asn118, Lys68, and Asp175, having values of ~0.9, ~0.7, ~1.1, and ~1.5, respectively. Figure  10B showed only the interactions between Y49 and the protein that occurred ≥ 30% of the simulation time. Figure 10C displayed the total specific interactions between ligand Y49 and the protein (top plot), whilst the bottom panel demonstrated the protein residues that interacted with the ligand at each time point. As mentioned earlier, Val116 made > 1 interaction with the ligand, which was represented by the dark orange color in the plot throughout the trajectory. Other residues: Lys86, Glu114, and Asp175, also made specific interactions with the ligand. It was also noticed that the reference ligand, Y49, as well as the compounds 19 and 31, all interacted with Asn118 through H-bonding; however, the contact points are different. While the residue interacted with Y49 at the purine carbonyl oxygen ( Figure 10B), it acted as an H-bond donor to the OH at the end of the aliphatic chain of 19 as well as to the OH at the fused ring system of 19 and 31 (Figures 11B and 12B). The reason for the different contact points of Asn118 with the compounds was probably attributed to the 3D conformation of the compounds inside the binding pocket, which is affected by the nature of the substitutions on the nucleus. Each compound assumed a pose that had the lowest possible free energy when interacting with the pocket's residues. Therefore, that pose with the different substitutions from one compound to the other created the unique binding interactions with the pocket amino acids.
30% of the simulation time in the selected trajectory. (C) A timeline representation of CK2-Y49 interactions presented in (A). The top panel presents the total number of specific interactions the protein made with the ligand over the course of the trajectory. The bottom panel presents the residues interacting with the ligand in each trajectory frame. The dark orange color indicates more than one specific interaction is made between some residues and the ligand. Figure 11 shows the amino acid residues of the protein binding pocket that interacted with compound 19. The fused ring system of 19 made similar interactions with the pocket residues as Y49, where Val116 and Lys86 interacted through the H-bond with the ring system ( Figure 11B). As previously seen in the molecular surface display (Figure 4A), the extended aliphatic chain occupied a distant pocket and created new interaction points between the two OH groups at the end of the chain and Asp120, where they occurred > 85% during the simulation ( Figure 11B). This interaction was not present with Y49 ( Figure  3B,D). It might be safe to assume that the enhanced binding affinity and stability of the complex were due to this new interaction with Asp120, which can also be inferred from the RMSD plot ( Figure 8A). The top panel presents the total number of specific interactions the protein made with the ligand over the course of the trajectory. The bottom panel presents residues interacting with the ligand in each trajectory frame. The dark orange color indicates that more than one specific interaction is made between some residues and the ligand.
Compound 31 also created new interacting points with amino acids in the main binding pocket: the chiral OH of the fused ring system made a strong H-bonding with His160 and Asn161 through bridging water molecules with a value of ~1.18 and 0.6, respectively. An enhanced interaction with Arg47 was observed as well (~0.98) ( Figure 12B). The tetrahydropyran ring at the end of the aliphatic chain of 31 also extended along the protein surface but did not occupy the distant pocket, as did compound 19 ( Figure 5A). The reason for that might be the fact that the chain in compound 31 is 6-carbon shorter than that of 19, so the group could not reach the second pocket. Another explanation could be the large size of the substituted tetrahydropyran hindered the group from occupying the pocket. Additionally, there is a high probability that the fluctuation observed in the RMSD of 31 towards the end of the simulation time ( Figure 9A) might be a result of the inability of this group to bind to the second pocket. The L-RMSF (ligand-root-mean-square fluctuation) represents how the atoms of the ligand interact with the protein and the changes in the positions of the ligand atoms. As seen in the L-RMSF plot for compound 31 (Figure 13), the positions of atoms 29-45 were dramatically changed because of the free rotation around the aliphatic chain, which in turn decreased the interaction between this part of the molecule with the protein and was reflected by > 3 Å fluctuation in the RMSF plot. The time-depended representation of the CK2α1-31 interactions showed that residues Arg47, His160, and Asp175 were the ones making specific interactions with the ligand, as indicated by the darker color in the plot ( Figure 12C).   It was also noticed that the reference ligand, Y49, as well as the compounds 19 and 31, all interacted with Asn118 through H-bonding; however, the contact points are different. While the residue interacted with Y49 at the purine carbonyl oxygen ( Figure 10B), it acted as an H-bond donor to the OH at the end of the aliphatic chain of 19 as well as to the OH at the fused ring system of 19 and 31 ( Figures 11B and 12B). The reason for the different contact points of Asn118 with the compounds was probably attributed to the 3D conformation of the compounds inside the binding pocket, which is affected by the nature of the substitutions on the nucleus. Each compound assumed a pose that had the lowest possible free energy when interacting with the pocket's residues. Therefore, that pose with the different substitutions from one compound to the other created the unique binding interactions with the pocket amino acids.

Conclusions
CK2 was related to many human illnesses, not only cancer, but also multiple sclerosis, cardiac hypertrophy, neurodegenerative and inflammatory disorders, cystic fibrosis, and virus infections [68]. It is noteworthy that the CK2 role is best recognized and investigated in cancer, where CK2 is almost positively upregulated, resulting in tumor progression because of its role in regulating nearly all the essential processes for developing apoptosis suppression [5,69]. It was reported that cancer cells rely on CK2 high levels compared to normal cells, which supports that the CK2 inhibitors can have a crucial contribution to cancer therapy development [5]. Recently, several compounds have been discovered and optimized via rational drug design approaches. Various structure-based drug design (SBDD) tools have been utilized for CK2 drug discovery for predicting possible compound and target interactions and their affinity. Various classes of natural metabolites, such as anthraquinones, benzoimidazoles, coumarins, pyrazolotriazines, and flavonoids, are recognized as CK2 inhibitors [70,71]. Fungal phenalenones are a fascinating class of fungal metabolites with diverse bioactivities that could be lead metabolites for drug discovery. With the aid of computational methods, i.e., ADMET, docking, and MD simulation, compounds 19 and 31 were identified as promising drug-like phenalenone derivatives that have better binding interactions and protein stability in a simulated aqueous

Conclusions
CK2 was related to many human illnesses, not only cancer, but also multiple sclerosis, cardiac hypertrophy, neurodegenerative and inflammatory disorders, cystic fibrosis, and virus infections [68]. It is noteworthy that the CK2 role is best recognized and investigated in cancer, where CK2 is almost positively upregulated, resulting in tumor progression because of its role in regulating nearly all the essential processes for developing apoptosis suppression [5,69]. It was reported that cancer cells rely on CK2 high levels compared to normal cells, which supports that the CK2 inhibitors can have a crucial contribution to cancer therapy development [5]. Recently, several compounds have been discovered and optimized via rational drug design approaches. Various structure-based drug design (SBDD) tools have been utilized for CK2 drug discovery for predicting possible compound and target interactions and their affinity. Various classes of natural metabolites, such as anthraquinones, benzoimidazoles, coumarins, pyrazolotriazines, and flavonoids, are recognized as CK2 inhibitors [70,71]. Fungal phenalenones are a fascinating class of fungal metabolites with diverse bioactivities that could be lead metabolites for drug discovery. With the aid of computational methods, i.e., ADMET, docking, and MD simulation, compounds 19 and 31 were identified as promising drug-like phenalenone derivatives that have better binding interactions and protein stability in a simulated aqueous physiological environment. The current work highlights the usefulness of these metabolites as lead for anticancer discovery. One of the important issues that require attentiveness is that several mechanistic studies are directed to the in silico methods because they provide information that cannot be obtained by other models and are less time-consuming. However, in vivo and in vitro investigations are warranted to strengthen the findings of in silico studies and provide opportunities for observing other mechanisms of the anticancer potential of these metabolites.

Target Prediction
The webserver SuperPred is a knowledge-based method that uses machine learning models for ATC code and target prediction of investigated compounds [57]. The machine learning model uses logistic regression and Morgan fingerprints of length 2048. The drugs approved by the WHO are classified by a drug classification system that connects the drugs' chemical properties and therapeutic properties and indications, where each classification is given a code called an Anatomical Therapeutic Chemical (ATC) code. Therefore, if a drug has more than one therapeutic indication, it is given an ATC code for each indication. The WHO has 6300 approved drugs that are linked to over 600,000 targets. Based on the hypothesis that compounds that have similar physiochemical properties exhibit similar biological effects, the webserver translates a user-defined compound into a structural fingerprint and compares this fingerprint to that of the WHO-approved drugs. When similarity is found, the webserver predicts the ATC code, the possible therapeutic target(s), and the putative therapeutic indication(s) for that compound. In other words, if an investigated compound is structurally similar to an approved drug, the compound is predicted to have biological activity on all possible targets of that drug. After targets are predicted, a probability score and a model accuracy score are reported. The probability represents the chance that the investigated compound will bind to a specific predicted target. The model accuracy reflects the performance accuracy of the used machine-learning model when predicting that specific target for the compound since the model performance differs between targets [57,72]. The targets and ATC codes for a library of investigated compounds were predicted using the SuperPred tool. The compounds that did not have the common target(s) of interest were excluded from further analysis. The ones sharing the common target(s) were advanced for the docking and further studies.

Ligand Preparation
Phenalenone derivatives were processed and prepared for docking using Schrodinger's LigPrep tool [40]. The 2D structures were converted to 3D and energy-minimized using OPLS3 force-field. After adding hydrogens, all possible ionization states and tautomeric forms were created at pH of 7.0 ± 0.2 by Epik; desalt option was also chosen. H-bonds were optimized by predicting the pKa of ionizable groups suing PROPKA [73].

Protein Preparation
CK2 crystal structure (PDB: 7BU4) was prepared using the Protein Preparation Wizard, added hydrogens to residues, changed covalent bonds to metal ions to zero-order, and created disulfide bonds. Water molecules > 5 Å from protein residues were deleted. Using Epik, the protonation state of residues was generated, and the formal charge on metal ions was adjusted. After removing the extra protein subunits of multi-subunit proteins and additional ligands, processing of the protein was refined by predicting the pKa of ionizable residues using PROPKA [73], and water molecules > 3 Å (not involved in water bridge) were removed. Finally, restrained minimization of the protein was applied using the OPLS4 force field.

Grid Generation and Docking
A grid box was generated around the co-crystallized ligand Y49 in the protein crystal structure (PDB: 7BU4) binding site using Glide's Receptor-Grid-Generation tool [62]. Inside this box is where the docking of the phenalenone compounds was performed. The nonpolar atoms were set for the VdW radii scaling factor by 1.0, and the partial charge cutoff was 0.25. Docking was then performed by the Schrodinger suite "Ligand Docking" tool [62,74]. The selected docking protocol was standard precision (SP), and the ligand sampling method was flexible. All other settings were default. Re-docking of the co-crystallized ligand (PDB: Y49) was performed to evaluate the docking method and docking of the investigated phenalenones followed.

ADMET Properties Prediction
The processed compounds were subjected to ADMET prediction using the QikPropmodule of the Schrodinger suite [63]. The descriptors: molecular weight (mol_MW), drug-likeness (#Stars), dipole moment (dipole), total solvent accessible surface area (SASA), number of hydrogen bond donors and acceptors (donorHB and acceptHB), predicted octanol-water partitioning (QPlogPo/w), predicted aqueous solubility (QPlogS), estimated binding to human serum albumin (QPlogKhsa), number of the possible metabolites (# metab), predicted blood-brain partitioning (QPlogBB), percentage of human oral absorption, predicted IC 50 for inhibiting HERG-K + channels (QPogHERG), central nervous system activity (CNS), and number of reactive functional groups present (#rtvFG), were predicted for these derivatives. The predicted values are compared to the recommended range derived from values determined/observed for 95% of known drugs.

MD Simulation
MD simulations were performed using Desmond software in the Schrodinger suite [65,66]. The protein-ligand complexes of interest were retrieved from the docking results where the force field was OPLS4. The complexes were tuned through the "System-Builder" tool to generate the solvated system for simulation. The solvent model was set as TIP3P, the selected box shape was orthorhombic, and the box dimensions were 10 Å. Na ions were added to neutralize the system. The simulation parameters were set up in the Molecular Dynamic tool, where the protein-ligand complexes were evaluated at pH 7.0 ± 0.2 over the 100 ns simulation time. The ensemble class was set as NPT in order to maintain the temperature and pressure constant during the run at 300 K and 1.01325 bar, respectively. After running the MD simulation, the generated results were analyzed.