3D-QSAR, molecular docking, ADMET, simulation dynamic, and retrosynthesis studies on new styrylquinolines derivatives against breast cancer

The sybyl program was used to draw the structure of the compounds [31], and the geometry of these compounds was optimized via a tripos force field with a convergence threshold of 0.01 kcal/mol [32]. Using the alignment technique, the minimized structures were employed as the main conformations to create a trustworthy 3D-QSAR model [33–35]. Compound 22 was used as a template to align all 43 compounds, as shown in Figure 2.

Figure 2

Alignment of molecules.

Comparative molecular similarity indices analysis (CoMSIA) contours were used to construct 3D-QSAR models to explore the relationship between 3D structure and the biological activity of the designed molecules [33,36,37]. CoMSIA contour used steric, electrostatic, hydrogen bond donor (HBD), hydrogen bond acceptor (HBA), and hydrophobic fields to enhance activity [38]. The energy cutoff was set at 30 kcal/mol, and the column filtering value was set at 2.0 kcal/mol, but the rest parameters of the partial least square (PLS) analysis were launched by default. There are parameters to consider while choosing a better model, such as the high values of correlation coefficient R ² and the coefficient of cross-validation correlation Q ² (R ² > 0.60 and Q ² > 0.50) [39,40]. Other parameters are taken into consideration such as the optimum number of components (N), the standard error of estimate (SEE), and F-test value (F). Furthermore, based on a test set, r ext 2 {r}_{\text{ext}}^{2} parameter was utilized to determine the best predictive model [41].

The Y-randomization method was used to check the correlation between the molecule’s structures and their pIC₅₀. The computation of the two parameters q Yrand 2 {q}_{\text{Yrand}}^{2} and r Yrand 2 {r}_{\text{Yrand}}^{2} after each random change in the pIC₅₀ is sufficient to assess the model’s dependability [38,42].

The energies of the boundary orbits [43], the lowest unoccupied molecular orbitals (LUMO), and the highest occupied molecular orbitals (HOMO) (E _LUMO, and E _HOMO) for the proposed compounds were calculated using Gaussian 09 software [44] with DFT and conjugated with the B3LYP/6-31G base [45,46]. E _HOMO and E _LUMO are used to calculate the following properties; the energy gap ∆E _gap = (E _LUMO − E _HOMO), chemical hardness (η = (∆E _gap)/2), the chemical potential μ is the first derivative of total energy [47], it is calculated by (μ = (E _LUMO + E _HOMO)/2), electronegativity (χ = −μ) explains the possibility and character of bond polarization [48], the chemical softness (S = 1/η), and the global electrophilicity (ω = (μ ₂/2η)) [49].

