**1. 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–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

**Citation:** Ibrahim, S.R.M.; Bagalagel, A.A.; Diri, R.M.; Noor, A.O.; Bakhsh, H.T.; Muhammad, Y.A.; Mohamed, G.A.; Omar, A.M. Exploring the Activity of Fungal Phenalenone Derivatives as Potential CK2 Inhibitors Using Computational Methods. *J. Fungi* **2022**, *8*, 443. https://doi.org/10.3390/jof8050443

Academic Editors: Tao Feng and Frank Surup

Received: 18 March 2022 Accepted: 21 April 2022 Published: 24 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

composed of two catalytic alpha proteins, CK2α and/or CK2α', that differ only in their C-termini [9–11], and two regulatory CK2β proteins that are responsible for substrate specificity [11–13]. It has minimal substrate specificity; therefore, it is able to phosphorylate a large number of proteins involved in multiple signal transduction pathways [14]. This enzyme is ubiquitously expressed [7] and constitutively active in cells; hence, its activity is not relied on extracellular stimuli [5,15,16]. Many reports have found that the overexpression of CK2 in many cancer types leads to inhibiting apoptosis and promoting cancer cell survival [5,8]. It is clear that its downregulation revered cancer progression and enhanced apoptosis [17], suggesting that CK2 can be considered a valuable target for anticancer agents [8,18].

Natural metabolites and their derivatives have long been established as a wealthy source for the discovery of novel anticancer drugs [19]. It was reported that ~5% of the anticancer drugs approved by the FDA are unmodified natural metabolites, and ~50% of the approved drugs are either chemically-modified natural metabolites or synthesized relied on natural metabolites structures [20].

It is noteworthy that different metabolites belonging to various classes, such as flavonoids, coumarins, anthraquinones, pyrazolotriazines, polyhalogenated benzimidazoles, and benzotriazoles, have been known as CK2 inhibitors [21]. Additionally, CIGB-300 and CX-4945 are CK2 inhibitors that are already in human trials as anticancer drugs. CX-4945 (Silmitasertib) has been designated by the FDA for treating cholangiocarcinoma, and many clinical investigations are ongoing with it versus different types of cancers [22]. Furthermore, CIGB-300 is under investigation for cervical cancers [23].

Fungi are eukaryotic micro-organisms that inhabit nearly all kinds of environments in nature, where they play a fundamental role in maintaining eco-balance [24–28]. Fungi have drawn immense research attention since their cultivation, isolation, characterization, and purification demonstrated the existence of a vast number of metabolites with unique and diversified chemical classes and bioactivities [29–34]. Many reports revealed the characterization of diverse classes of fungal metabolites with CK2 inhibitory potential [35–37].

Phenalenones are a unique class of secondary metabolites, having a fused three-ring system [38]. They are produced by higher plants and fungi [39]. They are known as phytoalexins, which are biosynthesized by the plant in response to an external threat, such as mechanical injury or pathogenic infections [40]. Fungal phenalenones have immense structural diversity, such as hetero- and homo-dimerization and high degrees of nitrogenation and oxygenation, as well as being complexed with metals, incorporating with additional carbon frameworks, or isoprene unit by the formation of either a linear ether (e.g., **8** and **12**) or a trimethyl-hydrofuran (e.g., **3**–**7**, **9**–**11**, **13**, and **14**) moiety. It has been reported to possess various bioactivities, including antimicrobial, anti-plasmodial, anticancer, antidiabetic, antioxidant, and cytotoxic effects [40]. These metabolites exhibit cytotoxic capacities towards various cancer cell lines; however, there are limited or lacking studies that explore the mechanism of their anticancer properties. It is noteworthy that there is no available report on their CK2 inhibitory potential. The molecular docking-based virtual screening process, along with increasing the data about the structure of CK2 alone or in complex with various inhibitors, could become of particular relevance in the discovery of new lead compounds as CK2 inhibitors [41,42]. In continuation of our interest in exploring the bioactivities of fungal metabolites, the present work focuses on the in silico assessment of CK2 inhibitory capacity of 33 phenalenone derivatives reported from various fungal sources (Table 1, Figures 1 and 2).

**Figure 1.** Structures of phenalenone derivatives **1**–**20**. **Figure 1.** Structures of phenalenone derivatives **<sup>1</sup>**–**20**.

