A Molecular Docking Study Reveals That Short Peptides Induce Conformational Changes in the Structure of Human Tubulin Isotypes αβI, αβII, αβIII and αβIV

Microtubules are cylindrical protein polymers assembled in the cytoplasm of all eukaryotic cells by polymerization of aβ tubulin dimers, which are involved in cell division, migration, signaling, and intracellular traffic. These functions make them essential in the proliferation of cancerous cells and metastases. Tubulin has been the molecular target of many anticancer drugs because of its crucial role in the cell proliferation process. By developing drug resistance, tumor cells severely limit the successful outcomes of cancer chemotherapy. Hence, overcoming drug resistance motivates the design of new anticancer therapeutics. Here, we retrieve short peptides obtained from the data repository of antimicrobial peptides (DRAMP) and report on the computational screening of their predicted tertiary structures for the ability to inhibit tubulin polymerization using multiple combinatorial docking programs, namely PATCHDOCK, FIREDOCK, and ClusPro. The interaction visualizations show that all the best peptides from the docking analysis bind to the interface residues of the tubulin isoforms αβl, αβll, αβlll, and αβlV, respectively. The docking studies were further confirmed by a molecular dynamics simulation, in which the computed root-mean-square deviation (RMSD), and root-mean-square fluctuation (RMSF), verified the stable nature of the peptide–tubulin complexes. Physiochemical toxicity and allergenicity studies were also performed. This present study suggests that these identified anticancer peptide molecules might destabilize the tubulin polymerization process and hence can be suitable candidates for novel drug development. It is concluded that wet-lab experiments are needed to validate these findings.


Introduction
Microtubules are cylindrical protein polymers assembled in the cytoplasm of all eukaryotic cells by polymerizing αß-tubulin heterodimers, which are involved in multiple biological functions such as cell division, migration, signaling, and intracellular traffic. These functions, especially the formation of mitotic spindles, make them important participants in the initiation and proliferation of cancerous cells and their subsequent metastases [1,2]. They exhibit an irregular temporal pattern of assembly and disassembly, which has been termed dynamic instability [3]. In addition, altered tubulin polymerization is detrimental to the formation of the tumor vascular system during the process of angiogenesis, which is a hallmark of cancer. Hence, unsurprisingly, the microtubular cytoskeleton and its building blocks, tubulin dimers, have become one of the key targets in cancer chemotherapy [4]. It

Peptide Screening and Preparation
We explored the DRAMP (data repository of antimicrobial peptides) database in this work [15]. DRAMP is an open-access and manually curated database containing diverse annotations of AMP, including information on their sequences, structures, biological activities, physicochemical, patent, clinical, and references of the peptides. DRAMP currently contains 22,407 entries, 6032 of which are general AMPs (having both natural and synthetic AMPs), 16,110 patented AMPs, 77 AMPs in drug development (preclinical or clinical stage), and 188 stapled antimicrobial peptides belonging to specific AMPs. To expand the scope of AMPs' design, DRAMP also holds 5909 candidate AMPs screened by some platforms whose antibacterial activities have not been assayed yet. In this study, we screened 826 peptides derived from plant sources. Anticancer peptides were extracted from the database. The peptides that did not have accurate information regarding their activity or function in the database were discarded. The peptides with more than 50 amino acids in length were also removed as being unlikely to hold promise for drug development. In addition, anticancer peptides with nonentity amino acids in their sequences were eliminated. The PEP-FOLD 2 webserver was used further to predict the information structure of 13 sequence peptides [16]. This web server visualizes the peptide structure using a hidden Markov model suboptimal sampling algorithm [17]. It uses a coarse-grained energy force field without including conformational entropy.

Homology Modeling and Protein Preparation
The sequences of human α and β tubulin isotypes were downloaded from the universal protein resource (UniProt) [18] (Q71U36: α-tubulin, P02773, Q13885, Q13509 and, P04350: βl-βlV-tubulin proteins). The 3D structure of bovine tubulin (Q2HJ86 and Q6B856: αβ-tubulin chain (Bos taurus)) was downloaded from the Protein Data Bank (PDB ID: 1SA0) as the template structure. The complex includes two tubulin αβ heterodimers, while the native ligand, namely, colchicine, is bound to β subunits at the interface with the α subunit. The homology model of the human tubulin was built with SWISS-MODEL [19]. The native ligand in the template structure was deleted before homology modeling was performed. The stereochemical quality of the template (referred to as tubulin 1SA0 hereafter) and the various tubulin isotypes were evaluated using PROCHECK [18]. Subsequently, the Verify-3D [20] and ERRAT [21] algorithms were employed to check the reliability of the generated models. The Ramachandran plots for all the isotypes are shown in Figure S1. All the structures resulting from the modeling αβl, αβll, αβlll, and αβlV were further subjected to protein preprocessing, optimization, and minimization using the protein preparation in the Schrodinger suite [22].