In recent decades, extensive research has been conducted on using tubulins as a biological target in cancer treatment, mainly due to their crucial role in the mitosis process [50]. Various tubulin inhibitors, specifically αβ-tubulin, have been investigated, revealing four main classes of compounds that bind to specific and corresponding sites: laulimalide, taxane, Vinca alkaloid, and colchicine. For docking with potent molecules, the protein tubulin complex with colchicine (PDB ID: 4O2B) was selected as the ideal target protein [51]. Downloading the protein tubulin (4O2B) from the Protein Data Bank (PDB) is the first step in molecular docking by link (https://www.rcsb.org). After that, it is necessary to eliminate water molecules to prepare receptors in their fundamental state. This step is followed by the usage of Surflex-Dock module in Sybyl 2.0, and in the final step, it is necessary to use Discovery and Pymol software [52,53] to determine the types of interactions [54,55].

The pharmacokinetics properties of absorption, distribution, metabolism, excretion, and toxicity are abbreviated in the word ADMET [56,30,38]; the SWISSADME web server has been programmed to determine these properties of the proposed derivatives (M1–M4) and 22; SWISSADME is a free web application on a website (http://www.swissadme.ch), and this step requires in the first step the drawing of these molecules, then the transformation of these structures in the form of codes (SMILES format), pkCSM (https://biosig.lab.uq.edu.au/pkcsm) webservers used to provide a set of molecules (SMILES format) [30,38], and in the second step, we put the codes (SMILES format) in the prediction of ADMET properties [54,57,58].

MD simulations during 100 ns were applied to get an overview of the stability of the [Tubulin–PredM1; Tubulin–PredM2] protein–ligand interaction [59]. The dynamic simulation run was performed based on the best docking position. In order to neutralize the system, Cl⁻/Na⁺ ions are introduced in this stage. The steepest descent (5,000 steps) was applied to minimize system energy. To check if the system has reached the pressure (1 bar) and temperature (300 K), balancing is applied under an NVT ensemble (constant number of particles N, constant system volume V and constant temperature T) at 298 K. A MD simulation with an NPT ensemble (constant pressure P) was followed to relax and refine the configuration at 1 atm, to stabilize the pressure and temperature of the complexes. Each complex is simulated at a temperature of 300 K for 125 ps, with the use of 40 and 400 kJ mol⁻¹ nm⁻² positional restraints on the side chains and the backbone, respectively. In addition, the NPT assembly was used in the production of 100 ns of simulation. The production and equilibration run were obtained during the simulation time by GROMACS 5.0.7 [60], and the complexes [Tubulin–PredM1; Tubulin–PredM2] are subjected to a CHARMM36 force field [61]. The most accentuated descent approach (5,000 cycles) was used for initial minimization by the simulation annealing method. Energy minimization for the entire system is performed at 50,000 energy minimization steps. The Ewald particle mesh method deals with electrostatic interactions with a physiological condition (0.9% NaCl, pH 7.4, 300 K) [62]. After the dynamic simulation duration (100 ns) is complete, the stability of [Tubulin–PredM1; Tubulin–PredM2] complexes was analyzed by root mean square fluctuation and deviation (RMSF and RMSD), number of H bond, solvent accessible surface area (SASA), and radius of gyration (R _g). The g-MMPBSA program was used to estimate the binding free energies (ΔG _Bind) of the screened complexes from the molecular mechanics of the Poisson–Boltzmann surface (MM-PBSA).

Generally, the binding free energy can be expressed as Δ G binding = Δ G complex − ( Δ G protein + Δ G ligand ) , \Delta {G}_{{\rm{binding}}}=\Delta {G}_{{\rm{complex}}}\left-(\Delta {G}_{{\rm{protein}}}\left+\Delta {G}_{{\rm{ligand}}}), where Δ G complex {\rm{where}}\Delta {G}_{{\rm{complex}}} represents the total free energy of the complex, Δ G protein \Delta {G}_{{\rm{protein}}} represents the total free energies of the isolated protein, and Δ G ligand \Delta {G}_{{\rm{ligand}}} represents the total free energy of ligands in a solvent. To decompose the binding energy, at first Δ E MM \Delta {E}_{{\rm{MM}}} , Δ G polar , \Delta {G}_{{\rm{polar}}},\hspace{1em} and Δ G non − polar \Delta {G}_{{\rm{non}}-{\rm{polar}}} were separately calculated for each residue.

Retrosynthesis is an essential technique for the organic synthesis of drugs. It facilitates the synthesis of candidate drug molecules by reducing the cost and time of the synthesis process [63]. The current work (spaya.ai) database was used to facilitate the synthesis of compounds M1 and M2 as drug candidates against breast cancer.

The Domain of Applicability determines the range in which a compound can be confidently predicted. This analysis aspect is explicitly requested [64]. It should be noted that the QSAR-built model is not intended to be used outside of this domain [65]. There are multiple ways to determine the scope of the quantitative structure-property relationship model, including the “leverage” method. The standardized residuals of the dependent variable vary depending on the leverages, which is the basis of this method. Suppose a compound has leverage that exceeds the threshold h* = 3(k + 1)/N (where N is the number of molecules in the learning set, and k is the number of descriptors). In that case, that compound is considered adequate for the developed model.

The accuracy of the 3D-QSAR model underwent both internal and external validation. CoMSIA produced the most favorable results with an R ² value of 0.97, an optimized component of 5, a Q ² value of 0.84, an SEE value of 0.09, and an F-test value of 224.65. Additionally, the SEE value was low at 0.09, demonstrating the created model's strong internal predictability. Hydrophobic, electrostatic, steric, HBD, and HBA contributions were 20, 9, 6, 11, and 54%, respectively. The external validation ( r ext 2 {r}_{\text{ext}}^{2} ) of CoMSIA reached 0.91. The HBA field has a remarkable role in the CoMSIA model, followed by the hydrophobic field. The statistical results of 3D-QSAR PLS analysis are summarized in Table 2, and the experimental and predicted pIC₅₀ of molecules are listed in Table 3. The values of experiments and predictions are closer. Figure 3 shows a regular distribution of styrylquinoline derivatives that target human breast cancer activity values as tubulin polymerization inhibitors.

Table 2

Model	R 2	Q 2	F	S CV	r ext 2 {r}_{\text{ext}}^{2}	N	Fraction
							Ster	Acc	Elec	Hyd	Don
CoMSIA	0.97	0.84	224.65	0.09	0.91	5	0.06	0.54	0.09	0.20	0.11

Figure 3

Graphical representation of the observed and predicted activity of training and test sets for CoMSIA models.

The green area of the steric field shown in Figure 4a shows a steric contribution to potency, while the yellow zone decreases potency. A green contour region is visible near two (O–CH₃) groups of R ₂ substituent, suggesting a preference for small bulky groups that explains the greater biological activity of compound 17 (pIC₅₀ = 5.72) than compounds 16 (pIC₅₀ = 5.47) and 15 (pIC₅₀ = 5.15). Figure 4b shows an electrostatic contour map; the blue contour represents the electronegative charge disfavoring, and the red contour represents the electropositive group disfavoring, respectively. The presence of blue near two (O–CH₃) of R ₂ and red near the other group (O–CH₃) explains why compound 29 (pIC₅₀ = 5.65) is more active than compound 28 (pIC₅₀ = 5.35) and compound 27 (pIC₅₀ = 5.27). Based on electrostatic and steric fields percents (9 and 6%, respectively), HBD, the hydrophobic fields, and HBA in CoMSIA were examined, as shown in Figure 5a, b, and c, respectively.

Figure 4

Contour maps of CoMSIA analysis: (a) steric field, and (b) electrostatic field.

Figure 5

Contour maps of CoMSIA: (a) HBD field, (b) hydrophobic field, and (c) HBA field.

HBD (Figure 5a) has a lower stake (11%) than HBA and hydrophobic fields, making them suitable for design in other fields. The cyan and purple in the HBD field indicate favorable and unfavorable conditions for increased activity, respectively. Where H-bond donor groups cause increased activity, as shown by a cyan area around the three –O–CH₃ groups of substituent R ₂, which explains the higher activity of compound 11 (pIC₅₀ = 5.19) compared to compound 8 (pIC₅₀ = 5.16), and the higher activity of compound 18 (pIC₅₀ = 5.36) compared to compound 13 (pIC₅₀ = 4.85). Figure 5b demonstrates the hydrophobic field; the white and yellow colors represent the hydrophobic parties decrease and increase activity. The R ₁ and R ₂ substituents are covered in yellow, meaning water repellency is preferred for increasing activity at these locations. For example, the higher activity of compound 34 (pIC₅₀ = 5.26) compared to compound 8 (pIC₅₀ = 5.16). While the magenta and red colors are shown in Figure 5c of the HBA contour map, they represent HBA groups that can increase and decrease biological activity, respectively. Two substituents, R ₁ and R ₂, with (O–CH₃) covered by magenta, show that this type of HBA group is preferred for increasing activity at these locations, which also explains the higher bioactivity of compound 40 (pIC₅₀ = 5.75) compared to the bioactivity of compound 21 (pIC₅₀ = 5.73).

It provides a further understanding of the relationship between the pIC₅₀ and the structures of the molecules. The value of pIC₅₀ changed randomly to determine the q Yrand 2 {q}_{\text{Yrand}}^{2} and r Yrand 2 {r}_{\text{Yrand}}^{2} values after each Y-randomization test; the results are listed in Table 4, where the lower values indicate that the correlation is not random.

Table 4

Iteration	CoMSIA		Iteration	CoMSIA
	q Yrand 2 {q}_{\text{Yrand}}^{2}	r Yrand 2 {r}_{\text{Yrand}}^{2}		q Yrand 2 {q}_{\text{Yrand}}^{2}	r Yrand 2 {r}_{\text{Yrand}}^{2}
1	0.234	0.367	4	0.345	0.356
2	0.453	0.543	5	−0.052	0.121
3	−0.234	0.053	6	0.123	0.321

To design new compounds with greater activity than the most active molecule [22], the 3D-QSAR model uses the most important fields to investigate suitable substituents of styrylquinoline derivatives as tubulin polymerization inhibitors against breast cancer cells. The most contributed fields (HBA) of four newly designed styrylquinoline compounds (M1–M4) have the most tubulin polymerization inhibitory activity, and the pIC₅₀ value of compound 22 is lower than the pIC₅₀ of the newly designed compounds. Table 5 shows the structures and calculated pIC₅₀ of the newly designed compounds.

Table 5

N°	Compounds	Predicted pIC50	N°	Compounds	Predicted pIC50
22		5.80	M3		6.42
M1		6.76	M4		6.14
M2		6.65

Molecular frontier orbitals can donate and accept electrons. HOMO and LUMO represent the capacity of a compound to donate and withdraw electrons from proposed molecules (M1–M4) electrons, respectively. The unit of energy was measured in the atomic units (au) and then converted to eV (1 au = 27.2113 eV). Table 6 presents the properties computed for proposed compounds in this article, such as hardness η, electronegativity χ, chemical potential µ, energy gap ∆E _gap, and electrophilicity ω. The energy gap (∆E _gap) is an important property of designed compounds, giving essential information about these compounds [33]. The smaller the gap energy, the more easily the compound can be excited. The results outlined in Table 6 demonstrate that M1 and M2 have energy gaps of 3.1382 and 3.0604 eV, respectively, indicating that these compounds are less stimulated. This information can be valuable when considering the properties of these compounds for potential applications. The electronegativity and hardness of the proposed compound M1 were found to be 3.9577 and 1.5610 eV, respectively. The LUMO and HOMO energies of the less excited molecule M1 were calculated as −2.3886 and −5.5268 eV. In contrast, the LUMO and HOMO energies of the designed molecule M2 were calculated as −2.1840 and −5.2444 eV, respectively.

Table 6

N°	Global properties (eV)
	E LUMO	E HOMO	µ	η \hspace{0.25em}\eta	χ	ω \omega	∆E gap
M1	−2.3886	−5.5268	−3.9577	1.5610	3.9577	4.7436	3.1382
M2	−2.1840	−5.2444	−3.7142	1.5302	3.7142	4.5076	3.0604
M3	−2.4756	−5.2085	−3.8420	1.3664	3.8420	5.4014	2.7329
M4	−2.4337	−5.2403	−3.8370	1.4033	3.8370	5.2456	2.8066
22	−2.1739	−5.5301	−3.8520	1.6781	3.8520	4.4210	3.3562

αβ-Tubulin dimers have four distinct interaction sites, including the colchicine-binding site (CBS). This makes CBS a promising target for developing new tubulin modulators. The tubulin co-crystal ligand was redocked into the binding pocket before the docking investigation to ensure that the docking approach was reliable, as shown in Figure 6. The redocked conformation of 4O2B was extremely consistent with co-crystallized ligand conformation with an RMSD of 1.388Å (<2.0Å). That indicates that the docking method was reliable and might be used in future research. Molecular docking was used to examine how the newly designed compounds (M1, M2, M3, and M4) and compound 22 bind to the target tubulin–colchicine complex (PDB ID: 4O2B). The results of the 2D view are listed in Table 7, and the interaction type and total score are listed in Table 8.

Figure 6

Re-docking pose with the RMSD value of 1.388 Å (green = original, orange = docked).

Table 7

N°	2D view	N°	2D view
22		M1
STRM2		M3
M4

The results show that all compounds have favorable interactions, with a total score between 3.45 and 6.53. Compound M1 is the most stable because it has the highest total score of 6.53. The interaction types of M1 with the tubulin protein active site consist of two conventional hydrogen bonds: Asn167 and Ala354, six bonds mixed between carbon H-bonds and Pi-DHB: (Asn258), Tyr202, Ala317, Thr353, Lys352, five Pi sigma/amide-Pi stacked: Leu242, Val238, Leu255, Lys254, and six Pi-alkyl types: Leu252, Ala250, Leu255(4.94), Ala354. Compound M2 has a total score of 5.74 with six different interaction types with the tubulin protein active site consisting of carbon H-bond and Pi-DHB: Val238, Gly237, Asn258, Lys352, one of amide-Pi stacked Lys254, one of Pi-cation Lys254, and five of Pi-alkyl Leu255, Cys241, Ile318, Ala250, Lys352. Moreover, compound M3 has a total score of 4.43 and has four carbon–hydrogen bonds: Gln247, Asn249, Asn350, VaL238, one of amide-Pi stacked Leu248; one of Pi-cation type Met259; and three of Pi-Alkyl Ala316, Gys241, Leu248. The binding mode of compound M4 revealed that it bonds inside the binding pocket with a total score of 3.67 and has four types of carbon–hydrogen bond interaction Thr353, Gln247, Asp251, one Pi-Sigma (Leu255), one Pi-cation (Lys254), and three types of pi-alkyl interaction Ala250, Lys352, and Leu 248. Compound 22 has a total energy of 3.45 with three types of carbon–hydrogen bond interactions: Asn349, Asn350, one interaction of type Pi-sigma Leu255, one interaction of type Pi-cation Lys254, and three interactions of type Pi-Alkyl Ala250, Lys352.

These compounds are more stable than the reference compound 22, which has a total score of 3.45. It has two interactions of type carbon–hydrogen bond interaction (Asn349, Asn350), one interaction of type Pi-sigma (Asn258), and two mixed interactions between Pi-alkyl and alkyl (Leu248, Cys241). Additionally, molecular docking studies indicated that compounds M1 and M2 have a strong affinity for the CBS of tubulin and could serve as lead molecules for future pharmaceuticals.

Computational pharmacology (ADMET) is a rapidly growing field that includes developing software to integrate and capture medical and biological data worldwide. It is important to study the ADMET properties of compounds, which can be done using in silico techniques, to find and create new medicines. Table 9 shows the ADMET prediction. The results show that all the compounds have an absorbance value of almost 90%, which means the human intestine absorbs them well. The blood–brain barrier (BBB) protects the brain and maintains good health. Furthermore, it blocks the flow of hazardous chemicals, preventing them from easily entering the brain; all designed compounds have a log BB value <−1 [66,67], and they cross the BBB poorly. The second parameter is the volume of distribution (VDss), which measures the intensity of drug diffusion in the body; the values are not accepted if (VDss) <−0.15 and are considered low, while compounds with VDss >0.45 are accepted and have a high distribution. All designed compounds (M1, M2, M3, and M4) have VDss >0.45.

Drug detoxification and elimination is an essential stage in metabolism; many cytochromes are involved, including 2D6, 3A4, 1A2, 2C19, 2C9, 2D6, and 3A4. All newly designed compounds (M1, M2, M3, and M4) were shown to be CYP3A4 inhibitors and substrates, but they are inhibitors only of 1A2 and 2C19 in this research. Constant clearance describes the relationship between drug concentration in the body and drug removal rate; when the clearance value is higher, the medication is considered safe. All designed compounds have a higher clearance value and are not toxic, making them the best-designed drug treatment against human breast cancer cells.

Under normal pressure and temperature conditions, the MD simulations were run for 100 ns for the best-newly designed compounds M1 and M2 with tubulin protein active site (PredM1 and PredM2) complexes on the “C-alpha position.” They demonstrated the best stability at molecular docking, and the results during dynamic simulation showed that the ligands (PredM1 and PredM2) remain bound to the protein (tubulin (pocket. Calculations of RMSD, RMSF, R _g, hydrogen bonding, SASA, and average center of the mass distance between protein–ligand and binding free energy (MMPBSA) were carried out to assess the stability of each structure.

All the proposed complexes [Tubulin–PredM1; Tubulin–PredM2] showed good stability in 100 ns simulation.

All complex configurations’ stability was studied using the RMSD during 100 ns MD simulation. Figure 7 shows RMSD for (PredM1 and PredM2); the RMSD of PredM1 ligand reaches its maximum value but is less than 40 ns during the first 30 s; after 40 ns, the graphs of M1 ligands and the protein need clarification during simulation. In the pred M2, generally, the whole chart remains stable during simulations (100 ns), and the RMSD values for the complex reach 2 Å. Generally, backbone RMSD stabilizes after 40 ns for the (Tubulin–PredM1) complex and after 10 ns for the (Tubulin–PredM2) complex. The average values of the backbone RMSD for M1 and M2 are 1.9 ± 0.3 and 1.3 ± 0.17 Å, respectively. Generally, all deviations of proteins, ligands, and complexes are less than 4 Å; this indicates the proteins, ligands, and complexes are stable in the majority time of simulation dynamics.

Figure 7

RMSD graphs of the complex, protein (Tubulin), and ligands (M1, M2) according to the simulation.

RMSF was utilized to evaluate the average fluctuation of the overall structure for each amino acid in the protein or to identify a flexible protein region. Figure 8 shows the RMSF graphs for the two complexes (Tubulin–PredM1, Tubulin–PredM2). In the graph of Pred M2, five peaks appear at the residue numbers 50, 180, 280, 380, and 440, while in the complex Pred M1, a peak appears at the 380 residues. As a result, the (Tubulin–PredM1) complex is more adaptable and of higher quality than the (Tubulin–PredM2) complex. The fluctuation intensity remains below 3.0 Å for both complexes except for some residues representing loops or turns in the protein (Tubulin). The (Tubulin–PredM2) complex has lower RMSF values than (Tubulin–PredM1), which indicates a better effect the ligand has on stabilizing protein (Tubulin) residues than M1. Generally, the RMSF of all complexes (Tubulin–PredM1, Tubulin–PredM2) are significantly similar, and all graphs of RMSF characteristics show less fluctuation, which indicates excellent stability of the two complexes. In addition, the appearance of the same amino acids in molecular docking and dynamic simulations as Asn167, Tyr202, Leu242, Asn258, Ala317, Lys352, Thr353, and Ala354 complex (Tubulin–PredM1) and Gly237, Val238, Leu255, Asn258, and Ile318 for (Tubulin–PredM2). The molecular simulation dynamics results indicate that all docking molecular results are validated.

Figure 8

Flexibilities of the protein backbones examined by RMSF values.

The total number of hydrogen bonds and average center-of-mass distance formed between the ligand (M1, M2) and protein (Tubulin) are shown in Figure 9. For the (Tubulin–PredM1) complex, the maximum number of hydrogen bonds reaches three at 25 and 85 ns, respectively. After that, this number reaches two hydrogen bonds at 6, 8, 18, 22–24, 62, 71–73, 78, and 98 ns. In the rest of the simulation, this number reaches 1. Generally, ligand M1 has H-bond gaps throughout most of the simulation time, even though the ligand is stable in the binding site. For the (Tubulin–PredM2) complex, the maximum number of hydrogen bonds reaches four at 16, 18, 58, and 98 ns. After that, this number reaches three, two, and one hydrogen bonds for the rest of the simulation. Ligand M2 forms a stable hydrogen bond network with the protein (Tubulin) over the simulation time.

Figure 9

(a) and (b) Hydrogen bonds of complex (Tubulin–PredM1) and (Tubulin–PredM2), respectively.

R _g is linked to the conformational state to show if a structure is folded and compact and has a stable conformation. The smaller the R _g value, the less flexible the proteins and the more stable the complexes. According to Figure 10, low R _g values for PredM1 (A) and PredM2 (B) led us to believe they were stable. R _g in the simulation was used to determine the compactness of the structure. The consistency throughout the simulation and the lower degree of fluctuation during the simulation indicate that the system is more rigid and compact. Both complexes (Tubulin–PredM1) and (Tubulin–PredM2) have the lowest R _g values, 21.13 and 20.97 Å, respectively. The R _g values varied between 21.13 and 21.43 Å for (Tubulin–PredM1) complexes and between 20.97 and 21.35 Å for (Tubulin–PredM2) complexes. In conclusion, all complexes showed relatively similar consistent R _g and compactness values, indicating they have good stability.

Figure 10

(Tubulin–PredM1) (a) and (Tubulin–PredM1) (b) R _g of the complexes during 100 ns MD simulation.

The magnitude of conformational changes that the water solvent can access was predicted using the graphs of SASA of the complexes (Tubulin–PredM1; Tubulin–PredM2). Hence, the interaction between these complexes and the solvent during the 100 ns MD simulation was analyzed using the SASA, as shown in Figure 11a. Figure 11a shows an SASA value varying between 1,825 and 1,875 Ų for Tubulin–PredM1 in the first 50 ns of simulations. However, after 50 ns, the solvent-accessible surface area varied between 1,850 and 1,925 Ų, indicating that the surface reached in the second period and remained stable. On the other hand, Figure 11b revealed values between 1,800 and 1,850 Ų in the first 38 ns. After the simulation interval [38–65 ns], the SASA value increases, reaching 1,825 and 1,875 Ų. In the rest of the simulation [66–100 ns], this value slowly decreases and varies between 1,800 and 1,850 Ų. The solvent-accessible zone reached by the [Tubulin–PredM1] complex is substantially similar to the zone reached by the [Tubulin–PredM1] complex; this value is higher and is very important, and SASA values for these two complexes (Tubulin–PredM1; Tubulin–PredM2) remain relatively stable, which indicates no significant changes during a 100 ns simulation.

Figure 11

SASA for 100 ns of the simulations. (a) and (b) represent Tubulin–PredM1 and Tubulin–PredM2, respectively.

MM/PBSA determines how stable newly designed compounds are for binding to tubulin proteins, using the fastest force field to calculate the free energy. The different types of energies calculated by the MM-PBSA method during 100 ns of simulation are listed in Table 10. The average binding free energies of (Tubulin–PredM1) and (Tubulin–PredM2) complexes calculated were found to be −101.775 ± 26.749 and −108.355 ± 77.540 kJ/mol, respectively. The average ΔE _MM for (Tubulin–PredM1) and (Tubulin–PredM2) complexes were found to be −181.36 ± 79.654 and −238.691 ± 62.879 kJ mol⁻¹, respectively. The average values of ΔG _Polar for the (Tubulin–PredM1) and (Tubulin–PredM2) complexes were −18.464 ± 3.421 kJ/mol and −22.222 ± 3.487 kJ/mol, respectively. The calculated average ΔG _SASA for (Tubulin–PredM1) and (Tubulin–PredM2) complexes were −18.464 ± 3.421 and 22.222 ± 3.487 kJ/mol, respectively. From the MM-PBSA result, we notice that all complexes showed almost closer values; this indicates good and excellent stability of our complexes.

Based on the retrosynthesis methods, the Spaya platform (https://www.spaya.ai) facilitates the synthesis of drug candidates with a high total score. Each route is a collection of similar synthetic steps documented in the literature (Scheme 1 and Table 11).

Scheme 1

The different stages of the prepared compound M1.

Table 11

Figure 12 shows the standardized prediction errors as a function of the lever values (hi) of the molecules in training, test sets, and the proposed M1 and M2 molecules.

Figure 12

Williams plot to evaluate the applicability domain of the model (h* = 0.53 and residual limits = ±2).

All compounds have levers hi < h ^∗ = 0.53, and all the tailings including M1 and M2 are found to be in the range of (2SD) (horizontal lines); this explains is a reliable model.