Integrative Multi-omics Analysis of Buti Huatan Tang in Chronic Obstructive Pulmonary Disease

Summary

Multi-omics profiling of Buti Huatan Tang (BTHTT) in chronic obstructive pulmonary disease reveals its potential to modulate the IL-1R2/IL-1β-thyroid hormone metabolic axis. This work provides a framework for further mechanistic validation of the bioactive components of BTHTT, establishing a foundation for future causal intervention studies required to confirm this axis.

Abstract

This study utilized a multi-omics and computational biology framework to investigate the therapeutic potential of the Traditional Chinese Medicine (TCM) formula Buti Huatan Tang (BTHTT) against chronic obstructive pulmonary disease (COPD). Significant physiological improvements were observed in a rat model following BTHTT intervention. Histological analysis showed a reversal of lung pathological damage, while biochemical assays, and transcriptomics confirmed the normalization of IL-1β and IL-1R2 levels. Additionally, metabolic profiling revealed that BTHTT corrected disruptions in T3 and T4 thyroid hormone levels. A negative correlation was observed between the IL-1β/IL-1R2 axis and these thyroid hormones, indicating that their regulation is associated with the formula’s therapeutic effect. Beyond direct measurements, machine learning algorithms identified ten COPD signature genes from clinical databases. Pathway enrichment analysis suggests that BTHTT may act through cytokine-cytokine-receptor interactions and thyroid hormone synthesis pathways. Furthermore, while 283 components were identified in vivo, compounds such as tanshinone IIA and cryptotanshinone are currently considered candidate active substances. Their role as primary drivers is supported by a model in which they stably bind to IL-1R2; this inference is based on molecular docking and molecular dynamics (MD) simulations rather than direct experimental isolation. Overall, the data support a model in which BTHTT exerts a multi-target effect on COPD by modulating inflammation and metabolic homeostasis. This integrated approach provides a refined scientific basis for the clinical application of BTHTT and highlights specific pathways for future experimental validation.

Introduction

Chronic obstructive pulmonary disease (COPD), one of the most common respiratory diseases in clinical practice, has become the third leading cause of death globally due to its extremely high disability rate and recurrent acute exacerbations¹. COPD was a significant factor leading to the loss of labor capacity in the population and a sharp increase in family care costs, with social and economic losses far exceeding those of most chronic diseases². In the “Healthy China Initiative 2030,” COPD has been listed as a key disease for prevention and treatment. Due to an incomplete understanding of the mechanisms underlying COPD, the combined use of bronchodilators and corticosteroids remained the primary treatment. While this approach provided relatively obvious symptom relief, it cannot prevent the progressive decline in lung function. Moreover, the long-term use of corticosteroids inevitably increases the incidence of adverse reactions, further affecting the long-term efficacy of the drugs and the patients' tolerance³. Studies have shown that even with triple therapy, 28%–31% of patients still experience acute exacerbations or worsening of symptoms. Therefore, the development of effective drugs for COPD remains a key and primary issue in current COPD research⁴.

To address these needs, classic TCM formulas such as Buti Huatan Tang (BTHTT), Shengmai Yin, and Bu Fei Tang have been developed. Among them, BTHTT was particularly notable for its significant effects in tonifying the lung, boosting energy, and resolving phlegm to relieve cough^5,6,7,8. BTHTT was widely used in the clinical treatment of COPD, not only for its remarkable efficacy and minimal side effects, but also for its potential to reverse the progressive decline in lung function. It held great promise in clinical and social contexts and had the potential for “secondary development” as a major drug, with significant scientific and economic value^9,10,11,12. However, due to the lack of systematic, long-term scientific research, studies on BTHTT have generally suffered from an unclear understanding of its component basis and unidentified core targets of action. Fundamentally, this was due to significant deficiencies in research on the main mechanism and potential therapeutic substance basis of BTHTT.

In response to this characteristic, integrating multiple omics technologies such as metabolomics and transcriptomics, as well as multi-algorithm mining methods, focusing on the main correlations between characteristic key genes and key metabolic pathways in the body's microenvironment^13,14,15 cannot only explore the essence of treatment but also evaluate the synergistic efficacy of the active components. A particular focus is placed on the thyroid hormone axis, as emerging research indicates that thyroid dysfunction is a frequent comorbidity in COPD, often correlating with increased systemic inflammation and impaired respiratory muscle function. Modulating this axis may represent a novel pathway for restoring metabolic homeostasis in pulmonary disease^16,17.

In summary, this study adopted a three-pronged approach combining multi-omics integration with computational biology, based on the premise of syndrome and formula correspondence. On the one hand, it utilized TCM serum pharmacochemistry techniques to identify the blood-borne components of BTHTT under effective conditions. On the other hand, focusing on key characteristic genes of COPD, it used comprehensive tissue metabolomics to discover disease biomarkers from endogenous small molecules and to construct a precise evaluation of the overall effects of the formula. Meanwhile, combined with multi-algorithm mining and transcriptomic identification analysis, it identified the key characteristic genes regulated by BTHTT under effective conditions and conducted correlation analysis between the key characteristic genes and key metabolic pathways. The potential therapeutic substance basis was screened by molecular docking and molecular dynamics analysis.

Based on our multi-omics framework and the pharmacological profile of BTHTT, this study proposes the following testable predictions. Antiinflammatory modulation: BTHTT intervention will significantly reduce systemic and pulmonary inflammation markers, specifically by normalizing the IL-1β/IL-1R2 expression axis to mitigate tissue damage. Metabolic restoration: BTHTT will alter signals in the thyroid hormone synthesis pathway, leading to a measurable correction of T3/T4 levels and a subsequent improvement in the metabolic state of the lung microenvironment. Metabolic Restoration: BTHTT will alter signals in the thyroid hormone synthesis pathway, leading to a measurable correction of T3/T4 levels and a subsequent improvement in the metabolic state of the lung microenvironment.

Protocol

NOTE: The reagents and the equipment used in this study are listed in the Table of Materials.

Clinical data acquisition and preprocessing
To investigate the molecular signatures of COPD, transcriptomic datasets were retrieved from the Gene Expression Omnibus (GEO) database. A total of five datasets were selected: GSE11784, GSE12472, GSE16972, GSE38974, and GSE222965. The raw data and platform annotation files were downloaded for probe-to-gene mapping. When multiple probes targeted a single gene, the maximum expression value was retained. The resulting gene expression matrices (rows as genes, columns as samples) were merged into a single discovery dataset. To account for technical variations across different microarray platforms and study cohorts, batch-effect correction was performed using the ComBat algorithm from the sva R package. The efficacy of the correction was validated via Principal Component Analysis (PCA) plots. Following batch correction, the merged discovery dataset was utilized for differential analysis. Differentially expressed genes (DEGs) between COPD patients and healthy controls were identified using the limma package. The significance thresholds were set at |logFC| ≥ 1 and p ≤ 0.05.

