Unleashing Nature’s potential: a computational approach to discovering novel VEGFR-2 inhibitors from African natural compound using virtual screening, ADMET analysis, molecular dynamics, and MMPBSA calculations

One of the characteristic features of cancer is angiogenesis, the process by which new, aberrant blood vessels are formed from pre-existing blood vessels. The process of angiogenesis begins when VEGF binds to its receptor, the VEGF receptor (VEGFR). The formation of new blood vessels provides nutrients that can promote the growth of cancer cells. When it comes to new blood vessel formation, VEGFR2 is a critical player. Therefore, inhibiting VEGFR2 is an effective way to target angiogenesis in cancer treatment. The aim of our research was to find new VEGFR-2 inhibitors by performing a virtual screening of 13313 from African natural compounds using different in silico techniques. Using molecular docking calculations and ADMET properties, we identified four compounds that exhibited a binding affinity ranging from −11.0 kcal/mol to −11.5 Kcal/mol when bound to VEGFR-2. These four compounds were further analyzed with 100 ns simulations to determine their stability and binding energy using the MM-PBSA method. After comparing the compounds with Regorafenib, a drug approved for anti-angiogenesis treatment, it was found that all the candidates (EANPDB 252, NANPDB 4577, and NANPDB 4580), with the exception of EANPDB 76, could target VEGFR-2 similarly effectively to Regorafenib. Therefore, we recommend three of these agents for anti-angiogenesis treatment because they are likely to deactivate VEGFR-2 and thus inhibit angiogenesis. However, it should be noted that the safety and suitability of these agents for clinical use needs further investigation, as the computer-assisted study did not include in vitro or in vivo experiments.

Introduction

Angiogenesis is a critical pathogenic process in many disease conditions (Niu and Chen, 2010). For instance, in cancer, angiogenesis is essential for the development and growth of the cancer cells (Nieves et al., 2009). Normally, the cancer cells relied on the formation of new blood vessels as the source of oxygen and nutrients for the developing tumor (Nishida et al., 2006). The process of angiogenesis begins when Vascular Endothelial Growth Factor (VEGF) binds to its receptors particularly the VEGFR-2 receptor (Nieves et al., 2009). Specifically, the overexpression of VEGF is associated with autophosphorylation of VEGFR-2 receptor in malignancy (Nieves et al., 2009). The VEGF and other growth factors produced by the tumor results in the production of new blood vessels which allows the cells to grow exponentially (Rajabi and Mousa, 2017). Under the influence of VEGF, the vasculature formed are abnormal resulting into abnormal conditions. Hence, the role of this growth factor in vasculature makes it a target for cancer treatment (Carmeliet, 2005).

The main goal of inhibiting tumor angiogenesis is to deprive cancer cells of the nutrients and oxygen they need to grow (Rajabi and Mousa, 2017). The therapeutic value of VEGFR-targeted cancer therapy is described by a large body of clinical evidence. When the VEGFR-2 pathway is blocked, it has a significant anti-angiogenic effect on human cancer (Ansari et al., 2022). Currently, drugs designed to block VEGF and its receptor are approved for the treatment of cancer and ocular diseases. These drugs include anti-VEGF antibodies such as bevacizumab, ranibizumab, and pegaptanib, and VEGFR inhibitors such as sunitinib, sorafenib, Regorafenib, and pazopanib (Carmeliet and Jain, 2011; Yoo and Kwon, 2013; Schmieder et al., 2014). Despite their efficacy, these drugs have limited efficacy and may develop resistance over time. For example, bevacizumab is known to cause severe inflammation in the eye (Sharma et al., 2012), while sunitinib can cause various side effects, such as thrombocytopenia and hypertension, in the treatment of metastatic renal cell carcinoma (Motzer et al., 2006).

Traditional medicine has been used for thousands of years, with a close relationship between medicinal and natural products. Extracts from different parts of medicinal plants were traditionally used in the treatment of different ailments (Rizvi et al., 2022). The awareness of the healing capacity of medicinal plants led to the discovery of phytochemicals (Thakur et al., 2020). The advent of technology led to the isolation of these chemicals that were believed to play a role in the plant healing process (Sasidharan et al., 2011). These natural compounds were reported to have several biological activities and serve as leads for the discovery of new drug candidates. For instance, some of the compounds were reported to exhibit anticancer (Prakash et al., 2013), antidiabetic (Ivorra et al., 1989; Li et al., 2004), antioxidant (Lourenço et al., 2019), anti-inflammatory (Ola et al., 2020), antitumor activities among others. Hence, natural products serve as the richest source of new classes of molecules for biological research. These molecules will continue to play a pivotal role in ethnopharmacology (Beutler, 2009). It is believed that many of these compounds are yet to be discovered. Also, the biological and pharmacological activities of the known compounds are not fully explored (Dias et al., 2012).

Although biochemical and cellular assays can evaluate more synthetic and natural compounds, their empirical screening is limited (McInnes, 2007), making computational methods important in drug discovery. Virtual screening allows searching for hits from chemical collections acquired for biological activity assays. By filtering vast virtual libraries, virtual screening based on molecular docking with a target protein of empirically known structure is becoming a standard approach for identifying promising lead molecules in drug discovery projects (Jain, 2004; Shoichet, 2004). To better understand the physical basis of interactions between protein receptors and their small molecule inhibitors, molecular dynamics simulation has become a fundamental method in drug discovery. The stability of receptor-ligand complexes can be evaluated by molecular dynamics simulations. In addition, this simulation considers the role played by individual amino acid residues.

Due to the limitations of currently available drugs, the development of small molecule inhibitors that can block VEGFR signaling is an attractive strategy to prevent angiogenesis. This approach is less likely to cause long-term toxicity (Byrne et al., 2005). In this study, we aim to identify potential VEGFR-2 inhibitors from natural sources using computational methods. This could help in finding alternative drug candidates that could be effective against the process of angiogenesis. The approach employed includes the collection of plant-based naturally occurring African compounds, followed by virtual screening of compounds using molecular docking approach against VEGFR-2 receptor to identify potential lead compounds. Finally, the lead compounds were subjected to molecular dynamics simulations and MMPBSA calculations to identify their stability Figure 1.

