Receptor-based pharmacophore modeling, molecular docking, synthesis and biological evaluation of novel VEGFR-2, FGFR-1, and BRAF multi-kinase inhibitors

A receptor-based pharmacophore model describing the binding features required for the multi-kinase inhibition of the target kinases (VEGFR-2, FGFR-1, and BRAF) were constructed and validated. It showed a good overall quality in discriminating between the active and the inactive in a compiled test set compounds with F1 score of 0.502 and Mathew’s correlation coefficient of 0.513. It described the ligand binding to the hinge region Cys or Ala, the glutamate residue of the Glu-Lys αC helix conserved pair, the DFG motif Asp at the activation loop, and the allosteric back pocket next to the ATP binding site. Moreover, excluded volumes were used to define the steric extent of the binding sites. The application of the developed pharmacophore model in virtual screening of an in-house scaffold dataset resulted in the identification of a benzimidazole-based scaffold as a promising hit within the dataset. Compounds 8a-u were designed through structural optimization of the hit benzimidazole-based scaffold through (un)substituted aryl substitution on 2 and 5 positions of the benzimidazole ring. Molecular docking simulations and ADME properties predictions confirmed the promising characteristics of the designed compounds in terms of binding affinity and pharmacokinetic properties, respectively. The designed compounds 8a-u were synthesized, and they demonstrated moderate to potent VEGFR-2 inhibitory activity at 10 µM. Compound 8u exhibited a potent inhibitory activity against the target kinases (VEGFR-2, FGFR-1, and BRAF) with IC50 values of 0.93, 3.74, 0.25 µM, respectively. The benzimidazole derivatives 8a-u were all selected by the NCI (USA) to conduct their anti-proliferation screening. Compounds 8a and 8d resulted in a potent mean growth inhibition % (GI%) of 97.73% and 92.51%, respectively. Whereas compounds 8h, 8j, 8k, 8o, 8q, 8r, and 8u showed a mean GI% > 100% (lethal effect). The most potent compounds on the NCI panel of 60 different cancer cell lines were progressed further to NCI five-dose testing. The benzimidazole derivatives 8a, 8d, 8h, 8j, 8k, 8o, 8q, 8r and 8u exhibited potent anticancer activity on the tested cell lines reaching sub-micromolar range. Moreover, 8u was found to induce cell cycle arrest of MCF-7 cell line at the G2/M phase and accumulating cells at the sub-G1 phase as a result of cell apoptosis. Supplementary Information The online version contains supplementary material available at 10.1186/s13065-024-01135-0.