**Figure 2. Figure 2.** Structures of phenalenone derivatives Structures of phenalenone derivatives **2121**––**3333**.

.

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.


#### **2. Results and Discussion**

### *2.1. 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.

**Table 2.** Prediction of target probability and model accuracy for phenalenone derivatives against CK2 using SuperPred target prediction webserver.


\* The probability of the test compound binding to a specific target, as determined by the respective target machine learning model. \*\* The accuracy of the performance of the prediction model displaying the 10-fold cross-validation score of the respective logistic regression model, as the model performance varied between different targets.

#### *2.2. 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-7*H*-purin-2 yl)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–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-7*H*-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,6*E*,10*E*,14*E*)-2,6,10,14-tetramethylhexadeca-6,10,14 triene-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.

**Figure 3.** Re-docking of the co-crystallized ligand to validate the docking method. The figure showed the 2D view of the binding interactions of the reference ligand **Y49** complexed with CK2α1; (**A**) after minimization of the crystal structure 7BU4, and (**B**) after re-docking of ligand **Y49** into the CK2α1 crystal structure. (**C**) 3D structure of the re-docked **Y49** (pink) superimposed on the co-crystallized Y49 (green). (**D**) Molecular surface display with electrostatic potential color scheme for CK2α1 complexed with ligand **Y49** after re-docking. After docking validation, docking the 3D structures of the > 100 phenalenones that **Figure 3.** Re-docking of the co-crystallized ligand to validate the docking method. The figure showed the 2D view of the binding interactions of the reference ligand **Y49** complexed with CK2α1; (**A**) after minimization of the crystal structure 7BU4, and (**B**) after re-docking of ligand **Y49** into the CK2α1 crystal structure. (**C**) 3D structure of the re-docked **Y49** (pink) superimposed on the co-crystallized Y49 (green). (**D**) Molecular surface display with electrostatic potential color scheme for CK2α1 complexed with ligand **Y49** after re-docking.

> 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.


**Table 3.** In silico docking results of phenalenone derivatives with CK2α1 (PDB: 7BU4).

**Figure 4.** CK2α1 in complex with **19**. (**A**) Molecular surface display with electrostatic potential color scheme for CK2α1- **19** complex and the close-up view presented. (**B**) 3D presumed binding mode of **19** in the binding site of CK2 (PDB:7BU4). Compound **19** is displayed as green ball-and-sticks. The amino acids in the binding site are shown as grey sticks, and hydrogen bonds are represented in yellow dotted lines. (**C**) 2D depiction of the ligand-protein interactions. **Figure 4.** CK2α1 in complex with **19**. (**A**) Molecular surface display with electrostatic potential color scheme for CK2α1-**19** complex and the close-up view presented. (**B**) 3D presumed binding mode of **19** in the binding site of CK2 (PDB:7BU4). Compound **19** is displayed as green ball-and-sticks. The amino acids in the binding site are shown as grey sticks, and hydrogen bonds are represented in yellow dotted lines. (**C**) 2D depiction of the ligand-protein interactions.

**Figure 5.** CK2α1 in complex with **31**. (**A**) Molecular surface display with electrostatic potential color scheme for CK2α1- **31** complex and the close-up view presented. (**B**) 3D presumed binding mode of **31** in the binding site of CK2 (PDB:7BU4). Compound **31** is displayed as green ball-and-sticks. The amino acids in the binding site are shown as grey sticks, and hydrogen bonds are represented in yellow dotted lines. (**C**) 2D depiction of the ligand-protein interactions. **Figure 5.** CK2α1 in complex with **31**. (**A**) Molecular surface display with electrostatic potential color scheme for CK2α1-**31** complex and the close-up view presented. (**B**) 3D presumed binding mode of **31** in the binding site of CK2 (PDB:7BU4). Compound **31** is displayed as green ball-and-sticks. The amino acids in the binding site are shown as grey sticks, and hydrogen bonds are represented in yellow dotted lines. (**C**) 2D depiction of the ligand-protein interactions.

**Figure 6.** Molecular display with electrostatic potential color scheme for CK2α1 (PDB: 7B4U) showing **19** (orange) and **31** (cyan) superimposed on **Y49** (green) inside the binding pocket. **Figure 6.** Molecular display with electrostatic potential color scheme for CK2α1 (PDB: 7B4U) showing **19** (orange) and **31** (cyan) superimposed on **Y49** (green) inside the binding pocket.