Molecular Docking
After cleaning the proteins using the protein preparation tool in the Schrodinger software suite, the output was used for the protein-peptide molecular docking. PATCH-DOCK [23] and ClusPro [24] web servers were used to analyze their interactions. FIRE-DOCK (Fast Interaction Refinement) webserver was used for further refinement [25]. The top 10 model solutions were determined and analyzed. The PATCHDOCK web server with default parameters was used to predict the best conformations. Since a rigid approach was utilized to obtain the docking solutions and because during protein−protein interactions, both side chains and backbones might change their conformation. The top 10 solutions were subjected to the FIREDOCK webserver to refine the interaction of protein-protein complexes resulting from molecular docking [25]. FIREDOCK performs side chain optimization and rigid body minimization to provide more extensive refinement. Subsequently, scoring and ranking identified the near-native advanced solutions. The final selection of the best-docked complexes was based on the global energy of the bound, predicted complexes after the refinement.

Molecular Dynamics
The molecular dynamics (MD) simulations were evaluated using the Desmond simulation package embedded in the Schroedinger suite [26]. Before the MD simulations were performed, the molecular system, including the protein, water molecules, and ions, was built. The water molecules were described using the TIP3P (transferable intermolecular potential with 3 points) model in an orthorhombic cubic box. The boundary condition box volume was computed according to the complex type and counter ions, including Na + and Cl − , to neutralize the system. The NPT ensemble with temperature and pressure of 300 K and 1 bar, respectively, was used in all the analyses for all the selected docked peptides in the αβl-αβlV tubulin receptor. The force field of OPLS_2005 was applied. Simulation progress was verified stepwise every 50 ps. The NPT assembly was launched following the simulation process, which covers the production of 100 ns. The frames were assembled and examined using the simulation interaction diagram that helped determine the fluctuations.

Post Molecular Dynamics MM-GBSA
The Molecular Mechanics-Generalized Born and Surface Area Continuum Solvation (MM/GBSA) was calculated for each MD trajectory using the thermal MM/GBSA script [27]. This was run via the Python command line. The average binding energy was calculated for 20 snapshots from the overall trajectory of 100 ns.

Physicochemical, Allergenicity, and Toxicity Prediction
Physicochemical, allergenicity, and toxicity assessments of the best peptides from the docked analysis were evaluated. The physiochemical properties of the best peptides were evaluated using ProtParam tools [28], which allow predicting the relevant properties, such as molecular weight, net charge at pH 7, peptide properties, stability, and charge. The AllerTop server was used to evaluate the nonallergenic nature of the peptides [29].

