Twenty-four chemical descriptors and the biological activity of 32 1,2,4-triazine-3(2H)-one derivative as microtubule assembly inhibitors were used in the QSAR analysis (Table S1). The most effective correlation models were determined using the multi-linear regression technique using a set of non-collinear descriptors. Models that do not meet the OCDE standards and Golbraikh and Tropsha’s requirements were disqualified39. The best linear model was represented by the following equation, which included the five chemical descriptors χ, TE, NHD, LogS, and I.MLR model for pIC50 breast cancer
$$\text{pIC}50 =-10.12+1.64\times\upchi -5.69{ 10}^{-05}\times \text{TE}+1.48 \times \text{NHD}-1.37\times \text{LogS}-0.36\times \text{I}$$
(18)
Since TE, LogS, and I all have negative signs in the model equation, lowering their values may raise the value of pIC50. The pIC50 value may be improved by raising X and NHD since they both have positive signs. Based on the MLR equation, it is evident that the values of the studied activity (pIC50) are linearly correlated with the five selected descriptors. As shown in Table 1, the greatest value of the coefficient of determination (R2 > 0.6), Furthermore, the created model is predictive, as shown by the high-value R2test. The variables were selected based on statistical significance, and we excluded variables that had high multicollinearity, as shown by the lower value of MSE and the strong statistical significance of the confidence level F (Fisher’s parameter). Additionally, the robustness of the developed MLR model is shown by Q2’s value (more than 0.5).
Table 1 MLR model’s parameters.The MLR results include normalized descriptor coefficients and the correlation between observed and predicted actions (Fig. 1)40. The experimentally determined and theoretical pIC50 values displayed in Fig. 2 exhibit a strong correlation. The model consisted of five descriptors (Table 2) based on the normalization diagram of coefficients (Fig. 1), with NHD, LogS, TE, and X being the most significant.Figure 1Contribution of the descriptors in model.Figure 2The correlation between the observed and predicted activities.Table 2 Chemical descriptors, observed and predicted activities, residuals using MLR model.There is no multicollinearity among the specified descriptors, and the MLR model has strong stability, as shown by the variance inflation factor (VIF) of less than 10 for all of the created model’s selected descriptors41. The VIF values of the derived model are shown in Table 3. Therefore, it is evident that all five of the model’s descriptors have a strong correlation with pIC50. Table 2 displays the observed and predicted activity values, while their correlational relationship is illustrated in Fig. 1.
Table 3 The selected descriptors VIF.The VIF values were within an acceptable range, suggesting that multicollinearity among the selected descriptors is not a significant concern in this specific model42.Validation of the QSAR model focusing on the reliability of two-dimensional chemical descriptorsY-Randomization testA randomization test is used, and 100 models for MLR have been generated, to prevent random correlation and check the MLR model that has been built. The low values of, \({Q}_{cv \left(LOO\right)(Rand)}^{2}\) and, \({R}_{Rand}^{2}\) found for the MLR model are shown in Table S3. The findings show that the model developed is not the product of random correlation43.Golbraikh and Tropsha criteriaThe MLR model’s outcomes and Golbraikh and Tropsha’s parameters were compared. The findings in Table 4 demonstrate that the MLR model adheres to the standards set forth by Golbraikh and Tropsha. The model proposed is able to predict with high performance for new compounds.
Table 4 Comparison of model parameters (MLR) with Golbraikh and Tropsha criteria.Applicability domainBy plotting the leverage effect values (hi) vs the residual values, the application domain of the verified model was identified. The distribution of normalized residual values, the leverage level values, and the threshold leverage effect (h* = 0.56), with \({h}^{*}=3\times \frac{k+1}{n}\); k = 5; N = 32, as well as the distribution of normalized residual values and the leverage level values have calculated as shown in following William’s plot (Fig. 3). Based on the levers (h* > 0.56), none of the compounds was found outside the AD of the developed model.Figure 3Williams plot of standardized residual versus leverage for the MLR model (with: h* = 0.56 and residual limits =  ± 2.5); training samples in the black color and test samples in the red color.Using MatlabV2021a software, we define the Williams plot type’s application domain in this study.Design of new compoundsThe suggested model demonstrates how adjustments to the QSAR model’s appeared descriptors might improve activity. The finished model’s descriptors, on the other hand, may be used to interpret the data. Finding the connection between the descriptor values and the researched compounds’ structural specifics is simple.Absolute electronegativity (χ)This descriptor represents the tendency of a molecule to attract electrons. In our QSAR model, a higher value of χ corresponds to increased electron-withdrawing capacity, which enhances the interaction of the molecule with the target protein, thereby increasing the pIC50 values. This suggests that molecules with higher electronegativity are more effective in disrupting the Tubulin dynamics critical for cancer cell viability.Total energy (TE)The total energy of a molecule reflects its stability; in our study, a lower TE is associated with higher biological activity. Molecules with lower total energy are more stable and can more readily interact with the target site on the Tubulin protein, which could explain their enhanced inhibitory effects on breast cancer cell lines.Number of hydrogen bond donors (NHD)Hydrogen bond donors play a crucial role in drug-target interactions by forming hydrogen bonds with the target protein, stabilizing the drug-protein complex. Our findings indicate that an increase in NHD enhances pIC50, suggesting that additional hydrogen bonds contribute positively to the binding affinity and thus the effectiveness of the inhibition.Water solubility (LogS)Solubility is critical for the bioavailability of a drug. In our model, compounds with better solubility (higher LogS values) generally exhibited lower pIC50, indicating that while solubility is crucial for overall drug performance, excessively soluble compounds may be too hydrophilic to efficiently permeate cell membranes and reach intracellular targets.Shape coefficient (I)This descriptor quantifies the three-dimensional shape and bulk of the molecule, influencing how well the molecule fits into the target binding site. Our results suggest that a smaller shape coefficient, indicating a more compact molecular structure, is beneficial for interacting with the specific contours of the Tubulin binding site, thereby enhancing inhibitory activity.As a result, we decided to search for novel compounds with specific alterations in the structural characteristics of the compounds under study. The impact of descriptors on the biological activity and structural characteristics of the most active compounds led to the suggestion of several compounds. ChemOffice software was used to sketch and optimize the structures of the suggested compounds, and ChemOffice and Gaussian software was used to determine the effective descriptors.The structures of the proposed compounds and their corresponding predicted pIC50 values are summarized in Tables S4 and 5. The data indicate that increases in the values of X and NHD enhance the anti-cancer activity of the dataset compounds, whereas increases in the values of TE, LogS, and I diminish this activity. The addition of elements with greater electronegativity than carbon (e.g., F, Cl, O, and N) and an increase in molecular size lead to an increase in the value of X and a decrease in the value of the descriptor TE. Based on this information, modifications were made to the 1,2,4-triazin-3(2H)-one ring in the template compound 5i. This phase involved the design of 28 compounds and the calculation of their leverage values to identify outliers. Equation (18) was utilized to determine the pIC50 values. Tables S4 and 5 present the structures, values of chemical descriptors, and pIC50 values of these compounds.
Table 5 The values of chemical descriptors, the predicted activities using the MLR model.Drug‑likeness assessment and ADMET predictionsThe physicochemical characteristics of the substances are often linked to certain filter variations when assessing how drug-like they are. So, via the ADMETlab 2.0 Web server, pertinent physicochemical parameters (Table S5) are produced. According to the radar charts, the physicochemical characteristics of the fourteen compounds fall between the upper (brown) and lower (red) bounds (Table S5). The fourteen compounds (Table S6) that passed the Lipinski rule of five were further evaluated using the SwissADME Web server utilizing various drug-likeness filter rules, including the Ghose filter rule, Veber’s rule, and Egan’s rule. The results are reported in Table S6.Using Lipinski’s rule of five as a filter, the SwissADME web server was used to theoretically forecast the drug-like characteristics of the examined compounds (Table S6), including the best-hit compounds. Any small molecule that fails to meet more than one of the aforementioned requirements may have problems with bioavailability. The novel compounds’ pharmacokinetic analysis is shown in Table S6, and All of the compounds pass the Lipinski rule of five tests, including the criterion for molecular weight, which stipulates a maximum of 500 g/mol. Furthermore, all of the best hit compounds had one hydrogen bond acceptors, which was within the accepted range. For the number of hydrogen bond donors, between 3 and 8 is also within the accepted range. Additionally, the estimated LogP value never went beyond or below the acceptable bounds. The best-hit compounds were determined to be drug-like based on these filtering requirements by not exceeding the established threshold values (or by not having more than one violation of the filtering conditions used). They were all further determined to be orally bioavailable based on the value of their bioavailability score (all having 0.56)26. Given that these substances complied with all of the filtering requirements, it was generally anticipated that there would be no problems with their bioavailability.The examined compounds, including the top hit compounds, were theoretically evaluated for their ADMET/pharmacokinetic characteristics using the online web server pkCSM (Table S7). Human intestinal absorption of the best-identified hit compounds was all found to vary between 63.568 and 96.568%. These small molecules can be absorbed inside the human gut if their values of intestinal absorption in humans are larger than the minimum suggested rate of 30% given for the assessment of this attribute. The accepted threshold value for the central nervous system (CNS) permeability is > − 2 to < − 3, while that for the blood–brain barrier (BBB) permeability is > − 0.3 to < − 1. The BBB permeability for these small molecules was found to be − 1 except for compounds Pred9, Pred10, Pred11, Pred15, Pred27, and Pred28, which was > − 1. This indicates that all of the compounds, except for compounds Pred9, Pred10, Pred11, Pred15, Pred27, and Pred28, only partially infiltrate or permeate the BBB. Except for compounds Pred2, Pred9, Pred10, Pred11, and Pred27, all of the compounds had a CNS permeability value of < − 3, indicating that they only partially infiltrate or permeate the CNS.The cytochrome (CYP) is recognized as being crucial for the body’s enzymatic breakdown and small molecule metabolism. As a result, it is important to consider how these small molecules are metabolized and processed by the human body. The body’s mechanisms for breaking down and metabolizing small compounds include CYP1A2, 2C9, 2C19, 2D6, and 3A4, with CYP3A4 playing the most important role (a good small molecule is expected to be both a substrate and inhibitor of CYP3A4)26. The most popular hits were all CYP3A4 substrates and inhibitors. This led to more confirmation that the body can break down these small compounds. Excretion/total clearance, which specifies the connection/relationship between the rate of these small molecules’ elimination and concentration inside the body, is another crucial consideration. The best-hit compounds showed a higher excretion value and were tested within the acceptable range for a medication. The comprehensive toxicity and pharmacokinetic evaluation detailed in Table S7 reveals a mixed safety profile for the best-identified hit drugs. While none of the compounds displayed mutagenicity, suggesting a low risk of genetic toxicity which bolsters their overall safety profiles, several compounds, specifically Pred1, Pred2, Pred3, Pred4, Pred10, Pred12, Pred14, and Pred15, exhibited potential hepatotoxic effects. These findings necessitate further in-depth investigations, including in vitro and in vivo studies, to thoroughly understand the implications for liver health and to mitigate any adverse effects. Moreover, only Pred11 was flagged for potential carcinogenic risks, indicating a need for extended research and evaluation through long-term animal studies to ascertain the full extent of this risk. Despite these concerns, the broad ADMET properties of these compounds were typical, aligning with expected pharmacokinetic profiles and confirming that none were overtly harmful.Molecular docking studyAn Intel (R) core (TM) i5 10th Gen processor-equipped microcomputer was used to carry out the computations for molecular docking. The operating system for all apps was Windows 10, 64-bit version 2020. The software AutoDockTools version 1.5.7 was used to carry out the molecular docking33, this employs the Genetic algorithm-based trajectory modeling technique. Additionally, we used Biovia Discovery Studio version 2021 software44 to display the many interactions that have been created between the ligands and the target active site.Molecular docking has become an increasingly popular approach for the development of new drugs, in part because of the favorable time and pecuniary costs of in silico drug screening compared with traditional laboratory experiments. In this study, the fourteen designed compounds were docked, using AutoDock 1.5.7 software, with the crystal structure 1SA0 of the target protein’s binding site31. To perform docking by re-docking the co-crystallized ligand at the Tubulin enzyme Fig. 4 represented the superimposed view of docked conformation45,46,47,48,49 and the co-crystallized ligand and the RMSD value is 1.084 Å. The RMSD (Root Mean Square Distance) of the docked co-crystallized ligand was within the reliable range of 2 Å. The binding affinity of the designed compounds with the receptors 1SA0 was reported in Table 6.Figure 4Re-docking pose with an RMSD value of 1.084 Å (Green = native ligand, blue = docked ligand).Table 6 Results of the binding affinity of the most stable conformation.Microtubule protein A and B chains: structural and functional insightsMicrotubules, essential components of the cytoskeleton, are composed of α-Tubulin (A chain) and β-Tubulin (B chain) that form heterodimers. α-Tubulin forms the stable minus end of the microtubule and harbors a non-exchangeable GTP molecule, providing structural stability. In contrast, β-Tubulin constitutes the dynamic plus end, featuring an exchangeable GTP site crucial for microtubule dynamics through GTP hydrolysis post-polymerization. This differential GTP activity between the chains underpins the dynamic instability vital for cellular processes such as mitosis and intracellular transport. Given β-Tubulin’s pivotal role in these processes and its drug target accessibility, our study focuses on inhibitors that target β-Tubulin to disrupt microtubule dynamics, offering potential therapeutic strategies against rapidly dividing cancer cells.All the designed compounds exhibited docking score values between − 8.5 and − 9.6 kcal/mol, having lower binding energies than the co-crystallized ligand -7.9 kcal/mol. The compounds Pred3, Pred5, Pred12, Pred13 and Pred28 showed the best binding score with the protein Tubulin. The best conformation from Autodock resulted with good binding affinity of compound Pred28 with 1SA0 protein as   − 9.6 kcal/mol, based on Table S8 and Fig. S1, has six conventional hydrogen bonds interactions with Lys (B: 254), Gly (A: 146), Gly (A: 142), Thr (A: 145), Gln (A: 11), Asn (A: 101), and H-donor, two carbon hydrogen bonds with Lys (B: 254) and Gln (A: 11), two halogen type fluorine bonds with Glu (A: 183) and Glu (A: 142), hydrophobic interactions with Lys (B: 254), Ala (A: 12), Ile (A: 171), Pro (A: 173), and Leu (B: 248) residues.The compound Pred13 formed six conventional hydrogen bonds interactions with Gly (A: 144), Asn (A: 101), Asp (A: 69), Gly (A: 142), Lys (B: 254), and Thr (A: 145), two carbon hydrogen bonds with Lys (B: 254) and Gln (A: 11), one halogen bond with Glu (A: 183), it exhibited also hydrophobic interactions with Tyr (A: 224), Lys (B: 254), Leu (B: 248), Ala (A: 12), Ile (A: 171), and Pro (A: 173) residues.The compound Pred12 formed five conventional hydrogen bonds interactions with Gly (A: 146), Asn (A: 101), Thr (A: 145), Gln (A: 11), and Ala (A: 12), one carbon hydrogen bond with Gln (A: 11), hydrophobic interactions with Ala (A: 12), Tyr (A: 224), Leu (B: 248), Ala (A: 12), and Ile (A: 171) residues.The compound Pred5 formed five conventional hydrogen bonds interactions with Gly (A: 146), Asn (A: 101), Tyr (A: 224), Ser (A: 140), and Lys (B: 254), one carbon hydrogen bond Pro (A: 173), one halogen bond with Asn (A: 206), it exhibited also hydrophobic interactions with Gly(A: 143), Gly (A: 10), and Ala (A: 12) residues.The compound Pred3 formed six conventional hydrogen bonds interactions with Gly (A: 146), Asn (A: 101), Tyr (A: 224), Glu (A: 183), Lys (B: 254), and Ser (A: 140), three carbon hydrogen bonds with Pro (A: 173), Glu (A: 71), and Asp (A: 98), two halogen bonds with Asn (A: 206) and Tyr (A: 172), hydrophobic interactions with Gly (A: 143), Gly (A: 10), and Ala (A: 12) residues.The compound co-crystallized ligand formed two carbon hydrogen bonds interactions with Ser (A: 178) and Val (B: 315), one Pi-Sulfur bond with Met (B: 259), hydrophobic interactions with Ala (B: 316), Ala (A: 180), Lys (B: 352), Leu (B: 248), Ala (A: 180), and Lys (B: 352) residues.A significant aspect of our study lies in correlating the molecular descriptors used in the QSAR model with the binding modes observed in the docking studies. For example, the descriptor ‘χ’ (absolute electronegativity) correlates with the affinity of the ligand for the negatively charged regions of the protein active site. On the other hand, ‘LogS’ (solubility log) is reflective of the hydrophobic interactions that are crucial for the stability of the ligand–protein complex. Thus, each descriptor serves not only as a statistical variable but also as a biologically relevant parameter that can be linked to specific aspects of ligand–protein interactions.We fully acknowledge that while computational methods like QSAR modeling and molecular docking are valuable tools for predicting the biological activity of compounds, they are not a substitute for experimental validation. In an ideal scenario, the next step would involve the chemical synthesis and experimental testing of the predicted 1,2,4-triazine-3(2H)-one derivatives. Due to limitations in scope and resources, this was not possible within the confines of the present study. However, we strongly advocate for such experimental validation in future work to confirm the predicted activities of these promising compounds.Molecular dynamic simulationIn this study, Molecular Dynamics Simulations (MDS) were conducted for 100 ns to investigate the binding interactions, structural dynamics, and flexibility of Tubulin when bound to the top three hits (Pred3, Pred13, and Pred28). Key metrics, including Root Mean Square Deviation (RMSD), Root Mean Square Fluctuation (RMSF), and Radius of Gyration (Rg), were calculated from the trajectory data collected over the 100-ns period50.RMSD calculationThe root mean square deviation (RMSD) is the most important indicator of a biomolecular system’s structural stability. The computer-aided drug design community believes that a lower RMSD value signifies a more stable system. In contrast, complexes with larger RMSD values are less likely to be stable51. The RMSD graph in Fig. 5A shows that, except for the complex formed by Pred3-Tubulin, all the remaining protein–ligand entities exhibit lower fluctuations in their spectra. This implies that their conformational dynamics have been minimally perturbed during the simulation. The average RMSD values for Pred3, Pred13, Pred28, and Colchicine (with their receptors) are calculated to be 0.65 nm, 0.34 nm, 0.29 nm, and 0.34 nm, respectively. This indicates that Pred28, with its lower RMSD value, is most similar to the standard Colchicine drug and may represent the most stable drug candidate among the promising inhibitors. This finding aligns with the docking simulation results, where Pred28 exhibited the best binding affinity (Table 6). However, the RMSD measurements are insufficient to ensure the stability of the system because they do not account for fluctuations in certain sections of the structure. Therefore, an RMSF diagram was also constructed to evaluate the atomic variations of each residue along the trajectory.Figure 5The results of the molecular dynamics study: (A) Time evolution of the backbone of three best hits inhibitors and the standard (Colchicine); (B); RMSF spectra of three promising inhibitors and the standard (Colchicine); (C) The comparative Radius of gyration values for the target protein with three best hits inhibitors and the standard (Colchicine); (D) The comparative hydrogen bonds for the target protein with the reference drug (Colchicine) and three best hits.RMSF calculationFor a comprehensive description of residue stability in MD simulation, RMSF measurement is the gold standard52. In the simulation, this approach assesses protein flexibility and can identify residues with high or low flexibility. Root-mean-square-fluctuations (RMSF) were calculated for each complex to reveal which residues changed the most during the MD simulation47,53. Higher RMSF values for biomolecular system corresponds to a lower residual stability and vice versa. As demonstrated in Fig. 5B, except the Pred3, all the remaining proposed compound exhibit a very similar RMSF pattern to the standard drug. The average RMSF values are 0.31 nm for Pred3, 0.16 nm for Pred13, 0.15 nm for Pred28, and 0.15 nm for Colchicine. Notably, Pred28, which has the lowest average RMSF value, indicates it may be a better inhibitor compared to the other compounds.Radius of gyration analysisWe investigated how the binding of various ligands influenced the overall compactness of the protein’s structure48. To achieve this, we calculated the radius of gyration (Rg) as a function of time. Ligands with higher Rg values tend to be more flexible and, consequently, less stable. Conversely, lower Rg values indicate a denser and more compact conformation38. Among the complexes examined (Fig. 5C), the Pred3 exhibited the largest radius of gyration, while compound Pred28 had the smallest, indicating it is the most stable. The average Rg values for protein-Reference drug complexe and protein-designed molecule complexes (including Pred03, Pred13, and Pred28) were 3.11 nm, 3.35, 3.07 nm, and 3.06 nm, respectively. This suggests that complex Pred28 is the most the most compact biomolecular system.H-BOND analysisIn a water environment, hydrogen bonds and their relative strength play a crucial role in protein–ligand interactions, especially when the mechanism involves hydrolysis, where water is essential for breaking down a chemical54. These hydrogen bonds form when an electronegative atom from a hydrogen-bond acceptor comes into contact with a hydrogen atom directly bonded to an electronegative atom from a hydrogen-bond donor55.In this study, we examined the intermolecular hydrogen bonds formed between Tubulin and the selected compounds (Pred3, Pred13, Pred28, and Colchicine). The results, presented in Fig. 5D, show that Pred28 and Colchicine have the most significant H-bond spectra among all the simulated compounds, with average H-bonds of 1.75 and 1.94, respectively. In contrast, Pred3 and Pred13 formed fewer intermolecular H-bonds with their receptor, averaging 1.37 and 1.02 H-bonds, respectively. This finding aligns with earlier RMSD and RMSF analyses, further supporting that Pred28 is the best compound, as the H-bonding spectrum results corroborate this conclusion.
QSAR, ADMET, molecular docking, and dynamics studies of 1,2,4-triazine-3(2H)-one derivatives as tubulin inhibitors for breast cancer therapy