#### *2.3. In Silico ADMET Properties of Selected Ligands 2.3. In Silico ADMET Properties of Selected Ligands*

ble 4.

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 Ta-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.


**Table 4.** Predicted in silico ADMET properties of the phenalenone derivatives.


**Table 4.** *Cont.*


and vice versa; Dipole: computed dipole moment; SASA: Total solvent accessible surface area; DonorHB: estimated number H+ to be donated in HB; AcceptHB: estimated number H+ to be accepted in HB; QLogPo/w: predicted octanol/water partition coefficient; QPlogS: Predicted aqueous solubility; QPlogKhsa: Prediction of binding to human serum albumin; #Metab: number of possible metabolic reactions; QPlogBB: Predicted brain/blood partition coefficient; % Human Oral Absorption: Predicted human oral absorption on 0 to 100% scale; QPlogHERG: Predicted IC50 value for blockage of HERG K+ channels; CNS: Predicted central nervous system activity; #RtvFG: Number of reactive functional groups.

#### *2.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.

*2.4. Molecular Dynamic (MD) Simulation* 

energy.

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 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

**Figure 7.** (**A**) The RMSD plot was obtained for the reference **Y49** complexed with CK2α1 protein (PDB-ID:7BU4). The simulation time (100 ns) reaffirmed the stability of the complex with no significant changes in the protein structure. (**B**) Stability of the secondary structure of CK2α1 protein (PDB ID: 7BU4) was evaluated by monitoring its SSE distribution (top plot), SSE composition (middle plot), and SEE assignment (bottom plot) over the 100 ns of MD simulation when complexed with **Y49. Figure 7.** (**A**) The RMSD plot was obtained for the reference **Y49** complexed with CK2α1 protein (PDB-ID:7BU4). The simulation time (100 ns) reaffirmed the stability of the complex with no significant changes in the protein structure. (**B**) Stability of the secondary structure of CK2α1 protein (PDB ID: 7BU4) was evaluated by monitoring its SSE distribution (top plot), SSE composition (middle plot), and SEE assignment (bottom plot) over the 100 ns of MD simulation when complexed with **Y49**.

**Figure 8.** (**A**) The RMSD plot was obtained for compound **19** complexed with CK2α1 protein (PDB-ID:7BU4). The simulation time (100 ns) reaffirmed the stability of the complex with no significant changes in the protein structure. (**B**) Stability of the secondary structure of CK2α1 protein (PDB ID: 7BU4) was evaluated by monitoring its SSE distribution (top plot), SSE composition (middle plot), and SEE assignment (bottom plot) over the 100 ns of MD simulation when complexed with **19**. **Figure 8.** (**A**) The RMSD plot was obtained for compound **19** complexed with CK2α1 protein (PDB-ID:7BU4). The simulation time (100 ns) reaffirmed the stability of the complex with no significant changes in the protein structure. (**B**) Stability of the secondary structure of CK2α1 protein (PDB ID: 7BU4) was evaluated by monitoring its SSE distribution (top plot), SSE composition (middle plot), and SEE assignment (bottom plot) over the 100 ns of MD simulation when complexed with **19**.

**Figure 9.** (**A**) Obtained RMSD plot for compound **31** complexed with CK2α1 protein (PDB-ID:7BU4). The simulation time (100 ns) reaffirmed the stability of the complex with no significant changes in the protein structure. (**B**) Stability of the secondary structure of CK2α1 protein (PDB ID: 7BU4) was evaluated by monitoring its SSE distribution (top plot), SSE composition (middle plot), and SEE assignment (bottom plot) over the 100 ns of MD simulation when complexed with **31. Figure 9.** (**A**) Obtained RMSD plot for compound **31** complexed with CK2α1 protein (PDB-ID:7BU4). The simulation time (100 ns) reaffirmed the stability of the complex with no significant changes in the protein structure. (**B**) Stability of the secondary structure of CK2α1 protein (PDB ID: 7BU4) was evaluated by monitoring its SSE distribution (top plot), SSE composition (middle plot), and SEE assignment (bottom plot) over the 100 ns of MD simulation when complexed with **31**.

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). *J. Fungi* **2022**, *8*, x FOR PEER REVIEW 19 of 30 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

> 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). 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.

