Unraveling the Therapeutic Mechanism of Saussurea involucrata against Rheumatoid Arthritis: A Network Pharmacology and Molecular Modeling-Based Investigation

Rheumatoid arthritis (RA) is a chronic autoimmune disease with a global prevalence of approximately 0.46%, causing significant impairments in patients’ quality of life and an economic burden. Saussurea involucrata (SI) has long been used in traditional medicine to treat RA, but its underlying mechanism remains unclear. This study utilized network pharmacology and molecular docking to explore the potential pharmacological effects of bioactive compounds in SI on RA. A total of 27 active compounds were identified, along with 665 corresponding targets. Additionally, 593 disease-related targets were obtained from multiple databases, with 119 common targets shared with SI. The high-ranking targets mainly belong to the MAPK family and NF-κB pathway, including MAPK14, MAPK1, RELA, TNF, and MAPK8, all of which are associated with inflammation and joint destruction in RA. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed significant pathways related to IL-17 signaling, Th17 cell differentiation, and osteoclast differentiation. Molecular docking and dynamic simulations demonstrated strong interactions between several flavonoids and RA-related targets. Xuelianlactone, Involucratin, and Flazin exhibit outstanding binding efficacy with targets such as MAPK1, MAPK8, and TNF. These findings provide valuable insights into the therapeutic potential of SI for RA and offer directions for further drug development.


Introduction
Rheumatoid arthritis (RA), a chronic autoimmune disease, exhibited a global prevalence of 0.46% between 1980 and 2018, with significant regional and country-specific variations [1].In mainland China, where the prevalence of RA was 0.42%, similar to the world average [2], the number of RA patients reached 5 million, with an average onset age of 45 years.Among all patients, 82% experienced moderate to severe disease, and patients often have complications [3], highlighting an ongoing and concerning situation.With a female-to-male ratio of about 4:1 in China [4], RA is more prevalent among women, and this overrepresentation is likely influenced by genetic (X-linked) factors and hormonal aspects, even though the exact reasons remain unclear [5].The occurrence of RA is influenced by a combination of genetic and environmental factors, including smoking, infection, alcohol intake, gender, and age [6][7][8].RA patients may suffer from various symptoms, including joint destruction, deformity, disability, and even death.In some cases, organs or systems such as the heart, kidneys, lungs, eyes, skin, digestive system, and nervous system may also be affected, leading to the development of various syndromes [9].RA causes physical and emotional pain to patients and shortens their life expectancy.Additionally, the teratogenic and disabling effects of RA contribute to a loss of labor capacity within the population, resulting in substantial economic losses [10].
There is no cure for RA at present, but there are some pharmacological and nonpharmacological treatments that can relieve symptoms, delay disease progression, and improve quality of life [11].The main purpose of drug treatment is to reduce inflammation, suppress excessive immune system response, and reduce joint damage.The currently used drugs include non-steroidal anti-inflammatory drugs (NSAIDs), disease-modifying anti-rheumatic drugs (DMARDs), glucocorticoids (GC), TNF inhibitors, and polymer colloids [12,13].These drugs may cause gastrointestinal side effects such as nausea, vomiting, or abdominal pain, and adverse reactions like stomatitis or mouth sores.Moreover, they can lead to liver toxicity (hepatic dysfunction), hematological disorders, bone density reduction, insomnia, and depression.In more severe cases, patients may experience infections, cardiovascular events, cancer, or even death [14,15].Therefore, the development of alternative drugs for RA treatment is necessary.Traditional Chinese Medicine and other ethnic medicines have a long history of clinical application and are considered valuable sources for natural drug development.These medicines offer higher safety profiles and have recently garnered significant attention in the field of RA drug research.
Saussurea, a medicinal plant, was first documented in the ancient Tibetan book Smandpyad Zla-ba'i-rgyal-po [16] in China.Among its species, Saussurea involucrata (SI), grown in Xinjiang, China, is particularly noteworthy and widely used in Uyghur Medicine, known as "Tage leylishi" [17].According to the first part of the Chinese Pharmacopoeia (2020 edition), this herb is utilized in both Traditional Chinese Medicine and Uyghur Medicine to address conditions such as RA, joint pain, irregular menstruation, and excessive leucorrhea [18].According to the current literature, SI is found to contain a variety of compounds, including phenolic acids, flavonoids, lignans, phenylpropanoids, steroids, coumarins, sesquiterpenes, ceramides, glycosides, and polysaccharides [19,20].Ongoing studies primarily concentrate on exploring the potential therapeutic effects of its flavonoids [21,22].
In China, Saussurea has found extensive clinical use in the treatment of RA, available in various dosage forms such as injections, capsules, and more [23].However, limited pharmacological studies have been conducted, with the majority of existing research focused on clinical efficacy verification.This study aims to employ network pharmacology methods to analyze the mechanism of SI in treating RA and verify the binding between its bioactive components and target proteins through molecular modeling.The findings will serve as a basis for subsequent drug development endeavors.

