Computational Investigation of 1, 3, 4 Oxadiazole Derivatives as Lead Inhibitors of VEGFR 2 in Comparison with EGFR: Density Functional Theory, Molecular Docking and Molecular Dynamics Simulation Studies

Abstract

Vascular endothelial growth factor (VEGF) is an angiogenic factor involved in tumor growth and metastasis. Gremlin has been proposed as a novel therapeutic pathway for the treatment of renal inflammatory diseases, acting via VEGFR 2 receptor. To date, most FDA-approved tyrosine kinase (TK) inhibitors have been reported as dual inhibitors of EGFR and VEGFR 2. The aim of the present study was to find the potent and selective inhibitor of VEGFR 2 specifically for the treatment of renal cancer. Fourteen previously identified anti-inflammatory compounds i.e., 1, 3, 4 oxadiazoles derivatives by our own group were selected for their anti-cancer potential, targeting the tyrosine kinase (TK) domain of VEGFR2 and EGFR. A detailed virtual screening-based study was designed viz density functional theory (DFT) study to find the compounds’ stability and reactivity, molecular docking for estimating binding affinity, SeeSAR analysis and molecular dynamic simulations to confirm protein ligand complex stability and ADMET properties to find the pharmacokinetic profile of all compounds. The DFT results suggested that among all the derivatives, the 7g, 7j, and 7l were chemically reactive and stable derivatives. The optimized structures obtained from the DFTs were further selected for molecular docking, and the results suggested that 7g, 7j and 7l derivatives as the best inhibitors of VEGFR 2 with binding energy values −46.32, −48.89 and −45.01 kJ/mol. The Estimated inhibition constant (IC₅₀) of hit compound 7j (0.009 µM) and simulation studies of its complexes confirms its high potency and best inhibitor of VEGFR2. All the derivatives were also docked with EGFR, where they showed weak binding energies and poor interactions, important compound 7g, 7j and 7i exhibited binding energy of −31.01, −33.23 and −34.19 kJ/mol respectively. Furthermore, the anticancer potential of the derivatives was confirmed by cell viability (MTT) assay using breast cancer and cervical cancer cell lines. At the end, the results of ADMET studies confirmed these derivatives as drug like candidates. Conclusively, the current study suggested substituted oxadiazoles as the potential anticancer compounds which exhibited more selectivity towards VEGFR2 in comparison to EGFR. Therefore, the identified lead molecules can be used for the synthesis of more potent derivatives of VEGFR2, along with extensive in vitro and in vivo experiments, that can be used to treat various cancers, especially renal cancers, and to prevent angiogenesis due to aberrant expression of VEGFR2.

1. Introduction

Cancer is one of the leading causes of mortality worldwide according to the World Health Organization (WHO) [1]. Among various cancer types, renal cancer is the seventh most common life-threatening disorder, causing 140,000 deaths annually [1,2]. The risk factors for renal cancer include smoking, alcohol intake, obesity, renal damage, and poor diet. Among various renal disorders, chronic kidney disease (CKD) is a fatal condition that affects between 5% and 7% of the world’s population and is a significant predictor of end-stage renal disease, cardiovascular morbidity, and ultimately death [3]. Regardless of the underlying cause, nearly all renal disorders result in permanent loss of kidney function due to gradual and irreversible nephron loss and diminished regeneration ability. The currently available chemotherapeutic agents only slow the progression of disease, necessitating the development of novel therapeutic drug approaches [4].

Tyrosine kinases are an important member of the enzyme family that catalyzes phosphorylation of tyrosine residues in protein by using ATP, which is a crucial step for normal cell development and homeostasis [5]. Potential anti-neoplastic therapeutic targets of tyrosine kinase that have been examined include epidermal growth factor receptor (EGF-R) and vascular endothelial growth factor receptor (VEGF-R). The epidermal growth factor receptor (EGFR) belongs to the tyrosine kinase family that is overexpressed in various cancers, such as breast cancer, prostate cancer and colorectal cancer, and plays an important role in cancer proliferation, growth and angiogenesis [6]. EGFR inhibitors have been developed in recent years that bind competitively to the active site of tyrosine kinase, thus blocking their uptake of ATP [7]. The gefitinib and erlotinib have been developed as EGFR inhibitors and many other drugs are undergoing clinical trials [8]. Another enzyme i.e., vascular endothelial growth factor receptors (VEGFRs) play an essential role by mediating angiogenesis, cell proliferation and blood vessel formation [9]. The VEGF receptors consist of three types, VEGFR1 (FLT1), VEGFR2 (KDR/FLK1) and VEGFR3 (FLT4). Among these, VEGFR2 is overexpressed in renal cancer, breast cancer, liver cancer and its activation by VEGFR signaling pathway is involved in tumor angiogenesis [10]. Recent advances in drug discovery have developed drugs that block VEGFR2 signaling pathways and thus inhibit angiogenesis [11], but these inhibitors because various adverse effects, such as nausea, skin rashes, and hypertension due to their non-selectivity [12]. The activation pathway of VEGFR and EGFR involves angiogenesis and cell proliferation [13], and is mentioned in Scheme 1 below.

Among various heterocyclic compounds, oxadiazoles are five-membered ring structures having one oxygen and two nitrogen atoms. Many oxadiazole derivatives have been reported to have anticancer potential, and three potent FDA-approved oxadiazole drugs with reported anti-cancerous properties are given below in Figure 1 [14,15].

In this study, we used structure-based virtual screening techniques to investigate the anticancer potential of 14 oxadiazole derivatives which have been reported earlier for their anti-inflammatory activity [16]. The purpose of the current study was to test the inhibitory potential of 1, 3, 4 oxadiazole derivatives against VEGFR2 for the treatment of renal cancer. The chemical structures of 14 oxadiazole derivatives are given in Figure 2. The IUPAC names of all derivatives are given in the Supplementary Materials (Table S1).

To measure the anti-cancer potential of our selected compounds (7a–n), a comprehensive multicomponent in silico study was designed. The study was comprised of density functional theory calculations and molecular docking studies, followed by molecular simulation studies. The pharmacokinetic profile of all the compounds was determined by calculating ADMET properties.

2. Materials and Methods

The whole procedure for the synthesis of selected derivatives, including the essential conditions, has been published in our earlier work [16]. Briefly, the study comprised of density functional theory studies with a focus on global reactivity descriptors i.e., hardness and softness of compounds, ionization potential, electron donating and accepting power. Molecular docking studies and molecular dynamic simulations studies were performed. The pharmacokinetic profile of all the compounds was then determined by ADMET software.

2.1. Density Functional Theory

The density functional theory is an effective method to obtain the optimized structures of compounds and for HOMO LUMO analysis. The Becke-3-Parameter-Lee-Yang-Parr (B3LYP) theory and SVP basis set was used for high accuracy in vibrational spectra both in gas and solvent phase calculations with the Gaussian 09 program [17]. The Gauss View 6 program [18] was used to obtain HOMO and LOMO structures, and to estimate the energy of frontier molecular orbitals.

2.2. Molecular Docking

2.2.1. Selection of Protein Targets

In order to execute the virtual screening process crystalline structure of VEGFR 2 and EGFR protein was downloaded from the RCSB protein data bank with the PDB ID: 3VHE for VEGFR 2 and 2GS6 for EGFR [19,20].

2.2.2. Softwares Required

The following softwares; MGLtools, binary files of Autodock4 and Autogrid4 [21], BIOVIA’s Discovery Studio Visualizer [22], SeeSAR [23], chemdraw ultra [24] and chemdraw 3D pro [25] were used for molecular docking.

2.2.3. Preparation of Protein

The downloaded protein from RCSB was further processed for executing autodock studies. All the heteroatoms along with co-crystal ligand and solvent molecules were removed from the protein molecule in BIOVIA’s discovery studio visualizer. The pure protein structure obtained was then prepared for docking with the help of autodock inbuilt tools by adding necessary polar hydrogen and Kollman charges to each atom and was saved in pdbqt format [26].

2.2.4. Preparation of Ligand and Molecular Docking

The structures of all the compounds were drawn in chemdraw ultra by using their IUPAC names and then the energy minimization for all the structures was carried out using chem 3D pro. The compound structures were saved in SDF format. These structures were then converted to autodock acceptable format i.e., pdbqt through openbabel GUI software. Structure based virtual screening was then performed using autodock4 software. The grid box dimensions were kept same as for co-crystal ligand to generate grid parameter file i.e., x: y: z = −28.96 Å, −4.140 Å, −14.515 Å for VEGFR 2 and grid box dimensions for EGFR was x: y: z =142.128 Å, 26.78 Å, 52.91 Å. For the generation of docking parameter file, Lamarckian genetic algorithm (LGA) and in house force field named Autodock4Zn were used. The no. of poses was set to 100 and population size was 300 to ensure maximum validity and reliability of the scoring function was estimated from extensive literature review [27]. The prepared ligand library was then docked into the active site of both proteins separately [21]. The docking of all the compounds was reconfirmed by SeeSAR Analysis [23].

