The RNA-sequencing datasets related to PTB were collected from the publicly available Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) of the National Center for Biotechnology Information (NCBI) (Barrett et al., 2012). The GEO database stands as the premier openly accessible repository housing a vast array of high-throughput molecular abundance data, with a primary emphasis on gene expression studies. Within its extensive collection, it encompasses a diverse spectrum of data types, spanning from traditional microarray-based evaluations of mRNA, genomic DNA, and protein levels to non-array methodologies like serial analysis of gene expression (SAGE) and mass spectrometry for comprehensive proteomic analyses. In our investigation, we procured two distinct RNA sequence datasets specifically focusing on females afflicted with preterm birth (PTB), with one dataset including pertinent data on preeclampsia (PE) status, sourced from the GEO accession GSE148241 (Yang et al., 2020). and the other with PE + IUGR (preeclampsia and intrauterine fetal growth restriction) from the GEO accession GSE114691 (Awamleh et al., 2019). The first RNA sequence dataset, GSE148241, encompassed RNA-seq data obtained from a total of 41 placenta tissue samples. Which included the placentae of 9 individuals diagnosed with early-onset severe preeclampsia (EOSE), in comparison with those of 32 healthy controls, utilizing the Illumina HiSeq 2,500 platform. Moving to the second RNA sequence dataset, GSE114691, high-throughput expression profiling was conducted using the Illumina HiSeq 2000 platform. This dataset, scrutinized in our investigation, entailed RNA-seq data from 79 samples. Of our main interest were the samples from intrauterine growth restriction (IUGR) and early-onset preeclampsia (PE), the pregnancies complicated by both conditions. The placental specimens sourced from four distinct patient cohorts, including PE + IUGR. 20 samples with PE_IUGR were included in the analysis, while a control group comprising 21 subjects served as reference.

To determine the differential expression of the genes in the above datasets, differential expression analysis was performed using the DESeq2 library in the R programming language (Love et al., 2014). The library operates by augmenting the robustness and comprehensibility of differential expression estimations via the implementation of shrinkage estimation techniques for dispersions and fold changes. In this analytical framework, differential expression analysis was governed by stringent criteria, requiring a minimum log2 fold change threshold of >1 for upregulated genes and <-1 for downregulated genes. Additionally, a significance threshold of p-value >0.05 was uniformly applied across both the selected Pre-Eclampsia (PE) and PE + Intrauterine Growth Restriction (IUGR) datasets, ensuring rigorous statistical scrutiny and reproducibility in the identified gene expression alterations.

Functional enrichment analysis of the set of DEGs was performed using the gseGO and gseKEGG methods of the clusterprofiler in the R programming language. This approach allows for the extraction of semantic pathways and Gene Ontology (GO) terms in which the list of DEGs is found to be enriched. The clusterProfiler was used to construct the GO and KEGG (Kyoto Encyclopedia of Genes and Genomes) datasets for multiple model organisms. The primary purpose of this approach was to facilitate the comparison of functional profiles across distinct conditions at a single level, such as various treatment groups (Wu et al., 2021). The parameters were set as follows: minimum gene set size, 3; maximum set size, 800; and p-value < 0.05.

We particularly focused on the upregulated DEGs from both datasets side by side; these DEGs were further analyzed to reveal the PPI network using STRING (Von Mering et al., 2003) and the plug-in in Cytoscape 3.9.1 (Kohl et al., 2011; Saito et al., 2012). Cluster (subnetwork) analysis was performed using MCODE, which allows for the recognition of local subnetworks from the larger network that have the most significance based upon the higher score and strength of the nodes and edges connecting them. These gene clusters were analyzed to closely study the processes and pathways in which the genes were enriched using ClueGO (Bindea et al., 2009) in Cytoscape.

Hub genes were defined as genes with a high correlation in candidate modules. These genes were basically defined by the number of associations (i.e., node degree); highly associated genes were denoted as hub genes. The hub genes were identified using CytoHubba in Cytoscape (Chin et al., 2014). This approach employs twelve distinct methods, including MCC, EPC, MNC, DMNC, eccentricity, radiality, stress, betweenness, closeness, clustering coefficient, bottleneck, and degree, to identify key genes within biological networks. These methods assess various network properties such as connectivity, centrality, and structural importance to pinpoint genes that play pivotal roles in regulating network behavior. Genes that were identified by five or more of these methods were classified as hub genes, signifying their significance in the network topology and their potential importance in biological processes under investigation.