Bioactive Components of SI and Prediction of Target Genes
To ensure a comprehensive collection of the bioactive components of SI, this study conducted searches across multiple databases.The databases used included the Traditional Chinese Medicine Systems Pharmacology Database (TCMSP, tcmspw.com/tcmsp.php,accessed on 27 February 2023) [24], the Chinese Medicine and Bioactive components Database of the Shanghai Institute of Organic Chemistry (organchem.csdb.cn/scdb/default.htm, accessed on 5 April 2023) [25], and SymMap (www.symmap.org,accessed on 12 April 2023) [26].The keyword "Saussurea involucrata" was employed for these searches, yielding relevant results.Furthermore, additional active ingredients were supplemented by searching the literature in the China National Knowledge Infrastructure(CNKI) [19,27,28].
The obtained components were screened based on parameters such as oral bioavailability and drug-likeness to ensure their potential for metabolism within the human body.For the query results from the SymMap database, components with an oral bioavailability (OB) threshold of ≥30% were selected [26].These selected components, along with the components obtained from literature searches, were used to query the PubChem database (pubchem.ncbi.nlm.nih.gov,accessed on 12 April 2023) using the pubchempy Python toolkit, acquiring the Simplified Molecular Input Line Entry System (SMILES) notations.These notations were then imported into the SwissADME database (www.swissadme.ch,accessed on 16 April 2023) and subjected to screening based on criteria such as a gastrointestinal absorption (GI) rating of "High" and drug-likeness having two or more items marked as "Yes" [29].Subsequently, the screened components were combined with the components directly obtained from the TCMSP database and underwent a further screening within TCMSP according to the standards of an oral bioavailability (OB) threshold of ≥30%, a drug-likeness (DL) threshold of ≥0.18, and a Caco-2 permeability threshold of ≥−0.4 [30,31].Finally, all the screened components were compiled, and any duplicate entries were removed.
For the components queried from the TCMSP, their targets were also obtained from this database, and their gene names were queried from the Uniprot database (www.uniprot.org,accessed on 27 February 2023).For other components, their SMILES were input into SwissTargetPrediction (swisstargetprediction.ch, accessed on 16 April 2023) with the species attribute set to "Homo sapiens" for target prediction, and targets with a probability greater than 0 were collected [32].Finally, all component targets were summarized and deduplicated.

Target Genes of RA
Disease targets were collected through searches across multiple databases.The DisGeNet database (www.disgenet.org,accessed on 6 April 2023) was queried using the keyword "rheumatoid arthritis", and disease-related target genes were obtained by applying a screening criterion of Score ≥ 0.3 [33].Similarly, the GeneCards database (www.genecards.org,accessed on 27 February 2023) was queried, and data with a Relevance score > 10 were selected [34].Furthermore, the OMIM database (Online Mendelian Inheritance in Man, www.omim.org/search/advanced/geneMap,accessed on 13 April 2023) and the TTD database (Therapeutic Target Database, db.idrblab.net/ttd,accessed on 14 April 2023) were employed to query targets associated with RA [35].Subsequently, all the retrieved results were compiled, summarized, and deduplicated to create the disease target library.

Protein-Protein Interaction (PPI) Network Analysis
Components exert their effects on the disease by targeting specific genes.The target genes for both the components and the disease were inputted into the Bioinformatics website (www.bioinformatics.com.cn,accessed on 17 May 2023) to generate a Venn diagram [36], obtaining the intersection of targets.
In order to elucidate the relationships among target proteins and further refine key targets, the intersecting target genes were inputted into the STRING database (string-db.org,accessed on 17 May 2023) with a confidence level set at >0.9, and the species selected as H. sapiens (human) [37].Based on this, the PPI network was constructed.The results were then imported into Cytoscape 3.9.1 software, of which, the CytoNCA plugin was utilized to calculate various centrality metrics for each node, including betweenness centrality (BC), closeness centrality (CC), degree centrality (DC), eigenvector centrality (EC), local average connectivity-based method centrality (LAC), and network centrality (NC) [38].The median values of each centrality metric were employed as the screening criteria, and target proteins surpassing or equaling these median values across all indicators were identified as key targets for this study.

Gene Ontology and Pathway Analysis
The key targets were uploaded to the Bioinformatics website for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis.The results were subjected to screening based on the criterion of p.adjust < 0.05 [39], and subsequently plotted using the Count, p.adjust, and GeneRatio values.As for the GO analysis results, they were sorted and visualized separately for the cellular component (CC), molecular function (MF), and biological process (BP) categories.

Network Visualization
Information such as component-target, disease-target, drug-component, pathwaytarget, function-target, etc., were summarized and imported into Cytoscape for network visualization, to construct their relationship network.Then, network analysis was performed using the network analysis function [37], and according to the obtained degree, 10 components and 10 targets with a higher ranking were selected as core components and core targets.

