Structure based drug design and machine learning approaches for identifying natural inhibitors against the human αβIII tubulin isotype

Homology modeling of human αβ_III tubulin isotype

Homology model of the human β_III tubulin isotype was constructed using the α_IBβ_IIB tubulin heterodimer (PDB ID: 1JFF) as template through Modeller 10.2⁹². Here, 100 models of the β_III tubulin isotype were generated, and the best model was selected based on the DOPE (Discrete Optimized Protein Energy) score (−51589.82422). DOPE is a statistical potential optimized for model assessment, which calculates an “energy” score for a given protein model. Lower DOPE scores generally indicate higher-quality models. The stereo-chemical quality of the β_III tubulin was assessed using the Ramachandran plot through PROCHECK⁹³ and results are shown in Supplementary Fig. 1. The Ramachandran plot revealed that 88.1% of residues are in the most favoured region, 9.6% residues in the additional allowed region, 1.5% residues in the generously allowed region, and 0.8% residues in the disallowed region. This indicates that the model quality is suitable for further molecular modelling study, as total of 97.7% of residues are in the most favoured and additional allowed region. Further, the complete model of αβ_III tubulin isotype was prepared by keeping the α_IB-tubulin constant and changing only the β_III tubulin isotype from crystal structure 1JFF.pdb, similar to an earlier study³². Hereafter, the homology modeled α_IBβ_III tubulin heterodimer is referred to as αβ_III tubulin isotype. This homology model of the human αβ_III tubulin heterodimer was then utilized for structure-based virtual screening.

Structure-based virtual screening

Structure-based virtual screening is a computational method for identifying compounds with high binding affinity by analysing their binding energy and mode of interactions with a target protein⁹⁴. Here, high-throughput virtual screening was performed by targeting the Taxol site of αβ_III tubulin isotype through AutoDock Vina⁹⁵ and compounds were filtered based on the binding affinities. A total of 1000 hit compounds were further selected with binding energy higher than − 10 kcal/mol. These hit compounds were further subjected to a machine learning approach, to identify the potential active compounds.

Identification of active compounds using machine learning

The present study employed the supervised ML approach, a widely used and highly effective method for separating active and inactive molecules from a test dataset. The virtual screening process, followed by hit identification, resulted in a total of 1000 hit compounds, which were considered as a test dataset. While, for the training dataset ‘Taxol site’ targeted and ‘non-Taxol site’ targeted drugs were considered as active and inactive datasets, respectively as shown in Supplementary Tables 1 and 2.

A total of 32 compounds were identified as active, and 26 compounds were identified as inactive for training dataset (Supplementary Tables 1 and 2). Next, DUD-E database was used to generate decoys for training dataset. Training dataset comprising 3030 decoys, delineated into 1580 active and 1450 inactive compounds. Next, chemical descriptor properties of the of test and training dataset were generated using the PaDEL Descriptor software⁵⁴and further used to build a classification algorithm to identify chemical activity. PaDEL descriptor uses the SMILES file format to calculate the descriptor properties also known as features. Initially, the prepared dataset underwent pre-processing where a binary ‘Label’ column was created based on compound ‘Activity’, subsequently eliminating rows with missing ‘pIC50’ values and metadata columns were identified⁹⁶. Subsequently, features containing zero or null values were excluded from the dataset to enhance computational efficiency. To ensure stratified sampling and preserve target class proportions, the remaining dataset was divided into training (75%) and testing (25%). We used the machine learning features of PyCaret⁹⁷ and entailed setting up a baseline for feature selection, comparing models using Area under Curve (AUC)^98,99shown in Table 1. Here, the prediction model was developed by using the following algorithms: AdaBoost, extreme gradient boosting (XGBoost), Light gradient boosting machine (LightGBM), CatBoost classifier (catboost), Gradient boosting classifier (gbc), Logistic regression (LR), Decision tree classifier (DT), Random forest (RF), Linear discriminant analysis (LDA), Extra trees classifier (ET), K Neighbors Classifier (KNN), naive Bayes (NB), Quadratic Discriminant Analysis (QDA), Dummy Classifier, SVM- Linear kernel, and Ridge classifier¹⁰⁰.