Introduction
Protein kinases are a class of phosphotransferases that play a fundamental role in the regulation of different cellular processes such as cellular survival, growth, proliferation, migration, and apoptosis [1].More than 30% of cellular proteins are phosphorylated by protein kinases.Protein kinases catalyse the transfer of a gamma phosphate group from ATP to an acceptor amino acid (Tyrosine, serine, or threonine) in a substrate protein [1].Therefore, protein kinases are categorized mainly into two main categories, the tyrosine kinases such as vascular endothelial growth factor receptor (VEGFR), fibroblast growth factor receptor (FGFR) and epidermal growth factor receptor (EGFR) and the serine/threonine kinases such as rapidly accelerated fibrosarcoma (RAF) kinases [2,3].In cancer, several protein kinases are dysregulated resulting in the uncontrolled growth, survival, and metastasis of tumour cells [2,4,5].Hence, targeting protein kinases has received a remarkable attention in recent years for the discovery of new targeted chemotherapeutic agents for cancer treatment [6][7][8].Based on the fact that cancer is regulated by multiple pathways that can compensate for one another when a single pathway is blocked, targeting multiple kinases is a more efficient strategy than targeting a single kinase [9].Moreover, multi-kinase inhibition has numerous advantages, such as increasing potency due to its synergistic effect, reducing probable polypharmacy toxicity, avoiding pharmacokinetics incompatibilities, and enhancing selectivity [9,10].
Angiogenesis, the formation of new blood vessels, plays an essential role in the growth and metastasis of tumour cells [11].Hence, targeting protein kinases that initiate and sustain the angiogenic process is a prominent approach in cancer treatment [12].Vascular endothelial growth factor (VEGF) and its receptors (VEGFRs) is a tyrosine kinase system that plays a curial role in angiogenesis both in the physiological as well as pathological conditions [13][14][15].In comparison to healthy tissues, VEGFR-2, in particular, is overexpressed in various types of cancer such as malignant melanoma, breast cancer, hepatocellular carcinoma, colon cancer, etc., [16,17].In addition, hFGFR family is a group of four isoforms; FGFR-1 to FGFR-4 that are expressed on the cell membrane and participate in various vital physiological and pathological processes, such as proliferation, differentiation, cell migration, survival, as well as angiogenesis [18].Binding of FGFR to its growth factor (FGF) results in its dimerization and phosphorylation of its intracellular kinase domain resulting in the initiation of a series of downstream signalling pathways.FGFRs overexpression has been reported in different types of solid tumours, for instance, FGFR-1 is amplified in breast and non-small cell lung cancers [19].Hence, it is believed that FGFRs inhibition by small molecules that competitively bind to the ATP binding pocket is an attractive tactic for the design of novel targeted anticancer agents [20].
Moreover, when the main pro-angiogenic factors VEGF and FGF bind to their target receptors, they result in an activation of the mitogen-activated protein kinase (MAPK) signalling pathway [21].Downstream signalling of this pathway leads to activation of RAS proteins which in turn causes subsequent activation of RAF kinases [22,23].RAF kinases are an intracellular serine/threonine family mediating several transcriptional factors leading to cell growth, survival, and proliferation [22,23].Among the RAF family, BRAF is the most sensitive isoform to activation and mutation [24].
The X-ray crystallographic structures of kinases such as VEGFR-2, FGFR-1 and BRAF demonstrated that their kinase domain comprises a smaller N-terminal lobe, larger C-terminal lobe, and an in-between ATP binding region which can be partitioned further into front pocket (front cleft or hinge region), gate area, and back cleft (allosteric back pocket).At the beginning of the C-terminal lobe there is an activation loop (A-loop) that is characterized by a highly conserved aspartatephenylalanine-glycine (DFG) motif.Based on the 3D orientation of the DFG motif, the A-loop can exist in different conformations resulting in the existence of the protein kinase in its active (DFG-in) or inactive (DFGout) conformations.In the catalytic cycle, the protein kinase switches between both open and closed conformations [25,26].Analysis of the binding modes of the co-crystallized protein kinase inhibitors (PKIs) at their target proteins demonstrated that PKIs can be classified according to their binding modes into six different types [26][27][28][29][30].Among them, type II inhibitors are regarded as promising ones performing their antagonistic activity on the inactive (DFG-out) conformation accommodating into the hinge region, the gate area and extend further to the less conservative back pocket enhancing their affinity, selectivity, and residence time [31].For example, sorafenib (I) (Fig. 1) is a VEGFR-2 (PDB ID: 4ASD) and BRAF (PDB ID: 1UWH) inhibitor in which the picolinamide moiety (coloured red in Fig. 1) occupies the front pocket and performs hydrogen bonding interaction with Cys919 (VEGFR-2)/Cys531 (BRAF).The ureido moiety (coloured pink in Fig. 1) extends through the gate area and forms by its NH group a hydrogen bond with the carboxylate group of αC-helix Glu885 (VEGFR-2)/Glu593 (BRAF), furthermore, the oxygen atom of the ureido moiety forms a hydrogen bond with the N-H group of DFG's Asp1046 (VEGFR-2)/Asp500 (BRAF).In addition, sorafenib (I) extends into the hydrophobic back pocket by a disubstituted phenyl group (coloured blue in Fig. 1) achieving multiple hydrophobic interactions with the surrounding residues [32,33].
In recent years, significant advancements have been achieved in PKIs exploration, mostly attributed to the utilization of computational techniques [44,45].These methods have proven instrumental in delivering valuable insights into diverse protein kinase structures and inhibitors [44,45].In computer-aided drug design (CADD), two primary strategies are commonly employed; structure-based drug design (SBDD) and ligand-based drug design (LBDD).These approaches allow researchers to predict and optimize the properties and activities of molecules even before they are synthesized and tested in the laboratory [46].
Due to the continuous resistance development by cancer cells on one hand and the more satisfactory effect and less drawbacks achieved by the concurrent targeting of multiple protein kinases on the other hand [47,48], there is a continuous demand for developing small molecules that target more than one protein kinase simultaneously (multi-kinase inhibitors).
Encouraged by these facts, the ultimate goal of the current investigation is to extract the common pharmacophoric features required for achieving multi-kinase inhibition of the target kinases; VEGFR-2, FGFR-1, and BRAF.The generated pharmacophore model will be then used to virtually screen a set of in-house synthetically feasible tailored diverse scaffolds to select those satisfy the pharmacophoric features of the target multi-kinase activity.Scaffolds satisfy the constructed pharmacophore model will be structurally optimized to enhance their target kinase binding.Molecular docking will be then used to confirm the ability of the designed derivatives to perform the essential interactions with the three target kinases.Afterwards, the in silico promising derivatives will be synthesized and evaluated for their biochemical inhibitory activity against the target kinases as well as for their cytotoxic activity on several cancer cell lines.Finally, the most potent candidate will be further analyzed for its effect on cell cycle progression and apoptosis induction.

