Abstract Sanbi Decoction (SBD) has demonstrated promising therapeutic potential in osteoarthritis (OA) treatment, yet its precise mechanisms remain unclear. This research combined computational and experimental approaches, including bioinformatics analysis, network pharmacology, molecular docking, molecular dynamics simulations, and laboratory validation, to investigate the mechanisms of action of SBD. A total of 114 active compounds and 113 intersecting targets were identified through TCMSP and multiple screening strategies. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses revealed that these targets are primarily involved in key signaling pathways, including the AGE-RAGE signaling pathway, IL-17 signaling pathway, and TNF signaling pathway. Among the active components, Shinflavanone, Gancaonin L, Xambioona, Phaseol, Gancaonin O, and Licoisoflavanone exhibited strong binding affinity and structural stability with core targets, as validated by molecular docking and molecular dynamics simulations. Experimental results confirmed that SBD alleviates oxidative stress, reduces inflammation, and protects cartilage by inhibiting the AGE-RAGE/JNK pathway. These findings highlight SBD’s potential as a promising therapeutic agent for OA treatment. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-025-99055-z. Keywords: Sanbi Decoction, Osteoarthritis, Network pharmacology, Molecular docking, Molecular dynamics simulations. Subject terms: Gene ontology, Gene regulatory networks, Virtual drug screening Introduction Osteoarthritis (OA) is a prevalent degenerative joint disease characterized by dysfunction, pain, and stiffness^[36]1. Synovitis is a key pathological feature of OA, which severely impacts patients’ quality of life^[37]2. According to the Global Burden of Disease study, OA is a major public health concern, with 607 million people affected worldwide in 2021^[38]3. The dual physical and psychological burden of OA imposes significant challenges on families and society. Traditional Chinese Medicine (TCM) has demonstrated unique advantages in the treatment of OA due to its efficacy and minimal side effects, contributing significantly to global health^[39]4,[40]5. TCM approaches to treating OA are multifaceted, focusing on promoting blood circulation, nourishing the liver and kidneys, and improving joint function^[41]6,[42]7. Sanbi Decoction (SBD) is a traditional herbal formulation created by Dr. Li Baochao, a distinguished practitioner of Traditional Chinese Medicine (TCM) from Luoyang, China. This formula consists of eight herbs: Radix Saposhnikoviae (Fangfeng), Notopterygium Root (Qianghuo), Radix Gentianae Macrophyllae (Qinjiao), Coix Seed (Yiyiren), Angelica Sinensis (Danggui), Radix Aconiti (Chuanwu), Radix Aconiti Kusnezoffii (Caowu), and Glycyrrhizae Radix (Gancao). SBD is known for its effects in dispelling wind and dampness and strengthening tendons and bones. While its components or their compounds have shown potential for treating OA, the specific mechanisms and targets remain unclear, necessitating more comprehensive insights. Embracing novel approaches and fresh insights is essential for gaining a more comprehensive understanding of the anti-osteoarthritis mechanisms of Traditional Chinese Medicine formulations. In recent years, with the rapid development of bioinformatics, network pharmacology combined with computational simulations has emerged as a promising approach, laying the foundation for comprehensively understanding the mechanisms of TCM in disease treatment^[43]8–[44]10. As a computational technique, molecular docking predicts the binding affinity and interaction patterns between small-molecule ligands from Traditional Chinese Medicine and protein receptors. This approach offers valuable insights into drug-target interactions and supports the pharmacological assessment of both TCM formulations and individual herbs^[45]11. However, molecular docking lacks detailed assessments of interaction stability and dynamics. Molecular dynamics (MD) simulations play a critical role in bioinformatics, particularly in analyzing ligand-receptor interactions^[46]12. These simulations provide in-depth insights into the dynamic behavior of biomolecules, uncovering key information about binding mechanisms, stability, and conformational changes. For instance, MD studies on G protein-coupled receptors (GPCRs) have revealed important insights into their conformational dynamics, essential for understanding their role as drug targets^[47]13. Similarly, MD simulations have been pivotal in uncovering hidden binding sites and predicting drug-binding poses when targeting small GTPases, crucial for developing novel inhibitors^[48]14. The integration of network pharmacology, molecular docking, MD simulations, and experimental validation enhances the reliability of results^[49]15,[50]16. Therefore, employing these methodologies to identify the active components and targets of SBD holds great potential for elucidating its anti-OA mechanisms (Fig. [51]1). Fig. 1. [52]Fig. 1 [53]Open in a new tab Flow chart of the study strategy. Materials and methods Acquisition and organization of GEO data Gene expression profile datasets [54]GSE55235 and [55]GSE55457 were retrieved from the GEO (Gene Expression Omnibus, [56]https://www.ncbi.nlm.nih.gov/geo/) database using “osteoarthritis” as the keyword. These datasets encompass microarray data from 50 samples. Data analysis using Perl and R Background correction and matrix normalization were performed in the R environment (version 4.4.1: [57]https://www.r-project.org/) based on the characteristics of the data samples. The Limma R package (version 3.62.1, [58]https://bioconductor.org/packages/limma/) was utilized to analyze differential gene expression in the microarray data. Significant differentially expressed genes (DEGs) were identified using filtering criteria of an adjusted p-value < 0.05 and an absolute value of log fold change (| log FC |) ≥ 1. Volcano plots for the microarray data were generated using the ggplot2 package (version 3.5.1, [59]https://ggplot2.tidyverse.org/), and heatmaps of the gene chips were created using the pheatmap package (version 1.0.12, [60]https://cran.r-project.org/package=pheatmap). Selection of SBD active components and targets The herbal composition of SBD was input into the TCMSP (Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform, [61]https://www.91tcmsp.com/#/database). Active components were filtered using thresholds of oral bioavailability (OB) ≥ 30% and drug-likeness (DL) ≥ 0.18^[62]17,[63]18, along with their corresponding targets. The UniProt (Universal Protein Resource, [64]https://www.uniprot.org/) was used for target validation by confirming gene names and UniProt IDs. Non-human genes were excluded to ensure the precise identification of potential targets of SBD active components. Identification of potential disease targets For osteoarthritis-related research, target information was collected from GeneCards (GeneCards^® Human Gene Database, [65]https://www.genecards.org/), DisGeNET (Disease-Gene Association Network, [66]https://www.disgenet.org/), the Online Mendelian Inheritance in Man database (OMIM: [67]https://www.omim.org/), and the Therapeutic Target Database (TTD: [68]https://db.idrblab.net/ttd/). These targets were merged and deduplicated using the UniProt database. The results were then combined with differentially expressed genes (DEGs) from the GEO database, with duplicate entries removed to establish a comprehensive disease target dataset. Construction of protein-protein interaction network and core target screening The intersection of potential targets for SBD and OA was identified using the Venn platform, generating a Venn diagram. A protein-protein interaction (PPI) network for intersecting genes was constructed using the STRING (Search Tool for the Retrieval of Interacting Genes/Proteins, [69]https://cn.string-db.org/), with “Homo sapiens” as the target species and a confidence threshold of “high confidence” (≥ 0.700)^[70]19,[71]20. Finally, Cytoscape (version 3.10.2, [72]https://cytoscape.org/) and the CytoHubba plugin were employed to analyze three key metrics: degree centrality (Degree), maximal clique centrality (MCC), and network centrality (MNC). The six highest-ranked genes were selected to establish the core target network. Construction of the drug-active component-common target network A network diagram was generated using Cytoscape to depict the interactions among drugs, active components, and shared targets. This visualization enhances the understanding of compound-target relationships and their potential mechanisms of action. GO and KEGG enrichment analysis GO (Gene Ontology, [73]http://geneontology.org/) and KEGG (Kyoto Encyclopedia of Genes and Genomes, [74]https://www.genome.jp/kegg/) enrichment analyses were conducted in the R environment. The org.Hs.eg.db package (version 3.2.0, [75]https://bioconductor.org/packages/org.Hs.eg.db/) was employed for gene annotation, while the clusterProfiler package (version 4.14.4, [76]https://bioconductor.org/packages/clusterProfiler/) facilitated gene ID conversion and pathway enrichment analysis of the identified common targets^[77]21. These methods accurately identified GO terms related to biological processes, molecular functions, and cellular components associated with SBD. KEGG pathway enrichment analysis was also performed^[78]22. The results were visualized using the ggplot2 package, offering a clear presentation of enriched GO terms and KEGG pathways to support further research. Target-KEGG network construction A network illustrating interactions between common targets and KEGG pathways was constructed using Cytoscape. The “Network Analyzer” module was used to analyze network properties. Targets meeting three critical indicators—Degree, MCC, and MNC—were identified as core targets. Molecular docking simulation The full names of key target proteins were verified using the UniProt database, and their corresponding three-dimensional (3D) structures were obtained from the PDB (Protein Data Bank, [79]https://www.rcsb.org/) database. The selected protein structures included IL1B (PDB ID: 8RYS), IL6 (PDB ID: 1IL6), TNF (PDB ID: 1TNF), JNK (PDB ID: 4QTD), and STAT3 (PDB ID: 6NJS). Using PyMOL (version 3.1, [80]https://pymol.org/), water molecules were removed from the protein structures. The three-dimensional (3D) molecular structures of key active compounds were retrieved from the PubChem and TCMSP databases. Using OpenBabel (version 2.3.2, [81]https://openbabel.org/), these structures were converted from SDF format to mol2 format for molecular docking analysis. Hydrogen atoms were added, and charges were calculated with AutoDockTools (version 1.5.7, [82]http://autodock.scripps.edu/resources/adt) to prepare the proteins for molecular docking^[83]23. After target annotation, docking was conducted using AutoDockTools^[84]24. Binding energy heatmaps were generated using the ggplot2 package in R to illustrate the docking affinities. To validate the RMSD values, ensure that the RMSD is less than 2.0 Å to confirm the reliability of the docking results. Two-dimensional interaction diagrams were generated using Discovery Studio Visualizer (version 19.1, ​[85]https://discover.3ds.com/discovery-studio-visualizer-download)^[86 ]25,[87]26. Finally, 3D images were generated using PyMOL. Molecular dynamics simulation Based on the molecular docking results from AutoDock, the top-ranked model was selected as the starting point for molecular dynamics (MD) simulations. The Amber99sb force field was used to characterize proteins, while the GAFF force field was applied to ligands. Topology files for each ligand were generated using the Sobtop tool (v.3.1, [88]https://www.sobtop.com/). MD simulations were conducted using GROMACS (version 2020.6, [89]http://www.gromacs.org/)^[90]27, with each complex solvated in an SPC water model and neutralized with sodium and chloride ions. Energy minimization was conducted to optimize the system, followed by a gradual temperature increase to 300 K over 200 ns. Subsequently, a 2-ns NPT equilibration was performed at 1 bar to stabilize the system. Finally, a 200-ns production molecular dynamics simulation was executed under controlled conditions of 300 K and 1 bar, with a timestep of 2 fs. The cutoff distance for van der Waals and short-range electrostatic interactions was set at 1.0 nm. The LINCS algorithm was used to constrain hydrogen bonds, with V-rescale and Berendsen methods employed to maintain stable temperature and pressure conditions. During the MD simulation, several parameters were analyzed to assess the conformational stability of the protein-ligand complex, including root mean square deviation (RMSD), root mean square fluctuation (RMSF), solvent-accessible surface area (SASA), radius of gyration (Rg), hydrogen bonds (HB), and principal component analysis (PCA). The MM/PBSA (Molecular Mechanics/Poisson-Boltzmann Surface Area) method was applied to calculate the binding free energy, providing insights into the thermodynamic stability of the complex. Drug preparation and cell culture The eight herbal components of SBD were ground into a fine powder in proportion and dissolved in distilled water at a concentration of 1 g/mL, The mixture was then centrifuged at 5000 rpm for 20 min. The supernatant was collected and filtered using a 0.22 μm microporous membrane to obtain the 1 g/mL crude extract, which was stored at − 80 °C. Murine macrophage-like RAW264.7 cells (Catalog No. JW-CL-0582) were purchased from Shanghai Jiwei Biotechnology Co., Ltd. The cells were cultured in DMEM high-glucose medium supplemented with 10% fetal bovine serum (FBS) under standard conditions of 37 °C with 5% CO₂. CCK-8 cell viability assay Cell viability was assessed using the Cell Counting Kit-8 (CCK-8, Catalog No. CK04-500). Briefly, RAW264.7 cells were seeded into 96-well plates at a density of 7.5 × 10³ cells/well and incubated for 24 h. The cells were then treated with SBD at concentrations of 12.5, 25, 50, 100, 200, and 400 µg/mL for 24 h. After treatment, 10µL of CCK-8 reagent was added to each well, followed by incubation at 37 °C for 1 h. The optical density (OD) was measured at 450 nm using a microplate reader. Cell viability was calculated using the standard formula, and the experimental data were visualized using GraphPad Prism 9.5. Real-time quantitative PCR (RT-qPCR) analysis of gene expression RAW264.7 cells were seeded into 6-well plates and incubated for 24 h. Cells in the control group were treated with phosphate-buffered saline (PBS), while the OA group and SBD group were stimulated with LPS (100 ng/mL) for 2 h. Subsequently, the SBD group was treated with SBD for 24 h. After treatment, the culture medium was discarded, and the cells were washed with PBS. Total RNA was extracted using an RNA extraction kit, followed by reverse transcription into cDNA using a reverse transcription kit. qPCR was performed using cDNA as a template on the LightCycler^®96 system with the following primer sequences: IL1-β: Forward: 5′-TGCCACCTTTTGACAGTGATG-3′. Reverse: 5′-ATGTGCTGCTGCGAGATTTG-3′. IL-6: Forward: 5′-GGAGCCCACCAAGAACGATAG-3′. Reverse: 5′-GTGAAGTAGGGAAGGCCGTG-3′. TNF: Forward: 5′-CCCTCACACTCAGATCATCTTCT-3′. Reverse: 5′-GCTACGACGTGGGCTACAG-3′. STAT3: Forward: 5′-CCGCAGCTTGGGCTGGAAGA-3′. Reverse: 5′-CAGGGCCGGGCTGTGGTAGT-3′. JNK: Forward: 5′-GTGGAATCAAGCACCTTCACT-3′. Reverse: 5′-TCCTCGCCAGTCCAAAATCAA-3′. GAPDH (Internal control): Forward: 5′-TGGATTTGGACGCATTGGTC-3′. Reverse: 5′-TTTGCACTGGTACGTGTTGAT-3′. All primers were synthesized by Nanning GenSys Biotechnology Co., Ltd. The 2−∆∆Ct method was used to calculate the fold change in gene expression, normalized to GAPDH. Experimental data were visualized using GraphPad Prism 9.5 software. Detection of ROS scavenging ability using the DCFH probe RAW264.7 cells were seeded into 12-well plates and incubated for 24 h. The grouping and treatment methods were the same as in Section “[91]Real-time quantitative PCR (RT-qPCR) analysis of gene expression”. After treatment, the cells were washed with phosphate-buffered saline (PBS) and incubated with the DCFH-DA fluorescent probe. Fluorescence imaging was performed using a FITC channel, and fluorescence intensity was quantified using ImageJ software, normalized, and visualized using GraphPad Prism 9.5 software. Immunofluorescence staining RAW264.7 cells were seeded into 12-well plates and incubated for 24 h. The grouping and treatment methods were the same as in Section “[92]Real-time quantitative PCR (RT-qPCR) analysis of gene expression”. After treatment, the cells were washed with PBS, fixed with 4% paraformaldehyde, permeabilized with Triton X-100, and blocked with a blocking solution. Cells were then incubated with JNK Polyclonal antibody (Cat No. 51153-1-AP), followed by incubation with a FITC-labeled secondary antibody (BA1105) for fluorescence detection. Finally, nuclei were counterstained with DAPI. Fluorescence images were captured under a fluorescence microscope, and fluorescence intensity was quantified using ImageJ software, normalized, and visualized using GraphPad Prism 9.5 software. Western blot analysis RAW264.7 cells were seeded into 6-well plates and incubated for 24 h under standard conditions. Cell grouping and treatment were performed as described in Section “[93]Real-time quantitative PCR (RT-qPCR) analysis of gene expression”. Following treatment, cells were washed with PBS, and total protein was extracted using RIPA lysis buffer supplemented with PMSF and phosphatase inhibitors. The protein lysates were separated via SDS-PAGE, transferred onto a PVDF membrane, and subsequently blocked. The membrane was then incubated with a JNK polyclonal antibody (Cat No. 51153-1-AP), followed by an HRP-conjugated goat anti-rabbit secondary antibody (Cat No. SA00001-2) after PBS washing. Protein bands were detected using an enhanced chemiluminescence (ECL) system, and their intensities were quantified with ImageJ software. The data were then normalized and graphically presented using GraphPad Prism 9.5. Results Identification of differentially expressed genes (DEGs) After data preprocessing, 432 differentially expressed genes (DEGs) were identified in osteoarthritis (OA) tissues (Fig. [94]2A). Among them, 199 genes were downregulated, and 233 genes were upregulated. The top 50 most significantly upregulated and downregulated genes were selected for further analysis (Fig. [95]2B). Fig. 2. [96]Fig. 2 [97]Open in a new tab (A) Volcano plot and (B) Dendrogram for DEGs. Identification of active components and their corresponding targets in SBD By integrating data from the Traditional Chinese Medicine Systems Pharmacology (TCMSP) database and eliminating duplicates, 114 potentially active compounds were identified across different herbal sources: Aconite (2 compounds), Angelica sinensis (2 compounds), Saposhnikovia divaricata (18 compounds), Coix lacryma-jobi (6 compounds), Glycyrrhiza uralensis (87 compounds), Notopterygium incisum (13 compounds), Gelsemium elegans (1 compound), and Gentiana macrophylla (1 compound). Utilizing the UniProt database for subsequent verification and refinement, we ultimately identified targets from Aconite (4 targets), Angelica sinensis (42 targets), Saposhnikovia divaricata (64 targets), Coix lacryma-jobi (27 targets), Glycyrrhiza uralensis (199 targets), Notopterygium incisum (37 targets), Gelsemium elegans (19 targets), and Gentiana macrophylla (27 targets). For additional details, please consult the supplementary materials. Potential OA-related target genes Using databases such as GeneCards, DISGENET, OMIM, and TTD, 2,046 target genes related to OA were identified after removing duplicates (Fig. [98]3A). Fig. 3. [99]Fig. 3 [100]Open in a new tab (A) Venn diagram of “OA” targets; (B) Venn diagram of “SBD-OA” targets; (C) PPI network diagram of “SBD-OA” targets; (D) Degree algorithm for screening results; (E) MCC algorithm for screening results; (F) MNC algorithm for screening results. Construction of PPI network Venn diagram analysis identified 113 overlapping targets between SBD and OA, highlighting potential therapeutic targets (Fig. [101]3B). These targets were imported into the STRING database with a high-confidence threshold of 0.700 for PPI network construction (Fig. [102]3C). Analysis using Cytoscape (v3.10.2) revealed 107 genes and 930 edges in the network. The CytoHubba plugin identified key genes based on Degree, MCC, and MNC metrics (Figs. [103]3D–F). Drug-active ingredient-common target network To clarify the interactions between SBD active compounds and their targets, a drug-active ingredient-target network was established (Fig. [104]4). This network comprised 237 nodes, including 1 formula name, 1 disease name, 8 herbal sources, 114 active compounds, and 113 target genes, interconnected by 948 edges. Fig. 4. [105]Fig. 4 [106]Open in a new tab “SBD-ingredient-target” network diagram. GO and KEGG pathway enrichment analysis To gain a clear understanding of the functions of the intersecting genes, GO (Gene Ontology) functional enrichment and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses were conducted on 113 intersecting targets using R software (v4.4.1) with a significance threshold of P < 0.05. GO analysis revealed that the biological processes (BP) were predominantly associated with antioxidant responses, TNF-mediated signaling pathways, cellular responses to bacterial components and lipopolysaccharides, as well as inflammation-related mechanisms, including cell proliferation and the regulation of apoptotic signaling pathways. The cellular component (CC) was mainly regions such as membrane rafts, membrane microdomains, vesicular lumens, secretory granule lumens, and the outer mitochondrial membrane; the molecular function (MF) was primarily related to ligand activity, DNA-binding transcription factor binding, protein kinase regulatory activity, and receptor binding functions. The top 10 terms for biological processes (BP), cellular components (CC), and molecular functions (MF) were visualized using multiple chart formats with the ggplot2 package (Fig. [107]5A). KEGG pathway analysis identified 175 significantly enriched pathways (P < 0.05), with key pathways including lipid metabolism and atherosclerosis, the RAGE signaling pathway in AGE-related diabetic complications, the IL-17 signaling pathway, and the TNF signaling pathway. The top 30 pathways were visualized using the ggplot2 package (Fig. [108]5B). The roles of the AGE-RAGE signaling pathway and IL-17 signaling pathway are particularly significant in the regulation of inflammation and oxidative stress by SBD in OA (Fig. [109]5C,D). Fig. 5. [110]Fig. 5 [111]Open in a new tab (A) Lollipop Plot for GO enrichment analysis; (B) Lollipop Plot for KEGG enrichment analysis; (C) AGE-RAGE signaling pathway diagram; (D) IL-17 signaling pathway diagram. Target-KEGG network construction To clarify the interconnections between SBD targets and pathways in OA treatment, we established a drug-active compound-common target-KEGG network. This network consisted of 267 nodes, including 4 core targets, 109 additional targets, 30 pathways, and 1,608 edges (Fig. [112]6). Among them, the 4 core targets, which hold significant positions in Degree, MCC, and MNC analyses, were selected for subsequent molecular docking studies. Fig. 6. [113]Fig. 6 [114]Open in a new tab Target-KEGG network diagram illustration. Note: Large circular nodes of various colors represent drugs, small circular nodes represent active ingredients, purple octagonal nodes represent KEGG pathways, blue circular nodes represent core targets, and green rectangular nodes represent intersecting targets. Molecular docking The core active compounds (Shinflavanone, Gancaonin L, Xambioona, Phaseol, Gancaonin O, and Licoisoflavanone) were used as ligands, while key therapeutic proteins (IL6, TNF, STAT3, IL1B, and JNK) were chosen based on their hydrogen bonding interactions and binding energies. The results were displayed as binding energy heatmaps using ggplot2 (Fig. [115]7A). Notably, JNK acts as an upstream regulator of IL6, TNF, and IL1B. Fig. 7. [116]Fig. 7 [117]Open in a new tab (A) Molecular Docking Heatmap: Docking Results of Core Targets and Binding Energy (kcal·mol⁻¹); (B) 3D and 2D images of Xambioona_IL6; (C) 3D and 2D images of Gancaonin L_TNF; (D) 3D and 2D images of Shinflavanone_STAT3; (E) 3D and 2D images of Shinflavanone_JNK; (F) 3D and 2D images of Shinflavanone_IL1B. Shinflavanone showed the lowest binding energy with IL1B, JNK, and STAT3, mainly interacting through hydrogen bonds (VAL A:41; ASN A:114) and Van der Waals forces (CYS A:251; ASP A:334) in the binding pocket. Xambioona had the strongest affinity for IL6, stabilizing the protein-ligand complex via hydrogen bonding (VAL A:97) and Van der Waals interactions (LEU A:93; ASN A:145). Gancaonin L exhibited the lowest binding energy with TNF, forming a stable complex primarily through Van der Waals forces (SER B:99) (see Supplementary Material S4 for details). The best-performing candidates for each complex were visualized using PyMOL and Discovery Studio (Fig. [118]7B–E). Additionally, three independent docking repetitions produced consistent and stable conformations, validating the reliability of molecular docking performed with AutoDockTools. MD simulation Based on the molecular docking results, 200-ns MD simulations were performed for the following protein-ligand complexes: Xambioona_IL6, Gancaonin L_TNF, Shinflavanone_STAT3, Shinflavanone_JNK, and Shinflavanone_IL1B. The stability of the trajectory files was assessed using multiple metrics, including root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), radius of gyration (Rg), solvent-accessible surface area (SASA), hydrogen bond dynamics (HB), and principal component analysis-free energy landscape (PCA). As shown in Fig. [119]8A, the RMSD variations for all ligands remained below 0.1 nm, indicating high stability of the ligands throughout the simulation. The RMSD profile of the Shinflavanone_IL1B complex exhibited noticeable fluctuations, suggesting structural and conformational changes due to protein folding. In contrast, the other complexes displayed relatively stable RMSD profiles throughout the simulation. Fig. 8. [120]Fig. 8 [121]Open in a new tab (A) RMSD; (B) RMSF; (C) Rg; (D) SASA; (E) HB; (F) PCA plots during molecular dynamics simulation. Note: From left to right are Xambioona_IL6, Gancaonin L_TNF, Shinflavanone_STAT3, Shinflavanone_IL1B, and Shinflavanone_JNK, PC1 representing RMSD and PC2 representing RG. As shown in Fig. [122]8B–D, the Shinflavanone_IL1B complex demonstrated increased flexibility in two regions (31–38 and 48–56), while the Shinflavanone_STAT3 complex exhibited similar flexibility in the 170–220 region. The Xambioona_IL6 complex displayed flexibility in three regions (42–70, 71–81, and 130–170). The Rg and SASA analyses revealed low fluctuation levels, indicating greater overall structural stability of the complexes. HB analysis revealed fewer than four hydrogen bonds between the molecules, suggesting that the interactions were primarily driven by van der Waals forces, which was further corroborated by the subsequent MM-PBSA analysis (Fig. [123]8E). PCA plots during molecular dynamics simulation revealed the dominant conformational changes and dynamic behaviors of the protein throughout the simulation (Fig. [124]8F). The MM-PBSA results, used to quantify the binding free energy of the complexes (see Supplementary Material S5 for details), confirmed that all ligand-protein complexes exhibited favorable binding affinity. Effects of SBD on cell viability As shown in Fig. [125]9A, CCK-8 assay results indicate that SBD treatment does not exhibit significant cytotoxicity at concentrations up to 100 µg/mL, However, at 200 µg/mL, a noticeable decrease in cell viability was observed. In the concentration range of 12.5–100 µg/mL, cell viability remained stable with a slight proliferative effect. Therefore, 100 µg/mL was selected as the non-cytotoxic concentration for subsequent experiments, ensuring both safety and therapeutic efficacy. Fig. 9. [126]Fig. 9 [127]Open in a new tab (A) Cell viability test of SBD treated with different concentrations; (B) Relative expression of IL1-β; (C) Relative expression of IL-6; (D) Relative expression of STAT3; (E) Relative expression of TNF; (F) Relative expression of JNK. Effects of SBD on inflammatory mediators and ROS QRT-PCR analysis revealed that LPS stimulation significantly induced the expression of inflammatory cytokines, while SBD treatment markedly downregulated the mRNA levels of JNK, STAT3, IL1B, IL6, and TNF (Fig. [128]9B–F). As shown in Fig. [129]10A and B, LPS stimulation led to a significant increase in ROS levels, with a noticeable enhancement in fluorescence intensity. However, SBD treatment significantly reduced fluorescence intensity and ROS levels, indicating its potential anti-inflammatory and antioxidant effects. Fig. 10. [130]Fig. 10 [131]Open in a new tab (A) ROS clearance was measured using a DCFH-DA probe; (B) Relative fluorescence intensity. Effects of SBD on JNK axis As illustrated in Fig. [132]11A–D, immunofluorescence, and Western blot analysis demonstrated that LPS stimulation significantly increased JNK levels, leading to elevated fluorescence intensity and gene expression. After SBD treatment, both fluorescence intensity and JNK were significantly reduced, suggesting that SBD effectively inhibits JNK activation. Fig. 11. [133]Fig. 11 [134]Open in a new tab (A, B) IF staining and fluorescence intensity analysis of JNK; (C, D) JNK protein expression was determined by western blot test. Discussion Osteoarthritis (OA) is a common orthopedic disease that imposes a substantial burden on society. Its prevalence increases with age, making OA an emerging global public health issue. Studies indicate that Traditional Chinese Medicine (TCM) shows potential in treating various chronic diseases and tumors^[135]28–[136]31. We employed a network pharmacology-based computational approach, combined with molecular docking, molecular dynamics simulation, and experimental validation, to explore the potential mechanisms of SBD in the treatment of OA. Through an analysis of multiple databases, including GEO, we identified 114 active chemical components and 113 intersecting targets, screening six active compounds with excellent binding affinities to core targets. Shinflavanone is an effective bone resorption inhibitor, suggesting its potential to reverse OA pathology^[137]32. Gancaonin L and Gancaonin O, flavonoid compounds extracted from licorice, possess immunomodulatory and other biological activities, making them promising candidates for OA treatment. Xambioona, a natural plant-derived compound, and its therapeutic mechanism in OA remain underexplored. Phaseol, a bioactive compound derived from legumes, holds antioxidant, anti-inflammatory, and antidiabetic potential. Licorisoflavanone, an isoflavonoid compound, may have a bone-repairing effect. By utilizing Cytoscape software and PPI network analysis, we identified four core targets from the compound-drug-component-disease-target network. The inflammatory cytokine Interleukin-1 beta(IL1β) induces an inflammatory phenotype characterized by increased RANKL and matrix metalloproteinases (MMPs)^[138]33. Additionally, IL-1β enhances the production of inflammatory cytokines like Tumor Necrosis Factor (TNF) and Interleukin-6(IL-6), exacerbating OA-related inflammation^[139]34. IL-1B inhibitors significantly reduce OA progression by preserving cartilage integrity and minimizing lesion size^[140]35. Elevated IL-6 levels are associated with increased pain, stiffness, and radiographic severity^[141]36. Excessive IL-6 production negatively affects the synthesis of cartilage matrix proteins, accelerating cartilage degradation^[142]37. Targeting IL-6 can alleviate its harmful effects on cartilage metabolism^[143]38. TNF influences various cellular processes that lead to cartilage degeneration and joint pathology^[144]39. TNF plays a critical role in OA by activating signaling pathways such as NF-kappa B, MAPK, and PI3K/Akt, promoting inflammatory cell recruitment and chondrocyte degradation^[145]40. TNF inhibitors have been shown to have the potential for inflammation relief and cartilage protection^[146]41. Activation of STAT3 is associated with the expression of IL-1β, IL-6, and TNF^[147]42,[148]43. The signal transducer and activator of transcription 3(STAT3) is involved in regulating chondrocyte apoptosis^[149]44,[150]45. Therefore, simultaneously targeting these targets offers significant potential for the treatment of OA. In the Gene Ontology (GO) analysis, Biological Processes (BP) were significantly associated with inflammation and cellular stress, including antioxidant responses, TNF-mediated signaling, cellular responses to bacterial components and lipopolysaccharides, as well as the modulation of cell proliferation and apoptotic pathways. Cellular Components (CC) were significantly enriched in regions such as lipid rafts, membrane microdomains, vesicle lumen, secretory granule lumen, and mitochondrial outer membrane, indicating that these cellular components play crucial roles in inflammatory signaling and functional regulation. Molecular Functions (MF) involve ligand activity, DNA-binding transcription factor binding, protein kinase regulatory activity, and receptor binding, which are likely associated with signal transduction and transcriptional regulation. GO analysis revealed the multidimensional roles of SBD in cellular functions and processes. In the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, SBD may be multi-pathway synergistic in OA treatment. Key inflammation- and immune-related pathways, including the AGE-RAGE pathway, TNF signaling pathway, IL-17 signaling pathway, and Toll-like receptor signaling pathway, indicate that SBD may exert therapeutic effects on OA by regulating these signaling cascades. Reactive oxygen species (ROS) is a key driver of OA progression, playing a crucial role in inflammation and cartilage degradation in the pathology of OA^[151]46. The AGE-RAGE axis enhances ROS generation while simultaneously inhibiting antioxidant enzyme activity, thereby weakening cellular defenses against oxidative stress^[152]47. Increased oxidative stress leads to DNA damage and mitochondrial dysfunction in chondrocytes, exacerbating cartilage degeneration and inflammation^[153]48. The JNK axis is a key component in the AGE-RAGE signaling pathway. JNK is considered a precursor of NF-κB. It not only exacerbates the inflammatory cascade but also promotes apoptosis by regulating the STAT3^[154]49. Continued activation of the JNK signaling pathway is implicated in cartilage degeneration and chronic synovial inflammation, making it a significant therapeutic target for OA^[155]50,[156]51. Inhibiting JNK expression helps alleviate the persistent cartilage degeneration and chronic inflammation in OA^[157]52. Inhibiting STAT3 activation is crucial in the progression of OA and cartilage damage^[158]53. IL-17 signaling pathway is closely associated with OA progression. The crosstalk between the IL-17 pathway and the AGE-RAGE pathway may involve JNK, which activates NF-κB, triggering inflammation and cartilage degradation, thereby further activating the IL-17 pathway^[159]54. Another study suggested the IL-17 axis as a novel therapeutic strategy for OA^[160]55. In OA animal models, SPHK2 knockdown inhibited the IL-17 pathway and alleviated cartilage damage and synovial inflammation^[161]56. In rat chondrocytes, TNF has been shown to increase senescence markers and glycosaminoglycan (GAG) degradation^[162]57. Clinical studies have demonstrated the potential of TNF inhibitors in treating arthritis^[163]58. Research indicates that IL-17 A and TNF synergistically induce pro-inflammatory mediators in synovial fibroblasts, accelerating OA progression^[164]59. Overall, enrichment analysis suggests that SBD may effectively treat OA through multi-gene, multi-pathway mechanisms. Using three algorithms (Degree Centrality, Maximum Clique Centrality, and Mean Nearest Neighbor Centrality), the top four core (IL1B, STAT3, IL6, TNF) genes were identified. Furthermore, within the AGE-RAGE pathway, JNK is a common upstream target for the four core targets and ranks highly. Therefore, this study specifically explores the relationship between the JNK axis in the AGE-RAGE pathway and OA. Molecular docking simulations suggest that the six key compounds (Shinflavanone, Gancaonin L, Xambioona, Phaseol, Gancaonin O, and Licoisoflavanone) exhibit significant binding activity with the five target proteins (IL1B, STAT3, IL6, TNF, and JNK) under investigation. Shinflavanone exhibited the highest docking score with JNK, demonstrating stable hydrogen bonds and strong van der Waals interactions. Shinflavanone_IL1B, Xambioona_IL6, Shinflavanone_STAT3, and Gancaonin L_TNF also showed favorable interactions, primarily characterized by strong van der Waals forces and hydrophobic interactions. Molecular dynamics simulations further validated this finding by demonstrating the dynamic stability of these interactions in solution. Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuation (RMSF) analyses indicated that all ligand-protein complexes exhibited relatively stable conformations, with Shinflavanone_JNK showing particularly stable behavior. The radius of gyration (Rg) and solvent-accessible surface area (SASA) analysis demonstrated the stability of the binding interface throughout the simulation. The consistent number of hydrogen bonds (HB) maintained throughout the simulation also supports the stability of the complex over time. The MM-PBSA method can quantitatively analyze the binding energy of ligand-protein complexes, offering unique advantages. MM-PBSA analysis further validated the dynamic stability of these complexes. In conclusion, all complexes exhibited significant stability, providing strong theoretical support and data evidence for the development of new therapeutic approaches. To validate the effects of SBD, we conducted in vitro experimental validation. The results indicate that SBD effectively scavenges ROS and inhibits the expression of IL-1β, IL-6, TNF-α, STAT3, and JNK. This significantly highlights the potential of SBD in alleviating oxidative stress, exerting anti-inflammatory effects, and inhibiting cartilage degradation. In addition, both immunofluorescence (IF) and Western blot (WB) experiments have demonstrated the importance of the JNK axis in the AGE-RAGE signaling pathway. Compared to conventional anti-inflammatory therapies such as NSAIDs and biologics, SBD may exert a broader regulatory effect by simultaneously targeting multiple inflammatory mediators. In conclusion, Experimental validation further supports the hypothesis that SBD alleviates OA by regulating the JNK axis within the AGE-RAGE signaling pathway, thereby disrupting the pathological loop in which JNK activation induces ROS production, which in turn amplifies AGE-RAGE signaling, exacerbating inflammation and tissue damage. By breaking this vicious cycle, SBD may mitigate inflammatory responses and preserve joint integrity. Future studies should explore its potential in combination with standard therapies to enhance efficacy and reduce adverse effects, offering a novel multi-targeted strategy for OA management. Conclusion This study comprehensively explored the potential of SBD in the treatment of OA through network pharmacology, molecular docking, MD simulations, and experimental validation. The results indicate that SBD may exert therapeutic effects on key inflammatory targets in OA through its key components, including Shinflavanone, Gancaonin L, Xambioona, Phaseol, Gancaonin O, and Licoisoflavanone. Moreover, we validated that SBD alleviates oxidative stress, exhibits anti-inflammatory effects, and protects cartilage through the JNK pathway. Overall, our findings confirm the reliability of computational techniques and highlight the potential of SBD as a therapeutic agent for OA (Fig. [165]12). Fig. 12. [166]Fig. 12 [167]Open in a new tab Validation mechanism of SBD in the treatment of OA. Electronic supplementary material Below is the link to the electronic supplementary material. [168]Supplementary Material 1^ (17.7KB, xlsx) [169]Supplementary Material 2^ (42.7KB, xlsx) [170]Supplementary Material 3^ (14.3KB, xlsx) [171]Supplementary Material 4^ (45.9KB, xlsx) [172]Supplementary Material 5^ (11.2KB, xlsx) [173]Supplementary Material 6^ (130.5KB, pdf) Acknowledgements