Results and Discussion
We retrieved plant peptides from the DRAMP database. The 13 peptides shown in Table 1 and Figure 1 were extracted from the set of plant-based peptides. Figure 2 displays the respective structures of the peptides. Our 3D predicted structures may not represent the native conformation in the plant, but they do represent at least one version of what might be created in the wet lab experiments. Which are linear cyclotides and they are aimed to be anticancer lead compounds.      The Ramachandran plots obtained by PROCHECK found 85.1% of residues in the most favored regions, 11.4% in the additional allowed regions, and 2.3% in generously allowed regions for αβI (see Table S1 and Figure S1). For αβII, 83.4% of residues in most favored regions, 13.3% in additional allowed regions, and 1.1% in generously allowed regions were found. In the αβIII, 83.7% of residues in the most favored regions, 12.9% in additional allowed regions, and 2.4% in generously allowed regions were observed. On the other hand, in αβIV tubulin, 84.1% of residues are found in most favored regions, 12.7% in additional allowed regions, and 2.3% in generously allowed regions. The results from the overall quality factor obtained by the ERRAT agreed with the acceptable range (>50 for a high-quality model) [30]. The overall quality factor obtained by the ERRAT tool for the isotypes αβI-αβIV was found to be 91.32%, 80.86%, 82.68%, and 90.7%, respectively. Subsequently, VERIFY3D was used to confirm that at least 80%, 97.56%, 96.52%, and 97.79% of the amino acids of all the analyzed isotypes achieved an average 3D/1D score higher than 0.2. In particular, most residues are found in the most favored regions; therefore, the built model was deemed reliable.
The identification of the binding site is significant for elucidating the actual binding mechanisms as well as the interaction between a drug and a protein. To characterize this, we used molecular docking and MD simulations, which are standard, generally accepted theoretical prediction approaches for examining peptide binding sites in the various protein receptors. ClusPro was used to reexamine the peptides with the bestdocked results from the FIREDOCK. The ClusPro score of docked peptides in αβIII and αβIV tubulin corroborates the result of the FIREDOCK. However, the results of αβI and αβII-tubulin obtained from CluPro and FIREDOCK do not corroborate it. Further, the peptides with the lowest global energy from the FIREDOCK ( Table 2) were utilized for further analysis. The MD simulation analysis was carried out on the complexes in order to answer many relevant inquiries, such as the stability and accuracy of the binding mode and the ligand stability over the simulation period. Most of the interactions predicted in the docking analysis of the lead peptides were stable within the binding site. They interacted appropriately with different target regions by forming hydrogen bonds, as well as hydrophobic and electrostatic connections directly with protein side chains.
bic interactions at αTyr222, βLeu42, and αAla18. The DRAM00789 peptide and αβ-tubulin receptor complex displayed ten hydrogen bonds at αGln279, βMet321, βMet323, αThr219 βAsp355, βGln245, αGly360, βSer322, and βASP355 and one hydrophobic bond a βArg320 (see Figure 2b, Table 3). Peptide binding to the receptor may cause the targe molecule to inhibit its interactions. All these results confirmed that the four peptides could effectively repress tubulin interactions by occupying the interface of the αβI-tubulin di mer.    The molecular dynamics simulations of the active peptides bound to the αβ tubulin heterodimer were carried out to confirm the structural rigidity and to validate the docking outcomes for the complexes. RMSD values of C-alpha atoms (Cα) have been studied to understand structural inflexibility. The results demonstrated that the DRAMP00776, DRAMP0782, DRAMP00783, and DRAM00789 complexes showed initial RMSD increases due to instability. The increase in RMSD was higher for the DRAM00789 peptide complex than for the other complexes. This result corroborates the molecular docking result. The trend in the RMSD was maintained between 40 and 65 ns, after which instability occurred. Meanwhile, DRAMP00776, DRAMP0782, and DRAMP00783 maintained stability from 70 ns to the end of the simulation period. Although the DRAMP00789 complex had a higher RMSD trend during the initial phase, the RMSD profile decreased with an increase in the simulation time (from 40 ns to 100 ns, see Figure 4). The average RMSD profile of the complexes formed with DRAMP00776, DRAMP00782, and DRAMP00783 was between 2.5 Å and 2.8 Å, respectively. The overall RMSD estimate for DRAMP00779 was 3.5 Å, and fluctuations were observed throughout the simulation time, indicating overall structural instability because higher RMSD profiles for Cα atoms are indicators of low stability ( Figure 4A). DRAMP00776 formed hydrogen bonds with the residues, αLys19, αAsn226, αArg276, αGln15, αGly71, αSer75, αGly79, αPro80, αGly223, and αGlu22, respectively. Compared to the interacting residues before MD simulations, DRAMP00783 peptides shift from the initial binding position. DRAMP00783 interacts with αGln15, βAsp355, βGln245, αAsp74, and βAsp355 via hydrogen bonding, which is accompanied by three hydrophobic interactions via βMet323, βVal353, and αTyr222 (Table 2). It was also observed that the higher the hydrophobic interaction, the better the binding energy. RMSF plots provide information on the flexible regions of the MD simulated structures. Flexible regions display a higher RMSF value, while constrained regions display a low RMSF value. There were no significant fluctuations in amino acid residues in the β-monomer after binding the peptides. The calculated RMSF for DRAMP00776, DRAMP00782, and DRAMP00783 was analyzed, and the peptides exhibited lower flexibility, below the 4 Å range, as shown in Figure 4B. due to instability. The increase in RMSD was higher for the DRAM00789 peptide complex than for the other complexes. This result corroborates the molecular docking result. The trend in the RMSD was maintained between 40 and 65 ns, after which instability occurred. Meanwhile, DRAMP00776, DRAMP0782, and DRAMP00783 maintained stability from 70 ns to the end of the simulation period. Although the DRAMP00789 complex had a higher RMSD trend during the initial phase, the RMSD profile decreased with an increase in the simulation time (from 40 ns to 100 ns, see Figure 4). The average RMSD profile of the complexes formed with DRAMP00776, DRAMP00782, and DRAMP00783 was between 2.5 Å and 2.8 Å, respectively. The overall RMSD estimate for DRAMP00779 was 3.5 Å, and fluctuations were observed throughout the simulation time, indicating overall structural instability because higher RMSD profiles for Cα atoms are indicators of low stability (Figure 4A). DRAMP00776 formed hydrogen bonds with the residues, αLys19, αAsn226, αArg276, αGln15, αGly71, αSer75, αGly79, αPro80, αGly223, and αGlu22, respectively. Compared to the interacting residues before MD simulations, DRAMP00783 peptides shift from the initial binding position. DRAMP00783 interacts with αGln15, βAsp355, βGln245, αAsp74, and βAsp355 via hydrogen bonding, which is accompanied by three hydrophobic interactions via βMet323, βVal353, and αTyr222 (Table 2). It was also observed that the higher the hydrophobic interaction, the better the binding energy. RMSF plots provide information on the flexible regions of the MD simulated structures. Flexible regions display a higher RMSF value, while constrained regions display a low RMSF value. There were no significant fluctuations in amino acid residues in the β-monomer after binding the peptides. The calculated RMSF for DRAMP00776, DRAMP00782, and DRAMP00783 was analyzed, and the peptides exhibited lower flexibility, below the 4 Å range, as shown in Figure 4B. (A)