The 3D model of the identified Hub proteins was investigated with the aid of protein-BLAST (Basic Local Alignment Search Tool) for identifying the structural similarity of the desired proteins in the PDB repository (Stover and Cavalcanti, 2017). In the context of incomplete structural data available in the PDB, indicating a gap in the structure, we employed an alternative homology modeling technique to construct a homology model for the Hub protein common to both datasets. To evaluate the completeness or maximum coverage of the 3D model, we utilized the SWISS-MODEL web server for homology modeling. (Robin et al., 2024). Following the generation of the respective 3D homology models, structural validation was conducted using verify-3D, ERRAT, and PROCHECK-based assessments in Saves6.0.(Dym et al., 2012).

In addition, a 3D compound library comprising bioactive flavonoids reported from various sources and from various classes was retrieved from the PubChem compound database (https://pubchem.ncbi.nlm.nih.gov). The 3D structures of all the reported compounds were retrieved, and following a visual examination of each 3D compound, the structures were initially downloaded in SDF (structure data format) and subsequently translated to PDB (protein database) using Open Babel software (Yoshikawa and Hutchison, 2019). All the 3D compounds were subjected to Pfizer’s rule of five for evaluating drug-likeness criteria (Tinworth and Young, 2020). Following the evaluation, only 3D bioactive flavonoids that satisfied all five of Pfizer’s Rule of Five physicochemical requirements were chosen because any one infraction was treated as a failure or an elimination aspect. Only eligible compounds were included and further subjected to docking analysis. According to the criteria previously reported by Azmi et al., in 2023, we selected different PTB drugs for reference and comparison (Azmi et al., 2023).

PyRx is a platform that automates and integrates various tools for molecular docking simulations, streamlining the process of virtual screening for small-molecule library screening. In our study, we utilized PyRx to prepare input files for docking by converting them to the PDBQT file format compatible with AutoDock Vina. After comprehensive stereo-chemical investigations of the selected 3D model, the protein was docked to a library of 61 bioactive flavonoid compounds and 5 PTB drugs in PyRx (Ounthaisong and Tangyuenyongwatana, 2017). The grid parameters used for this analysis were X: 20.0265, Y: 35.5676, and Z: 61.2142. The complexes with the most favorable binding energies and RMSD (root mean square deviation) values equal to zero were selected for further pharmacokinetic studies using pkCSM (Pires et al., 2015).

Desmond, a software application from Schrödinger LLC, was used to run 100 nanoseconds of molecular dynamics (MD) simulations. For this purpose, rigid binding of the chosen chemical compound against the target protein was calculated via receptor–ligand docking via MD simulation (Bowers et al., 2006; Ferreira et al., 2015). By using Newton’s classical equation of motion, MD simulation investigations were carried out to predict the ligand binding status in the physiological environment (Hildebrand et al., 2019; Rasheed et al., 2021). Using Maestro’s Protein Preparation Wizard, the chosen proteins and ligands were optimized and reduced. There was no longer any steric conflict, poor contact, or distorted geometry. The systems were constructed using the System Builder tool, and TIP3P (Intermolecular Interaction Potential 3 Points Transferable), an orthorhombic box with the OPLS_2005 force field, was utilized as the solvent model (Shivakumar et al., 2010; Malik et al., 2023). Throughout the simulation period, a 300 K temperature and 1 atm pressure were utilized to simulate physiological conditions, while counterions were added to neutralize the models, and 0.15 M sodium chloride was added. Trajectories were saved after every 100 picoseconds (ps) for inspection, and the stability of the protein–ligand complex was determined by measuring the root mean square deviation (RMSD) over time. Using the Bio3D package of R, principal component analysis (PCA) and a dynamic cross-correlation matrix (DCCM) were examined; however, PCA and DCCM were computed using an R language script (Malik et al., 2023; Ullah et al., 2023).

In generalized molecular mechanics, during MD simulations of the above protein‒ligand complex, the Born surface area (MM-GBSA) module of Prime was utilized to calculate the binding free energy (Gbind) of the docked complex. The binding free energy was calculated using the OPLS 2005 force field, VSGB solvent model and rotamer search methods. After the MD run, the MD trajectory frames were selected at intervals of 10 ns (Ullah et al., 2023). Δ G bind = G complex   –   G protein + G ligand

In this equation, ΔG_bind = binding free energy, G_complex = free energy of the complex, G_protein = free energy of the target protein, and G_ligand = free energy of the ligand.

The thresholds of < −1 and >1 upon log2(fold change) and p-value <0.05 yielded 270 upregulated genes and 119 downregulated genes, respectively, in the case of PE. Some of the top genes included FSTL3, SPAG4, FLT1, LOC105370135, TMEM45A, HTRA4, SLC6A8, and LEP, all of which were overexpressed (Table 1; Supplementary Table S1). However, under the PE + IUGR condition, 262 genes were upregulated, 525 were downregulated, and the top genes included SLAMF1, QPCT, NTRK2, PNCK, ENG, HEXB, MYO7B, LOC102723566, PROCR, CLDN9 and FSTL3 (Table 2; Supplementary Table S2). We also found that none of the genes were upregulated; these genes were common in both datasets, except for PTGS2, and were highly expressed in both datasets. CA1 was the only gene commonly downregulated in both datasets. Interestingly, 20 of the genes whose expression was suppressed in the PE dataset were overexpressed in the PE + IUGR dataset; these genes included RALYL, LOC101928957, OGN, MYO1A, CLDN16, ZNF385B, LOC105377546, LOC107986064, SLAMF1, LOC107985906, RNASE11, SELE, LRRN3, NRAD1, EWSAT1, TRIL, SCN4B, SPP1, ITIH3, and CNTN6 (Figure 2; Figure 3). Figures 4A, B indicates the comparative status (with relation to controls) of both datasets (i.e., PE and PE + IUGR)

Gene set enrichment analysis revealed the molecular functions, biological processes, cellular components, and pathways in which our set of DEGs was enriched under both conditions, that is, placenta with PE and PE + IUGR. The set of PE + IUGR DEGs was mostly enriched in GO terms such as response to hypoxia, regulation of organelle organization, regulation of receptor internalization, cell recognition, synaptic transmission, glutamatergic, positive regulation of type II interferon production, and the JNK cascade and in KEGG pathways such as the apelin signaling pathway and regulation of the actin cytoskeleton. Specifically, enrichment in terms related to “response to hypoxia” could indicate that these DEGs are involved in the cellular response to reduced oxygen levels, which is often associated with complications during pregnancy, including PTB and IUGR. Furthermore, terms such as “regulation of organelle organization” and “regulation of receptor internalization” may imply that these DEGs are involved in maintaining cellular structure and signaling processes, which are critical for fetal development. Additionally, the presence of terms related to the immune response (“positive regulation of type II interferon production”) and cellular signaling pathways (“JNK cascade”) suggested potential mechanisms underlying PTB and IUGR, possibly related to immune regulation and cell signaling abnormalities (Figure 5).

The upregulated KEGG pathways in the case of PE DEGs were related to neuroactive and prolactin signaling pathways, with a set of 6 and 3 genes involved, respectively. The neurodegenerative pathway was downregulated, with a net enrichment score (NES) of −1.82. The top enriched GO terms were steroid biosynthetic process, regulation of endothelial cell proliferation, cellular response to endogenous stimulus, and phosphatidylinositol 3-kinase signaling (Figure 6).

The construction of a PPI network for further analysis of DEGs serves as a highly significant step, as it shows how the proteins encoded by the set of DEGs interact with each other to carry out cellular processes and how disturbing their expression can be problematic in maintaining the balance. CytoHubba analysis of the STRING network of proteins from both datasets helped unveil the hub genes based on the 12 methods. We identified 8 hub genes involved in PE, namely, PTGS2, ENG, KIT, MME, CGA, GAPDH, GPX3 and P4HA1, while the network of the PE + IUGR gene set included 9 hub genes, namely, PTGS2, FGF7, FGF10, IL10, SPP1, MPO, THBS1, CYBB and PF4 (Table 3). MCODE revealed 4 clusters (cluster 1 with a score >3.50) for the PE condition (Table 4) and 7 cluster modules (4 clusters with a score >3.50) for the PE + IUGR condition (Table 4). These clusters help analyze the local subnetworks from the larger networks so that nodes of the highest significance can be identified.

Cluster 1 from the PE dataset included 10 edges and 5 nodes and was particularly significant due to its functional enrichment related to the interaction between chorionic gonadotropin alpha subunit (CGA) and chorionic gonadotropin beta subunit 3 (CGB3) to form human chorionic gonadotropin (hCG) Figures 71A, 1B. When analyzing the significant clusters from the PE + IUGR dataset, we found that cluster 1 proteins were involved in processes such as negative regulation of the response to wounding, the fibroblast growth factor receptor signaling pathway, positive regulation of keratinocyte migration, and lung fibrosis [Figures 72A, 2B]. Cluster 2, containing 6 nodes and 14 edges, was enriched in the sodium channel complex, and Ankyrins linked voltage-gated Na+ and K+ channels to spectrin and L1 [Figures 73A, 3B]. Cluster 3 included 11 nodes, and 25 were associated with processes related to defense mechanisms, such as neutrophil chemotaxis and lymphocyte chemotaxis [Figures 74A, 4B]. Cluster 4 was linked to cellular processes, such as miRNA involvement in the immune response in sepsis, SARS-CoV-2 innate immunity evasion, and negative regulation of T-cell proliferation (Figures 75A, 5B)). The geneMANIA network of the hub genes from both conditions revealed local interactions with other genes in a functional network based on co-expression, shared protein domains, co-localization, and genetic interactions (Figure 8).