FIGURE 1

FIGURE 1. Workflow illustrating the screening process used to identify novel VEGFR-2 inhibitors.

Material and methods

Database preparation

We performed virtual screenings with compounds contained in the African Natural Products Database (ANPDB) (http://african-compounds.org/nanpdb/) and South African Natural Compounds Database (SANCDB) (https://sancdb.rubi.ru.ac.za/). It is important to note that the ANPDB contains information on compounds from different parts of Africa (Ntie-Kang et al., 2017; Simoben et al., 2020), while the SANCDB is specifically for natural compounds collected from South African (Beutler, 2009; Ola et al., 2020). Altogether, a total of 14332 natural compounds were retrieved from the databases. After removal of duplicate, the compounds reduced to 13313 (Ntie-Kang et al., 2017; Simoben et al., 2020). Thereafter, the compounds were downloaded in “Sdf” format protonated and subjected to energy minimization (Aljuaid et al., 2022). Subsequently, the compounds were converted to pdbqt format using the python script (mk_prepare_ligand.py found in MGLTools 1.5.6.) in preparation for molecular docking studies.

VEGFR-2 preparation

The crystal structures of the protein were retrieved from the Protein Data Bank with the respective PDB codes: 4ASD and 4AGC. A multiple sequence alignment demonstrated that both 4AGC and 4ASD exhibit 100% sequence similarity and are well aligned when superimposed. Both structures share conserved residues within their active sites and the docking results of Regorafenib, an FDA-approved drug, with both 4AGC and 4ASD showed similar interactions. Due to the fewer missing residues in 4ASD, we selected it for our subsequent studies. To address the missing residues, we employed “Build Structure” from Chimera’s modeling tools to rectify the protein structure. Prior to the docking analysis, we prepared the VEGFR-2 protein by removing water molecules, ions, and the ligand (Tabti et al., 2023). The protonation of the amino acids was analyzed on the H++ server, and the hydrogen atoms were added to the whole structure.

Molecular docking

Virtual screening of retrieved compounds against VEGFR-2 was conducted using AutoDock Vina (Abdelkader et al., 2022). Before commencing the docking studies, MGLTools 1.5.6 was used to convert the VEGFR-2 structure from its native format, pdb, to the docking-ready format, pdbqt by adding the polar hydrogen atoms and Gasteiger charges to the protein chain (Baammi et al., 2022). Except for the exhaustiveness parameter, which was set to 100, all default values were retained. The exhaustiveness plays a role in controlling thoroughness of the search space exploration during molecular docking. Although the default value of exhaustiveness is 8, increasing it to 100 led to an increase in the reproducibility of our docking results (Devaurs et al., 2019). The grid box spacing was set to 0.375 Å, the center to (−24.611 Å, −0.388 Å, −10.929 Å), and the lattice size to 20 Å × 20 Å × 20 Å. Nine poses were constructed for each protein-ligand complex based on docking affinity. The discovery Studio Viewer was used to display and analyze the docking results to find the important interactions between the ligands and the protein binding site (Baammi et al., 2023). In addition, the co-crystalline ligand was re-docked as an inhibitor of VEGFR-2 using the above parameters and values following by comparing the RMSD (root-mean-square deviation) of the heavy atoms between the docked pose and the crystallographic pose of the ligand (Xu and Meroueh, 2016). The nine best ranked poses were visualized using Discovery Studio Visualizer version 17.2 and PyMol version 1.1. The best docked complexes with the lowest docking score and the most favorable interactions were used for molecular dynamics simulations (Hosseini and Amanlou, 2020).

Pharmacokinetic (ADME) and toxicological predictions

The determination of compounds ADMET properties has a crucial impact in drug design and discovery (Guan et al., 2019). Information on these properties helps in understanding the pharmacokinetics and pharmacodynamics of a drug candidate (Walker, 2004). Of course, 60% of lead compounds fail during screening in the drug discovery pipeline due to their unacceptable ADMET properties. Hence, early prediction of these properties would result in substantial cost savings within the field of drug research (Mandlik et al., 2016). Some of the parameters evaluated for the ADMET includes adsorption, distribution, metabolism, excretion, and toxicity (En-nahli et al., 2022). To accurately predict drug efficacy, penetration through the blood-brain barrier, absorption in the human intestine, CNS permeability, inhibition of cytochrome P450 2D6 and 3A4, hepatotoxicity, and AMES toxicity are also considered (Lagorce et al., 2017). For this study, the ADMET parameters for the four lead compounds were predicted using the pkCSM web server (Pires et al., 2015). These compounds were further subjected to MD simulation studies.

Molecular dynamics simulation

The behavior of the selected ligands was investigated together with VEGFR-2 over a period of 100 ns using the GROMACS software package version 2019.3 (Park et al., 2020) on a high-performance cluster (POWEREDGE C6420, CRC-STACKHPC, XEON PLATNIUM 8276L 28C 2.2GHZ, MELLANOX INFINIBAND HDR100)). The CHARMM27 force field was used for the protein (Lindahl et al., 2010), and the topology for the ligands was generated using the Swissparam server (Zoete et al., 2011). Prior to neutralization in the system with counterions, each complex was resolved in a dodecahedral box (1.0 nm) using the TIP3P water model. The steepest descent method was used to achieve both the minimum energy and maximum force, with Fmax set to 1,000 kJ/mol/nm (Al-Khafaji and Taskin Tok, 2020). To equilibrate the system at 300 K and 1 bar, two 100 ps simulations were performed in rapid succession using canonical NVT and isobaric NPT ensembles. Subsequently, 100 ns molecular dynamics simulations were performed for each molecule. The output trajectories were generated, and the data files were analyzed to better understand the behavior of the protein.