Peptide-αβII-Tubulin Receptor Complex
The binding region surrounding the peptides consists of residues within a 4 Å range from the peptide atoms (Figure 5a). Subsequently, two peptides were found to have potent binding energy, namely, DRAMP00779 and DRAMP00788, respectively, after docking the extracted peptides to the αβII-tubulin dimer. For the DRAMP00779 complex, the residues αArg221, αGlu279, αGlu77, αThr82, αArg229, αThr225, αThr82, αTyr83, αAla19, αGly365 were detailed in the α-tubulin and βGln245, βGlu45 in the β-tubulin subunit (Figure 5b). While residues αArg221, αThr225, and αTyr224 were detailed in the α-tubulin βLeu246, βAsp41, and βAsp41 in the β-tubulin subunit for the DRAMP00788 complex ( Figure 5b, Table 4). The trajectories obtained from a set of MD simulations correspond to conformational changes of the complex of DRAMP00779 and DRAMP00788, and the corresponding data for individual biological systems were analyzed. The difference in the RMSD of the liganded and the unliganded tubulin heterodimer is not so dramatic ( Figure  6A). The structural motions of the residues in each tubulin subunit achieve additional sta-

Peptide-αβII-Tubulin Receptor Complex
The binding region surrounding the peptides consists of residues within a 4 Å range from the peptide atoms (Figure 5a). Subsequently, two peptides were found to have potent binding energy, namely, DRAMP00779 and DRAMP00788, respectively, after docking the extracted peptides to the αβII-tubulin dimer. For the DRAMP00779 complex, the residues αArg221, αGlu279, αGlu77, αThr82, αArg229, αThr225, αThr82, αTyr83, αAla19, αGly365 were detailed in the α-tubulin and βGln245, βGlu45 in the β-tubulin subunit (Figure 5b). While residues αArg221, αThr225, and αTyr224 were detailed in the α-tubulin βLeu246, βAsp41, and βAsp41 in the β-tubulin subunit for the DRAMP00788 complex (Figure 5b, Table 4). The trajectories obtained from a set of MD simulations correspond to conformational changes of the complex of DRAMP00779 and DRAMP00788, and the corresponding data for individual biological systems were analyzed. The difference in the RMSD of the liganded and the unliganded tubulin heterodimer is not so dramatic ( Figure 6A). The structural motions of the residues in each tubulin subunit achieve additional stability upon binding the peptides. The primary residues involved in the interaction of the DRAMP00779 complex after 100 ns of simulation time include αGly81, αArg221, βMet323, αGly81, αThr80, βAsp355, αArg221, αGly365, βMet321 (conventional and carbon-hydrogen bonds). Hydrophobic and electrostatic interactions also stabilized the complex. While DRAMP00788 interacted with residues, αAsp76, αGlu77, αArg221, αGly365, αAsp76, αGlu77, αArg221, αThr80, αGly81, βMet323, βLys324, βMet323, βAsp355, βLeu42, βMet321, and βLys324. Furthermore, the peptide formed distinct interactions such as electrostatic, hydrophobic, and hydrogen bonds with the intermediate domain residues of the αβ tubulin dimer (206-381). The domain opens with helices H6 and H7, a long loop, and helix H8 at the longitudinal interface sandwiched between the monomers. Many areas of the outer surface of tubulin are negatively charged and can attract hydrogen ions [31]. These regions play a role in tubulin-tubulin interactions and tubulin interactions with motor proteins such as kinesin. The MD results demonstrate that DRAMP00779 interacts electrostatically more strongly than DRAMP00788. To explore the atomic fluctuations at the binding sites in detail, we determined the RMSF value of each amino acid residue at the αβ tubulin binding site ( Figure 6B). The residues fluctuated at the N-terminal domain (amino acids 1 to~205). The atomic fluctuation profile shows a low fluctuation of the DRAMP00779 and DRAMP00788-tubulin complex in the α-tubulin compared to the β-tubulin monomer.  [31]. These regions play a role in tubulin-tubulin interactions and tubulin interactions with motor proteins such as kinesin. The MD results demonstrate that DRAMP00779 interacts electrostatically more strongly than DRAMP00788. To explore the atomic fluctuations at the binding sites in detail, we determined the RMSF value of each amino acid residue at the αβ tubulin binding site ( Figure 6B). The residues fluctuated at the N-terminal domain (amino acids 1 to ~205). The atomic fluctuation profile shows a low fluctuation of the DRAMP00779 and DRAMP00788-tubulin complex in the α-tubulin compared to the βtubulin monomer.

