Screening of plant-related small molecules is considered an initial step in understanding the multi-target effect of medicinal plants against diseases. In the current study, database and literature search were carried out for obtaining the main active constituents of plants. In the literature and databases, the search was confined to “O. heteracanthum,” “A. ehrenbergiana,” “O. vaillantii,” “C. rotundus,” “C. carandas,” “C. spinarum,” and “C. sinensis.” PubChem (https://pubchem.ncbi.nlm.nih.gov/) and Google Scholar were considered for the extensive literature survey, while the phytochemical databases, including KNApSAcK (http://www.knapsackfamily.com/KNApSAcK/) (Shinbo et al., 2006) and IMPPAT (https://cb.imsc.res.in/imppat/home) (Mohanraj et al., 2018), were searched for obtaining the active compounds of selected plants. After retrieval, the drug-like potential of active compounds was analyzed by predicting their oral bioavailability (OB) and drug-likeness (DL) values. In pharmacology, OB denotes the fraction of small molecules that successfully enter the systemic circulation after an oral dose, thereby enabling them to exert their desired pharmacological effects (Pathak and Raghuvanshi, 2015). A minimum OB of 30% is generally employed as a benchmark criterion for screening potential drug candidates since it indicates that a considerable portion of the oral usage of small molecules is capable of being absorbed and reaching systemic circulation. Conversely, small molecules with OB < 30% are likely to exhibit limited efficacy due to poor absorption. Thus, a higher OB value is positively correlated with drug effectiveness and is a key determinant of drug development and optimization. In a similar vein, DL analysis is employed to assess the potential of a molecule to serve as an oral drug by evaluating its bioavailability using qualitative measures (Jia et al., 2020). Following that, the DL and OB values of compounds were obtained using Molsoft (https://molsoft.com/mprop/) (James et al., 2015) and SwissADME software (https://www.swissadme.ch/) (Daina et al., 2017), and only those compounds having OB and DL values greater than 30% and 18% were considered the finalist compounds. Furthermore, the molecular weight and two-dimensional (2D) structure of predicted phytoconstituents were collected from PubChem (Kim et al., 2019) and Molinspiration (https://molinspiration.com/) (Molinspiration, 2011).

Although DL and OB are key elements in evaluating the drug-like potential of a compound, it is important to note that these parameters do not serve as the sole determinants of a compound’s effectiveness for further development. ADMET properties, which encompass absorption, distribution, metabolism, excretion, and toxicity, are equally important factors that must be evaluated to assess a compound’s safety and efficacy (Rehman et al., 2022). The evaluation of ADMET properties is conducted during the initial phases of drug development, which substantially plays an essential role in predicting the potential safety, effectiveness, and pharmacokinetics of a drug. To this end, the SwissADME server (https://www.swissadme.ch/) (Daina et al., 2017) and the ProTox II tool (http://tox.charite.de/protox_II) (Banerjee et al., 2018) were used to assess the ADMET properties of small biologically active molecules. Small molecules demonstrating excellent absorption, favorable solubility properties, and low toxicity were selected for further analysis.

To determine the integrative efficacy of the compounds derived from local plants in Saudi Arabia, we used two different platforms, namely, STITCH (http://stitch.embl.de/) (Kuhn et al., 2007) and SwissTargetPrediction (Gfeller et al., 2014) databases, with a species limitation of “Homo sapiens.” Furthermore, the SMILES representation of the chosen ingredients was used as input for the SwissTargetPrediction (https://www.swisstargetprediction.ch/) platform. Specifically, only proteins with a probability > 0.7 were considered putative targets of Saudi plants. For the STITCH database, targets exhibiting a combined score >0.7 were considered the significant targets of particular compounds.

Microarray data analysis is a technique used for the identification of genes that are differentially expressed among different conditions (Noor et al., 2021; Sufyan et al., 2021). This involves comparing the gene expression levels of samples from disease and control groups and calculating statistical measures to determine which genes are significantly differentially expressed. In this study, gene expression data were acquired from the NCBI–GEO database (https://www.ncbi.nlm.nih.gov/) (Barrett et al., 2005). The search term “Ovarian Cancer” was used to retrieve four microarray datasets, namely, GSE54388, GSE69428, GSE36668, and GSE40595. The selection of datasets was based on their inclusion of both affected and normal ovaries; the dataset must have >3 samples. After data selection, the gene expression data were processed through the limma package of R (Ritchie et al., 2015) for the screening of differentially expressed genes (DEGs) as disease-specific targets for ovarian cancer.

Upon identifying putative targets associated with both ovarian cancer and local plants, a Venn diagram was created to determine overlapping proteins. Subsequently, these identified genes were regarded as potential targets of indigenous Saudi Arabian plants, holding promise as potential biomarkers for intercepting the pathophysiological processes associated with ovarian cancer. Additionally, Cytoscape version 3.8 (https://cytoscape.org/) (Smoot et al., 2011) was used to construct an active ingredient–target network based on these overlapping genes. Within the compound–target network, the nodes indicated the small molecules along with their corresponding targets, while gray dotted lines depicted the possible interactions among these nodes. Finally, the NetworkAnalyzer plugin was used to evaluate the connectivity of targets exhibited by the chemical constituents among networks.

After identifying the common target, our study categorizes these targets based on their Gene Ontology (GO). GO covers three aspects of biology: cellular components (CCs), molecular functions (MFs), and biological processes (BPs). The BP term describes the protein’s role in various biological activities, and the CC term signifies the precise intracellular localization of a protein. Concurrently, the MF records the specific molecular tasks and interactions in which the protein or gene is involved. To achieve this, the current study used the DAVID database (https://david.ncifcrf.gov/tools.jsp) (Dennis et al., 2003) for GO and Kyoto Encyclopedia for Genes and Genomes (KEGG) pathway analyses. However, to filter-out those with significant relevance, a statistical cutoff of p-value<0.05 was implemented.

The common targets were then analyzed using the STRING database (https://string-db.org/) (von Mering et al., 2003) to produce a protein–protein interaction (PPI) network. PPIs exhibit high specificity, adaptability, and versatility and are, therefore, of considerable significance. PPI network common genes were then analyzed using Cytoscape version 3.8 (https://cytoscape.org/) (Smoot et al., 2011) to identify hub genes. Hub genes exhibit a high number of interactions with other proteins and are considered essential elements of the network as they contribute to its integrity and stability. Furthermore, central or ‘hub’ genes hold significance in a range of disease-linked pathways and biological processes, suggesting their pivotal role in diverse cellular functions. In the current study, we examined the topology of the PPI network to identify ‘hub’ genes. Metrics such as degree, closeness, and betweenness centrality can facilitate this process. Specifically, we employed the degree methods available in cytoHubba to discern these hub genes.

The effects of Saudi Arabian flora on ovarian cancer were investigated by constructing compound–target and target protein–pathway networks using Cytoscape version 3.8 (https://cytoscape.org/) (Smoot et al., 2011). By merging these networks, a final compound–target–pathway network was subsequently established. The networks comprise various nodes and edges, where nodes represent pathways, compounds, and targets relevant to the disease, while edges indicate the interactions among these nodes. The compound–target–pathway network offers valuable insights regarding the synergistic activity of phytoconstituents in ovarian cancer treatment.

Using advanced molecular docking methodologies, intricate interactions between small molecules and anticipated protein structures were analyzed. This rigorous analytical procedure illuminated potential pharmacological combinations, suggesting potential synergistic responses for precision therapeutics. Initially, the protein structures were obtained from the RCSB PDB database, which were then carefully prepared for docking. Missing atoms, residues, and side chains were addressed using the protein preparation tools within Chimera, ensuring a complete and accurate representation of the protein’s conformation. Water molecules and non-essential ions were removed, and hydrogen atoms were added to the protein structures to reflect physiological pH conditions. After this refinement, the proteins were subjected to an energy minimization process using the steepest descent and conjugate gradient methods, terminating at an energy convergence threshold of 0.01 kcal/mol to attain stable conformations. Subsequent docking studies were conducted using AutoDock Vina 1.1.2 implemented within the PyRx 0.8 interface (https://sourceforge.net/projects/pyrx/) (Dallakyan and Olson, 2015). The SMILES strings for the phytoconstituents were sourced from the PubChem database and converted to 3D structures using Open Babel, which is integrated into the PyRx platform. A comprehensive energy minimization was configured for a total of 2000 steps to ensure the reliability of the ligand conformations prior to the docking simulations. Later, the active site residues of the target proteins were uncovered with the assistance of the CASTp tool (http://sts.bioe.uic.edu/castp/index.html?201l) (Dundas et al., 2006). PyRx 0.8 was employed for target docking, which facilitated the calculation of binding affinities between the small molecules and target proteins. In our study, docking scores were used as indicators of binding affinity between the identified phytochemicals and target proteins. A docking score of less than −5.00 kcal/mol was considered to indicate strong binding affinity, while a score of less than −7.00 kcal/mol was interpreted as signifying very strong binding affinity. These thresholds were based on established computational chemistry conventions and were used to prioritize compounds for further analysis. The docked complexes were then visualized employing Discovery Studio (https://discover.3ds.com/discovery-studio-visualizer-download) (Biovia, 2017), PyMOL (https://pymol.org/2/) (Seeliger and de Groot, 2010), and ChimeraX (https://www.cgl.ucsf.edu/chimerax/) (Pettersen et al., 2021) programs for a better understanding of the binding interactions.

All-atom molecular dynamics (MD) simulation is a computational approach that precisely represents each atom and bond within a system, facilitating an in-depth analysis of molecular dynamics (Sahoo et al., 2020). This method calculates the movements of every atom in a system, which are determined by the interactions amongst them, represented by interatomic potentials (Sahoo et al., 2022). In this research, to execute the MD simulations of the ultimate complexes, we implemented the OPLS-AA/L force field using GROMACS 2018 software. The three-dimensional configurations of the protein served as the initial structural inputs for the simulations. To enhance their quality and suitability for the simulations, additional optimization was performed using Dock Prep (Pettersen et al., 2004). The parameterization of the active ingredients was executed using the SwissParam webserver (Zoete et al., 2011). Subsequently, MD simulations were conducted for a duration of 20 nanoseconds (ns), following a methodology similar to previous studies (Alamri, 2020). Various general MD simulation parameters were assessed for each complex, including root mean square deviation (RMSD), radius of gyration (Rog), and root mean square fluctuation (RMSF) (Needle et al., 2015). These parameters provide important information about the stability, conformational changes, and flexibility of the complexes during the simulation period.

The screening of active compounds represents a critical initial step in the identification of promising phytoconstituents with potential therapeutic efficacy against ovarian cancer. After searching and filtering, four compounds, namely, hispidulin, apigenin, luteolin, and beta-sitosterol were obtained from O. heteracanthum; two compounds, namely, stigmasterol and ascorbic acid were obtained from A. ehrenbergiana; five compounds namely, octopamine, 3,4-dihydroxybenzoic acid, cyperene, rotundine A, and rotundine B were obtained from O. vaillantii; six compounds namely, orientin, quercetin, kaempferol, pungenin, rotundine B, and cyperene were obtained from C. rotundus; nine compounds namely, carissic acid, alpha1-sitosterol, citric acid, d-tartaric acid, ascorbic acid, beta-ionone, anethole, ursolic acid, and beta-sitosterol were obtained from C. carandas; six compounds, namely, carissic acid, kaempferol, quercetin, digitoxigenin, ursolic acid, and beta-sitosterol were obtained from C. spinarum; and four compounds, namely, (−)-epicatechin gallate, (+)-catechin, quercetin, and kaempferol were obtained from C. sinensis. These 36 compounds were considered potential phytoconstituents of local Saudi plants as they fulfill the criteria of molecular weight <500 g/mol, DL ≥ 0.18, and OB ≥ 0.30 (Table 1).

The chosen compounds were subjected to ADME analysis, and 12 compounds, namely, hispidulin, stigmasterol, ascorbic acid, octopamine, cyperene, kaempferol, pungenin, citric acid, d-tartaric acid, beta-sitosterol, (−)-epicatechin gallate, and (+)-catechin, were identified as the main active constituents, with local good GI absorption and limited permeability of BBB, as presented in Table 2. Within the scope of ADMET analysis, hepatotoxicity refers to the capacity of compounds to inflict harm on the liver, potentially resulting in liver malfunction or complete failure. Notably, all selected compounds were found to have inactive hepatotoxicity in the current study. On the other hand, carcinogenicity characterizes the potential to induce cancer, while mutagenicity denotes its capacity to provoke DNA alterations that could lead to developmental anomalies, mutations in cancer, or many others. Compared to our results, the final active ingredients demonstrated themselves to be non-mutagenic and non-carcinogenic. In sum, these outcomes further support that the native plants of Saudi Arabia have drug-like potential that could be integral to preventing and treating diseases.

Analysis of microarray data offers crucial insights into the molecular processes associated with disease progression and assists in the identification of possible targets for therapeutic intervention. Four gene expression datasets, namely, GSE54388, GSE69428, GSE36668, and GSE40595, were obtained from the NCBI–GEO database. Subsequent processing for the detection of DEGs was carried out using the limma package (Table 3, Supplementary Material S4: Supplementary Tables S1–S3). Within the limma package, the LogFC >1.0 and p-value <0.05 criteria were set for the screening of DEGs. From the GSE54388 dataset, a total of 4,687 DEGs were obtained (2,457 upregulated and 2,230 downregulated), 2,084 DEGs (1,020 upregulated and 1,064 downregulated) from GSE69428, 6,924 DEGs (3,147 downregulated and 3,777 upregulated) from GSE36668, and 2,840 DEGs (1,423 upregulated and 1,417 downregulated) from the GSE40595 dataset (Figure 1, Supplementary Material S5). The identified DEGs were subsequently used as therapeutic targets for ovarian cancer.

Following the prediction of disease-associated targets, a Venn tool was employed to uncover shared genes between plant-associated targets and disease-related targets. SwissTargetPrediction and STITCH databases yielded 469 potential targets of active compounds (Supplementary Material S1: Supplementary Table S1). Finally, 200 targets were found to be common among both plants and diseases (Supplementary Material S2: Supplementary Table S1). Compounds targeting these 200 genes were then imported into Cytoscape. This network suggests that these 200 identified targets could potentially work in synergy when these indigenous Saudi plants are used as an anticancer remedy.

Analyses for GO and pathway enrichment of common 200 target genes were carried out employing the DAVID tool to discern their biological characteristics (Supplementary Material S3: Supplementary Tables S1–S4). A total of 388 BP, 57 CC, and 109 MF terms met the p-value criterion of <0.05. Analysis of the top GO terms suggested that the shared genes were primarily associated with the positive regulation of the MAPK, ERK1, and ERK2 cascade, positive regulation of vasoconstriction, intracellular signal transduction, regulation of cell proliferation, positive regulation of the apoptotic process, Hsp90 protein binding, estradiol 17-beta-dehydrogenase activity, estrogen response element binding, steroid binding, and steroid hormone receptor activity (Figure 2). A total of 103 KEGG pathways were identified, indicating that the intersecting targets were predominantly involved in the mechanisms of growth hormone action, secretion, and synthesis, pathways in cancer, the estrogen signaling pathway, pI3K-Akt and ErbB signaling pathways, progesterone-mediated oocyte maturation, endometrial cancer, and many others (Figure 3).

Constructing a PPI network is a fundamental step in network pharmacology, offering a structure to comprehend the intricate interactions between biological molecules and their impacts on disease pathways. Interactions among 200 common genes were predicted using the STRING database. This network comprised a total of 196 genes, forming 1,426 interactions with each other. One protein was found to interact with many proteins; therefore, the top 10 proteins with greater connectivity than other nodes were selected as hub genes. These genes are AKT1 (91), VEGFA (73), HSP90AA1 (63), JUN (63), HIF1A (60), PTGS2 (530, MAPK1 (53), PPARA (38), PIK3R1 (36), and AR (36) (Figure 4, Supplementary Material S6: Supplementary Table S1).

For gaining a deep understanding of the underlying action mechanisms of indigenous plants in ovarian cancer, an integrated “active ingredient–target–pathway” network was constructed using the GO and KEGG pathway databases, and non-ovarian cancer-associated pathways were excluded from the analysis (Figure 5). Subsequently, a comprehensive analysis of the PPI network and active ingredient–target–pathway network was performed at the systematic level to identify top-ranked proteins, including AKT1 and VEGFA, which were further subjected to molecular docking analysis. AKT1 and VEGFA were selected as most of the compounds targeted these genes, and these genes are mainly linked to ovarian cancer-related pathways, including the estrogen signaling pathway, progesterone-mediated oocyte maturation, the prolactin signaling pathway, Rap1 and ErbB signaling pathways, colorectal cancer, growth hormone action, secretion, and synthesis, and endometrial cancer.

After network analysis, two proteins, AKT1 and VEGFA, were chosen for docking analysis, as these proteins were found to be targeted by multiple compounds and also linked to disease-relevant pathways. All selected compounds were docked against AKT1 and VEGFA, allowing for the assessment of their binding affinity, interaction stability, and free energy within the target protein’s active site. These proteins were docked with selected compounds; however, the top five compounds based on their binding affinity were selected. In the case of the AKT1 protein, the top five compounds with maximum binding affinities are ascorbic acid (−11.1227), citric acid (−10.2178), kaempferol (−11.4719), stigmasterol (−10.0388), and d-tartaric acid (−10.9147) (Figure 6). In the case of the AKT1–ascorbic acid complex, binding affinity contributed to H-bonds with the Ser A:240; Leu A:239; Arg A:243; Arg C:2; and Glu A:341 residues. In terms of citric acid, AKT1 formed H-bonds with Glu A:341 and Gly A:345; for kaempferol, AKT1 demonstrated H-bond interaction with Glu A:341; for stigmasterol, the H-bond interaction was found with Ser A:240 and Arg C:2; and d-tartaric acid showed H-bond interaction with Arg C:2; Glu A:341; Arg A:346; Gly C:1; and Leu A:239 residues (Table 4).

In the case of the VEGFA protein, (+)-catechin (−15.6498), citric acid (−14.7792), (−)-epicatechin gallate (−15.1182), hispidulin (−17.3714), and octopamine (−17.517) were found to have more binding affinity than other active ingredients (Figure 7). In the case of the VEGFA-(+)–catechin complex, binding affinity contributed to H-bonds with Asn B:158; Ser A:179; Ser A:176; Lys A:143; and Ser A:177; for citric acid, VEGFA has H-bond interactions with Ser A:177; Ser A:179; Leu A:178; Gln B:160; Ser B:176 Ser A:177; Ser A:179; Leu A:178; Gln B:160; and Ser B:176; for (−)-epicatechin gallate, VEGFA demonstrated H-bond interactions with Asn A:158; Thr B:180; Gln B:160; Lys A:143; and Asp A:144; for hispidulin, H-bond interactions were found with Asn A:158; Gln B:160; Ser A:179; and Ser B:176; and octopamine showed H-bond interactions with Ser A:177; Ser A:179; Ser B:176; and Gln A:171 residues (Table 5).

All-atom Molecular Dynamics (MD) simulations were conducted to analyze the relationship between the AKT1 and VEGFA proteins and to assess the stability of the active compounds. The simulation was run at 20 ns using GROMACS software. For each protein, the top two docked complexes based on their binding energies were chosen for MD simulation for evaluating their RMSD and RMSF values, which assist in analyzing conformational changes and interactions among small molecules and target proteins. RMSD offers an understanding of the extent of deviation that a set of atoms (complexes, compounds, and proteins) experience from their original reference structure. On the other hand, RMSF calculates the temporal progression of the average discrepancy for each residue from its initial position in the minimized structures. The RMSD analysis revealed that selected complexes stabilized within a 20-ns timeframe. The average RMSD of the backbone atoms of AKT1 with kaempferol and ascorbic acid–protein systems was found to be ∼0.175 and ∼0.23Å, respectively (Figure 8A), while for VEGFA, the average RMSD for hispidulin and octopamine was found to be 0.45 and 0.5 Å, respectively (Figure 9A). On the other hand, the RMSD of ligand atoms for AKT1 with ascorbic acid exhibits an average RMSD value of ∼0.04, with a sharp increase in RMSD (∼0.12 Å) observed after 4–13 ns of simulation (Figure 8B), while for kaempferol, the RMSD falls within the range of ∼0.08 to ∼0.10 Å. In the case of the VEGFA protein, the average RMSD value for octopamine and hispidulin was ∼0.14 and ∼0.12 Å (Figure 9B). These findings suggest that octopamine was more stable at 20 ns in the case of the VEGFA protein, while for AKT1, kaempferol was found to be more stable. This was further substantiated by analyzing RMSF versus the residue number of AKT1, which showed that in comparison to AKT1–kaempferol, the AKT1–ascorbic acid complex exhibited greater fluctuations in backbone residues (Figure 8C). However, in terms of the VEGFA complex, hispidulin indicated greater fluctuations in backbone residues than octopamine (Figure 9C).

The radius of gyration (Rg) serves as a crucial metric to evaluate changes in the compactness of a docked complex. For AKT1, ascorbic acid exhibited more pronounced fluctuations than kaempferol (Figure 8D). In the context of VEGFA, notably, the Rg values for VEGFA–hispidulin exhibited a narrower distribution, denoting a more consistent structural compactness relative to the VEGFA–octopamine complex, which displayed a broader range of Rg values, suggestive of increased flexibility and structural variability. This differential behavior in Rg implies that the binding of hispidulin may confer a greater degree of conformational restraint upon the VEGFA protein, potentially influencing the stability and functional specificity of the ligand–protein interaction (Figure 9D). Furthermore, insights derived from the evaluation of Rg, RMSF, and RMSD metrics indicate that all complexes maintained consistent stability throughout the simulation period. Therefore, it could be concluded that these small molecules hold the potential to be used as an inhibitor for AKT1 and VEGFA: however, the finding needs further in vivo and in vitro validation.