Virtual screening (VS)

The virtual screening procedure was described in great details in methods and Fig. 3. In summary, the ligand-based virtual screening relies on the structure of the reference ligand of Eg5 protein (PDB: 3ZCW). A library of drug-like ZINC database was screened to obtain lead compounds. Subsequently, 600 compounds were obtained, then the physicochemical properties and molecular docking studies were carried out and finally 80 compounds were selected. Further, visual inspections were applied and based on the binding mode within the target pocket, the binding score, the H-bond, the hydrophobic and aromatic interactions of the potential compounds with an allosteric pocket of Eg5 protein. Six compounds, named ZINC43068613 (compound 1), ZINC1300905 (compound 2), ZINC587060476 (compound 3), ZINC2639178 (compound 4), ZINC06444857 (compound 5) and ZINC43011180 (compound 6), were selected for further analysis (Fig. 4, Table 1). Interestingly, compound 1 and 6 were identified by Lahue [21], as Eg5 inhibitors that target α2/L5/α3 pocket, however in our study both compounds were identified as α4/α6/L11 inhibitors.

Molecular docking

The re-docked reference ligand showed a similar binding mode to co-crystalized ligand (Fig. 5), which suggests the docking protocol is efficient and solid. The docking results are summarized in Table 2, and Fig. 6. We can conclude that the shortlisted compounds are mimicking the same binding mode as the reference ligand of Eg5 protein. All the ligands occupied the back hydrophobic sub-pocket, that formed of Ile288. In addition, all molecules showed hydrophobic interactions with the key residues Leu292, Leu293 and Ile299. Furthermore, the key residues Tyr104 and Tyr352 formed pi-pi interactions with the six selected hits. All the six compounds formed H-bonds with the allosteric pocket residues. In further detail, compound 1 and compound 2 formed hydrogen bonds with the key residue Arg355. The C=O of compound 1 formed H-bond with NH-Arg355, while the C=O group of compound 2 formed two H-bonds with OH-Thr300 and NH-Arg355 at distances of 2.4 Å and 2.3 Å respectively. An additional H-bond was observed between compound 2 and Tyr352. Both compounds were embedded in the front and back hydrophobic pockets (Fig. 6). Interestingly compound 1 has a similar structure to the reference ligand, while the fluoro atom in position 3 of methoxybenzyl-amine group of the reference compound showed closer contact to Leu292 and Leu293 than compound 1. The methoxy group of compound 3 and compound 4 formed H-bonds with the allosteric pocket. Methoxy group of compound 3 formed an H-bond with NH-Asn271 at distance 2.9 Å, while the methoxy group of compound 4 formed an H-bond with OH-Tyr104 at a distance of 2.9 Å. Compound 4 formed another H-bond with NH-Arg355 at a distance of 3.4 Å. Both compounds 5 and 6 formed two H-bonds, compound 5 formed an H-bond between N-imidazole and NH-Asn289 at a distance of 2.5 Å and the second H-bond formed between N of purine and OH-Thr300 at a distance of 2.6 Å. While compound 6 formed an H-bond between the oxygen atom of SO2 and NH-Asn289 at a distance of 2.4 Å and the second H-bond formed between OH and NH-Arg297 at a distance of 2.5 Å. Both compounds showed hydrophobic interactions with the residues Leu292, Leu293 and Leu299.

Fig. 5

Superimposition of the original co-crystalized ligand benzimidazole (pink) and re-docked pose (yellow) in Eg5 allosteric site

Table 2 Docking results and interactions of hits compoundFig. 6

Binding mode A compound 1, B compound 2, C compound 3. Right Panel represented the 3D interactions of compounds and the pocket’s residues (purple, stick). H-bond represented as black dotted-line.Left panel represented the selected compounds (brown, stick) compared to the reference compound (magenta, stick) within the binding surface (Orange; hydrophobic, white; neutral; blue; hydrophilic). Binding mode D compound 4, E compound 5, F compound 6. Right Panel represented the 3D interactions of compounds and the pocket’s residues (purple, stick). H-bond represented as black dotted-line. Left panel represented the selected compounds (brown, stick) compared to the reference compound (magenta, stick) within the binding surface (Orange; hydrophobic, white; neutral; blue; hydrophilic)

Physiochemical and ADME properties

