In Silico and In Vitro Studies for Benzimidazole Anthelmintics Repurposing as VEGFR-2 Antagonists: Novel Mebendazole-Loaded Mixed Micelles with Enhanced Dissolution and Anticancer Activity.

Cancer is a leading cause of death worldwide and its incidence is unfortunately anticipated to rise in the next years. On the other hand, vascular endothelial growth factor receptor 2 (VEGFR-2) is highly expressed in tumor-associated endothelial cells, where it affects tumor-promoting angiogenesis. Therefore, VEGFR-2 is considered one of the most promising therapeutic targets for cancer treatment. Furthermore, some FDA-approved benzimidazole anthelmintics have already shown potential anticancer activities. Therefore, repurposing them against VEGFR-2 can provide a rapid and effective alternative that can be implicated safely for cancer treatment. Hence, 13 benzimidazole anthelmintic drugs were subjected to molecular docking against the VEGFR-2 receptor. Among the tested compounds, fenbendazole (FBZ, 1), mebendazole (MBZ, 2), and albendazole (ABZ, 3) were proposed as potential VEGFR-2 antagonists. Furthermore, molecular dynamics simulations were carried out at 200 ns, giving more information on their thermodynamic and dynamic properties. Besides, the anticancer activity of the aforementioned drugs was tested in vitro against three different cancer cell lines, including liver cancer (HUH7), lung cancer (A549), and breast cancer (MCF7) cell lines. The results depicted potential cytotoxic activity especially against both HUH7 and A549 cell lines. Furthermore, to improve the aqueous solubility of MBZ, it was formulated in the form of mixed micelles (MMs) which showed an enhanced drug release with better promising cytotoxicity results compared to the crude MBZ. Finally, an in vitro quantification for VEGFR-2 concentration in treated HUH7 cells has been conducted based on the enzyme-linked immunosorbent assay. The results disclosed that FBZ, MBZ, and ABZ significantly (p < 0.001) reduced the concentration of VEGFR-2, while the lowest inhibition was achieved in MBZ-loaded MMs, which was even much better than the reference drug sorafenib. Collectively, the investigated benzimidazole anthelmintics could be encountered as lead compounds for further structural modifications and thus better anticancer activity, and that was accomplished through studying their structure-activity relationships.

Introduction

Cancer is the second leading cause of mortality in the world, killing about 8 million people each year. Furthermore, it is a pity that cancer incidence is anticipated to rise by more than 50% in the next years.[1−4] In comparison to healthy cells, cancer cells frequently require more oxygen and nutrients, but it is far from cancer’s only requirement. Therefore, researchers are paying efforts and attempts to discover and develop new therapies for many cancer types.[5−9]

Furthermore, in adults, angiogenesis is essential during tissue growth, repair, and pregnancy. In addition, angiogenesis is a basic underlying process in the pathogenesis of several human diseases, including cancer.[10,11] Angiogenesis is a fundamental step in the turning of a benign tumor into a malignant one, where tumor masses are infiltrated by new blood vessels furnishing them with oxygen and nutrients to promote tumor growth and metastasis.[12] Since its discovery in 1983, vascular endothelial growth factor (VEGF) has been considered the most important regulator of blood vessel formation and is a key mediator of neovascularization in cancer.[10] Besides, VEGF can promote endothelial cell proliferation and motility. Hence, endothelial cells are abundantly represented in the malignant tissue.[5,13,14] Cancer cells are not capable of growing nor metastasizing when they cannot express VEGF. Notably, VEGFs are overexpressed in several types of cancers including hepatocellular carcinoma, colorectal cancer, and breast cancer.[12]

The biological duties of the VEGFs are mediated by a family of protein tyrosine kinase receptors called VEGF receptor tyrosine kinases (VEGFR-TKs) including VEGFR-1, VEGFR-2, and VEGFR-3.[15] VEGFR-1 and VEGFR-2 both enhance angiogenesis, whereas VEGFR-3 stimulation can induce lymphangiogenesis. Although VEGFR-1 has been shown to influence the function of VEGFR-2, VEGFR-2 has been revealed to mediate nearly all known VEGF cellular responses.[12]

When VEGF binds to VEGFR, it causes a conformational change in the receptor, which is followed by dimerization and phosphorylation of tyrosine residues. VEGF signaling through VEGFR-2 has been shown to play a key role in tumor angiogenesis regulation.[16,17] As a result, VEGF/VEGFR-2 signaling represents a promising therapeutic target in cancer treatment.[18−20] Therefore, inhibiting VEGFR-2 or downregulating its signaling is a pivotal strategy for developing novel drugs for a variety of human angiogenesis-dependent cancers.[21]

Additionally, concerning their common pharmacophoric features shared, VEGFR-2 antagonists showed four main structural features.[18,22,23] These features are as follows: (i) a flat heteroaromatic ring system containing at least one H-bond acceptor, (ii) a central aryl ring (spacer), (iii) a pharmacophore moiety contains two groups of an H-bond acceptor and an H-bond donor (e.g., amide or urea), and (iv) a terminal hydrophobic group.[12]

On the other hand, drugs that have been clinically approved or experimentally evaluated for diseases other than cancer but have been revealed to have unexpected cytotoxicity against malignant cells could be treated as good repurposed anticancer candidates.[24] Furthermore, several approved benzimidazole anthelmintic drugs were repurposed as anticancer agents.[25−29] It is worth mentioning that fenbendazole possesses an efficient antiproliferative activity. It was introduced as a novel microtubule interfering drug with antineoplastic action that could be investigated as a promising anticancer agent.[24] Besides, mebendazole was chosen for cancer cell suppression studies based on its pharmacokinetic properties.[30] In addition, the phosphorylation of Bcl-2, which results in dosage and time-dependent intrinsic apoptotic action in melanoma cells driven by microtubule depolymerization, is one of the mebendazole’s anticancer mechanisms.[30] Moreover, albendazole is a well-known FDA-approved anthelmintic medicine that is also cytotoxic to healthy cells and has been found as an anticancer agent. It was proved to have an extensive effect against human paclitaxel-resistant ovarian cancer cells.[31]

Accordingly, concerning the pivotal role of VEGFR-2 antagonism for cancer treatment and applying the repurposing approach (Figure ), our perspective in this article was targeting VEGFR-2 through virtual screening on a small library of some FDA-approved benzimidazole anthelmintics (1–13) which were found to have a great structural similarity to the co-crystallized benzimidazole urea inhibitor (14) of VEGFR-2 (PDB ID: 2OH4(32)), as depicted in Figures and 3. Also, we managed to enhance the solubility, dissolution, and the in vitro anticancer activity of one of the promising benzimidazole anthelmintics (MBZ) by preparing it in the form of mixed micelles (MMs) incorporating vitamin E.

Figure 1

Graphical representation for the repurposing of benzimidazole anthelmintic drugs as VEGFR-2 antagonists.

Figure 2

Chemical structures of FDA-approved benzimidazole-based anthelmintic agents; fenbendazole 1, mebendazole 2, albendazole 3, ricobendazole 4, cyclobendazole 5, oxibendazole 6, oxfendazole 7, dribendazole 8, parbendazole 9, bendazole 10, tiabendazole 11, triclabendazole 12, flubendazole 13, as well as the crystalline benzimidazole-urea inhibitor (14) as potent VEGFR-2 inhibitor.

Figure 3

Crystal structure of VEGFR-2 (PDB: 2OH4) with a benzimidazole urea inhibitor shows its essential binding interactions.

Accordingly, we carried out molecular docking[33] of the selected ligands on the 3D crystal structure of the mentioned VEGFR-2. Furthermore, molecular dynamics (MD) simulations were performed for the most promising members of the screened anthelmintics to confirm the obtained docking results and obtain more insights regarding the drug’s thermodynamic properties at the target receptor.[34,35]

Then, our proposed in silico studies were complemented through different in vitro studies on the best-screened candidates of benzimidazole anthelmintic drugs for their cytotoxicity and VEGFR-2 inhibiting potentials as well.

On the other hand, MBZ suffers from poor aqueous solubility which has been a major negative impact on its oral bioavailability. Therefore, a high dose of MBZ is administered to achieve a proper therapeutic effect, causing many adverse effects. High doses of MBZ cause anemia and liver damage.[36] Moreover, studies have shown evidence of teratogenic effects of MBZ in rats and mice. To overcome poor aqueous solubility and enhance the bioavailability of MBZ at a lower dose, MMs were prepared to contain a hydrophobic block that enables the poorly soluble drugs to be incorporated as an internal core and a hydrophilic block as a surrounding shell improving its activity.[37]

Rationale of the Work

Compounds with the benzimidazole carbamide scaffold have been introduced as potent inhibitors of the VEGFR-2 and endothelium-specific receptor tyrosine kinase (TIE-2).[32] The respective prototypical compound of such chemical class was albendazole being identified through initial high-throughput screening efforts and originally targeted as an anthelmintic agent.[38] The X-ray crystallographic structure of the top active benzimidazole-urea compound in complex with VEGFR-2 showed favored ligand anchoring at the ATP-specific site mediating polar contacts with Cys919 of the hinge region (Figure ). While the methoxycarbamido functional group of the crystalline ligand extended toward the solvent front region, the rest of the skeleton showed deep anchoring toward the back pocket region of the receptor with its urea moiety being sandwiched between Asp1046 and Glu885 at the Asp-Phe-Gly (DFG) motif of the VEGFR-2 activation loop. The compound exhibited potent in vitro activity against both VEGFR-2 and TIE-2 with IC50 values down to single-digit nanomolar concentrations.[32]