Visualization

The BIOVIA discovery studio visualizer was used for the analysis of ligand and protein interactions. All the 2D and 3D conformations were generated by adding autodock result files and protein pdbqt files in BIOVIA’s Discovery Studio. The various bonding and non-bonding interactions of the ligand and active pocket were identified accordingly.

Validation

The docking protocol was validated on the basis of the RMSD value and by re-docking the co-crystal ligand into the active pocket of protein. Only poses with the RMSD value of docking and experimental value of ligand less than 2.0 Å were accepted [28].

2.2.5. Molecular Dynamics Simulations

Desmond from Schrödinger LLC suite, developed by Bowers and his team [29], was used to simulate protein ligand complex for 100 nanoseconds. In the case of protein and ligand complexes, docking investigations were the first step in preparing them for molecular dynamics modeling. The binding status of a ligand into the active pocket of protein was calculated using molecular docking studies under static conditions. MD simulations calculate the motion of atoms over the period of time by applying Newton’s classical equation of motion into the molecule’s binding position in a protein’s active site [30]. By employing MD simulations, it was possible to anticipate the ligand binding status in extreme physiological conditions [31,32].

The Protein Preparation Wizard or Maestro software was used to preprocess the protein–ligand complex, which includes optimization and energy minimization of the protein–ligand complex structure. All the files were created with the help of the System Builder tool. Because of its simplicity, TIP3P was chosen as a solvent model with an orthorhombic box (Transferable Intermolecular Interaction Potential 3 Points). The OPLS-AA 2005 force field was utilized to execute the MD simulation [33]. Counter ions were added to the models to make them behave in a neutral manner i.e., 0.15 M sodium chloride (NaCl) was added to mimic physiological circumstances. For the whole simulation, the NPT ensemble with a temperature of 300 K and a pressure of 1 atm was employed. Prior to the simulation, the models were made more flexible. After every 100 ps, trajectories were saved for inspection, and the simulation’s stability was confirmed by calculating the root mean square deviation (RMSD) of the protein and ligand over time.

2.2.6. Cell Viability Assay

In order to estimate the anti-cancer potential, all the derivatives were tested against a human cervical cancer cell line (HeLa) and human breast cancer cell line (MCF-7), by using the already reported method of Mosmann (1983) and Nik and Otto (1990) [34,35] as discussed in our previously published article [36]. The experiment was carried out in 96-well flat-bottom plates with 90 μL of medium and 10 × 10⁴ cells seeded into each well. 100 μL of test chemical solution was added to each well, and the plate was incubated for 24 h at 37 °C and 5% CO₂. The positive and negative control wells were seeded with 10 μL of standard drug (cisplatin) and 100 μL of cells medium (without chemical). After that, each well was pipetted with 10 μL of MTT reagent and incubated for 4 h at 37 °C. Then, 100 μL of 10% sodium dodecyl sulphate solution was added and the mixture was maintained at room temperature for 30 min with intermittent shaking. Finally, the optical density was computed. All tests were done in triplicate, and the findings were given as percent growth inhibition values, as previously described [36].

2.2.7. ADMET Properties

Predicting ADMET properties is a crucial step in the drug development process. The online web server ADMETlab 2.0 was used to calculate the physicochemical parameters, medicinal properties, absorption, distribution, metabolism and toxicity of derivatives. It is an accurate and reliable online platform for predicting ADMET properties [37].

3. Results and Discussion

3.1. Synthesis of 1, 3, 4-Oxadiazole Amide Derivatives

