Structure-based molecular screening and dynamic simulation of phytocompounds targeting VEGFR-2: a novel therapeutic approach for papillary thyroid carcinoma

Papillary thyroid carcinoma (PTC) is the most prevalent type of thyroid cancer, with aggressive variants presenting major therapeutic challenges. Vascular endothelial growth factor receptor-2 (VEGFR-2) is a key regulator of tumor angiogenesis and is highly expressed in PTC, making it a promising target for therapeutic intervention. This highlights the potential of VEGFR-2 inhibition as an effective strategy for managing PTC. In this study, we employed virtual drug screening, molecular dynamics simulations, and binding free energy calculations to identify potential VEGFR-2 inhibitors from the African natural product database (AfroDb). Our virtual drug screening identified three lead compounds SA_0090, 17.3.1.7.8 and BMC_0005 with a docking scores of −9.04 kcal/mol, −8.96 kcal/mol, and −8.33 kcal/mol respectively, surpassing the control compound (−8.39 kcal/mol). Molecular dynamics simulation analysis confirmed the dynamic stability, structural compactness, and minimal residual fluctuations of the 17.3.1.7.8 and BMC_0005 compounds-VEGFR2 complexes. The binding free energy calculations further supported the strong interactions, with values recorded as −60.3861 ± 0.39 kcal/mol for the control, −52.2732 ± 0.37 kcal/mol for SA_0090, −52.7797 ± 0.62 kcal/mol for 17.3.1.7.8, and −61.476 ± 0.59 kcal/mol for BMC_0005. Additionally, the selected compounds exhibited highly favorable ADMET properties, including optimal water solubility, efficient gastrointestinal absorption, and a non-hepatotoxic profile, all aligning with Lipinski’s rule of five. In conclusion, these findings highlight 17.3.1.7.8 and BMC_0005 compounds as compelling candidates for VEGFR-2 inhibition, offering a promising therapeutic avenue for papillary thyroid carcinoma, warranting further in vitro and in vivo validation for potential therapeutic use.

Introduction

Thyroid carcinoma is the most prevalent endocrine malignancy, with Papillary Thyroid Carcinoma (PTC) being the most common subtype, representing 80%–85% of thyroid cancers in both adults and children (Limaiem et al., 2024). Originating from the follicular cells of the thyroid, PTC is typically well-differentiated and has a low, stable mortality rate. However, certain aggressive variants and advanced stages pose a greater risk, often leading to recurrence, metastasis, and poor outcomes (Gordon et al., 2022) (Mao et al., 2025). Notably, PTC with multifocality and early lymph node metastasis shows a recurrence rate of up to 35%, and the 10-year survival rate for advanced cases drops below 50% (Tuttle et al., 2017). Given PTC’s high incidence and the aggressive behavior of some forms, understanding the risk factors, molecular mechanisms, and metastatic processes is essential, as current treatment approaches remain insufficient. The exact cause of papillary thyroid carcinoma (PTC) is unknown, but risk factors include genetic syndromes (e.g., familial adenomatous polyposis), radiation exposure, and iodine imbalance (Zhang and Xu, 2024). Environmental and hormonal factors, such as endocrine disruptors and prolonged estrogen exposure, also play a role. PTC may be asymptomatic or present with a painless thyroid nodule, voice changes, or difficulty swallowing (Harahap and Jung, 2024).

Metastasis is most frequent in Hurthle cell tumors (33%) and is also common in medullary and anaplastic thyroid cancers. Overall, distant metastasis occurs in about 4% of thyroid tumors at diagnosis, while skin metastasis from papillary thyroid cancer is rare (<1%) (Somoza et al., 2013). Vasculogenesis, angiogenesis, tumorigenesis, and inflammation involve multiple physiological and pathological processes driven by factors like fibroblast growth factor, VEGF, HGF, and interleukin-6 (Ghalehbandi et al., 2023). Increased angiogenesis supports tumor growth by stimulating new capillary formation from existing blood vessels (Modi and Kulkarni, 2020). Tyrosine kinases (TKs), particularly VEGFR-2, are key regulators of this process and are overexpressed in various cancers. VEGFR-2 activation triggers signaling pathways that enhance cell survival, proliferation, and growth (Sana et al., 2020). Previous studies have revealed that VEGFR-2 is highly expressed in both nodular hyperplasia and papillary thyroid carcinoma (PTC). Its expression is particularly pronounced in PTC, where it is accompanied by the co-expression of VEGF and VEGFR-1. This suggests that targeting VEGFR could be a valuable strategy for managing papillary thyroid carcinoma (Gogiashvili et al., 2024).

The VEGFR-2 gene, located at chromosome 4q11-12, encodes a receptor tyrosine kinase composed of 1,356 amino acids. This receptor exists in three forms: a non-glycosylated version (150 kD), an intermediate glycosylated form (200 kD), and a fully mature glycosylated receptor (230 kD). Notably, only the mature glycosylated VEGFR-2 is capable of initiating intracellular signal transduction (Shah et al., 2025). VEGFR2 (Vascular Endothelial Growth Factor Receptor 2) plays a central role in cancer angiogenesis by mediating the effects of VEGF-A. When VEGF-A binds to VEGFR2, it activates several downstream signaling pathways (like PI3K/AKT, MAPK, and Src pathways) that promote endothelial cell survival, proliferation, migration, vascular permeability, and capillary formation. This leads to the development of an abnormal, leaky, and disorganized tumor vasculature, which supports tumor growth and metastasis (Ghalehbandi et al., 2023). VEGFR-2 tyrosine kinase inhibitors completely blocked VEGF-induced angiogenesis and significantly reduced bFGF-induced angiogenesis in both in vivo and in vitro models. In endothelial cell invasion assays, these inhibitors suppressed VEGF- and bFGF-driven invasion by 100% and about 90%, respectively (Tille et al., 2001; Sarabipour et al., 2016).

VEGFR-2 is overactive in cancer cells but not in normal cells, making it a prime target for selective cancer therapies. Blocking VEGFR-2 activation with inhibitors prevents tumor angiogenesis without affecting healthy tissues. Several FDA-approved VEGFR-2 inhibitors have been developed to treat various cancers by restricting blood vessel growth. However, these drugs often cause significant side effects, leading researchers to explore new small molecules with better efficacy and fewer adverse effects (Claesson-Welsh and Welsh, 2013). Computational chemistry has become a key tool in drug discovery, aiding in the design, optimization, and ADMET evaluation of potential VEGFR-2 inhibitors (Suleimen et al., 2021) (El-Adl et al., 2021) (Alanazi et al., 2021). Natural products (NPs) represent a vast reservoir of chemically diverse and biologically active molecules, many of which have served as essential drugs or lead compounds for treating various diseases (Shahrajabian et al., 2022) (Mehrbod et al., 2021). For centuries, traditional medicine has relied on these naturally derived compounds, demonstrating their therapeutic potential. Compared to synthetic drugs, natural compounds often exhibit greater selectivity, reduced toxicity, and cost-effectiveness, making them attractive candidates for drug development (Aware et al., 2022). African natural products offer unique chemical diversity and structural novelty, often guided by traditional medicinal knowledge. They show strong potential against infectious and drug-resistant diseases and are underrepresented in global libraries, making them a valuable source for novel drug discovery. In this study, we employed virtual drug screening, molecular dynamics simulations, and binding free energy calculations to identify potential VEGFR-2 inhibitors from the African natural product database. Inhibiting VEGFR-2 could serve as an effective therapeutic strategy for managing papillary thyroid carcinoma by suppressing tumor angiogenesis.

Methodology

Retrieval and preparation of crystal structure of VEGFR2