The predicted physicochemical properties of the selected hits, were calculated using DataWarrior and the Lipinski rule of five [57]. The co-crystalized ligand was used as a reference in these calculations, the results were summarised in Table 3. The values of physicochemical parameters of the 6 hits are H-bond donors ≤ 5, the H-bond acceptors ≤ 10, logPw/o ≤ 5, the rotatable bond ≤ 10, total polar surface area is ≤ 140 A [2], and the molecular weights ≤ 500 g/mol. The results indicate that the selected hits follow Lipinski rule, expect compound 6 as it showed a molecular weight of 503.49 g/mol which is slightly more than the Lipinski limit, but it is still in the accepted range. The log S of the 6 hits were calculated and compounds (1, 2, 4, 5, and 6) showed moderate solubility. However, compound 3 showed slightly poor solubility, but so did the reference compound. These results revealed that 5 hits are exhibiting drug-like properties and are absorbed well by the biological system. The pharmacokinetic (absorption, distribution, metabolism, excretion) properties of the six hits and reference compound were investigated (Table 4). The first five compounds showed high absorption from the intestine similar to the reference compound, while compound 6 showed low absorption. All the hits do not penetrate the blood brain barrier. Cytochrome P450 is one of the important detoxification enzymes in humans; many drugs are activated or inhibited by cytochrome P450. It was found that compounds 1 and 3 are predicted to be deactivated by CYP2D6 similar to the reference compound. Compounds 2, 4, 5, and 6 are predicted to be non-inhibitors for CYP2D6. Regarding excretion compounds 1, 3, 4, and 6 are not substrates to renal OCT2 enzyme, like the reference compound, while compound 2 and 5 are substrates to the renal OCT2 enzyme.

Table 3 Predicted the physicochemical properties of the selected compoundsTable 4 ADME (Absorption, Distribution, Metabolism, and Excretion) properties of selected compoundsToxicity studies

The toxicity prediction study of the hits and reference compound are summarised in Table 5, the hits were expected to have a non-toxicity profile. The hits and the reference compound are predicted to be non–mutagenic. In addition, the hits were predicted to be non-carcinogenic based on the FDA rodent carcinogenicity model. The oral rate acute toxicity LD50 values of all the compounds, except compound 4, are around 2.42–2.78 mol/kg. This is similar to the reference compound’s rate of 2.50 mol/kg, while the LD50 of compound 4 is 3.09 mol/kg. The six compounds showed rat chronic lowest observed adverse effect level (LOAEL) values that ranged from − 0.02 to 1.51 g/kg, which is extremely close to the reference compound 1.10 g/kg. Interestingly, the hits are non-irritant to the skin.

Table 5 Predicted toxicity of selected hitsBiological studies of selected hitsMT ATPase assay

The Eg5 ATPase activity of the six hits at 10 µM was investigated in comparison to Monastrol (10 µM), as it acted as a positive control compound. The results are shown in Fig. 7A and indicate that compound 5 showed significant anti MT ATPase activity. Compound 5 showed 50% inhibition activity, in comparison to Monastrol with 60% inhibition activity. Similarly compounds 1, 3 and 4 caused 40, 35 and 30% inhibition activity respectively. Compounds 2 and 6 showed the lowest inhibitory activity at around 25%. The results revealed that compound 5 showed the most promising inhibition activity, which prompted the measurement of its IC50 against Eg5 ATPase; it was found that the IC50 of compound 5 is 2.37 ± 0.15 µM. Although some compounds showed better IC50 than compound 5, however cellular resistance was developed for these compounds, which makes compound 5 a promising lead compound with new scaffolding.

Fig. 7

A Represented enzymatic Eg5 –ATPase inhibition activity of the 6 hits. Each result is a mean of 3 replicate samples and values are represented as % inhibition (± standard deviation). Immunofluorescence assay of monastrol and compound 5, the Hela cells were treated for 24 h with B DMSO as a negative control, C Monastrol as a positive control, D compound 5 at 10 µM then fixed and stained with anti-α-tubulin antibody (green) and with DAPI for DNA (blue) to visualize the microtubules

Immunofluorescence assay and inhibition of mitotic spindle formation

Further investigation of compound 5 was performed, which consisted of an immunofluorescence assay that was used to analyse the mechanism of compound 5 on the organisation of the microtubules into mitotic spindles during divisions of the cells. The cells were treated with the negative control (10 µM), which showed typical bipolar spindles (Fig. 7B), while Monastrol and compound 5 (10 µM) treated cells formed monopolar spindle profiles (Fig. 7C, D). These results revealed that compound 5 caused tubulin assembly distortion with irregular morphology, showing a typical mitotic arrest similar to Monastrol.