PTGS2 was found to be common between both conditions as a highly expressed gene and was among the top hub genes, suggesting that it can be used as a promising target in combating both conditions related to PTB, that is, PE and PE + IUGR. The prostaglandin G/H synthase 2 precursor (PTGS2) consists of 604 amino acid residues. The protein 3D structure is available in the Protein Data Bank under accession number 5f19; however, to fill these gaps, five 3D homology models were generated (Table 5). The selected 3D model of PTGS2 had 91.5% of the residues lying in the most favored regions according to the Ramachandran plot. Verification of 3D showed that 76.04% of the residues had an average 3D-1D score >0.1, while ERRAT showed a 3D model quality score of 95.94%.

Molecular docking analysis of the flavonoid compounds and PTB drugs revealed significant binding affinities for the PTGS2 protein. Overall, 59 out of the 66 flavonoid compounds exhibited good binding affinities in the range of −7.5 to −9.1 kcal/mol, with a root mean square deviation (RMSD) of zero (Supplementary Table S3). Similarly, the PTB drugs allylestrenol, ritodrine, and retinosban also showed good affinities within the favored range, while terbutaline and hydroxyprogesterone caproate had binding energies higher than that of the target protein PTGS2 (Supplementary Table S3). The flavonoid compounds with the most favorable binding energies were flavoxate (−9.1 kcal/mol), disomin (−9 kcal/mol), 7,8-dihydroxyflavone (−8.9 kcal/mol), 3-hydroxyflavone (−8.9 kcal/mol), gossypetin (−8.8 kcal/mol) and rhamnetin (−8.8 kcal/mol). Some of the common interacting residues in the target protein were HIS374, ALA188, TRP373, VAL102, LEU517, ALA513, and HIS372 (Table 6). In addition, the PTB drugs, such as Ritodrine (−8.1 kcal/mol), Hydroxyprogesterone caproate (−6.7 kcal/mol), Retosiban (−8.1 kcal/mol), Terbutaline (−6.8 kcal/mol) and Allylestrenol (−7.8 kcal/mol), had less favorable binding energies (Supplementary Table S3; Figure 9).