Molecular Docking
The core targets were searched in the RCSB PDB database (RCSB Protein Data Bank, www.rcsb.org,accessed on 18 May 2023) to retrieve structures with species as H. sapiens and a resolution of less than or equal to 2.5 Å.The structures were selected based on parameters such as Gene Name, Resolution (Å), ligand, R-Value, pH, Experimental Method, Publication Year, etc. [40].Priority was given to structures containing similar ligands, with particular attention paid to structures with higher resolution, pH conditions closely resembling the human body environment, and more recent publication years.For each target, one structure was chosen and downloaded.
The downloaded protein structure was refined using Swiss-PdbViewer 4.1.0[41], remove water and ligands with Pymol, and then imported into AutoDockTools for receptor preparation [42], which involved operations such as hydrogen addition and charge assignment.The structures of the core components were obtained from the PubChem database, and initial conformations for each compound were established by energy minimization using the MM2 force field in Chem3D 14.0 software [42].After converting these files into appropriate formats, they were imported into AutoDockTools for ligand preparation.
To determine the position and size of the Grid Box, priority was given to the original ligand positions within the structures, taking into account the target position descriptions from the Uniprot database.In the case of structures without their own ligands, the Deepsite tool (www.playmolecule.com/deepsite,accessed on 28 May 2023) [43] on the PlayMolecule website was utilized for docking pocket prediction.The position and size of the grid box were determined based on the pocket scores and the 3D structure observations.Finally, AutoDock Vina 1.1.2[42] was employed to perform molecular docking on a total of 100 receptor-ligand combinations.Each combination yielded 10 conformations along with binding energy results.