Molecular dynamic simulations

Molecular dynamics (MD) provide deep insights into ligand-receptor interactions and a wide understanding of the ligand-receptor conformational change over the time. MD simulations consider the flexibility of the protein, which is not taken into account by molecular docking programs, which results in more reliable results. In this study molecular dynamics were performed 3 times for 200 ns with different initial velocities. The top docked poses of compound 5 and reference compound complexes with the allosteric site of Eg5 protein were selected, and the Apo-protein for molecular dynamic simulations using GROMACS software.

Root mean square deviations (RMSDs) and root mean square fluctuations (RMSF)

In order to investigate the stability of compound 5—Eg5 complex and analyse any conformational changes, the Cα root-mean-square deviations (RMSDs) of inhibitor complex was calculated and compared with the reference ligand complex and Apo-protein. The average values of RMSD of the complexes are represented in Table S1. The average RMSD of the Apo-protein, compound 5-complex and reference compound—complex showed values 3.5 ± 0.2 Å, 3.0 ± 0.3 Å and 3.1 ± 0.3 Å. It could be observed from Fig. 12A, B that the compound 5 complex reached equilibrium at 20 ns, with an average RMSD value around 3 Å, the complex continued smoothly until the end of the simulation. Compound 5 remained in a similar position to the original pose in runs 1, 2 and 3. Interestingly, RMSD of compound 5-complex is similar to RMSD of Apo-protein during most of the simulations, the RMSD plot of compound 5 complex and Apo-protein overlap in most of the simulations (Fig. 8A). The RMSD of the reference compound complex reached equilibrium at 20 ns with a value around 3 Å. RMSD plot continued smoothly and no significant changes occurred during the remaining time period of the trajectory (Fig. 8B). These results indicate the stability of the compound 5 complex and its similarity to the reference ligand complex and Apo-protein.

Fig. 8

A RMSD of Cα of compound 5-complex (red) and Cα of Apo-protein (black). B RMSD of Cα of native ligand complex (red) and Apo-protein (black). C RMSD of compound 5 (red), native ligand (yellow). D RMSF of compound 5-complex (red), native ligand –complex (yellow) and Apo-protein (black)

In order to understand the binding mode of compound 5 and its stability within the allosteric pocket, the RMSD of the compound 5 alone was calculated relative to the starting docked pose, and compared with the RMSD of the reference compound (Fig. 8C). It could be seen that RMSD of compound 5 reached equilibrium at 20 ns and the plot continued in a straightforward manner until the end. The average RMSD value was 2 Å. On the other hand, the reference compound reached the equilibrium at 10 ns with RMSD value 1.5 Å, at 70 ns the RMSD decreased sharply to around 0.5 Å and continued like that until the end of the 200 ns. These results reveal compound 5 adopted one stable conformation over the dynamics. We investigated the conformations of the reference ligand at RMSD 1.5 Å and 0.5 Å. It was found that the fluoro-methoxy phenyl group adopted two conformations at 1.5 Å while at 0.5 Å; the fluoro-methxy group adopted the same conformation (Figure S1). The overall results are that compound 5-complex is stable over the dynamics and adopts one conformation during the entire trajectory.

Root-mean-square fluctuations (RMSF) of Cα atoms allow direct insight into the structural fluctuations and the degree of the residues flexibility in the protein. RMSF values of Cα atoms of compound 5-complex, reference compound-complex and Apo-protein were computed (Fig. 8D). The average RMSF (Table S1) of the Apo-protein, compound 5-complex and reference compound were 0.33 ± 0.4 Å, 0.28 ± 0.30 Å and 0.30 ± 0.3 Å. A close observation of RMSF plot reveals that the residues of the three systems displayed a low degree of flexibility with an average value 0.3 Å. The highest two peaks with RMSF were around 0.8 Å, the first one is between Asp43-Glu48 which is a loop, and the second one is between residues Val264-Asn271, which represents another loop. The three systems show nearly the same RMSF profile with low fluctuations in allosteric site α4 (292–297) and α6 (348–355). Overall, the three systems show lower residue fluctuations and lower RMSF values, than the Apo-protein residues. Most of the fluctuations are in the loops. These small ranges of RMSFs demonstrate that compound 5 is capable of forming suitable and stable interactions with the allosteric pocket during the dynamics. These results are in accordance with the findings from the RMSD results.

Solvent accessible surface area (SASA) and radius of gyration (Rg) analyses