Peptides-αβIII-Tubulin Dimer Complex
The DRAMP00776 and DRAMP00781 peptides bind to the αβlll-tubulin dimer with binding energies of −43.38 and −32.08 kcal/mol, respectively, which represents a strong binding affinity of the peptides to the αβlll-tubulin dimer ( Table 2 and Figure 7a). The docking result gives an idea about the principal binding sites of DRAMP00776 and DRAMP00781 on the α-tubulin monomer surface, to a lesser extent on the β-tubulin monomer surface. Figure 7b and Table 5 indicate that DRAMP00776 expressly interacts with αGlu22, αThr82, αArg229, αAsn18, αGlu22, αThr82, and αTyr224 on α-tubulin, and it also interacts with residues, βGly244, βLeu246, βPro243, on β-tubulin. While DRAMP00781 binds to residues αThr225, αArg229, αGlu22, αGlu22, αThr225, αGln11, αPro364, and αTyr224 on the α-tubulin surface, and residues βGly244, and βLeu246 on the β-tubulin surface. The αβ-tubulin interface exhibits hydrogen bonding with DRAMP00776 and DRAMP00781, within a 4 Å range. The docking results suggest that the hydrogen bonds are vital for binding between DRAMP00776 and DRAMP00781 and the αβlll-tubulin dimer. In contrast, hydrophobic and electrostatic interactions also provide a minimal contribution to stabilizing the complexes. To investigate the conformational fluctuations, the aggregate RMSD deviations of tubulin atomic coordinates were considered during the interaction of DRAMP00776 and DRAMP00781 with the αβlll-tubulin dimer ( Figure 8A). A specific ligand should interact with the H6-H7 loop to efficiently inhibit the switch of nucleotide in tubulin (T216-Y224) [32,33]. Meanwhile, the H6-H7 loop favors hydrophobic interactions; thus, DRAMP00776 formed a hydrophobic interaction with αTyr224, while DRAMP00781 formed hydrophobic interactions with the βPro243 residue in the H6-H7 loop. Both peptides interact with residues at the H7-H8 loop (βGln245 and βGly244) via hydrogen bonding. As mentioned above, both peptides have a similar binding landscape and exhibit similar RMSD values of~3.2 Å. The calculation of RMSF is an exceptional tool to determine local protein mobility. However, the rate of hydrogen bond formation may be directly implicated in the flexibility of peptides in the binding site, as shown in the RMSF graph ( Figure 8B). contribution to stabilizing the complexes. To investigate the conformational fluctuations, the aggregate RMSD deviations of tubulin atomic coordinates were considered during the interaction of DRAMP00776 and DRAMP00781 with the αβlll-tubulin dimer ( Figure 8A). A specific ligand should interact with the H6-H7 loop to efficiently inhibit the switch of nucleotide in tubulin (T216-Y224) [32,33]. Meanwhile, the H6-H7 loop favors hydrophobic interactions; thus, DRAMP00776 formed a hydrophobic interaction with αTyr224, while DRAMP00781 formed hydrophobic interactions with the βPro243 residue in the H6-H7 loop. Both peptides interact with residues at the H7-H8 loop (βGln245 and βGly244) via hydrogen bonding. As mentioned above, both peptides have a similar binding landscape and exhibit similar RMSD values of ~3.2 Å. The calculation of RMSF is an exceptional tool to determine local protein mobility. However, the rate of hydrogen bond formation may be directly implicated in the flexibility of peptides in the binding site, as shown in the RMSF graph ( Figure 8B).