MM-PBSA binding energy calculation

The Molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) was used to compute the binding free energies of the complexes (Kumari et al., 2014). The calculations were conducted using the g_mmpbsa script tool (reference (Simoben et al., 2020)), which employs an approach based on the average of two energy values: the solvation energy and the potential energy in a vacuum.

$∆\mathrm{E}\left(\text{MM}-\text{PBSA}\right)=∆{\mathrm{E}}_{\text{MM}}+∆{\mathrm{G}}_{\text{solvation}.}(1)$

In Equation 1, E_MM represents the potential energy in a vacuum, while Gsolvation corresponds to the free solvation energy. The molecular mechanical energy (E_MM) is determined by considering the contributions of the electrostatic component (E_ele) and the van der Waals interaction (E_vdW). The solvation energy is calculated using the polar solvation energy (G_pol) and the non-polar solvation energy (G_nonpol). The polar solvation energy (G_pol) is determined using the Poisson-Boltzmann equation (PB), while the non-polar solvation energy (G_nonpol) is evaluated based on the solvent-accessible surface area (SASA).

Result and discussion

Virtual screening and molecular docking

Computer-aided drug discovery (CADD) relies heavily on molecular docking as a fundamental technique. This method uses computer models to evaluate the interaction between a receptor and the multiple compounds in a virtual environment by using computational algorithms that calculate factors such as binding affinity, energy, and molecular interactions (Yu and Mackerell, 2017). This approach allows researchers to select promising lead molecules for further investigation. To discover potential compounds against VEGFR-2, a molecular docking simulation of the binding pocket of VEGFR-2 was performed using 13363 molecules from African plants. Using the docking method, all compounds were predicted, ranked from lowest to highest docking score, and compared to that of the co-crystallized ligand (−11.0 kcal/mol). Of the 13363 compounds with structural differences, only four compounds had a binding affinity ≤ −11.0 Kcal/mol as Regorafenib in terms of binding energy. Further in-depth molecular studies of the interactions of these ligands with VEGFR-2 are shown in Table 1. It should be noted that the best compounds are those with the lowest docking score and the most beneficial interactions that could fit well into the binding site (Hosseini and Amanlou, 2020).

TABLE 1

TABLE 1. The analyzed binding affinity for the best six African compounds against the VEGFR-2 protein.

The most effective ligand for docking with VEGFR-2 was Naringenin 7-p-coumaroylglucoside (NANPDB 4577) with a binding affinity of −11.5 kcal/mol. There is a strong interaction with VEGFR-2, as evidenced by the presence of four hydrogen bonds between Glu885, Ser884, Ala881, and Lys868, and several hydrophobic, electrostatic, and other interactions between Leu840, Leu889, Asp814, and Cys817. These molecular interactions of NANPDB 4577 with VEGFR-2 are shown in Figure 2. For Naringenin 7-p-coumaroyl glucoside isolated from Phlomis aurea, there is no evidence yet of therapeutic use in cancer therapy. Further investigation of the medicinal potential of this compound is recommended in view of our results.

FIGURE 2

FIGURE 2. 3D and 2D binding interactions of the selected compounds.

Other compounds, Lettowifuraquinone (EANPDB 76) and Abyquinone B (EANPDB 252), were best docked to VEGFR-2 with a binding affinity of −11.2 kcal/mol. They showed strong interactions with VEGFR-2. EANPDB 76 exhibits two hydrogen bonds with Val914, Asp1046 and several hydrophobic, electrostatic, and other interactions via Phe918, Ala866, Leu1035, Val899, Leu889, Cys1024, Ile888, Val848, Phe1047, Leu840, Ly868, Glu885, and Cys1045. However, EANPDB 252 exhibits three hydrogen bonds with Asp1046, His1026, Ala881, and several hydrophobic, electrostatic, and other interactions via Val898, Leu1019, Ile1044, Gly1048, and Cys817. The molecular interactions of EANPDB 76 and EANPDB 252 with VEGFR-2 are shown in Figure 2. These two compounds, Lettowifuraquinone and Abyquinone B, belong to the class of organic compounds known as quinones, which form an important class of cytotoxins against cancer, tumor, antimicrobials, and antiparasitic effects (Campos-Xolalpa et al., 2021). NANPDB 4580 (chrysoeriol 7-p-coumaroylglucoside) has the same binding affinity as the reference (−11.0 kcal/mol) with four hydrogen bonds with Arg842, Cys919, Asn923, and Ile1044 and eight hydrophobic bonds with Arg1051, Ala866, Val848, Val916, Val899, His1026, Val916, and Leu889, and one electrostatic bond with Cys1045. This molecule belongs to the class of Flavonoids that have been widely used as anticancer (Zhao et al., 2019), antimicrobial, antiviral, antiangiogenic (Camero et al., 2018; Zhao et al., 2018), antimalarial, antioxidants, neuroprotective, antitumor, and anti-proliferative agents (Patel et al., 2018).

Drug metabolism and toxicity profiling

Poor pharmacokinetic profiles and toxicity concerns are often cited as reasons for excluding lead compounds from preclinical and clinical investigations (Singh, 2006). Therefore, when lead compounds are discovered, it would be advisable to use in silico methods to predict the potential toxicity and pharmacokinetic properties of the hit compounds (Salifu et al., 2023). Therefore, the top four compounds from the virtual screening were subjected to evaluation of their absorption, distribution, metabolism, excretion, and toxicity. The ADMET properties of all compounds are listed in Table 2.

TABLE 2

TABLE 2. Absorption, distribution, metabolism, excretion, and toxicity (ADMET) prediction properties using the pkCSM server.