To identify the most robust characteristic genes, three independent machine learning algorithms were integrated. To ensure the reliability of the models and prevent data leakage, the feature selection process was nested within the cross-validation loops where applicable, and the discovery dataset was strictly separated from the independent validation sets. LASSO model was applied to the DEGs using the glmnet package. We employed 10-fold cross-validation to determine the optimal penalty parameter. The optimal penalty parameter value corresponding to the minimum cross-validation error was selected as the threshold to identify core feature genes. SVM was utilized to rank genes based on their discriminative power. A 10-fold cross-validation strategy was implemented to identify the point of minimum generalization error, thereby determining the optimal number of characteristic genes. The randomForest package was used to rank DEGs based on their Mean Decrease Accuracy and Gini index. Genes with the highest importance scores were selected as disease-related features.

The intersection of the features identified by LASSO, SVM, and randomForest was taken to define the final core characteristic genes. The diagnostic performance of these genes was evaluated using Receiver Operating Characteristic (ROC) curve analysis. The Area Under the Curve (AUC) and associated 95% Confidence Intervals (CIs) were calculated using the pROC package. A gene was considered to have high diagnostic value if the AUC > 0.70. Finally, the expression levels and diagnostic accuracy of these genes were further validated in an independent validation set to ensure the generalizability of the findings.