In our previously reported study [16], the target compounds 7a–n were synthesized in a series of steps. First, 5-chlorophenyl-1, 3-dimethyl-4-oxadiazol-2-thiol (3) and 2-bromo-N-[aryl/aralkyl]propionamide (5a–n) were produced separately. The following steps were used to produce compound 3: The compound 4 chlorobenzoic acid (a) was first refluxed in ethanol under the catalysis of concentrated H₂SO₄ to produce the corresponding ethyl ester (1), which was then reacted with hydrated hydrazine in methanol to produce carbohydrazide (2), which was then refluxed with carbon disulfide in the presence of ethanolic potassium hydroxide to yield the cyclized product 5 (4 chlorobenzoic acid (3). The other aryl/aralkyl amine precursors (5a–n) of the target compounds were synthesized by combining the aryl/aralkyl amines (4a–n) with 2 bromopropionyl bromide in basic conditions (pH 9–10). The required compounds (7a–n) were successfully synthesized in high yield by treating the precursor 3 with the electrophiles 2-bromo N [aryl/aralkyl] propionamide (5a–n). Infrared, nuclear magnetic resonance (NMR), and 13C NMR spectroscopy, as well as electron ionization mass spectrometry (EIMS) and high-resolution electron ionization mass spectrometry (HREIMS) [16] were used to determine the structures of substances (Figure 3).

3.2. Density Functional Theory Calculations (DFTs)

All the DFT Calculations were carried using Gaussian-09 via DFT35 with B3LYP36 Hybrid GGA functional theory and SVP basis set. The Table 1 shows optimized geometry parameters in gaseous and ethanol phase.

Results of DFT Studies: Global and Local Descriptors

The table of global reactivity descriptors reveals that, according to the ionization potential data, 7j and 7n require less energy to remove an electron from the ground state, whereas 7b required the highest energy indicating the highly reactive nature of 7j. However, the value of this attribute for 7j and 7b was 0.2300 and 0.24068 eV, as estimated with SVP basis set. The compound 7j has the highest value of LUMO and HOMO energy gap representing the highest stability of this derivative. The electronegativity value indicated that 7b has the maximum ability to attract electrons in both phases due to negligible effect of solvent. This compound would be energetically favorable for a nucleophilic attack. The findings of all parameters are significantly different for each of the compounds, as shown in Table 2.

The electron-accepting and electron-donating power of the various oxadiazoles as well as their tendency to donate or accept a small amount of charge required for chemical interactions, were determined. Chemical hardness of any compound is an important parameter to determine its reactivity potential and compounds with higher values of hardness are less prone to be reactive. Among 14 oxadiazoles, most of the derivatives exhibited almost similar hardness values and are within the acceptable range of reactivity. The values of all other parameters for all the compounds are given below in Table 3.

The optimized structure of top 3 potent compounds having Frontier molecular orbitals (HOMO & LUMO) along with energy gap is presented in Figure 4, while the other compounds are given in Supplementary File (Figure S3).

3.3. Molecular Docking

The crystalline structure of VEGFR 2 indicated that it is separated into two lobes by its folded form. The phosphotransfer catalysis takes place in an inter lobar gap between the two lobes. When the N-terminal lobe (residues 820–920) is folded, one helix and one twisted beta sheet structures are formed (aC). The b structure is made up of five antiparallel strands (b1–b5), three of which are severely coiled and curl over the other two strands (b4–b5), while the other two strands are straight. The larger C-terminal domain (residues 921–1168) has two antiparallel b strands (b7–b8) near the N-terminal b sheet at the top of the C-terminal domain. At the bottom of the C-terminal domain is the N-terminal b sheet. The C-terminal domain’s basic structure is made up of seven a helices (aD, aE, aE-F, aF, aG, aH, and aI). The glycine-rich nucleotide binding loop (residues 841–846), the catalytic loop (residues 1026–1033), and the activation loop (residues 1046–1075) are all functionally essential loops for the proper functioning of VEGFR 2.

The active pocket of VEGFR 2 consists of the following sequence of amino acids i.e., Asp814, Cys817, Leu840, Val848, Ala866, Lys868, Ala881, Leu882, Ser884, Glu885, Ile888, Leu889, Ile892, Val898, Val899, Val914, Val916, Glu917, Phe918, Cys919, Lys920, Phe921, Gly922, Asn923, Leu1019, Cys1024, Ile1025, His1026, Arg1027, Leu1035, Ile1044, Cys1045, Asp1046, Phe1047, Gly1048, and Leu1049 [38] and the active pocket of EGFR protein consists of the following sequence of amino acid Leu694, Gly695, Ser696, Phe699, Val702, Ala719, Ile720, Lys721, Leu723, Lys730, Ala731, Lys733, Glu734, Ile735, Asp737, Glu738, Val741, Met742, Cys751, Leu764, Ile765, Thr766, Arg812, Asp813, Arg817, Asn818, Leu820, Thr830, Asp831, Phe832, Gly833, Leu834, Ala835, Lys836, Tyr845, Ala847, Glu848, Gly849, Gly850, Ys851, Val852, Pro853, Arg865, and Tyr867. The Tyr6 is only amino acid from 2^nd chain that is involved in active pocket activation [39].

3.3.1. Structure Activity Relationship (SAR) of 1, 3, 4-Oxadiazole Amide Derivatives

The binding energies of all the derivatives were determined against VEGFR2 and EGFR. All the derivatives showed strong interactions within the active pocket of VEGFR2 and exhibited less binding score against EGFR. The differential analysis of molecular docking was based on the fact that any potent derivative with strong interactions has minimum binding energy [26]. The structure activity relationship was studied on the basis of binding interactions, binding scores and predicted inhibitory concentration values obtained during docking analysis (Table 4).

Among all the derivatives, the parent compound i.e., 7a exhibited high binding energy value and less estimated inhibitory constant value among all the derivatives. The compound 7j (−48.89 kJ/mol) which was found as the best inhibitor has methyl substitution at 3 and 5 position of the benzene ring. When the structure and activity of compound 7j was compared with other derivatives it was found that the substitution of methyl group at meta position is responsible for the inhibitory potential of this compound. This effect can be seen in case of compound 7g (−46.32 kJ/mol) which was found as the second most potent inhibitor of VEGFR2. In this compound the substitution of methyl was done at ortho and meta position i.e., 2nd and 5th position of phenyl ring. It can be suggested that the shifting of one methyl from 3rd to 2nd position resulted in the slight difference in inhibitory potential of compounds. Both compounds 7j and 7g exhibited equipotent Estimated inhibitory concentration values i.e., 0.0009 μM. The methyl group is electron donating group which imparts a positive mesomeric effect (+M) reasonable for inhibitory potential. The effect of substitution was compared between those derivatives in which ortho substitution was carried out i.e., derivative 7b, 7d, 7e and 7f. The compounds 7b and 7d have methyl substitution at ortho position which was found less favorable for the ligand-protein interaction as well as Estimated inhibitory values were found less potent. On the other hand, compound 7e and 7f exhibited better results. The detailed analysis suggested that the compound having substitution at ortho position alone (either at position 2 or position 4) or di substituted (2 and 4 position) were found less potent whereas, the compound with di-substitution at one ortho and one meta showed better result suggesting that in these derivatives. the substitution with meta group resulted in improved binding interaction as well as Estimated inhibitory constant value. An interesting behavior was observed when the substitution at ortho position was done with the bulkier group i.e., ethyl and ethoxy as in case of 7k and 7l (ethyl substitution); 7m and 7n (ethoxy substitution). All these compounds were found less potent as compared to derivatives with less electronegative CH₃ group.

The structure activity relationship of these derivatives was also studied against EGFR, and it was observed that these compounds showed less potent behavior against EGFR with higher binding energy and less inhibitory potential as depicted by their estimated inhibitory concentration values. However, the pattern of the substitution effect was same as observed in VEGFR2. The compounds having meta substitution were found more potent than para and other substitution.

3.3.2. Binding Interaction Studies

Binding Interactions of 1, 3, 4 Oxadiazoles with VEGFR2

In order to study the binding interactions of these derivatives, first of all the co-crystal ligand was re-docked within the active pocket of the selected protein (Figure 5). Here, compound 1-{2-fluoro-4-[(5-methyl-5H-pyrrolo [3, 2-d] pyrimidin-4-yl) oxy] phenyl}-3-[3-(trifluoromethyl) phenyl] urea was the co-crystal ligand for our target protein. The re-docking analysis suggested that this compound formed 3 covalent bonds by the three fluorine atom attached to the methyl phenyl ring with three amino acid residues i.e., His1024, Ile1044 and Cys1045, respectively. This ring was also involved in forming one sigma bond with Leu889 and another sigma bond was formed by 6 membered pyrimidine rings with Leu1035 amino acid residue. This derivative formed 5 conventional hydrogen bonds with different amino acid residues. The two amine groups of phenyl urea ring formed two hydrogen bonds with the same amino acid residue, Glu885, the carbonyl oxygen formed one hydrogen bond with Asp1046, the nitrogen atom of pyrimidine ring formed one hydrogen bond with Cys919 and one hydrogen bond was formed between fluorine atom and His1026 amino acid residue. The methyl group attached to the pyrole ring was involved in making one carbon hydrogen bond with Leu840. This Leu840 also formed one sigma bond with pyrole ring and one weak alkyl bond and one pi-sigma bond with pyrimidine ring. Besides these bonds, single pi-pi T shaped bond (phe1047) was formed with benzene ring present adjacent to pyrimidine ring. The other with various weak alkyl and pi-alkyl bonds were formed with Val898, Ile1044, His1026, Lys868, Val916 and Val848.

The most potent derivative i.e., 7j formed different types of bonding and non-bonding interactions with Phe1047, Val916, Ile1044, Val848, Leu1035, Ala866, Phe918 and Leu840. This derivative formed 2 pi-sigma bonds, 2 hydrogen bonds, 3 pi-alkyl interactions and one pi-anion interactions. Regarding 2 sigma bonds, one sigma bond was formed between xylene ring and Leu889 and one sigma bond was formed between 1-chloro-4-methylbenzene ring and Leu1035. This 1-chloro-4-methylbenzene was also involved in making three pi-alkyl interactions with Ala866, Cys919 and Phe918 amino acid residues. Among two hydrogen bonds, one was formed between N-methylformamide group and Glu885, and second hydrogen bond was formed by the oxygen with Asp1046. Moreover, among 3 pi-alkyl interactions, one pi-alkyl interaction was formed by oxadiazole ring with Val848. This oxadiazole ring was also involved in forming one pi-anion interaction with Cys1045. The amino acids involved in non-bonding interactions were: Val916, Val848 and Leu840 amino acid residues.

The second most potent derivative i.e., 7g formed different types of bonding and non-bonding interactions with Leu1035, Leu889, Cys1045, Glu885, Val899, Lys868, Asp1046, Ile892, Phe1047, Val916, Ile1044, Val848, Leu1035, Ala866, Phe918 and Leu840. This derivative formed 3 pi-sigma bonds, 2 hydrogen bonds, 8 pi-alkyl interactions, one pi-anion and one pi-cation interactions. The chlorobenzene ring formed one pi-sigma interaction with Leu889 and one pi-alkyl interaction with Ile892. The oxadiazole ring was involved in forming Hydrogen bond interaction with ASP1046, pi-anion interaction with Glu885, Lys868 and pi-alkyl interaction with Val899. Moreover, N-propylacetamide formed strong hydrogen bond with Lys721. The toluene ring showed three pi-alkyl interactions Ala866, Cys919 and Val848 with amino acid residues. The toluene ring was also involved in forming one pi-sigma interaction with Leu1035.

The ligand interaction of compound 7l with the active pocket rediues showed that the chlorobenzene ring formed pi-alkyl interactions with Ala866, Leu840, Cys919, Phe918 amino acids and also formed one pi-sigma interaction with Leu1035. Moreover, the oxadiazole ring was involved in forming pi-anion interaction with Cys1045 and one Pi-alkyl interaction with Val916. The toluene ring formed one pi-sigma interaction with Leu889 and two pi-alkyl interaction with His1026 and Leu1019 amino acid residues. The 3D and 2D binding modes of most potent compounds 7j, 7g and 7l are mentioned in Figure 6 and Figure 7 while 2D, 3D interactions of other compounds are given in Supplementary File (Figures S4 and S5).

Binding Interactions of 1, 3, 4 Oxadiazoles with EGFR

Following the same method as with VEGFR protein, the co-crystal ligand i.e., analogue of ATP (Thiophosphoric acid o-((adenosyl-phospho) phosphor)-s-acetamidyl-diester) of EGFR protein was re-docked within the active pocket (Figure 7). The various amino acid residues which were involved in bonding and non-bonding interactions are: Ala698, Phe699, Asp831, Asn818, Arg817, Glu738, Phe832, Met769 and Lys721. The analysis of docking results suggested that 11 conventional hydrogen bonds were formed between co-crystal ligand and different amino acid residues of active pocket. The two-hydroxyl group of pentose sugar formed four hydrogen bonds with Asp831, Asn818, Arg817 amino acid residues, one hydrogen bond was observed between oxygen atom and Lys721. One conventional hydrogen bond and one carbon hydrogen bond was formed between Ala698 and nitrogen atom of adenine ring. In the same way, Phosphate group adjacent to adenosine moiety formed two hydrogen bonds with Glu738 and Phe832, respectively. Moreover, acetamides group formed two hydrogen bonds with Met769 amino acid residue. The co-crystal ligand showed one pi-sigma bond by adenine ring with Phe799 and also exhibited two cationic bonds by electropositive Phosphorus atom with Asp831 and Glu738, amino acid residues.

According to the docking results 7j is the most potent compound of the series. The compound 7j formed different bonding and non-bonding interactions with following amino acids i.e., Ala719, Leu768, Met769, Val702, Leu820, Thr766, Lys721, Leu764, Glu738, Thr830 and Asp831. The most common bonding interactions included hydrogen bonding, pi-sigma bonds, weak cationic bonds and pi-alkyl bonds. The docking results of 7j suggested the presence of 3 hydrogen bonds, 1 pi-sigma bond, 2 cation bonds and 5 pi-alkyl bonds. The amino acid residues Thr766, Lys721 and Thr830 formed 3 hydrogen bonds with nitrogen atom of propanamide group, oxygen atom of oxadiazole ring and 1st nitrogen atom of oxadiazole ring. The pi-sigma bond was formed by the xylene ring and Leu820. The two cationic bonds were formed by the Glu738 and Asp831 with the oxadiazole ring and adjacent phenyl ring. In addition to these 5 weak pi-alkyl bonds were also formed by xylene ring and 3 amino acid residues i.e., Ala719, Leu768, and Met769.

The docking interaction of 2nd best compound 7g involved the following amino acids residues i.e., Thr766, Leu764, Lys721, Ala719, Leu694, Leu768, Gly772, Met769, Leu820, Thr830 and Asp831. The analysis of docking results revealed the presence of one conventional hydrogen bond, one carbon hydrogen bond, two pi-sigma bonds and seven pi-alkyl bonds. The carbon hydrogen bond was present between Thr766 and xylene ring while conventional hydrogen bond was formed by the Lys721 and oxygen atom of amide group. One pi-sigma bond was formed by amino acid Leu820 and oxadiazole ring and other pi-sigma bond Leu694 and phenyl ring. The weak pi-alkyl bonds were formed by 5 amino acid residues i.e., Leu764, Lys721, Ala719, Leu768 and Leu820.

The sequence of amino acid residues involved in the docking interaction of 7l are Asp831, Thr830, Leu820, Met769, Leu694, Ala719, Leu768, Gln767, Met742, Lys721 and Thr766. The bonding and non-bonding interactions of 7l included 2 conventional hydrogen bonds, 1 carbon hydrogen bond, 1 pi sigma bond, 4 pi-alkyl bonds and 1 sulphur bond. The only carbon hydrogen bond was formed by the Asp831 and ethyl benzene ring. In the same way one conventional hydrogen bond was formed by Thr830 and Lys721 with the hydrogen and oxygen atoms of propanamide group. The only pi-sigma bond was formed by Leu820 and Phenyl ring and the sulphur bond was formed by the sulphur atom of Met742 and oxadiazole ring. The pi-alkyl bonds were formed by Leu768, Leu694, Ala719 and Leu820 with the ethylbenzene ring and methyl group attached to phenyl ring. The 3D interactions of the most potent derivatives are shown in Figure 6 while 2D and 3D interactions of other derivatives are given in the Supplementary File (Figures S6 and S7).

3.4. SeeSAR Analysis

The SeeSAR study was used to compute the HYDE (Hydrogen bonding and Dehydration) score using BiosolveIT’s SeeSAR [23]. The computation of desolvation value, hydrophobic interaction, and hydrogen bonds established between the ligand and active pocket are all part of the HYDE analysis of a molecule. HYDE also aids in anticipating the specific region of the complex that experiences favorable and unfavorable ligand receptor binding and is represented by green and red coronas around the specific atom. The larger coronas gave higher HYDE value representing greater involvement in binding affinity. All the atoms (dark green sphere) that contributed favorably towards binding affinity and their individual HYDE values for the best molecule are highlighted for the most potent compounds. The role of atoms in each ligand must be determined to anticipate overall binding affinity or interactions with other molecules. It is evident from the results that specific region of each ligand contained few atoms that contribute positively toward the binding affinity. The HYDE value of all potent compounds is greater for VEGFR 2 than EGFR (Figure 8).

3.5. Molecular Dynamics Simulations

The protein-ligand complexes of VEGFR and EGFR with the most favorable docking pose of 7j were simulated in aqueous conditions for 100 ns using Desmond by Schrodinger. The structural behavior of both proteins and their complexes was determined by Root Mean Square Deviation (RMSD), Root Mean Square Fluctuation (RMSF), Solvent accessible surface area (SASA) and Radius of gyration (Rg) along with other MD simulation parameters. The RMSD graph provides the detail of any structural fluctuation occurred because of ligand-protein interaction. For the analysis of MD simulation results the RMSD plot of both proteins and their complexes was generated and compared to evaluate the stability of protein complex. Figure 8 shows the RMSD graphs of both proteins along with their complexes. The RMSD graph of VEGFR complex showed some initial fluctuation and then became stable after 30 ns at RMSD value of 1.8 while the graph of EGFR showed that the protein is relatively unstable, and its 7j-complex became stable after 50 ns with RMSD value of 2.5. The data clearly shows that VEGFR 2 and its complex was found more stable in aqueous medium as compared to EGFR which augment the docking results (Figure 9).

Further analysis of MD simulation was carried out by generating RMSF plots, and the results are presented in Figure 10. Any fluctuation in the amino acid residues of C and N terminal lobe in target protein can be identified by RMSF graphs. The average RMSF value for C_α chain with majority of amino acid residues of VEGFR was found to be less than 1 Å and for the residues of EGFR protein was below 1.2 Å suggesting higher stability of VEGFR protein. The very few residues of the protein VEGFR, on the other hand, experienced small variations, which could be related to the hanging position. The RMSF data further established the stability of VEGFR proteins and its complex in aqueous conditions. The RMSF value of both ligands exhibited some fluctuations which showed their dynamical shift at their binding site in respective proteins. (Further 2 replicates of 50 ns for each protein has been added in the Supplementary File, i.e., Figures S8 and S9).

The radius of gyration (Rg), solvent accessible surface area (SASA), MSD value and intra-molecular hydrogen bonding of both systems was calculated by analyzing MD simulation data. The mass-weighted RMS distance of atomic cluster from their common center of mass is known as Radius of gyration (Rg). In MD simulation studies Rg is an important factor in finding the compactness of protein [40]. Rg value of VEGFR 2 was found to be uniform throughout the simulation length, demonstrating their stability with an average value of 5.2 and for EGFR it was 4.0. The Rg of both proteins remained constant during the simulation, indicating that the systems did not experience any notable conformational changes [41]. SASA is another factor used in MD simulations to assess protein stability [42]. Throughout the simulation period, SASA of VEGFR 2 and EGFR were remained stable. The average SASA of VEGFR 2 and EGFR was 55Å and 160Å, respectively. Figure 11 summarizes different graphs for both proteins.

Moreover, the alpha-helices and beta-strands were considered as secondary structure elements during the simulation (SSE). The distribution of SSE across the protein structure by residue index is given in Figure 12. The figure below shows the SSE composition for each trajectory frame during the simulation, SSE assignment for each residue over time. The residue index of VEGFR ranges from 60–70, 140–155 and 230 to 300 belongs to beta strand. The % SSE graph of VEGFR complex shows the least fluctuation and is comparatively more stable than EGFR.

At the end, the protein interactions with the ligand were also observed during the simulation process. These developed linkages were identified by four types of protein-ligand interactions i.e., hydrogen bonds, hydrophobic interactions, ionic interactions, and water bridges. The ‘Simulation Interactions Diagram’ panel in Maestro was used to study the subtypes of each interaction type. Over the course of the trajectory, the stacked bar charts were standardized: for example, a value of 0.7 indicated that the specific interaction has been maintained for 70% of the simulation duration. Due to several interactions with the same subtype of ligand, values above 1.0 were observed which is also feasible.

Interestingly, the hydrogen bonds observed during the study exhibited a considerable influence on drug selectivity, metabolism, adsorption, and hydrogen-bonding properties that is important in drug designing [43,44]. Here, the four types of hydrogen bonds included: backbone acceptors, backbone donors, side-chain acceptors, and side-chain donors that exists between a protein and a ligand (Figure 13).

The majority of significant ligand–protein interactions detected by MD simulations are water bridges and hydrogen bonds, as shown in Figure 12. In terms of hydrogen bonding, ASP-1046 and CYS-919 are the most significant residues for 7j-VEGFR 2, while MET-769 is the most critical amino acid for 7j-EGFR complex. Over the course of the trajectory, the stacked bar charts were standardized: for example, a value of 1.0 indicates that the specific interaction was maintained for 100% of the simulation time. Many protein residues may have several interactions with the ligand of the same subtype, therefore values above 1.0 are also attainable. Figure 14 depicts individual ligand atom interactions with both protein residues. Only those interactions that occur in the chosen trajectory (100.0 ns) for more than 30.0 percent of the simulation duration are mentioned.

3.6. Cell Viability Assay

To support the in silico studies, the preliminary screening of the most potent derivatives was carried out using in vitro cell viability assay i.e., MTT assay (Table 5).

The human cervical cancer cells (HeLa) and MCF-7 cells were treated with 100 µM concentration for 24 h and 48 h. The concentrations were selected based on the predicted inhibitory values obtained during docking studies (as shown in correlation plot Figure 15).

A time dependent linear response of cell death was observed at single dose concentration. The compounds showed good results and caused % growth reduction but the derivatives 7c, 7g, 7h, 7j, 7l, 7k and 7m showed maximum cell death justifying the computational studies where these compounds were found to be the best inhibitors of VEGFR-2 and also EGFR. The single concentration cisplatin was used as positive control. The results were calculated by comparing with the total activity control (without inhibitor) i.e., un-treated cells. The % viability graph was generated via graph pad Prism software and is given in the Figure 16.

3.7. ADMET Properties

The ADMET (absorption, distribution, metabolism, excretion, and toxicity) properties of all the VEGFR 2 inhibitors were calculated to compute its appropriateness as a drug molecule. Other physicochemical properties explored are mentioned in our previous paper [45].

The physiochemical attributes of a clinically active substance have a substantial impact on its ability to become a therapeutically effective medicine. Compounds with MW < 500, nON < 10, nNH < 5, nRot < 10, and cLogP < 5 are considered orally bio-available and have a good ADMET profile [37]. The most powerful molecule, 7j, has an excellent physicochemical profile, according to the findings. The physicochemical properties of compounds are mentioned in Table 6.

According to ADMET results of absorption and distribution, the compound 7k has maximum volume of distribution. Other properties like blood–brain barrier, placenta barrier, human intestinal absorption, plasma protein binding, P-gp inhibitor and P-gp substrate indicated that 7a, 7j and 7g are pharmacokinetically the most suitable. The values of other compounds for all these properties are given in Table 7.

The metabolism and excretion profile of all the 1,3,4 oxadiazoles suggested that most of the potent derivatives are the best target for liver enzymes i.e., CYP1A2 and CYP2C19 while all other liver enzymes are best complexed with different compounds. The value of 7m, 7c and 7j is highest for CYP2C9, CYP2D6 and CYP3A4. The ADMET results suggested that the compound 7h is excreted slowly from the body and the small half-life of 7i is responsible for its shortest stay in the body (Table 8).

The medicinal properties of a compound explained the potential of a candidate molecule to become a drug. The data showed that all the 14 compounds are in agreement with Lipinski rule of 5 which included molecular weight, donors and acceptors of H bonds, and partitioning coefficients (log P). The synthetic accessibility score of all compounds was found below 3.0, indicating the ease of synthesizing them (Table 9).

The toxicity profile of any molecule is a paramount criterion for its selection regarding further studies. The drugs with high toxicity cannot be considered for human trials which is indicated by −sign while +sign indicates acceptable results. The results of our study showed that 7e and 7i are most toxic and 7a and 7n are least toxic among the given oxadiazoles derivatives while all the compounds are equally dangerous for eyes. The data showed that 7g, 7j and 7n are least toxic. All the values are given in Table 10.

As a result, we can conclude that our molecules, particularly the most potent derivatives 7g, 7j and 7l exhibit drug-like properties. Further research is required to investigate toxicity profiles of given compounds.

4. Conclusions

The findings of current study convincingly suggested that selected 1,3,4 oxadiazole derivatives are potent inhibitors of VEGFR 2 compared to EGFR and can be used in the treatment of cancers where selective VEGFR2 inhibitors are required i.e., renal cancer and to mitigate angiogenesis in other cancers. According to the results of the molecular docking investigations, density functional theory, cell viability assays and MD simulation, 7g, 7j and 7l are the most potent compounds and should be further explored for the development of new molecules.

Author Contributions

Conceptualization, S.A.E.; Data curation, M.S.B. and F.S.; Formal analysis, M.S.B., S.Z., A.T.A. and N.A. (Nojood Altwaijry); Methodology, A.T.A., F.S. and H.I.U.; Project administration, S.A.E., N.A. (Naveed Akhtar) and T.A.W.; Resources, S.Z., N.A. (Naveed Akhtar), T.A.W. and H.M.A.; Software, S.Z. and H.M.A.; Supervision, S.A.E. and T.A.W.; Validation, N.R. and A.T.A.; Visualization, A.T.A. and H.I.U.; Writing—original draft, M.S.B., S.A.E., S.Z. and F.S.; Writing—review & editing, S.A.E. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The required data will be available upon reasonable request. All the relevant is included in the Supplementary File.

Acknowledgments

The authors extend their appreciation to researchers supporting project number (RSP-2021/357) King Saud University, Riyadh Saudi Arabia for funding this research.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures, Scheme and Tables

Scheme 1. Activation pathway of enzymes VEGFR and EGFR. In the above figure, SOS is representing the son of sevenless, which is an RAS-activating factor. Flt 1 is a VEGFR receptor as Flt1 is a gene encoding VEGFR receptor, ERK (Extracellular signal-regulated kinase), Nrp1 is a cell surface protein and Akt is a serine/threonine kinase. Activation mechanism is denoted by thick and thin arrows, whereas double vertical lines are representing EGFR receptor [12,13].

Figure 1. Already reported FDA approved oxadiazoles [14,15].

Figure 2. Chemical structures of the selected 1, 3, 4 oxadiazole derivatives [16].

Figure 3. Schematic Diagram for the synthesis of 1, 3, 4 oxadiazoles [16].

Figure 4. Optimized structures along with different energy transitions for 7g, 7j and 7l. The transition of energy is depicted in green and red color around the atoms of evaluated compounds. The green color indicates the negative periods, while the red color represents the positive periods. The optimized structures of selected 1, 3, 4 oxadiazoles derivatives are given in (Supplementary File Figures S1–S3).

Figure 5. 3D diagram showing the binding interactions of 7g and 7j and 7l within the active pocket of Protein VEGFR2.

Figure 6. 2D binding interactions of 7g, 7j and 7l in best conformation with (A): VEGFR2 and (B): EGFR.

Figure 7. 3D diagram showing the binding interactions of 7g, 7j and 7l within the active pocket of Protein EGFR.

Figure 8. SeeSAR analysis of the most potent compounds with both proteins: (A) VEGFR; and (B) EGFR.

Figure 9. Root mean square deviation (RMSD) of the C-alpha atoms of 7j-VEGFR 2 Complex and 7j-EGFR complex with time.

Figure 10. Residue wise Root Mean Square Fluctuation (RMSF) of 7j-VEGFR 2 and 7j-EGFR complex with time.

Figure 11. Radius of gyration (Rg), Solvent accessible surface area (SASA) and Intra-molecular hydrogen bonds (IntraHB) graphs of VEGFR 2 and EGFR.

Figure 12. Residue index throughout the protein structure of complexes 7j-VEGFR 2 and 7j-EGFR showing alpha helices (Red) and beta-strands (Blue).

Figure 13. Protein-ligand contact histograms 7j-VEGFR 2 and 7j-EGFR showing various interactions.

Figure 14. Ligand atom interactions with the amino acid residues of target proteins. (7j-VEGFR 2 and 7j-EGFR).

Figure 15. Correlation plot between docking score (kcal/mol) and experimental potency.

Figure 16. The depiction of % growth reduction of HeLa and MCF-7 cell lines by 1, 3, 4oxadiazoles derivatives. The results were measured for 24 h and 48 h.

Table 1

Optimization energy and dipole moment of selected compounds in gaseous and solvent phase.

Codes	Gas		Ethanol
	Optimization Energy (Hatree)	Dipole Moment (Debye)	Optimization Energy (Hatree)	Dipole Moment(Debye)
7a	−1828.415	1.982	−1828.431	4.061
7b	−1867.705	2.210	−1867.721	4.005
7c	−1867.707	2.124	−1867.722	4.471
7d	−1867.707	2.258	−1867.722	4.394
7e	−1906.995	2.371	−1907.011	4.058
7f	−1906.996	2.565	−1907.012	4.316
7g	−1906.996	2.608	−1907.012	5.233
7h	−1906.997	2.404	−1907.013	4.948
7i	−1906.999	2.427	−1907.014	4.821
7j	−1906.999	2.205	−1907.014	4.724
7k	−1906.994	2.903	−1907.008	4.694
7l	−1906.995	2.308	−1907.011	4.211
7m	−1982.149	2.869	−1982.165	4.163
7n	−1982.151	2.933	−1982.166	4.787

Table 2

Various Global descriptors and their calculated values for selected compounds.

Code	EHOMO(eV)	ELUMO(eV)	∆Egap(eV)	PotentialIonization I (eV)	Affinity A (eV)	ElectronDonating Power(ω-)	ElectronAccepting Power (ω+)	Electrophilicity (Δω±)
7a	−0.2370	−0.0731	0.163	0.23707	0.0731	0.234	0.079	0.314
7b	−0.2406	−0.0713	0.169	0.24068	0.0713	0.232	0.076	0.308
7c	−0.2337	−0.0727	0.161	0.07275	0.23379	0.079	0.232	0.312
7d	−0.2289	−0.0723	0.156	0.22892	0.07236	0.230	0.079	0.309
7e	−0.2392	−0.0703	0.168	0.23929	0.07037	0.230	0.075	0.305
7f	−0.2378	−0.0706	0.167	0.23787	0.07069	0.230	0.076	0.306
7g	−0.2386	−0.0707	0.167	0.23866	0.07077	0.230	0.076	0.306
7h	−0.2403	−0.0712	0.169	0.24036	0.07125	0.232	0.076	0.308
7i	−0.2255	−0.0719	0.153	0.07197	0.22553	0.079	0.228	0.307
7j	−0.2320	−0.0721	0.177	0.23009	0.07217	0.230	0.079	0.309
7k	−0.2396	−0.0707	0.168	0.23966	0.07073	0.231	0.076	0.306
7l	−0.2294	−0.0723	0.157	0.22947	0.07231	0.230	0.036	0.199
7m	−0.2312	−0.0689	0.162	0.23125	0.06894	0.177	0.044	0.221
7n	−0.2138	−0.0716	0.142	0.21381	0.07161	0.185	0.053	0.238

Table 3

Chemical properties including softness and hardness values of selected compounds.

Code	Hardness(η)	Softness(ζ)	Electronegativity(χ)	Chemical Potential(μ)	Electrophilicity Index (ω)
7a	0.082	6.099	0.155	−0.155	0.147
7b	0.085	5.904	0.156	−0.156	0.144
7c	−0.081	−6.210	0.153	−0.153	−0.146
7d	0.078	6.387	0.151	−0.151	0.145
7e	0.084	5.920	0.155	−0.155	0.142
7f	0.084	5.982	0.154	−0.154	0.142
7g	0.084	5.956	0.155	−0.155	0.143
7h	0.085	5.913	0.156	−0.156	0.144
7i	−0.077	−6.512	0.149	−0.149	−0.144
7j	0.085	6.332	0.151	−0.151	0.145
7k	0.084	5.920	0.155	−0.155	0.143
7l	0.079	6.363	0.151	−0.151	0.145
7m	0.081	6.161	0.150	−0.150	0.139
7n	0.071	7.032	0.143	−0.143	0.143

Table 4

Docking score of selected compounds by considering bound and unbound states of the ligand along with their predicted inhibition constant (pIC₅₀).

Code	Vascular Endothelial Growth Factor (VEGFR)		Endothelial Growth Factor Receptor (EGFR)		Selectivity for VEGFR2
	Docking Score (kJ/mole)	Predicted Inhibition Constant (µM)	Docking Score (kJ/mole)	Predicted Inhibition Constant (µM)
7a	−41.79	0.047	−32.31	2.20	46.80
7b	−42.51	0.035	−32.43	2.05	58.57
7c	−44.39	0.016	−32.93	1.67	104.37
7d	−43.38	0.024	−33.18	1.50	62.5
7e	−44.18	0.018	−34.65	0.084	4.66
7f	−40.93	3.77	−31.64	2.80	0.74
7g	−46.32	0.009	−31.01	3.65	405.55
7h	−44.54	0.014	−31.93	2.51	179.28
7i	−44.78	0.017	−33.02	0.718	42.23
7j	−48.89	0.009	−33.23	1.50	166.66
7k	−43.51	0.023	−31.72	2.73	118.69
7l	−45.01	0.012	−34.19	1.01	84.16
7m	−42.63	0.033	−27.29	6.45	195.45
7n	−42.92	0.029	−32.68	1.86	64.13
VEGFR2 *	−51.49	0.0009	-	-	-
EGFR *	-	-	−23.65	71.18	-

Selectivity index defined as IC [EGFR]/IC [VEGF]. * control ligands.

Table 5

Values for the % Growth reduction after 24 and 48 h for HeLa and MCF-7 cells.

Code	% Growth Reduction of HeLa Cells		%Growth Reduction MCF-7 Cells
	After 24 h (µM)	After 48 h (µM)	After 24 h (µM)	After 48 h (µM)
7a	30.1 ± 1.23	55.6 ± 2.81	10.3 ± 0.88	49.3 ± 1.66
7b	50.1 ± 1.98	75.2 ± 3.11	54.2 ± 1.78	65.1 ± 2.88
7c	35.8 ± 0.89	51.2 ± 1.98	66.1 ± 2.17	75.2 ± 1.67
7d	60.2 ± 2.89	87.4 ± 2.34	58.3 ± 1.78	83.1 ± 3.11
7e	56.8 ± 1.67	79.7 ± 2.99	69.3 ± 0.76	78.9 ± 1.67
7f	66.3 ± 4.11	72.1 ± 1.78	71.2 ± 2.11	89.9 ± 3.11
7g	70.6 ± 2.98	98.1 ± 3.01	66.1 ± 1.56	96.7 ± 2.89
7h	51.4 ± 2.88	76.9 ± 1.89	28.9 ± 1.67	67.5 ± 1.88
7i	62.6 ± 2.66	82.5 ± 2.04	54.7 ± 2.09	78.9 ± 2.88
7j	67.9 ± 2.78	97.8 ± 2.98	77.2 ± 1.22	89.6 ± 1.98
7k	45.1 ± 2.31	57.9 ± 2.88	81.2 ± 2.12	94.3 ± 2.87
7l	56.1 ± 2.11	89.9 ± 2.65	88.8 ± 1.33	97.5 ± 3.01
7m	43.1 ± 1.99	73.1 ± 2.11	62.1 ± 1.11	88.9 ± 2.11
7n	38.9 ± 0.88	66.9 ± 1.56	57.3 ± 1.58	89.3 ± 1.67
Cisplatin	79.2 ± 2.44	89.1 ± 2.44	85.2 ± 2.11	98.2 ± 1.34

Table 6

Various physicochemical properties of selected compounds.

Physicochemical Properties
	Molecular Weight	Density	nHA	nHD	TPSA	LogS	LogP	LogD
7a	359.05	1.064	5	1	68.02	−5.279	4.237	4.0
7b	373.07	1.052	5	1	68.02	−5.194	4.442	4.11
7c	373.07	1.052	5	1	68.02	−5.572	4.69	4.176
7d	373.07	1.052	5	1	68.02	−5.582	4.697	4.2
7e	387.08	1.04	5	1	68.02	−5.371	4.9114	4.356
7f	387.08	1.04	5	1	68.02	−5.484	4.868	4.319
7g	387.08	1.04	5	1	68.02	−5.462	4.914	4.343
7h	387.08	1.04	5	1	68.02	−5.039	4.333	4.287
7i	387.08	1.04	5	1	68.02	−5.717	5.178	4.426
7j	387.08	1.04	5	1	68.02	−5.865	5.135	4.265
7k	387.08	1.04	5	1	68.02	−5.398	4.812	4.499
7l	387.08	1.04	5	1	68.02	−5.708	5.031	4.532
7m	403.08	1.058	6	1	77.25	−5.38	4.439	4.269
7n	403.08	1.058	6	1	77.25	−5.622	4.642	4.288

Table 7

Absorption and distribution properties of selected compounds.

Absorbtion & Distribution Properties
	Volume of Distribution (vd) Liters	Human IntestinalAbsorption (hia)	Caco-2 Permeability (log)	Blood Brain Barrier (bbb) & Blood-Placenta Barrier (bpb)	Plasma Protein Binding (ppb) %	pgp-Inhibitor	p-Glycoprotein Substrate (pgp-Substrate)	Mdck Permeability (cm/s)
7a	2.52	0.008	−4.668	0.218	99.08	0.087	0.001	1.3 × 10−5
7b	2.442	0.006	−4.584	0.176	99.53	0.008	0.001	1.3 × 10−5
7c	2.398	0.006	−4.604	0.185	99.73	0.057	0.003	1.2 × 10−5
7d	2.505	0.006	−4.589	0.165	99.52	0.069	0.002	1.2 × 10−5
7e	2.243	0.006	−4.581	0.187	96.94	0.004	0.004	1.3 × 10−5
7f	2.155	0.005	−4.513	0.172	98.12	0.045	0.004	1.2 × 10−5
7g	2.179	0.005	−4.508	0.161	100	0.056	0.005	1.2 × 10−5
7h	2.817	0.005	−4.548	0.159	99.99	0.007	0.002	1.3 × 10−5
7i	2.418	0.006	−4.611	0.16	100	0.015	0.008	1.2 × 10−5
7j	1.701	0.005	−4.557	0.13	100	0.33	0.015	1.1 × 10−5
7k	2.964	0.005	−4.539	0.167	99.86	0.057	0.001	1.3 × 10−5
7l	2.846	0.006	−4.542	0.166	99.94	0.341	0.002	1.2 × 10−5
7m	2.629	0.005	−4.513	0.077	99.89	0.717	0.001	1.3 × 10−5
7n	2.854	0.006	−4.539	0.545	100	0.627	0.001	1.2 × 10−5

Table 8

Metabolism and Excretion values of selected compounds.

Metabolism							Excretion
	CYP1A2 Inhibitor	CYP2C19 Inhibitor	CYP2C9 Inhibitor	CYP2D6 Inhibitor	CYP3A4 Inhibitor	CL (mL/min.)	T1/2 (h)
7a	0.927	0.962	0.886	0.527	0.378	3.452	0.136
7b	0.824	0.956	0.9	0.477	0.704	2.956	0.135
7c	0.875	0.954	0.892	0.678	0.809	3.738	0.112
7d	0.758	0.936	0.861	0.63	0.55	3.486	0.099
7e	0.758	0.949	0.911	0.584	0.833	3.26	0.112
7f	0.632	0.936	0.886	0.509	0.808	3.263	0.126
7g	0.743	0.944	0.903	0.597	0.847	3.257	0.103
7h	0.687	0.952	0.902	0.429	0.914	2.551	0.166
7i	0.78	0.939	0.902	0.675	0.801	3.802	0.08
7j	0.791	0.951	0.908	0.69	0.837	4.106	0.13
7k	0.869	0.951	0.919	0.588	0.644	2.945	0.125
7l	0.768	0.93	0.875	0.72	0.467	3.433	0.094
7m	0.776	0.958	0.922	0.454	0.595	2.768	0.129
7n	0.607	0.925	0.809	0.532	0.334	3.279	0.081

Table 9

Medicinal properties and toxicity profile of selected compounds.

Medicinal Properties			Toxicity
	Synthetic Accessibility Score	Lipinski Rule	AMES Toxicity	Carcinogenicity	Eye Corrosion	Eye Irritation	Respiratory Toxicity
7a	2.331	None	0.010	0.817	0.003	0.023	0.965
7b	2.393	None	0.022	0.838	0.003	0.027	0.954
7c	2.426	None	0.014	0.846	0.003	0.018	0.952
7d	2.379	None	0.015	0.852	0.003	0.019	0.942
7e	2.49	None	0.025	0.859	0.003	0.027	0.933
7f	2.465	None	0.016	0.866	0.003	0.022	0.933
7g	2.477	None	0.016	0.857	0.003	0.02	0.927
7h	2.521	None	0.023	0.865	0.003	0.018	0.944
7i	2.471	None	0.012	0.881	0.003	0.021	0.926
7j	2.547	None	0.01	0.857	0.003	0.018	0.947
7k	2.484	None	0.019	0.805	0.003	0.019	0.922
7l	2.428	None	0.013	0.846	0.003	0.017	0.89
7m	2.476	None	0.019	0.783	0.003	0.022	0.872
7n	2.435	None	0.010	0.836	0.003	0.016	0.867

Table 10

Various Toxicological parameters of selected compounds.

TOX21 Pathway
	NR-AR		NR-AR-LBD		NR-ER		Antioxidant Response Element
	Result	Probability	Result	Probability	Result	Probability	Result	Probability
7a	+	0.01	+	0.156	-	0.751	+	0.921
7b	+	0.013	+	0.076	+	0.675	+	0.91
7c	+	0.011	+	0.035	-	0.745	+	0.916
7d	+	0.011	+	0.042	-	0.78	+	0.924
7e	+	0.021	+	0.038	+	0.61	-	0.912
7f	+	0.015	+	0.025	+	0.693	+	0.906
7g	+	0.015	+	0.024	+	0.682	+	0.906
7h	+	0.023	+	0.024	-	0.511	+	0.88
7i	+	0.02	+	0.03	-	0.732	+	0.917
7j	+	0.011	+	0.014	-	0.745	+	0.906
7k	+	0.012	+	0.058	-	0.695	+	0.908
7l	+	0.011	+	0.04	-	0.801	+	0.926
7m	+	0.012	+	0.271	-	0.65	+	0.921
7n	+	0.004	+	0.167	+	0.805	+	0.936

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/biom12111612/s1, Figure S1: Optimized structures of the oxadiazoles derivatives (7a–7h) in gas phase; Figure S2: Optimized structures of the oxadiazoles derivatives (7i–7n) in gas phase; Figure S3: HOMO and LUMO structures of all the oxadiazoles derivatives in gas phase; Figure S4: 3D binding interactions of all compounds with VEGFR; Figure S5: 2D binding interactions of all compounds with VEGFR; Figure S6: 3D binding interactions of all compounds with EGFR; Figure S7: 2D binding interactions of all compounds with EGFR; Figure S8: Simulation replicate of Vegfr-7j Complex; Figure S9: Simulation replicate of Egfr-7j Complex; Table S1: Chemical structures of oxadiazoles with IUPAC names.

References

1. Bray, F.; Laversanne, M.; Weiderpass, E.; Soerjomataram, I. The ever-increasing importance of cancer as a leading cause of premature death worldwide. Cancer; 2021; 127, pp. 3029-3030. [DOI: https://dx.doi.org/10.1002/cncr.33587] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34086348]

2. Brown, M.L.; Lipscomb, J.; Snyder, C. The burden of illness of cancer: Economic cost and quality of life. Annu. Rev. Public Health; 2001; 22, 91. [DOI: https://dx.doi.org/10.1146/annurev.publhealth.22.1.91] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11274513]

3. Huang, J.; Leung, D.K.W.; Chan, E.O.T.; Lok, V.; Leung, S.; Wong, I.; Wong, M.C. A global trend analysis of kidney cancer incidence and mortality and their associations with smoking, alcohol consumption, and metabolic syndrome. Eur. Urol. Focus; 2022; 8, pp. 200-209. [DOI: https://dx.doi.org/10.1016/j.euf.2020.12.020]

4. Murai, M.; Oya, M. Renal cell carcinoma: Etiology, incidence and epidemiology. Curr. Opin. Urol.; 2004; 14, pp. 229-233. [DOI: https://dx.doi.org/10.1097/01.mou.0000135078.04721.f5] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15205579]

5. Roskoski, R., Jr. Src protein-tyrosine kinase structure and regulation. Biochem. Biophys. Res. Commun.; 2004; 324, pp. 1155-1164. [DOI: https://dx.doi.org/10.1016/j.bbrc.2004.09.171] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15504335]

6. Normanno, N.; De Luca, A.; Bianco, C.; Strizzi, L.; Mancino, M.; Maiello, M.R.; Salomon, D.S. Epidermal growth factor receptor (EGFR) signaling in cancer. Gene; 2006; 366, pp. 2-16. [DOI: https://dx.doi.org/10.1016/j.gene.2005.10.018]

7. Raymond, E.; Faivre, S.; Armand, J.P. Epidermal growth factor receptor tyrosine kinase as a target for anticancer therapy. Drugs; 2000; 60, pp. 15-23. [DOI: https://dx.doi.org/10.2165/00003495-200060001-00002]

8. Xu, Y.; Liu, H.; Chen, J.; Zhou, Q. Acquired resistance of lung adenocarcinoma to EGFR-tyrosine kinase inhibitors gefitinib and erlotinib. Cancer Biol. Ther.; 2010; 9, pp. 572-582. [DOI: https://dx.doi.org/10.4161/cbt.9.8.11881]

9. Veikkola, T.; Karkkainen, M.; Claesson-Welsh, L.; Alitalo, K. Regulation of angiogenesis via vascular endothelial growth factor receptors. Cancer Res.; 2000; 60, pp. 203-212.

10. Wise, L.M.; Veikkola, T.; Mercer, A.; Savory, L.J.; Fleming, S.B.; Caesar, C.; Vitali, A.; Makinen, T.; Alitalo, K.; Stacker, S.A. Vascular endothelial growth factor (VEGF)-like protein from orf virus NZ2 binds to VEGFR2 and neuropilin-1. Proc. Natl. Acad. Sci. USA; 1999; 96, pp. 3071-3076. [DOI: https://dx.doi.org/10.1073/pnas.96.6.3071]

11. Wang, Z.; Wang, N.; Han, S.; Wang, D.; Mo, S.; Yu, L.; Chen, J. Dietary compound isoliquiritigenin inhibits breast cancer neoangiogenesis via VEGF/VEGFR-2 signaling pathway. PLoS ONE; 2013; 8, e68566. [DOI: https://dx.doi.org/10.1371/journal.pone.0068566]

12. Joensuu, H.; Trent, J.C.; Reichardt, P. Practical management of tyrosine kinase inhibitor-associated side effects in GIST. Cancer Treat. Rev.; 2011; 37, pp. 75-88. [DOI: https://dx.doi.org/10.1016/j.ctrv.2010.04.008] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20570050]

13. Lichtenberger, B.M.; Tan, P.K.; Niederleithner, H.; Ferrara, N.; Petzelbauer, P.; Sibilia, M. Autocrine VEGF Signaling Synergizes with EGFR in Tumor Cells to Promote Epithelial Cancer Development. Cell; 2010; 140, pp. 268-279. [DOI: https://dx.doi.org/10.1016/j.cell.2009.12.046]

14. Siwach, A.; Verma, P.K. Therapeutic potential of oxadiazole or furadiazole containing compounds. BMC Chem.; 2020; 14, pp. 1-40. [DOI: https://dx.doi.org/10.1186/s13065-020-00721-2] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33372629]

15. Dokla, E.M.; Fang, C.-S.; Abouzid, K.A.; Chen, C.S. 1,2,4-Oxadiazole derivatives targeting EGFR and c-Met degradation in TKI resistant NSCLC. Eur. J. Med. Chem.; 2019; 182, 111607. [DOI: https://dx.doi.org/10.1016/j.ejmech.2019.111607]

16. Iftikhar, M.; Saleem, M.; Riaz, N.; Ahmed, I.; Rahman, J.; Ashraf, M.; Htar, T.T. A novel five-step synthetic route to 1, 3, 4-oxadiazole derivatives with potent α-glucosidase inhibitory potential and their in-silico studies. Arch. Der Pharm.; 2019; 352, 1900095. [DOI: https://dx.doi.org/10.1002/ardp.201900095] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31544284]

17. Frisch, M.-J.; Trucks, G.-W.; Schlegel, H.-B.; Scuseria, G.-E.; Robb, M.-A.; Cheeseman, J.-R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G.-A. et al. Gaussian 09; Gaussian, Inc.: Wallingford, CT, USA, 2009.

18. Dennington, R.; Keith, T.A.; Millam, J.M. GaussView V 6.1; Semichem Inc.: Shawnee Mission, KS, USA, 2016.

19. Available online: https//www.rcsb.org/3VHE (accessed on 1 July 2022).

20. Available online: https://www.rcsb.org/2GS6 (accessed on 1 July 2022).

21. Morris, G.M.; Huey, R.; Lindstrom, W.; Sanner, M.F.; Belew, R.K.; Goodsell, D.S.; Olson, A.J. AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility. J. Comput. Chem.; 2009; 30, pp. 2785-2791. [DOI: https://dx.doi.org/10.1002/jcc.21256] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19399780]

22. Accelrys Software Inc. Discovery Studio Visualizer, 2; Accelrys Software Inc.: San Diego, CA, USA, 2005.

23. BioSolveIT-SeeSAR. Available online: https://www.biosolveit.de/SeeSAR/ (accessed on 15 July 2022).

24. CambridgeSoft. ChemDraw Ultra 12.0 0 (Copyright) 1986 to 2009; CambridgeSoft Corp.: Cambridge, MA, USA, 2009.

25. CambridgeSoft. Chem 3D Pro 12.0 (Copyright) 1986 to 2009; CambridgeSoft Corp.: Cambridge, MA, USA, 2009.

26. Ferreira, L.G.; Dos Santos, R.N.; Oliva, G.; Andricopulo, A.D. Molecular Docking and Structure-Based Drug Design Strategies. Molecules; 2015; 20, pp. 13384-13421. [DOI: https://dx.doi.org/10.3390/molecules200713384] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26205061]

27. Zentgraf, M.; Steuber, H.; Koch, C.; La Motta, C.; Sartini, S.; Sotriffer, C.A.; Klebe, G. How reliable are current docking approaches for structure-based drug design? Lessons from aldose reductase. Angew. Chem. Int. Ed.; 2007; 46, pp. 3575-3578. [DOI: https://dx.doi.org/10.1002/anie.200603625]

28. Yusuf, D.; Davis, A.M.; Kleywegt, G.J.; Schmitt, S. An alternative method for the evaluation of docking performance: RSR vs. RMSD. J. Chem. Inf. Model.; 2008; 48, pp. 1411-1422. [DOI: https://dx.doi.org/10.1021/ci800084x]

29. Bowers, K.J.; Chow, D.E.; Xu, H.; Dror, R.O.; Eastwood, M.P.; Gregersen, B.A.; Shaw, D.E. Scalable algorithms for molecular dynamics simulations on commodity clusters. Proceedings of the SC’06: 2006 ACM/IEEE Conference on Supercomputing; Tampa, FL, USA, 11–17 November 2006; 43.

30. Hildebrand, P.W.; Rose, A.S.; Tiemann, J.K. Bringing Molecular Dynamics Simulation Data into View. Trends Biochem. Sci.; 2019; 44, pp. 902-913. [DOI: https://dx.doi.org/10.1016/j.tibs.2019.06.004]

31. Rasheed, M.; Iqbal, M.; Saddick, S.; Ali, I.; Khan, F.; Kanwal, S.; Ahmed, D.; Ibrahim, M.; Afzal, U.; Awais, M. Identification of Lead Compounds against Scm (fms10) in Enterococcus faecium Using Computer Aided Drug Designing. Life; 2021; 11, 77. [DOI: https://dx.doi.org/10.3390/life11020077] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33494233]

32. Shivakumar, D.; Williams, J.; Wu, Y.; Damm, W.; Shelley, J.; Sherman, W. Prediction of Absolute Solvation Free Energies using Molecular Dynamics Free Energy Perturbation and the OPLS Force Field. J. Chem. Theory Comput.; 2010; 6, pp. 1509-1519. [DOI: https://dx.doi.org/10.1021/ct900587b] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26615687]

33. Anderson, P.M.; Wilson, M.R. Developing a force field for simulation of poly (ethylene oxide) based upon ab initio calculations of 1,2-dimethoxyethane. Mol. Phys.; 2005; 103, pp. 89-97. [DOI: https://dx.doi.org/10.1080/00268970412331293811]

34. Mosmann, T. Rapid colorimetric assay for cellular growth and survival: Application to proliferation and cytotoxicity assays. J. Immunol. Methods; 1983; 65, pp. 55-63. [DOI: https://dx.doi.org/10.1016/0022-1759(83)90303-4]

35. Nikš, M.; Otto, M. Towards an optimized MTT assay. J. Immunol. Methods; 1990; 130, pp. 149-151. [DOI: https://dx.doi.org/10.1016/0022-1759(90)90309-J]

36. Iqbal, J.; Ejaz, S.A.; Saeed, A.; Al-Rashida, M. Detailed investigation of anticancer activity of sulfamoyl benz (sulfon) amides and 1H-pyrazol-4-yl benzamides: An experimental and computational study. Eur. J. Pharmacol.; 2018; 5, pp. 11-24. [DOI: https://dx.doi.org/10.1016/j.ejphar.2018.05.011]

37. ADMETlab 2.0. Available online: https://admetmesh.scbdd.com/ (accessed on 5 October 2022).

38. McTigue, M.A.; Wickersham, J.A.; Pinko, C.; Showalter, R.E.; Parast, C.V.; Tempczyk-Russell, A.; Appelt, K. Crystal structure of the kinase domain of human vascular endothelial growth factor receptor 2: A key enzyme in angiogenesis. Structure; 1999; 7, pp. 319-330. [DOI: https://dx.doi.org/10.1016/S0969-2126(99)80042-2]

39. Zhang, X.; Gureasko, J.; Shen, K.; Cole, P.A.; Kuriyan, J. An allosteric mechanism for activation of the kinase domain of epidermal growth factor receptor. Cell; 2006; 125, pp. 1137-1149. [DOI: https://dx.doi.org/10.1016/j.cell.2006.05.013]

40. Sargolzaei, M. Effect of nelfinavir stereoisomers on coronavirus main protease: Molecular docking, molecular dynamics simulation and MM/GBSA study. J. Mol. Graph. Model.; 2021; 103, 107803. [DOI: https://dx.doi.org/10.1016/j.jmgm.2020.107803]

41. Barbosa, L.R.; Ortore, M.G.; Spinozzi, F.; Mariani, P.; Bernstorff, S.; Itri, R. The importance of protein-protein interactions on the pH-induced conformational changes of bovine serum albumin: A small-angle X-ray scattering study. Biophys. J.; 2010; 98, pp. 147-157. [DOI: https://dx.doi.org/10.1016/j.bpj.2009.09.056] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20085727]

42. Badieyan, S.; Bevan, D.R.; Zhang, C. Study and design of stability in GH5 cellulases. Biotechnol. Bioeng.; 2012; 109, pp. 31-44. [DOI: https://dx.doi.org/10.1002/bit.23280] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21809329]

43. Davis, A.M.; Teague, S.J. Hydrogen bonding, hydrophobic interactions, and failure of the rigid receptor hypothesis. Angew. Chem. Int. Ed.; 1999; 38, pp. 736-749. [DOI: https://dx.doi.org/10.1002/(SICI)1521-3773(19990315)38:6<736::AID-ANIE736>3.0.CO;2-R]

44. Rajagopal, S.; Vishveshwara, S. Short hydrogen bonds in proteins. FEBS J.; 2005; 272, pp. 1819-1832. [DOI: https://dx.doi.org/10.1111/j.1742-4658.2005.04604.x] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15819878]

45. Aziz, M.; Ejaz, S.A.; Tamam, N.; Siddique, F.; Riaz, N.; Qais, F.A.; Iqbal, J. Identification of potent inhibitors of NEK7 protein using a comprehensive computational approach. Sci. Rep.; 2022; 12, 6404. [DOI: https://dx.doi.org/10.1038/s41598-022-10253-5]