Owing to the great conserved structure of the ATP binding site of several human kinases’ receptors[39] as well as sharing the 2-methoxycarbamido-1H-benzimidazole core skeleton as with different anthelmintic agents, it was highly rationalized to explore the VEGFR-2 inhibition activity of market-released benzimidazole-based anthelmintic drugs acting as promising anticancer agents. The latter has driven us to determine the efficiency of FDA-approved benzimidazole-based anthelmintic drugs (Figure ) against VEGFR-2 using in silico and in vitro approaches, while the preliminary data obtained from this study would guide further optimization within the future based on structure–activity relationships (SARs) studied being attaining. Adopting this drug repurposing approach possesses the advent of assuring medical safety because these FDA-approved benzimidazole-based anthelmintic drugs have already been tested in animal models, undergone all the essential clinical trials, and are already infrastructured toward manufacture at large scales.[40,41]

Results and Discussion

Docking Studies

Molecular docking stimulation of fenbendazole 1, mebendazole 2, albendazole 3, ricobendazole 4, cyclobendazole 5, oxibendazole 6, oxfendazole 7, dribendazole 8, parbendazole 9, bendazole 10, tiabendazole 11, triclabendazole 12, flubendazole 13, and the benzimidazole urea inhibitor 14 into the active site of VEGFR-2 receptor was done. They got stabilized at the VEGFR-2 binding site by variable several electrostatic bonds. The order of strength of binding was as follows: the benzimidazole urea inhibitor (14, docked) > fenbendazole (1) > mebendazole (2) > albendazole (3) > ricobendazole 4 > cyclobendazole 5 > oxibendazole 6 > oxfendazole 7 > dribendazole 8 > parbendazole 9 > bendazole 10 > tiabendazole 11 > triclabendazole 12 > flubendazole 13.

Most drugs showed nearly binding modes similar to the benzimidazole urea co-crystallized inhibitor. Many poses were obtained with better binding modes and interactions inside the receptor pocket. The poses with the most acceptable scores (related to the stability of the pose) and rmsd_refine values (related to the closeness of the selected pose to the original ligand position inside the receptor pocket) were selected. Results of scores and different interactions with the amino acids of protein pocket are depicted in Table .

Table 1

Receptor Interactions and Binding Energies of the Tested Anthelmintic Drugs (1–13) and the Docked Benzimidazole–Urea Inhibitor (14) at the Binding Site of the VEGFR-2 Receptor

no.	anthelmintic drug	Sa	rmsdb	amino acid bonds	distance (A)
1	fenbendazole	–8.18	1.08	Leu838/H2O bridged H bond	2.74
				Leu838/H2O bridged H bond	2.89
				Cys917/H-donor	2.96
				Cys917/H-acceptor	3.12
				Leu838/pi-H	3.98
				Val846/pi-H	4.69
2	mebendazole	–8.12	0.93	Leu838/H2O bridged H bond	2.68
				Leu838/H2O bridged H bond	2.93
				Cys917/H-donor	2.97
				Cys917/H-acceptor	3.08
				Leu838/pi-H	4.03
3	albendazole	–7.91	1.15	Leu838/H2O bridged H bond	2.74
				Leu838/H2O bridged H bond	2.90
				Cys917/H-donor	3.01
				Cys917/H-acceptor	3.12
				Leu838/pi-H	4.01
4	ricobendazole	–7.85	1.31	Leu838/H2O bridged H bond	2.77
				Cys917/H-donor	2.89
				Leu838/H2O bridged H bond	2.99
				Cys917/H-acceptor	3.03
				Leu838/pi-H	4.09
5	cyclobendazole	–7.60	1.66	Leu838/H2O bridged H bond	2.83
				Cys917/H-donor	2.86
				Leu838/H2O bridged H bond	2.97
				Cys917/H-acceptor	3.02
				Leu838/pi-H	4.06
6	oxibendazole	–7.58	1.05	Cys917/H-donor	2.80
				Leu838/H2O bridged H bond	2.92
				Cys917/H-acceptor	3.09
				Leu838/pi-H	4.03
7	oxfendazole	–7.54	1.69	Leu838/H2O bridged H bond	2.75
				Cys917/H-acceptor	3.63
				Leu838/pi-H	4.04
8	dribendazole	–7.53	1.31	Cys917/H-donor	2.85
				Leu838/H2O bridged H bond	2.91
				Cys917/H-acceptor	3.06
				Leu838/pi-H	4.03
				Val846/pi-H	4.82
9	parbendazole	–7.46	1.06	Cys917/H-donor	2.80
				Leu838/H2O bridged H bond	2.92
				Cys917/H-acceptor	3.05
				Leu838/pi-H	4.06
				Val846/pi-H	4.72
10	bendazole	–6.84	0.70
11	tiabendazole	–6.62	0.59	Leu838/H2O bridged H bond	2.82
				Leu838/H2O bridged H bond	2.92
				Cys917/H-acceptor	3.30
				Leu838/pi-H	3.97
12	triclabendazole	–6.40	1.18	Leu838/H2O bridged H bond	2.83
				Cys917/H-acceptor	3.22
				Leu838/pi-H	3.89
13	flubendazole	–5.94	1.05	Leu838/H2O bridged H bond	2.72
				Leu838/H2O bridged H bond	2.88
				Cys917/H-donor	2.96
				Cys917/H-acceptor	3.09
				Leu838/pi-H	4.06
				Val846/pi-H	4.78
14	benzimidazole urea inhibitor (GIG)	–11.15	1.39	Cys917/H-donor	2.82
				Glu883/H-donor	2.86
				Glu883/H-donor	2.89
				Leu838/H2O bridged H bond	2.92
				Asp1044/H-acceptor	2.94
				Cys917/H-acceptor	3.10
				Leu838/pi-H	4.02

S: score of a compound inside the binding pocket of protein (kcal/mol).

rmsd/refine: root-mean-square deviations among heavy atoms of crystallized structure (prerefinement) and those of the obtained binding mode (postrefinement).

By analyzing docking results of the selected anthelmintics, it was found that most of the selected compounds manifested very close binding scores and modes compared to the co-crystallized inhibitor (GIG) at the VEGFR-2 target receptor. FBZ (1), MBZ (2), and ABZ (3) were found to have the best binding affinities and modes against VEGFR-2 with binding scores of −8.18, −8.12, and −7.91 kcal/mol, respectively (Table ). These energy values were close to that of the docked GIG inhibitor (binding energy = −11.15 kcal/mol). When we put eyes on the binding interactions of FBZ (1) with the VEGFR-2 target receptor, we could reveal that it showed two water bridged hydrogen bonds with Leu838, interacted with Cys917 through both hydrogen bond donor and acceptor, and bound to Leu838 and Val846 via pi–H bonds. However, regarding the binding interactions of MBZ (2) and ABZ (3), both revealed that they are bound to Leu838 through two water bridged hydrogen bonds, interacted with Cys917 through hydrogen bond donor and acceptor, and bound to Leu838 via a pi–H bond as well.

The detailed binding modes of the docked benzimidazole urea inhibitor (GIG, 14) and all of the tested anthelmintics (1–13) are presented in Table . Nevertheless, the 3D-protein positioning, as well as 3D-binding interactions of the best selected three anthelmintics (1, 2, and 3), is described in Table .

Table 2

3D View of Binding Interactions and the 3D Positioning of the Best-Docked Benzimidazole Anthelmintic Drugs (1, 2, and 3) within the VEGFR-2 Receptor Pocket (PDB: 2OH4) Compared to the Benzimidazole Urea Inhibitor (Docked, 14)

MD Simulations

For gaining an in-depth understanding of the relative stability and thermodynamic behavior of the assigned ligand–protein complexes across a specific time, a MD investigation was conducted. The latter computational tool is considered particularly beneficial for investigating the conformational spaces of ligand–protein complexes which are considered more advantageous over other different computational techniques including mechanics-energy minimization and molecular docking approaches concerning the sole analysis of static images.[42] Exhibiting relevant ligand–protein docking interactions, the top-docked binding modes of the benzimidazole-based anthelmintic compounds (FBZ, MBZ, and ABZ) in addition to the crystalline reference inhibitor (GIG), at the VEGFR-2 active site, were enrolled within a 200 ns explicit MD simulation for comparative dynamic investigation.

Analyzing the Stability of Protein–Ligand Complexes

Under 200 ns MD simulation, all investigated anthelmintic compounds showed significant thermodynamic stability at the hVEGFR-2 active site which was confirmed via monitoring the root-mean standard-deviation (rmsd) trajectories. Typically, the rmsd determines the deviation of the molecules in relation to a designated reference/original structure. This analytical parameter would indicate ligand–protein stability and confirm the validity of the MD simulation protocol being adopted. High rmsd values are indicative of significant conformational changes and target instability.[43] However, high complex rmsds are correlative for minimal ligand–protein affinity having the ligand unable to be maintained at the protein’s pocket throughout the MD simulation interval.[44]