In the process of training, we used a 5-fold cross-validation approach to estimate the models. From the dataset, the precision, recall, AUC, accuracy, MCC, TT, Kappa, and F-score were estimated. Here, Table 1 shows that several models achieved very high AUC scores, including: AdaBoost Classifier, XGBoost, LightGBM. However, AdaBoost model exhibits high values for recall, precision, F1-score, Kappa, and MCC, indicating strong overall performance. Therefore, by evaluating the highest value of these indices the best prediction model, AdaBoost, was selected. Through finalization and tuning, the AdaBoost model was further refined. Visualization methods including confusion matrices, AUC plots, and feature importance analysis were performed⁹⁸. Additionally, Fig. 2 displays the receiver operating characteristic (ROC) curve for the AdaBoost model. The ROC plot (Fig. 2) shows the true positive rate (y-axis) and the false positive rate (x-axis) for each compound’s threshold between 0 and 1. Further, the ROC curve of each model was generated, and it is given in Fig. 2. Here, all curves have an AUC of 1.00, indicating that the classifier is perfectly distinguishing between the two classes. The above data undoubtedly explains the significant efficiency of each model.

The probability scores for positive and negative class predictions were visualized using a scatter plot to illustrate the model’s classification performance (Supplementary Fig. 2). The plot shows high accuracy, with most predictions aligning with the true labels, indicated by a horizontal line at the top. This suggests strong model performance and high confidence in both positive and negative class predictions. However, based on the predicted probability scores for the positive class, we selected 20 compounds exhibiting the highest activity.

Pass prediction

The PASS activity of selected 20 active compounds (see Supplementary Table 1) from the machine learning approach was further evaluated using the Way2Drug server⁶⁸. Out of 20 active compounds, ZINC03847075, ZINC12889138, ZINC08952577, and ZINC08952607 show microtubule formation inhibitors, beta-tubulin antagonistic, tubulin inhibitor, and anticancer, as well as cytostatic activity as shown in Table 2. Selected compounds are listed in Table 2, along with their PubChem IDs for further reference. This shows that the selected drug compounds have a tubulin binding activity and hence they were further considered for the ADME-T study.

ADME-T properties

The ADME and toxicity properties of ZINC03847075, ZINC12889138, ZINC08952577, and ZINC08952607 compounds were assessed using Swiss-ADME⁶⁹ and pkCSM⁷⁰ web server respectively, as shown in the Table 3. The ZINC03847075, ZINC12889138, ZINC08952577, and ZINC08952607 have molecular weight of 343.42, 487.54, 472.57, and 448.6, and iLogP value of 3.93, 4.48, 4.26, and 3.93 respectively Water solubility criteria such as Esol, Ali and Silico-IT shows that ZINC03847075 and ZINC08952577 compounds are poorly soluble, and ZINC08952607 and ZINC12889138 compounds are moderately soluble as shown in Table 3.

Lipinski rule, specifying < 5 hydrogen bond donors, < 10 hydrogen bond acceptors, showed only 1 violation for two compounds, which is acceptable. Bioavailability was predicted as 0.55, indicating neutrality. Only ZINC08952607 compound exhibited PAINS (Pan-Assay Interference Compounds) interference. Synthetic availability ranged between 1 (very easy) to 10 (difficult), with all compounds scoring < 6 as shown in Table 3. The toxicity analysis demonstrated that only ZINC03847075 compound displayed AMES mutagenesis, all the drug compounds show maximum tolerated dose (human) less than 0.5 and minnow toxicity < 0; and only ZINC08952607 compound showed hepatotoxicity. The selected compounds show appropriate pharmacokinetics activity which is helpful in further evaluations¹⁰¹.

Binding mode of αβ_III tubulin isotype with ZINC compounds using Docking