**Figure 10.** *Cont*.

**Figure 10.** (**A**) Stacked bar graph for CK2 interactions with reference **Y49** throughout the simulation. (**B**) Schematic diagram shows the detailed 2D atomic interactions of **Y49** with CK2 that occurred > **Figure 10.** (**A**) Stacked bar graph for CK2 interactions with reference **Y49** throughout the simulation. (**B**) Schematic diagram shows the detailed 2D atomic interactions of **Y49** with CK2 that occurred > 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.

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. 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).

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

*J. Fungi* **2022**, *8*, x FOR PEER REVIEW 21 of 30

specific interaction is made between some residues and the ligand.

**Figure 11.** *Cont*.

**Figure 11.** (**A**) CK2 interactions with compound **19** throughout the simulation. (**B**) Schematic diagram shows the detailed 2D atomic interactions of **19** with CK2 that occurred > 30% of the simulation time in the selected trajectory. (**C**) A timeline representation of CK2-**19** 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 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 bind-**Figure 11.** (**A**) CK2 interactions with compound **19** throughout the simulation. (**B**) Schematic diagram shows the detailed 2D atomic interactions of **19** with CK2 that occurred > 30% of the simulation time in the selected trajectory. (**C**) A timeline representation of CK2-**19** 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 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. *J. Fungi* **2022**, *8*, x FOR PEER REVIEW 23 of 30 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).

**Figure 12.** *Cont*.

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 indi-

cated by the darker color in the plot (Figure 12C).

**Figure 12.** (**A**) CK2 interactions with compound **31** throughout the simulation. (**B**) Schematic diagram shows the detailed 2D atomic interactions of **31** with CK2 that occurred > 30% of the simulation time in the selected trajectory. (**C**) A timeline representation of CK2-**31** 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 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 12.** (**A**) CK2 interactions with compound **31** throughout the simulation. (**B**) Schematic diagram shows the detailed 2D atomic interactions of **31** with CK2 that occurred > 30% of the simulation time in the selected trajectory. (**C**) A timeline representation of CK2-**31** 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 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 13.** Ligand RMSF shows the fluctuations of **31** broken down by atom as represented by the compound's 2D structure. It provides ideas on how ligand atoms interact with the protein and their entropic role in the binding event. The 'Fit Ligand on Protein' line presents the ligand fluctuations, with respect to the protein. The CK2α1-**31** complex was first aligned on the protein backbone, and then the ligand RMSF was measured on the atoms of the ligand. It was also noticed that the reference ligand, **Y49,** as well as the compounds **19** and **Figure 13.** Ligand RMSF shows the fluctuations of **31** broken down by atom as represented by the compound's 2D structure. It provides ideas on how ligand atoms interact with the protein and their entropic role in the binding event. The 'Fit Ligand on Protein' line presents the ligand fluctuations, with respect to the protein. The CK2α1-**31** complex was first aligned on the protein backbone, and then the ligand RMSF was measured on the atoms of the ligand.

#### **31,** all interacted with Asn118 through H-bonding; however, the contact points are differ-**3. Conclusions**

ent. 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. **3. 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 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.

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

#### **4. Materials and Methods**

#### *4.1. 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.

#### *4.2. Preparation of PDB Structures*

#### 4.2.1. 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].

### 4.2.2. 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.

#### *4.3. 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 cocrystallized ligand (PDB: **Y49**) was performed to evaluate the docking method and docking of the investigated phenalenones followed.

#### *4.4. 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<sup>50</sup> for inhibiting HERG-K<sup>+</sup> 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.

#### *4.5. 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.

**Author Contributions:** Conceptualization, G.A.M., A.M.O., and S.R.M.I.; methodology, G.A.M., A.M.O., S.R.M.I., and Y.A.M.; validation, A.A.B., R.M.D., A.O.N., and H.T.B.; resources, A.A.B., R.M.D., A.O.N., and H.T.B.; data curation, G.A.M., A.M.O., S.R.M.I., and Y.A.M.; software, A.M.O. and Y.A.M.; writing—original draft preparation, A.M.O., S.R.M.I., and Y.A.M.; writing—review and editing, A.A.B., R.M.D., A.O.N., and H.T.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Journal of Fungi* Editorial Office E-mail: jof@mdpi.com www.mdpi.com/journal/jof

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-0365-6143-1