In reference to the backbone, the obtained VEGFR-2 protein’s rmsds showed overall typical MD simulation behaviors (Figure A). At the beginning of MD runs and over the initial 20 ns frames, the protein backbone rmsds elevated due to the release of constraining applied during the prior minimization and equilibration stages. After the initial 20 ns, steady protein rmsds were depicted over more than half the MD run (>150 ns), except for limited fluctuations being illustrated at FBZ-associated protein across the 50–60 ns time-frames. Nearly all simulated proteins leveled off at almost similar rmsds throughout the equilibration plateau and until the MD simulation ends. Comparable rmsd trajectories were depicted for ABZ-, MBZ-, and GIG-bounded VEGFR-2 protein, following their respective equilibration with average values of 3.34 ± 0.18, 3.59 ± 0.37, and 3.96 ± 0.23 Å, respectively. However, slightly higher values were assigned for FBZ-bound protein (5.12 ± 0.31 Å), being correlated to its depicted limited fluctuations as well as late equilibration following the 50 ns of the MD simulation timeframes. The ABZ-bound VEGFR-2 protein managed to achieve the steadiest rmsd trajectories, illustrating minimal standard deviation (SD) following the attained equilibrium. All described VEGFR-2 protein thermodynamic behaviors emphasize the successful target protein convergence across the designated MD simulation timeframe. Additionally, the above-described protein rmsds inferred the adequacy of the prior system minimization/relaxation and thermal equilibration as well as the adopted 200 ns MD runs that required no further extensions.

Figure 4

rmsd trajectory analysis of the examined anthelmintic compounds and reference inhibitor in bound with VEGFR-2 target across the 200 ns explicit MD runs. (A) Protein’s backbone-rmsds; (B) ligand–protein complex backbone-rmsds; and (C) only ligand backbone-rmsds (Å), along MD timeframe (ns).

To insure the maintenance of the simulated ligand at the VEGFR-2 ATP-binding site, the combined ligand–protein complex rmsd deviations were monitored having the protein backbones as their reference frame (Figure B). All simulated VEGFR-2/ligand complexes successfully attained their respective thermodynamic stability state showing backbone rmsd plateau despite limited fluctuations. Although differential backbone rmsds were shown across the initial MD frames, each of the ABZ, MBZ, and GIG models was capable to converge across the second 100 ns timeline, achieving the final rmsd value at 3.80 ± 0.33 Å. Nonetheless, the latter dynamic behaviors were nontypical for the FBZ–protein complex because the depicted model showed two equilibration plateau; an initial equilibration around 20–40 ns and a latter one beyond 50 ns until the end of the MD simulation run. Such a dynamic behavior suggested significant FBZ conformation change beyond 50 ns. Nevertheless, the FBZ ligand itself was suggested to be confined within the target’s binding site because a slight complex rmsd shift was depicted between both equilibria (3.87 ± 0.23 → 4.49 ± 0.18 Å) and such attained values were not higher than those depicted for the protein’s rmsds. On the other hand, compounds ABZ, MBZ, and crystalline reference GIG achieved earlier equilibrations of steadier complex rmsd trajectories and lower comparative average values (3.30 ± 0.29, 3.62 ± 0.29, and 3.97 ± 0.30 Å, respectively). All the above findings highlight the significant ligand’s retainment of the investigated anthelmintic agents within the protein active site, which was highly comparable to the crystalline GIG potent VEGFR-2 inhibitor.

Further analysis was proceeded through monitoring the only ligand’s rmsd tones in relation to the reference backbone frame of the protein which was considered as a relevant stability indicator of ligand’s confinement within the pocket as well as simulated protein convergence (Figure C). Interestingly, both MBZ and GIG showed the steadiest and lowest average ligand rmsd trajectories across the whole MD simulation runs (1.52 ± 0.29 and 1.68 ± 0.22 Å, respectively) inferring limited ligand conformational changes within the VEGFR-2 binding site. Concerning the ABZ compound, limited fluctuations were depicted at the initial MD simulation frames where subsequently the rmsds reached equilibration maintaining the initial ligand’s conformation/orientation across the same timeframes where the respective VEGFR-2 protein showed successful convergence (beyond the 20 ns). Notably, the average ABZ ligand rmsd value following equilibration (1.29 ± 0.29 Å) was nearly comparable to those of both MBZ and GIG suggesting further stability of these ligand-docked poses within the target site. In concordance with the above FBZ–protein rmsd trajectories, the ligand showed stable trajectories until the 50 ns time frames where after that the rmsds were elevated and attained stable trajectories until the end of the MD simulation. Such a dynamic behavior ensures the significant change of the FBZ compound attaining a second conformation/orientation following the 50 ns simulation run. Notably, the rmsd tones of each VEGFR-2 protein did not exceed the 1.5-fold of those of their corresponding ligand in complex, which further confirms the successful convergence of the ligand–protein complexes and infers the suitability of 200 ns MD simulation runs requiring no further extension.

The ligand–protein global stability was further investigated through monitoring both the radii of gyration (Rg) and solvent-accessible surface area (SASA) trajectories of the complex entities along with the whole MD timeframe. Typically, the estimated radii of gyration of the investigated complexes permitted the exploration of the complex rigidity and compactness where this stability parameter accounts for the complex’s mass-weighted root-mean-square distance relative to its common mass center. In these regards, minimal Rgs with achieved plateau around the average value would be correlated to the sustained stability/compactness of an investigated complex.[45] On similar bases, decreased SASA trajectories confer ligand–protein structural shrinkage due to the influence of solvent-surface charges, which would lead to more conformational compactness and stability. The latter has been correlated to the SASA calculation which estimates the molecular surface area being assessable to solvent molecules, providing a quantitative measurement of the complex-solvent interaction.[46]

Herein, the steadiest Rg trajectories were assigned for the crystalline ligand, GIG, in complex with VEGFR-2, showing an average value of 20.65 ± 0.13 Å (Figure A). Concerning both FBZ and MBZ, the complexes seemed to be expanding at the initial MD simulation frames, with the MBZ–protein complex being at higher Rg tones (max values 21.08 and 21.59 Å, respectively). However, both systems achieved compactness and significant contraction following the 80 ns and until the MD ends. The ABZ–protein complex, on the other hand, showed initial compactness for more than half the MD runs (20.68 ± 0.13 Å), while after that the complex seemed to be expanding until the next 40 ns where then it finally attained significant compactness with the lowest Rg trajectories among all investigated molecules (20.16 Å). Interestingly, the four investigated ligand–protein complexes converged around similar Rg values (20.68 ± 0.19 Å) at 200 ns, confirming the comparable compactness and stability profiles of the three anthelmintic models as well as the crystalline potent inhibitor. Significant intra-/intermolecular interactions among the investigated complexes were further confirmed through the calculated complex SASA tones (Figure B). Following the 60 ns MD timeframe, comparable complex SASA trajectories were assigned for the ABZ– and MBZ–protein complexes as well as the crystalline inhibitor (172.82 ± 3.46, 173.07 ± 2.77, and 172.19 ± 3.20 nm2, respectively), conferring preferential ligand confinement within the VEGFR-2 binding site. Notably, the FBZ–protein complex illustrated significantly reduced SASA tones (165.78 ± 3.23 nm2) following 50 ns and until the end of the simulation run, suggesting better ligand–protein interactions, particularly for the second ligand’s conformation, because binding is a solvent-substitution process.

Figure 5

Global stability profiles of obtained Rg and SASA tones of the examined anthelmintic compounds and reference inhibitor complexed with VEGFR-2 target across the 200 ns explicit MD run. (A) Complex Rg (Å); (B) complex SASA tones (nm2), along MD timeframe (ns).

Because the presented rmsd analysis highlights the significant ligand-target stability for the examined ligands, it was beneficial to further investigate the local protein flexibility and how this could be contributed to the ligand–protein binding. The fluctuation of VEGFR-2’s residues was monitored by estimating the RMS-fluctuation (RMSF) stability validation parameter being able to highlight the residue-wise contribution within the target protein stability. Typically, RMSF provides a valuable evaluation of the target’s residue dynamic behavior represented as both fluctuation and flexibility, through estimating the average deviation of each protein’s amino acid in relation to its respective reference position across time.[47] Within the presented article, the difference RMSF (ΔRMSF) was a better estimation of the protein local flexibility being the RMSF difference for each ligand-bound protein relative to the VEGFR-2 apo state (ΔRMSF = apo RMSF – holo RMSF). Adopting a ΔRMSF cutoff value of 0.30 Å was relevant for estimating the significant alterations within the protein’s structural movements, meaning that amino acids with ΔRMSF above 0.30 were considered of limited mobility.[48] Investigating the RMSF trajectories essentially execute for a trajectory region considered stable. Based on the above protein’s rmsd analysis (Figure A), the VEGFR-2 proteins targets were of significant conformational stability along the 200 ns MD simulations for all systems despite the limited fluctuations for the FBZ system. Therefore, the backbone RMSF calculations were reasoned to be estimated across the whole MD simulation trajectories.