SASA is the surface area accessible to the water molecule; SASA calculations can help to investigate the conformational dynamics of the compound 5—Eg5 protein. The average SASA values (Fig. 9A) for compound 5—protein, reference compound—complex and the Apo—proteins are 18,200Å2, 18,200Å2 and 18,300Å2. It could be seen that the SASA plots for compound 5 and reference compound complexes have overlapped in all the simulations. The SASA for both complexes started high with value 19,500, then reduced gradually until 50 ns to be 18,200 Å2 and continued smoothly until the end of trajectories. The SASA of Apo-protein started with 18,500 Å2 then reduced to 18,300 Å2, suddenly at 80 ns the SASA value of Apo-protein increased again to 19,000 Å2 until 110 ns and decreased again to 18,300Å2. The SASA results reveal the stability of compound 5—complex is similar to the reference compound, and that there were less inner residues interacting with the surroundings. The SASA fluctuations of the Apo-protein are slightly higher than the compound 5 complex which indicates that there are more inner residues in contact with the solvent, and binding compound 5 to the allosteric site reduced this contact and stabilized the Eg5 protein.

Fig. 9

A SASA of compound 5-complex (red), native ligand –complex (orange) and Apo protein (black). B Rg of compound 5-complex (red), native ligand –complex (orange) and Apo protein (black)

Radius of gyration (Rg) is another parameter which is linked with the tertiary structure and the general conformational state defining our understanding of compactness and the folding of proteins. The Rg of the three complexes were calculated and the results are represented in Fig. 9B. The three systems showed the same Rg average 20.09 Å. A fluctuation of 0.1 Å, around the average values were observed for compound 5 and reference complexes (Fig. 9B) at 70 ns and immediately after 20 ns the graph returned to its average value. These sudden deviations may be attributed to protein’s packing. All three graphs attained equilibrium around their average value, thus suggesting that the complex relaxed throughout the dynamics and confirms that compound 5-complex remained compact, and the folding of the Eg5 protein was maintained.

H-bonds and ligand–protein contact

The MD trajectory of compound 5 was analysed to gain more information about the stability of the ligand within the binding pocket and intermolecular interactions (H-bonds, hydrophobic, pi-pi interactions) over the simulation. It was found that compound 5 occupied the same position within the binding pocket from the starting point till the end (Fig. 10A). The distance from the centre of mass of compound 5 to the centre of mass of the allosteric pocket was around 2 Å (Fig. 10B) and the plot was smooth and stable over the 200 ns, these results confirm the stability of the inhibitor within the pocket. Investigations about the intermolecular interactions of compound 5 over the dynamics showed that compound 5 formed four H-bonds with the binding site (Fig. 10C) with the interactions percentage being 97%. The most stable and continuous H-bond was with Thr300 which persists over the simulation with occupancy of 95%. Furthermore, at 20 ns two H-bonds formed between compound 5 and residues Gly296 and Arg355 with percentages of 70% and 50% respectively, but when Arg355 or Gly296 moved away from compound 5, these two H-bonds were broken. At the same time, another two H-bonds formed between compound 5 and residue Asn289 and Tyr352 (Fig. 10D). It could be observed compound 5 was able to maintain the initial binding mode of the simulation through forming strong polar interactions with residues Asn271, Asp279, Asn287 and Arg297. Besides the pi-pi interactions with Tyr104 and Tyr352 (Fig. 10D). Forming continuous H-bonds with Asn289, Gly296, Tyr352 and Arg355 with the average percentage more than 50%, this suggests that the H-bonds may be the main reason for the stabilization of compound 5. In particular, the stable H-bond formed with Thr300 was almost preserved in throughout the entire dynamics. Additional hydrophobic interactions were established with Ile288, Leu292, Leu293, and Ile332, as they had high occupancy which was greater than 60% (Fig. 10D) so fortified the stability of the ligand over the whole trajectory.

Fig. 10

A Snapshot of binding mode of compound 5 at the initial of the simulation (green) and at the end of 200 ns. B The Distance from the centre of mass of Eg5 allosteric site to the centre of mass of the compound 5. C No of H-bonds of compound 5 with the allosteric pocket’s residues. D Protein–ligand contact mapping for Eg5 protein with compound 5

Principle component analysis (PCA)