Molecular docking was performed to identify the least binding energy conformation of ZINC03847075, ZINC12889138, ZINC08952577, and ZINC08952607 with αβ_III tubulin isotype using AutoDock 4.2.3⁷². The least binding energy conformation of ZINC03847075, ZINC12889138, ZINC08952577, and ZINC08952607 was found to be −8.52, −9.71, −10.80, and − 10.57 kcal/mol, respectively (Table 4). The analysis of αβ_III-ZINC03847075 complex shows that the ZINC03847075 is stabilized by Leu273 (5.46Å), Leu361 (4.61Å), Ala231 (4.74Å), Leu361 (4.87Å), Val23 (5.43Å), Ala231 (4.60Å) and Pro358 (4.71Å) as shown in Fig. 3A; Table 4. Next, the analysis of αβ_III-ZINC12889138 complex showed that ZINC12889138 is stabilized by only the non-bonded interactions of the alkyl and π-alkyl type with Ala23 (4.16 Å), Leu215 (4.26 Å), Leu284 (3.93 Å), Pro358 (4.97 Å), Leu361 (4.41 Å), Phe270 (4.92 Å), and Leu361 (5.42 Å) as shown in Fig. 3B; Table 4. The analysis of αβ_III tubulin-ZINC08952577 complex shows that ZINC08952577 is stabilized by conventional hydrogen bonding interactions with Arg359 (3.16 Å), and Gly360 (2.85 Å), and CH type of bonding interaction with Thr274 (3.79 Å). Also, Phe270 and Leu361 form a π-donor and π-sigma type of bonding, Ala231 and Pro358 form alkyl, and Phe270, Leu361, and Leu284 form a π-alkyl type of interactions as shown in Fig. 3C; Table 4. Finally, the analysis of αβ_III-ZINC08952607 complex showed that the ZINC08952607 is stabilized by only the CH type of interactions with Ala275 (3.37 Å) and Thr274 (3.47 Å) as shown in Fig. 3D; Table 4. While Leu215, Leu217, Ala231, Pro358, and Leu361 forms alkyl-type of interactions whereas Leu215, Leu217, Phe270, and Ala275, form a π-alkyl type of non-bonded interactions as listed in the Table 4.

The analysis docking complexes revealed that only the non-bonded type of interactions stabilizes the αβ_III-ZINC12889138 and αβ_III-ZINC03847075 complexes. While, the αβ_III-ZINC08952577 complex is stabilized by both the type of interactions such as bonded and non-bonded interactions. Similarly, the αβ_III-ZINC08952607 complex is stabilized by the CH and non-bonded interactions as shown in Fig. 3; Table 4. Hence, to refine the binding mode and interaction, molecular dynamics simulation was employed and discussed below.

Molecular dynamic simulation

To investigate the refined mode of interaction, stability and affinity of αβ_III tubulin isotype with natural drug compounds, MD simulations were performed using Gromacs2021.5⁷⁴. The least binding energy docked complex of αβ_III-ZINC03847075, αβ_III-ZINC12889138, αβ_III-ZINC08952577, and αβ_III-ZINC08952607 were used as a starting conformation for MD simulations (Fig, 3).

The stability of the MD simulation trajectories was examined by calculating the root mean square deviations (RMSD) of the C_α backbone atoms of αβ_III tubulin isotype. RMSD plot revealed that all αβ_III tubulin-drug complexes reached equilibrium after 200ns and remained stable, with fluctuations between 0.20 ± 0.40 nm as shown in Fig. 4A and structural dynamics of all the complexes are shown in Supplementary Movie 1–5. Furthermore, RMSD plot reveals that αβ_III-ZINC12889138 complex has less RMSD fluctuations (below 0.3 nm) and is found to be highly stable compared to αβ_III tubulin and αβ_III tubulin-drug complexes as shown in Fig. 4A. Additionally, RMSD of the drug compounds were plotted to visualize its conformational changes in bound state over time (Supplementary Fig. 3). All drug compounds exhibited deviation patterns, with ZINC12889138 and ZINC03847075 showing more stable fluctuations compared to the others, as illustrated in Supplementary Fig. 3. Additionally, the RMSD analysis revealed that ZINC12889138 and ZINC03847075 maintained lower fluctuations (below 0.3 nm), reflecting more stable binding throughout the simulation. In contrast, ZINC08952577 and ZINC08952607 showed persistent deviations, indicating relatively unstable interactions.