Throughout the ΔRMSF analysis, the free terminal residues showed a typical fluctuation pattern with the highest negative ΔRMSF values in comparison to the core residues (Figure ). The latter could be reasoned for the lower extent of intramolecular interactions among the terminal residues as compared to those of core amino acids which have been considered ideal behavior throughout MD simulation runs. Notably, patterns of high fluctuation were illustrated for all ligand-bound VEGFR-2 amino acids near the N-terminal relative to those settled near the carboxy end (average −3.90 ± 0.51 vs 0.12 ± 1.43 Å). Almost all residue ranges within the C-lobe and N-lobe of the kinase domain illustrated significant immobility depicting positive ΔRMSF values. Interestingly, the residue ranges Ser875–Leu880 of the N-lobe and Lys1053–Arg1059 at C-lobe showed the highest immobility profiles with ΔRMSF up to 4.13 and 4.74 Å, respectively. Such dynamic behavior confers significant influence of ligand’s binding upon the stability of these residue ranges (particularly at C-lobe) or in other terms the pivotal role of these residues for the ligand stability at the VEGFR-2’s active site. This came in great agreement with several reported studies investigating the potential antiangiogenic VEGFR-2 inhibition activity of naturally occurring metabolites and chemical library deposits.[49−54] However, it is worth noting that both residue ranges were proven to possess relatively conserved hydrogen bond interactions among the constituting residues as well as with the binding inhibitors.[55] On the contrary, the residue range Ala943–Glu993 was of the most flexible pattern (ΔRMSF down to the highest negative values −9.92 Å) inferring the negligible contribution of such residues within the ligand–protein interactions. This was also consistent with the reported data by where the residue range of Thr940–Glu989 possessed an insignificant impact upon the catalytic activity of the VEGFR-2 protein.[56]

Figure 6

Difference RMSF analysis across VEGFR-2 residues bound to the investigated anthelmintic compounds and reference ligand across the 200 ns MD runs. The protein’s backbone-ΔRMSFs were determined considering independent 200 ns MD run for holo VEGFR-2 states (in complex with investigated ligands or crystalline reference inhibitor, GIG) against the unliganded/apo state (PDB ID: 1VR2). Trajectories of ΔRMSF are illustrated as functional residue numbers (residues 814-up to-1169).

Concerning the comparative local stability of ligand-bound VEGFR-2 proteins, trends of high negative ΔRMSF values were depicted for FBZ-bound amino acids as compared to those of other investigated anthelmintic agents as well as reference inhibitors. These trends were recognized along with several VEGFR-2 residue ranges while being most noted for the flexible Ala943–Glu993 range as well as Asp1062–Lys1068 residue range at the C-lobe domain. On the other hand, the highest residue-wise stability and rigidity profiles were assigned for the ABZ being comparable to that of the crystalline potent inhibitor, which further confirms the relevant stability of such systems as being previously discussed within the rmsd findings.

Further investigation of ligand-VEGFR-2 interactions was proceeded through comparative studying the furnished ΔRMSF values and specific flexibility of the pocket’s canonical lining amino acids. Interestingly, several key pocket residues, as well as vicinal residues, showed significant inflexibility having ΔRMSF values beyond the threshold index 0.30 Å (Table ). The hinge region residues depicted moderate immobility profiles with both Cys917 and Lys918 represented as the most rigid ones among such residue range, exhibiting comparable ΔRMSF values (0.76 ± 0.10 Å) across the four investigated ligands. The latter infers the important role of these two hinge residues for anchoring the polar carbamoyl of anthelmintic agents as well as the urea scaffold of the crystalline ligand through hydrogen bond interactions. Both the hydrophobic back-pocket and ATP-association cleft showed moderate rigidity profiles for all ligands. Such a significant inflexibility pattern highlights the preferential and deep anchoring of all anthelmintic compounds within the binding site of the kinase domain. Despite being short-length ligands, the anthelmintic agents managed to achieve a relevant residue-wise stability profile for the deep hydrophobic back pocket being comparable to that of the crystalline potent inhibitor. Such important observation further confirms the potential inhibition activity of these anthelmintic agents.

Table 3

Estimated ΔRMSFa Values for Ligand-VEGFR-2 Proteins across the 200 ns MD Runs

canonical subsites	comprising residues	FBZ	MBZ	ABZ	GIG
hinge region	Val914	0.35	0.33	0.34	0.438
	Glu915	0.46	0.45	0.45	0.50
	Phe916	0.47	0.45	0.45	0.44
	Cys917	0.72	0.65	0.66	0.66
	Lys918	0.91	0.81	0.84	0.80
ATP-binding cleft	Leu838	1.34	1.19	1.32	1.27
	Ala864	1.15	0.98	1.07	1.05
	Val865	0.78	0.63	0.76	0.68
	Lys866	0.70	0.60	0.67	0.61
	Glu883	0.78	0.85	0.92	1.13
	Gly920	1.15	1.01	1.15	1.16
	Leu1033	0.64	0.58	0.71	0.67
hydrophobic deep back pocket	Ile886	0.39	0.52	0.54	0.77
	Leu887	1.04	0.99	1.01	1.16
	Ile890	0.96	0.87	0.82	0.99
	Val896	0.62	0.53	0.66	0.68
	Val897	0.45	0.28	0.47	0.46
	Val912	0.34	0.28	0.16	0.50
	Leu1017	0.36	0.31	0.47	0.41
	Cys1043	0.59	0.54	0.68	0.66
	Ile1042	0.82	0.19	0.77	0.82
glycine-rich region	Gly841	1.28	1.11	1.23	1.26
	Ala842	1.31	1.19	1.26	1.31
	Phe843	1.49	1.42	1.44	1.49
	Gly844	1.14	1.01	1.09	1.11
	Gln845	1.39	1.04	1.41	1.23
	Val846	2.11	1.88	2.11	2.03
catalytic loop	Cys1022	0.48	0.46	0.65	0.69
	Ile1023	0.50	0.38	0.61	0.70
	His1024	0.33	0.28	0.38	0.41
	Arg1025	0.33	0.43	0.55	0.48
	Asp1026	0.30	0.31	0.34	0.31
activation loop DFG motif (Asp–Phe–Gly) and vicinal residues	Cys1043	0.42	0.09	0.45	0.43
	Asp1044	0.27	0.06	0.26	0.22
	Phe1045	0.45	0.23	0.42	0.49
	Gly1046	0.29	–0.01	0.29	0.29
	Leu1047	0.46	0.08	0.44	0.54
	Ala1048	0.99	0.13	0.90	1.07
	Arg1049	0.69	–0.39	0.66	0.74

ΔRMSF values were determined for each ligand-associated VEGFR-2 protein in relation to unliganded/apo state (PDB ID: 1VR2). ΔRMSF > 0.30 Å immobility threshold was set, where amino acids with significant inflexibility/immobility profiles showed values above this designated cutoff.

Interestingly, the glycine-rich region exhibited the highest immobility profiles with ΔRMSF values up to 2.03 ± 0.11 Å being assigned for the Val846 hydrophobic residue. The latter could be reasoned for the conserved hydrogen bonding among the region comprising residues as well as the close proximity of Val846 toward the ATP-binding site inferring its pivotal role for the stability of the investigated ligand and crystalline inhibitor.[39] Moving toward the significant loops near the canonical substrate binding site, higher flexibility trends were assigned for residues of the DFG motif of the activation segment rather than those of the catalytic loop. Notably, the ΔRMSF values for the activation loop DFG motif were lower than those of any of the binding site subsites. Residues including Asp1044 and Gly1046 as well as Arg1049 only for the FBZ–protein model were below the 0.3 Å threshold or even at negative values inferring their insignificant role in stabilizing the binding ligands. Interestingly, significantly higher flexibility was assigned for activation loop residues in complex with FBZ which may be due to the fact that the ligand exhibited conformational shift throughout the MD simulation. Based on the ΔRMSF furnished stability results across different pocket subsites, comparable dynamic behaviors were depicted for the three anthelmintic agents in relation to the crystalline inhibitor. Despite adopting altered conformation/orientation, the FBZ-bound protein residues exhibited almost comparable inflexibility and stability profiles to those of ABZ, MBZ, and the crystalline inhibitor (GIG). All of which are in good agreement with the above-described MD behaviors illustrated through the analysis of Rg, SASA, and rmsd findings. However, examining the differential conformations of investigated ligand-VEGFR-2 complexes as well as time-evolution of ligand-VEGFR-2 binding interactions at selected frames would provide valuable information regarding the nature and conformationally directed ligands’ affinity toward the VEGFR-2 binding site.

Conformational and Intermolecular Binding Analysis

Analysis of key conformational alterations across the MD simulation timeframe was performed through examining the ligand–protein models at trajectories of the start and final timeframes. For each ligand-VEGFR-2 model, frames at 0 and 200 ns were extracted and minimized to a gradient of 1 × 10–3 kcal mol–1 A–2 using the MOE-package. Stable ligand-target binding profiles were depicted for almost all investigated anthelmintic agents and reference inhibitors. Interestingly, both ABZ and MBZ showed deeper anchoring into the VEGFR-2 pocket at the MD simulation end (Figure A,B). However, the initial polar interaction between both ligand’s benzimidazole carbamoyl scaffold and the main chain of Cys917 (CO and NH) were maintained at overall frequencies of 75.51 and 77.34% for ABZ and MBZ, respectively, across the MD simulation. Notably, the occupancy of the Cys917 (C=O) mainchain hydrogen bond pair is always around twofold higher than that of the Cyst917 (NH) mainchain functionality. Several hydrophobic pocket residues, including Leu838, Val846, Val897, Phe916, Leu1033, and Phe1045, maintained the sandwich-like conformation around the linear hydrophobic skeleton of both ABZ and MBZ until 200 ns MD runs. This was observed because these residues kept short distances from both ligand’s structure along the 200 ns simulation run as being presented within the given heatmap (Figure A,B lower panels). On the other hand, few ionizable residues such as Lys866, Glu883, and Asp1044 might have imposed unfavored repulsion forces with the ligand’s hydrophobic terminal. This was particularly obvious with MBZ where the terminal benzoyl of MBZ showed 45° rotation from its initial docked position to minimize steric hindrances as well as unfavored contacts with Lys866 and Glu883 (Figure B).