The optimal chemical complex of the protein‒ligand complex, i.e., PTGS2 with flavoxate, was simulated using molecular dynamics for 100 ns. Desmond’s simulated trajectories were analyzed. The root-mean-square deviation (RMSD) and root-mean-square fluctuation (RMSF) were determined via MD trajectory analysis.

Figure 101A, 1B shows the time-dependent changes in the RMSD values for the C-alpha atoms in the flavoxate-bound PTGS2 protein. The RMSD plot showed that the protein complex with flavoxate ligands was stabilized throughout the indicated time interval. Once the simulation begins, the RMSD remains between ∼1.5 and 2.5 Å for the remainder of the run. The flavoxate-bound PTGS2 protein complex remained stable over the entire time interval, according to the RMSD values. Throughout the simulation, the ligand fit-to-protein RMSD of the flavoxate-bound PTGS2 protein complex was more stable, that is, it ranged between ∼2.0 and 4.0 Å for the remainder of the run. However, some variations were observed after 60 ns but remained in the range with the least RMSD variation.

Figures 102A, 2B shows the RMSF values for the flavoxate-bound PTGS2 protein complex. Based on the MD trajectories, we found that the overall residues of the protein had stabilized peaks, that is, ∼2.5 Å in the N- and C-terminal zones. The stability of ligand binding to the protein was demonstrated by the low RMSF values of the binding site residues. The flavoxate-bound PTGS2 protein complex included HIS374, VAL281, LEU280, ILE394, GLN189, ALA188, TRP373, and HIS193, which exhibited a greater interaction fraction during the 100 ns analysis. Figures 103B also shows the radius of gyration and analysis of the protein‒ligand interaction area. The PCA eigenvalue was plotted against the percentage of variance, which showed some varying areas that were displayed in three separate sections. Variations in PC1, PC2, and PC3 accounted for 26.88%, 13.05%, and 6.56%, respectively (Figure 103A).

