In Silico Exploration of Microtubule Agent Griseofulvin and Its Derivatives Interactions with Different Human β-Tubulin Isotypes

Tubulin isotypes are known to regulate microtubule stability and dynamics, as well as to play a role in the development of resistance to microtubule-targeted cancer drugs. Griseofulvin is known to disrupt cell microtubule dynamics and cause cell death in cancer cells through binding to tubulin protein at the taxol site. However, the detailed binding mode involved molecular interactions, and binding affinities with different human β-tubulin isotypes are not well understood. Here, the binding affinities of human β-tubulin isotypes with griseofulvin and its derivatives were investigated using molecular docking, molecular dynamics simulation, and binding energy calculations. Multiple sequence analysis shows that the amino acid sequences are different in the griseofulvin binding pocket of βI isotypes. However, no differences were observed at the griseofulvin binding pocket of other β-tubulin isotypes. Our molecular docking results show the favorable interaction and significant affinity of griseofulvin and its derivatives toward human β-tubulin isotypes. Further, molecular dynamics simulation results show the structural stability of most β-tubulin isotypes upon binding to the G1 derivative. Taxol is an effective drug in breast cancer, but resistance to it is known. Modern anticancer treatments use a combination of multiple drugs to alleviate the problem of cancer cells resistance to chemotherapy. Our study provides a significant understanding of the involved molecular interactions of griseofulvin and its derivatives with β-tubulin isotypes, which may help to design potent griseofulvin analogues for specific tubulin isotypes in multidrug-resistance cancer cells in future.


Introduction
Griseofulvin is the secondary metabolite of several fungal species. Since its first isolation from Penicillium griseofulvum in 1939 [1], it has been used to treat dermatophyte infections [2]. The biosynthetic griseofulvin gene cluster including 13 gsf genes responsible for producing griseofulvin was first identified in Penicillium aetiopicum through shotgun sequencing and genome mining over a decade ago [3]. In addition to Penicillium, other genera of ascomycetes, including Aspergillus alliaceus, Xylaria flabelliformis, and Memoniella echinata, have shown the ability to produce griseofulvin as they harbor the gsf gene cluster in their genomes [4][5][6]. However, our recent study showed that only 7 out of 13 gsf genes might be essential for griseofulvin biosynthesis [7].
Mechanistically, griseofulvin stabilizes microtubule dynamics and mediates microtubule abnormalities such as misaligned chromosomes and multipolar spindles. These abnormalities suppress cell division dramatically and lead to a cell-death response [8]. Additionally, increased p53 accumulation in the nucleus of these cells revealed that they were

Virtual Screening and Molecular Docking Studies
In this study, PyRx was used for virtual screening of ligand candidates to determine the binding modes and affinities toward TUBB5. The PyRx results contain binding energies and RMSDs values for each binding mode of the 464 analogues are shown in Supplementary File S4. To further investigate the ability of binding griseofulvin and its derivatives to tubulin, the binding mode and interactions of griseofulvin and the most active derivatives with the human β-tubulin were examined by molecular docking studies using AutoDock. The 3D structure and docking mode of griseofulvin at the human β-tubulin active site are shown in Figure 1.