The collective motion of compound 5—Eg5 protein, reference compound—Eg5 protein and Apo-protein was computed from the trajectories using PCA method. The method based on the constructions of the diagonal covariance matrix from Cα atoms of the Eg5 protein, which captures the global motion of the atoms through eigenvectors or the principal components (PCs) and eigenvalues. The eigenvectors explain the global direction of motion of the atoms, while the eigenvalues represent the atomic contribution of motion in MD trajectories of each system. The reduction in the subspace size was determined using the scree plot, the distribution of the eigenvectors versus eigenvalues. Figure 11A demonstrates a sharp fall in the slope at the fifth PC. The first eigenvector accounted for 78.9% of the overall variance, the first three eigenvectors together accounted for roughly 92% of the total variance. The first three eigenvectors were selected to calculate the reduced subspace.

Fig. 11

A The changes in the eigenvalues with increasing the eigenvectors. B The projection of each trajectory on the first two eigenvectors. C The projection of each trajectory on the first and third eigenvectors. D The projection of each trajectory on the second and third eigenvectors

Bi-dimensional projection studies

The projecting of each trajectory of the three systems was represented in Figs. 11B–D. The conformational subspace of the systems was evaluated using the first three eigenvectors of the total Cα. The PCA study can help in understanding the dynamic behaviour of compound 5—complex and compare that with dynamics behaviour of the reference complex and Apo-protein. The PCA of the first and second eigenvectors revealed that the conformational clusters of compound 5-complex are well defined and covered a minimal number of subspaces (Figs. 11B) in comparison to the Apo-protein which occupied large subspaces. Moreover, the trajectories of compound 5 and the reference compound overlapped in most of the dynamics with different average structures. Similarly, the PC between the eigenvectors one and three (Fig. 11C) and between eigenvectors two and three (Fig. 11D) showed nearly the same atomic motions. These results indicate that compound 5 and reference compound complexes cover the minimum number of subspace, while the Apo-protein showed large atomic motions and conformational changes. The trajectories of compound 5—Eg5 and reference compound—Eg5 overlapped in most of the PCA analyses and occupied a small subspace range. These finding indicate that compound 5 showed restricted subspaces in a complex with Eg5 protein, leading to a well-defined internal motion behaviour, which stabilized the complex better than Apo-protein.

MM/GBSA binding free energy

The MM/GBSA method was used to calculate the free binding energy of compound 5 and the reference compound. The components of the binding free energy were represented in Table 6. Compound 5 has a binding energy slightly higher than the reference. The free binding energy for compound 5—complex was (− 46.2 ± 0.2 kcal/mol) compared to (− 44.1 ± 0.2 kcal/mol) for the reference complex). The most favourable interaction was Van der Waals interaction with an average around − 63.1 ± 0.22 kcal/mol for compound 5—Eg5 and an average around − 57.8 ± 0.12 kcal/mol for the reference complex. The electrostatic interaction for compound 5—complex was a favourable interaction with a value of − 21.0 ± 0.12 kcal/mol, which was positive value for the reference compound 12.4 ± 0.1 kcal/mol. Both systems showed the same non-polar interactions of − 6.8 ± 0.1 kcal/mol (Table 6).

Table 6 Binding free energies and its components, compound 5, reference compound complexes

Further, per-residue decomposition analysis of compound-5 complex was computed to determine the contributions of residues to the binding free energy (Fig. 12). It was found that Tyr104 and Tyr352 contributed with an average energy of around − 5 kcal/mol and − 3 kcal /mol and this can be explained by the face-face aromatic interactions with compound 5, and H-bond with Tyr352. Furthermore, Asn289, Gly296 and Thr300 contributed well to the binding affinity; the three residues formed H-bonds in high incidences, more than 50%. The hydrophobic residues Leu292, Leu293 and Leu295 formed strong hydrophobic interactions with compound 5 with proportion more than 90%. In addition, the two residues Asn271 and Arg277 formed polar interactions with compound 5, around 60% occurrence and contributed significantly to the binding free energy. It could be concluded that the key residues Tyr104, Asn289, Leu292, Leu293, Gly296, Arg297, Thr300 and Arg355 which formed the allosteric site (α4/α6/L11) contributed to the binding free energy of compound 5 with an of average − 5, − 3, − 4, − 5, − 4, − 3, − 4 and − 6 kcal/mol. Interestingly compound 5 showed a slightly more favourable binding energy than that for the reference compound, which may be attributed to the continuous H-bonds with the key residues Asn289, Gly296, Tyr352 and Arg355 and the pi-pi interactions with Tyr104, which resulted in significant contribution ≥ − 3 kcal/mol to the binding free energy of compound 5.

Fig. 12

MM-GBSA free energy decomposition of EG5 residues