Figure 7

Conformational and intermolecular distance analysis of simulated ligand–protein complexes. (A) FBZ; (B) MBZ; (C) ABZ; and (D) GIG. Upper panels are overlaid snapshots of the ligand-VEGFR-2 complexes at 0 and 200 ns of MD runs. The VEGFR-2 proteins are illustrated in red and green 3D representation (cartoon) relative to the last and initially extracted frames, respectively. Both ligands (represented as sticks) and polar hydrogen bond interactions (dashed lines) are colored in correspondence to their respective extracted frames. Lower panels represent the heatmap representation of the time evolution of intermolecular distances between binding ligand and protein residues functionalities during the whole MD simulation. Polar interaction (hydrogen bonding) distances were measured between the designated interacting donor-H···acceptor, while the distances of hydrophobic interactions were measured from the nearest interacting atom of the ligand to the Cα of a particular residue. The values of the intermolecular distances were conditionally formatted through a color scale from 0 Å (green) and up to 10 Å (white) using the Microsoft EXCEL spreadsheets.

The third anthelmintic compound, FBZ, showed a significant conformational shift for its terminal aromatic moiety, yet kept the rest of the compound without relevant conformational/positional changes (Figure C). The latter could explain the second equilibration plateau beyond the 50 ns at the previously described rmsd tones (Figure ). Interestingly, the initial close proximity of the phenyl substitution at Lys866 might have caused such depicted orientation flip, allowing a more favored position where the aromatic group was at an ideal T-shaped π–π stacking with Phe1045. However, such a new orientation allowed the phenyl substitution to face another ionizable residue, Arg1049, yet stabilization imposed by Phe1045 might overcompensate any depicted ring-residue repulsions. Finally, the great superimposition of the FBZ’s benzimidazole carbamoyl scaffold at initial and 200 ns allowed stronger hydrogen bond pairing with Cys917 main chain with slightly higher frequencies (80.17%) as compared to other anthelmintic compounds. Moving toward the crystalline inhibitor, minimal conformational changes were depicted for the ligand’s core skeleton, while the terminal aromatic moiety was the most being altered (Figure D). The hydrogen pair with Cys917 was also maintained, yet at a much higher frequency (91.61%) with the preferentiality for the carbonyl mainchain functionality of this hinge residue. Another polar hydrogen bond pairing was deduced through analyzing the running trajectories including urea carbonyl: Asp1044 mainchain (53.20%), urea NH: Glu883 sidechain (3.35%), and fluoro substitutions: His1024 NH sidechain (0.25%) hydrogen bond pairs. The heat map for ligand-residue intermolecular distance showed the much more close-range hydrogen-donor acceptor distances for Asp1044 polar interactions over those of Glu883 and His1024 across the simulation where both weak bonds were lost by the end of the MD simulation. Finally, the crystalline ligand showed close range hydrophobic distances with more pocket lining nonpolar residues which further stabilize the ligand at the target pocket.

Free Binding Energy Calculations

The calculation of the binding free energy was performed to understand the ligand-VEGFR-2 interaction nature, investigate comparative ligand-directed pocket affinities, and obtain more insights regarding individual residue-wise energy contributions.[57] In these regards, the molecular mechanics Poisson/Boltzmann’s surface area (MM/PBSA) calculations were used to estimate the binding-free energies, where lower negative values confer less ligand’s affinity toward the target binding site.[58] Accuracy of the MM/PBSA calculation is considered comparable to the free-energy perturbation methods, yet with the advantage of lower computational expenditures.[58] Both single trajectory approach and only SASA model (ΔGTotal = ΔGMolecular Mechanics + ΔGPolar + ΔGAPolar) were adopted. Extracted and saved representative frames along the entire 200 ns MD runs were used for estimating each term of the binding energy as well as their average values. The latter was reasoned because ligand–protein complexes rmsds rapidly achieved the convergence/equilibration plateau after few initial MD trajectories (Figure B).

To our delight, all simulated anthelmintic agents illustrated significant binding-free energies and in turn affinities toward the VEGFR-2’s binding site (Table ). Comparable binding-free energies at significant negative values were assigned for the three anthelmintic agents; −105.28 ± 11.82, −99.52 ± 23.18, and −96.28 ± 39.37 kJ/mol for FBZ, MBZ, and ABZ, respectively. Interestingly, the latter binding free energy pattern came in great concordance with the preliminary docking investigation showing preferential higher docking scores for FBZ, followed by that of MBZ and then ABZ. Similarly, higher binding-free energy was depicted for the crystalline reference inhibitor being at almost 1.5-fold those of the investigated ligands. This preferential binding affinity highlights the significant impact of the fluorinated aromatic carbamide tail substitution in guiding the ligand-VEGFR-2 binding and mediating preferential ligand anchoring at the deep target pocket. Considering the cited potent inhibitory bioactivity of the crystalline compound against VEGFR-2 and serving as a relevant antiangiogenic agent,[32] the nearly matchable binding affinity data for these anthelmintic agents imply their promising activities toward this same protein target.

Table 4

Total Free Binding Energies (ΔGTotal binding ± SD) and Dissected Energy Terms Regarding Investigated Anthelmintic Agents and Reference Inhibitor at the VEGFR-2 Active Site

	ligand–protein complex
energy (kJ/mol ± SD)	FBZ	MBZ	ABZ	GIG
ΔGvan der Waal	–158.88 ± 10.57	–164.78 ± 35.98	–155.90 ± 16.17	–247.06 ± 49.17
ΔGElectrostatic	–23.68 ± 35.00	–29.06 ± 44.78	–23.72 ± 21.09	–36.04 ± 52.00
ΔGSolvation; Polar	93.36 ± 47.82	110.30 ± 38.65	98.82 ± 38.73	160.76 ± 49.30
ΔGSolvation; SASA	–16.08 ± 1.91	–15.95 ± 2.46	–15.48 ± 1.63	–25.07 ± 4.78
ΔGTotal binding	–105.28 ± 11.82	–99.52 ± 23.18	–96.28 ± 39.37	–147.41 ± 19.77

The furnished free binding energies illustrated the dominance of van der Waal potentials within the binding energy calculations of each simulated ligand. Higher electrostatic and lipophilic energy contributions were depicted for the crystalline reference over those of the anthelmintic agents due to possessing profound aromatic characteristics and hydrogen bond features (high numbers of hydrogen bong acceptors) at its terminal substitutions. This was highly reasoned because the crystallized ligand depicted a higher hydrogen bond number/frequency with VEGFR-2 pocket residues than any of the investigated anthelmintic agents. It is worth mentioning that the total nonpolar interactions (ΔGvan der Waal plus ΔGSASA) were higher for the crystalline inhibitor as compared to the anthelmintic compounds, which could be correlated to the large surface area of the VEGFR-2 pocket. Interestingly, several reported data have considered the pocket of VEGFR-2 as deep and more hydrophobic than other comparable targets possessing conserved nonpolar residues lining the target’s pocket.[32,39,56] Being with a large surface area and hydrophobic, the VEGFR-2 pocket would favor higher hydrophobic interactions with GIG as this ligand can attain more extended structural conformations at the target’s binding site. On the other hand, deep anchoring of anthelmintic agents at the target pocket while attaining relevant could be achieved through furnishing favored strong hydrophilic binding interactions with VEGFR-2 pocket key amino acids. This was clearly presented throughout the above illustrated intermolecular bonding as well as conformational analyses.

All furnished data confirm the important role of hydrophobic/polar combined functional groups at the terminal fluorinated aromatic carbamide tail substitutions for anchoring the reference ligand at the VEGFR-2 pocket. However, these particular polar functional moieties might serve as double-bladed influencers upon ligand-target binding. This could be reasoned that such functional groups impose higher ΔGsolvation for the GIG–protein system, compromising the ligand’s anchoring as the binding process is considered a solvent-substitution one. Therefore, these anthelmintic compounds could be optimized through future developmental approaches via introducing ionizable groups with relevant lipophilic characters (e.g., tetrazole ring) has been considered significant for ΔGsolvation reduction, ligand-VEGFR-2 binding extension, and in turn potential target inhibition.

Free binding energy decomposition via g_mmpbsa scripts for identifying the key amino acids enrolled in the furnished free binding energies.[58] Interestingly, almost comparable residue-wise energy contributions were depicted for three anthelmintic agents except for limited differences regarding the magnitude of energy contribution (Figure A–C). This was expected as these benzimidazole agents showed close values of the dissected energy terms and total free binding energies. The significant residues showing favored contribution (high negative values) within the ligand–protein binding energy include the following: Val914, Phe916, and Cys917 (hinge region); Leu838, Ala864, Gly920, and Leu1033 (ATP-binding cleft); Val897 and Cys1043 (hydrophobic back pocket); Val846 of the glycine-rich region; and finally Cys1043, Phe1045, and Leu1047 at the activation loop. On the other hand, the crystalline ligand GIG shared many of the residue-wise energy contributions being comparable to those of the anthelmintic compounds, however, of significantly higher magnitudes (Figure D). Moreover, additional residues such as Ile886, Leu887, Ile890, Leu1017, and His1024 were depicted as significant for the binding of the crystalline ligand at the protein pocket.