Virtual Screening and Molecular Docking Studies
In this study, PyRx was used for virtual screening of ligand candidates to determine the binding modes and affinities toward TUBB5. The PyRx results contain binding energies and RMSDs values for each binding mode of the 464 analogues are shown in Supplementary File S4. To further investigate the ability of binding griseofulvin and its derivatives to tubulin, the binding mode and interactions of griseofulvin and the most active derivatives with the human β-tubulin were examined by molecular docking studies using AutoDock. The 3D structure and docking mode of griseofulvin at the human β-tubulin active site are shown in Figure 1. The five derivatives with the lowest binding energy and the highest number of interactions in terms of hydrogen bonds and hydrophobic interactions were selected for further analysis. Table 1 shows two-dimensional (2D) structures of the best-docked ligand molecules, and computed binding energies towards the tubulin protein. The number of rotatable bonds for each molecule was defined by AutoDock. It is expected that molecules with a greater number of rotatable bonds will require a greater number of energy evaluations to converge to an energy minimum, due to a greater number of degrees of freedom and conformational states [30]. Moreover, it is known that increasing the number of rotatable bonds makes it more difficult for Autodock to find the correct docking conformations [31]. Based on our computational predictions, griseofulvin derivatives G1-G5, with the number of rotatable bonds ranging from 3 to 11, show lower binding energies ranging from −8.8 to −9.1 kcal mol −1 compared to griseofulvin with a binding energy of −6.6 kcal mol −1 with only 3 rotatable bonds. The structure analysis of selected derivatives suggests that functional groups containing oxygen, such as methoxy in G1 and aromatic ether group in derivative G2, provide a better interaction and higher affinity toward receptors. Moreover, structures containing benzene, such as dichlorobenzene in G3, fluoromethyl benzene in G4, and acetate ester group in derivative G5, could enhance the interaction between ligands and proteins.
Even though binding energy indicates how strong the interaction between ligand and β-tubulin protein could be, Moriguchi octanol-water partition coefficient (MlogP) values were calculated and considered for the drug-diffusion ability into the cells (Table 1). MlogP is a useful structure-based factor developed by Moriguchi et al. [32] to estimate and compare the distribution of drugs within the cells, organs, and the body [33]. As shown in Table 1, griseofulvin shows a lower partition coefficient (MlogP = 0.71) than other compounds, indicating a stronger ability of chemical to distribute through the body The five derivatives with the lowest binding energy and the highest number of interactions in terms of hydrogen bonds and hydrophobic interactions were selected for further analysis. Table 1 shows two-dimensional (2D) structures of the best-docked ligand molecules, and computed binding energies towards the tubulin protein. The number of rotatable bonds for each molecule was defined by AutoDock. It is expected that molecules with a greater number of rotatable bonds will require a greater number of energy evaluations to converge to an energy minimum, due to a greater number of degrees of freedom and conformational states [30]. Moreover, it is known that increasing the number of rotatable bonds makes it more difficult for Autodock to find the correct docking conformations [31]. Based on our computational predictions, griseofulvin derivatives G1-G5, with the number of rotatable bonds ranging from 3 to 11, show lower binding energies ranging from −8.8 to −9.1 kcal mol −1 compared to griseofulvin with a binding energy of −6.6 kcal mol −1 with only 3 rotatable bonds. The structure analysis of selected derivatives suggests that functional groups containing oxygen, such as methoxy in G1 and aromatic ether group in derivative G2, provide a better interaction and higher affinity toward receptors. Moreover, structures containing benzene, such as dichlorobenzene in G3, fluoromethyl benzene in G4, and acetate ester group in derivative G5, could enhance the interaction between ligands and proteins. and into the cells. However, the MlogP values of griseofulvin derivatives G1-G5 are in good agreement, ranging from 1.96 to 2.98 (Table 1). and into the cells. However, the MlogP values of griseofulvin derivatives G1-G5 are in good agreement, ranging from 1.96 to 2.98 (Table 1).   As a result, the five derivatives described above, as well as griseofulvin and taxol were then docked into the griseofulvin binding site of nine distinct β-tubulin isotypes. Taxol, an antimicrotubule agent with remarkable antimitotic activity, was used as the control for comparing the binding energies and interactions with tubulin isotypes in this study. Interestingly, griseofulvin-binding sites overlap with the taxol binding site at the beta-tubulin H6-H7 loop [27]. Docking studies revealed that griseofulvin and its derivatives have different binding conformations and energies, depending on the residue composition variations in and around the binding pocket of different β-tubulin isotypes (Table −8.8 1.96 Even though binding energy indicates how strong the interaction between ligand and β-tubulin protein could be, Moriguchi octanol-water partition coefficient (MlogP) values were calculated and considered for the drug-diffusion ability into the cells (Table 1). MlogP is a useful structure-based factor developed by Moriguchi et al. [32] to estimate and compare the distribution of drugs within the cells, organs, and the body [33]. As shown in Table 1, griseofulvin shows a lower partition coefficient (MlogP = 0.71) than other compounds, indicating a stronger ability of chemical to distribute through the body and into the cells. However, the MlogP values of griseofulvin derivatives G1-G5 are in good agreement, ranging from 1.96 to 2.98 ( Table 1).
As a result, the five derivatives described above, as well as griseofulvin and taxol were then docked into the griseofulvin binding site of nine distinct β-tubulin isotypes. Taxol, an antimicrotubule agent with remarkable antimitotic activity, was used as the control for comparing the binding energies and interactions with tubulin isotypes in this study. Interestingly, griseofulvin-binding sites overlap with the taxol binding site at the betatubulin H6-H7 loop [27]. Docking studies revealed that griseofulvin and its derivatives have different binding conformations and energies, depending on the residue composition variations in and around the binding pocket of different β-tubulin isotypes ( Table 2).  As shown in the heat map above, the binding energies of five griseofulvin derivatives with binding energies −7.29 to −10.3 kcal mol −1 are more negative than taxol and griseofulvin, with binding energies ranging from −6.29 to −8.11 and −6.8 to −7.38 kcal mol −1 , respectively. The studies exhibited that compounds G1 and G3 are good examples of lower binding energies, greater than −8.09 kcal mol −1 towards all β-tubulin isotypes, while for compound G2, interaction with βIIb, βIVa, βIVb, and βV are the strongest binding tubulin isotypes with binding energies −10.3, −9.88, −9.68, and −9.63 kcal mol −1 , respectively. Compound G5 had the best binding energies when docked to isoforms βV, βIVa, and βIII dominates. However, other derivatives docked well in the griseofulvin-binding site with a binding free energy of at least −7.36 kcal mol −1 .
The analysis of hydrogen bonds and hydrophobic interactions of griseofulvin derivatives docked with β-tubulin isotypes shows differences in the binding pocket of different β-tubulin, which are listed in Table 3. Worthy of note are griseofulvin and the hydrogen bonds formation of its derivatives with the residues Thr274; and Gln279 of β-tubulin isotypes IIa, IIb, III, IVa, IVb, and VIII. More hydrogen-bond interaction can be observed between griseofulvin and derivatives with residues His227, Arg276, Gln279, and Ser230. H-bond formation with the βI tubulin isotype shows differences at the residues level. Notably, the residues Lys19, Gln276, and Gln279 of the βI tubulin isotype play a predominant role in establishing hydrogen bonds. It is noticeable that the maximum hydrogen-bond distance is 3.5 Å, which corresponds to strong H bonds. Many more residues with significant interactions existed and were involved in hydrophobic interaction at the binding site of β-tubulin isotypes and ligands ( Table 3). The results of interaction analysis showed the dominance of hydrophobic interactions to hydrogen bonding in this study. Table 3. Binding energy as well as interactions detail of ligands with different human β-tubulin isotypes after molecular docking.
As shown in the heat map above, the binding energies of five griseofulvin derivatives with binding energies −7.29 to −10.3 kcal mol −1 are more negative than taxol and griseofulvin, with binding energies ranging from −6.29 to −8.11 and −6.8 to −7.38 kcal mol −1 , respectively. The studies exhibited that compounds G1 and G3 are good examples of lower binding energies, greater than −8.09 kcal mol −1 towards all β-tubulin isotypes, while for compound G2, interaction with βIIb, βIVa, βIVb, and βV are the strongest binding tubulin isotypes with binding energies −10.3, −9.88, −9.68, and −9.63 kcal mol −1 , respectively. Compound G5 had the best binding energies when docked to isoforms βV, βIVa, and βIII dominates. However, other derivatives docked well in the griseofulvin-binding site with a binding free energy of at least −7.36 kcal mol −1 .
The analysis of hydrogen bonds and hydrophobic interactions of griseofulvin derivatives docked with β-tubulin isotypes shows differences in the binding pocket of different β-tubulin, which are listed in Table 3. Worthy of note are griseofulvin and the hydrogen bonds formation of its derivatives with the residues Thr274; and Gln279 of β-tubulin isotypes IIa, IIb, III, IVa, IVb, and VIII. More hydrogen-bond interaction can be observed between griseofulvin and derivatives with residues His227, Arg276, Gln279, and Ser230. Hbond formation with the βI tubulin isotype shows differences at the residues level. Notably, the residues Lys19, Gln276, and Gln279 of the βI tubulin isotype play a predominant role in establishing hydrogen bonds. It is noticeable that the maximum hydrogen-bond distance is 3.5 Å, which corresponds to strong H bonds. Many more residues with significant interactions existed and were involved in hydrophobic interaction at the binding site of β-tubulin isotypes and ligands ( Table 3). The results of interaction analysis showed the dominance of hydrophobic interactions to hydrogen bonding in this study.   We then investigated the residue composition changes in the griseofulvin binding pocket of different β-tubulin isotypes. The multiple sequence analysis study demonstrates residue variations in human β-tubulin isotypes at different locations ( Figure 2). Griseofulvin shows differences in binding energy ranging from −6.8 to −7.38 kcal mol −1 with respect to the residue composition variations in and around the binding pocket of different β-tubulin isotypes. The griseofulvin binding pocket of βI has five residue changes, including Arg276-Gln, Lys362-Ser, Val23-Met, Ala231-Leu, and Asp26-Glu, compared to other β-tubulin isotypes. Here, Leu215, Lue217, and Cys211 interact hydrophobically with the methoxy group of griseofulvin. Afterwards, an analysis of the docking complex of βI-tubulin isotype-griseofulvin shows that methoxy and carbonyl groups of griseofulvin form hydrogen bonding interactions with residues Gln276 and Gln279, respectively. Thus, these results suggest that the residue composition variation in and around the griseofulvin binding pocket results in differences in binding energy, conformation, and hydrogen bonding interactions among the different β-tubulin isotype-griseofulvin complexes.
Phylogenetic reconstruction with PhyML based on aligned amino acid sequences of α and β tubulin proteins ( Figure 3) shows that βI (TUBB1) resulted from the earliest gene duplication within the β tubulin clade. This explains the sequence divergence of TUBB1 from the rest of the β tubulin proteins in Figure 2. The earliest lineage within the α tubulin clade is TUBAL3 (Figure 3). TUBB1 and TUBAL3 have the highest isoelectric point (pI = 4.93 and 5.98, respectively) among all the tubulin proteins in Figure 3, whereas the rapidly evolving TUBB8B has the smallest pI (=4.59) and features the smallest number of lysine residues. There is a tendency for more recently derived paralogues to have low pI values, consistent in both α and β clade (Figure 3), although the biological significance of this pI change is not obvious. All the tubulin proteins are expected to be negatively charged at normal cytoplasmic pH near 7 because their pI values are all smaller than 7. methoxy group of griseofulvin. Afterwards, an analysis of the docking complex of βI-tubulin isotype-griseofulvin shows that methoxy and carbonyl groups of griseofulvin form hydrogen bonding interactions with residues Gln276 and Gln279, respectively. Thus, these results suggest that the residue composition variation in and around the griseofulvin binding pocket results in differences in binding energy, conformation, and hydrogen bonding interactions among the different β-tubulin isotype-griseofulvin complexes.  clade is TUBAL3 (Figure 3). TUBB1 and TUBAL3 have the highest isoelectric point (pI = 4.93 and 5.98, respectively) among all the tubulin proteins in Figure 3, whereas the rapidly evolving TUBB8B has the smallest pI (=4.59) and features the smallest number of lysine residues. There is a tendency for more recently derived paralogues to have low pI values, consistent in both α and β clade (Figure 3), although the biological significance of this pI change is not obvious. All the tubulin proteins are expected to be negatively charged at normal cytoplasmic pH near 7 because their pI values are all smaller than 7. TUBA1A   TUBB4B   TUBB8B   TUBAL3   TUBB3   TUBB4A   TUBB1   TUBB   TUBA8   TUBA3E   TUBB2A   TUBA4A   TUBA1C   TUBA3C_TUBA3D   TUBA1B   TUBB2B   1  The study was continued for compound G1 (CID 25171849) with the best binding energies and best interactions, in terms of hydrogenic and hydrophobic bonds towards all the β-tubulin isotypes for a better understanding of the effect of the β-tubulin isotypes interactions with a griseofulvin derivative. The structural model of the binding mode and conformation of compound G1 at the binding site of nine distinct β-tubulin isotypes are shown in Table 4. Table 4. The 3D-and 2D-binding conformations and active residues involved in interactions between compound G1 and nine distinct β-tubulin isotypes.

Β-Tubulin Isotype
Binding Energy (kcal/mol) 3D Interactions 2D Interactions The study was continued for compound G1 (CID 25171849) with the best binding energies and best interactions, in terms of hydrogenic and hydrophobic bonds towards all the β-tubulin isotypes for a better understanding of the effect of the β-tubulin isotypes interactions with a griseofulvin derivative. The structural model of the binding mode and conformation of compound G1 at the binding site of nine distinct β-tubulin isotypes are shown in Table 4. Table 4. The 3D-and 2D-binding conformations and active residues involved in interactions between compound G1 and nine distinct β-tubulin isotypes.

Molecular Dynamic (MD) Simulations
Molecular docking results suggested that compound G1 (CID 25171849) serves as the best interacting compound with all the β-tubulin isotypes and, therefore, it was consid-

Molecular Dynamic (MD) Simulations
Molecular docking results suggested that compound G1 (CID 25171849) serves as the best interacting compound with all the β-tubulin isotypes and, therefore, it was considered for the MD analysis. The cytotoxicity of this compound against the breast cancer cell

Molecular Dynamic (MD) Simulations
Molecular docking results suggested that compound G1 (CID 25171849) serves as the best interacting compound with all the β-tubulin isotypes and, therefore, it was considered for the MD analysis. The cytotoxicity of this compound against the breast cancer cell line MDA-MB-231 has been shown by MTT colorimetric assay in a previous study [34].

Molecular Dynamic (MD) Simulations
Molecular docking results suggested that compound G1 (CID 25171849) serves as the best interacting compound with all the β-tubulin isotypes and, therefore, it was considered for the MD analysis. The cytotoxicity of this compound against the breast cancer cell line MDA-MB-231 has been shown by MTT colorimetric assay in a previous study [34]. The compound G1 showed the highest number of hydrogenic and hydrophobic interactions with the lowest free binding energies ranging from −8.65 to −9.32 kcal mol −1 (Table 5), corresponding to the highest affinity towards β-tubulin isotypes. To investigate the protein structure stability and conformational changes upon binding ligands, we performed molecular dynamic (MD) simulations over the β-tubulin-G1 docked complexes for 200 ns. The trajectory files from MD simulations were analyzed for root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), radius of gyration (Rg), solvent-accessible surface area (SASA), and molecular mechanics Poisson-Boltzmann surface area (MM-PBSA). The root-mean-square deviations (RMSD) were calculated to examine the stability of the molecular dynamics simulation. The RMSD analysis complexes suggest that most β-tubulin-G1 complexes reached their equilibrium conformation after a time period of 100 ns and then, retained their stability during the simulation ( Figure 4A). However, βIIb still seems to be fluctuating even after 180 ns from the RMSD graph in Figure 4A. It is notable that unfavorable interaction with His227 of G1 at tubulin βI resolved in the first few ns of the MD simulation; and then, a new conformation is stabilized, as also suggested by the high RMSD at this isotype. The average RMSD values were found to be the lowest for the G1 complexes with the βII (0.29 Å), followed by βIVa (0.32 Å) and βIII (0.37 Å) isotypes. No significant variations were observed in the average values of RMSFs, demonstrating the consistent complex exposure for a favorable conformational cavity of the β-tubulin proteins ( Table 5). The average Rg values show that the packing of β-tubulin upon ligand binding is not considerably different from β-tubulin protein alone during the last 20 ns of MD simulations (Table 5). From the SASA plot in Figure 4D, we can see that the SASA values of almost all the β-tubulin-G1 complexes are greater overall than the β-tubulin without ligand, indicating the increases of solvent-accessible surface area and stability of the protein. The highest SASA value was observed in βIII (200.1 Å 2 ), followed by βV (193.7Å 2 ) and βIV (191.6 Å 2 ), revealing the protein volume expansion and low variation throughout the simulation time (Table 5).
proteins ( Table 5). The average Rg values show that the packing of β-tubulin upon ligand binding is not considerably different from β-tubulin protein alone during the last 20 ns of MD simulations (Table 5). From the SASA plot in Figure 4D, we can see that the SASA values of almost all the β-tubulin-G1 complexes are greater overall than the β-tubulin without ligand, indicating the increases of solvent-accessible surface area and stability of the protein. The highest SASA value was observed in βIII (200.1 Å 2 ), followed by βV (193.7Å 2 ) and βIV (191.6 Å 2 ), revealing the protein volume expansion and low variation throughout the simulation time (Table 5).  Table 5. The average of RMSD, RMSF, Rg, and SASA for nine distinct β-tubulin isotypes during the last 20 ns of MD simulations. MM/PBSA was originally defined by Kollman et al. [35], which can be used as a reliable and practical approach to predicting binding affinities. Using the last 30 ns steady RMSD trajectories, 1000 snapshots were taken at regular intervals to study the binding free energy upon binding, which is one of the most common techniques for calculating MM/PBSA was originally defined by Kollman et al. [35], which can be used as a reliable and practical approach to predicting binding affinities. Using the last 30 ns steady RMSD trajectories, 1000 snapshots were taken at regular intervals to study the binding free energy upon binding, which is one of the most common techniques for calculating interaction energies of biomolecular complexes. Table 6 shows the calculated average values of the various binding energies, such as the van der Waals, electrostatic, polar solvation, and SASA energies. Our hit compound G1 showed the average free binding energy of −34.6 kcal mol −1 , which is favorable for binding to β-tubulin isotypes. These findings confirm the docking analysis, indicating that compound G1 has a strong binding affinity towards β-tubulin isotypes, as seen in Table 6.