pKCSM is an online tool that we used as a reference for several different variables in our study. According to the data from the pKCSM server, a higher HIA value was obtained for compounds EANPDB 76, EANPDB 252, NANPDB 4577, and NANPDB 4580, suggesting that these compounds may be more absorbed in the intestine after oral administration. To determine whether a compound is mutagenic, the researchers used the AMES toxicity assay. All four compounds tested (EANPDB 76, EANPDB 252, NANPDB 4577, and NANPDB 4580) proved negative in the AMES toxicity test. According to hepatotoxicity prediction research, the new chemical poses no risk to liver function. The group of isoenzymes that includes cytochrome P450 (CYP) is the key factor in pharmacokinetics. It plays an important role in the metabolic processes of numerous substances, such as drugs, bile acids, steroid hormones, fatty acids and carcinogens. Some of the substances tested may be able to prevent certain cytochrome P450 isoforms from functioning properly. These isoforms are important enzymes involved in the metabolism of various drugs and can cause negative interactions and side effects when combined with other drugs. In silico data suggests that all tested compounds (EANPDB 76, EANPDB 252, NANPDB 4577, and NANPDB 4580) are neither substrates nor inhibitors of CYP2D6. However, EANPDB 76 was found to be an inhibitor and substrate of CYP3A4.

Molecular dynamics simulation analysis

The virtual and ADMET screening methods were used to choose the top four compounds from a library containing 13336 molecules. Further, molecular dynamics (MD) simulations were conducted for a duration of 100 ns to better understand the binding mechanism and dynamic behavior of these compounds in complex with VEGFR-2 (Figure 3).

FIGURE 3

FIGURE 3. The results of the molecular dynamics study: (A) Time evolution of the backbone of the selected complex; (B) The comparative RMSF values for the selected compounds; (C) The comparative Radius of gyration values for the selected compounds; (D–H) The comparative hydrogen bonds and pairs within 0.35 nm for the target protein with the Regorafenib, EANPDB 76, EANPDB 252, NANPDB 4577, NANPDB 4580 respectively.

Stability of the docked complexes

Molecular dynamics simulation is a method that could provide a real-time assessment of receptor-ligand interactions. Unlike molecular dynamics simulation, which can show how ligand binding evolves over time in terms of conformational changes, molecular docking can only represent a single position of the ligand-protein complex (Almehmadi et al., 2022). Therefore, we simulated the complexes of VEGFR-2-Regorafenib, EANPDB 76, EANPDB 252, NANPDB 4577, and NANPDB 4580 for 100 ns using molecular dynamics to learn more about the binding interactions between these molecules. These compounds were identified as the top molecules based on the docking results. In the first 30 ns, the four possible VEGFR-2 inhibitors (EANPDB 76, EANPDB 252, NANPDB 4577, and NANPDB 4580) and the reference ligand (Regorafenib) exhibited backbone RMSD values between 0.1 and 0.4 nm (Figure 3A), suggesting that these molecules would adapt to a new conformation within the binding pocket (Mahmood et al., 2022). Subsequently, all systems except EANPDB 76 clearly reached the plateau after 30 ns, having shown a deflection of 0.19 nm and 0.15 nm at 65 ns in the case of EANPDB 252 and NANDB 4580, respectively. Compared with the Regorafenib-VEGFR-2 system, the complexes EANPDB 252, NANPDB 4577, and NANPDB 4580-VEGFR-2 demonstrated better stability. This suggests that these three derivatives may have a higher binding strength with VEGFR-2 for VEGFR-2 than Regorafenib.

Root mean square fluctuation (RMSF) analysis

The root mean square fluctuation (RMSF) of the molecular dynamics trajectory of a protein is used to characterize changes in residue flexibility (Tabti et al., 2022). Figure 3B shows that the RMSF profiles of the four systems were very similar. Consistent with previous findings, the Lys941-Met1016 residues were the most flexible, suggesting that they play a minor role in receptor-ligand interactions. Among these residues, the region containing Thr940-Glu989 residues were reported not to significantly affect VEGFR-2 catalytic activity (McTigue et al., 1999). Moreover, the Stability around Ala866-Lys941 and Met1016-Arg1050 suggested that these residues were important for binding of the potential inhibitors to VEGFR-2. Both areas represented stable hydrogen interaction zones (Zhang et al., 2013). They contained important residues such as Cys917, Cys919, and Asp1046. These residues formed hydrogen bonds with the inhibitors and reduced the flexibility of the corresponding domains. In this study, we discovered that the three natural derivatives formed hydrogen bonds with Asp1046, Cys919, or Lys868, leading to identical RMSF profiles. Indeed, our findings showed that EANPDB 252, NANPDB 4577, and NANPDB 4580 were more stable than the Regorafenib-VEGFR-2 complex, whereas the EANPDB 76-VEGFR-2 complex was more flexible.

Radius of gyration analysis

We examined how the compactness of the protein’s structure is altered by binding to different compounds. For this purpose, we calculated the radius of gyration (Rg) as a function of time (Figure 3C). Consisten with the results of RMSD analysis of the protein backbone, EANPDB 252, NANPDB 4577, and NANPDB 4580 exhibited a higher degree of protein compactness than Regorafenib. In addition, the systems containing FDA-approved inhibitors (Regorafenib) showed a relatively stable protein compactness profile. The VEGFR-2-NANPDB 4577 system, which was the most compact of all systems, had an average value of 1.98 nm, whereas the VEGFR-2- NANPDB 4580, VEGFR-2 - EANPDB 252 complexes had an average value of 2.00 nm. The VEGFR-2-ANDB-76 system had an exceptional profile (2.02 nm).

Hydrogen bonds analysis