Figure 8

Binding-free energy/residue decomposition illustrating the protein residue contribution at ligand–protein complex ΔGTotalbinding calculation. (A) FBZ; (B) MBZ; (C) ABZ; and (D) GIG.

The nature of the top energy-contributing residues is mostly hydrophobic which further emphasizes the predominance of the van der Waals potentials for binding the investigated ligands deep into the target binding site. However, the strong polar contacts with the hinge Cys917 provided a significant role in Coulomb’s electrostatic potentials for enforced ligand–protein biding. On the contrary, several other VEGFR-2 pocket-lining amino acids depicted positive energy contribution values with all ligands, which inferred repulsion forces and unfavored ligand-VEGFR-2 affinity. Among these energy-positive contributing residues are the ionizable Lys866, Glu883, Glu915, Asp1044, and Arg1049 amino acids. The furnished residue-wise energy contributions were consistent with ΔRMSF results, showing significant immobility and rigidity for all high energy-negative contributing residues while depicting higher fluctuations for the energy-positive contributors inferring significant repulsion and instability. Notably, the number of the energy-positive contributing residues was less in the case of GIG as compared to any of the anthelmintic agents where Lys866 and Arg1049 lacked unfavored energy contributions. Based on the previous conformational analysis, three of the repulsive residues (Lys866, Glu883, and Asp1044) showed close proximity toward the terminal aromatic substitutions of both FBZ and MBZ, the terminal propyl in ABZ, or the central aromatic ring of GIG ligands. That is why the substitution of the anthelmintics’ terminal structure with polar functionality (e.g., carbamide or guanidine moiety as close analogue) was suggested to overcompensate the repulsive-derived binding energy contributions through providing preferential polar binding interactions with Glu883 and/or Asp1044. Concerning the other side of the simulated ligands, the hinge Glu915 residues impose repulsive force against the aromatic phenyl ring within the core benzimidazole scaffold of all simulated compounds.

Certain residue-wise energy contributions were characteristic for the FBZ and GIG simulated systems. Unlike the other investigated anthelmintic agents, FBZ lacked the energy-positive contribution by Lys866, while expressed nonprecedential repulsive forces by Arg1049 as well as attractive favored binding interactions with Val865. The latter was reasoned because FBZ adopted a second different conformation following 50 ns simulation where its terminal aromatic tail became at close proximity toward Arg1049 rather than Lys866. Moving toward GIG, the ATP-binding specific residue, Glu883, imposed much higher positive energy contributions against the ligand-binding as compared to anthelmintic agents. Based on the conformational analysis, the Glu883 negatively charged residues were at additional closeness toward the fluorinated ring substitution. On the contrary, the unfavored positive-energy contribution by Asp1044 was much reduced in the case of GIG because being best oriented toward the polar carbamide rather than the hydrophobic parts of the ligand. Based on the above differential residue-wise contribution, it became highly reasoned that the fluorinated aromatic carbamide terminal substitution could impose a double-bladed impact upon the preferential anchoring of GIG within the VEGFR-2 pocket, which redeems further structure-based modification and optimization.

In brief, the MD findings were in great concordance with docking ordering values for the examined anthelmintic compounds and co-crystallized ligand. The free binding energy furnished through MM/PBSA was GIG > FBZ > MBZ > ABZ being at comparable ordering as obtained by the docking simulation findings. This may be reasoned for the strong polar interactions mediated by the ligand’s benzimidazole carbamoyl scaffold at the VEGFR-2 hinge region, while only the ligand’s terminal aromatic functionality would provide the differential ligand-target binding interactions. The same findings were consistent through protein ΔRMSF and ligand-target hydrogen bond analysis of the MD simulation trajectories. Nevertheless, the MD simulation provided one of the most interesting findings that have been highlighted, where FBZ managed to attain a second orientation/conformation beyond the 50 ns MD run. FBZ showed a significant conformational shift for its terminal aromatic moiety, yet kept the rest of the compound without relevant conformational/position changes (Figure A). Despite such conformational shift, FBZ managed to attain stabilized orientation within the VEGFR-2 binding site and top calculated free binding energy values across the investigated anthelmintic agents. This was suggested for a more favored position attained by the aromatic group where it formed an ideal T-shaped π–π stacking with Phe1045 overcompensating any depicted ring-residue unfavored electrostatic interactions. The latter observation further confirms the fact that the ligand’s terminal aromatic functionality would provide the differential ligand-target binding interactions. This was obvious with the superior free-energy of binding for co-crystallized ligand (GIG) which was correlated to the significant role of the polar/hydrophobic combined functionalities of its terminal fluorinated aromatic carbamide tail substitution for anchoring within the VEGFR-2 binding site.

Characterization of the Prepared MBZ-Loaded MMs

Obviously, the solubilization of a hydrophobic drug can be enhanced by its structure in biological media, and thus, the drug bioavailability and stability could be improved. The nanoscopic size of micelles (<100 nm) lets them be readily carried to the systemic circulation. Besides, MM possesses drug-loading efficiency and higher stability compared to micelles with single components.[37] Incorporation of vitamin-E in nanocarriers has been proved to improve the hydrophobicity of the drug delivery system, enhance the solubility of the loaded poorly soluble drugs, increase the biocompatibility of the polymeric drug carriers, and improve the anticancer potential of the anticancer agents by reversing the cellular drug resistance via simultaneous administration. In addition to being a powerful antioxidant, vitamin E demonstrated its anticancer potential by inducing apoptosis in various cancer cell lines.[59]

Particle Size, PDI, and Zeta Potential

The particle size of MBZ-loaded MMs was found to be 110.8 ± 1.1 nm, as shown in Figure . The polydispersity index (PDI) is a measure of the distribution homogeneity of molecular size in a given sample. Low PDI indicates that the BBV dispersion is homogeneous in nature.[60] The PDI was 0.17 ± 0.014, indicating good size homogeneity of the micellar solution.

Figure 9

(A) Size and size distribution of the micelles, (B) zeta potential of the MBZ-loaded MMs.

MMs with ZP > 30 mV were considered electrophoretically stable. MBZ-loaded MMs showed a negative ZP value of −22.8 ± 1.3 mV, and the negative charge is due to the presence of sodium deoxycholate (SDC) in the micelles.[61]

Transmission Electron Microscopy

To ensure the morphological structure of the prepared MBZ-loaded MMs, they were inspected using negative stain transmission electron microscopy (TEM), as shown in Figure . The shape of the prepared MMs was almost spherical and uniform in shape.

Figure 10

TEM micrograph of MBZ-loaded MMs. Scale bar in 200 (A) and 100 nm (B).

Entrapment Efficiency and Drug Loading

The entrapment efficiency of MBZ in the prepared MMs was 25 ± 2%. The drug loading capacity expressed as % of MBZ from the weight of the micelles forming materials was 1%.

In Vitro Dissolution of MBZ and MBZ-Loaded MMs

The cumulative percent of MBZ released, as a function of time, from prepared lyophilized MMs in comparison to plain drug powder is illustrated in Figure . Dissolution testing was performed in 0.1 N HCl (USP dissolution media of MBZ). Dissolution is well known as an important indicator of absorption and eventual bioavailability.[62]

Figure 11

Dissolution profiles of MBZ drug and MBZ-loaded MMs in 0.1 N HCl medium.

Interestingly, lyophilized MBZ-loaded MMs showed almost 100% (97 ± 0.02) release of MBZ after 2 h, while MBZ plain powder showed only 34 ± 1.3%. The preparation of MBZ in the form of MMs considerably improved its dissolution ability; consequently, this will increase the absorption of MBZ from GIT and provide a greater benefit in the biological system.

In Vitro Cytotoxic Activity

Based on the previously discussed in silico studies, three anthelmintic drugs, namely, FBZ (1), MBZ (2), and ABZ (3), together with MBZ-loaded MMs, and plain MMs were tested for their cytotoxic activity against three cancer cell lines, including liver cancer (HUH7) cell line, lung cancer (A549) cell line, and breast cancer (MCF7) cell lines. Different concentrations of drugs were prepared in dimethyl sulfoxide (DMSO) ranging from (1, 10, 100, and 1000 μg/mL), and incubated with the cancer cells for 48 h, the cell viability was tested using 3-[4,5-methylthiazol-2-yl]-2,5-diphenyl-tetrazolium bromide (MTT). Tested drugs showed potent cytotoxic activity against the liver cancer cell line (Figure ), while MBZ-loaded MMs was the best cytotoxic agent, with IC50 (6.54 μg/mL), which was even better than both plain MMs and MBZ free form (IC50 75.19 and 32.42 μg/mL, respectively).

Figure 12

Cytotoxic effect of the tested compounds against different cancer cell lines. (A) HUH7% cell viability upon treatment with a series of tested compounds, (B) A549% cell viability upon treatment with a series of tested compounds. (C) MCF7% cell viability upon treatment with a series of tested compounds. Concentrations were used starting from 1 to 1000 μM for 48 h, and the cytotoxicity effect was detected by the MTT assay (n = 3).

