Abstract Rhizomes of Drynaria quercifolia have long been traditionally used to manage rheumatic pain. However, there is limited research supporting this traditional practice and insufficient evidence demonstrating the molecular mechanisms of action of plant-derived bioactives in rheumatoid arthritis (RA). The current study aims to identify the effective components in Drynaria quercifolia methanol rhizome extract (DME) and their probable pharmacological mechanisms in alleviating Rheumatoid Arthritis (RA) using network-pharmacology, molecular docking, molecular-dynamics simulations, and gene expression-based validation. Gas chromatography–mass spectrometry (GC-MS) based screening identified 41 volatile phytocomponents from DME having drug-like potentiality. Network pharmacology-based screening revealed 117 therapeutic targets for RA of which 11 have been identified as core targets. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis indicated that key target genes were mostly enriched in the inflammatory response associated with multiple signalling pathways. Molecular docking and molecular dynamics studies revealed that key target proteins like serine/threonine-protein kinase (AKT1), peroxisome proliferator-activated receptor alpha (PPARA), and peroxisome proliferator-activated receptor gamma (PPARG), exhibited strong binding affinity and stable interactions with multiple phytocomponents present in DME. For experimental verification FCA (Freund’s complete adjuvant)-induced chronic arthritis model employed for further molecular investigation. Quantitative reverse transcription polymerase chain reaction (qRT-PCR) results validated that DME significantly (p ≤ 0.05) regulate the expression of key identified target genes AKT and PPARG in experimental RA model. Moreover, this study further confirmed that DME significantly (p ≤ 0.05) downregulated pro-inflammatory mediators like COX-2, IL-6 and TNF-α at gene and protein levels and also normalized (p ≤ 0.05) different oxidative stress parameters in both the low and high dose groups of DME-treated arthritic animals. In conclusion, the network-based in silico approach indicated that the phytocomponents present in DME probably act in a synergistic way to modulate key identified targets associated with RA, which was further validated by experimental studies. Therefore, DME could be a potential alternative in immunomodulatory therapies to combat RA and related chronic inflammatory conditions. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-025-87461-2. Keywords: Drynaria Quercifolia, Rheumatoid arthritis, Network pharmacology, Inflammation, Cytokine, Molecular docking, Gene expression Subject terms: Computational biology and bioinformatics, Drug discovery Introduction Inflammation is the body’s natural defence mechanism, triggered by pathogen invasion, tissue damage, or chemical disturbances^[42]1. It initiates a self-limiting healing process, but imbalances lead to chronic inflammation, causing cellular and tissue damage and it can aggravate into several chronic inflammatory disorders, including rheumatoid arthritis (RA)^[43]2–[44]4. Rheumatoid arthritis (RA) is the prevalent form of chronic inflammatory arthritis, affecting 1% of the population, characterized by synovial hyperplasia and infiltration of inflammatory cells in the synovial tissues^[45]5,[46]6. Its progression leads to chronic disability and gradual deterioration of cartilage and bone^[47]7,[48]8. Anti-inflammation is a key part in treating RA patients because it is the primary cause of their clinical symptoms^[49]9,[50]10. Conventional corticosteroids, disease-modifying anti-rheumatic medications (DMARDs), such as methotrexate and hydroxychloroquine, non-steroidal antirheumatic drugs (NSAIDs), such as indomethacin, ibuprofen, aspirin, etc., and immunosuppressive agents like TNF-α and IL-6 inhibitors, are presently included in the treatment regime^[51]9,[52]11–[53]15. However, these drugs have significant side-effects, which can include nephrotoxicity, gastrointestinal discomfort, hematological abnormalities, and an elevated risk of cardiovascular diseases^[54]16–[55]19. As a result alternative therapeutic approaches are gaining considerable attention due to their potential for safe and effective management of arthritis^[56]20–[57]23. Drynaria quercifolia (L.) J. Sm, also known as Oak leaf fern or “Pankhiraj” locally in West Bengal (India), is a medicinal pteridophyte, found abundantly in the tropical and subtropical nations like India, Bangladesh, Pakistan, North America, and Africa^[58]24. Traditionally, it is used as “Ashwakatri” in Ayurvedic medicine, with its cooked-rhizome decoction providing anti-pyretic properties, and tuberculosis cure^[59]25. The entire plant is used to treat chest and skin conditions, appetite loss, and rheumatic illness by tribes of Tamil Nadu’s Eastern Ghats^[60]27. The Koch tribe from the Cooch Behar district of northern West Bengal, India, as well as the sub-Himalayan tribal communities like Mech, Toto, and Rabha tribes, use the fronds part of this plant to restore bone fractures and strengthen and heal damaged bones^[61]28,[62]29. In several in vitro and in vivo studies, the efficacy of Drynaria quercifolia rhizome’s effectiveness has already been proven as it has antimicrobial, thrombolytic, analgesic, membrane stabilizing, anti-inflammatory, anti-arthritic, and antioxidant potentials^[63]26,[64]30–[65]33. However, due to the complexity of compounds and their multiple associated targets, it remains challenging to comprehensively elucidate the mechanisms of action of traditional herbal medicines at a holistic level using conventional research methodologies. In the post-genomic era network pharmacology has become a widely used analytical technique that incorporates multi-disciplinary approaches from systems biology, biochemistry, and bioinformatics to rebuild network pharmacology^[66]34–[67]36. This approach evaluates complex compound-protein/gene-disease interactions to create effective remedies^[68]37. Network pharmacology studies focus on multiple goals and diseases, identifying the factors for intercepting therapeutics, such as toxicants, signaling pathways, and disease-ameliorating capabilities of the phytocomponents^[69]37,[70]40–[71]42. Although Drynaria quercifolia’s anti-inflammatory action has been confirmed previously^[72]33, detailed molecular information regarding the plant’s anti-arthritic properties has not yet been elucidated. Therefore, the current study for the first time, explored the molecular basis of Drynaria quercifolia rhizome methanol extract’s (DME) role in RA management through network pharmacology-based experimental validation. The study identified compounds, targets, and protein-protein interaction (PPI) interaction networks, providing a holistic view of DME’s action. The study also investigated the interaction between active ingredients and core target genes, and validated the predicted targets in Freund’s Complete Adjuvant (FCA) induced Wister rat model through mRNA expression and protein level analyses. The workflow of study design is provided in Fig. [73]1. Fig. 1. [74]Fig. 1 [75]Open in a new tab Study workflow diagram of Drynaria quercifolia methanol rhizome extract (DME) against Rheumatoid arthritis (RA). The Drynaria rhizome photograph is captured by D. Modak and was modified with the help of freely distributed Adobe express online tool ([76]https://www.adobe.com/express). Figure was generated using, Microsoft PowerPoint 2021 (Version 2302, Build 16130.20306) 64-bit. Results Identification of active compounds and targets All the identified phytocomponents of DME showed drug-like potentiality based on Lipinski’s rule of Five (Supplementary Table [77]S1). Among them, 41 components reported 792 target proteins from Swiss target prediction (STP) and 561 from the Similarity Ensemble approach (SEA) (Supplementary Table [78]S2) databases, respectively. Venn diagram showed that 280 target proteins were overlapping between the STP and SEA databases (Supplementary Fig. [79]S1a). By combining the Online Mendelian Inheritance in Man (OMIM) and Gene-Disease Association (DisGeNet) databases, a total of 2787 distinct target proteins, linked with RA, were identified (Supplementary Table S3). Finally, among the 280 targets shared between STP and SEA, we identified 117 overlapping common target sets (Supplementary Table S4) that were also present in the OMIM + DisGeNet datasets (Supplementary Fig. [80]S1b). Protein-protein interactions from 117 overlapping target protein PPI network represented proteins as nodes and their interactions as edges. PPI (Fig. [81]2a) network analysis revealed that 114 out of 117 overlapping target proteins closely interacted with each other. However, 3 target proteins (namely Peptidyl-glycine alpha-amidating monooxygenase or PAM, Glutathione S-Transferase Kappa 1 or GSTK1, and Heparanase or HPSE) did not correlate with others in the network. There was a total of 117 nodes (representing active proteins) and 674 edges (representing the interactions between the active proteins and other proteins) in the PPI network, which were significantly more than the expected number of 218 edges (p < 1.0e^− 16). Most of the research on the PPI network was done by finding sub-networks, or modules, that had specific topological characteristics^[82]43. Topology refers to the arrangement of nodes and edges in a network. By analyzing network’s topology, we were able to identify important sub-structures within the network. Results of the network topology analysis indicated that network density was 0.105 and network heterogeneity was 0.967 and the path length was 2.427. The average betweenness centrality of nodes was 0.0124, and there were 33 nodes larger than the average betweenness centrality. Analysis of network topology indicated that the average node degree was 11.5, with 41 nodes exhibiting a higher degree of the average node values. The AKT1 target protein showed PPI’s highest node degree interactions (62) and was regarded as a central target protein (Supplementary Table S6). The core PPI sub-network was evaluated according to the highest degree of centrality (DC) and betweenness centrality (BC) using the Cytoscape plugin CytoNCA^[83]44. The core PPI sub-network comprised of 11 key target proteins (Fig. [84]2b) (Supplementary Table [85]S5). Among the 11 key target proteins, CXCL8 had the highest betweenness value (72.32), followed by AKT1 (60.19), ACE (59.43), ESR1 (52.12), PPARA (43.79), PPARG (38.65), VEGFA (32.73), MMP9 (30.44), TLR4 (30.22), CASP3 (29.85), and PTGS2 (27.33). Fig. 2. [86]Fig. 2 [87]Open in a new tab The PPI network of 117 common target analyses and the Herb-compound-target-pathway network. (a) The node represents the target proteins and the edges represent the protein-protein interactions of 117 common target proteins. (b) The core PPI subnetwork represents the 11 key target proteins. Details of abbreviation of 117 final target proteins and core target proteins are provided in Supplementary Table S4 and Supplementary Table S5 respectively. (c) Herb-compound-target-pathway network. The green colour represents herb, the blue colour represents phytocomponents, the red colour represents disease, brown represents the pathways that were related to the disease, and yellow represents the targets of RA that were directly modulated by the phytocomponents. Detail abbreviations of phytocomponents (blue blocks) and RA target genes (yellow blocks) and are provided in Supplementary Table [88]S1 and Supplementary Table S6 respectively. The image is generated in Cytoscape 3.9.1. ‘Herb-compound-target-pathway network’ analysis Herb-compound-target-pathway (H-C-T-P) network (Fig. [89]2c) illustrates the interrelationship between the compounds and related protein-targets along with the pathways associated with RA. There was a total of 413 nodes and 1891 edges in the network. In this network, the red ellipsoid represented RA, and the green hexagon represented the plant extract. 41 blue hexagon nodes represented the phytocomponents (represented by PubChem IDs), 249 brown circular nodes represented the metabolic pathways related to RA (represented by KO numbers), and 121 yellow hexagon nodes represented the probable protein-targets that phytocomponents could modulate. The edges indicated that nodes can interact with each other. Figure [90]2c illustrated that certain targets were influenced by several phytocomponents, whereas individual phytocomponent might interact with several targets. The higher the number of connected nodes indicated the more crucial roles of the disease targets or the phytocomponents within this network. The lollipop plot (Supplementary Fig. [91]S2), based on the count of interacting phytocomponents and target proteins, indicated that a total of 23 phytocomponents of DME might have interactions with PPARA, which was the highest in number. Another crucial inflammatory marker protein, PTGS2 or Cyclo-oxygenase (COX)-2, showed the second highest interactions and interacted with 21 phytocomponents of DME (Supplementary Fig. S2a). Other inflammatory markers like AKT1, CXCL8, AKT1, IL6ST, MMP-2, and MMP-9 are also modulated by several phytocomponents (Supplementary Table S6). Subsequently, we identified the principal phytocomponents from the H-C-T-P network with multiple targets. Components like 1-phenanthrene-carboxylic acid, 7-ethyl-1,2,3,4,4a,4b,5,6,7,9,10,10a-dodecahydro-1,4a,7-trimethyl, methyl ester, 5-(7a-Isopropenyl-4,5-dimethyl-octahydroinden-4-yl)-3-methyl-pent-2-en- 1-ol, and Beta-d-Ribopyranoside, methyl, 3-acetate were found to have the highest number of interactions with key inflammatory markers like MMP2, PTGES, PTGS2, MMP8, PPARG, ACE, PPARA, IL6ST and MMP9 etc. (Supplementary Fig. S2b). Other phytocomponents like 1,3,4,5-Tetrahydroxycyclo-hexanecarboxylic acid, n-Hexadecanoic acid, Octadecanoic acid, Ergost-5-en-3-ol, (3 beta 24r), Alpha-D-Glucopyranoside, methyl were also found to be significantly active ingredients responsible for modulating the RA pathways (Supplementary Table S7). GO enrichment analysis and KEGG pathway analysis We performed a GO enrichment analysis to elucidate the relevant biological functions of target proteins. The analysis revealed 110 BP, 29 CC, and 39 MF GO entries. The top 20 significantly enriched terms in each GO category (Supplementary Table S8) are shown in Fig. [92]3a. The “positive regulation of transcription” had the highest gene counts (gene counts 28) in the BP categories, followed by “inflammatory response” (gene counts 24). “Cytoplasm”, “plasma membrane”, and “cytosol” were the primary enriched CC categories having 60, 51, and 49 gene counts, respectively. While “zinc ion binding” and “enzyme binding” were the top two ranked categories in MF, having 31 and 21 gene counts, respectively. Fig. 3. [93]Fig. 3 [94]Open in a new tab GO enrichment, Kyoto Encyclopaedia of Genes and Genomes (KEGG) enrichment analysis analysis of 117 targets (p value < 0.05 and FDR < 0.05) and Volcano plot of differential gene expression patten between RA and healthy clinical samples. (a) The top 20 per cent enriched GO terms in “biological process” (BP), “cellular component” (CC), and “molecular function” (MF) are shown. The x-axis represents the GO terms: green for BP, brown for CC, and blue for MF. The y-axis represents the number of genes in each GO term. (b) In KEGG enrichment analysis, the top 23 KEGG pathways are listed by bubble chart according to the parameters “p value < 0.05 and FDR < 0.05”. Each bubble represents different KEGG pathways on the vertical axis, and the proportion of the genes is represented on the horizontal axis. The size of each bubble indicates the number of genes enriched in each KEGG pathway. The larger the bubble, the greater the number of genes involved in the pathway. The colour of each bubble represents the adjusted p-value for each KEGG path. The corrected p-value increases as the bubble become more red. The Figure is generated in SRPlot online toolkit ([95]http://www.bioinformatics.com.cn/srplot). (c) Volcano plot of differential gene expression shows upregulated genes are red in colour, whereas downregulated genes are depicted as blue in colour. KEGG pathway analysis was further carried out to understand the relationship between target proteins and biological pathways. In total, 23 terms (p < 0.05, FDR < 0.05) were obtained from KEGG analysis (Supplementary Table S9). The most significant KEGG targets for chronic inflammation included the neutrophil extracellular trap formation (hsa04613, gene counts 12), arachidonic acid metabolism (hsa00590, gene counts 10), Ras (hsa04014 gene counts 10), phospholipase D (hsa04072, gene counts 8), HIF-1 (hsa04066, gene counts 7), and VEGF signalling (hsa04370, gene counts 6) pathways. In addition to that, we noticed that several other pathways were also significantly enriched. This includes pathways in cancer (hsa05200, gene counts 31), micro RNAs in cancer (hsa05206, gene counts 12), Coronavirus disease-COVID-19 (hsa05171, gene counts 10), estrogen signalling pathway (hsa04915, gene counts 8), AGE-RAGE signalling pathway in diabetic complications (hsa04933, gene counts 7), and regulation of lipolysis in adipocytes (hsa04923, gene counts 5). The bubble plot represented the top 23 KEGG entries selected according to the –logP value (Fig. [96]3b). GEO analysis and identification of differentially expressed genes (DEGs) in human RA synovial tissue samples A differential gene expression study was conducted using a publicly accessible dataset derived from human RA related synovial tissue. The DEGs were identified through fold change analysis (Limma package) with the cutoff criteria (log2FC) > 2 and a p-value < 0.05. A total of 197 DEGs with were identified with 84 upregulated genes and 113 downregulated genes (Supplementary Table S10). The volcano plot (Fig. [97]3c) illustrates the relationship between statistical significance (p- value) and the extent of change (fold change). The plot analysis indicates an upregulation in the expression of concerned genes, including TLR4, PTGS2, VEGFA, MMP9, and ESR1. Although upregulation of CXCL8 was not noticed in this dataset but a linked gene, CXCL2, was identified as upregulated. Conversely, CASP3 exhibited downregulation in this analysis, along with genes like HLA-DQA1 and HLA-DQB1. Molecular docking analysis Molecular docking studies were performed to elucidated the interactions between the active ingredient of DME and RA-related key target proteins at the molecular level. A total of 41 compounds were selected for docking on the 12-target protein under the docking procedure. The best-docked two-dimentional (2D) conformation between the active ingredient and target proteins is represented in Fig. [98]4. The detailed three-dimentional (3D) conformation, binding energies, interacting residues between the active ingredient and target proteins are represented in Supplementary Fig. S3. Docking results of all phytocomponents binding to the key target proteins have been represented in a heatmap in Supplementary Fig. S4. Generally, a low docking score indicates higher binding-interaction between the ligand and the protein. Among the 41 components and the target proteins, 9,12-Octadecadienoic acid (Z, Z)-, 2-hydroxy-1-(hydroxymethyl) ethyl ester showed the highest binding affinity (-10.329 Kcal/Mol) with the target protein PPARA (pdb id-1I7G) and formed four hydrogen (H) bonds at Phenylalanine (PHE) 273, Histidine with hydrogen on the epsilon nitrogen (HIE) 440, Tyrosine (TYR) 464, TYR 314 residues covering the active sites and several hydrophobic interactions (Supplementary Table S11). The same ligand also showed the highest binding energy (-9.845 Kcal/Mol) when docked with PPARG (pdb id-1FM9) and formed three H-bond with Histidine (HIS) 323, Serine (SER) 289, HIE 449 residues and several hydrophobic, and polar interactions (Supplementary Table S11). Fig. 4. [99]Fig. 4 [100]Open in a new tab Molecular docking 2D interaction between best docked active compounds and key targets proteins. (a) 1,3,4,5-tetrahydroxycyclohexanecarboxylic acid binds with ACE; (b) 1,3,4,5-tetrahydroxycyclohexanecarboxylic acid binds with AKT1; (c) 1,3,4,5-tetrahydroxycyclohexanecarboxylic acid binds with CASP3; (d) Ergost-5-en-3-ol, (3.beta.,24r) binds with COX-2; (e) 5-(7a-Isopropenyl-4,5-dimethyl-octahydroinden-4-yl)-3-methyl-pent-2-en- 1-ol binds with ESR1; (f) 1,3,4,5-tetrahydroxycyclohexanecarboxylic acid binds with IL-6; (g) Beta.-d-Ribopyranoside, methyl, 3-acetate bind with MMP-9; (h) 9,12-Octadecadienoic acid (Z, Z)-, 2-hydroxy-1-(hydroxymethyl)ethyl ester binds with PPARA; (i) 9,12-Octadecadienoic acid (Z, Z)-, 2-hydroxy-1-(hydroxymethyl)ethyl ester binds with PPARG; (j) 1,3,4,5-tetrahydroxycyclohexanecarboxylic acid binds with TLR4; (k) Catechol binds with TNF-α; and (l) Catechol binds with VEGFA. All 2D figures are generated in Glide (Maestro Version 12.5.139, Schrödinger). Phytocompounds like 5-(7a-Isopropenyl-4,5-dimethyl-octahydroinden-4-yl)-3-methyl-pent-2-en- 1-ol showed the highest binding affinity (-10.061 Kcal/Mol) with ES1 (pdb id-1A52) and formed two H-bonds at the residues Glutamic acid (GLU) 353 and Arginine (ARG) 394. Docking results depicted that the active ingredient 1,3,4,5-Tetrahydroxycyclo-hexanecarboxylic acid was the most prominent phytocomponent. This ligand showed highest binding affinities with target proteins like ACE (− 7.451 Kcal/Mol; pdb id-1O86), AKT1 (− 7.283 Kcal/Mol; pdb id-1H10), CASP3 (− 7.096 Kcal/Mol; pdb id-1CP3), IL-6 (− 5.091 Kcal/Mol; pdb id-1ALU), and with TLR-4 (− 4.556 Kcal/Mol; pdb id-3FXI) with several H-bonds and hydrophobic interactions (Supplementary Table S11). Ergost-5-en-3-ol showed the highest binding affinity (− 7.178 Kcal/Mol) and formed one H-bond with the target protein COX-2 (pdb id- 5IKQ) at Threonine (THR) 212 residue. It also formed a number of hydrophobic interactions that contribute to its binding energy thus reflecting its COX-2 inhibition activity (Supplementary Table S11). Catechol docked on TNF-α (pdb id- 2AZ5) and VEGFA (pdb id- 1CZ8) showed the highest binding affinity (− 4.3 Kcal/Mol and − 4.159 Kcal/Mol respectively) and formed H-bond and several hydrophobic interactions. Similarly, compound like Beta-d-Ribopyranoside, methyl, 3-acetate showed highest binding energy (-6.228 Kcal/Mol) with three H-bonds at Proline (PRO) 421, Neutral Glutamic acid (GLH) 402, Leucine (LEU) 188 residues formation with target protein MMP-9 (pdb id -1GKC). Molecular dynamics analysis To investigate the intricate molecular interactions, molecular dynamics simulations were conducted utilizing the Schrodinger Desmond package over a 100 ns time scale. This facilitated a thorough analysis of the dynamic characteristics of the selected target proteins interacting with the optimally docked phytocomponents. The average distance between a group of atoms in each frame compared to the reference frame was measured using RMSD (root mean square deviation) to estimate the conformational change from the reference structure. The RMSD plot of AKT1 bound the ligand 1,3,4,5-Tetrahydroxycyclohexanecarboxylic acid (PubChem Id: 1064) presented in Fig. [101]5a reveals the protein and ligand RMSD values were very stable throughout the 100 ns, showing a good interaction. The ligand-bound protein showed significant stability in the major part of the simulation; however, it was slightly unstable around 83 ns of simulation, and after that, again, it stabilised. The RMSD plot (Fig. [102]5b) of PPARA binding to the ligand 9,12-Octadecadienoic acid (Z, Z)-, 2-hydroxy-1-(hydroxymethyl) ethyl ester (PubChem Id: 5365676) established a stabilisation between 12 ns and 30 ns of simulation, and then 36 ns to 80 ns of MD simulation. Fig. 5. [103]Fig. 5 [104]Open in a new tab Root mean square deviation (RMSD) and root mean square fluctuation (RMSF) of the protein-ligand complex. (a) RMSD plot of AKT1 bounded with ligand 1,3,4,5-Tetrahydroxycyclohexanecarboxylic acid (PubChem Id: 1064); (b) RMSD plot of PPARA bounded with ligand 9,12-Octadecadienoic acid (Z, Z), 2-hydroxy-1-(hydroxymethyl) ethyl ester (PubChem Id: 5365676); (c) RMSD plot of PPARAG bounded with ligand 9,12-Octadecadienoic acid (Z, Z)-, 2-hydroxy-1-(hydroxymethyl) ethyl ester (PubChem Id: 5365676); (d) visualisation of the average RMSF plot for AKT; (e) visualisation of the average RMSF plot for PPARA; (f) visualisation of the average RMSF plot for PPRAG. Figures are generated using Desmond Molecular Dynamics System, Schrödinger Release 2020-3. However, the RMSD plot (Fig. [105]5c) of PPARG bound with its ligand 9,12-Octadecadienoic acid (Z, Z)-, 2-hydroxy-1-(hydroxymethyl) ethyl ester (PubChem Id: 5365676) showed instability from 0 ns to 22 ns of simulation. After 23 ns, it became stabilized up to 39 ns. Even though there was a small fluctuation between 40 and 63 ns, the complexes stabilized after that point. In addition, root mean square fluctuations (RMSF) was evaluated and plotted to equate the flexibility of each residue in the ligand-protein complexes. RMSF is useful for characterising local changes along the protein chain. Loop areas (= N and C terminal zones) showed significant peaks, whereas the structural sections (alpha helices and beta strands) typically exhibited fewer variations. The average RMSF values were maintained constantly for all the complexes; the same is shown in Fig. [106]5d-f. Throughout the simulation, the hydrogen bond interaction, hydrophobic bond interactions, ionic bond interactions, and water bridges between the bioactive components and the protein active site residues were also analyzed and depicted in a bar chart (Supplementary Fig. S5a-c). The most significant interactions throughout the simulation in the case of the AKT docked complex were hydrogen bonds in ARG_86 and ASN_54. In the case of simulation of the PPARA complex, HIS_440 and TYR_464 were the most vital H-bond, and PHE_318 and Methionine (MET)_325 were the most critical hydrophobic interaction. In PPARG docked complex, H-bonds in SER_289, HIS_323, TYR_327, TYR_473, and hydrophobic bonds in PHE_282, LEU_330, and Isoleucine (ILE)_341 are some of the important interactions that were noticed. A timeline representation of the interactions and contacts were also studied (Supplementary Fig. S5 d-f). In timeline representations, the number of distinct overall interactions that the protein had with the ligand during the trajectory was displayed in the top panel. The residues that interacted with the ligand in each trajectory frame are displayed in the bottom panel. Moreover, the energy-wise distribution through every 10 ns time intervals of protein-ligand complexes is provided in Additional files (Supplementary Fig. S6, S7 and S8). Results of gene expression patterns of inflammatory biomarkers in FCA-induced chronic RA model After the 28 days feeding experiment, the mRNA expression of several inflammatory biomarkers were assayed. The target genes were selected from our PPI network results and from our docking studies. Following the 2^−ΔΔCt method, it was found that AKT expression was elevated in arthritic control group (3.49 ± 0.51 fold) when compared to vehicle control group. However, treatment with DME high-dose could significantly (p < 0.05) reduce the elevated AKT expression level (1.66 ± 0.47 fold) (Supplementary Table S12) compared to arthritic control group (Fig. [107]6a). A considerable decrease in PPARG and Iκβ expression was observed in the arthritic control group (2.25 ± 0.23 fold and 0.77 ± 0.91 fold respectively) when compared to the vehicle control group. However, treatment with DME increased the PPARG mRNA expressions in low dose treated group (3.63 ± 1.68 fold) and significantly (p < 0.001) increased the PPARG mRNA expressions in high dose treated group (10.4 ± 0.02 fold) (Fig. [108]6b) when compared to the arthritic control group. In the high-dose treated group, Iκβ mRNA expressions also significantly (p < 0.05) increased (2.96 ± 0.38 fold) when compared to arthritic control group (Fig. [109]6c). An elevated expression (expressed as fold change) of COX-2 gene products was noticed in the arthritic control group (1.47 ± 0.14 fold) in comparison to the vehicle control group. However, both the low (0.99 ± 0.07 fold, p < 0.01) and the high (0.38 ± 0.05 fold, p < 0.001) dose groups of DME significantly decreased the elevated COX-2 expression when compared to arthritic control group (Fig. [110]6d). An exaggerated IL-6 expression occurred in the arthritic control group (3.48 ± 0.17 fold). Treatment of arthritic rat with DME both at low (3.20 ± 0.25 fold) and high dose (2.23 ± 0.43 fold) reduced the expression of IL-6 significantly as compared with arthritic control group (Fig. [111]6e). It was found that mRNA expression of TNF-α was raised in the arthritic control group (6.98 ± 0.41 fold) (Supplementary Table S12). However, treatment with DME ameliorated the elevated TNF-α in both the low (3.08 ± 1.99 fold) and high (0.36 ± 0.66 fold) dose groups significantly when compared to the arthritic control group (Fig. [112]6f). Fig. 6. [113]Fig. 6 [114]Open in a new tab Effect of DME on mRNA expressions of various inflammatory markers in experimental arthritic rat groups. Where a represents AKT mRNA fold change, b represents PPARG mRNA fold change, c represents IKβ mRNA fold change, and d represents the COX-2 mRNA fold change, e represents IL-6 mRNA fold change, and f represents TNF-α mRNA fold change in experimental rat groups. All data were represented as mean ± SD and analysed by one-way ANOVA followed by Dunnett’s post hoc test compared with arthritic control group, where *** indicates p ≤ 0.001, ** indicates p ≤ 0.01, and * indicates p ≤ 0.05. Effects of DME on serum concentration of COX-2, TNF α, and IL-6 The effects of DME on the serum concentration of pro-inflammatory markers are presented in Fig. [115]7 and Supplementary Table S13. In arthritic control group, the COX-2 levels (11.6 ± 1.04 ng/ml) showed significant upregulation than that in vehicle control group (3.13 ± 0.46 ng/ml). However, treatment with DME lowered the elevated COX-2 level in both the low dose (6.03 ± 0.287 ng/ml; ) and high dose (5.08 ± 0.505 ng/ml) treated groups significantly when compared to the arthritic control group (Fig. [116]7a). Simultaneously, the arthritic control group had a significantly high TNF-α level (444 ± 28.7 pg/ml). Treatment with DME with both the low dose (235 ± 27.2 pg/ml) and high doses (163 ± 23.6 pg/ml) substantially lowered this high level (Fig. [117]7b). Additionally, the arthritic control group’s IL-6 level was considerably higher at 97.6 ± 17 pg/ml than it was in the vehicle control group (22.6 ± 8.54 pg/ml). However, in both the low dose (67.6 ± 5.84 pg/ml) and high dose (45.3 ± 7.89 pg/ml) DME-treated experimental animals, the DME significantly normalized the elevated IL-6 level (Fig. [118]7c). Fig. 7. [119]Fig. 7 [120]Open in a new tab Effect of DME on serum inflammatory biomarkers and different oxidative stress parameters of liver in different experimental rat groups after treatment schedule. Where, (a) serum COX-2 level; (b) serum TNF-α level; (c) serum IL-6 level; (d) liver SOD activity; (e) liver GSH level; (f) liver MDA level, and (g) liver CAT activity. All data is presented as mean ± SD of six animals and analyzed by One-way ANOVA followed by Dunnett’s post hoc test where ***indicates p ≤ 0.001, and **indicates p ≤ 0.01 as compared to arthritic control group. Effects of DME on oxidative stress parameters The results of oxidative stress triggered by FCA and the alleviating impact of both the dose groups of DME are presented in Fig. [121]7d-g and in Supplementary Table S14. Compared to a vehicle control group (SOD: 10.2 ± 0.30 U/gm tissue, CAT: 18.6 ± 0.97 U/min/gm tissue), there was a significant decrease in SOD (5.90 ± 0.86 U/mg tissue) and CAT activity (9.90 ± 0.45 U/mg tissue) in the arthritic control group as shown in Fig. [122]7d and g respectively. However, DME low dose substantially restored the liver SOD activity (6.71 ± 0.61 U/gm tissue), and significantly normalized the CAT level (12.9 ± 0.79 U/min/gm tissue). Similarly, the high dose of DME significantly normalized the SOD (8.70 ± 0.55 U/gm tissue) and CAT (16.4 ± 0.64 U/min/gm tissue) activity in experimental rat groups. At the same time, the liver GSH level showed a significant (p < 0.05) decrease in the arthritic control group when compared to the vehicle control group (Fig. [123]7e). However, both the low and high dose of DME treatment groups significantly normalized the liver GSH level. Additionally, the arthritic control group’s liver had a considerably higher amount of MDA (48.8 ± 5.89 nmol/gm tissue) than that in the vehicle control group (28.1 ± 2.53 nmol/gm tissue). As shown in Fig. [124]7f, treatment with low and high DME dosages significantly decreased the MDA level in arthritic rats. Discussion In the present network pharmacological analysis, database-based screening of biological targets of 42 volatile active components from DME revealed final set of 117 target genes. Results of the GO enrichment and KEGG analysis of these target genes suggested that they are mainly involved in the regulation of inflammatory response, apoptotic process, regulation of cell proliferation, neutrophil extracellular trap formation, arachidonic acid metabolism, phospholipase D, HIF-1, and VEGF signalling pathways. Multiple signalling pathways like PI3K-Akt, TNF, HIF-1, and VEGF signalling pathways have been found to play a vital role in the therapeutic mechanism of DME against RA. These pathways are triggered by pro-inflammatory cytokines, resulting in immune-mediated inflammation that further leads to an aggravation of RA. In the arachidonic acid metabolism pathway, the PTGS2 gene is overexpressed by the modulation of TNF-α and nuclear factor κ beta (NF-κB) signalling pathways during the inflammatory cascades. Overexpression of the COX-2 enzyme, encoded by PTGS2, in the synovial tissues of RA patients causes the infiltration of inflammatory cells and synovial hyperplasia, resulting in joint degradation^[125]45,[126]46. Uncontrolled proliferation of synovial fibroblasts is characteristic of the pathology of RA. RA synovial fibroblasts (RASF) often exhibit hypoxia, and one of the main regulators of this is a hypoxia-inducible factor (HIF)-1^[127]5,[128]47,[129]48. Along with pro-inflammatory cytokines, synovial cellular infiltration and invasive behaviour are triggered by VEGF, which is a critical regulator of angiogenesis in RASFs under hypoxic conditions^[130]47. Researchers recently discovered that blocking HIF-1 with RNA interference (RNAi) lowers the expression of inflammatory cytokines like IL-6, IL-1, and TNF in peripheral blood serum and synovial cell culture in a collagen-induced arthritis Wistar rat model^[131]49. Therefore, inhibiting the expression of PTGS2 and HIF signalling pathways in synovial cells might be therapeutically beneficial in reducing the inflammatory response in RA. Network pharmacology analysis also highlighted major phytoconstituents of DME having therapeutic effects such as 1,3,4,5-Tetrahydroxycyclo-hexanecarboxylic acid, n-Hexadecanoic acid, Octadecanoic acid, and Ergost-5-en-3-ol, (3 beta 24r) which have multiple targets in multiple pathways. This result supports the findings of Wink et al. (2015) which indicated that the activity of a medicinal plant was due to synergistic interactions of several phytoconstituents^[132]50. Our H-C-T-P network analysis established that DME may intervene in the inflammatory pathways through “multi-component-single target” or “multi-component multi-target” interactions to minimize the release of pro-inflammatory mediators and inflammatory cytokines. This finding further justifies the usage of DME as traditional medicine in the treatment of the chronic inflammatory disease like RA. PPI network analysis revealed AKT1 as a central target protein in RA. It is well known that RA patients’ synovial tissues have abnormally active PI3K/AKT signaling pathways^[133]51,[134]52. Additionally, activation of the PI3K/AKT signaling pathway results in the aberrant production of various downstream effector molecules that have anti-apoptotic and pro-survival effects. AKT also activates the NF-κB pathway, and upon its activation, NF-κB additionally activates a vast spectrum of downstream genes and effectors molecules like COX-2 that further contribute to synovial proliferation and inflammation^[135]53. A previous study showed that curcumin, a plant-derived tetra-terpenoid obtained from Curcuma longa, can inhibit the activation of PI3K/AKT pathways by suppressing the insulin like growth factor or IGF1 (the upstream cytokine of PI3K/AKT pathway) in RA fibroblast-like synoviocyte (RA-FLS) cell line^[136]54. Inhibition of the PI3K/AKT also leads to decrease in TNF-α, IL-6, MMP-2, and MMP-9 production, which suppresses cell proliferation and invasion, and promote apoptosis in cultured RA-FLS^[137]54. From our H-C-T-P network and docking results it was found that the compounds like Gamma-Tocopherol, Hexadecanoic acid, 2-hydroxy-1-(hydroxymethyl) ethyl ester, Benzenepropanoic acid, 4-hydroxy-, methyl ester, 1,3,4,5-Tetrahydroxycyclo-hexanecarboxylic acid can interact with the AKT1 and can inhibit it. Several studies also reported the modulatory role of PI3K/AKT pathway and anti-inflammatory role by these compounds^[138]55–[139]57. Thus, it is plausible to infer from these previous data that the inhibition of AKT1 might significantly attenuate the PI3K/AKT signalling-mediated chronic inflammation progression during RA. Apart from AKT1, the PPI analysis displayed that 10 signalling pathways were directly associated with RA development, suggesting that DME could modulate these pathways against chronic inflammation like RA. The peroxisome proliferator-activated receptors (PPARs) are a group of nuclear receptor proteins that promote ligand-dependent transcription of target genes which can regulate inflammation^[140]58. Most of the anti-inflammatory effects of the PPARs act by inhibiting NF-κB and activator protein (AP)-1, which ultimately repress the expression of downstream genes like PTGS2, inducible nitric oxide synthase (iNOS), IL-6, and IL-12 that are involved in the inflammatory responses^[141]58. PPARα has also been reported to control the duration and magnitude of the inflammatory response through its ability to induce the expression of genes encoding proteins that are involved in the catabolism of pro-inflammatory lipid mediators^[142]59. It was reported that compounds like squalene, benzeneacetic acid, and endogenous fatty acid like n-hexadecanoic acid, and octadecanoic acid act as a potential PPARA agonist in both in vivo and in vitro condition^[143]60–[144]62. Our results indicate that 9,12-Octadecadienoic acid (Z, Z), 2-hydroxy-1-(hydroxymethyl) ethyl ester present in DME possibly interact with PPARA and PPARG in a significant way to control inflammation. Even before the first clinical signs and symptoms appear in RA, activation of the innate system leads to an increase of effector molecules involved in RA-FLS aggression. Specific TLR2, TLR3, and TLR4 are highly expressed in the synovial fluid of RA patients, and their stimulation results in the activation of NF-κB pathway that leads to an increase in the production and upregulation of pro-inflammatory cytokines (such as TNF-α, IL-1, and IL-12), Type I interferon (Type I IFN), and several MMPs within the disease joint^[145]45,[146]46. Our differential expression analyses revealed that TLR4, PTGS2, VEGFA, MMP-9, and ESR1 were notably upregulated in the synovial samples of RA patients. Activated NF-κB regulates genes that contribute to inflammation, including TNF-α, IL-6, IL-8, iNOS, and PTGS2 (also known as COX-2)^[147]63. In the inflamed synovium of RA, unregulated inflammatory condition promotes synovial hyperplasia and angiogenesis^[148]5. ACE (Angiotensin-converting enzyme) is increased in RA-FLS that catalyses the formation of angiotensin II from its inactive precursor, angiotensin I, and then raises the angiotensin II concentration in RA patients that leads to synovial hypoxia by hypoxia-inducible factor (HIF-1α)^[149]64. Studies reported that induction of HIF-1α with TLR signalling also enhances the production of inflammatory cytokines (IL-6, IL-8, TNF-α), MMPs (MMP- 1, MMP-3, MMP-9), and VEGF, thereby escalating the inflammatory state within the RA synovium^[150]47. A positive-feedback regulation between HIF-1α and VEGF triggers angiogenesis in hypoxic conditions^[151]47. Several reports have demonstrated that the synergistic effects of plant extracts can modulate and downregulate inflammatory cytokines like iNOS, TNF-α, IL-6, and even COX-2 in animal models^[152]20,[153]21,[154]65. The apoptotic proteins that trigger the aberrant proliferation of RA synovial fibroblasts (RASFs) can be inhibited by RASFs^[155]5,[156]66. The apoptosis mechanism in osteoclasts of RA joint is regulated by the Bcl-2 protein family^[157]67. They directly contribute to the release of cytochrome c from the mitochondria, which triggers the activation of caspase (CASP)-3 and caspase-9^[158]67. NF-κB, however, controls the overexpression of Bcl proteins in RA and this inhibition could potentially prevent RASFs and macrophages for undergoing apoptosis^[159]68. These changes in RA synovium have been linked to altered apoptotic response of synovial and inflammatory cells^[160]5. Our differential expression data also shows the down regulation of CASP-3 in the synovial samples of RA patients, confirming further the accuracy of the key target proteins identified using PPI network analysis. In summary, by interacting with these targets, active ingredients in DME may contribute to the treatment of RA. To further document the potential effect of the active compounds on the key targets, in silico molecular docking analysis followed by MD simulation was performed. Docking results suggest that 1,3,4,5-Tetrahydroxycyclo-hexanecarboxylic acid showed the best ligand efficacy with interactions with TLR4 (-4.556 kcal/mol), AKT1(-7.283 kcal/mol), ACE (-7.451 kcal/mol), IL-6 (-5.091 kcal/mol), and CASP3 (-7.096 kcal/mol) proteins. Previous studies reported that TLR4 activates AKT1, which leads to the activation of NF-κB^[161]69. NF-κB, a transcription factor that promotes chronic inflammation, is inhibited by derivatives of 1,3,4,5-tetrahy-droxy-cyclo-hexanecarboxylic^[162]70. So, it is possible that 1,3,4,5-tetrahy-droxycy-clo-hexanecarboxylic acid suppressed the TLR4 mediated AKT1 activation to regulate NF-κB activation. Another phytocompound, 9,12-Octadecadienoic acid (Z, Z)-, 2-hydroxy-1-(hydroxymethyl) ethyl ester, also showed a strong binding affinity with PPARA (-10.329 kcal/mol) and PPARG (-9.845 kcal/mol. Most of the anti-inflammatory effects of PPARs result from their inhibition of the NF-κB signaling pathway, which also inhibits the production of various genes implicated in the inflammatory response, including iNOS, IL-6, and IL-12^[163]58. Additionally, inhibition of NF-κB could simultaneously lead to a decrease in COX-2, TNF-α production, and Caspase (CASP)-3 activation, decreasing the level of pro-inflammatory cytokines and promoting apoptosis in the inflamed synovium^[164]63. The protein-ligand docked complexes of AKT1, PPARA, and PPARG were subjected to MD simulation to assess the structural alteration of the complexes under dynamic circumstances. From our simulation study, AKT1 formed 13 hydrogen bonds and among those, Arg 86 and Lys 14 were also indicated in our docking studies. Similarly, PPARA formed 23 hydrogen bonds during the simulation studies, among them three hydrogen bonds (Tyr 314, His 440, and Tyr 464) were similar to docking studies. PPARG also formed 23 hydrogen bonds in the simulation study, and the hydrogen bonds of PPARG were also in line with docking analysis results. All three protein-ligand complexes developed good stability in the target protein’s active site pocket during the course of the simulation, indicating their modulatory actions. The RMSD values were calculated to analyze the structural alterations of the protein-ligand complexes, and both ligand compounds revealed minimal differences. In addition, all three protein-ligand complexes also had reduced RMSF values, indicating the presence of a secondary structure, such as a helix, which denotes the stability of bound structures. Our in silico-based approach identified the key targets, key genes, and key pathways that could be modulated by volatile compounds of DME during RA amelioration. The impact of DME on controlling the expression of crucial gene targets, as identified by network-based investigations, has been further confirmed through experimental validation using qRT-PCR. Our previous study showed that DME can attenuate inflammatory paw-edema, soft tissue swelling and bone erosion in experimental animal groups via modulating the different hematological, biochemical parameters when orally fed to arthritic rats^[165]33. Although stimuli like FCA are thought to activate the NF-κB pathway, they also increase the level of reactive oxygen species (ROS) and inflammatory cytokines, which in turn trigger immune cells to release more inflammatory cytokines and enzymes to aggravate arthritis^[166]71. Furthermore, excessive formation of ROS induces oxidative stress, hence activating pro-inflammatory pathways, including NF-κB and mitogen-activated protein kinase (MAPK), which aggravate inflammation and tissue injury^[167]72,[168]73. According to several reports, PI3K/AKT pathway activation upregulates genes that support vascular permeability, thrombogenicity, and inflammatory cascades by activating the NF-κB pathway^[169]54,[170]74. Moreover, several recent studies on network pharmacology, molecular docking, and in vitro tests have validated that key identified targets from plant extracts like Hedyotis diffusa can modulates TNF-α, IL6, and AKT1, primarily via the PI3K/AKT signalling pathway, establishing a basis for its clinical use in RA treatment^[171]36,[172]37. Meanwhile, previous studies on adjuvant-induced arthritis model have shown that, activation of PPARG expression is linked with the anti-inflammatory activity by inhibiting Iκβ phosphorylation and thus inhibiting NF-κB nuclear translocation^[173]58,[174]75. Similarly, network pharmacology analysis revealed the potential of Gaultheria trichophylla extract and Fangji huangqi decoction in ameliorating RA through identification of the key targets and signalling pathways involved in their anti-RA actions through experimental in vitro and in vivo studies^[175]76,[176]77. Thus our gene expression studies depicted that, in experimental animal models, the active ingredients in DME can inhibit the core target protein AKT and through modulating the expression of PPARG, it can ultimately limit the production of Iκβ. The anti-inflammatory characteristics of PPARA and PPARG also result from their ability to prevent NF-κB translocation in the nucleus, which in turn suppresses the expression of downstream inflammatory response-related genes such COX-2, IL-6, and TNF-α inside the RA synovium (Fig. [177]8). This current study demonstrated that oral administration of DME significantly downregulated the pro-inflammatory cytokines (TNF-α and IL-6) and key inflammatory enzymes like COX-2 mRNA expression and protein level expression in arthritic rats, suggesting the plant’s immunomodulatory properties. Moreover, upon treatment with DME, it significantly restored the oxidative stress parameters. Thus, it is plausible that DME’s reduction of oxidative stress serves as a principal mechanisms for the inhibition of COX-2 and other inflammatory cytokine gene expressions during RA condition. Fig. 8. [178]Fig. 8 [179]Open in a new tab Probable schematic pathway of the mechanism of DME in ameliorating rheumatoid arthritis. TLR toll-like receptor, PI3K phosphoinositide-3-kinase, AKT serine/threonine-protein kinase, PPARA/PPARG peroxisome proliferator-activated receptors Alpha/Gamma, IκB inhibitor of κ beta, NK-κB nuclear factor κ beta, TNF-α tumor necrosis factor-alpha, IL-6 interleukin-6. Figure generated using, Microsoft Word 2021 (Version 2302, Build 16130.20306) 64-bit. This study combined network pharmacology and molecular structure-based approach to elucidate the molecular mechanism of DME and the results suggest that key target proteins related to anti-inflammation can be effectively modulated by multiple active components of DME. One of the probable mechanisms of DME in treatment of RA may be related to the inhibition of AKT1 which may significantly attenuate the PI3K/AKT signalling-mediated chronic inflammation progression during RA (Fig. [180]8). PPARA and PPARG are the other two significant factors which may be modulated by DME. Additionally, mRNA and protein level analysis from our in vivo study supported the anti-inflammatory anti-rheumatoid properties of DME. Methods Collection of plant and identification of phyto-components Fresh rhizomes of the naturally growing Drynaria quercifolia (L.) J. Sm. (Class Polypodiopsida, Order Polypodiales, and Family Polypodiaceae) were collected and authenticated^[181]33. Volatile phyto-compounds identified through GC-MS-based profiling from Drynaria quercifolia rhizome methanol extract (DME) was considered in the current investigation. The details of the procedure can be accessed from our previously published work^[182]33. Out of 47 components, after removal of a few repeated ones, the molecular formula and canonical SMILE (Simplified Molecular-Input Line-Entry System) of 41 identified components were retrieved from NCBI PubChem database ([183]https://pubchem.ncbi.nlm.nih.gov/) and their respective PubChem ID and 2D chemical structures in .SDF (structure data file) format were also obtained^[184]78,[185]79. The detailed name of GC-MS identified components and their PubChem IDs are provided in Supplementary Table [186]S1. Active phyto-component screening Only a small percentage of the phyto-components found in traditional solvent extracts of plant specimens can have medicinal benefits. The physicochemical characteristics of the phyto-components were examined to screen out active components having drug-like potentiality. The identification of drug-likeness qualities of phyto-components is based on Lipinski’s Rule of Five and was checked using the QikProp module of Maestro (QikProp, Maestro Version 12.5.139, Schrödinger, LLC, New York, USA)^[187]80. Target Prediction by database screening The SMILE format of selected bioactive phyto-components confirmed by Lipinski’s rule was further used for potential gene target identification^[188]80. Two public cheminformatics platforms were used for this purposes, Swiss Target Prediction (STP, [189]http://www.swisstargetprediction.ch/ 2022-03-16) and Similarity Ensemble Approach (SEA, [190]https://sea.bkslab.org/ 2022-03-22) implemented in “Homo sapiens” mode^[191]81,[192]82. The overlapping targets between STP and SEA were further selected. Simultaneously, Online Mendelian Inheritance in Man database (OMIM, [193]https://www.ncbi.nlm.nih.gov/omim 2022-03-27) and Gene-Disease Association database (DisGenet, [194]https://www.disgenet.org/ 2022-03-29) were used to identify the target proteins associated with Rheumatoid Arthritis (RA) and after eliminating the repeated targets, we then acquired all RA targets and merged into a single disease-specific protein target set^[195]83,[196]84. The final common set of targets between the selected phytocomponent-related targets and RA-disease-associated targets was visualized on the Venn diagram plotter web tool ([197]https://bioinformatics.psb.ugent.be/webtools/Venn/ 2022-03-29). Construction of PPI network and extraction of the core PPI sub-networks The final target protein set was utilized to construct protein-protein interactions (PPI) on Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) ([198]https://string-db.org/) and defined the gene source as “Homo sapiens”^[199]85. The confidence level was adjusted to 0.4 to determine the relationship between protein interactions. The result of the PPI network was exported to Cytoscape 3.9.1 software for further visualization^[200]86. The network nodes represent proteins, while the edges represent protein-protein associations. The core PPI subnetwork was extracted using CytoNCA, a Cytoscape plugin for network centrality analysis^[201]44,[202]86. First, genes with the highest degree centrality (DC) values in the top 30% were selected for subnetwork construction using CytoNCA^[203]44. Then, genes with the top 30% highest value of betweenness centrality (BC) in the subnetwork were identified as key target genes and formed the core PPI network. Construction of Herb-Compound-Target-Pathway network The herb-compound-target-pathway network was constructed using Cytoscape 3.9.1, using the merged network approach to gather a broad overview of the interactions of phytocompounds with final target proteins along with the disease-associated pathways^[204]86. Gene enrichment analysis To determine the biological functions of target proteins and RA-related pathways Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was performed using The Database for Annotation, Visualization, and Integrated Discovery (DAVID) tools version 6.8 ([205]https://david.ncifcrf.gov/).^[206]87,[207]88 The select identifier was limited to “OFFICIAL GENE SYMBOL,” the list type was limited to “Gene List,” and the species was limited to “Homo sapiens” p-value < 0.05 and FDR < 0.05 were set as the parameters for further visualization. GO enrichment analysis interprets the biological function of target genes in terms of gene function. GO has three ontologies: molecular function (MF), cellular component (CC), and biological process (BP). KEGG pathway enrichment analysis identified significantly enriched metabolic pathways involving targets compared with the whole genome background^[208]89,[209]90. The SRPlot online toolkit ([210]http://www.bioinformatics.com.cn/srplot 2022-04-02) was used to represent the findings of the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment studies^[211]42. GEO analysis and identification of differentially expressed genes (DEGs) The gene expression profiles based on microarray experiments associated with RA were obtained from the Gene Expression Omnibus (GEO) database^[212]91. Following the screening, microarray data submitted by Woetzel D et al. (GEO accession number [213]GSE55235) conducted on the experimental platform [214]GPL96 (Affymetrix Human Genome U133A Array) was chosen for the subsequent differential gene expression (DEG) analysis^[215]92. The dataset consisted of 30 synovial tissue samples, comprising 10 from healthy joints, 10 from osteoarthritic joints, and 10 from rheumatoid arthritis-affected joints^[216]92. For the DEG analysis, only samples related to healthy and rheumatoid arthritis patients were further considered. The objective was to examine whether core targets determined from protein-protein interaction (PPI) network exhibited differential gene expression between healthy and rheumatoid arthritis human data set. To normalize the expression datasets, the robust multichip averaging (RMA) method was applied^[217]93. The R limma package was utilized to identify differentially expressed genes (DEGs)^[218]94. Finally, a volcano plot was generated to visually represent the differential expression of genes. Molecular docking study Core target proteins identified from PPI network along with two previously reported biomarkers, namely IL-6 and TNF-α, were selected for molecular docking analysis. IL- 6 and TNF-α were selected based on the pivotal role they play in chronic inflammatory diseases like RA^[219]20,[220]22,[221]95. The 3D structures of proteins were retrieved in pdb format from Protein Data Bank^[222]96. The proteins were prepared using Schrödinger (Maestro Version 12.5.139, Schrödinger, LLC, New York), utilizing the force field Optimized Potentials for Liquid Simulations (OPLS)-3. A total of 41 phytochemicals identified in GC-MS analysis were considered as ligands for the docking studies. NCBI PubChem ([223]https://pubchem.ncbi.nlm.nih.gov) database was used to retrieve the structures of a selected phytochemical in SDF format^[224]79. The ligands were prepared using the force field OPLS3 in LigPrep (Maestro Version 12.5.139, Schrödinger, LLC, New York). For executing molecular docking analysis, we first generated receptor grids for each protein in Glide (Maestro Version 12.5.139, Schrödinger, LLC, New York). The grids were generated with specific coordinates for the ligands and are provided in Supplementary Table S15. The binding affinities in terms of docking scores for each ligand were recorded using Glide’s Extra-precision (XP) docking models. Molecular dynamics simulation study Docked complexes of AKT (AKT Serine/Threonine Kinase), PPARA (Peroxisome Proliferator Activated Receptor Alpha), and PPARG (Peroxisome Proliferator Activated Receptor Gamma) were further subjected to MD simulation. Desmond (Schrödinger Release 2020-3: Desmond Molecular Dynamics System, D. E. Shaw Research, New York, NY, 2020) software was used for MD simulation of best scoring receptor-ligand docked complexes to assess their dynamic behavior. Protein Preparation Wizard of Maestro (Schrödinger LLC, New York, NY, USA) was utilized to pre-process the protein-ligand complex. Initially, water molecules and other cofactors were removed from the crystal structure, followed by the addition of hydrogen atoms and missing side-chain residues. The simulation orthorhombic box was generated using the System Builder module and solvated with a TIP3P (Intermolecular Interaction Potential 3 Points Transferable) water model with a periodic boundary condition for a 10 Å buffer region (10 Å×10 Å×10 Å) in OPLS 2005 force^[225]97,[226]98. A salt concentration of 0.15 M was added, and the system’s overall charge was neutralized by adding chloride and sodium ions^[227]99. The system was minimized using the steepest descent technique for 1000 iterations. After equilibration, the unrestrained production phase ran under isothermal–isobaric (NPT) ensemble (number of atoms, pressure, and temperature were kept constant) at 300 K temperature and 1.01325 bar pressure^[228]100. Finally, 100 ns MD simulation was performed. The system’s stability was assessed using the Root Mean Square Deviation (RMSD), Root Mean Square Fluctuations (RMSF), and hydrogen bond analysis by using the MD trajectory. In vivo anti-arthritic study Animal maintenance We procured Wistar albino male rats from M/s Chakraborty Enterprise, Kolkata (erstwhile Calcutta), India (Regd. No. 1443/PO/Bt/s/11/CPCSEA), which is an authorized animal dealer. The rats were 8–12 weeks old and weighed 120 ± 10 g. All animals were kept in well-ventilated polypropylene cages (Tarson, India) with rice husk utilized as bedding materials. All the experimental rats were kept at the animal house of the departmental facility, in a controlled environment with a 12-hour cycle of light and darkness and a constant room temperature of 25 °C ± 3 °C. All the experimental animals received unrestrained access to food and water ad libitum. Before the study, animals had a two-week acclimatization period. The experimental procedures were carried out from 2022 to 2023 in strict compliance with the ethical guidelines approved by the Institutional Animal Ethical Committee. The ARRIVE (Animal Research: Reporting of In Vivo Experiments) guidelines were followed in reporting the study^[229]101. Extract preparation and dose selection According to a previously described procedure, methanol extract of the D. quercifolia rhizome methanol extract (DME) was prepared^[230]33. The DME was stored at 4 °C in an airtight container for future experimental usage. The oral feeding-dose was chosen for in vivo investigations based on previously published acute oral toxicity tests^[231]33. Experimental design The adjuvant-induced RA model established in this study was similar to that described in previously published protocols^[232]20,[233]33,[234]102. A total number of twenty-four (24) Wistar albino male rats (120 ± 10 g) were randomly divided into four groups of six rats each as under: Group 1 (Arthritic control): Arthritic rats received orally normal saline (p.o.) only, from day 0 to day 28. Group 2 (Low dose group of DME): Rats were treated with DME (250 mg/kg b.w.). The therapy was commenced on day 0 and continued till the 28th day. Group 3 (High dose group of DME): Rats were doses with DME (500 mg/kg b.w.). The therapy was commenced on day 0 and continued till the 28th day. Group 4 (Vehicle control): 0.5% of DMSO was administered (p. o.) in each healthy rat of the vehicle group. The duration of treatment was similar to other experimental groups. A 0.1 ml of FCA was injected intradermally into the right hind-paw sub-plantar region of each rat to generate arthritis. All animals in the groups received an injection of FCA (Sigma-Aldrich, India) on day zero, except the animals in the vehicle control group. On day 14 of the experimentation, the animals were given a booster dosage of 0.1 ml of FCA. The course of treatment started on the first day of arthritic induction and continued through the day 28th. According to the experimental groups, DME doses were administered once daily orally during the treatment schedule. After the treatment schedule, all animals were anesthetized with sodium pentobarbital (60 mg/kg; i.p) and euthanized by cervical decapitation. Thereafter, blood samples were obtained through the cardiac puncture into EDTA coated vials and serum samples were collected for biochemical analysis as per our previously published protocol^[235]33. Quantitative real‑time polymerase chain reaction (qRT‑PCR) Total RNA was isolated using Trizol LS (Invitrogen, USA) following the manufacturer’s protocol. Before preparing the cDNAs, all the RNA samples were quantified in nanodrop spectrophotometer (SPECTROstar Nano, BMG LABTECH, Germany) and subjected to DNase (New England Biolabs, USA) treatment following the manufacturer’s protocol. The total RNA from the animals of each group was pooled together and cDNA was prepared using reverse transcriptase enzyme (Thermo Fischer, USA; Hi-cDNA Synthesis Kit, India), random hexamers, and dNTP (Sigma Aldrich, USA) with some modifications following a published protocol^[236]20. The prepared cDNA was used to quantify the expressions of AKT, PPARG, COX-2, IL-6, TNF-α, and IKβ gene in a fold-change manner^[237]20. For, PPARG, COX-2, TNF-α and Iĸβ, the gene expression study was done through real time quantitative PCR (qRT-PCR) in Roche LightCycler 96 (Roche, Switzerland) machine using FastStart Universal Sybr Green Master dye (Roche Diagnostics, USA). For AKT and IL-6 gene, the whole study was performed in Bio-Rad CFX96 Touch Real Time Detection system (Bio-Rad, USA). The specific primers of all the genes were selected from existing research papers after validating through primer blast and nucleotide blast to ensure best binding affinities. An internal control housekeeping gene GAPDH was also included for estimating the relative expression of genes using 2^−ΔΔCT method^[238]103. The details list of primers is provided in Supplementary Table S16. Estimation of serum level COX-2, TNFα, and IL-6 COX-2 was estimated by ELISA kits (provided by MyBioSource, USA), and TNF-α, and IL-6 were estimated by ELISA kits (provided by RayBiotech, USA) following the manufactures’ protocol. Estimation of oxidative stress biomarkers in liver After sacrifice, the liver tissue was removed immediately from all experimental animal groups and 10% w/v liver homogenate was prepared followed by standard protocol^[239]104. From this homogenate the activities of superoxide dismutase (SOD), reduced glutathione (GSH), malondialdehyde (MDA) content and catalase (CAT) were spectrophotometrically determined following established protocols^[240]104,[241]105. Statistical analysis The results were expressed as mean ± standard deviation (S.D.). All data were analyzed using analysis of variance (ANOVA) followed by Dunnett’s multiple comparisons test. Values of p ≤ 0.05 were considered to indicate a statistical difference. All the statistical analyses were performed using Graph Pad Prism Version 7.00 for Windows (GraphPad Software Inc., San Diego, USA). Conclusion In summary, this study investigated the therapeutic potential and mechanism of action of DME on RA through a combined approach involving network pharmacology, molecular docking and molecular dynamics study and experimental validation through qRT-PCR and ELISA. Network-based investigation reveals that the DME may influence different targets and multiple pathways associated with RA. The phytocompounds probably act in a synergistic mechanism to modulate key targets of inflammatory pathways associated with rheumatoid arthritic conditions. Our in vivo findings indicated that DME inhibits inflammatory cytokines production in RA via PI3K/Akt, Iκβ- NF-κB and PPARG signaling pathways. Thus, through the inhibition the expression of pro-inflammatory mediators such as COX-2, TNF-α, and IL-6, both at the gene and protein levels, DME administration also inhibits the progression of arthritis. Thereby, the findings of this current study offered a novel approach to further explore the potential molecular pathways through which the rhizome methanol extract of D. quercifolia demonstrates its anti-RA and anti-inflammatory properties. However, it is essential to perform a thorough screening and identification of the key phytochemicals with anti-inflammatory properties, including their isolation, purification, and characterization. Moreover, a comprehensive examination of the impact of phytocomponents on the gene-protein expression levels of additional inflammatory and toxicological biomarkers is necessary to refine and strengthen the findings of this study. Electronic supplementary material Below is the link to the electronic supplementary material. [242]Supplementary Material 1^ (9.7MB, docx) Acknowledgements