Both the affinity for the ligand and its stability depends on the H-bonds that take place between the protein and the ligand, and it increases in direct proportion to the number of H-bonds present in the system (Chen et al., 2016). The Gromacs gmx_hbond module was used to estimate the number of H-bonds formed between the protein and its ligands during 100 ns. The results of this estimation are shown in Figures 3D–G. Regorafenib was found to form an average of 2.01 hydrogen bonds and 3.23 bond pairs within 0.35 nm of the active pocket of VEGFR-2. Similarly, NANPDB 4577 was linked to VEGFR-2 in the binding site via an average of 3.97 hydrogen bonds, while the average number of pairs within 0.35 nm was 4.15, followed by VEGFR-2 and NANPDB 4580, VEGFR-2 and NANPDB 252 with an average number of hydrogen bonds of 1.5 and 1, and an average number of pairs within 0.35 nm of 2.1 and 3.2, respectively. However, for VEGFR-2/NANPDB 4580, the average number of hydrogen bonds was 2.48 and the average number of pairs within 0.35 nm was 3.45. Although the latter two compounds formed fewer hydrogen bonds compared to the reference compound (Regorafenib), they still created stable interactions in the active site of VEGFR-2 by forming hydrogen bonds with crucial residues. The four derivative compounds (EANPDB 252, NANPDB 4577, and NANPDB 4580) were found to be equivalent to Regorafenib and were able to effectively target VEGFR-2 against angiogenesis, as shown by the results of molecular dynamics simulations, despite minor differences between these inhibitors.

MM-PBSA binding free energy analysis

The binding free energy calculations were performed using the Python script MmPbSaStat.py from the package g_mmpbsa to calculate the average binding free energy for the entire trajectory and to further investigate the interactions of the selected complexes (Jayaraman et al., 2021) (Table 3). Using the output files from g_mmpbsa, this script determines the average free energy of binding and the associated standard deviation/error.

TABLE 3

TABLE 3. MMPBSA calculations of the predicted binding free energy for all complexes.

Binding energy indicates how much energy is released when a bond is formed or when a ligand and a protein interact. Several types of energy were considered when calculating binding energy, namely Van der Waal energy, Electrostatic Energy, Solvent Accessible Surface Energy (SASA), and the Polar solvation energy, which was the only one excluded from the analysis (Bhardwaj et al., 2020). All other energy types had a positive effect on the interaction between VEGFR-2 and various molecules. Of all the identified compounds, NANPDB 4577 showed the lowest binding free energy (−161.152 kJ/mol). This value represents the energy required to bind two molecules together. Regorafenib, a newly developed drug, had the second lowest binding free energy (−148.036 kJ/mol). The connections between VEGFR-2 and three putative anti-angiogenic compounds were found to be mediated mainly by Van der Waals interactions (Evdw) rather than electrostatic interactions (Eele). All three derivatives contributed as strongly to Van der Waals forces as regorafenib, but their electrostatic contributions were much weaker in the case of NANPDB 4577 and EANDB 252. It is important to mention that the MMPBSA energetics has two major issues related to the intrinsic dielectric constant and the determination of the entropic term efficiently and accurately. Based on this, other methods such as experimental assays can be conducted in order to validate the accuracy and reliability of the computational predictions.

Conclusion

In the present study, a virtual screening method was used to search for potential anticancer drugs against VEGFR-2 among 13313 natural products in the African database. Using molecular docking calculations, and ADMET properties, 4 compounds were identified that, when bound to VEFGR-2, had a binding affinity −11.0 kcal/mol and favorable physiochemical and pharmacokinetic properties. These 4 compounds were then subjected to 100 ns long MD simulations to determine their stability and their binding energy using the MM-PBSA method. Compared to Regorafenib as an approved drug, root means square deviations, root means square fluctuation, radius of gyration, hydrogen bonds, and binding-free energy analysis showed that, with the exception of EANPDB 76, all candidatess, including EANPDB 252, NANPDB 4577, and NANPDB 4580, which are derived from natural sources, can effectively target VEGFR-2, similar to Regorafenib. Therefore, we suggest three substances that can be used in anti-angiogenesis treatment. These compounds are expected to deactivate VEGFR-2 and thus reduce the process of angiogenesis. However, it should be noted that their suitability for clinical use and safety need to be further explored, because in vitro and in vivo experiments were not performed in our study.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.

Author contributions

RD, AA, and SB designed the study; the initial draft of the article was done by SB; Formal analysis and editing of the manuscript were done by all authors. Supervision was done by RD and AA. All authors contributed to the article and approved the submitted version.

Acknowledgments