Tubulin Isotype Expression in Cancer
Clinical observations showed altered expression of tubulin isotypes in a range of cancers. These β-tubulins are encoded by multiple genes that are expressed tissue-specifically; e.g., βI is expressed in bone marrow and lymphoid tissue, βII in brain, βIII in neuronal and testicular cells, βIV in neuronal and glial cells, βVI in the erythroid cells and platelets, and βVIII is observed in testicular cells (Table 4). High expression of several β-tubulin isotypes, including I, II, III, IVa, and V tubulins, has been linked to aggressive clinical behavior, chemotherapy drug resistance, and poor patient outcome in many cancers, including lung, breast, ovarian, and gastric cancers. β-tubulin class III has been extensively investigated because of its role in resistance to antimitotic drugs. For example, it has been demonstrated that paclitaxel-resistant lung cancer cells and ovarian tumors have higher levels of β-tubulin classes III and IVa isotypes [36]. Tubulin isotypes βIVa and βIVb accounted for more than half of the β-tubulin in breast and colon cancer cells [37]. Detailed information on β-tubulin isotypes in cancer cells is provided in Table 7.

Ligand and Receptor Preparation
The 3D SDF file formats for griseofulvin (compound CID 441140), taxol (compound CID 36314), as well as 464 griseofulvin analogues, were retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). To have a comprehensive study, all the available data for griseofulvin derivatives from PubChem database were included in this study. The biological activity of some griseofulvin derivatives such as antifungal and anticancer activities was reported in previous studies; however, there is limited information on the biological activity of most derivatives recorded in PubChem. The CIDs and more details about 464 griseofulvin derivatives are listed in Supplementary File S1. The pdbqt ligand file was then preprocessed for docking by adding partial charges, atom types, and assigning the torsional flexibility of atoms to the ligand; in addition, their optimization was generated using Open Babel version 2.4.1.
Since there is no crystal structure available for some β-tubulin isotypes, homology modeling was performed for protein sequences of human Tubulin beta-5 chain (TUBB5) using MODELLERv18 [45]. This sequence was retrieved from UniProt (ID P07437) and then aligned against Protein Data Bank deposited sequences using BLAST. A protein of known tubulin structure (PDB ID: 6I2I) with the highest sequence identity to the target protein was chosen as a template for homology modeling. The template and the target file were used as input to MODELLER to align the sequence and to build the model structure. The best model was chosen based on the discrete optimized potential energy (DOPE) score. The energy minimization of the model structure was carried out using Chimera 1.15 software, and the quality of the protein structure was evaluated using PROCHECK by Ramachandran plot [46]. Evaluation of the model using Ramachandran plot showed that 87.1% of residues lie in favored regions, followed by 11.6% in allowed regions, indicating a good quality model (Supplementary File S2). In addition, the Ramachandran plot of the human Tubulin beta-5 chain (TUBB5) protein retrieved from AlphaFold was evaluated and the results showed that the number of residues present in the favoured region is 94.6%, and the allowed region is 4.9% (Supplementary File S3). This finding shows an agreement between the results of the structure generated with homology modeling and AlphaFold. Therefore, the nine different human β-tubulin isotypes proteins were retrieved from AlphaFold database for further analysis in this study [47]. The AlphaFold IDs of these β-tubulin sequences are as follows: βI (Q9H4B7), βIIa (Q13885), βIIb (Q9BVA1), βIII (Q13509), βIVa (P04350), βIVb (P68371), βV (P07437), βVI (Q9BUF5), and βVIII (Q3ZCM7).