Peptide-αβIV-Tubulin Receptor Complex
The DRAMP00776 and DRAMP00784 bind to the αβlV-tubulin receptor ( Table 2 and Figure 9a) with binding energies of −42.27 and −48.97 kcal/mol, respectively, which shows a strong binding affinity of the peptides with the αβlV-tubulin receptor compared to their counterpart. The visualization of the docking result details the binding sites of DRAMP00784 and DRAMP00776 on αβIV-tubulin. DRAMP00776 and DRAMP00784 display a similar extent for the α and β-tubulin monomers. Figure 9b and Table 6 show that DRAMP00776 specifically interacts with αGlu113, and αLys96, at α-tubulin, and it also interacts with residues, βGlu158, βAsp161, βPro160, and βGlu125, in β-tubulin. By comparison, DRAMP00784 forms attractive charge interactions with residues αGln31, αThr82, αAsn228, αArg229, βGly244, and αGlu22, respectively. This is accompanied by four hydro-gen bonds with αThr82, αThr225, αGln15, βGly244, and βGln245 (Figure 9b, Table 6). The αβ-tubulin interface exhibits hydrogen bonding with DRAMP00784 and DRAMP00776 within a 4 Å range. The stability of the DRAMP00784 and DRAMP00776 peptides with the proteins was observed for the complete length of the simulation. The RMSD analysis shows that the protein complexes made with DRAMP00784 and DRAMP00776 were stabilized after 70 ns of the simulation time, and this was maintained for DRAMP00776 until the completion of the simulation ( Figure 10A). Meanwhile, the RMSD of DRAMP00784 displayed a deviation between 80 ns and 90 ns and continued in a stable conformation until the completion of the simulation. The two peptides maintained contact with the proteins by hydrogen bond interactions with residues αThr109, βPhe159, βPro160, αThr94, αGly95, βAsp128, βCys127 αGly95, βGlu125, and βAsp128 for DRAMP00776 and αGlu77, αAsn228, αGln15, αThr73, αAsn18, and βAsp41 for DRAMP00784. The RMSF values for each amino acid residue in the protein backbone are depicted in Figure 10B. The peaks represent the fluctuation of every amino acid residue throughout the simulation. This means that higher RMSF values represent greater flexibility of residues, while lower RMSF values reflect lower flexibility of residues and better system stability. The RMSFs of the α and β-domains are depicted separately. A slight amount of fluctuation in residues present at the active site indicates minimal conformational change, implying that the reported lead compound was consistently bound to the target protein.   depicted in Figure 10B. The peaks represent the fluctuation of every amino acid residue throughout the simulation. This means that higher RMSF values represent greater flexibility of residues, while lower RMSF values reflect lower flexibility of residues and better system stability. The RMSFs of the α and β-domains are depicted separately. A slight amount of fluctuation in residues present at the active site indicates minimal conformational change, implying that the reported lead compound was consistently bound to the target protein.