The authors acknowledge the African Supercomputing Center at Mohamed VI Polytechnic University for the supercomputing resources (https://ascc.um6p.ma/) made available for conducting the research reported in this paper.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Abdelkader, A., Elzemrany, A. A., El-Nadi, M., Elsabbagh, S. A., Shehata, M. A., Eldehna, W. M., et al. (2022). In-silico targeting of SARS-CoV-2 NSP6 for drug and natural products repurposing. Virology 573, 96–110. doi:10.1016/J.VIROL.2022.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Al-Khafaji, K., and Taskin Tok, T. (2020). Molecular dynamics simulation, free energy landscape and binding free energy computations in exploration the anti-invasive activity of amygdalin against metastasis. Comput. Methods Programs Biomed. 195, 105660. doi:10.1016/J.CMPB.2020.105660

PubMed Abstract | CrossRef Full Text | Google Scholar

Aljuaid, A., Salam, A., Almehmadi, M., Baammi, S., Alshabrmi, F. M., Allahyani, M., et al. (2022). Structural Homology-based drug repurposing approach for targeting NSP12 SARS-CoV-2. Molecules 27, 7732. doi:10.3390/MOLECULES27227732

PubMed Abstract | CrossRef Full Text | Google Scholar

Almehmadi, M., Allahyani, M., Aljuaid, A., Muhammad, T., Baammi, S., et al. (2022). Identification of human immunodeficiency virus-1 E protein-targeting lead compounds by pharmacophore based screening. Article Saudi Med. J. 43, 1324–1332. doi:10.15537/smj.2022.43.12.20220599

PubMed Abstract | CrossRef Full Text | Google Scholar

Ansari, M. J., Bokov, D., Markov, A., Jalil, A. T., Shalaby, M. N., Suksatan, W., et al. (2022). Cancer combination therapies by angiogenesis inhibitors; a comprehensive review. Cell Commun. Signal. 20, 49–23. doi:10.1186/s12964-022-00838-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Baammi, S., Daoud, R., and El Allali, A. (2022). Assessing the effect of a series of mutations on the dynamic behavior of phosphite dehydrogenase using molecular docking, molecular dynamics and quantum mechanics/molecular mechanics simulations. J. Biomol. Struct. Dyn. 41, 4154–4166. doi:10.1080/07391102.2022.2064912

PubMed Abstract | CrossRef Full Text | Google Scholar

Baammi, S., Daoud, R., and El Allali, A. (2023). In silico protein engineering shows that novel mutations affecting NAD+ binding sites may improve phosphite dehydrogenase stability and activity. Sci. Rep. 13, 1878. doi:10.1038/S41598-023-28246-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Beutler, J. A. (2009). Natural products as a Foundation for drug discovery. Curr. Protoc. Pharmacol. 46, 9.11.1–9. doi:10.1002/0471141755.PH0911S46

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhardwaj, V. K., Singh, R., Sharma, J., Rajendran, V., Purohit, R., and Kumar, S. (2020). Identification of bioactive molecules from tea plant as SARS-CoV-2 main protease inhibitors. J. Biomol. Struct. Dyn. 39, 3449–3458. doi:10.1080/07391102.2020.1766572

PubMed Abstract | CrossRef Full Text | Google Scholar

Byrne, A. M., Bouchier-Hayes, D. J., and Harmey, J. H. (2005). Angiogenic and cell survival functions of vascular endothelial growth factor (VEGF). J. Cell Mol. Med. 9, 777–794. doi:10.1111/J.1582-4934.2005.TB00379.X

PubMed Abstract | CrossRef Full Text | Google Scholar

Camero, C. M., Germanò, M. P., Rapisarda, A., D’Angelo, V., Amira, S., Benchikh, F., et al. (2018). Anti-angiogenic activity of iridoids from Galium tunetanum. Rev. Bras. Farmacogn. 28, 374–377. doi:10.1016/J.BJP.2018.03.010

CrossRef Full Text | Google Scholar

Campos-Xolalpa, N., Pérez-Ramos, J., Esquivel-Campos, A., Pérez-González, C., Sánchez-Pérez, L., Pérez-Gutiérrez, S., et al. (2021). Cytotoxic and antimicrobial activities of quinones isolated from different organism. Cytotox. - New Insights into Toxic. Assess. doi:10.5772/INTECHOPEN.95598

CrossRef Full Text | Google Scholar

Carmeliet, P., and Jain, R. K. (2011). Molecular mechanisms and clinical applications of angiogenesis. Nature 473, 298–307. doi:10.1038/NATURE10144

PubMed Abstract | CrossRef Full Text | Google Scholar

Carmeliet, P. (2005). VEGF as a key mediator of angiogenesis in cancer. Oncology 69 (3), 4–10. doi:10.1159/000088478

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, D., Oezguen, N., Urvil, P., Ferguson, C., Dann, S. M., and Savidge, T. C. (2016). Regulation of protein-ligand binding affinity by hydrogen bond pairing. Sci. Adv. 2, e1501240. doi:10.1126/SCIADV.1501240

PubMed Abstract | CrossRef Full Text | Google Scholar

Devaurs, D., Antunes, D. A., Hall-Swan, S., Mitchell, N., Moll, M., Lizée, G., et al. (2019). Using parallelized incremental meta-docking can solve the conformational sampling issue when docking large ligands to proteins. BMC Mol. Cell Biol. 20, 42. doi:10.1186/S12860-019-0218-Z

PubMed Abstract | CrossRef Full Text | Google Scholar

Dias, D. A., Urban, S., and Roessner, U. (2012). A Historical overview of natural products in drug discovery. Metabolites 2, 303–336. doi:10.3390/METABO2020303

PubMed Abstract | CrossRef Full Text | Google Scholar

En-nahli, F., Baammi, S., Hajji, H., Alaqarbeh, M., Lakhlifi, T., and Bouachrine, M. (2022). High-throughput virtual screening approach of natural compounds as target inhibitors of plasmepsin-II. J. Biomol. Struct. Dyn., 1–11. doi:10.1080/07391102.2022.2152871

CrossRef Full Text | Google Scholar

Guan, L., Yang, H., Cai, Y., Sun, L., Di, P., Li, W., et al. (2019). ADMET-Score – A comprehensive scoring function for evaluation of chemical drug-likeness. Medchemcomm 10, 148–157. doi:10.1039/C8MD00472B

PubMed Abstract | CrossRef Full Text | Google Scholar

Hosseini, F. S., and Amanlou, M. (2020). Anti-HCV and anti-malaria agent, potential candidates to repurpose for coronavirus infection: virtual screening, molecular docking, and molecular dynamics simulation study. Life Sci. 258, 118205. doi:10.1016/J.LFS.2020.118205

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivorra, M. D., Payá, M., and Villar, A. (1989). A review of natural products and plants as potential antidiabetic drugs. J. Ethnopharmacol. 27, 243–275. doi:10.1016/0378-8741(89)90001-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Jain, A. N. (2004). Virtual screening in lead discovery and optimization. Curr. Opin. Drug Discov. Devel 7, 396–403. (accessed May 10, 2023). Available at: https://europepmc.org/article/med/15338948.

PubMed Abstract | Google Scholar

Jayaraman, M., Loganathan, L., Muthusamy, K., and Ramadas, K. (2021). Virtual screening assisted discovery of novel natural products to inhibit the catalytic mechanism of Mycobacterium tuberculosis InhA. J. Mol. Liq. 335, 116204. doi:10.1016/J.MOLLIQ.2021.116204

CrossRef Full Text | Google Scholar

Kumari, R., Kumar, R., and Lynn, A. (2014). g_mmpbsa—a GROMACS tool for high-throughput MM-PBSA calculations. J. Chem. Inf. Model 54, 1951–1962. doi:10.1021/ci500020m

PubMed Abstract | CrossRef Full Text | Google Scholar

Lagorce, D., Douguet, D., Miteva, M. A., and Villoutreix, B. O. (2017). Computational analysis of calculated physicochemical and ADMET properties of protein-protein interaction inhibitors. Sci. Rep. 7, 46277–46315. doi:10.1038/srep46277

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, W. L., Zheng, H. C., Bukuru, J., and De Kimpe, N. (2004). Natural medicines used in the traditional Chinese medical system for therapy of diabetes mellitus. J. Ethnopharmacol. 92, 1–21. doi:10.1016/j.jep.2003.12.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindahl, E., Bjelkmar, P., Larsson, P., Cuendet, M. A., and Hess, B. (2010). Implementation of the CHARMM force field in GROMACS: analysis of protein stability effects from Correction Maps, virtual interaction sites, and water models. J. Chem. Theory Comput. 6, 459–466. doi:10.1021/CT900549R

PubMed Abstract | CrossRef Full Text | Google Scholar

Lourenço, S. C., Moldão-Martins, M., and Alves, V. D. (2019). Antioxidants of natural plant Origins: from sources to Food Industry applications. Molecules 24, 4132. doi:10.3390/MOLECULES24224132

PubMed Abstract | CrossRef Full Text | Google Scholar

Mahmood, R. A., Hasan, A., Rahmatullah, M., Paul, A. K., Jahan, R., Jannat, K., et al. (2022). Solanaceae Family phytochemicals as inhibitors of 3C-like protease of SARS-CoV-2: an in silico analysis. Molecules 27, 4739. doi:10.3390/MOLECULES27154739

PubMed Abstract | CrossRef Full Text | Google Scholar

Mandlik, V., Bejugam, P. R., and Singh, S. (2016). Application of Artificial neural Networks in modern drug discovery. Artif. Neural Netw. Drug Des. Deliv. Dispos., 123–139. doi:10.1016/B978-0-12-801559-9.00006-5

CrossRef Full Text | Google Scholar

McInnes, C. (2007). Virtual screening strategies in drug discovery. Curr. Opin. Chem. Biol. 11, 494–502. doi:10.1016/J.CBPA.2007.08.033

PubMed Abstract | CrossRef Full Text | Google Scholar

McTigue, M. A., Wickersham, J. A., Pinko, C., Showalter, R. E., Parast V, C., Tempczyk-Russell, A., et al. (1999). Crystal structure of the kinase domain of human vascular endothelial growth factor receptor 2: A key enzyme in angiogenesis. Structure 7, 319–330. doi:10.1016/S0969-2126(99)80042-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Motzer, R. J., Rini, B. I., Bukowski, R. M., Curti, B. D., George, D. J., Hudes, G. R., et al. (2006). Sunitinib in patients with metastatic renal cell carcinoma. JAMA 295, 2516–2524. doi:10.1001/jama.295.21.2516

PubMed Abstract | CrossRef Full Text | Google Scholar

Nieves, B. J., D’Amore, P. A., and Bryan, B. A. (2009). The function of vascular endothelial growth factor. Biofactors 35, 332–337. doi:10.1002/BIOF.46

PubMed Abstract | CrossRef Full Text | Google Scholar

Nishida, N., Yano, H., Nishida, T., Kamura, T., and Kojiro, M. (2006). Angiogenesis in cancer. Vasc. Health Risk Manag. 2, 213–219. doi:10.2147/VHRM.2006.2.3.213

PubMed Abstract | CrossRef Full Text | Google Scholar

Niu, G., and Chen, X. (2010). Vascular endothelial growth factor as an anti-angiogenic target for cancer therapy. Curr. Drug Targets 11, 1000–1017. doi:10.2174/138945010791591395

PubMed Abstract | CrossRef Full Text | Google Scholar

Ntie-Kang, F., Telukunta, K. K., Döring, K., Simoben, C. V., Moumbock, A. F. A., Malange, Y. I., et al. (2017). Nanpdb: A resource for natural products from Northern African sources. J. Nat. Prod. 80, 2067–2076. doi:10.1021/acs.jnatprod.7b00283

PubMed Abstract | CrossRef Full Text | Google Scholar

Olajide, O. A., and Sarker, S. D. (2020). Anti-inflammatory natural products. Annu. Rep. Med. Chem. 55, 153–177. doi:10.1016/BS.ARMC.2020.02.002

CrossRef Full Text | Google Scholar

Park, J. Y., Lee, Y., Lee, H. J., Kwon, Y. S., and Chun, W. (2020). In silico screening of GABA aminotransferase inhibitors from the constituents of Valeriana officinalis by molecular docking and molecular dynamics simulation study. J. Mol. Model 26, 228–313. doi:10.1007/s00894-020-04495-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Patel, K., Kumar, V., Rahman, M., Verma, A., and Patel, D. K. (2018). New insights into the medicinal importance, physiological functions and bioanalytical aspects of an important bioactive compound of foods ‘Hyperin’: health benefits of the past, the present, the future. Beni Suef Univ. J. Basic Appl. Sci. 7, 31–42. doi:10.1016/J.BJBAS.2017.05.009

CrossRef Full Text | Google Scholar

Pires, D. E. V., Blundell, T. L., and Ascher, D. B. (2015). pkCSM: predicting small-molecule pharmacokinetic and toxicity properties using Graph-based Signatures. J. Med. Chem. 58, 4066–4072. doi:10.1021/ACS.JMEDCHEM.5B00104

PubMed Abstract | CrossRef Full Text | Google Scholar

Prakash, O., Kumar, A., Kumar, P., and Ajeet, A. (2013). Anticancer potential of plants and natural products: A review. Am. J. Pharmacol. Sci. 1, 104–115. doi:10.12691/AJPS-1-6-1

CrossRef Full Text | Google Scholar

Rajabi, M., and Mousa, S. A. (2017). The role of angiogenesis in cancer treatment. Biomedicines 5, 34. doi:10.3390/BIOMEDICINES5020034

PubMed Abstract | CrossRef Full Text | Google Scholar

Rizvi, S. A. A., Einstein, G. P., Tulp, O. L., Sainvil, F., and Branly, R. (2022). Introduction to traditional medicine and their role in prevention and treatment of Emerging and Re-Emerging diseases. Biomolecules 12, 1442. doi:10.3390/BIOM12101442

PubMed Abstract | CrossRef Full Text | Google Scholar

Salifu, E. Y., Abugri, J., Rashid, I. A., Osei, F., and Ayariga, J. A. (2023). In silico identification of potential inhibitors of acyl carrier protein reductase and acetyl CoA carboxylase of Plasmodium falciparum in antimalarial therapy. Front. Drug Discov. 3, 1. doi:10.3389/FDDSV.2023.1087008

CrossRef Full Text | Google Scholar

Sasidharan, S., Chen, Y., Saravanan, D., Sundram, K. M., and Yoga Latha, L. (2011). Extraction, isolation and Characterization of bioactive compounds from plants’ Extracts. Afr. J. Traditional, Complementary, Altern. Med. 8, 1–10. doi:10.4314/ajtcam.v8i1.60483

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmieder, R., Hoffmann, J., Becker, M., Bhargava, A., Müller, T., Kahmann, N., et al. (2014). Regorafenib (BAY 73-4506): antitumor and antimetastatic activities in preclinical models of colorectal cancer. J. Int. Du Cancer 135, 1487–1496. doi:10.1002/IJC.28669

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, S., Johnson, D., Abouammoh, M., Hollands, S., and Brissette, A. (2012). Rate of serious adverse effects in a series of bevacizumab and ranibizumab injections. Can. J. Ophthalmol. 47, 275–279. doi:10.1016/J.JCJO.2012.03.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Shoichet, B. K. (2004). Virtual screening of chemical libraries. Nature 432, 862–865. doi:10.1038/NATURE03197

PubMed Abstract | CrossRef Full Text | Google Scholar

Simoben, C. V., Qaseem, A., Moumbock, A. F. A., Telukunta, K. K., Günther, S., Sippl, W., et al. (2020). Pharmacoinformatic investigation of medicinal plants from East Africa. Mol. Inf. 39, 2000163. doi:10.1002/MINF.202000163

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, S. (2006). Preclinical pharmacokinetics: an approach towards safer and efficacious drugs. Curr. Drug Metab. 7, 165–182. doi:10.2174/138920006775541552

PubMed Abstract | CrossRef Full Text | Google Scholar

Tabti, K., Baammi, S., ElMchichi, L., Sbai, A., Maghat, H., Bouachrine, M., et al. (2022). Computational investigation of pyrrolidin derivatives as novel GPX4/MDM2–p53 inhibitors using 2D/3D-QSAR, ADME/toxicity, molecular docking, molecular dynamics simulations, and MM-GBSA free energy. Struct. Chem. 33, 1019–1039. doi:10.1007/s11224-022-01903-5

CrossRef Full Text | Google Scholar

Tabti, K., Baammi, S., Sbai, A., Maghat, H., Lakhlifi, T., and Bouachrine, M. (2023). Molecular modeling study of pyrrolidine derivatives as novel myeloid cell leukemia-1 inhibitors through combined 3D-QSAR, molecular docking, ADME/Tox and MD simulation techniques. ADME/Tox MD Simul. Tech., 1–17. doi:10.1080/07391102.2023.2183032

CrossRef Full Text | Google Scholar

Thakur, M., Singh, K., and Khedkar, R. (2020). Phytochemicals: extraction process, safety assessment, toxicological evaluations, and regulatory issues. Funct. Preserv. Prop. Phytochemicals, 341–361. doi:10.1016/B978-0-12-818593-3.00011-7

CrossRef Full Text | Google Scholar

Walker, D. K. (2004). The use of pharmacokinetic and pharmacodynamic data in the assessment of drug safety in early drug development. Br. J. Clin. Pharmacol. 58, 601–608. doi:10.1111/J.1365-2125.2004.02194.X

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, D., and Meroueh, S. O. (2016). Effect of binding pose and modeled structures on SVMGen and GlideScore Enrichment of chemical libraries. J. Chem. Inf. Model 56, 1139–1151. doi:10.1021/ACS.JCIM.5B00709

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoo, S. Y., and Kwon, S. M. (2013). Angiogenesis and its therapeutic opportunities. Mediat. Inflamm. 2013, 127170. doi:10.1155/2013/127170

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, W., and Mackerell, A. D. (2017). Computer-aided drug design methods. Methods Mol. Biol. 1520, 85–106. doi:10.1007/978-1-4939-6634-9_5

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Yang, S., Jiao, Y., Liu, H., Yuan, H., Lu, S., et al. (2013). An integrated virtual screening approach for VEGFR-2 inhibitors. J. Chem. Inf. Model 53, 3163–3177. doi:10.1021/CI400429G

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, K., Yuan, Y., Lin, B., Miao, Z., Li, Z., Guo, Q., et al. (2018). LW-215, a newly synthesized flavonoid, exhibits potent anti-angiogenic activity in vitro and in vivo. Gene 642, 533–541. doi:10.1016/J.GENE.2017.11.065

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, L., Yuan, X., Wang, J., Feng, Y., Ji, F., Li, Z., et al. (2019). A review on flavones targeting serine/threonine protein kinases for potential anticancer drugs. Bioorg Med. Chem. 27, 677–685. doi:10.1016/J.BMC.2019.01.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Zoete, V., Cuendet, M. A., Grosdidier, A., and Michielin, O. (2011). SwissParam: A fast force field generation tool for small organic molecules. J. Comput. Chem. 32, 2359–2368. doi:10.1002/JCC.21816

PubMed Abstract | CrossRef Full Text | Google Scholar