Moreover, to understand the impact of drug on the protein conformation, root mean square fluctuations (RMSF) of backbone C_α atoms were calculated. Regions with higher RMSF values correspond to greater flexibility, whereas lower RMSF values indicate more rigid or stable regions. As shown in Fig. 4B, regions in proximity to the Taxol binding site-specifically the M-loop (residues 272–287) and the H6-H7 loop (residues 210–240)-exhibit notable changes in dynamics. The RMSF analysis indicates that the M-loop becomes increasingly flexible upon drug binding to β_III tubulin. Similarly, the H6-H7 loop also showed enhanced dynamic behaviour. In addition to these, other regions such as the T5 loop (residues 170–180) and a segment of helix H11 (residues 405–426) display a significant increase in flexibility. The region comprises the H6 helix, T7 loop, H7 helix, S7 sheet, and M-loop, which collectively play a key role in accommodating drug compounds at the Taxol binding site (Supplementary Movie 1–5). The analysis of RMSF plot revealed that β_III-ZINC12889138 complex show lower RMSF value (specifically in the region; 225–375 amino acid) as compared to other β_III and drug complexes. Moreover, ZINC12889138 and ZINC08952577 compound has been found to be stable at the ‘Taxol site’ of βIII tubulin and showing the significant interactions with the M-loop (Supplementary Movie 3–4), compared to ZINC03847075 and ZINC08952607 (Supplementary Movies 2 & 5). Hence, to further understand the effect of drug on the protein’s structural stability and conformation the radius of gyration (Rg) and solvent accessible surface area (SASA) were calculated (Fig. 4C and D).