Protein–ligand interaction analysis revealed the active and dominant involvement of hydrophobic interactions with the top key amino acid residues. In addition, water bridge and hydrogen bond formations were observed (Figure 104A). Secondary structural features, such as alpha helices and beta strands, were tracked throughout the simulation (SSE). The percentage of helices was 33.89, the percentage of strands was 4.29%, and the percentage of secondary structure elements (SSEs) was 38.18% (Figure 104B). The ligand properties also summarized the different conformational evolution patterns of flavoxate throughout the simulation trajectory (from zero to 100 ns) (Figure 104C).

The detailed ADMET analysis of the six top-ranked docked candidate flavonoid compounds showed that these compounds were safe for use as drugs (Supplementary Table S4). Small molecules are generally preferred for drug development because they are easier to handle and tend to exhibit better bioavailability. The flavonoid “3-hydroxyflavone” has the lowest molecular weight (238.24), making it a favorable candidate for this purpose. LogP is a measure of the lipophilicity of a molecule and can affect its ability to cross cell membranes. Low logP values are generally preferred. Gossypetin had the lowest LogP (1.6936), indicating good hydrophilicity. The H-acceptor and H-donor parameters are related to the ability of the molecule to form hydrogen bonds. Generally, molecules with moderate numbers of H-acceptors and H-donors are preferred. Diosmin had the highest number of H-acceptors and H-donors. A higher water solubility is desirable for drug candidates because it can enhance their bioavailability. Flavoxate has relatively good water solubility. Caco2 permeability is a measure of a molecule’s ability to pass through intestinal cells. Higher values were generally preferred. Flavoxate has a high Caco2 permeability. Molecules that are not P-glycoprotein substrates or inhibitors may have advantages in terms of drug development. Neither “gossypetin” nor “rhamnetin” are substrates or inhibitors. The volume of distribution in humans can indicate the extent of the distribution of a molecule in the body. Values close to zero were preferable. “3-Hydroxyflavone” has a relatively low VDss value. The unbound fraction parameter is related to the fraction of the drug that is not bound to proteins in the bloodstream. A higher unbound fraction can enhance the efficacy of a drug. “Diosmin” has a relatively high fraction of unbound sites. The blood‒brain barrier (BBB) and central nervous system (CNS) permeability are essential for drugs to target the brain. Molecules that can cross the BBB and exhibit CNS permeability are preferable. None of the molecules in the dataset exhibited strong BBB or CNS permeability. A lower total clearance suggested a longer duration of action of the drug. “7,8-Dihydroxyflavone” had the lowest total clearance value. Molecules with no AMES toxicity are preferred because they are less likely to be mutagenic. None of the tested molecules exhibited AMES toxicity. Molecules that do not inhibit the hERG channel are preferred because hERG inhibition can lead to cardiac side effects. Maximum tolerated dose and toxicity values, as well as acute and chronic toxicity data, should be evaluated to ensure the safety of the molecule as a drug. “Flavoxate” has a relatively high “Maximum Tolerated Dose (humans)". Overall, based on the properties evaluated, “flavoxate” appears to have several favorable characteristics for further analysis as an optimal candidate compound (Supplementary Table S4).