Testing anticancer activity against lung cancer showed that FBZ, MBZ, and ABZ had promising antilung cancer activity with IC50 (5.74, 22.45, and 32.98 μg/mL, respectively). Free MBZ showed slightly higher cytotoxic activity compared to MBZ-loaded MMs (31.93 μg/mL).

From another point, free FBZ, MBZ, and ABZ showed low cytotoxic activity against breast cancer cell line (MCF7), with IC50 values (216.38, 948.46, and 854.33 μg/mL, for FBZ, MBZ, and ABZ respectively). These results prove the high affinity of chosen drugs for VEGFR-2, as VEGFR-2 was noted to be poorly expressed in the MCF-7 cells both in vitro and in vivo.[63] Interestingly, MBZ-loaded MMs showed higher anticancer activity against MCF7 with IC50 (26.09 mg/mL), compared to plain MMs (240.99 μg/mL) and free MBZ (948.46 μg/mL) (Table ).

Table 5

IC50 of the Tested Compounds on Different Cell Lines

	IC50 (μg/mL)
compound	HUH7	A549	MCF7
FBZ (1)	28.23729	5.74467	216.38202
MBZ (2)	32.42383	22.45513	948.46979
ABZ (3)	21.1508	32.98596	854.33276
MBZ-loaded MMs	6.54465	31.9354	26.09982
plain MMs	75.19966	472.62847	240.99073
5FU	14.8919	19.6689	ND
Taxol	6.68	37.72706	25.61472
Dox.	4.95317	6.35791	11.40258

However, the enhanced cytotoxic activity of plain MMs incorporating vitamin E only or MBZ-loaded MMs is attributed to higher endocytosis by the cells and lower efflux of MBZ by the P-group pumps, allowing the MBZ to remain inside the cells for a longer duration and resulting in higher cytotoxicity (26.1 μg/mL).[59,64,65] Furthermore, the enhanced cytotoxicity of MMs may be due to the presence of F127 and vitamin E, which are known to be inhibitors of the P-group and reduce the efflux of drugs.[59,66] Vitamin E also shows both in vitro and in vivo anticancer activities against various cancer cells due to its apoptosis-inducing properties.[67]

Cell-Based VEGFR-2 Assay

To further validate the previously investigated in silico study, in vitro quantification for VEGFR-2 concentration in treated HUH7 cells with the drug under investigation has been conducted based on enzyme-linked immunosorbent assay. The results disclosed that all tested drugs significantly (p < 0.001) reduce the concentration of VEGFR2 from (3255 ± 312 pg/mL) for control to (1787 ± 37.3 pg/mL) for FBZ, (2041 ± 30.1 pg/mL) for MBZ, and (1213 ± 21.5 pg/mL) for ABZ, while the lowest inhibition for VEGFR-2 was achieved in MBZ-loaded MMs (860.8 ± 312 pg/mL), which was even much better than the reference drug, sorafenib (1073 ± 41.1) (Figure ) (Table ).

Figure 13

Schematic diagram illustrating the concentration of VEGFR-2 in HUH7 cells treated with the IC50 of tested drugs after 48 h using ELISA technique (n = 3), *p < 0.05, ***p < 0.001 compared to control.

Table 6

Concentration of VEGFR-2 in HUH7 Cells Treated with the IC50 of Tested Drugs after 48 h

compound	VEGFR-2 (pg/mL)
control	3255 ± 181
FBZ (1)	1787 ± 64.6
MBZ (2)	2041 ± 52.1
ABZ (3)	1213 ± 37.2
MBZ-loaded MMs	860.8 ± 57.1
sorafenib	1073 ± 41.1

SAR Studies

Therefore, based on the stabilities and binding scores of benzimidazole anthelmintic drugs to VEGFR-2, we could put eyes on the SAR that showed very promising and interesting activities against VEGFR-2, and hence the best expected antineoplastic activity. Moreover, concerning the anticancer activity of benzimidazole derivatives, the best affinity of anthelmintic drugs toward VEGFR-2 was attained when benzimidazole scaffold substituted at position 6 with phenyl sulfanyl group, FBZ (1), with benzoyl group, MBZ (2), or with propylthio group, ABZ (3). Furthermore, the studied SAR revealed that better activity could be accomplished when benzimidazole scaffold is substituted at position 2 with carbamate than any other group. Besides, it was noticed that substitution of benzimidazole scaffold at position 6 showed better activity against VEGFR-2 than unsubstituted ones (Figure ).

Figure 14

SAR study of the tested benzimidazole anthelmintics as VEGFR-2 antagonists.

Conclusions

Thirteen benzimidazole anthelmintic drugs were subjected to molecular docking studies as promising anticancer agents targeting VEGFR-2. The tested drugs showed comparable binding mode toward the ATP binding pocket of the VEGFR-2 receptor as compared to the potent co-crystallized benzimidazole urea inhibitor especially for FBZ (1), MBZ (2), and ABZ (3). Moreover, MD simulations revealed the great stability and binding affinity of the aforementioned investigated anthelmintic agents toward the VEGFR-2 active site. Achieving comparable free-binding energy and thermodynamic behavior along the 200 ns simulation time further confirms the suitability of these benzimidazole-based anthelmintic agents for repurposing approach targeting cancer’s angiogenesis pathway. Also, the enhanced release of MBZ from the MMs in USP dissolution media of MBZ ensured better oral availability after administration. Thus, allowing the administration of lower doses of MBZ consequently decreasing its side effects. Hence, these three promising drugs besides the MBZ-loaded MMs were subjected to further in vitro screening using three different cancer cell lines, including HUH7, A549, and MCF7 cell lines to confirm their anticancer activities. Accordingly, IC50 values for FBZ were 28.23, 5.74, and 216.38, for MBZ were 32.42, 22.45, and 948.46, for ABZ were 21.15, 32.98, and 854.33, and for MBZ-loaded MMs showed better promising in vitro cytotoxicity values of 6.54, 31.93, and 26.09 concerning HUH7, A549, and MCF7 cell lines, respectively. The MBZ-loaded MMs showed increased in vitro cytotoxicity in MCF7, HUH7, and neutral effect on A549 compared to free FBZ, MBZ, and ABZ. In conclusion, MMs could be used for enhancing the bioavailability and anticancer activity of MBZ. Finally, an in vitro quantification for VEGFR-2 concentration in treated HUH7 cells has been conducted based on an enzyme-linked immunosorbent assay. All tested drugs significantly (p < 0.001) reduced the concentration of VEGFR-2 from (3255 ± 312 pg/mL) for control to (1787 ± 37.3 pg/mL), (2041 ± 30.1 pg/mL), and (1213 ± 21.5 pg/mL) for FBZ, MBZ, and ABZ, respectively, while the lowest inhibition was achieved in MBZ-loaded MMs (860.8 ± 312 pg/mL), which was even much better than the reference drug sorafenib (1073 ± 41.1). Hence, MM nanoformulations could be used for the delivery of MBZ in cancer treatment. Additionally, the investigated benzimidazole anthelmintics can be treated as lead compounds for further structural modifications. That was accomplished after shedding light on the SARs improving their anticancer activity.

Materials and Methods

Materials

FBZ, MBZ, ABZ, and vitamin E were kindly gifted from Adwia, Pharaonic, and EVA pharmaceutical companies, Egypt. Pluronic F127 was procured from Sigma-Aldrich. SDC was purchased from BASF Co. (Florham Park, New Jersey, USA). Acetonitrile HPLC grade and ethanol 96% HPLC grade were purchased from CARLO ERBA Reagents (France). Ammonium acetate and hydrochloric acid were obtained from Merck (Darmstadt, Germany). Deionized (DI) water from the Ultrapure (type 1) water system (Direct-Q3 UV) was used for the preparation of all buffer and water-based solutions.

Docking Studies

Molecular docking studies of 13 benzimidazole anthelmintic drugs (1–13) at the VEGFR-2 receptor (PDB: 2OH4)[32] using MOE 2019.0102 drug design software[68,69] were performed to evaluate their predicted affinity as potent VEGFR-2 antagonists compared to the co-crystallized benzimidazole urea inhibitor (GIG, 14).

Preparation of the Tested Benzimidazole Anthelmintic Drugs

The chemical structures of the tested benzimidazole anthelmintic drugs were downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). They were converted into their 3D forms and checked for their chemical structures and the formal charges on atoms to be prepared for docking.[70−74] Then, energy minimization and automatic calculation of the partial charges were also done as previously described.[75−79] Finally, they were imported with the isolated co-crystallized benzimidazole urea inhibitor in the same database and saved as an MDB file before docking calculations with the target VEGFR-2 receptor.

Target VEGFR-2 Receptor Optimization