Molecular modeling study
For the intended study, a common 3D multi-kinase pharmacophore model for type II kinase inhibitors of the target kinases (VEGFR-2, FGFR-1, and BRAF) was constructed using the receptor-based pharmacophore technique.The generated pharmacophore was then validated for its ability to discriminate between active and inactive compounds of the different kinases of interest using a pre-compiled test set of active inhibitors and inactive decoys.Next, the validated pharmacophore was used to virtually screen several in-house datasets with diverse scaffolds and the most promising scaffold was used as a starting point to develop several optimized derivatives with potential synthetic feasibility which should be, by design, multi-kinase inhibitors of the target kinases.Finally, molecular docking simulations were used to investigate the ability of the designed compounds to bind to the active sites of the target kinases accomplishing the key interactions responsible for the kinase inhibitory activity.Promising compounds were passed to the next chemical synthesis step.

Common 3D multi-kinase receptor-based pharmacophore model generation
X-ray crystallographic structures Several X-ray crystal structures for the target kinases (VEGFR-2, FGFR-1, and BRAF) are available in the protein data bank [49].For the current work, we are focusing on the design of type II kinase inhibitors due to their reported superiority [20].Thus, representative structures were selected that are cocrystalized with potent structurally diverse type II kinase inhibitors which bind to the inactive DFG-out kinase conformation occupying the front cleft (hinge region), the gate area and extend beyond the gatekeeper into the hydrophobic allosteric back cleft [50].Different kinase structures co-crystalized with the same ligand were also preferred.Hence, the X-ray crystallographic structures of VEGFR-2 (PDB ID: 3VHE, 3VNT, and 3VO3), FGFR-1 (PDB ID: 4V01 and 3RHX), and BRAF (PDB ID: 4DBN and 6B8U) were downloaded from the protein data bank [49].The protein structures were then prepared and aligned using their alpha carbons (See Additional file 1: Section 1: computational studies for further details).
Manual 3D receptor-based pharmacophore models generation Using the aligned prepared protein structures, several manual 3D pharmacophores were generated to describe the common inhibitors' interactions.The main common ligand-target interactions involve H-bonding interaction with the hinge region Cys919, Ala564, and Cys532 in VEGFR-2, FGFR-1, and BRAF, respectively, H-bonding with DFG Asp1046, Asp641, and Asp594 in VEGFR-2, FGFR-1, and BRAF, respectively, and H-bonding with αC-helix Glu885, Glu531, and Glu501 in VEGFR-Fig.2 Benzimidazole-based multi-kinase inhibitors and their interactions in their targets' kinase domain 2, FGFR-1, and BRAF, respectively.These interactions were described by hydrogen bond acceptor, donor, and acceptor features, respectively, with their corresponding projected features (Site features).In addition to hydrophobic interactions with the hydrophobic allosteric back pocket in each protein structure which were described by a broader hydrophobic pharmacophoric feature.In addition, excluded volumes were employed to define the binding sites' steric extent.The different 3D pharmacophores obtained are the result of several combinations of the different pharmacophoric features (in terms of their number and volume) giving a set of 17 combinations (See Additional file 1: Section 1: computational studies for further details).