Post-MM-GBSA Analysis
The post-MM-GBSA analysis of the free binding energy calculation was carried out by processing and analyzing the frame generated during molecular dynamics with a 10step sampling size using the thermal MM-GBSA script in the Schrödinger suite. These are generally based on simulations of the molecular dynamics of the receptor-ligand complex. Consequently, they are intermediate in precision and computational effort between empirical scoring and strict alchemical disturbance methods [34]. These estimates have been applied to various systems with varying degrees of success. The free energies of binding have been enhanced after the post-MM-GBSA for all the analyzed compounds. The estimated values were negative for all the active docked complexes, of which DRAMP00776 showed the most negative MM/GBSA value (−109.02 kcal/mol), followed by DRAMP00783 (−87.04 kcal/mol) and DRAMP00789 (−64.01 kcal/mol) for the αβI-tubulin complex. DRAMP00789 peptides have low binding energy, which may have led to the instability of the peptides during the MD simulations, as represented in the RMSD and RMSF results. The results show that compound DRAMP00776 has a substantial binding energy based on the docked result and a significant post-MM-GBSA compared to the other compounds. These results indicate that the inhibition activity of DRAMP00776 against αβI-tubulin is likely to be more substantial compared to the other compounds. DRAMP00779 has a binding energy of −93.66 kcal/mol, whereas DRAMP00788 has a binding energy of −80.72 kcal/mol, which directly reveals the complex's stability for αβIItubulin. Notably, the high binding free energy was apparently due to strong electrostatic interactions, which demonstrate significant contributions of residues βLys324, αGlu77, αAsp76, βAsp355, and βAsp39, which were lacking in the binding of DRAMP00788. DRAMP00776 (−107.20 kcal/mol) exhibits superior predicted binding energy to the αβIIItubulin dimer compared to the DRAMP00781 peptide (−106.12 kcal/mol). Moreover, both peptides show good stability results in this study, which corroborate with their ClusPro and the global energy value obtained in the docking analysis. Nevertheless, this provides initial evidence that these peptides represent potential candidates for the development of novel therapeutic applications to cancer cell inhibition, which renders them suitable for further examination. The binding energy calculation of the peptides-αβIV-tubulin complex shows that the average binding free energies of the DRAMP00776 and DRAMP00784 docked complexes are −102.97 and −104.90 kcal/mol, respectively. The resulting binding energy difference could be due to hydrophobic interactions between DRAMP00784 and the αTyr224 residue (H6-H7 loop), which are lacking in the DRAMP00776 docked complex. The MM-GBSA data suggest that the intermolecular electrostatic, hydrophobic, and hydrogen bonding interactions are vital in the binding of the peptides to the respective tubulin receptor sites. The more pronounced stabilization observed in the evaluated peptides appears to correlate with their binding energy. Importantly, the activity or inactivity of the peptides that inhibit tubulin polymerization corresponds to their tubulin-binding ability as assessed by the molecular dynamics simulation results [11].