Using the PubChem (https://pubchem.ncbi.nlm.nih.gov/) and drugbank (https://go.drugbank.com/), a total of 12 specific drug profiles were obtained, such as LABAs, LAMAs, SABAs, SAMAs, Fluticasone propionate, Budesonide, Beclomethasone, Fluticasone, Salmeterol, Umeclidinium, Vilanterol, Theophylline. The compiled targets were calibrated using the Uniprot database (https://www.uniprot.org/), during which non-human genes were removed, and invalid duplicate targets were deleted to obtain standardized gene names. By entering the keyword “chronic obstructive pulmonary disease”, “COPD” in the GeneCards (https://www.genecards.org/), OMIM (https://www.omim.org/), and DisGENET (https://www.disgenet.org/) databases, disease-related targets were retrieved. All targets from the three databases were consolidated into an Excel file, duplicate genes were removed, and the data was calibrated using the Uniprot database to obtain the final disease target gene information.

Machine learning model development
A multi-algorithmic machine learning framework was constructed by sequentially utilizing R packages glmnet, e1071, and randomForest. Specifically, LASSO regression was performed for penalty-based dimensionality reduction, Support Vector Machine (SVM) analysis was employed to evaluate validation errors based on sample grouping, and Random Forest was used to filter features according to importance scores. This process yielded corresponding diagnostic visualizations, including cross-validation curves and gene importance bubble plots. Following model construction, a Venn diagram analysis was performed on the gene sets identified by these multiple algorithms to extract the overlapping "intersection" feature genes, thereby enhancing the reliability of the potential biomarkers. The expression matrix of these intersection genes was then extracted to visualize inter-group expression differences via violin plots. Finally, ROC curves were generated through iterative loops for each gene to calculate the Area Under the Curve (AUC), validating their diagnostic value as candidate biomarkers.

Animal model induction
The experimental protocol was approved by the Animal Ethics Committee of Guizhou Medical University (2303411) and followed the ARRIVE guidelines and animal welfare regulations. A total of 72 male Sprague-Dawley (SD) rats of SPF grade (230 ± 20 g). After 1 week of adaptive housing under standard conditions (25 ± 1 ℃, 50 ± 5% humidity, 12 h light-dark cycle), the rats were randomly divided into the Control group (C group) (n = 24) and the Model group (M group) (n = 48).

The M group was subjected to dual-factor modeling: Intermittent cigarette smoke exposure (9 weeks, 6 days per week, 3 research-grade cigarettes per day divided into 2 sessions, 30 min per session) and intratracheal instillation of LPS (200 µg per instillation, once every other week)^18,19. The C group received an equivalent volume of normal saline. After successful modeling, the groups were randomly divided into the M group (n = 12), the high-dose BTHTT group (H group) (High, 1× clinical dose, n = 12), and the low-dose BTHTT group (L group) (Low, 1/2× clinical dose, n = 12), with continuous gavage intervention for 2 weeks. The C and M groups were given distilled water synchronously.

Body weight and spontaneous activity parameters were recorded weekly. At the end of modeling and treatment, lung tissue, serum, and serum and bronchoalveolar lavage fluid (BALF) were collected. Levels of inflammatory markers (i.e., TNF-α, IL-1β, IL-6, IL-8, OPN, and MCP-1) were measured using ELISA according to the kit instructions. Tissue pathology and multi-omics analyses were also performed.

Treatment administration
Buti Huatan Tang (BTHTT) consists of nine traditional Chinese medicinal herbs: Astragalus membranaceus (Fisch.) Bge. (Astragali Radix, 15 g), Pseudostellaria heterophylla (Miq.) Pax (Pseudostellariae Radix, 15 g), Cinnamomum cassia (L.) J. Presl (Cinnamomi Cortex, 15 g), Angelica sinensis (Oliv.) Diels (Angelicae Sinensis Radix, 10 g), Salvia miltiorrhiza Bge. (Salviae Miltiorrhizae Radix et Rhizoma, 15 g), Perilla frutescens (L.) Britt. (Perillae Folium, 10 g), Raphanus sativus L. (Raphani Semen, 10 g), Lepidium apetalum Willd. (Descurainiae Semen, 10 g), and Mahonia fortunei (Lindl.) Fedde (Mahoniae Folium, 10 g). The herbal materials were soaked in water at 10x their combined weight for 30 min and subsequently decocted for 1 h. The decoction was filtered, and the filtrate was collected and divided into three equal portions for oral administration.

Based on clinical safety and efficacy, the standard dose of BTHTT for adults was 1.57 g∙kg^-1∙day^-1 according to the clinical medication guidelines. Considering the dose conversion factor for rats, which indicated that the standard rat drug dose was 6.3× the human standard dose, the corresponding H group for rats was set at 9.9 g∙kg-1∙day-1, and the L group was set at half this dose (4.95 g∙kg^-1∙day^-1). Given the final concentrated volume of BTHTT was 50 mL, the administered volume for the H group in rats was approximately 1.6 mL, and for the L group, approximately 0.8 mL. The drug was administered once daily via oral gavage.

Tissue and BALF collection
At the experimental endpoints (Week 9 and Week 11), rats were anesthetized via intraperitoneal injection of 5% sodium pentobarbital (1 mL/100 g). Blood was collected from the portal vein, allowed to stand for 30 min at room temperature, and then centrifuged at 13,000 × g for 15 min at 4 ℃. The supernatant was stored at -80 ℃. Lung tissues were rapidly frozen in liquid nitrogen and stored at -80 ℃. BALF was collected by three consecutive lavages with cold PBS.

Transcriptomics analysis
Total RNA was extracted from lung tissue and quality-controlled using (A₂₆₀/A₂₈₀ > 1.8) (RIN ≥ 7.0. RNA sequencing libraries were constructed: mRNA was enriched with oligo(dT), fragmented, and used to synthesize double-stranded cDNA, which was then ligated to adapter oligonucleotides and amplified by PCR. Libraries were quantified and quality-checked before sequencing. Including content distribution inspection, FPKM density distribution analysis of each sample, and overall quality assessment analysis of RNA-seq²⁰.

Raw transcriptomics data generated by the sequencing platform were processed using Perl scripts to remove adapter sequences and low-quality reads (reads with Q ≤ 25 bases accounting for > 60% or N rate > 5%). Clean reads were obtained after this filtering process. The clean reads were aligned to the reference genome using HISAT2, and gene expression levels were quantified to calculate FPKM values. Differential expression analysis was performed (with screening criteria of |log2FC| > 1 and p-adj < 0.05). Transcription factor annotation was based on Animal TFDB or Pfam/DBD databases, matching gene IDs and protein domain information²¹.

An orthologous gene mapping strategy was employed to ensure methodological rigor for cross-species validation. This strategy involved retrieving orthologs in rats for core human genes (e.g., SRPX2, IL-1R2, TFF3) using the NCBI HomoloGene and Ensembl BioMart databases. The selection was restricted to gene pairs exhibiting a clear "one-to-one" mapping relationship and high protein sequence identity. In instances where multiple candidates were present, preference was given to orthologous pairs certified by the HGNC. To ensure high detection precision, specific RT-qPCR primers were designed based on the mRNA sequences of the identified rat orthologs. Validation criteria were defined by consistency in expression direction and functional verification. Consistency in Expression Direction: In a smoke-induced rat COPD model, RT-qPCR showed that the expression trends of target genes in rat lung tissue were fully consistent with those observed in human GEO clinical datasets. Genes exhibiting the same polarity of change across species were considered to possess conservation as disease biomarkers. Pathological Verification: Upon expression consistency, correlation analysis was performed to verify the involvement of these genes in the pathological evolution of COPD.

Metabolomics workflow
Lung tissue (20–50 mg) was homogenized in prechilled methanol-acetonitrile-water (2:2:1, v/v) and sonicated in an ice bath. The homogenate was centrifuged 13,000 × g for 20 min at 4 ℃. The supernatant was vacuum-concentrated, re-dissolved in acetonitrile-water (1:1, v/v), and filtered through a 0.22 µm membrane for LC-MS analysis²².

Chromatographic separation was achieved using an amide-based UPLC column (1.7 µm, 2.1 mm × 100 mm). Column temperature was maintained at 25 ℃. The mobile phase consisted of (A) water containing 25 mM ammonium acetate and 25 mM ammonia, and (B) acetonitrile. Flow rate was set at 0.5 mL/min, and the injection volume was 2 µL. Gradient elution program was as follows: 0–0.5 min, 95% B; 0.5–7 min, linear decrease of B from 95% to 65%; 7–8 min, linear decrease of B from 65% to 40%; 8–9 min, B held at 40%; 9–9.1 min, linear increase of B from 40% to 95%; 9.1–12 min, B held at 95%. During the entire analysis, samples were kept at 4 ℃ in the autosampler. To ensure system stability and the reliability of the experimental data, samples were analyzed in a random sequence, with quality control (QC) samples interspersed in the queue. Mass spectrometry analysis was performed using an ultra-high-performance liquid chromatography (UHPLC) system coupled to a mass spectrometer. Samples were ionized using electrospray ionization (ESI) in both positive and negative ion modes. ESI source and MS settings were as follows: nebulizer gas (Gas 1) was set to 50, auxiliary gas (Gas 2) to 2, ion source temperature to 350 ℃, and spray voltage (ISVF) to 3,500 V in positive ion mode and 2,800 V in negative ion mode. The mass range for MS1 was set from 70 to 1,200 Da, with a resolution of 60,000 and a scan accumulation time of 100 ms. For MS2, data-dependent acquisition (DDA) with stepped collision energy was used. The mass range for MS2 was also set from 70 to 1,200 Da, with a resolution of 60,000 and a scan accumulation time of 100 ms. The dynamic exclusion time was set at 4 s.

Raw metabolomics data were converted to the mzXML format and then processed for peak alignment, retention time correction, and peak area extraction. Data preprocessing workflow included the following steps: First, ion peaks with a missing rate > 50% were removed. Second, the remaining missing values were imputed using the KNN algorithm. Third, metabolic features with a relative standard deviation (RSD) >50% were discarded. The quality of the experimental data was assessed using principal component analysis (PCA) and clustering of QC samples. Subsequent analyses included univariate statistics (e.g., t-tests), multivariate statistics (PLS-DA), differential metabolite screening (VIP > 1 and p < 0.05), and KEGG pathway enrichment analysis (hypergeometric test)^23,24.

Establishment of analytical methods for in vivo and in vitro components
Preparation of BTHTT samples for in vitro testing: BTHTT was extracted by decoction in water (2 x 30 min), concentrated to 1.1–1.2 g/mL, and then freeze-dried. Before analysis, 600 µL of the freeze-dried powder solution was mixed with 400 µL of methanol, re-dissolved in 40% methanol, and centrifuged to collect the supernatant.

Preparation of BTHTT samples for in vivo testing: Serum was deproteinized by mixing with methanol (1:1) and precipitating at −20 ℃ for 30 min, followed by centrifugation at for 20 min. The supernatant was vacuum-dried and re-dissolved in 40% methanol to obtain the final sample. For the preparation of blank serum + BTHTT samples, an appropriate amount of blank serum was spiked with the in vitro BTHTT supernatant, and the remaining steps were performed as described.

Samples were separated using a UHPLC system equipped with a reverse-phase UPLC column (2.1 mm × 100 mm, 1.8 µm). The column temperature was maintained at 35 °C, and the flow rate was set at 0.3 mL/min. The mobile phase consisted of (A) 0.1% formic acid in water and (B) 0.1% formic acid in acetonitrile. Gradient elution was performed as shown in Table 1.

A mass spectrometer was used for the acquisition of both MS1 and MS2 spectra. The mass spectrometer was coupled with the UHPLC system and operated in both positive and negative ESI modes. ESI parameters were as follows: spray voltage 3,800 V (ESI+) / 3500 V (ESI-), sheath gas pressure 45 arb, auxiliary gas pressure 20 arb, ion transfer tube temperature 320 °C, and vaporizer temperature 350 °C. The detection mode was set to full scan/data-dependent MS2 (Full-MS/dd-MS2) with resolutions of 60,000 for MS1 and 15,000 for MS2. The top 10 MS1 ions were selected for MS/MS fragmentation with stepped normalized collision energies of 20, 40, and 60. Mass range for MS1 was set from 90 to 1,300 Da.

For in vivo analysis, including blank group samples, dosed group samples, and blank group + BTHTT samples, 6 µL of each sample was precisely injected. For in vitro analysis of BTHTT, 2 µL of the sample was injected. Each batch of blank and dosed group samples was injected once, while the blank group + BTHTT samples were injected in triplicate, and the BTHTT samples were injected in quintuplicate.

Data in the mzXML format were processed and compounds were identified based on a local high-resolution commercial TCM mass spectrometry database. The criteria for identification were set as follows: a mass error of <25 ppm for MS1 and a match score > 0.7 for MS2 (where the Score reflected the similarity of fragment ions, with ≥0.7 being a reliable threshold)^25,26. Statistical analysis included compound counting and classification (e.g., flavonoids, alkaloids), which was performed in conjunction with annotations from the mass spectrometry database²⁷.

Molecular docking and MD simulation
To investigate potential binding modes between the identified characteristic proteins and their corresponding ligands, in silico molecular docking was performed. The three-dimensional structures of the small molecules were retrieved from the PubChem database, and their geometric configurations optimized. The crystal structures of the target proteins were obtained from the RCSB Protein Data Bank (PDB). Using PyMOL, water molecules and heteroatoms were removed, and co-crystallized ligands were extracted to define the active site coordinates. Hydrogen atoms were added, and Gasteiger charges were assigned using software. Docking simulations were executed, generating 15 independent conformations per run. The conformation with the lowest binding energy was selected for further analysis. To rigorously characterize the non-covalent interactions, the receptor-ligand complexes were analyzed using the Protein-Ligand Interaction Profiler (PLIP). NOTE: It is important to emphasize that these docking results provide structural support for potential molecular interactions and serve as a foundation for further dynamic refinement; however, they do not constitute standalone evidence of biological efficacy.

To evaluate the stability and conformational evolution of the predicted protein-ligand complexes under physiologically relevant conditions, molecular dynamics simulations were performed using the GROMACS software package. Topology files for both the proteins and ligands were generated based on the GROMOS96 43a1 force field. Each complex was positioned at the center of a dodecahedral box, maintaining a minimum distance of 1.0 nm from the box edges, and solvated using the SPC water model. To ensure electrical neutrality, sodium or chloride ions were added to the system as required. Energy minimization was performed using the steepest descent algorithm until the maximum force was less than 1,000.0 kJ∙mol^-1∙nm^-1. The system was then equilibrated in two stages: first, an NVT ensemble was employed to heat the system to 300 K over 100 ps using a V-rescale thermostat; second, an NPT ensemble was used to stabilize the pressure at 1 bar over 100 ps using a Parrinello-Rahman barostat. Production simulations were carried out for a total duration of 10 ns with a time step of 2 fs. Long-range electrostatic interactions were calculated using the Particle Mesh Ewald (PME) method, while short-range van der Waals and electrostatic interactions were managed with a cutoff radius of 1.2 nm. To ensure the reliability of the simulations, three independent runs were performed where feasible. The stability of the complexes was quantitatively assessed by calculating the Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuation (RMSF) of the protein backbone atoms relative to the initial structure. The attainment of a plateau in the RMSD profile was utilized as the primary criterion for system equilibration and structural stability.

General statistical analysis
In this experiment, group calculations were performed using t-tests or one-way analysis of variance (ANOVA).

Results

Prediction of key characteristic genes in COPD
Following batch-effect correction of clinical data from multiple microarray platforms (Supplemental File 1–Supplemental Figure S1), a differential expression analysis was performed. 34 candidate genes (32 upregulated and 2 downregulated) were identified, including SRPX2, CLDN10, and TFF3 (Figure 1A).

To refine these candidates into a robust predictive signature, three machine learning algorithms were employed in shaping a joint model. For LASSO regression, using the λ value corresponding to the minimum cross-validation error, 16 features were selected (Figure 1B,C). For SVM, 30 characteristic genes were determined at the point of minimum screening error (Figure 1D,E), and 21 genes were identified based on their importance scores in random Forest (Figure 1F,G). By intersecting these results, 10 core characteristic genes (SRPX2, CLDN10, TFF3, HRASLS2, ALDH3A1, MUC1, MUC4, SERPINF1, IL-1R2, and FOSB) were identified (Figure 1H). These genes are considered predictive biomarkers for COPD diagnosis.

The diagnostic utility of the 10-gene signature was assessed using ROC analysis. In the discovery phase, all genes exhibited high diagnostic accuracy with Area Under the Curve (AUC) values > 0.7 (Figure 1J). To evaluate model generalizability and performance uncertainty, the biomarkers were tested against an independent validation set. Differential expression was significantly maintained for IL-1R2, SRPX2, and TFF3 (Figure 1K). Notably, FOSB, IL-1R2, MUC4, SRPX2, and TFF3 sustained AUC values > 0.7, with IL-1R2 achieving a peak AUC of 0.816 (Figure 1L). These results suggest that while these genes are strong candidates for clinical stratification, their functional roles in COPD pathogenesis remain to be elucidated through future mechanistic studies.

To explore potential therapeutic intersections, we analyzed 12 drug profiles based on clinical literature and PubChem. A total of 280 potential therapeutic targets were identified. Simultaneously, 8,722 COPD-related targets were retrieved from GeneCards, OMIM, and DisGeNET. We identified 253 overlapping targets between the drug profiles and COPD-related genes (Figure 2A). It should be noted that these "DRUG-COPD" interaction targets are derived from indirect therapeutic drug-target mapping and represent a generalized intersection rather than drug-specific experimental binding. Protein-protein interaction (PPI) analysis (Confidence > 0.7) and topological screening (Relevance Score > 150) identified TERT, BMPR2, PKHD1, SERPINA1, and IL1B as central nodes (Figure 2B,C).

By integrating the machine learning-derived predictive biomarkers with the therapeutic network analysis, the IL-1R2/IL-1β axis was identified as a critical focal point for the progression and potential treatment of COPD.

Pathological and phenotypic manifestations of the COPD rat model
The COPD rat model was constructed by intratracheal instillation of LPS combined with exposure to cigarette smoke, and the model was evaluated multidimensionally from behavioral, physiological, and biochemical indicators, as well as histopathological aspects. Compared with the control group (C), the model group (M) rats exhibited a significant decrease in body weight (^**p < 0.01) (Figure 3A), and the number of activities in the M group significantly decreased compared with the premodeling (^**p < 0.01) (Figure 3B).

BALF tests indicated that the levels of pro-inflammatory factors (TNF-α, IL-6, IL-8, MCP-1, and OPN) were significantly elevated in the M group (Figure 3C). Histopathological observations revealed extensive inflammatory infiltration in the lung tissue of the M group, characterized by thickening of the alveolar walls with granulocyte infiltration (Figure 3D), lymphocyte aggregation around the bronchioles, and focal macrophage infiltration. Typical pathological changes included hydropic degeneration of the bronchial epithelium (pale cytoplasmic swelling), abnormal mucus secretion in the lumen, and proliferation of epithelioid cells forming cystic structures containing necrotic debris. No significant pathological changes were observed in the lung tissue of the C group. These results were highly consistent with the clinical characteristics and pathological mechanisms of human COPD. So the animal model with typical pathological features of COPD was successfully established in this study.

Transcriptomic validation of key characteristic genes in COPD
Based on the above results, further transcriptomic research was conducted on the COPD rat model to identify relevant, predictive results for key characteristic genes. As shown in Supplemental File 1–Supplemental Table S1 and Supplemental File 1—Supplemental Figure S2, the sequencing data and quality control met the experimental requirements. Compared with the C group, a total of 173 differentially expressed genes (DEGs) (p-adj < 0.05; |log2(fold change)| > 1) were screened in the M group, including 134 up-regulated genes (including IL-1β) and 39 down-regulated genes (including IL-1R2). The results are shown in Figure 4A and Supplemental File 1—Supplemental Table S2.

Subsequently, functional annotation analysis was performed on the differentially expressed genes. The KEGG pathway annotation results of the DEGs, shown in Figure 4B, revealed the strongest association with the cytokine-cytokine-receptor interaction pathway (p-adj < 0.01), where the IL-1R2 and IL-1β changing were probably one of the primary reasons for this strong association (Figure 4D). Further construction of an interaction network for the DEGs is shown in Figure 4C. Network analysis indicated that based on a centrality measure of degree ≥ 20, IL-1β was identified as a key hub node in the interaction network. Therefore, integrated transcriptomic analysis results further confirmed the accuracy of IL-1R2/IL-1β as key characteristic genes for COPD, and the cytokine-cytokine-receptor interaction pathway is presumed as a key signaling pathway in COPD.

Organizational metabolomics and identification of key metabolic pathways in COPD rat model
The total ion current (TIC) chromatograms of the QC samples are shown in Supplemental File 1—Supplemental Figure S3A,B. The response intensity and retention times of the chromatographic peaks largely overlapped, with correlation coefficients between QC samples exceeding 0.9, indicating good experimental reproducibility (Supplemental File 1-–Supplemental Figure S3C,D). The proportion of peaks with relative standard deviation (RSD) ≤ 30% in the QC samples accounted for over 70% of the total number of peaks in the QC samples, confirming the stability of the analytical system and the suitability of the data for further analysis (Supplemental File 1-–Supplemental Figure S3E,F).

PCA analysis revealed distinct clustering within groups and clear separation between the C and M groups, indicating successful induction of metabolic disturbances by the COPD model (Figure 5A). The OPLS-DA model effectively distinguished the M group from the blank group (Figure 5B), with Q2 values of 0.659 and 0.757 in the positive and negative ion modes, respectively (Q2 > 0.5), demonstrating the model's stability and reliability. Permutation tests (Figure 5C,D) showed that the R2 and Q2 values of the random models decreased with increasing permutation retention, confirming the absence of overfitting and the robustness of the original model. A total of 29 differential metabolites were identified based on VIP > 1 and p < 0.05 (Supplemental File 1–Supplemental Table S3), as shown in Figure 5E,F.

Pearson correlation analysis was performed on the identified small-molecule metabolites to explore metabolic relationships and regulatory interactions during biological state changes. Metabolites with correlation coefficients greater than 0.8 were presumed as key markers involved in synthesis and transformation. As shown in Figure 6A,B, glutathione (oxidized), phosphocholine, and L-palmitoylcarnitine were presumed as key markers according to the average correlation coefficient greater than 0.8. In addition, topological analysis of differential metabolic pathways was performed using the MetPA analysis method. Sixteen major metabolic pathways were identified by MetPA analysis, including tyrosine metabolism, one-carbon pool by folate, and pyrimidine metabolism (Figure 6C). The pathway impact analysis showed that phenylalanine, tyrosine, and tryptophan biosynthesis had an impact value greater than 0.5 in the COPD rat model. To avoid the loss of important signals due to filtering thresholds and to capture coordinated changes across the metabolite network, major metabolic pathways were further identified using metabolite set enrichment analysis (MSEA) (Figure 6D). Among the enriched pathways, thyroid hormone synthesis emerged as one of the most crucial pathways.

Tyrosine serves as a direct substrate for the synthesis of the thyroid hormones triiodothyronine (T3) and thyroxine (T4), while phenylalanine acts as a precursor that is converted to tyrosine by phenylalanine hydroxylase, thereby contributing to thyroid hormone biosynthesis. Based on the combined results of MetPA and MSEA analyses, thyroid hormone synthesis is therefore inferred to be a key metabolic pathway involved in COPD.

Amelioration of COPD metabolic dysregulation by BTHTT via the phenylalanine-tyrosine-thyroid hormone axis
Continued tissue metabolomics analysis showed that during BTHTT treatment, metabolic markers in the treatment groups exhibited a trend toward normalization. In the H group, 16 metabolic markers were significantly reversed (p < 0.05, p < 0.01), including the key markers glutathione (oxidized), phosphocholine, and L-palmitoylcarnitine (p < 0.05, p < 0.01), as shown in Figure 6E. In contrast, only seven metabolic markers were significantly reversed in the L group, with no significant reversal observed for the key markers (Supplemental File 1-–Supplemental Table S3). These results further confirmed the superior regulatory capacity of the high-dose BTHTT treatment group. The H group significantly regulated 12 metabolic pathways by MetPA analysis, especially phenylalanine, tyrosine, and tryptophan biosynthesis (impact > 0.5) (Figure 6F). Compared with the M group, MSEA analysis of the H group showed similar significant regulation of metabolic pathways, including thyroid hormone synthesis (Figure 6G), and it was particularly notable that thyroid hormone synthesis was a significant regulatory pathway shared by both MSEA analyses. In summary, based on the significant therapeutic effects of BTHTT (H group) on COPD, it markedly reversed disturbances in multiple biomarkers and metabolic pathways. Thyroid hormone synthesis could be one of the key metabolic pathways through which BTHTT exerts its therapeutic effects in COPD.

Amelioration of pathological and phenotypic manifestations in a COPD rat model by BTHTT
After 2 weeks of intervention with BTHTT (Figure 7), the body weight of rats in the M group increased slowly, whereas weight gain in the treatment groups, especially the high-dose group (H), was significant (p < 0.01), and the overall weight returned to the level of the control group. Spontaneous activity analysis showed a significant improvement in the number of movements in the treatment groups (p < 0.01). Serum and BALF tests indicated that BTHTT intervention significantly reduced the elevated levels of TNF-α, IL-6, IL-8, MCP-1, and OPN in the M group (p < 0.05, p < 0.01), and the levels of inflammatory factors in the H group were no longer significantly different from those in the C group. Histopathological examination revealed persistent lymphocyte infiltration around the bronchioles and thickening of the alveolar walls with granulocyte infiltration in the M group. The low-dose group (L) showed reduced inflammatory infiltration (local lymphocyte/neutrophil infiltration, slight alveolar thickening), and the H group exhibited near-complete repair of pathological damage (no significant inflammatory cell infiltration or alveolar structural abnormalities). In summary, BTHTT exerted a dose-dependent therapeutic effect on COPD, with the H group demonstrating the greatest capacity for pathological repair.

Amelioration of COPD metabolic dysregulation by BTHTT via the phenylalanine–tyrosine–thyroid hormone axis
Based on the aforementioned therapeutic outcomes, further transcriptomic analysis was conducted comparing the H group with the M group. Compared with the M group, a total of 224 differentially expressed genes were identified in the H group (Supplemental File 1-–Supplemental Table S4), including 70 downregulated genes and 154 upregulated genes. Among these, the H group regulated the key characteristic genes IL-1R2 and IL-1β (Figure 8A).

Further KEGG molecular functional annotation analysis of all regulated differentially expressed genes revealed that regulation in the BTHTT H group remained primarily focused on the cytokine–cytokine receptor interaction pathway, with regulation of IL-1R2/IL-1β continuing to play a dominant role (Figure 8B). In addition, regulation of multiple genes within the CC subfamily (CCR8, CCL4, CCL4L1, CCL4L2, CCL5) and the CXC subfamily (CXCL1, CXCL3, CXCL2, CXCL4, CXCL4L1, CXCR4) also showed significant changes, which may contribute to the therapeutic effects of BTHTT (Figure 8C).

Rebalancing of IL-1β/IL-1R2 and thyroid hormone synthesis in COPD by BTHTT
Thyroid hormone receptors are expressed across various pulmonary cell types, including alveolar epithelial cells, bronchial epithelial cells, and intrapulmonary immune cells. Consequently, the lung serves as a direct target for thyroid hormones, which can be locally delivered via bronchoalveolar lavage fluid (BALF).

Furthermore, lung tissue expresses type II deiodinase (D2), an enzyme responsible for intracellular conversion of circulating inactive thyroxine (T4) into its biologically active form, triiodothyronine (T3). Therefore, the levels of T3 and T4 within BALF do not merely reflect systemic infiltration from the blood but also indicate localized uptake, metabolism, and activation within lung tissue. During states of inflammation or injury, this localized metabolic profile may undergo significant alterations. Based on the above results, we further validated the changes in IL-1β, T3, and T4 in lung tissue using BALF (Figure 9A–C). In the detection of IL-1β, the M group maintained a high level, whereas all treatment groups showed significantly reduced IL-1β levels (p < 0.05, p < 0.01). The H group exhibited the most pronounced reduction, showing no significant difference compared with the C group. In the detection of T3 and T4, the M group showed significantly decreased levels, while all treatment groups led to substantial increases in T3 and T4 (p < 0.01), with no significant difference compared with the C group. Finally, immunohistochemical analysis of IL-1R2 in lung tissue was performed to observe protein expression changes across groups. Compared with the M group, the H group showed a higher expression level of IL-1R2-positive cells in rat lung tissue (Figure 9D).

The Spearman correlation coefficient method was used to confirm the pathway correlation between IL-1β/IL-1R2 and thyroid hormone synthesis (Table 2). IL-1β was significantly negatively correlated with T3 (ρ = −0.587, p = 0.002) and T4 (ρ = −0.622, p = 0.001), indicating that higher levels of inflammation were associated with lower levels of thyroid hormones. IL-1R2 showed a significant moderate positive correlation with T3 (ρ = 0.536, p = 0.006) and T4 (ρ = 0.567, p = 0.003), suggesting that increased receptor antagonism or regulation was associated with elevated thyroid hormone levels. A strong positive correlation was observed between T3 and T4 (ρ = 0.818, p < 0.001), consistent with their known physiological relationship.

These results further indicate that the transcriptional changes in IL-1β and IL-1R2 were confirmed, indicating abnormal gene expression. Abnormal levels of T3 and T4 also suggested dysregulation of thyroid hormone synthesis identified in metabolomics analysis. BTHTT treatment significantly altered these abnormal changes. Overall, the results suggest a potential negative correlation between the cytokine–cytokine receptor interaction pathway (mainly IL-1β/IL-1R2) and thyroid hormone synthesis (mainly T3/T4) in COPD. In summary, we propose a speculative mechanistic model involving the IL-1β–IL-1R2–thyroid hormone synthesis axis, which may be modulated by BTHTT.

Integrated UHPLC-MS and chemometrics for characterization of BTHTT's in vivo/in vitro components
Ultra-high-performance liquid chromatography coupled with a speculative mechanistic model here involving the IL-1β-–IL-1R2-–thyroid hormone synthesis axis, and which may be modulated by BTHTT, coefficients greater than 0.9, indicating stable and reliable data (Supplemental File 1-–Supplemental Figure S4). The base peak chromatograms (BPC) in both positive and negative ion modes are shown in Figure 10A,B. The serum samples from the H group and the in vitro test samples exhibited significant differences in the chromatograms. Moreover, clear distinctions were observed between these groups and the blank group, as well as the blank serum plus in vitro test samples.

The acquired data, including mass, isotope distribution, and MS/MS fragmentation information, were compared with the commercial standard Traditional Chinese Medicine (TCM) database. The results were further matched with public databases, such as GNPS²⁸, ReSpect²⁹, and Massbank³⁰ for compound identification and annotation. A total of 2,547 in vitro components of BTHTT were identified (1,649 in positive ion mode and 968 in negative ion mode) using a mass error threshold of < 25 ppm for MS1 and a matching score > 0.7 for MS2. NPClassifier analysis indicated that the predominant in vitro components were alkaloids (24%) and shikimate/phenylpropanoid derivatives (24%) (Figure 10C). Further analysis of BPC led to the selection of 38 high-abundance in vitro components of BTHTT (Figure 10D,E and Supplemental File 1—Supplemental Table S5), with shikimate/phenylpropanoid derivatives accounting for 51% (Figure 10F).

Based on the analysis and identification results of blank control serum and blank serum plus in vitro test samples of BTHTT, combined with the in vitro full component analysis results, a background subtraction algorithm in the chemometric module was applied to the H group to ultimately determine the in vivo components of BTHTT. A total of 283 in vivo components were identified (153 in positive ion mode and 132 in negative ion mode). NPClassifier chemical classification revealed that the predominant in vivo components were shikimate/phenylpropanoid derivatives (39%), alkaloids (20%), and terpenoids (16%) (Figure 10G).

In vivo components were further cross-referenced with public databases such as PubChem and ChemSpider. Using a mass error threshold of ppm ≤ ±5 and a matching score > 0.9 for MS2, 71 in vivo components of BTHTT were ultimately precisely identified (Supplemental File 1–Supplemental Table S6). Among them, seven high-abundance in vitro components that were also detected in vivo were identified: Vitamin B1, magnoflorine, cryptotanshinone, tanshinone IIA, 3,4-dihydroxyphenylacetic acid, salvianolic acid A, and gibberellin A4.

Prediction and analysis of potential therapeutic substance basis
To explore the structural basis of the therapeutic effects of BTHTT, we performed molecular docking of its in vivo components against the core protein target IL-1R2. Using molecular docking software via an AI supercomputing platform, we assessed the binding orientations of these compounds.

Based on a binding energy threshold of ≤-8 kcal/mol, 12 compounds were identified as having high structural compatibility with the IL-1R2 binding pocket (Supplemental File 1–Supplemental Table S7). Notably, Tanshinone IIA and Cryptotanshinone exhibited the lowest predicted binding energies. Furthermore, these two compounds are significant as they represent the primary high-content bioactive components identified in in vitro analysis, suggesting they are plausible candidates for further investigation.

To further evaluate the supportive computational evidence for these interactions, we conducted molecular dynamics (MD) simulations to assess the stability of the Tanshinone IIA–IL-1R2 and Cryptotanshinone–IL-1R2 complexes under simulated physiological conditions. As shown in Figure 11A,B, both systems reached their minimum potential energy within the first 300 ps. Following NVT equilibration, the systems maintained a stable temperature of approximately 300 K. NPT equilibration successfully stabilized the pressure at approximately 1 bar. While minor fluctuations were observed, the consistent system density confirmed effective pressure control. The Root Mean Square Deviation (RMSD) for both complexes reached a stable plateau after 17 ns, with fluctuations remaining within a narrow range of 0.2 nm. Throughout the 20 ns production run, the total energy of both systems remained steady with minimal variance.

In summary, the MD simulations indicate that Tanshinone IIA and Cryptotanshinone maintain a stable structural association with IL-1R2. These findings provide computational support for their potential role as active constituents, though subsequent in vitro and in vivo functional assays are required to establish their biological efficacy.

Despite the significant findings of this study, several limitations must be acknowledged. First, the relationship identified between inflammatory signaling and thyroid hormone levels is based on associative omics data. Further functional perturbation experiments are required to establish a definitive causal link.

Specifically, future studies should directly evaluate the hypothalamic-pituitary-thyroid (HPT) axis, alongside histological and functional assessments of the thyroid gland itself. Additionally, the identification of active compounds relied heavily on molecular docking and bioinformatic inference; therefore, targeted in vitro validation of specific compounds is necessary. While the animal model utilized simulates human Chronic Obstructive Pulmonary Disease (COPD), interspecies physiological differences may limit the direct clinical translation of these findings.

In summary, this study integrates multi-omics data and computational modeling to generate mechanistic hypotheses regarding the role of the IL-1R2/IL-1β axis and thyroid hormone metabolism in COPD. Although our results demonstrate the therapeutic potential of BTHTT and provide a paradigm for modern Traditional Chinese Medicine (TCM) research, further experimental validation is essential to confirm the proposed mechanisms.

Data Availability Statement:
Publicly available gene expression datasets used in this study were obtained from GEO under the following accession numbers: GSE11784, GSE12472, GSE16972, GSE38974, and GSE222965. Additionally, the original high-throughput sequencing data generated during the experimental phase of this study have been deposited in the NCBI BioProject database under accession number PRJNA1286104, accessible via the following link:https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1286104. All datasets are publicly available and meet the journal's data sharing requirements.

Figure 1: Machine-learning-based identification of COPD diagnostic biomarkers. (A) Volcano plot of DEGs. (B-G) Feature selection via (B,C) LASSO regression, (D,E) SVM-RFE, and (F,G) Random Forest. (H) Venn diagram showing 10 consensus genes. (I) Expression heatmap of key signatures. (J-L) ROC curve and expression analysis for (J) training and (K,L) validation cohorts. *p < 0.05, **p < 0.01, ***p < 0.001. Please click here to view a larger version of this figure.

Figure 2: Network pharmacology analysis. (A) Compound-target interactions. (B) PPI network. (C) Topological analysis of the PPI network. Please click here to view a larger version of this figure.

Figure 3: Validation of the COPD rat model. (A,B) Body weight and spontaneous activity (pre- vs. post-modeling). (C) Pro-inflammatory factors in serum and BALF. (D) H&E-stained lung sections (200×). Control shows normal architecture; Model shows (i) septal thickening, (ii) bronchial epithelial degeneration, (iii) mucin hypersecretion, and (iv) cyst-like structures. Scale bars = 50 µm. *p < 0.05, **p < 0.01. Please click here to view a larger version of this figure.

Figure 4: Transcriptomic profiling of COPD versus Control rats. (A) Volcano plot highlighting 39 downregulated and 134 upregulated genes. (B) KEGG pathway enrichment. (C) PPI network of DEGs identifying IL-1β as a central hub. (D) Cytokine-cytokine receptor interaction pathway (red: upregulated; blue: downregulated). Please click here to view a larger version of this figure.

Figure 5: Metabolic profiling of lung tissue in COPD rats. (A) PCA, (B) OPLS-DA, and (C,D) permutation tests in negative and positive ion modes. (E,F) Heatmaps of differential metabolites. Please click here to view a larger version of this figure.

Figure 6: Metabolic markers and pathway determination. (A,B) Correlation analysis in negative and positive modes. (C,D) MetPA and MSEA analysis identifying key pathways (e.g., phenylalanine, tyrosine, and tryptophan biosynthesis). (E) Heatmap of 29 metabolic markers following BTHTT treatment. (F,G) MetPA and MSEA analysis of the High-dose BTHTT group. Please click here to view a larger version of this figure.

Figure 7: Therapeutic effects of BTHTT on COPD rats. (A,B) Dose-dependent changes in body weight and spontaneous activity during treatment. (C) Pro-inflammatory factors in serum and BALF. (D) Representative H&E lung histopathology (200×). High-dose group shows near-normal architecture; Low-dose group shows localized infiltration. Scale bars = 50 µm. *p < 0.05, **p < 0.01, #p < 0.05 (vs. Control). Please click here to view a larger version of this figure.

Figure 8: Transcriptomic profiling of High-dose BTHTT versus Model group. (A) Expression levels of IL-1R2 and IL-1β. (B) KEGG enrichment. (C) Cytokine-cytokine receptor interaction pathway (red: upregulated; blue: downregulated). Please click here to view a larger version of this figure.

Figure 9: The changes in levels of IL-1β, T3, T4, and IL-1R2 in the lung tissue of rats. (A–C) Concentrations of IL-1β, T3, and T4 in BALF. (D) Immunohistochemical staining of IL-1R2 in lung tissue. *p < 0.05, **p < 0.01. Please click here to view a larger version of this figure.

Figure 10: Chemical characterization of BTHTT's in vivo/in vitro components. (A,B) BPC chromatograms of BTHTT in vitro and in vivo with positive and negative ion modes. (C) Npclassifier distribution of main chemical categories. (D,E) BPC of high-content components, peak numbers 15 and 16 correspond to Cryptotanshinone and Tanshinone IIA, respectively. (F) Classification of high-content components. (G) Classification of in vivo transitional components. Please click here to view a larger version of this figure.

Figure 11: Molecular dynamics simulations of tanshinone IIA and cryptotanshinone with IL-1R2. (A) Tanshinone IIA—IL-1R2, (B) Cryptotanshinone—IL-1R2. Please click here to view a larger version of this figure.

Figure 12: Schematic diagram of the mechanism of BTHTT-mediated thyroid hormone synthesis pathway in the treatment of COPD through IL-1R2/IL-1 β and phenylalanine, tyrosine, and trypsin biosynthesis. Please click here to view a larger version of this figure.

Time (min)	Mobile Phase A (%)	Mobile Phase B (%)
Initial	95	5
3	75	25
8.5	55	45
14	5	95
17	2	98
17.2	95	5
20	95	5

Table 1: Gradient elution method for in vivo and in vitro component analysis of BTHTT.

			IL-1β	IL-1R2	T3	T4
Spearman Rho	IL-1β	correlation coefficient	1	-0.176	-.587**	-.622**
		Significance	.	0.401	0.002	0.001
		(dual tailed)
		N	25	25	25	25
	IL-1R2	correlation coefficient	-0.176	1	.536**	.567**
		Significance	0.401	.	0.006	0.003
		(dual tailed)
		N	25	25	25	25
	T3	correlation coefficient	-.587**	.536**	1	.818**
		Significance	0.002	0.006	.	0
		(dual tailed)
		N	25	25	25	25
	T4	correlation coefficient	-.622**	.567**	.818**	1
		Significance	0.001	0.003	0	.
		(dual tailed)
		N	25	25	25	25

Table 2: Correlation analysis based on Spearman coefficient.
NOTE: **. At the 0.01 level (double tailed), the correlation is significant.

Supplemental File 1. Supplementary figures and tables providing additional quality control analyses, differential expression results, metabolomic profiling, and chemical characterization supporting the study.Please click here to download this file.

Discussion

The synthesis and maintenance of thyroid hormone levels are significant for COPD patients and extend far beyond mere endocrine function, profoundly influencing disease progression, nutritional status, pulmonary function, and overall prognosis^31,32,33. Particularly during acute exacerbations, COPD patients often exhibited a hypothyroidism-like state referred to as “low T3 syndrome”. It was reported that decreased thyroid hormone levels were positively correlated with the severity of COPD^34,35. First, T3 was a key hormone regulating basal metabolic rate and protein synthesis. A decline in T3 levels led to reduced muscle protein synthesis and enhanced catabolism, which was one of the important mechanisms behind muscle atrophy, weight loss, and malnutrition in COPD patients^36,37. Second, thyroid hormones were crucial for maintaining respiratory muscle function and driving respiration; insufficiency can cause respiratory muscle weakness, easy fatigue, and worsened dyspnea³⁸. Finally, thyroid hormones also affected lung tissue repair, surfactant production, and bronchial tone regulation, and their deficiency promoted the deterioration of pulmonary function³⁹.

COPD was a systemic inflammatory disease in which inflammatory factors interfere with the synthesis, release, or peripheral conversion of thyroid hormones, specifically the conversion of T4 to the more biologically active T3. IL-1β, as a critically important pro-inflammatory cytokine, not only directly inhibited thyroid hormone secretion and disrupted the hypothalamic-pituitary-thyroid (HPT) axis, but also impaired thyroid immune tolerance. This promoted the production of autoantibodies and the infiltration of inflammatory cells, leading to the destruction of thyroid follicles and impaired hormone synthesis^40,41. Furthermore, it was important to note that the key metabolic pathways of phenylalanine, tyrosine, and tryptophan biosynthesis were also closely related to thyroid hormone synthesis. Tyrosine was the direct precursor for thyroid hormone synthesis, and phenylalanine served as the precursor for tyrosine, being converted directly to tyrosine in the liver by the action of phenylalanine hydroxylase^42,43.

Building upon these preliminary findings from the literature and integrating them with the results of the current study, we propose a novel mechanistic hypothesis that BTHTT may exert therapeutic effects via the "IL-1R2/IL-1β–thyroid hormone synthesis pathway." In this hypothesized model, IL-1β acts as a critical pro-inflammatory cytokine that may disrupt thyroid hormone synthesis by interfering with the hypothalamic-pituitary-thyroid (HPT) axis or impairing thyroid immune tolerance. Our findings speculated that IL-1R2 functions as a "molecular lock" or decoy receptor for IL-1β, precisely inhibiting its pro-inflammatory activity. We hypothesize that BTHTT utilizes this signaling role to negatively regulate IL-1β, thereby potentially reducing the downstream impact on thyroid hormone synthesis.

In addition, through multi-algorithm collaborative mining and molecular docking, we also identified a potential pharmacodynamic base including Apigenin 7-glucuronide, Tanshinone IIA, and Cryptotanshinone, among others. Notably, Tanshinone IIA and Cryptotanshinone were identified as high-content blood-entry components. While molecular dynamics simulations showed excellent binding affinity within the hypothesized IL-1R2/IL-1β pathway, these results currently represent computational predictions that require further biological validation.

Despite the multiple findings reported in this study, several limitations must be acknowledged. First, the identified relationship between inflammatory signaling and thyroid hormone levels is based on associative multi-omics data. Consequently, further functional perturbation experiments are required to establish a definitive causal link. Future investigations should include a direct evaluation of the hypothalamic-pituitary-thyroid (HPT) axis, complemented by histological and functional assessments of the thyroid gland.

Furthermore, the identification of active compounds relied heavily on molecular docking and bioinformatic inference; these findings necessitate further in vitro validation targeting specific compounds. While the animal model employed successfully simulates human Chronic Obstructive Pulmonary Disease (COPD), inherent physiological differences between species may limit the direct clinical translation of these results. In summary, this study first integrated multi-omics data and computational modeling to propose a mechanistic hypothesis regarding the role of the IL-1R2/IL-1β axis and thyroid hormone metabolism in COPD (Figure 12). Second, our results demonstrate the therapeutic potential of BTHTT and provide a methodological paradigm for modern Traditional Chinese Medicine (TCM) research. Finally, this work establishes a foundation for subsequent experimental validation to confirm the proposed novel therapeutic mechanisms.

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Acetonitrile (chromatographic-grade)	Merck, Germany	1499230-935
Agilent 4150 bioanalyzer	Agilent Technologies, USA
Ammonium acetate	SIGMA, Germany	SRE0035-20230321
Angelica sinensis (Oliv.) Diels (Angelicae Sinensis Radix),	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240807
Astragalus membranaceus (Fisch.) Bge. (Astragali Radix),	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240601
ChemOffice		chemical structure drawing and molecular modeling software
Cinnamomum cassia (L.) J. Presl (Cinnamomi Cortex),	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240901
Eppendorf Centrifuge 5430 R	Eppendorf, Germany
Experimental cigarettes	Harbin Taihua, China	20230208
GraphPad Prism 8		statistical analysis and graphing software
Huangshan cigarettes		research-grade cigarettes for smoke exposure experiments
IL-1β ELISA kit	Shanghai Enzyme-Linked Biotechnology, China	E-EL-R0012c-20230522
Illumina Novaseq 6000 sequencing platform	Illumina, USA
Lepidium apetalum Willd. (Descurainiae Semen)	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240201
Lepidium sativum L. (Raphani Semen)	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240301
Lipopolysaccharide (LPS)	Shandong Huino Pharmaceutical, China	20230315
Mahonia fortunei (Lindl.) Fedde (Folium Mahoniae)	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240301
Methanol (chromatographic-grade)	Fisher, USA	DK745764
Nanodrop ND-2000 spectrophotometer	Thermo Scientific, USA
Perilla frutescens (L.) Britt. (Perillae Folium)	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240510
Pseudostellaria heterophylla (Miq.) Pax (Pseudostellariae Radix),	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20250201
Qubit 4.0 fluorometer	Thermo Fisher Scientific, USA
RNA library construction kit	ABclonal, USA	RK20301-20230415
SMINA		molecular docking software
Salvia miltiorrhiza Bge. (Salviae Miltiorrhizae Radix et Rhizoma)	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240902
Small animal spontaneous activity recorder (KW-ZF)	Nanjing Karlvin Biotechnology, China
Sodium pentobarbital	Shanghai Experimental Reagent Procurement Company, China	20230412
TNF-α ELISA kit	Shanghai Enzyme-Linked Biotechnology, China	E-EL-R2856c-20230522
Traditional Chinese Medicine MS database (Shanghai Applied Protein Technology Co., Ltd.)		commercial traditional Chinese medicine mass spectrometry spectral database
Trizol reagent	Magen, China	20230508
UPLC BEH Amide column		amide-based UHPLC analytical column
UPLC HSS T3 column		reversed-phase UHPLC analytical column (T3-type)
Vanquish Neo UHPLC-Orbitrap Exploris 480 MS system	Thermo Scientific, USA
Vanquish UHPLC-Q-Exactive HF-X MS system	Thermo Scientific, USA