The X-ray structure of the target VEGFR-2 receptor complex was downloaded from the Protein Data Bank (http://www.rcsb.org/, PDB: 2OH4, resolution of 2.05 Å).[32] It was prepared for docking studies as follows: adding hydrogen atoms with their standard 3D geometry, and checking for any errors in the atom’s connection and the type through automatic correction as described earlier.[80−83] The applied force field was Amber10: EHT, and the co-crystallized ligand site was selected to be the docking site to show the different interactions with it in the complex structure.[84−86]

Docking of the Tested Drugs to the Target VEGFR-2 Receptor Active Site

The prepared database containing both the 13 benzimidazole anthelmintics (1–13) together with the benzimidazole urea co-crystallized inhibitor (GIG, 14) was docked using the MOE 2019 suite.[68,87] The applied methodology was as follows:[71,88−90] the docking was initiated as a general process after loading the file of the target protein active site. The docking site was selected to be the ligand site, triangle matcher was chosen as the placement methodology, and London dG was selected as the scoring methodology.[91] The refinement methodology was selected as a rigid receptor and the scoring methodology was GBVI/WSA dG for the selection of the best poses.[11,86] The MDB file was generally docked automatically. After the end of the docking process, the obtained poses were carefully studied, and the best ones for each having the best interactions and scores with the protein pocket were selected.

Best docked models of promising anthelmintic agent leads, as well as reference inhibitor, in complex with VEGFR-2, represented the starting coordinates for 200 ns explicit MD runs via GROMACS-2019 package.[92−94] The automatic CHARMM-General ForceField program (https://cgenff.umaryland.edu/) was used for ligand parameterization and topology file generation.[95] The CHARMM36m forcefield was considered suitable for simulated proteins.[33,89,96] Models were solely solvated in TIP3P cubic 3D-box at periodic boundary conditions with 10 Å marginal distances between protein and 3D-box sides.[97] The VEGFR-2 amino acids were set at their respective standard ionization within physiological pH 7.0. The entire system net charge was neutralized via sufficient chloride and sodium ions introduced using the Monte-Carlo ion-placement method.[98]

Following construction, the systems were minimized through 5 ps under the steepest descent algorithm[99] for optimizing the system geometry.[80] Under a constant number of particles, volume, and temperature (NVT) ensembles (303.15 K), the minimized system was equilibrated for 100 ps under the Berendsen-temperature coupling approach.[100] A second equilibration was done at a constant number of particles, pressure, and temperature (NPT) ensembles (303.15 K/1 atm pressure) under the regulation of the Parrinello–Rahman barostat approach.[101] A 1000 KJ/mol.nm2 force constant was used for all heavy atoms restrainment and original protein folding maintenance throughout the minimization/equilibration stages. Finally, 200 ns MD runs were produced under NPT ensembles using particle mesh Ewald algorithm for long-range electrostatic interaction estimation.[102] Linear constraint LINCS method allowed modeling of all covalent bond lengths (including hydrogen atoms) at 2 fs integration time step.[103] The van der Waals and Coulomb nonbonded potentials were truncated at 10 Å via Verlet cutoff schemes.[104]

Analysis tools, such as rmsd, RMSF, Rg, and SASA were estimated using GROMACS analysis scripts. The ΔRMSF was determined for each simulated liganded/holo protein relative to VEGFR-2 unliganded/apo state (PDB ID: 1VR2; ΔRMSF is RMSFapo–holo). The apo VEGFR-2 was subjected to similar preparation, minimization, equilibration, and 200 ns production stages, yet without ligand preparations. Visual Molecular Dynamics V-1.9.3 (VMD; University of Illinois, Urbana-Champaign, USA) was used for ligand–protein intermolecular hydrogen-bond across the whole MD timeframes. Hydrogen bond (acceptor···H-donor) cutoff values for angles and distances were established at 20° and 3.0 Å, respectively.[48,105] Using GROMACS “g_mmpbsa” module, the ligand–protein free binding energies were estimated through MM/PBSA energy calculation applied on representative frames of whole MD run (200 ns). This calculation approach provided more insights concerning ligand–protein magnitude affinity, nature of the ligand–protein binding, and protein residue-associated contributions at the free binding energy calculations.[58] Schrödinger-Pymol V.2.0.6 was used for representing the ligand-VEGFR-2 within the depicted conformational analysis at specific frames.[106,107]

MMs Loaded Mebendazole Preparation

MBZ-loaded MMs were prepared by using thin-film hydration followed by the sonication method.[62,108] Briefly, 10 mM concentration of F127 and SDC at molar ratio (7:3) with vitamin E (10% w/w) and MBZ (5% w/w) were dissolved in ethanol. The solvent was removed under low pressure, in a rotary vacuum evaporator set at a temperature of 40 °C. Then, 10 mL of DI water was added to the thin film and sonicated in a water bath sonicator for 10 min at room temperature. The resultant micellar solution was centrifuged at 10,000 rpm for 10 min (Sigma 3K30, Germany) to remove the unentrapped MBZ. The mixed micellar solution was kept at −80 °C for 2 h and then lyophilized (Alpha 2–4, CHRIST, Osterode am Harz, Germany) for 18 h. Similarly, plain MMs (without MBZ) were also prepared.

Characterization of the Prepared MBZ-Loaded MMs

Determination of Particle Size and Zeta Potential

The average size of the micelles and their PDI (size distribution) were determined by the dynamic light scattering method (Malvern Instruments, Malvern, UK). Cell temperature was kept at 25 °C with a detection angle of 90°. All measurements were performed in triplicate after appropriate reconstitution and dilution using DI water. The zeta potential of the prepared MMs was measured using the same instrument. All measurements were carried out in triplicate and the data are represented as mean ± SD.

Transmission Electron Microscopy

The morphology of the lyophilized MBZ-loaded MMs was observed by TEM (JEM-2100, JEOL, Tokyo, Japan). The lyophilized powder was reconstituted with DI water and shook for the formation of MMs. A drop (10 μL) of the prepared micellar solution was submitted on a carbon-coated copper grid and left to dry until a thin film is formed. This film was negatively stained with phosphotungstic acid (2% w/v) aqueous solution for 5 min and then air-dried at room temperature for 10 min before viewing and photography by TEM.

Drug Loading and Encapsulation Efficiency

10 mg of lyophilized MBZ-loaded MMs was accurately weighed and 5 mL of ethanol was added and sonicated for 2–3 min, and the samples were measured by HPLC using an Agilent 1100 HPLC (Agilent Technologies, Santa Clara, CA) instrument. Chromatographic separation of the samples was attained by a μBondapak, C18, 150 × 4.6 mm, 5 μm column (Waters, USA). The mobile phase was composed of acetonitrile 25% and buffer 75% (7.5 g of ammonium acetate in 1000 mL of DI water). The mobile phase was set at a flow rate of 1.2 mL/min with UV detection at 250 nm. The injection volume was 20 μL.[109] A calibration curve in ethanol was constructed in the concentration range of 1.6–240 μg/mL (linearity R2 = 0.9988). Encapsulation efficiency (EE) and drug loading (DL) was calculated according to the following equations

In Vitro Drug Release Study

The dissolution test was performed by adding 5 mg of MBZ drug and an amount of lyophilized MMs powder equivalent to 5 mg of MBZ was introduced into the USP II dissolution apparatus (Hanson Research, SR8PLUS, USA) filled with 100 mL of 0.1 N HCl in minivessels, maintained at 37 ± 0.5 °C, and agitated with a paddle rotating at a speed of 75 rpm.[109] Five-milliliter samples were withdrawn periodically at predetermined time intervals of 5, 15, 30, 60, and 120 min and replaced instantly by equal amounts of dissolution medium to maintain the initial volume. The dissolved amounts of MBZ at each interval were determined by HPLC as previously mentioned. A calibration curve in 0.1 N HCl was constructed for the concentration range of 1.6–20 μg/mL (linearity R2 = 0.9906).[109]

In Vitro Cytotoxic Assay

Cell Culture and Maintenance

Human hepatoma (HUH7), lung small cell adenocarcinoma (A549), and breast cancer (MCF7) cell lines were propagated in Dulbecco’s modified Eagle’s medium. All mediums were supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin antibiotics (Seralab, UK). For growth maintenance, the cells were incubated at 37 °C in 5% humidified CO2.

Evaluation of Cell Proliferation by MTT Assay

After treatment with variable concentrations of the compounds, the HUH7, A549, and MCF7 viable cell percentages were evaluated by the MTT assay as reported previously,[102,110−113] with a small modification. In brief, the cell count and viability were evaluated by the trypan blue dye-based method, and then, cancer cells (1 × 104 cells/well) were seeded using a 96-well plate. The cancer cells were then kept overnight for attachment. On the next day, the medium was exchanged completely with a fresh one, and then different concentrations of the synthesized compounds (0, 0.1, 1, 10, and 100 μM) were examined on each cell line. Thereafter, cells were permitted to grow for 24 h, and then, 10 μL of the MTT (5 mg/mL) was added to each well for 4 h before completion of the incubation period. After the incubation accomplishment, the formazan crystals were dissolved by adding 100 μL of DMSO to each well and left for 20 min, and then, the 96 well plates were waved for 5 min to assure the dye homogeneity in the solution. After the reaction occurs, a BioTek microplate reader was used to measure the color development at 490 nm. Based on the MTT results attained, we have chosen the most outstanding cytotoxic compound to investigate its mode of action.

VEGFR-2 Cell-Based Quantification

HUH7 cells were treated with the IC50 values of the investigated drugs. After 48 h, cell extraction buffer was used to lyse treated and nontreated cells, and then, the standard diluent buffer was used to dilute the lysate formed over the range of the assay and then measured for vascular endothelial growth factor (ab213476—Human VEGF Receptor 2 SimpleStep ELISA Kit). All procedures for standard ELISA techniques were done according to the methodology described earlier.[114]