Radius of gyration (Rg) value has been utilized to get insight into the stability of a system by indicating the degree of compactness of a protein. Rg plot analysis revealed that αβ_III-ZINC12889138 and αβ_III-ZINC08952577 complex has less Rg value and more compact structure compared to αβ_III-tubulin, αβ_III-ZINC03847075 and αβ_III-ZINC08952607 complexes as shown in Fig. 4C. Similarly, αβ_III-ZINC12889138 (~ 300–325 nm² and αβ_III-ZINC08952577 (~ 300–340 nm² complexes had a decrease in surface area as shown in Fig. 4D, while αβ_III and αβ_III with other drug complexes showed higher SASA value (Fig. 4D).

Overall, MD analysis revealed that the ZINC12889138 and ZINC08952577 have profound effect on the αβ_III tubulin heterodimer as compared to other drug compounds. Furthermore, to explore effect of drugs on αβ_III tubulin heterodimers, principal component analysis (PCA), free energy landscape (FEL) and binding energy calculations were performed and are discussed in the below sections.

Principal component analysis (PCA)

PCA is an important method to find the total combined movements of the C_α atoms in a protein, which are represented by the Eigen Vectors (EV) s of the covariance matrix. In this study, PCA was performed using the Cartesian coordinates of the 500 ns MD simulation trajectory to assess the conformational dynamics of the αβ_III tubulin heterodimer and its complexes with selected compounds, as illustrated in Fig. 5.

The first two principal components (PC1 and PC2), which represent the dominant motions, were extracted and plotted. The αβ_III tubulin heterodimer complex with αβ_III-ZINC03847075 exhibits wider diversity (PC1 range ~ −4 to 6 nm; PC2 ~ −3 to 4 nm) of conformations during the simulation. Among the complexes, the αβ_III-ZINC08952577 (PC1 range ~ −4 to 6 nm; PC2 ~ −2 to 5.5 nm) and αβ_III-ZINC08952607 (PC1 range ~ −4 to 6 nm; PC2 ~ −3 to 3 nm) also show wider distribution, indicating significant conformational flexibility in these complexes. While αβ_III-ZINC12889138 complex (Fig. 5) exhibits the constricted clustering of data points along both PC1 and PC2 axes (PC1 range ~ −3 to 5 nm; PC2 ~ −3 to 2 nm), indicating limited structural fluctuations and enhanced rigidity. This suggests that the interaction with ZINC12889138 significantly restricts the dynamic movements of αβ_III tubulin, likely favouring a more stable conformation. Overall, the PCA analysis revealed that ZINC12889138 induces notable structural stabilization of αβ_III tubulin, as reflected by its minimal conformational variation in comparison to other ligand-bound forms. This rigidity may play a crucial role in its effectiveness as a potential inhibitory compound targeting the αβ_III tubulin isotype.

Free energy landscape

To explore the conformational dynamics and stability of αβIII tubulin in both its unbound form and when complexed with the ZINC03847075, ZINC12889138, ZINC08952577 and ZINC08952607, FELs and their corresponding contour maps were generated and analysed (Fig. 6). FELs allowed us to visualize the conformational state and identify potential energy barriers for transitions between the distinct structural configurations. Figure 6A and E illustrate the conformational ensembles of the various systems, where deep purple areas in the FEL indicate the most energetically favourable states (global minima) whereas other colors indicates the energetically unfavourable states (local minima). The unbound αβ_III tubulin (Fig. 6A) demonstrates a broad distribution of conformations, requiring more extensive sampling to reach its lowest energy state, which lies within the range of 0 to 1.740 kJ/mol. In comparison, the αβ_III-ZINC03847075 complex (Fig. 6B) reaches its energy minima between 0 and 1.820 kJ/mol, while αβ_III-ZINC08952577 (Fig. 6D) and αβ_III-ZINC08952607 (Fig. 6E) display minima in the range of 0 to 1.710 kJ/mol (Fig. 6D and E). These complexes exhibit a single dominant energy basin but show greater variability in their folding transitions. Notably, αβ_III-ZINC12889138 complex (Fig. 6C) shows a narrower conformational distribution, with a global energy minimum between 0 and 1.680 kJ/mol, indicating enhanced structural stability throughout the simulation. The FEL profile of this complex suggests that it maintains a more stable conformation over time compared to the other ligand-bound forms and maintain the folding transition. Hence, to gain a deeper understanding of intermolecular interactions such as the bonded and non-bonded interactions between αβ_III tubulin and the drug compounds, we further employed binding energy calculations using the MM-GBSA method.

Pair distribution functions (PDFs) analysis

Pair Distribution Functions (PDFs) provide a spatial distribution profile of drug molecules relative to a reference point during molecular dynamics simulations, offering crucial insights into molecular packing, interaction strength, and structural stability. As illustrated in Fig. 7, each drug compound exhibits a unique interaction profile. Notably, the PDF for the αβ_III-ZINC12889138 complex shows the highest and most distinct first peak, centered around ~ 7 Å with a g(r) value of ~ 4.2 as shown in (Fig. 7). This sharp and prominent peak reflects a strong local ordering and a high likelihood of interaction at this specific distance, suggesting a robust and well-defined structural arrangement. The breadth of the peak also implies a degree of conformational flexibility in the interaction range. Beyond this point, the PDF oscillates before stabilizing near unity at larger distances, indicating persistent, long-range interactions (Fig. 7).

In contrast, the αβ_III-ZINC08952577 and αβ_III-ZINC08952607 complexes display first peaks around ~ 6 Å with peak g(r) values near ~ 3.3. Although less intense than ZINC12889138, these peaks still reflect moderate interaction strengths (Fig. 7). The αβ_III-ZINC03847075 complex shows the weakest interaction, with a broad and relatively low-intensity first peak around ~ 6 Å and a peak height of approximately ~ 5. This indicates a lower probability of interaction and weaker structural ordering compared to the other systems. Overall, among the compounds studied, ZINC12889138 exhibits the strongest and most stable interactions with αβ_III tubulin, while ZINC03847075 demonstrates the weakest, as supported by its low and diffuse PDF profile (Fig. 7).

Binding energy calculations

The binding energy calculations were performed using the MM/GBSA approach through gmx_mmpbsa tool⁹⁰ to quantitatively clarify the energetics of the binding of αβ_III tubulin with ZINC03847075, ZINC12889138, ZINC08952577 and ZINC08952607. The estimated binding energies (ΔE_bind) of αβ_III tubulin with ZINC03847075, ZINC12889138, ZINC08952577, and ZINC08952607 was found to be −0.02, −44.88, −29.62, −28.38 kcal/mol, respectively (Table 5). Among the analyzed complexes, αβ_III tubulin exhibited the strongest binding affinity with the ZINC12889138 compound, while the lowest binding affinity was observed with ZINC03847075. (Table 5). The binding affinity of the αβ_III tubulin towards the selected compounds followed a descending order: ZINC12889138 > ZINC08952577 > ZINC08952607 > ZINC03847075. Furthermore, In the αβ_III– ZINC12889138, the van der Waals energy makes highest energy contribution, while in case of αβ_III-ZINC08952577 both electrostatic and van der Waals make a contribution in the binding as compared to other complexes (Table 5).

The solvation energy (ΔE_sol) plays a significant role in ligand binding by offsetting the gas-phase binding energy (ΔE_gas) to varying extents. For αβ_III-ZINC03847075, a significant reduction in the van der Waals interaction (−0.03 kcal/mol), a high ΔE_sol value (143.28 kcal/mol) which reduces the strong gas-phase interaction, resulted in less binding energy (ΔE_bind = −0.02 kcal/mol) while maintaining stability. In the case of αβ_III-ZINC12889138, minimal solvation effects (ΔE_sol = 7.48 kcal/mol) allow the gas-phase interactions to dominate, yielding a stronger binding energy (ΔE_bind = −44.88 kcal/mol) indicating hydrophobic interactions. For αβ_III-ZINC08952577, ΔE_sol (20.70 kcal/mol) partially offsets gas-phase energy, leading to balanced and moderate binding energy (ΔE_bind = −29.62 kcal/mol). In αβ_III-ZINC08952607, substantial solvation effects (ΔE_sol = 321.15 kcal/mol) nearly neutralize gas-phase interactions, resulting in a moderate binding energy (ΔE_bind = −28.38 kcal/mol) and adaptability within the solvent-protein environment. The net binding free energy, which is decided by the competition of Egas and Esol, is lowest for αβ_III-ZINC12889138 and αβ_III-ZINC08952577 as shown in Table 5. As a result, the binding free energy calculation confirms the findings of MD simulation investigations indicating that the ZINC12889138 compound has highest binding affinity with αβ_III-tubulin isotype compared to other drug compounds.

The per-residue energy decomposition analysis highlights the key binding site residues involved in stabilizing the interaction between αβ_III tubulin and various drug candidates. In the αβ_III tubulin–ZINC03847075 complex (Fig. 8A), Gly360 and Thr274 emerge as significant contributors to ligand binding, indicating their active role in the interaction. For the αβ_III tubulin–ZINC12889138 complex (Fig. 8B), residues such as Leu228, Thr232, Phe270, and Leu273 show the highest energy contributions. Additionally, Leu215, His227, Ala231, Met300, and Leu361 also participate in stabilizing the complex, reinforcing the involvement of a broader binding interface within the αβ_III tubulin isotype. In the αβ_III tubulin–ZINC08952577 complex (Fig. 8C), prominent energy contributions are observed from Phe270, Pro272, Thr274, Arg276, and Leu361, along with Ala231 and His227, suggesting a strong network of interactions with the ligand. The αβ_III tubulin–ZINC08952607 complex (Fig. 8D) features key contributions from residues Val23, Asp26, Ala231, Phe270, Arg359, Gly360, and Leu361, indicating their role in initiating and maintaining the binding process. Among all the complexes, the αβ_III tubulin–ZINC12889138 complex (Fig. 8B) demonstrates the most diverse and energetically favorable set of interactions, involving multiple residues throughout the binding pocket. This extensive engagement likely accounts for its enhanced binding affinity compared to the other drug-tubulin complexes illustrated in Fig. 8A and C, and 8D.