The crystal structure of VEGFR2 bound to a native benzimidazole-urea inhibitor (PDB ID: 2OH4) was retrieved from the Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB PDB) (https://www.rcsb.org/) (Suleman et al., 2025). To prepare the structure for further analysis, all water molecules were eliminated using PyMOL. Subsequently, hydrogen atoms were incorporated into the protein, and energy minimization was performed using Chimera to refine the overall structure and resolve any steric clashes (Zhang et al., 2021; Sayaf et al., 2023).

Virtual screening of natural products libraries against VEGFR2 protein

The AfroDb database (https://zinc12.docking.org/pbcs/afronp), comprising a diverse collection of bioactive natural products derived from African medicinal plants, was retrieved and formatted for compatibility (Ntie-Kang et al., 2013). To ensure drug-likeness and minimize toxicity risks, compounds violating Lipinski’s Rule of Five (R5) were filtered out using the FAF4drug online webserver (Sayaf et al., 2023). Before virtual screening with EasyDock Vina 2.0, ligand structures were converted to the.pdbqt format. Open Babel was employed to assign atomic charges, atom types, and prepare ligands. Receptor preparation involved generating grid maps in AutoGrid, defining the known active site residues (Glu883, Val914, Cys917, Asp1044, and Phe1045) (Elkaeed et al., 2022). Preliminary docking was conducted using the AUTODOCK4 algorithm with an exhaustiveness value of 16 for initial screening, followed by a more refined evaluation at exhaustiveness 64 to eliminate false positives. The top 54 compounds, with docking scores ranging from −7.601 to −8.839 kcal/mol, proceeded to IFD (induced fit docking) via AutoDockFR. This method incorporated receptor flexibility, covalent docking, and using force fields like Amber or CHARMM, along with force-field-based scoring functions. Afterword, the default IFD parameters were applied (Ravindranath et al., 2015). In the rigid docking phase (initial screening), the receptor was kept fixed, and only ligand flexibility was considered. However, In the flexible docking (IFD) stage, key amino acid residues around the binding site were allowed to move (side-chain flexibility), while the rest of the receptor remained relatively rigid. Scoring functions in AutoDockFR were used both before and after the flexible adjustment. Initially, standard AutoDock scoring functions estimated the binding poses. After incorporating receptor flexibility, the force-field-based scoring recalculated the binding energies, providing a refined and more accurate ranking of ligand binding affinities. Finally, the top four candidates were subjected to structural validation and molecular dynamics simulations. PyMOL and Schrödinger Maestro (academic version) were used for visual analysis to further assess binding interactions and stability (Khan et al., 2022).

Molecular dynamic simulation analysis of VEGFR2-ligands complexes

Molecular dynamic simulation of VEGFR2-ligands complexes was performed by using AMBER21 software known for its sophisticated computational techniques and offering great detail for the stepwise or holistic visualization of biomolecular systems. The tLeap module was used to make coordinates and topology files corresponding to each protein-ligand complex (Case et al., 2005; Salomon-Ferrer et al., 2013). The system was solvated in a TIP3P water box (14 Å) and counterions (Na + or Cl−) were added for neutralization. To generate the necessary topology and force field modification (frcmod) files for the ligands, the parameters were assigned according to GAFF2 force field by using the antechamber and parmchk2 tools. By alternating between the conjugate gradient and steepest descent approaches, energy minimization was carried out in phases, enabling the system to achieve a stable conformation while minimizing unfavorable steric conflicts. The next phase included increasing the temperature, which was kept at the setpoint using coupling algorithms like Langevin dynamics or Berendsen thermostat. For the production phase, each system underwent a 400-nanosecond simulation with NPT (constant pressure and temperature) or NVT (constant volume and temperature) ensemble. This stage enabled evaluation of the protein-ligand interactions’ dynamics over a period, which was necessary for understanding binding strength and conformational changes throughout the duration of the interactions (Salomon-Ferrer et al., 2013).

Post-simulation analysis of VEGFR2-ligands complexes

For the post-simulation analysis (residual fluctuation, dynamic stability, compactness, hydrogen bonds) of shortlisted compounds-VEGFR2 complexes, we used the CPPTRAJ and PTRAJ packages (Roe and Cheatham, 2013). At the residue level, we assessed the degree of flexibility using Root Mean Square Fluctuation (RMSF) analysis. Rather than studying the total movement of the complex, the RMSF technique focused on tracking the movement of specific residues throughout a given period of time. The RMSF values were obtained using the following equation:

$ThermalfactororB-factor=\left[\left(8\pi \ast \ast 2\right)/3\right]\left(msf\right)$

The structural compactness of the complexes over the simulation period was assessed by computing the radius of gyration (Rg) using the following mathematical formula:

${\mathit{R}}_{\mathit{g}\mathit{y}\mathit{r}}^2=\frac{1}{\mathit{M}}{\displaystyle \sum _{\mathit{i}=1}^{\mathit{N}}}{\mathit{m}}_{\mathit{i}}\left({\mathit{r}}_{\mathit{i}}-{\mathit{R}}^2\right)$

$\mathit{M}={\displaystyle \sum _{i=1}^N}{m}_i$

${R=N}^{-1}{\displaystyle \sum _{i=1}^N}{r}_i$

Additionally, to determine the stability of the complexes, the Root Mean Square Deviation (RMSD) was calculated. The RMSD values provided insight into the overall structural deviations throughout the simulation, using the mathematical expression below:

$RMSD=\sqrt{\frac{\sum d^{2}i=1}{N_{atoms}}}$

Post-simulation binding free energy calculation of VEGFR2-ligand complexes

To calculate the binding free energies of lead compounds-VEGFR2 complexes we used the MMPBSA. PY script (Sayaf et al., 2023). For the binding free energies calculation, we selected the last 1,000 frames of the MD simulation. Numerous studies use this computational method extensively to assess the Total Binding Energy (TBE) of various ligands (Sayaf et al., 2024) (Khan et al., 2025). The free energy of the ligand in its unbound solvated form (G_{ligand, solvated}), the receptor in its solvated state (G_{receptor, solvated}), and the completely solvated complex (G_{complex, solvated}) were all evaluated. The relationship between these energy components can be further represented through the following equation:

${∆G}_{bind}=G_{\left(complex,solvated\right)}-G_{\left(ligand,solvated\right)}-G_{\left(VEGFR2,solvated\right)}(1)$

To delve deeper into the specific energy contributions, we reformulated Equation as:

$G=E_{MolecularMechanics}-G_{solvated}-TS(2)$

To calculate the specific energy term, the formula was restructured as follows:

${∆G}_{bind}={∆E}_{MolecularMechanics}+{∆G}_{solvated}-∆TS={∆G}_{vaccum}+{∆G}_{solvated}(3)$

${∆E}_{MolecularMechanics}={∆E}_{\mathit{int}}+{∆E}_{electrostatic}+{∆E}_{vdW}(4)$

$∆G_{solvated}=∆G_{Generalizedborn}+∆G_{surfacearea}(5)$

$∆G_{surfacearea}=\gamma .SASA+b(6)$

Lipinski’s rule, and pharmacokinetics analysis

Lipinski’s Rule of Five, which specifies important characteristics for oral bioavailability, is crucial in drug design. These include molecular weight ≤500, hydrogen bond donors ≤5, hydrogen bond acceptors ≤10, and log P ≤5 (Pollastri, 2010). In drug discovery, Lipinski’s Rule of Five is a key rule that aids in predicting the oral bioavailability of possible therapeutic candidates. To check adherence to these computational criteria, we used SwissADME (http://www.swissadme.ch/), a web tool that predicts physicochemical parameters, pharmacokinetics, drug-likeness, and drugability in medicinal chemistry (Daina et al., 2017). Moreover, we conducted an ADMET (Absorption, Distribution, Metabolism, Excretion, and Toxicity) analysis using pkCSM (https://biosig.lab.uq.edu.au/pkcsm/) (Pires et al., 2015), which provided insights into critical pharmacokinetic parameters such as skin sensitivity, hepatotoxicity, solubility, blood-brain barrier permeability, and intestinal absorption. Our goal was to find lead phytocompounds with promising drug-like qualities so that they may be developed further in the pharmaceutical industry by using these computational approaches.

Results and discussion

Papillary Thyroid Carcinoma (PTC) is the most common type of thyroid cancer, where aggressive variants often lead to recurrence with a poor prognosis (Zhang and Xu, 2024). The pro-angiogenic factor, known as VEGFR-2, is overexpressed in a plethora of cancers, including PTC, where it aids in tumor growth. Suppression of tumor angiogenesis by targeting VEGFR-2 makes this approach highly promising (Dhar et al., 2022). FDA-approved inhibitors of VEGFR-2 do exist, but their associated side effects of other such treatment options have become unsatisfactory (Claesson-Welsh and Welsh, 2013). Natural products contain diverse structures and biologically active compounds that can be selectively applied with reduced side effects. Thus, we aim to design putative inhibitors of VEGFR2 using virtual docking, MD simulations, and BFE approaches. This study explored VEGFR-2 inhibitors from the African natural product database to develop effective treatments for PTC. The overall workflow of this study is shown in the Figure 1.

Figure 1

Figure 1. Overall work flow of the study, illustrating the virtual screening and validation of the selected compounds through MD simulation and BFE calculation.

Screening of AfroDb natural product database against VEGFR2 protein

The Virtual drug screening (VDS) is an innovative approach in current drug discovery because it utilizes natural phytocompounds as a prospective source of therapeutics in an economical and efficient manner (Giordano et al., 2022). The unique structural and bioactive characteristics of phytocompounds allow them to be quickly screened for target specific interactions through virtual drug screening techniques. This enables the rapid discovery of novel drug candidates that possess good pharmacokinetics and low toxicity, and decreases the need for time extensive and costly experimental screening. Moreover, VDS offers the possibility of repurposing phytochemicals for several diseases and supports the design of herbal remedies by scientifically confirming their medicinal value. VDS combines herbal products with technology which increases the ability to discover new drugs for cancer, infectious diseases, and other complicated health issues (Oliveira et al., 2023) (Lin et al., 2020). In this study, we conducted virtual screening of the AfroDB database to identify potential inhibitors targeting the active site (Glu883, Val914, Cys917, Asp1044, and Phe1045) of VEGFR2 protein as shown in the Figure 2b. AfroDB is a valuable repository of bioactive compounds derived from African medicinal plants. To refine our selection, we applied Lipinski’s Rule of Five, a widely accepted criterion for assessing drug-like properties, eliminating non-compliant molecules before proceeding with screening (Suleman et al., 2024a). Our computational approach utilized AutoDock Vina for molecular docking against VEGFR2. The initial dataset comprised 954 compounds, which was reduced to 743 after filtering for drug-likeness. These compounds underwent a multi-step screening process. In the primary docking stage, binding affinities ranged from −8.839 to 5.27 kcal/mol, leading us to shortlist 54 compounds with docking scores between −7.601 and −8.839 kcal/mol. These top candidates were then subjected to induced-fit docking, refining the selection further with binding affinities between −7.37 and −9.04 kcal/mol. On the basis of high docking score only three compounds such as 8-oxo-16-[(2R, 3S, 4S, 5S, 6R)-3,4,5-trihydroxy-6-(hydroxymethyl) tetrahydropyran-2-yl]oxy-hexadecanoic, [(1aR, 1bR, 2R, 5aR, 6S, 6aR)-1a-(hydroxymethyl)-2-(2S, 3R, 4R, 5S, 6R)-3,4,5-trihydroxy-6-(hydroxymethyl)te and (2S, 3R, 4R, 5S, 6S)-3,4,5-trihydroxy-6-[2-(4-hydroxyphenyl) ethoxy]tetrahydropyran-2-yl]methyl with a docking scores of −9.04 kcal/mol, −8.96 kcal/mol, and −8.33 kcal/mol were further processed for the bonding network and stability evaluation. The selected compounds with docking scores and binding residues are shown in the Table 1. Furthermore, to validate the docking protocol, a re-docking experiment was conducted using the reference ligand of VEGFR2. This process involved reintroducing the native ligand into the active site, successfully replicating its original binding conformation. The analysis confirmed the reliability of the docking approach by demonstrating a close match with the binding pattern observed in the downloaded ligand structure. The superimposed native ligand on over re-docked is shown in the Figure 2a.

Figure 2

Figure 2. Representation of VEGFR2 active site and validation of docking accuracy. (a) Represents the superimposition of native ligand (Green) over re-docked ligand (orange). (b) Represents the active sites of VEGFR2 and the shortlisted compounds.

Table 1

Table 1. List of lead compounds from African natural compounds and TCM database.

Bonding network analysis of lead compounds-VEGFR2 complexes

Binding network analysis of protein-drug complexes is essential for understanding the molecular interactions that drive drug efficacy, specificity, and stability. By examining non-covalent interactions such as hydrogen bonds, hydrophobic contacts, electrostatic forces, and van der Waals interactions, researchers can identify critical residues involved in binding and optimize drug candidates for improved affinity and selectivity (Klebe, 2025). The binding analysis of the control drug (benzimidazole-urea inhibitor) demonstrated a docking score of −8.39 kcal/mol. This compound formed two hydrogen bonds with Glu883 and Asp1044 residues within the active site (Figure 3a). In contrast, the interaction analysis of the shortlisted compound SA_0090 revealed a stronger binding affinity with a docking score of −9.04 kcal/mol. SA_0090 established five hydrogen bonds with key amino acid residues, including Cys917, Asn921, His1024, Asp1044, and Arg1049 (Figure 3b). Notably, two of these residues, Cys917 and Asp1044, have previously been identified as critical for ligand binding (Elkaeed et al., 2022), suggesting that SA_0090 interacts with essential regions of the active site. Additionally, the increased number of hydrogen bonds indicates more extensive molecular interactions, which may contribute to greater binding stability and specificity. In conclusion, the higher docking score and increased number of hydrogen bond interactions suggest that SA_0090 exhibits stronger and more specific binding to the target protein compared to the control drug.

Figure 3

Figure 3. Docking analysis of control and SA_0090-VEGFR2 complexes. (a) Represents the interaction (2D and 3D) of control and VEGFR2 protein (b) Represents the interaction (2D and 3D) of SA_0090 and VEGFR2 protein.

The docking analysis of the 17.3.1.7.8-VEGFR2 complex revealed a docking score of −8.96 kcal/mol, indicating a stronger binding affinity compared to the control drug (−8.39 kcal/mol). This compound formed four hydrogen bonds with key amino acid residues within the active site of the VEGFR2 protein, specifically Glu883, Lys866, Val912, and Asp1044. Notably, three of these residues Glu883, Val912, and Asp1044 were previously identified as critical drug targets, reinforcing the biological relevance of these interactions. The higher docking score and the presence of key hydrogen bonds suggest that the 17.3.1.7.8 compound exhibits enhanced binding affinity and specificity relative to the control drug (Figure 4a). In comparison, the BMC_0005-VEGFR2 complex showed a slightly lower docking score of −8.33 kcal/mol, forming five hydrogen bonds with essential residues, including Glu883, Lys866, Ile1023, Asp1044, and Phe1045 (Figure 4b). The interaction of BMC_0005 with these key amino acids suggests effective binding; however, the slightly lower docking score indicates a reduced binding affinity compared to the 17.3.1.7.8 compound. In conclusion, the compounds SA_0090 and 17.3.1.7.8 demonstrated stronger binding affinity and increased molecular interactions with VEGFR2 compared to the control drug. These findings suggest their potential as more effective VEGFR2 inhibitors.

Figure 4

Figure 4. Docking analysis of 17.3.1.7.8 and BMC_0005-VEGFR2 complexes. (a) Represents the interaction (2D and 3D) of 17.3.1.7.8 and VEGFR2 protein (b) Represents the interaction (2D and 3D) of BMC_0005 and VEGFR2 protein.

Lipinski’s rule of five evaluation for the selected compounds

Lipinski’s Rule of Five (Ro5) is a crucial guideline in drug discovery used to evaluate the drug-likeness and oral bioavailability of compounds (Nhlapho et al., 2024). It helps predict whether a molecule can be efficiently absorbed by the human body based on key physicochemical properties: molecular weight (≤500 Da), lipophilicity (LogP ≤5), hydrogen bond donors (≤5), and hydrogen bond acceptors (≤10). By applying Ro5 early in drug development, researchers can optimize compounds for better pharmacokinetics, reduce failures in later clinical trials, and enhance the efficiency of medicinal chemistry efforts (Nhlapho et al., 2024). Therefore, to check the physiochemical properties of our selected compounds we performed the Lipinski’s rule of five evaluation for each compound. As shown in the Table 2 all compounds except the control meet the molecular weight criterion (≤500 Da), with the control slightly exceeding it (503.412 Da). “SA_0090″and “BMC_0005″fully comply with Ro5, having zero violations, making them strong drug candidates. “17.3.1.7.8″has one violation due to an excess of HBAs (11 instead of ≤10), which may impact solubility. The control compound has two violations, exceeding the limits for both molecular weight and Log P (6.3354), making it more lipophilic and potentially less soluble. The bioavailability scores range from 0.11 to 0.55, with “SA_0090″showing the lowest (0.11) due to its high polarity, while “17.3.1.7.8″and “BMC_0005″have 0.55, suggesting better oral absorption.

Table 2

Table 2. Lipinski’s rule five analysis for all selected top hits.

ADMET properties (absorption, distribution, metabolism, excretion, toxicity) analysis of selected compounds

ADMET (Absorption, Distribution, Metabolism, Excretion, and Toxicity) evaluation is crucial in drug discovery as it ensures the selection of compounds with optimal pharmacokinetic and safety profiles. It helps enhance bioavailability, predict toxicity risks, and improve drug distribution, reducing failures in clinical trials. Early ADMET screening aids in eliminating unsuitable candidates, optimizing molecular structures, and ensuring compliance with regulatory guidelines, ultimately saving time and costs. By integrating computational and experimental ADMET assessments, researchers can develop safer and more effective drugs, increasing their chances of success in therapeutic applications. As shown in the Table 3 the ADMET evaluation of four compounds (control, SA_0090, 17.3.1.7.8, and BMC_0005) reveals distinct pharmacokinetic and toxicological profiles. In terms of absorption, all compounds have moderate water solubility (Log S values between −2.886 and −3.51), but their Caco-2 permeability varies significantly, with SA_0090 showing the lowest (−0.278) and 17.3.1.7.8 the highest (0.368). Caco-2 cells, derived from human colorectal adenocarcinoma, serve as a standard model for assessing the permeability of substances across the intestinal epithelium (Kus et al., 2023). Furthermore, the control exhibits the highest intestinal absorption (88.506%), while SA_0090 has the lowest (25.89%), with 17.3.1.7.8 and BMC_0005 showing moderate absorption (46.488% and 49.305%, respectively). In terms of distribution, none of the compounds cross the blood-brain barrier, and all localize in mitochondria. Volume of distribution (VDss) suggests SA_0090 has limited tissue distribution (−1.26), whereas BMC_0005 (0.554) distributes more widely. The volume of distribution (VD) measures how a drug disperses between plasma and tissues. If the VD value is lower than −0.15, the drug is more likely to remain in plasma rather than being distributed into tissues. Conversely, a VD value exceeding 0.45 indicates a broader distribution across tissues. Metabolically, none of the compounds act as CYP2D6 or CYP3A4 substrates, but the control inhibits CYP1A2 and CYP2C19, potentially leading to drug-drug interactions. Excretion data show SA_0090 has the highest clearance (1.915), meaning it is rapidly eliminated, whereas BMC_0005 has the lowest clearance (0.127), suggesting longer retention. None of the compounds are substrates for renal OCT2 transporters. Toxicity assessment indicates that all compounds are free from AMES toxicity, skin sensitization, carcinogenicity, and respiratory toxicity; however, the control is hepatotoxic, whereas the other three compounds are not. Overall, our selected compounds exhibit favorable pharmacokinetics properties making them promising alternatives to the control drug.

Table 3

Table 3. Evaluation of ADMET properties of selected compounds.

Post-simulation stability analysis of shortlisted compounds-VEGFR2 complexes

Post-simulation RMSD (Root Mean Square Deviation) analysis is crucial for evaluating the stability, conformational changes, and binding behavior of drug-protein interactions during molecular dynamics simulations. It helps determine whether the complex reaches equilibrium, with stable RMSD values indicating structural stability and significant fluctuations suggesting conformational shifts or ligand dissociation (Hassan et al., 2024). Analyzing the RMSD of the ligand within the binding pocket reveals how well the drug remains bound, while comparing RMSD across multiple drug candidates helps identify the most stable and effective compounds. Additionally, RMSD analysis can detect binding site rearrangements, ensure the reproducibility of simulations, and validate computational findings against experimental data, making it an essential tool in rational drug design and optimization (Pieroni et al., 2023). Therefore, to check the dynamic stability of shortlisted compounds-VEGFR2 complexes we calculated the RMSD over 400 ns simulation. As shown in the Figure 5 the post-simulation RMSD (Root Mean Square Deviation) analysis of VEGFR2 complexes with various drug candidates provides insights into the structural stability and dynamic behavior of each system over a 400 ns molecular dynamics (MD) simulation. The control VEGFR2 complex maintains a relatively stable RMSD between 2 and 3 Å over the entire 400 ns simulation (Figure 5a) which also validated by the superimposition of post-simulation retrieved structures (Figure 5b). The SA_0090-VEGFR2 complex exhibits minor variation among the analyzed complexes, with values fluctuating between 3 and 4 Å. This increased RMSD suggests that the complex undergoes substantial structural rearrangements throughout the simulation. The early stages of the simulation show a gradual increase in RMSD, indicating an initial conformational adjustment and accommodation for the ligand in the active site pocket (Figure 5c). However, the structural alignment of snapshots extracted from the simulation trajectories at 50 ns, 100 ns, 200 ns, 300 na and 400 ns demonstrated that the ligand remained stably positioned within the binding pocket throughout the simulation (Figure 5d). In contrast, the 17.3.1.7.8-VEGFR2 complex shows consistent RMSD values (∼3 Å) with minimal fluctuations, indicating stable binding and limited structural perturbation as compared to the control complex (Figure 5e). Additionally, the post-simulation retrieved structure further validated the stable binding of compounds in the active site throughout the simulation (Figure 5f). Similarly, the BMC_0005-VEGFR2 complex also demonstrates a relatively stable RMSD profile, maintaining values between 2 and 3 Å (Figure 5g). This stability suggests strong and consistent ligand binding, which preserves the protein’s structural integrity. The structural alignment in Figure 5h confirms these findings, showing limited conformational changes across the simulation timeframe. The ligand appears to remain anchored within the binding pocket. Overall, the RMSD data suggest that the 17.3.1.7.8 and BMC_0005 compounds form more stable VEGFR2 complexes, while SA_0090 leads to structural deviations.

Figure 5

Figure 5. RMSD trajectories analysis for the stability of shortlisted compounds-VEGFR2 complexes. (a,b) Represents the dynamic stability of control drug-VEGFR2 complex. (c,d) Represents the dynamic stability of SA_0090-VEGFR2 complex. (e,f) Represents the dynamic stability of 17.3.1.7.8-VEGFR2 complex. (g,h) Represents the dynamic stability of BMC_0005-VEGFR2 complex.

Post-simulation fluctuation analysis of shortlisted compounds-VEGFR2 complexes at residues level

RMSF (Root Mean Square Fluctuation) analysis is crucial in drug-protein interaction studies as it provides insights into the flexibility and dynamic behavior of protein residues during molecular dynamics simulations. It helps identify flexible and rigid regions, assess the stability of the drug-binding site, and reveal conformational changes induced by ligand binding. Reduced fluctuations in the binding site indicate stable drug binding, while increased fluctuations suggest weak or unstable interactions. This analysis is valuable for optimizing drug design by targeting dynamic hotspots and comparing the effects of different drug candidates on protein stability (Yang and Kar, 2024) (Abdullahi et al., 2024). Consequently, we calculated the RMSF to evaluate the flexibility of compounds-VEGFR2 complexes at residues level. As shown in the Figure 6 the RMSF analysis of the VEGFR2-ligand complexes highlights the flexibility of individual residues throughout the 400 ns molecular dynamics simulation. The RMSF values for the control, SA_0090, 17.3.1.7.8, and BMC_0005 complexes generally remain low (<2 Å) across most residues, indicating overall structural stability. However, a significant peak is observed around residue ∼190–210, suggesting increased flexibility in this loop region across all systems. This region is associated with the ligand-binding pocket, where interactions with the drug candidates induce structural adjustments. Among the complexes, the SA_0090 exhibits slightly higher fluctuations in this region compared to others (Figure 6a). The superimposition of retrieved post-simulation 3D structures revealed that the fluctuation was primarily localized to the loop regions. This indicates that the intrinsic flexibility of these loops likely drives the observed variations, potentially influencing the dynamic behavior of drug-protein interaction (Figures 6b–e). RMSF results further validated the RMSD data showing the stable binding of shortlisted compounds in the active site of VEGFR2 protein.

Figure 6

Figure 6. Residual fluctuation analysis of shortlisted compounds-VEGFR2 complexes by processing the RMSF trajectories. (a) Showing the fluctuation of each residues in control and shortlisted compounds-VEGFR2 complexes. (b–e) Showing the superimposition of 3D structures retrieved at different time point of simulation.

Structural compactness analysis of shortlisted compounds-VEGFR2 complexes

The radius of gyration (Rg) calculation is crucial for understanding drug-protein interactions as it provides insights into the conformational variations of the protein-ligand complex. Rg measures the distribution of atomic mass around the center of mass, indicating the compactness or flexibility of the system. A significant change in Rg values during molecular dynamics simulations suggests structural rearrangements, which can reflect the binding strength and stability of the drug-protein complex (Suleman et al., 2024a) (Khan et al., 2022). The Radius of Gyration (Rg) plots depict the compactness and structural stability of the VEGFR2 protein in complex with shortlisted compounds over a 400 ns molecular dynamics simulation (Figure 7). In the control-VEGFR2 system, the Rg fluctuates between 19.8 Å and 20.7 Å, indicating a relatively stable and compact structure throughout the simulation (Figure 7a). The SA_0090-VEGFR2 complex exhibits a slightly higher Rg range (20.1 Å to 21.0 Å) with increased fluctuations, suggesting a more flexible and less compact structure (Figure 7b). The 17.3.1.7.8-VEGFR2 complex maintains a relatively stable Rg (19.8 Å to 20.7 Å) with minimal variation, implying that this ligand stabilizes the VEGFR2 structure effectively (Figure 7c). Similarly, the BMC_0005-VEGFR2 complex shows slight fluctuations initially but stabilizes around 20.1 Å after 200 ns, indicating that this ligand maintains the structural integrity of VEGFR2 over time (Figure 7d). Overall, the 17.3.1.7.8 and BMC_0005 provide greater stability and compactness, suggesting stronger and more consistent interactions with the VEGFR2 protein as compared to the control and SA_0090.

Figure 7

Figure 7. Compactness analysis of shortlisted compounds-VEGFR2 complexes by calculating Rg. (a) Illustrating the structural compactness of control-VEGFR2 complex. (b) Illustrating the structural compactness of SA_0090-VEGFR2 complex. (c) Illustrating the structural compactness of 17.3.1.7.8-VEGFR2 complex. (d) Illustrating the structural compactness of BMC_0005-VEGFR2 complex.

Post-simulation hydrogen bonds analysis of shortlisted compounds-VEGFR2 complexes

Post-simulation hydrogen bond analysis is crucial for understanding drug-protein interactions as it provides insights into the stability, specificity, and binding affinity of the complex under dynamic conditions. It helps assess the persistence of key hydrogen bonds, which indicates interaction stability and supports binding affinity estimation. This analysis also validates molecular docking predictions by revealing whether interactions remain stable during molecular dynamics (MD) simulations (Suleman et al., 2021). Consequently, to check the binding stability of the shortlisted compounds and VEGFR2 protein we calculated the average post-simulation hydrogen bonds. The Figure 8 presents the post-simulation analysis of hydrogen bonds (H-bonds) over 400 ns for VEGFR2 in complex with different drug candidates and a control. Figure 8a represents the control-VEGFR2 complex, showing a relatively stable H-bond count with fluctuation between 80 and 120 ns. In contrast the SA_0090-VEGFR2, 17.3.1.7.8-VEGFR2 and BMC_0005-VEGFR2 complexes showed relatively similar pattern of hydrogen bonds. The average hydrogen bonds for the control, SA_0090-VEGFR2, 17.3.1.7.8-VEGFR2 and BMC_0005-VEGFR2 complexes were recorded to be 139.14, 134.57, 139.64 and 139.63 respectively (Figures 8a–d). Across all four systems, the H-bond count remains relatively stable, suggesting the stable binding of shortlisted compounds with the active site of the VEGFR2 protein.

Figure 8

Figure 8. Calculation of average hydrogen bonds in shortlisted compounds-VEGFR2 complexes. (a) Illustrates the number of hydrogen bonds in control-VEGFR2 complex. (b) Illustrates the number of hydrogen bonds in SA_0090-VEGFR2 complex. (c) Illustrates the number of hydrogen bonds in 17.3.1.7.8-VEGFR2 complex. (d) Illustrates the number of hydrogen bonds in BMC_0005-VEGFR2 complex.

Solvent Accessible Surface Area analysis of shortlisted compounds-VEGFR2 complexes

Solvent Accessible Surface Area (SASA) analysis is a key tool in computational biology and structural bioinformatics. It calculates how much of a molecule’s surface is exposed to surrounding solvent, such as water. This method is widely used to study biological macromolecules like proteins, DNA, and RNA. By analyzing exposed surface areas, researchers can better understand molecular interactions, locate ligand binding sites, examine protein-protein interfaces, and predict how molecules may behave in different environments (Suleman et al., 2024b). The Figure 9 presents Solvent Accessible Surface Area (SASA) over a 400 ns molecular dynamics (MD) simulation for VEGFR2 in complex with the control and three lead compounds: SA_0090 (b), 17.3.1.7.8 (c), and BMC_0005 (d). The Control-VEGFR2 complex displays moderate fluctuations in SASA values, ranging from approximately 15,000 to 16,500 Å², suggesting relatively stable but dynamic behavior (Figure 9a). In comparison, the SA_0090-VEGFR2 complex shows slightly higher and more variable SASA values, peaking near 17,500 Å², indicating enhanced solvent exposure, possibly due to ligand-induced conformational flexibility (Figure 9b). The 17.3.1.7.8-VEGFR2 complex shows the most consistent SASA values, largely maintaining within the 15,000–16,500 Å² range, suggesting stable ligand binding with minimal conformational disruption (Figure 9c). Conversely, the BMC_0005-VEGFR2 complex shows decreasing SASA values over time, indicating a potential compaction of the protein structure upon ligand binding (Figure 9d). Overall, the SASA trends suggest that different ligands influence VEGFR2 structural dynamics distinctively, with BMC_0005 potentially stabilizing a more compact conformation.

Figure 9

Figure 9. Surface area analysis of shortlisted compounds-VEGFR2 complexes. (a) Represents the SASA analysis of control-VEGFR2 complex. (b) Represents the SASA analysis of SA_0090-VEGFR2 complex. (c) Represents the SASA analysis of 17.3.1.7.8-VEGFR2 complex. (d) Represents the SASA analysis of BMC_0005-VEGFR2 complex.

Binding free energy calculation

MM/GBSA binding free energy is crucial for understanding drug-protein interactions as it provides an accurate estimation of binding affinity, helping to identify and prioritize potent drug candidates (Tuccinardi, 2021). By decomposing the total binding energy into van der Waals, electrostatics, and solvation components, it reveals the key forces driving molecular recognition and stability. This method aids in predicting the stability of drug-protein complexes during molecular dynamics simulations and guides rational drug design by highlighting areas for chemical modification to improve affinity. Additionally, MM/GBSA helps evaluate the effects of protein mutations on drug binding, offering insights into drug resistance. It complements experimental techniques by providing atomic-level details, accelerating drug discovery and optimizing lead compounds (Yau et al., 2024). Therefore, to check the binding strength of shortlisted compounds, we calculated the total binding free energy by using the MM/GBSA approach. The MM/GBSA analysis presented in Table 4 evaluates the binding free energies of the control and three shortlisted compounds (SA_0090, 17.3.1.7.8, and BMC_0005) by breaking down the energy contributions. The van der Waals (ΔEvdw) interactions remain relatively consistent across all compounds, with energies of −67.1847 ± 0.37 kcal/mol, −58.4277 ± 0.33 kcal/mol, −56.447 ± 0.45 kcal/mol and −66.7041 ± 0.32 kcal/mol for control, SA_0090, 17.3.1.7.8, and BMC_0005 respectively. However, electrostatic energy (ΔEele) is significantly more favorable in the shortlisted compounds, with SA_0090 showing the highest contribution (−109.9232 kcal/mol) compared to the control (−40.513 kcal/mol). However, the polar solvation energy (EGB) is also substantially higher for the shortlisted compounds, particularly for SA_0090 (124.6813 kcal/mol), suggesting greater desolvation penalties. The non-polar solvation energy (ESURF) remains similar across all samples, with slight variation, where SA_0090 has the largest contribution (−8.6036 kcal/mol). Despite the increased electrostatic attraction in the shortlisted compounds, the overall binding free energy (ΔG total) reveals that BMC_0005 (−61.476 kcal/mol) shows the most stable binding, as compared to the control (−60.3861 kcal/mol), while the recorded total binding free energies for 17.3.1.7.8 and SA_0090 were −52.7797 kcal/mol and −52.2732 kcal/mol, respectively. These results suggest that BMC_0005 may have the strongest binding affinity among the shortlisted compounds, while the higher solvation reduce the overall binding efficiency for 17.3.1.7.8 and SA_0090. Binding free results further validated the RMSD, RMSF and Rg data.

Table 4

Table 4. List of binding free energies calculated by using the MM/GBSA.

Conclusion

Papillary thyroid carcinoma (PTC) remains a significant clinical challenge, particularly in aggressive and metastatic cases where current therapeutic options are limited. Given the critical role of VEGFR-2 in tumor angiogenesis and its overexpression in PTC, targeting this receptor presents a promising avenue for therapeutic intervention. In this study, we employed a structure-based virtual screening approach, molecular dynamics simulations, and binding free energy calculations to identify potent VEGFR-2 inhibitors from the African natural product database. Our findings revealed three lead phytocompounds SA_0090 (−9.04 kcal/mol), 17.3.1.7.8 (−8.96 kcal/mol), and BMC_0005 (−8.33 kcal/mol) that demonstrated superior binding affinities compared to the control compound (-8.39 kcal/mol). Among them, 17.3.1.7.8 and BMC_0005 exhibited remarkable stability in molecular dynamics simulations, with minimal structural fluctuations and strong binding interactions. Additionally, ADMET analysis of 17.3.1.7.8 and BMC_0005 confirmed their favorable pharmacokinetic properties, including optimal solubility, gastrointestinal absorption, and non-hepatotoxicity, making them promising drug candidates. Notably, binding free energy calculations identified BMC_0005 as the most potent inhibitor, surpassing both the control and other shortlisted compounds in overall stability and affinity. These results underscore the potential of natural product-derived VEGFR-2 inhibitors as promising therapeutic agents for PTC. The limitations of this study is the lack of experimental validation and analysis of potential off target affects. Therefore, further in vitro and in vivo studies are essential to validate their biological efficacy and safety.

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

SW: Conceptualization, Formal Analysis, Methodology, Software, Visualization, Writing – original draft. LZ: Formal Analysis, Investigation, Methodology, Writing – review and editing. WZ: Data curation, Formal Analysis, Methodology, Validation, Writing – review and editing. XZ: Data curation, Formal Analysis, Investigation, Writing – review and editing. JM: Conceptualization, Data curation, Formal Analysis, Writing – review and editing. WX: Investigation, Methodology, Software, Supervision, Writing – review and editing. LY: Conceptualization, Project administration, Software, Supervision, Validation, Writing – review and editing.

Funding

The author(s) declare that no financial support was received for the research and/or publication of this article.

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.

Generative AI statement

The authors declare that no Generative AI was used in the creation of this manuscript.

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

Abdullahi, S. H., Moin, A. T., Uzairu, A., Umar, A. B., Ibrahim, M. T., Usman, M. T., et al. (2024). Molecular docking studies of some benzoxazole and benzothiazole derivatives as VEGFR-2 target inhibitors: in silico design, MD simulation, pharmacokinetics and DFT studies. Intell. Pharm. 2 (2), 232–250. doi:10.1016/j.ipha.2023.11.010

CrossRef Full Text | Google Scholar

Alanazi, M. M., Mahdy, H. A., Alsaif, N. A., Obaidullah, A. J., Alkahtani, H. M., Al-Mehizia, A. A., et al. (2021). New bis [(1, 2, 4) triazolo] (4, 3-a: 3′, 4′-c) quinoxaline derivatives as VEGFR-2 inhibitors and apoptosis inducers: design, synthesis, in silico studies, and anticancer evaluation. Bioorg. Chem. 112, 104949. doi:10.1016/j.bioorg.2021.104949

PubMed Abstract | CrossRef Full Text | Google Scholar

Aware, C. B., Patil, D. N., Suryawanshi, S. S., Mali, P. R., Rane, M. R., Gurav, R. G., et al. (2022). Natural bioactive products as promising therapeutics: a review of natural product-based drug development. South Afr. J. Bot. 151, 512–528. doi:10.1016/j.sajb.2022.05.028

CrossRef Full Text | Google Scholar

Case, D. A., Cheatham, T. E., Darden, T., Gohlke, H., Luo, R., Merz, K. M., et al. (2005). The Amber biomolecular simulation programs. J. Comput. Chem. 26 (16), 1668–1688. doi:10.1002/jcc.20290

PubMed Abstract | CrossRef Full Text | Google Scholar

Claesson-Welsh, L., and Welsh, M. (2013). VEGFA and tumour angiogenesis. J. Intern. Med. 273 (2), 114–127. doi:10.1111/joim.12019

PubMed Abstract | CrossRef Full Text | Google Scholar

Daina, A., Michielin, O., and Zoete, V. (2017). SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Sci. Rep. 7 (1), 42717. doi:10.1038/srep42717

PubMed Abstract | CrossRef Full Text | Google Scholar

Dhar, S., Sarkar, T., and Sa, G. (2022). Neoangiogenesis and immune-regulation: two armour of VEGF in the tumor microenvironment. J. Breast Cancer Res. 2 (1), 28–39. doi:10.46439/breastcancer.2.015

CrossRef Full Text | Google Scholar

El-Adl, K., Sakr, H. M., Yousef, R. G., Mehany, A. B. M., Metwaly, A. M., Elhendawy, M. A., et al. (2021). Discovery of new quinoxaline-2 (1H)-one-based anticancer agents targeting VEGFR-2 as inhibitors: design, synthesis, and anti-proliferative evaluation. Bioorg. Chem. 114, 105105. doi:10.1016/j.bioorg.2021.105105

PubMed Abstract | CrossRef Full Text | Google Scholar

Elkaeed, E. B., Yousef, R. G., Khalifa, M. M., Ibrahim, A., Mehany, A. B. M., Gobaara, I. M. M., et al. (2022). Discovery of new VEGFR-2 inhibitors: design, synthesis, anti-proliferative evaluation, docking, and MD simulation studies. Molecules 27 (19), 6203. doi:10.3390/molecules27196203

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghalehbandi, S., Yuzugulen, J., Pranjol, M. Z. I., and Pourgholami, M. H. (2023). The role of VEGF in cancer-induced angiogenesis and research progress of drugs targeting VEGF. Eur. J. Pharmacol. 949, 175586. doi:10.1016/j.ejphar.2023.175586

PubMed Abstract | CrossRef Full Text | Google Scholar

Giordano, D., Biancaniello, C., Argenio, M. A., and Facchiano, A. (2022). Drug design by pharmacophore and virtual screening approach. Pharmaceuticals 15 (5), 646. doi:10.3390/ph15050646

PubMed Abstract | CrossRef Full Text | Google Scholar

Gogiashvili, L. (2024). Expression of Vascular endothelial growth factor (VEGF) and CD34 in different thyroid disorders. Transl. Clin. Medicine-Georgian Med. J. 9 (1).

Google Scholar

Gordon, A. J., Dublin, J. C., Patel, E., Papazian, M., Chow, M. S., Persky, M. J., et al. (2022). American Thyroid Association guidelines and national trends in management of papillary thyroid carcinoma. JAMA Otolaryngology–Head Neck Surg. 148 (12), 1156–1163. doi:10.1001/jamaoto.2022.3360

PubMed Abstract | CrossRef Full Text | Google Scholar

Harahap, A. S., and Jung, C. K. (2024). Cytologic hallmarks and differential diagnosis of papillary thyroid carcinoma subtypes. J. Pathology Transl. Med. 58 (6), 265–282. doi:10.4132/jptm.2024.10.11

PubMed Abstract | CrossRef Full Text | Google Scholar

Hassan, A. M., Gattan, H. S., Faizo, A. A., Alruhaili, M. H., Alharbi, A. S., Bajrai, L. H., et al. (2024). Evaluating the binding potential and stability of drug-like compounds with the monkeypox virus VP39 protein using molecular dynamics simulations and free energy analysis. Pharmaceuticals 17 (12), 1617. doi:10.3390/ph17121617

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, A., Ali, S. S., Zahid, M. A., Abdelsalam, S. S., Albekairi, N., Al-Zoubi, R. M., et al. (2025). Exploring the dynamic interplay of deleterious variants on the RAF1–RAP1A binding in cancer: conformational analysis, binding free energy, and essential dynamics. Proteins Struct. Funct. Bioinforma. 93 (3), 684–701. doi:10.1002/prot.26759

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, A., Randhawa, A. W., Balouch, A. R., Mukhtar, N., Sayaf, A. M., Suleman, M., et al. (2022). Blocking key mutated hotspot residues in the RBD of the omicron variant (B. 1.1. 529) with medicinal compounds to disrupt the RBD-hACE2 complex using molecular screening and simulation approaches. RSC Adv. 12 (12), 7318–7327. doi:10.1039/d2ra00277a

PubMed Abstract | CrossRef Full Text | Google Scholar

Klebe, G. (2025). “Protein–ligand interactions as the basis for drug action,” in Drug design: from structure and mode-of-action to rational design concepts (Springer), 39–65.

Google Scholar

Kus, M., Ibragimow, I., and Piotrowska-Kempisty, H. (2023). Caco-2 cell line standardization with pharmaceutical requirements and in vitro model suitability for permeability assays. Pharmaceutics 15 (11), 2523. doi:10.3390/pharmaceutics15112523

PubMed Abstract | CrossRef Full Text | Google Scholar

Limaiem, F., Rehman, A., and Mazzoni, T. (2024). “Papillary thyroid carcinoma,” in StatPearls (Florida, United States: StatPearls Publishing).

Google Scholar

Lin, X., Li, X., and Lin, X. (2020). A review on applications of computational methods in drug screening and design. Molecules 25 (6), 1375. doi:10.3390/molecules25061375

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, Y., Ye, F., Jiang, Q., Liu, S., and Gong, Y. (2025). A visualization analysis of global research trends in targeted therapies for thyroid carcinoma (2013–2023). Medicine 104 (11), e41835. doi:10.1097/MD.0000000000041835

PubMed Abstract | CrossRef Full Text | Google Scholar

Mehrbod, P., Safari, H., Mollai, Z., Fotouhi, F., Mirfakhraei, Y., Entezari, H., et al. (2021). Potential antiviral effects of some native Iranian medicinal plants extracts and fractions against influenza A virus. BMC complementary Med. Ther. 21, 246–12. doi:10.1186/s12906-021-03423-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Modi, S. J., and Kulkarni, V. M. (2020). Discovery of VEGFR-2 inhibitors exerting significant anticancer activity against CD44+ and CD133+ cancer stem cells (CSCs): reversal of TGF-β induced epithelial-mesenchymal transition (EMT) in hepatocellular carcinoma. Eur. J. Med. Chem. 207, 112851. doi:10.1016/j.ejmech.2020.112851

PubMed Abstract | CrossRef Full Text | Google Scholar

Nhlapho, S., Nyathi, M., Ngwenya, B., Dube, T., Telukdarie, A., Munien, I., et al. (2024). Druggability of pharmaceutical compounds using Lipinski rules with machine learning. Sci. Pharm. 3 (4), 177–192. doi:10.58920/sciphar0304264

CrossRef Full Text | Google Scholar

Ntie-Kang, F., Zofou, D., Babiaka, S. B., Meudom, R., Scharfe, M., Lifongo, L. L., et al. (2013). AfroDb: a select highly potent and diverse natural product library from African medicinal plants. PloS one 8 (10), e78085. doi:10.1371/journal.pone.0078085

PubMed Abstract | CrossRef Full Text | Google Scholar

Oliveira, T. A. d., Silva, M., Maia, E., Silva, A., and Taranto, A. (2023). Virtual screening algorithms in drug discovery: a review focused on machine and deep learning methods. Drugs Drug Candidates 2 (2), 311–334. doi:10.3390/ddc2020017

CrossRef Full Text | Google Scholar

Pieroni, M., Madeddu, F., Di Martino, J., Arcieri, M., Parisi, V., Bottoni, P., et al. (2023). MD–ligand–receptor: a high-performance computing tool for characterizing ligand–receptor binding interactions in molecular dynamics trajectories. Int. J. Mol. Sci. 24 (14), 11671. doi:10.3390/ijms241411671

PubMed Abstract | CrossRef Full Text | Google Scholar

Pires, D., Blundell, T., and Ascher pkCSM, D. (2015). Predicting small-molecule pharmacokinetic and toxicity properties using graph-based signatures 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

Pollastri, M. P. (2010). Overview on the rule of five. Curr. Protoc. Pharmacol. 49 (1), Unit 9.12–9.12.8. doi:10.1002/0471141755.ph0912s49

PubMed Abstract | CrossRef Full Text | Google Scholar

Ravindranath, P. A., Forli, S., Goodsell, D. S., Olson, A. J., and Sanner, M. F. (2015). AutoDockFR: advances in protein-ligand docking with explicitly specified binding site flexibility. PLoS Comput. Biol. 11 (12), e1004586. doi:10.1371/journal.pcbi.1004586

PubMed Abstract | CrossRef Full Text | Google Scholar

Roe, D. R., and Cheatham, T. E. (2013). PTRAJ and CPPTRAJ: software for processing and analysis of molecular dynamics trajectory data. J. Chem. theory Comput. 9 (7), 3084–3095. doi:10.1021/ct400341p

PubMed Abstract | CrossRef Full Text | Google Scholar

Salomon-Ferrer, R., Case, D. A., and Walker, R. C. (2013). An overview of the Amber biomolecular simulation package. Wiley Interdiscip. Rev. Comput. Mol. Sci. 3 (2), 198–210. doi:10.1002/wcms.1121

CrossRef Full Text | Google Scholar

Salomon-Ferrer, R., Götz, A. W., Poole, D., Le Grand, S., and Walker, R. C. (2013). Routine microsecond molecular dynamics simulations with AMBER on GPUs. 2. Explicit solvent particle mesh Ewald. J. Chem. theory Comput. 9 (9), 3878–3888. doi:10.1021/ct400314y

PubMed Abstract | CrossRef Full Text | Google Scholar

Sana, S., Reddy, V. G., Bhandari, S., Reddy, T. S., Tokala, R., Sakla, A. P., et al. (2020). Exploration of carbamide derived pyrimidine-thioindole conjugates as potential VEGFR-2 inhibitors with anti-angiogenesis effect. Eur. J. Med. Chem. 200, 112457. doi:10.1016/j.ejmech.2020.112457

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarabipour, S., Ballmer-Hofer, K., and Hristova, K. (2016). VEGFR-2 conformational switch in response to ligand binding. Elife 5, e13876. doi:10.7554/eLife.13876

PubMed Abstract | CrossRef Full Text | Google Scholar

Sayaf, A. M., Kousar, K., Suleman, M., Albekairi, N. A., Alshammari, A., Mohammad, A., et al. (2024). Molecular exploration of natural and synthetic compounds databases for promising hypoxia inducible factor (HIF) Prolyl-4-hydroxylase domain (PHD) inhibitors using molecular simulation and free energy calculations. BMC Chem. 18 (1), 236. doi:10.1186/s13065-024-01347-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Sayaf, A. M., Ullah Khalid, S., Hameed, J. A., Alshammari, A., Khan, A., Mohammad, A., et al. (2023). Exploring the natural products chemical space through a molecular search to discover potential inhibitors that target the hypoxia-inducible factor (HIF) prolyl hydroxylase domain (PHD). Front. Pharmacol. 14, 1202128. doi:10.3389/fphar.2023.1202128

PubMed Abstract | CrossRef Full Text | Google Scholar

Shah, F. H., Nam, Y. S., Bang, J. Y., Hwang, I. S., Kim, D. H., Ki, M., et al. (2025). Targeting vascular endothelial growth receptor-2 (VEGFR-2): structural biology, functional insights, and therapeutic resistance. Archives Pharmacal Res. 48, 404–425. doi:10.1007/s12272-025-01545-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Shahrajabian, M. H., Sun, W., and Cheng, Q. (2022). The importance of flavonoids and phytochemicals of medicinal plants with antiviral activities. Mini-Reviews Org. Chem. 19 (3), 293–318. doi:10.2174/1570178618666210707161025

CrossRef Full Text | Google Scholar

Somoza, A. D., Bui, H., Samaan, S., Dhanda-Patil, R., and Mutasim, D. F. (2013). Cutaneous metastasis as the presenting sign of papillary thyroid carcinoma. J. Cutan. Pathology 40 (2), 274–278. doi:10.1111/cup.12044

PubMed Abstract | CrossRef Full Text | Google Scholar

Suleimen, Y. M., Metwaly, A. M., Mostafa, A. E., Elkaeed, E. B., Liu, H. W., Basnet, B. B., et al. (2021). Isolation, crystal structure, and in silico aromatase inhibition activity of ergosta-5, 22-dien-3β-ol from the fungus gyromitra esculenta. J. Chem. 2021 (1), 1–10. doi:10.1155/2021/5529786

CrossRef Full Text | Google Scholar

Suleman, M., Moltrasio, C., Tricarico, P. M., Marzano, A. V., and Crovella, S. (2024a). Natural compounds targeting thymic stromal lymphopoietin (TSLP): a promising therapeutic strategy for atopic dermatitis. Biomolecules 14 (12), 1521. doi:10.3390/biom14121521

PubMed Abstract | CrossRef Full Text | Google Scholar

Suleman, M., Sayaf, A. M., Aftab, S., Alissa, M., Alghamdi, A., Alghamdi, S. A., et al. (2025). Medicinal phytocompounds as potential inhibitors of p300-hif1α interaction: a structure-based screening and molecular dynamics simulation study. Pharmaceuticals 18 (4), 602. doi:10.3390/ph18040602

PubMed Abstract | CrossRef Full Text | Google Scholar

Suleman, M., Sayaf, A. M., Khan, A., Khan, S. A., Albekairi, N. A., Alshammari, A., et al. (2024b). Molecular screening of phytocompounds targeting the interface between influenza A NS1 and TRIM25 to enhance host immune responses. J. Infect. Public Health 17 (7), 102448. doi:10.1016/j.jiph.2024.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Suleman, M., Yousafi, Q., Ali, J., Ali, S. S., Hussain, Z., Ali, S., et al. (2021). Bioinformatics analysis of the differences in the binding profile of the wild-type and mutants of the SARS-CoV-2 spike protein variants with the ACE2 receptor. Comput. Biol. Med. 138, 104936. doi:10.1016/j.compbiomed.2021.104936

PubMed Abstract | CrossRef Full Text | Google Scholar

Tille, J.-C., Wood, J., Mandriota, S. J., Schnell, C., Ferrari, S., Mestan, J., et al. (2001). Vascular endothelial growth factor (VEGF) receptor-2 antagonists inhibit VEGF-and basic fibroblast growth factor-induced angiogenesis in vivo and in vitro. J. Pharmacol. Exp. Ther. 299 (3), 1073–1085. doi:10.1016/s0022-3565(24)29231-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Tuccinardi, T. (2021). What is the current value of MM/PBSA and MM/GBSA methods in drug discovery? Expert Opin. Drug Discov. 16 (11), 1233–1237. doi:10.1080/17460441.2021.1942836

PubMed Abstract | CrossRef Full Text | Google Scholar

Tuttle, R. M., Haugen, B., and Perrier, N. D. (2017). Updated American Joint Committee on cancer/tumor-node-metastasis staging system for differentiated and anaplastic thyroid cancer: what changed and why? NY, USA: Mary Ann Liebert, 751–756.

Google Scholar

Yang, S., and Kar, S. (2024). Protracted molecular dynamics and secondary structure introspection to identify dual-target inhibitors of Nipah virus exerting approved small molecules repurposing. Sci. Rep. 14 (1), 3696. doi:10.1038/s41598-024-54281-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Yau, M. Q., Liew, C. W. Y., Toh, J. H., and Loo, J. S. E. (2024). A head-to-head comparison of MM/PBSA and MM/GBSA in predicting binding affinities for the CB1 cannabinoid ligands. J. Mol. Model. 30 (11), 390. doi:10.1007/s00894-024-06189-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., and Xu, S. (2024). High aggressiveness of papillary thyroid cancer: from clinical evidence to regulatory cellular networks. Cell Death Discov. 10 (1), 378. doi:10.1038/s41420-024-02157-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, S., Krumberger, M., Morris, M. A., Parrocha, C. M. T., Kreutzer, A. G., and Nowick, J. S. (2021). Structure-based drug design of an inhibitor of the SARS-CoV-2 (COVID-19) main protease using free software: a tutorial for students and scientists. Eur. J. Med. Chem. 218, 113390. doi:10.1016/j.ejmech.2021.113390

PubMed Abstract | CrossRef Full Text | Google Scholar