Sequence Analysis of β-Tubulin Isotypes
To identify the various residues at the griseofulvin binding pocket, the sequences of nine different human β-tubulin isotypes were aligned with MAFFT with the slow, but accurate, G-INS-i option [48].

Virtual Screening and In Silico Molecular Docking
PyRx 0.8 [49] from MGL Tools (https://ccsb.scripps.edu/mgltools/) was used for the virtual screening of griseofulvin analogues library molecules against β-tubulin, using default settings. In this study, TUBB5 (AlphaFold ID P07437) was used as a receptor for screening using AutoDock vina in PyRx. A grid box of size 23 Å × 25 Å × 26 Å with a spacing of 0.375 Å and co-ordinates x = 200.01, y = 268.74, and z = 334.49 was prepared at the β-tubulin. The five lowest-scoring analogues of griseofulvin, griseofulvin, and taxol were then docked into the β-tubulin isotypes binding sites using the Lamarckian genetic algorithm in AutoDock 4 software [50]. The geometric center of the grid box was placed at the macromolecule with coordinates x = −19.276, y = 9.865, and z = −1.569. A grid box of size 70 Å × 64 Å × 70 Å with a spacing of 0.375 Å was prepared at the β-tubulin to cover the putative griseofulvin binding site. The ligand with the lowest binding-free energies and the best modes of interactions with different β-tubulin binding sites was considered for further analysis. The interactions between the receptors and the ligands were visualized using Discovery Studio Visualizer v.20 (https://discover.3ds.com/discovery-studio-visualizerdownload).

Molecular Dynamics Simulation
To validate the docking results, molecular dynamics (MDs) simulation was performed for the best derivative-docked complexes with nine different human β-tubulin isotypes; i.e., βI, βIIa, βIIb, βIII, βIVa, βIVb, αβV, βVI, and αβVIII, using the GROMACS-2022 software package [51] to carry out 200 ns simulations with an AMBER99SB force field. The Ante Chamber PYthon Parser interface (ACPYPE) was used to create the force field parameters for the ligand [52]. The TIP3P water model was selected for solvating complexes, and then sodium and chloride ions were added to neutralize the system. The conditions of periodic boundary (PBC) were used. Energy minimization was carried out at 1000 kJ/mol/nm. The system was equilibrated in the NVT and NPT ensembles for 1 ns. Sufficient amount of Na+ and Cl counter-ions were added to achieve the overall system charge neutrality and ionic strength of 0.15 M. To maintain the temperature and pressure at 310 K and 1.0 bar, the Nose-Hoover thermostat and the Berendsen barostat were used, respectively [53]. The trajectories were generated every 2 fs and saved every 2 ps. Post-MD analyses were performed; these included root-mean-square deviation (RMSD), root-mean-square fluctuations (RMSF), the radius of gyration (Rg), solvent-accessible surface area (SASA), and molecular mechanics Poisson-Boltzmann surface area (MM-PBSA).

Conclusions
The FDA-approved drug griseofulvin has shown the potential for suppressing cancer cell division and inducing cell death through interaction with the mitotic spindle microtubule [8], but its mechanism and interaction with tubulin molecules have remained poorly understood. In this study, sequence analysis, molecular docking, molecular dynamics simulations, and binding-free energy calculations were used to investigate the binding affinity of griseofulvin and its derivatives with nine human β-tubulin isotypes I, IIa, IIb, III, IVa, IVb, V, VI, and VIII. The human βI-tubulin isotype was found to have different residue compositions at the griseofulvin binding pocket, while no significant differences were found in other β-tubulin isotypes.
The results of molecular docking demonstrate that griseofulvin has a remarkable affinity toward β-tubulin isotypes, i.e., at the taxol binding site, as reported previously [27]. Griseofulvin exhibited different binding modes and binding energies for different β-tubulin isotypes; this could be due to changes in residue composition in and around the binding pocket of different β-tubulin isotypes. The residues Gln276 and Gln279 in βI-tubulin, and residue Thr274 in other β-tubulin isotypes are involved in the hydrogen bonding interactions with griseofulvin. Interestingly, our docking analysis showed that the binding energies of griseofulvin toward β-tubulin isotypes are close to the binding energies estimated for taxol. Further, the binding free-energy calculation revealed that derivatives G1-G5 have more negative binding energies ranging from −7.29 to −10.3 kcal/mol −1 and they could bind favorably to human β-tubulin isotypes even better than taxol with binding energies ranging from −6.29 to −8.11 kcal/mol −1 .
The structure analysis of selected derivatives revealed functional groups such as methoxy, and aromatic ether groups in derivatives G1 and G2 provided better interactions and higher affinities toward receptors. Furthermore, structures containing benzene in derivatives G3 and G4, and the acetate ester group in derivative G5 enhanced the interactions between ligands and receptors.
Molecular dynamics simulations were performed on β-tubulin isotype-G1 docked complexes, to further investigate the structural stability and conformational changes of β-tubulin isotypes upon binding to ligands. Furthermore, β-tubulin isotypes IVb, followed by IIb and IIa, showed lower (more negative) binding energies with −30.8, −25.9, and −25.3 kcal mol −1 corresponding to more favorable interactions. From MM-PBSA, it is evident that compound G1 has the most favorable binding to the βI-tubulin isotype with a free-binding energy of −34.6 kcal mol −1 . The residue changes at the binding site could be one of the reasons behind the more negative binding-free energy of the G1 derivative toward the βI-tubulin isotype. Similarly, a previous study investigated the differential binding affinity of an anti-mitotic agent colchicine analogue to different αβ-tubulin isotypes. They found the changed residues at the binding site of an αβIII isotype were responsible for affecting the binding affinity of a colchicine analogue [54]. Another interesting study explored the interaction of the microtubule-depolymerizing agent indanocine with different human αβ tubulin isotypes. Their findings suggested differences in amino acid sequences at the indanocine binding pockets of βI, βIIa, βIII, and βVI isotypes, which resulted in less binding-free energy toward indanocine [55].
One of the limitations of the present study is the using of AlphaFold β-tubulin isotypes protein structures to perform molecular docking analysis, as there were no available crystal structures for some β-tubulin isotypes. However, homology modeling was performed for one of these β-tubulin isotypes, and evaluation of the model using a Ramachandran plot showed a good agreement between the results of the structure generated with homology modeling and AlphaFold structure. The other limitation of the current study is that molecular dynamics simulations were performed only for the G1 derivative due to the time-consuming data analysis of molecular dynamics simulations and our limited computational resources, although other derivatives such as G2 showed good binding affinities and favorable interactions toward β-tubulin isotypes. Our current computational study provides a significant understanding of the involved molecular interactions of griseofulvin with human β-tubulin isotypes, which provides insight for designing potent griseofulvin analogues for specific tubulin isotypes. These superior analogues may improve the treatment of patients with advanced carcinoma that demonstrate tubulin isotype specificity, as well as in the development of innovative new cancer medicines in the future.