Physicochemical, Allergenicity, and Toxicity Prediction
The physicochemical properties of the peptide sequences of interest were generated using the Expasy web server tool. The following properties were examined: length, aliphatic index, instability, and molecular weight (Table 7). Peptides DRAMP00782 and DRAMP00789 were determined to be unstable, with an instability index of 50.50 and 50.46, respectively. These peptides will most likely also be found unstable in vitro because a value of the instability index above 40 is considered unstable. Thus, preventive procedures are to be taken to stabilize the unstable therapeutic peptides through suitable biochemical processes. In comparison, the remaining peptides are predicted to be stable. Further, the sizes of the peptides range from 378 to 434 atoms, with molecular weights between 2902 and 3199 Dalton. The theoretical isoelectric point (pI) denotes the respective pH of the peptides. The predicted aliphatic indices of the peptides were found to be 84.84, 46.90, 74.67, and 50.35 for peptides DRAMP00776, DRAMP00782, DRAMP00783, and DRAMP00789, respec-tively. The aliphatic index often indicates the relative volume of the aliphatic lateral chains (alanine, valine, isoleucine, and leucine). This is a positive factor in increasing the thermostability of globular proteins. In terms of residue charge, DRAMP00776, DRAMP00779, and DRAMP00783 tended to be charged negatively. AllerTOP was used for the in silico prediction of allergens based on the primary physicochemical properties of proteins. Meanwhile, the application uses the amino acid z-descriptors, ACC protein transformation, and k nearest neighbors clustering parameters. Most of the peptides are non-allergic and non-toxic, apart from the DRAMP00779 and DRAMP00783 peptides, which are identified as possible allergens. The ToxinPred results for all the active peptides showed that they were non-toxic compared to the mutated peptides. Interestingly, all the peptides belong to the cyclotide family, with most of them classified as plant defensins. Cyclotides are macrocyclic peptides with a knotted arrangement of three disulfide bonds formed from their six conserved cysteine residues. They contribute to their exceptional stability and natural functions as plant defense peptides. Cyclotides have many pharmaceutically relevant activities, especially their significance in drug design. Several synthetic cyclotides have also been made for applications in drug design. The first cyclotides were generated during the rise of natural products through their discovery as active compounds in studies that screened plant extracts for medicinal properties [35][36][37][38]. They were initially discovered because of their uterotonic activity by identifying Kalata B1 as an active agent from Oldenlandia affinis, thus used in African traditional medicine as a utertonic tea to quicken childbirth [39]. Circulin A and circulin B isolated from Chassalia parvifolia extract have been reported to act as anti-HIV agents [40]. Hence, further explorations of the best plant-based peptides analyzed in this work, as tubulin polymerization agents are hoped to be of great benefit. After the MD analysis, all the peptides were found to lack interactions with residue β45. This may affect the microtubule dynamics in the class βII-VIII isotypes [3]. Our results also conclude that the different β tubulin isotype interactions with the peptides are unique. Among the docked peptides, DRAMP00776 was found to be involved in strong interactions with αβI, αβIII, and αbIV tubulin proteins. However, it shows superior binding energy with αβI (−109.02 kcal/mol) and αβIII (−107.20 kcal/mol), compared to αβIV (−102.97 kcal/mol). Numerous preclinical studies have demonstrated that elevated levels of βIII-tubulin expression are linked with drug resistance in human cancer cell lines such as lung, ovary, prostate, and cancer [41]. On the other hand, the βI-tubulin isotype is the most highly expressed tubulin isoform in humans and the most common isotype found in cancerous cells, hence its importance as a target for inhibition in cancer cells. Chemical synthesis can readily produce the promising peptide for the respective tubulin isotypes since specific differences will result in an improved therapy protocol. Moreover, recombinant technologies can be employed to satisfy the multiple prerequisites, which are in place in the pharmaceutical sector. Peptides can be efficiently designed, functionalized, and modified to optimize their bioavailability, stability, specificity, and effectiveness to enable the peptide to fulfill clinical drug requirements [42,43].

Conclusions
An in silico analysis of the protein-peptide interactions by docking using tubulin dimers as targets was carried out in sequence with molecular dynamics. In this study, the four most important isotypes of β-tubulin were explored (βI-βIV). Reported peptides with potential anticancer properties can serve as potential candidates for developing new therapeutic options for destabilizing microtubules. Consequently, they can provide candidate structures for anticancer applications. Interfering with mitotic spindle formation is a time-tested strategy for cancer chemotherapy design and development. The molecular dynamics simulation was further explored to rule out false interactions and investigate the stability of the proteins when interacting with the selected peptides. In addition, allergenicity profiling confirmed the non-allergic properties of the peptide molecules selected in this current study. The subsequent data analysis has led to the identification of the most promising peptide molecules that bind to the residues in the tubulin dimer and thus inhibit its polymerization dynamics. The presence of minor variations of peptide binding in the protein structure of β-tubulin isoforms can be an initial step for developing a new medicinal product with high levels of selectivity and specificity for a tubulin isotype of choice. This can ultimately lead to the development of secondary lines of treatment for cancer cell types that are pharmaco-resistant after conventional chemotherapy fails due to somatic mutations, changed expression levels of target proteins, or increased expression of efflux pumps. We hope that this research will provide further insights and possibly a new focus for the future development of drugs and/or biologics that target critical regions within tubulin. We have aimed here to initiate investigations of a new generation of potential peptide-based therapeutics that would selectively and specifically target tubulin variants and thus help reduce the adverse side effects associated with a vast majority of the present cancer chemotherapy drugs.

Data Availability Statement:
The final dataset used for this manuscript can be requested from the corresponding author.