Molecular Dynamics (MD) Simulation
The well-combined target protein-ligand complex was selected from the molecular docking results, and MD simulations were performed using Gromacs 2023.2.For protein molecules, Charmm36 force field and TIP3P water model were chosen for processing [44].Ligand conformations with good docking affinity were imported into Avogadro 1.2.0 to add hydrogens and format conversion.Then, the Sobtop tool (Tian Lu, Sobtop, Version 1.0 (dev3.1),http://sobereva.com/soft/Sobtop,accessed on 4 September 2023) was used to generate topology files with the GAFF force field.The receptor-ligand complex was subjected to solvation, ion balance and energy minimization.In order to simulate the human body environment, the ion concentration was set at 0.145 M when adding ions, and the MD simulation of the system was carried out at a temperature of 310 K and a pressure of 1 bar for 50 ns.The cut-off value of the non-bonded interaction was set to 1.2 nm, and the long-distance electrostatic interaction was calculated [45].Energy calculations were carried out on complexes that exhibit favorable results in MD simulations using the gmx_mmpbsa.py[46] and gmx_mmpbsa script of gmxtool (Jicun Li, gmx_mmpbsa.bsh,DOI 10.5281/zenodo.6408973,accessed 15 September 2023).

Bioactive Components and Key Targets
After the above component-mining, a total of 27 components of SI were obtained (Table 1), corresponding to 665 non-redundant targets.For disease targets, a total of 593 disease-related targets were obtained from four databases.By analyzing the Venn diagram, 119 intersection targets were identified, representing the overlapping targets between the component targets and the disease targets (Figure 1).

PPI Network Analysis
Figure 2 presents the results of PPI analysis, revealing the interactions between target proteins.The significance of a target protein is determined by the number of interactions it has with other proteins.Based on the network analysis conducted in Cytoscape, a total of 33 key targets were identified (Table 2).These key targets play crucial roles in the network and are considered important in this study.

KEGG Pathway Enrichment Analysis
The key targets were input into the Bioinformatics website and, after screening, a total of 131 pathways with p.adjust < 0.05 were obtained.Among these pathways, the top 20 pathways with the smallest p.adjust values were selected and presented in Figure 3.The color in the plot transitions from green to red, indicating an increasing level of significance.Furthermore, the size of each circle corresponds to the number of genes enriched within the respective pathway.Larger circles signify a higher gene enrichment count.

PPI Network Analysis
Figure 2 presents the results of PPI analysis, revealing the interactions between targe proteins.The significance of a target protein is determined by the number of interactions it has with other proteins.Based on the network analysis conducted in Cytoscape, a tota of 33 key targets were identified (Table 2).These key targets play crucial roles in the net work and are considered important in this study.

PPI Network Analysis
Figure 2 presents the results of PPI analysis, revealing the interactions between target proteins.The significance of a target protein is determined by the number of interactions it has with other proteins.Based on the network analysis conducted in Cytoscape, a total of 33 key targets were identified (Table 2).These key targets play crucial roles in the network and are considered important in this study.

GO Enrichment Analysis
Based on the output results from the Bioinformatics website, the findings were categorized into three groups, as shown in Figure 4.In both the MF and BP analyses, multiple results met the criteria of p.adjust < 0.05.The top 10 results with the smallest p.adjust values were selected for plotting.As for the CC analysis, only four results met the requirements, and all of them were included in the plot.The plot exhibits a color gradient shifting from blue to red, indicating an increasing level of significance.Additionally, the length of each bar corresponds to the number of genes enriched within the specific category, with longer bars signifying a higher gene enrichment count.

Network Constructions and Analysis
Information such as bioactive components, key targets, related pathways, and the functions obtained above was formatted in a network diagram using Cytoscape software (Figure 5).In the diagram, green circles represent the bioactive components present in SI, which intersect with RA-related targets.The intersection results are depicted by blue squares.Among them, the inner two circles represent the 33 key targets identified through network analysis, while the innermost circle represents the 10 core targets selected for molecular docking.Yellow hexagons represent the top 20 pathways obtained from the key target analysis.At the bottom of the diagram, three distinct groups illustrate the results of the GO analysis.Purple, orange, and red represent the BP, CC, and MF analysis results, respectively.This comprehensive analysis unveils the intrinsic correlation between the bioactive components of SI and their potential in treating RA.

GO Enrichment Analysis
Based on the output results from the Bioinformatics website, the findings were categorized into three groups, as shown in Figure 4.In both the MF and BP analyses, multiple results met the criteria of p.adjust < 0.05.The top 10 results with the smallest p.adjust values were selected for plotting.As for the CC analysis, only four results met the requirements, and all of them were included in the plot.The plot exhibits a color gradient shifting from blue to red, indicating an increasing level of significance.Additionally, the length of each bar corresponds to the number of genes enriched within the specific category, with longer bars signifying a higher gene enrichment count.

GO Enrichment Analysis
Based on the output results from the Bioinformatics website, the findings were categorized into three groups, as shown in Figure 4.In both the MF and BP analyses, multiple results met the criteria of p.adjust < 0.05.The top 10 results with the smallest p.adjust values were selected for plotting.As for the CC analysis, only four results met the requirements, and all of them were included in the plot.The plot exhibits a color gradient shifting from blue to red, indicating an increasing level of significance.Additionally, the length of each bar corresponds to the number of genes enriched within the specific category, with longer bars signifying a higher gene enrichment count.

Molecular Docking
According to the network analysis results of Cytoscape, 10 compounds with est degree values, Quercetin, Luteolin, Acacetin, Xuelianlactone, Moslosooflav pidulin, Mosloflavone, Involucratin, 5,7-Dihydroxy-6,8-Dimethoxyflavone a were selected as molecular docking ligands (core components).Ten targets with est degree values, MAPK14, MAPK1, RELA, TNF, MAPK8, IL6, IL1B, CHUK, IK NFKBIA, were selected as the core targets for molecular docking.Nine protein were selected as the receptors, with CHUK and IKBKB targets located in diff ments of the same protein structure (Table 3).The results of the molecular doc visualized in the form of a heat map (Figure 6), where the color gradient shifting to blue indicates a decrease in affinity.Lower binding energy signifies a bette

Molecular Docking
According to the network analysis results of Cytoscape, 10 compounds with the highest degree values, Quercetin, Luteolin, Acacetin, Xuelianlactone, Moslosooflavone, Hispidulin, Mosloflavone, Involucratin, 5,7-Dihydroxy-6,8-Dimethoxyflavone and Flazin were selected as molecular docking ligands (core components).Ten targets with the highest degree values, MAPK14, MAPK1, RELA, TNF, MAPK8, IL6, IL1B, CHUK, IKBKB and NFKBIA, were selected as the core targets for molecular docking.Nine protein structures were selected as the receptors, with CHUK and IKBKB targets located in different segments of the same protein structure (Table 3).The results of the molecular docking were visualized in the form of a heat map (Figure 6), where the color gradient shifting from red to blue indicates a decrease in affinity.Lower binding energy signifies a better binding activity between the ligand and receptor protein.In general, affinity values > −5 kcal/mol indicate no predicted binding, affinity < −5 kcal/mol suggests moderate predicted binding, and affinity < −7 kcal/mol suggests strong predicted binding [47].In this section, all the selected ligands and receptors demonstrated excellent binding activity, as the lowest binding energy with the best ligand for each receptor was less than −7 kcal/mol.Notably, Xuelianlactone demonstrates the lowest affinity at −10.01 kcal/mol.In addition, the docking results of MAPK1, MAPK8, and TNF showed significantly lower binding affinities compared to the other targets (p < 0.05).However, the binding effect of IL6, ILB, CHUK, and IKBKB to each ligand is not outstanding.selected ligands and receptors demonstrated excellent binding activity, as the lowest binding energy with the best ligand for each receptor was less than −7 kcal/mol.Notably, Xuelianlactone demonstrates the lowest affinity at −10.01 kcal/mol.In addition, the docking results of MAPK1, MAPK8, and TNF showed significantly lower binding affinities compared to the other targets (p < 0.05).However, the binding effect of IL6, ILB, CHUK, and IKBKB to each ligand is not outstanding.For each target protein, the ligand with the lowest affinity was selected from the molecular docking results.Visualization was achieved using Biovia Discovery Studio 2021 software, and the docking results are presented in Table 4 and Figure 7.For each target protein, the ligand with the lowest affinity was selected from the molecular docking results.Visualization was achieved using Biovia Discovery Studio 2021 software, and the docking results are presented in Table 4 and Figure 7.

Molecular Dynamics Simulation
Based on the degree values obtained in network analysis and molecular docking results, the following study primarily focuses on the MAPK1 and TNF.For both targets, the top four compounds with the best docking results, as well as the original ligands (38Z and VGY) of the target proteins, were selected for MD simulations.In addition, for MAPK8, two compounds and the original ligand (38Z) were chosen for simulations.Considering the certainty of the target position, the simulation was conducted with NFKIBA as the representative target for the NF-κB pathway, focusing on the top two best-performing compounds based on the docking results.Table 5 and Figure 8 present the Root Mean Square Deviation (RMSD) values for each of these simulation combinations.

Molecular Dynamics Simulation
Based on the degree values obtained in network analysis and molecular docking results, the following study primarily focuses on the MAPK1 and TNF.For both targets, the top four compounds with the best docking results, as well as the original ligands (38Z and VGY) of the target proteins, were selected for MD simulations.In addition, for MAPK8, two compounds and the original ligand (38Z) were chosen for simulations.Considering the certainty of the target position, the simulation was conducted with NFKIBA as the In MD simulations, RMSD measures the average deviation or fluctuations in atomic positions in a molecular system from their initial positions, providing insight into the stability of a target protein-ligand complex.Generally, it is considered that once the RMSD value stabilizes and the fluctuation range is less than 0.2 nm, the system can be deemed stable [48].It can be observed that, for MAPK1, Involucratin and Luteolin reach equilibrium after about 13 ns, with smaller fluctuations than the protein's original ligand 38Z.The stabilized RMSD value for Involucratin is less than 0.4 nm, and both compounds exhibit fluctuations below 0.2 nm after stabilization.For TNF, the RMSD values for all analyzed compounds remain below 0.4 nm, with fluctuations below 0.2 nm for all except Mosloflavone.Among them, Xuelianlactone and Acacetin demonstrate the highest stability.Regarding MAPK8, the original ligand 38Z exhibits a low RMSD value, but shows significant fluctuations towards the end of the simulation.On the other hand, Flazin maintains an RMSD of 0.39, with fluctuations also within the 0.2 nm range.As for NFKIBA, both selected ligand Luteolin and Acacetin display higher overall fluctuations.In the last 5 ns, their RMSDs tend to stabilize, with fluctuations smaller than 0.2 nm.However, compared to other complexes, their stability remains poor, and thus, they cannot be considered to have achieved stable binding.Based on this, further analysis will be conducted of the following compounds: Involucratin and Luteolin for MAPK1, Xuelianlactone and Acacetin for TNF, and Flazin for MAPK8.representative target for the NF-κB pathway, focusing on the top two best-performing compounds based on the docking results.Table 5 and Figure 8 present the Root Mean Square Deviation (RMSD) values for each of these simulation combinations.In MD simulations, RMSD measures the average deviation or fluctuations in atomic positions in a molecular system from their initial positions, providing insight into the stability of a target protein-ligand complex.Generally, it is considered that once the RMSD value stabilizes and the fluctuation range is less than 0.2 nm, the system can be deemed Solvent-Accessible Surface Area (SASA) quantifies the surface area of a molecule that is accessible to solvent molecules and can indicate changes in a molecule's conformation or exposure during simulation, aiding in the assessment of binding interactions in a target protein-ligand complex.The average SASA values for complexes XT (black), AT (red), IM (blue), LM (green), and FM (purple) during MD simulation are as follows: 202.49± 3.72, 201.65 ± 3.17, 179.03 ± 3.25, 175.97 ± 2.59, and 179.9 ± 2.9 nm 2 .Throughout the entire simulation process, the SASA values of these compounds show small fluctuations, indicating their stable interaction with the protein complexes (Figure 9).tains an RMSD of 0.39, with fluctuations also within the 0.2 nm range.As for NFKIBA, both selected ligand Luteolin and Acacetin display higher overall fluctuations.In the last 5 ns, their RMSDs tend to stabilize, with fluctuations smaller than 0.2 nm.However, compared to other complexes, their stability remains poor, and thus, they cannot be considered to have achieved stable binding.Based on this, further analysis will be conducted of the following compounds: Involucratin and Luteolin for MAPK1, Xuelianlactone and Acacetin for TNF, and Flazin for MAPK8.
Solvent-Accessible Surface Area (SASA) quantifies the surface area of a molecule that is accessible to solvent molecules and can indicate changes in a molecule's conformation or exposure during simulation, aiding in the assessment of binding interactions in a target protein-ligand complex.The average SASA values for complexes XT (black), AT (red), IM (blue), LM (green), and FM (purple) during MD simulation are as follows: 202.49± 3.72, 201.65 ± 3.17, 179.03 ± 3.25, 175.97 ± 2.59, and 179.9 ± 2.9 nm 2 .Throughout the entire simulation process, the SASA values of these compounds show small fluctuations, indicating their stable interaction with the protein complexes (Figure 9).Molecular Mechanics-Poisson Boltzmann Surface Area )MMPBSA) is a method used to estimate the binding energy of a target protein-ligand complex by post-processing the MD trajectory.MMPBSA can characterize the strength and stability of the protein-ligand interaction by calculating the difference between the free energy of the bound state and the free energy of the unbound state.For the five complexes mentioned above, we selected Molecular Mechanics-Poisson Boltzmann Surface Area (MMPBSA) is a method used to estimate the binding energy of a target protein-ligand complex by post-processing the MD trajectory.MMPBSA can characterize the strength and stability of the protein-ligand interaction by calculating the difference between the free energy of the bound state and the free energy of the unbound state.For the five complexes mentioned above, we selected segments with stable RMSD values and extracted 10 ns trajectories for MMPBSA energy analysis.The total binding energy (∆G) is the sum of the change in enthalpy (∆H) and the change in entropy (-T∆S) [49].Using gmx_mmpbsa.py with a time interval of 10 ps, the total energy difference (∆G total ) between the ligand, receptor and their respective complexes in the stable phase of MD simulation was calculated, as shown in Figure 10.It can be seen from the figure that the ∆G total curves of each compound have small fluctuations, and the average is below −10 kcal/mol.Furthermore, by utilizing the gmx_mmpbsa.bshscript, the entropy change in each complex is calculated at a time interval of 1000 ps, allowing for the calculation of the total binding energy (∆G).Table 6 presents the contributions of the Van der Waals force (VDW), coulomb interaction energy (COU), molecular mechanics (MM), polar solvation energy (PB) and nonpolar solvation energy (SA).The ∆H calculated by gmx_mmpbsa.bshclosely aligns with the ∆G total computed using gmx_mmpbsa.py.The results from both methods exhibit a close resemblance, which enhances the credibility of the computed free energy values.It can be observed from the figure and the table that ∆G and ∆H for XT are significantly lower than that of the other ligands, consistent with the results observed in molecular docking.
tions of the Van der Waals force (VDW), coulomb interaction energy (COU), molecular mechanics (MM), polar solvation energy (PB) and nonpolar solvation energy (SA).The ΔH calculated by gmx_mmpbsa.bshclosely aligns with the ΔGtotal computed using gmx_mmpbsa.py.The results from both methods exhibit a close resemblance, which enhances the credibility of the computed free energy values.It can be observed from the figure and the table that ΔG and ΔH for XT are significantly lower than that of the other ligands, consistent with the results observed in molecular docking.This study investigated the binding free energy and residue contributions of active components of SI to the RA target using methods such as MD simulations and MMPBSA energy calculations.As shown in Figure 11, among the four complexes, Van der Waals forces have the most significant impact in promoting protein-ligand binding, while nonpolar solvation energy has a relatively weaker influence on binding affinity.On the other  This study investigated the binding free energy and residue contributions of active components of SI to the RA target using methods such as MD simulations and MMPBSA energy calculations.As shown in Figure 11, among the four complexes, Van der Waals forces have the most significant impact in promoting protein-ligand binding, while nonpolar solvation energy has a relatively weaker influence on binding affinity.On the other hand, polar solvation energy strongly impedes the formation of the complexes.Specifically, for the MAPK1 target, we compared the results of MD simulations for the IM and LM complexes and identified four residues (36, 39, 156, and 166) that played important roles in both complexes.In the LM complex, the sum of the ∆H values for these four residues accounted for 42.9% of the total ∆H.The maximum resistance to binding was observed from residues such as 54, 151, 153, and 167, but it was insufficient to offset the binding energy.For the TNF target, in both the XT and AT complexes, four residues (133, 135, 231, and 233) made significant contributions to the binding free energy.In the AT complex, the sum of the ∆H values for these four residues accounted for 43.0% of the total ∆H, while the maximum resistance to binding came from residues 174 and 195, although their impact was relatively small, with residue 174 contributing only 0.687 kJ/mol to ∆H.Compared to complexes IM and LM, complexes AT and XT exhibit higher binding energy contributions from residues and a lower binding resistance.This observation aligns with what is reflected in the RMSD curves.
complex, the sum of the ΔH values for these four residues accounted for 43.0% of the total ΔH, while the maximum resistance to binding came from residues 174 and 195, although their impact was relatively small, with residue 174 contributing only 0.687 kJ/mol to ΔH.Compared to complexes IM and LM, complexes AT and XT exhibit higher binding energy contributions from residues and a lower binding resistance.This observation aligns with what is reflected in the RMSD curves.The analysis of residue interactions with different ligands targeting the same receptor reveals crucial regions that facilitate complex binding.In future drug design efforts, it is essential to pay attention to the contributions of these residues.This can be achieved by introducing functional groups or optimizing their molecular spatial structures, aiming to retain residues with strong binding forces and reduce the impact of residues with significant resistance [50].

Discussion
In this study, we employed network pharmacology to analyze SI, resulting in the identification of 27 bioactive components and 119 targets that intersect with RA.Among these targets, 33 were identified as key targets.Additionally, we conducted KEGG and GO analyses to investigate the pathways and biological processes involved.The findings suggest a potential correlation between the components present in SI and the treatment of RA.The analysis of residue interactions with different ligands targeting the same receptor reveals crucial regions that facilitate complex binding.In future drug design efforts, it is essential to pay attention to the contributions of these residues.This can be achieved by introducing functional groups or optimizing their molecular spatial structures, aiming to retain residues with strong binding forces and reduce the impact of residues with significant resistance [50].

Discussion
In this study, we employed network pharmacology to analyze SI, resulting in the identification of 27 bioactive components and 119 targets that intersect with RA.Among these targets, 33 were identified as key targets.Additionally, we conducted KEGG and GO analyses to investigate the pathways and biological processes involved.The findings suggest a potential correlation between the components present in SI and the treatment of RA.
Among the core targets, ERK (MAPK1), JNK (MAPK8), and p38 MAPK (MAPK14) are well-defined families of MAPKs, activated by phosphorylation, and play significant roles in the inflammatory and destructive mechanisms observed in rheumatoid arthritis, including regulating pro-inflammatory cytokine production and mediating downstream signaling from IL-1, IL-17, and TNF-α receptors.These properties make them attractive therapeutic targets for rheumatic diseases [51,52].RELA, NFKBIA, CHUK, and IKBKB are associated with the NF-κB pathway.CHUK and IKBKB serve as different catalytic subunits of the IκB kinase (IKK) complex [53,54].Highly activated NF-κB induces the production of proinflammatory cytokines, including TNF-α, IL-1β, and IL-6, accelerating RA progression.The upregulation of these cytokines also triggers positive feedback, regulating NF-κB activation, and forming a vicious cycle that worsens RA development [52].Among these targets, MAPK targets exhibit a higher degree of centrality in network analysis.According to molecular docking results, both MAPK1 and MAPK8 demonstrate prominent binding efficiency with all core components.Therefore, MAPK targets play a crucial role in the process of treating RA with IS.Regarding NF-κB-related targets, their network analysis shows lower degree centrality rankings, and taking NFKIBA as an example, the results of molecular dynamics simulations fail to confirm its stable binding.Therefore, it can be inferred that the active components of SI have a weaker effect on these targets.
From a pathway perspective, the pathways exhibiting high gene enrichment and confidence levels include IL-17, TNF, TLR (Toll-like receptor), NLR (NOD-like receptor), and CLR (C-type lectin receptor) signaling pathways, as well as Th17, Th1, and Th2 cell differentiation pathways, among which the IL-17 and Th17 pathways are the most confident.Interactions between RA, Fibroblast-Like Synoviocytes (FLSs), and Th17 cells further contribute to tumorous FLS growth and joint pannus formation [55].Previous research has identified IL-17 as a potential therapeutic target for RA, leading to several trials being conducted [56,57].There is also considerable confidence in the osteoclast (OC) differentiation pathway, which plays a key role in the pathological bone destruction observed in RA.Pro-inflammatory cytokines, including TNF, IL-1 and IL-17, induce OC maturation and differentiation through different signaling pathways such as NF-κB [58].This highlights the dual therapeutic role of SI, not only targeting the inflammatory symptoms of RA but also reducing bone destruction.
The GO analysis primarily centers on responses to biotic stimuli, specifically lipopolysaccharides and other molecules of bacterial origin, known to trigger immune and inflammatory responses.Previous research has indicated that bacterial infection plays a significant role in contributing to RA pathogenesis [59].The pathways showing the highest enrichment, such as Coronavirus disease, Chagas disease, and Yersinia infection, further reinforce this notion, as they are involved in immune responses to viruses, parasites, and bacteria.In terms of molecular function analysis, the active components primarily influence functions related to binding activities, such as cytokine receptor binding, growth factor receptor binding, and phosphatase binding, while also affecting cytokine activity and receptor ligand activity.This suggests that these components may modulate processes like phosphorylation within the pathways by binding to specific targets, thereby inhibiting excessive immune responses and achieving therapeutic effects.Additionally, in the biological process analysis, certain processes are associated with the regulation of cell-cell adhesion.Adhesion molecules are important regulators of leukocyte recruitment into the synovial tissue [60], while Leukocytes, which are aberrantly recruited and activated, are characteristic features of rheumatoid arthritis (RA), as they migrate out of the blood, cross the endothelium, and infiltrate inflamed tissues, triggering a cascade of pathological processes, particularly in the early stages of the disease [61].
Based on the above analysis results, components Luteolin, Involucratin, Xuelianlactone, Flazin, Quercetin, and Acacetin, all belonging to the flavonoid class, exhibit favorable binding effects with the targets.Quercetin exerted a bone-protective effect by reducing MMPs, RANKL production, and osteoclast formation via the MAPKs and NF-κB pathways [62].Current animal model studies have demonstrated the effectiveness of quercetin in treating RA [63].Some studies suggest that Luteolin can inhibit the proliferation and affect the function of stimulated rat synovial fibroblasts [64].Furthermore, it has the ability to suppress the IL-1β-induced production of cytokines and MMPs, which are crucial in tissue degradation in rheumatoid synovium, through the activation of p38 MAPK, JNK, NF-κB, and AP-1 [65].Acacetin inhibits p38 and JNK phosphorylation and reduces MMP-1, MMP-3, and MMP-13 expression in IL-1β-induced FLSs, suggesting that Acacetin has antiarthritic effects in FLSs [66].In addition, it has been shown to inhibit the activation of NF-κB by stimulation with TNF-α [67].Moreover, this study reveals the significant roles of Involucratin, Xuelianlactone, and Flazin in the process of treating RA with SI.Among all complexes, Xuelianlactone exhibits the best binding efficiency with the TNF target, displaying the lowest binding energy in both molecular docking and MD simulations, along with good binding stability.Involucratin, on the other hand, stands out in the molecular docking results, showing the lowest binding energy among ligands binding to MAPK1, MAPK8, and RELA.Flazin also demonstrates strong binding abilities with MAPK8 in MD simulations.Considering the limited number of related studies conducted on these compounds to date, they represent potential new research avenues for the development of RA treatment drugs in the future.

Conclusions
In summary, the network pharmacology analysis suggests that the mechanism of SI in treating RA mainly involves its flavonoids, which target key proteins through MAPK and NF-κB-related pathways.By doing so, they effectively inhibit inflammation and mitigate bone destruction, leading to therapeutic effects.The binding efficiency of complexes like XT, AT, and FM has been confirmed through MD simulations.This suggests that the action of SI in treating RA is achieved through multiple components acting on various targets.However, it is essential to note that network pharmacology is primarily a data-driven and network-based research method.Therefore, these findings should be further validated and supported by experimental studies to establish their clinical relevance and potential as RA treatment options.

Figure 1 .
Figure 1.Common targets of SI and RA.

Figure 1 .
Figure 1.Common targets of SI and RA.

Figure 6 .
Figure 6.Heat map of lowest binding energy in molecular docking.

Figure 6 .
Figure 6.Heat map of lowest binding energy in molecular docking.

Figure 7 .
Figure 7. Molecular docking results in each target with the bioactive compound of SI: (a) Luteolin and MAPK14; (b) Involucratin and MAPK1; (c) Involucratin and RELA; (d) Xuelianlactone and TNF; (e) Involucratin and MAPK8; (f) Quercetin and IL6; (g) Flazin and IL1B; (h) Flazin and CHUK; (i) Quercetin and IKBKB; (j) Luteolin and NFKBIA.The images depict the binding positions of ligands within the protein, along with the surrounding residues, showcasing the interactions between the ligands and the residues.In this figure, pink, magenta, and purple dashed lines represent hydrophobic interactions, while the green dashed lines represent hydrogen bonds.

Figure 7 .
Figure 7. Molecular docking results in each target with the bioactive compound of SI: (a) Luteolin and MAPK14; (b) Involucratin and MAPK1; (c) Involucratin and RELA; (d) Xuelianlactone and TNF; (e) Involucratin and MAPK8; (f) Quercetin and IL6; (g) Flazin and IL1B; (h) Flazin and CHUK; (i) Quercetin and IKBKB; (j) Luteolin and NFKBIA.The images depict the binding positions of ligands within the protein, along with the surrounding residues, showcasing the interactions between the ligands and the residues.In this figure, pink, magenta, and purple dashed lines represent hydrophobic interactions, while the green dashed lines represent hydrogen bonds.

Figure 8 .
Figure 8.The fluctuation plot of the target protein-ligand complexes RMSD.(a) MAPK1, (b) TNF, (c) MAPK8, (d) NFKIBA.38Z is the original ligand for MAPK1 and MAPK8, while VGY is for TNF, and NFKIBA does not have an original ligand.

Figure 8 .
Figure 8.The fluctuation plot of the target protein-ligand complexes RMSD.(a) MAPK1, (b) TNF, (c) MAPK8, (d) NFKIBA.38Z is the original ligand for MAPK1 and MAPK8, while VGY is for TNF, and NFKIBA does not have an original ligand.

Figure 10 .
Figure 10.The fluctuation plot of the target protein-ligand complexes ΔH.

Figure 10 .
Figure 10.The fluctuation plot of the target protein-ligand complexes ∆H.

Table 1 .
Bioactive components of SI.

Table 2 .
List of Key Targets.

Table 4 .
Best protein-ligand complexes and their binding energies.

Table 4 .
Best protein-ligand complexes and their binding energies.

Table 5 .
The Average RMSD values for each target protein-ligand complex.

Table 5 .
The Average RMSD values for each target protein-ligand complex.

Table 6 .
The ΔG of each complex and the contributions of individual energy components(kcal/mol).Note: In this table, MM is the sum of VDW and COU, while ΔH is the sum of MM, PB, and SA.

Table 6 .
The ∆G of each complex and the contributions of individual energy components(kcal/mol).Note: In this table, MM is the sum of VDW and COU, while ∆H is the sum of MM, PB, and SA.