Pharmacophore model selection and validation
Selection of the best 3D pharmacophore model was carried out with the aid of a compiled test set of 2387 compounds (Table 1) (See Additional file 1: Section 1: computational studies for further details).
The 3D pharmacophore ability to discriminate between the test set active and inactive compounds was used to evaluate its quality which was assessed based on its collective results on the whole test set.For each 3D pharmacophore, the total number of true positives TP t , false positives FP t , true negatives TN t , and false negatives FN t were determined from its performance on each kinase test set (See Additional file 1, for further details) and were used in calculating the different assessment metrics; Sensitivity Se, specificity Sp, yield of actives Ya, enrichment E, accuracy Acc, discrimination ratio DR, F1 score F1 and Mathew's correlation coefficient MCC to evaluate the models' performance (Table 2 (Metric's values of the best performing pharmacophore model (Ph4-4) are shown in bold) and see Additional file 1: Section 1: computational studies for further details).
As can be seen in Table 2, Ph4-16 and Ph4-17 showed low sensitivity (0.479 and 0.425, respectively) meaning that they yielded a low number of true positives, however, they showed good specificity (0.952 and 0.986, respectively) and so could discard decoys and correctly consider them as inactive compounds, so these two models are biased towards inactive compounds and that is reflected in their low F1 score and MCC (0.3196 and 0.3100, respectively, for Ph4-16 and 0.4526 and 0.4375, respectively, for Ph4-17) (Table 2 and Fig. 3).On the contrary, models Ph4-12 to Ph4-15 showed good sensitivity (0.808, 0.849, 0.849 and 0.849, respectively) meaning that they yielded a high number of true positives, however, they showed low specificity (0.905, 0.773,  0.824 and 0.861, respectively) and so could not discard decoys properly and predict large number of decoys as active compounds, so these models are biased towards active compounds and this is reflected in their low MCC (0.3822, 0.2486, 0.2918 and 0.3328, respectively) (Table 2 and Fig. 3).Models Ph4-1 to Ph4-10 showed a balance between sensitivity and specificity with a sensitivity range of (0.726-0.767) and a specificity range of (0.919-0.964) indicating that the models are not biased towards either of actives or decoys (Table 2 and Fig. 3).Model Ph4-4 (Fig. 4) was selected as the best model because it showed the best overall performance on the test set.Ph4-4 selected 146 hits out of 2387 compounds of which 55 compounds were true positives (Se  2 and Fig. 3).
Figure 4 shows the selected 3D pharmacophore model, Ph4-4, its pharmacophoric features, and interfeature distances (in Å) in 3D space.Ph4-4 consists of 4 main features [F1-F4].Feature 1 (F1:Acc), a hydrogen bond acceptor, where ligands bind to the hinge region Cys or Ala, and the direction of this hydrogen bond acceptor lone pair is indicated by its projected site point feature (F5:Acc2).Feature 2 (F2:Don), a hydrogen bond donor, describing the feature required for binding to the glutamate residue of the Glu-Lys αC helix conserved pair and its projected site point feature (F6:Don2) indicates the direction of the hydrogen bond donor hydrogen.Feature 3 (F3:Acc), a hydrogen bond acceptor mapping where the ligands bind to the DFG motif Asp at the activation loop, in addition to its projected site point feature (F7:Acc2) indicating the direction of the hydrogen bond acceptor lone pair.Finally, the broadest feature (F4:Hyd) where ligands' hydrophobic moieties occupy the allosteric back pocket next to the ATP binding site.Moreover, thirty excluded volumes (Not shown in Fig. 4 for clarity) were also added to this pharmacophore to define the steric extent of the binding sites and to restrict the highly flexible compounds (if any) to bind in the desired conformations to the binding site, simulating the actual binding scenario.

Virtual screening
The pharmacophore model exhibited the best performance on the test set (Best discrimination between actives and inactive compounds), Ph4-4, was then used to screen an in-house dataset of diverse scaffolds.The benzimidazole scaffold 8a was selected by Ph4-4 as the most promising hit with the least RMSD from the assigned pharmacophore model features' centroids (RMSD = 0.979Å) (Fig. 5).

Hit optimization
The promising hit scaffold was then used as a starting point to develop several optimized derivatives 8a-u with potential chemical synthesis feasibility and probable good binding affinity which should be, by design, multi-kinase type II inhibitors for the target kinases (Fig. 6).
The design strategy took into consideration that the benzimidazole core would occupy the gate area of the target protein kinases and the imidazole moiety would be involved in hydrogen bonding with the key amino acids Glu885 and Asp1046 in VEGFR-2, Glu531 and Asp641 in FGFR-1, as well as Glu501 and Asp594 in BRAF.Thus, the 5-position of the benzimidazole moiety was substituted with different aryl groups which are decorated with hydroxy and methoxy groups at different positions to be involved in hydrogen bonding with the key amino acid Cys919 (VEGFR-2), Ala564 (FGFR-1), and Cys532 (BRAF) at the hinge region.Moreover, (un) substituted aryl groups were introduced at 2-position of the benzimidazole scaffold to occupy the allosteric hydrophobic back pocket to engage in hydrophobic interactions with the surrounding amino acid (Fig. 6).

Molecular docking simulation
Molecular docking is a well stablished technique for the investigation of the binding mode and binding affinity of drug-like molecules in their proposed biological targets [51][52][53][54][55].In the current study, to confirm and to study the binding characteristics of the designed compounds in the binding sites of the target kinases VEGFR-2, FGFR-1 and BRAF, molecular docking studies were performed using Molecular Operating Environment (MOE, 2022.02)software.The X-ray crystallographic structures of VEGFR-2 (PDB ID: 4ASD), FGFR-1 (PDB ID: 4V01) and BRAF (PDB ID: 5CT7) in their DGF-out inactive conformation were downloaded from the Protein Data Bank (PDB) [32,39,49,56].The downloaded protein structures are co-crystalized with a type II PK inhibitor, sorafenib (I), ponatinib, and RAF265, respectively.Molecular docking setup was initially validated by selfdocking of the co-crystalized ligands in the binding sites of their corresponding target kinases.These simulations successfully reproduced the binding pattern of the co-crystalized ligands in the target binding sites, VEGFR- The docked compounds showed promising binding patterns in VEGFR-2, FGFR-1 and BRAF interacting with the key amino acids in their binding pocket.The benzimidazole ring fits in the gate area stabilized via hydrogen bond interactions.By its imidazole ring, it interacts with the side chain carboxylate of Glu885, Glu531 and Glu501 of the αC helix in VEGFR-2, FGFR-1, and BRAF, respectively, and/or with Asp1046, Asp594,

ADME properties prediction
SwissADME online web tool [57][58][59] was used to predict the physicochemical and AMDE properties of the target compounds 8a-u.SwissADME showed that the newly synthesized compounds possess promising predicted physiochemical and pharmacokinetic properties.
All target compounds possess promising predicted physicochemical properties and moderate predicted aqueous solubility.Moreover, they complied with Lipinski's rule of 5 indicating that they are predicted to be orally bioavailable, and they possess a predicted SwissADME bioavailability score of 0.55 (See Additional file 1: Section 1: computational studies for further details).showing its interaction with the FGFR-1 active site Furthermore, as shown in SwissADME Boiled-Egg chart (Fig. 10), all target compounds showed high predicted GIT absorption with no predicted blood brain barrier (BBB) permeation and so devoid of CNS side effects.Moreover, Fig. 10 shows that all compounds are not p-glycoprotein (P-gp) substrates.
In summary, the designed benzimidazole derivatives 8a-u are predicted to be promising type II-like multikinase inhibitors in terms of binding affinity and pharmacokinetic properties and can be progressed further to chemical synthesis and biological evaluation.

Chemistry
As can be seen in Fig. 11, the target 2,5-disubstituted benzimidazole derivatives 8a-u were synthesized by the reaction of the benzaldehyde derivatives 1a-c with sodium metabisulfite to give the corresponding intermediates 2a-c which were condensed with 3,4-diaminobenzoic acid (3) to give the 2-aryl-benzimidzole-5-carboxylic acids 4a-c [34].Esterification of 4 was then carried out to afford the corresponding ethyl esters 5a-c [60].Hydrazinolysis of the benzimidazole esters 5a-c was carried out to yield the benzimidazole acid hydrazides 6a-c [61] which was followed by the reaction with different hydroxy and methoxybenzaldehyde derivatives 7a-g to afford the target benzimidazoles 8a-u in excellent yield.

Biochemical assay
VEGFR-2 inhibitory activity All the target 2,5-disubstituted benzimidazoles 8a-u were screened for their inhibitory activity on VEGFR-2 at 10 µM concentration and the % of inhibition was depicted in Table 3 using sorafenib (I) as a reference standard.

Multi-kinase inhibitory activity of 8u
In reference to the potent activity of 8u on VEGFR-2 (Table 3), it was fur-  4.

Antiproliferative activity
Antiproliferative activity on NCI cancer cell lines at 10 µM concentration All the target 2,5-disubstituted benzimidazoles 8a-u were screened by NCI (USA) for their ability to supress the growth of NCI 60 cancer cell lines at 10 micromolar concentration and the results were depicted in Table 5.

Antiproliferative activity on NCI cancer cell lines at five different concentrations
The 2,5-diaryl benzimidazole derivatives 8a, 8d, 8h, 8j, 8k, 8o, 8q, 8r, and 8u were selected by NCI to be further assayed for their growth inhibitory activity at five dose level and their GI 50 results were depicted in Table 6.The selected 2,5-diaryl benzimidazoles displayed potent inhibitory activity against the tested cell lines with GI 50 up to 0.12 µM.Close Fig. 11 Schematic pathway for the synthesis of the target 2,5-disubstituted benzimidazoles 8a-u examination showed that 8d displayed the most potent GI 50 against most of the cell lines with mean GI 50 of 2.62 µM.Meanwhile compounds 8a, 8h, 8j, 8k, 8o, 8q, 8r and 8u demonstrated potent inhibitory activity against most of the tested cell lines with mean GI 50 range of 2.98 to 7.98 µM.

Cell cycle analysis
Encouraged by the potent multi-kinase inhibitory activity of the 2,5-diaryl benzimidazole derivative 8u as well as its potent and broad spectrum of antiproliferative activity, it was selected as a representative to be examined for its influence on the cell cycle progression of MCF-7 cell line at its GI 50 concentration by flow cytometry analysis using propidium iodide (PI) stain.Comparison with breast cancer MCF-7 cells treated with DMSO as control was carried out, and the results were presented in Fig. 13 and Table 7.The % of MCF-7 cells that was accumulated in the G1 phase showed a decrease from 56.30% to 45.32% after treatment with 8u.On the other side, the accumulated % of cells in the G2/M phase increased from 20.60% to 25.57% (Table 7).These results emphasized that 8u arrests the MCF-7 cell line at the G2/M phase.Furthermore, an increase in the percent of cells accumulated in the sub-G1 phase, from 2.64% in the control to 4.40% in the treated cells was noticed as a result of cell apoptosis.

Apoptosis assay
The capability of 8u to enhance apoptosis of MCF-7 cell line at its GI 50 concentration was explored (Fig. 14).The presented results showed that the % of cells in the early apoptosis and late apoptosis phase increased from 0.25% and 2.36% to 1.46 and 4.54%, respectively, after treatment with 8u, which indicates that 8u induces cell apoptosis in MCF-7 cell line.

Conclusion
Several receptor-based pharmacophore models describing the binding features required for the multikinase inhibition of the target kinases (VEGFR-2, FGFR-1, and BRAF) were constructed based on the experimental binding mode and binding interactions of several inhibitors for these target kinases.Using a compiled test set of 73 active inhibitors for the target kinases as well as 2314 inactive decoys, (Ph4-4) was selected as the best model showing F1 score of 0.5023 and Mathew's correlation coefficient of 0.5131 indicating its good overall quality in discriminating between the active and the inactive compounds.Virtual screening of an in-house dataset of diverse scaffolds using the selected pharmacophore model yielded a benzimidazole-based scaffold as a promising hit among the dataset compounds     with RMSD of 0.979Å.Structural optimization of the hit benzimidazole-based scaffold through (un) substituted aryl substitution on 2 and 5 positions of the benzimidazole ring produced compounds 8a-u.Based on molecular docking simulations and ADME properties predictions, the optimization products were predicted to be promising type II-like multi-kinase inhibitors in terms of binding affinity and pharmacokinetic properties and can be progressed further to chemical synthesis and biological evaluation.

Molecular modeling study
The molecular modeling study was carried out using Molecular Operating Environment software (MOE 2022.02)according to the following steps:

Common 3D multi-kinase receptor-based pharmacophore model generation
Retrieving X-ray crystallographic structures The X-ray crystallographic structures of VEGFR-2 (PDB ID: 3VHE, 3VNT, and 3VO3), FGFR-1 (PDB ID: 4V01 and 3RHX), and BRAF (PDB ID: 4DBN and 6B8U) were downloaded from the protein data bank [49].MOE was used to prepare the retrieved protein structures (For further details see Additional file 1: Section 1: computational studies).Finally, correctness of ligands' structures and reported ligand interactions at the active site were further checked after the protonation step.The different prepared protein structures of VEGFR-2, FGFR-1 and BRAF were aligned and superposed using Align protocol in MOE using protein structures' αCs.Consequently, the co-crystalized  ligands were aligned in their bioactive conformation in the binding sites of the different protein structures.
Manual 3D receptor-based pharmacophore models generation Using the pharmacophore query editor in MOE, the aligned co-crystalized ligands were used to generate several manual 3D pharmacophores based on their com-  mon interactions with the target kinases' binding sites in the different protein structures.The assigned pharmacophoric features include recognition, shape, and projected site point features representing the main common features responsible for the inhibitors' binding to the kinase domain hotspots in the different proteins.Moreover, several excluded volumes (with different volumes and number) were included to define the steric extent of the binding sites (See Additional file 1: Section 1: computational studies for further details).

Pharmacophore model selection and validation
For pharmacophore model selection and validation, for each kinase of the target kinases, a test set was constructed from active inhibitors and inactive decoys with type II-kinase-inhibitor-like structures (Based on their visual inspection).The test set compounds for each kinase were compiled from the Directory of Useful Decoys-Enhanced (DUD-E) [63] and/or DEKOIS 2.0 [64] databases.A large decoys/actives ratio (≈30) was maintained to mimic the natural ratio in the chemical space between the active and inactive compounds.
Table 1 shows the distribution of the actives and decoys for each target kinase (See Additional file 1, Section 1: computational studies for the structure of the active inhibitors used in the test set for each target kinase).Conformational search was then carried out for the test set compounds using Stochastic search in MOE which generates conformations by randomly rotating all bonds (including ring bonds) and randomly inverting tetrahedral centres followed by an all-atom energy minimization.
The generated conformers were virtually screened using the different manually generated pharmacophore models to test their ability to discriminate between the active and inactive compounds in the compiled test sets.MOE pharmacophore search-algorithm begins with prefiltration of the conformers database by calculating the conformer similarity to the pharmacophore model with respect to feature type and distance; followed by a more computationally expensive alignment of the conformer atoms to the query feature points minimizing their deviation from each other using root mean square deviation (RMSD) as the fitness criteria for the alignment quality.
The 3D pharmacophore ability in discriminating between the test set active and inactive compounds was assessed based on its collective results on the whole test set.For each 3D pharmacophore, the total number of true positives TP t , false positives FP t , true negatives TN t , and false negatives FN t were determined (see Additional file 1, Section: computational studies).To assess the performance of the different generated pharmacophores, a set of assessment metrics were used to select and validate the best one.These metrics include sensitivity Se, specificity Sp, yield of actives Ya, enrichment E,  accuracy Acc, discrimination ratio DR, F1 score and Mathew's correlation coefficient MCC to evaluate the models' performance (For further details see Additional file 1, Sect.1: computational studies) [65].

Virtual screening
The pharmacophore model exhibited the best performance on the test set in discriminating between actives and inactive compounds was used to screen an in-house dataset of diverse scaffolds after performing conformational search on the in-house dataset using the same protocol used for the test set compounds (Vide supra).

Hit optimization
The promising scaffold was then used as a starting point to develop several optimized derivatives with potential synthetic feasibility which should be, by design, multikinase type II inhibitors for the target kinases.

Molecular docking study
Finally, molecular docking was used to investigate the ability of the designed compounds to bind to the binding sites of the target kinases accomplishing the key interactions responsible for the kinase inhibitory activity.The X-ray crystallographic structure of VEGFR-2, FGFR-1, and BRAF co-crystallized with Type II kinase inhibitors (PDB ID: 4ASD, 4V01 and 5CT7, respectively) were downloaded from the protein data bank [32,39,49,56] and were utilized to perform the molecular docking study.Details of the molecular docking procedures are discussed in the Additional file 1, Section 1: computational studies.

ADME properties prediction
SwissADME online web tool [57][58][59] was used to predict the physicochemical and AMDE properties of the target compounds 8a-u (See Additional file 1, Section 1: computational studies for further details).

General remarks
Chemicals, reagents and solvents were purchased from commercial suppliers.Chemical reactions were followed up by TLC using aluminium plates precoated with silica gel 60 F 245 (Merck).Uncorrected melting points were measured on a Stuart SMP30 melting point apparatus.Spectral and elemental analyses of 8a-u were recorded in the laboratory central services, National Research Centre, Cairo, Egypt and faculty of pharmacy, Cairo University.IR spectra (4000-400 cm −1 ) were detected on a Jasco FT/ IR 300 E Fourier transform infrared spectrophotometer.
The precipitated crude products 8a-u were filtered and were crystalized from EtOH to give the pure 2,5-diaryl benzimidazoles 8a-u. N

Screening of the inhibitory activity of 2,5-disubstituted benzimidazoles 8a-u on VEGFR-2
The VEGFR-2 inhibitory activity of 8a-u was investigated at 10 µM using VEGFR-2 assay kit (BPS Biosciences-San Diego-CA-US) according to the protocol provided by the manufacturer [68] and the % of inhibition was determined (For further details see Additional file 1: section II: practical results) [67].

Screening of the inhibitory activity of 8u on diverse kinases
The inhibitory activity of 8u was investigated at different concentrations using VEGFR-2, FGFR-1 and BRAF assay kits (BPS Biosciences-San Diego-CA-US) according to the protocol provided by the manufacturer (For further details see Additional file 1: section II: practical results).

Cell cycle analysis assay
The diaryl benzimidazole 8u was applied at its GI 50 concentration to MCF-7 cancer cells.After the cells were handled as previously described [69], the distribution of the cells at each stage of the cell cycle was analysed.(For further details see Additional file 1: section II: practical results).

Apoptosis assay
As previously reported, the Annexin V-FITC apoptosis detection kit (Abcam Inc., Cambridge, UK) in conjunction with two fluorescent channels flow cytometry was used to identify the populations of apoptosis and necrosis cells [69].(

Fig. 6
Fig.6 Suggested chemically feasible optimized derivatives 8a-u from the selected promising scaffold 8a

Fig. 7
Fig. 7 2D diagram (A) and 3D representation (B) of compound 8u showing its interaction with the VEGFR-2 active site

Fig. 8
Fig. 8 2D diagram (A) and 3D representation (B) of compound 8u showing its interaction with the FGFR-1 active site

Fig. 9
Fig. 9 2D diagram (A) and 3D representation (B) of compound 8u showing its interaction with the BRAF active site

Fig. 10
Fig. 10 SwissADME BOILED-Egg chart for the designed compounds 8a-u = 2.62-7.70µM).Meanwhile, incorporation of 3-methoxyphenyl group is favourable at the five position of 2-(4-methoxyphenylbenzimidazole) 8j and 2-(4-chlorophenylbenzimidazole) 8q (mean GI 50 = 5.24 and 5.40 µM, respectively).Interestingly, the introduction of 3,4,5-trimethoxyphenyl derivative is promising only in 2-(4-chlorophenylbenzimidazole) 8u (mean GI 50 = 7.98 µM).Encouraged by the potent activity of 8u on both the biochemical and cellular assays, it was further assessed for its effect on cell cycle and apoptosis of MCF-7 cell line.Interestingly, 8u was found to induce cell cycle arrest in MCF-7 cell line at the G2/M phase and accumulating cells at the sub-G1 phase as a result of cell apoptosis.

Fig. 12
Fig. 12 Structure activity relationship diagram showing the effect of the substitution pattern of 2-and 5-phenyl moieties of 8a-u on the antiproliferative activity

Fig. 13
Fig. 13 Effect of compound 8u on the phases of cell cycle of MCF-7 cells

Table 1
Distribution of the test set active inhibitors and inactive decoys for the target kinases VEGFR-2, FGFR-1,

Table 2
The collective assessment metrics of the generated pharmacophore models target kinases with predicted docking energy score ranges of − 14.89 to − 12.73 kcal/mol in VEGFR-2, − 14.18 to − 11.62 kcal/mol in FGFR-1, and − 13.65 to − 11.49 kcal/mol in BRAF, in comparison to the co-crystalized ligands docking score of − 15.19, − 17.00, and − 15.82 kcal/mol, respectively (See Additional file 1: Section 1: computational studies for further details).

Table 6
GI 50 of compounds 8a

Table 6
(continued) a Not detected
1H NMR (400 MHz; For further details see Additional file 1: section II: practical results).