Multi-Dimensional Characterization of Programmed Cell Death Patterns for Prognostic Stratification and Therapeutic Insights in Sepsis

Introduction

Sepsis, a life-threatening organ dysfunction caused by a dysregulated host response to infection, remains a pressing global health challenge, despite clinical advances.¹ Sepsis is associated with significant morbidity and mortality, with global report in 2020 estimates indicating approximately 48.9 million cases and 11 million deaths annually, accounting for nearly 20% of the overall deaths.² This varies considerably across regions, from 20–30% in developed countries to as high as 60% in resource-limited regions.³ From 1998 to 2023, United Kingdom sepsis-coded hospital admissions rose from 27.9 to 210.4 per 100,000 population, up to a 7.5-fold increase,⁴ and a total of 4,177,071 deaths were attributed to sepsis-related complications from 1999 to 2022 in the United Status,⁵ highlighting the growing burden of sepsis. The pathophysiology of sepsis is complex and involves excessive inflammation, immune suppression, coagulation abnormalities, and metabolic dysfunction, ultimately leading to multi-organ failure.⁶ This uncontrolled systemic inflammatory response, often accompanied by cytokine storms and endothelial damage, drives progression to septic shock and multiple organ dysfunction syndrome (MODS).⁷ Although early interventions, such as fluid resuscitation, antibiotics, and organ support have improved outcomes, mortality remains unacceptably high, especially in critically ill patients.^6,8 Conventional biomarkers, such as procalcitonin, C-reactive protein, and lactate, offer limited accuracy in predicting sepsis prognosis.⁹ Consequently, there is an urgent need to identify reliable diagnostic and prognostic biomarkers to guide personalized treatment strategies and improve outcomes in patients with sepsis.

Programmed cell death (PCD) is a genetically regulated process essential for maintaining cellular and internal environmental homeostasis,¹⁰ and could be classified into inflammatory and non-inflammatory forms according to their immunological consequences. Inflammatory PCD, such as pyroptosis,¹¹ necroptosis,¹² netotic cell death,¹³ immunogenic cell death (ICD),¹⁴ alkaliptosis,¹⁵ is typically characterized by the release of damage-associated molecular patterns (DAMPs), cytokines, or extracellular traps, thereby amplifying systemic inflammation and tissue injury in sepsis. By contrast, non-inflammatory PCD generally proceeds in a more immunologically silent and context-dependent manner, including apoptosis,¹⁶ ferroptosis,¹⁷ autophagy-dependent death,¹⁸ lysosome-dependent cell death,¹⁹ parthanatos,²⁰ mitochondrial permeability transition (MPT)-driven necrosis,²¹ entotic death,¹³ oxeiptosis,²² cuproptosis,²³ etc. Increasing evidence highlights extensive crosstalk and compensatory mechanisms among PCD pathways.²⁴ Compared with analyses focused on a single PCD type, prognostic models that integrate multiple PCD patterns provide a more comprehensive representation of disease pathophysiology and demonstrate improved performance in predicting patient outcomes and drug responsiveness.²⁵ Previous studies mainly examined the sepsis prognostic relevance of the individual PCD type,^26–28 however, a comprehensive understanding of how multiple PCD patterns collectively influence disease progression and patient outcomes remains lacking. These observations underscore the need for multi-targeted therapeutic strategies and integrated PCD-based prognostic models for effective sepsis management.

Traditional Chinese Medicine (TCM) has garnered increasing attention as a complementary approach to sepsis management because of its multitarget and multipathway therapeutic potential.²⁹ Simiao Yongan Decoction (SMYAD), a classic TCM formula composed of Lonicerae japonicae flos, Figwort root, Angelicae sinensis radix, and Licorice, is traditionally used to treat heat-toxin syndromes and suppurative infections.³⁰ Recent studies have shown that SMYAD possesses anti-inflammatory, immunomodulatory, antioxidative, and microcirculatory protective effects,³¹ which closely align with the pathophysiological mechanisms of sepsis. Experimental evidence suggests that SMYAD can attenuate systemic inflammation, inhibit excessive cytokine release, and protect against organ dysfunction,^32,33 indicating its promising potential as an adjunctive therapy in sepsis treatment. However, its mechanistic relevance to PCD in sepsis remains unclear.

In this study, we aimed to elucidate how different PCD patterns collectively shape prognosis and therapeutic response in sepsis. We conducted the first multi-dimensional investigation of 14 PCD patterns in sepsis, identified molecular subtypes with distinct immune and clinical characteristics, developed and validated a robust seven-gene prognostic model, and explored the interaction between SMYAD-derived compounds and PCD-related proteins. Collectively, this framework deepens the mechanistic understanding of PCD crosstalk in sepsis, and provides a translational basis for combining biomarker-guided precision medicine with traditional Chinese medicine approaches.

Materials and Methods

Data Source and Processing

Fourteen PCD pattern genes were obtained from a previous study,³⁴ including apoptosis, pyroptosis, necroptosis, ferroptosis, mitotic cell death, oxeiptosis, MPT-driven necrosis, parthanatos, cuproptosis, immunogenic cell death, alkaliptosis, lysosome-dependent cell death, entotic cell death, and autophagy. The microarray datasets for sepsis were downloaded from the Gene Expression Omnibus (GEO) database. GSE65682 included 802 healthy control (n = 42) and sepsis (n = 760) samples. In addition, other datasets were used to further verify the results: GSE48080 contains microarray expression patterns of peripheral blood mononuclear cells (PBMCs) from healthy volunteer (n = 3), survivor (n = 5) and non-survivor septic (n = 5) patients;³⁵ GSE28750 is a microarray cohort of PBMCs from healthy control (n = 20), sepsis (n = 10) and post-surgical (n = 11) patients.³⁶ GSE57065 contains patients based on the Simplified Acute Physiology Score II (SAPSII) score: SAPSII-high (n = 40) and SAPSII-low (n = 42) groups;³⁷ GSE26440 from normal control (n = 32), survivor (n = 81), and non-survivor septic shock (n = 17) patients;³⁸ GSE106878 comprised 84 septic shock patients treated with either Placebo or Hydrocortisone in the morning (Pre) and 24 h (Post) with different survivor status; GSE33119 of survival outcome from 20 patients before and after activated protein C (Xigis) treatment.

Correlation Analysis

Correlations among different PCD patterns and between immune cell infiltration and risk score were calculated and visualized using the corrplot package (version 0.5.9) in R. Significant correlations were defined as those with p < 0.05.

Survival Analysis

Univariate Cox proportional hazard regression analysis was performed to evaluate the prognostic significance of each PCD pattern in the sepsis cohort. Survival analysis was performed by survival package (version 3.8–3). Kaplan-Meier survival curves were generated to visualize the survival probability differences between groups, and the Log rank test was applied to assess the statistical significance of these differences. Hazard ratios (HR) and 95% confidence intervals (CI) were calculated to identify the PCD patterns associated with patient outcomes.

Consensus Clustering

Based on the activity scores of the PCD patterns with prognostic significance, we conducted consensus clustering using the ConsensusClusterPlus R package (version 1.70.0). The optimal number of clusters was determined using cumulative distribution function (CDF) curves and delta-area plots via Kvec and calcICL function. A principal component analysis (PCA) was performed by prcomp function to validate the separation of the identified clusters.

Gene Set Enrichment Analysis

Gene Set Enrichment Analysis (GSEA) was applied using the clusterProfiler package (version 4.14.6) in R. Hallmark gene sets from the Molecular Signatures Database (MSigDB) were used as the reference pathways. Statistical significance was set at a false discovery rate (FDR) < 0.05.

Differentially Expressed Genes Identification

Differentially expressed genes (DEGs) between the clusters were identified using the limma package (version 3.58.1). Gene expression values were first normalized and fitted to a linear model using lmFit function. Empirical Bayes moderation was applied via eBayes function to improve variance estimates across genes, and differential expression was assessed using moderated t-tests. P values were adjusted for multiple testing using the Benjamini-Hochberg false discovery rate (FDR < 0.05).

LASSO Regression

Prognosis-associated genes were first screened using univariate Cox regression analysis to identify genes significantly correlated with patient survival. Genes with p < 0.05 were retained for further analysis. Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression was conducted using the glmnet package (version 4.1–8) to construct a prognostic gene signature. LASSO applies an L1 penalty to the regression coefficients, shrinking some coefficients to zero, which effectively performs variable selection and reduces model overfitting. A ten-fold cross-validation was performed via cv.glmnet function to determine the optimal lambda value according to lambda.min parameter. Genes with non-zero coefficients at this lambda were retained and used to construct the final prognostic model. The coxph function from the survival package was then applied to fit a multivariate Cox proportional hazards model based on these selected genes, generating hazard ratios and coefficients for each gene. The risk score for each patient was calculated as a linear combination of the expression levels of the selected genes weighted by their Cox regression coefficients:

Patients were stratified into high- and low-risk groups based on the median risk score.

Receiver Operating Characteristic Curve

Time-dependent predictive performance of the prognostic model was evaluated using the timeROC R package (version 0.4). The area under the receiver operating characteristic (ROC) curve (AUC) was calculated at multiple time points (4, 7, 14, and 21 days) to assess the model’s predictive accuracy over time, providing a comprehensive assessment of its prognostic value.

Molecular Docking

The structures of prognostic target proteins were obtained from the Protein Data Bank (PDB), including MPO (PDB: 1CXP), ELANE (PDB: 5ABW), and IL1B (PDB: 1L2H). The active compounds of SMYAD were retrieved from the Traditional Chinese Medicine Systems Pharmacology (TCMSP) database, then obtained their structures from PubChem and ChemSpider database, including quercetin (PubChem CID: 5280343), oleic acid (PubChem CID: 445639), ursolic acid (PubChem CID: 64945), rutin (PubChem CID: 5280805), thymol (PubChem CID: 6989), and glycyrrhizin (ChemSpider ID: 14263). Protein structures were prepared using PyMOL (version 3.1) to remove water molecules and ligands. Molecular docking was performed using AutoDock (version 1.5.6). The docking grid was centered on the active site of each protein, with grid box dimensions set to cover the binding pocket. The Lamarckian Genetic Algorithm was applied with default parameters, and multiple docking runs were performed to ensure sampling of diverse binding poses. Docking results were ranked based on binding free energy, and the top binding conformations with the lowest energies were selected for further analysis. PyMOL was used to visualize binding poses, highlighting hydrogen bonds, hydrophobic interactions, and key residues involved in ligand recognition.

Single-Cell RNA Sequencing Analysis

Single-cell RNA sequencing (scRNA-seq) datasets (GSE167363³⁹ and GSE175453⁴⁰) were integrated via harmony method using the Seurat package (version 4.0.3) in R. Cells were remained with (1) detected genes in 200~6000; (2) read counts > 1000; (3) mitochondrial gene content < 20%. Highly variable genes were identified using SCTransform function, and PCA was performed on the scaled data, then cells were clustered using the FindClusters function. Uniform Manifold Approximation and Projection (UMAP) was used for visualization of cell clusters in two dimensions. Cells were annotated into major immune cell types based on the expression of canonical markers. PCD activity scores were estimated using AddModuleScore function for each cell type.

Multi-Organ Temporal RNA-Seq Data Analysis

Multi-organ and multi-time point RNA-seq data from a mouse sepsis model (GSE224127⁴¹) were analyzed to evaluate the spatiotemporal expression dynamics of prognostic genes. Tissues included bone marrow, peripheral blood, brain, heart, kidney, liver, skin, colon, small intestine, spleen, thymus, lymph nodes, and lungs, collected at 0, 0.25, 0.5, 1, 2, 3, and 5 d after lipopolysaccharide (LPS) injection. The mean gene expression was visualized to determine the temporal expression peaks and tissue-specific patterns. The spatial transcriptome (GSE248904⁴¹) across different organs in the control and LPS groups was subjected to the Seurat object.

Statistical Analysis

Data processing, statistical analyses and visualization were conducted using R software (version 4.1.2). Survival curve analyses (Kaplan-Meier curves) were performed via the Log rank test. Hazard ratios in the were calculated using univariate cox regression analyses. Correlations were calculated using the Pearson method. P value < 0.05 was considered statistically significant (* p < 0.05; ** p < 0.01; *** p < 0.001; **** p < 0.0001).

Results

Identification of Programmed Cell Death Clusters in Sepsis

The workflow of this study was shown as Figure 1A. To systematically delineate the landscape of PCD in sepsis, we quantified the scores of 14 PCD patterns in the sepsis cohort (GSE65682), including 42 healthy controls and 760 sepsis patients. The PCD patterns were significantly different between the healthy and sepsis groups (Figure S1). Correlation analysis revealed substantial interactions among several PCD modalities. Notably, pyroptosis showed a strong positive correlation with necroptosis (r = 0.70), and necroptosis was strongly correlated with immunogenic cell death (r = 0.70) (Figure 1B), suggesting potential crosstalk among these cell death pathways within the septic microenvironment.

Figure 1 Continued.

Figure 1 Identification of programmed cell death clusters in sepsis. (A) The workflow of this study. (B) Correlation network of 14 programmed cell death (PCD) patterns in sepsis. (C) Prognostic analysis of each PCD pattern using Kaplan-Meier survival curves. (D) Cumulative distribution function (CDF) curve for consensus clustering. (E) The area under curve (AUC) for consensus clustering. (F) Consensus matrix heatmap at k = 2. (G) Principal component analysis (PCA) plot demonstrating transcriptional separation between the two clusters. (H) Kaplan-Meier survival of two clusters. (I) Kaplan-Meier survival of two clusters according to gender and age.

Based on the prognostic impact of each PCD pattern, we stratified the cohort into high and low groups and identified six PCD patterns with significant prognostic characteristics (p < 0.05, Figure 1C). Ferroptosis and mitotic cell death are associated with increased mortality risk, whereas lysosome-dependent cell death, MPT-driven necrosis, necroptosis, and pyroptosis are associated with improved survival outcomes. Therefore, these six PCD patterns were selected for subsequent analyses.

We next performed consensus clustering based on PCD pattern scores to identify distinct molecular subtypes of sepsis. Using the CDF and consensus matrix, we determined k = 2 was the optimal number of clusters (Figure 1D and E), resulting in the classification of sepsis patients into Cluster1 and Cluster2 with high intra-cluster average consistency (Figure 1F). PCA demonstrated a clear transcriptional separation between the two clusters (Figure 1G). Kaplan-Meier survival analysis revealed that patients in Cluster 2 had significantly better survival probabilities than those in Cluster 1 (Figure 1H), indicating that Cluster 1 was associated with a higher mortality risk. Importantly, sex and age distributions were comparable between the two clusters, thus minimizing potential demographic confounders (Figure 1I).

Characteristics of the Two Clusters in Sepsis

To elucidate the immunological characteristics of the identified clusters, we applied a single-sample gene set enrichment analysis (ssGSEA) to quantify immune cell infiltration across 28 immune-related gene sets. As shown in Figure 2A, cluster 1 exhibited markedly lower infiltration of immune cells, particularly Activated CD4 T cells, type 17 T helper cells, Activated B cells, CD56dim natural killer cells, activated/immature dendritic cells, and mast cells. In contrast, Cluster2 demonstrated significantly higher levels of immune infiltration, characterized by Activated CD8 T cells, effector memory CD8 T cells, T follicular helper cells, type 1 helper cells, Regulatory T cells, Immature B cells, myeloid-derived suppressor cells, natural killer T cells, macrophages, eosinophils, and neutrophils, suggesting a potential relationship between immune infiltration and patient outcomes. Furthermore, we observed substantial differences in the expression of immune checkpoint genes between the two clusters (Figure 2B). Cluster 2 generally exhibited higher immune checkpoint expression, indicating an active, yet regulated immune microenvironment.

Figure 2 Immunological and functional characteristics of the two PCD clusters in sepsis. (A) Immune cell infiltration scores quantified by single-sample gene set enrichment analysis (ssGSEA). (B) Differential expression of immune checkpoint genes between Cluster 1 and Cluster 2. (C) Gene set variation analysis (GSVA) results of the hallmark pathways differences between clusters. (D) Distribution of candidate prognostic genes across the two clusters. ns not significant, * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001.

To further investigate the biological underpinnings of the two clusters, we constructed gene set variation analysis (GSVA) to identify dysregulated hallmark pathways (Figure 2C). Cluster 1 was enriched in pathways related to hypoxia, epithelial-mesenchymal transition (EMT), glycolysis, and reactive oxygen species (ROS) signaling, suggesting a metabolic and stress-associated phenotype. In contrast, cluster 2 was predominantly enriched in immune-related pathways, including IL6-JAK-STAT3 signaling, inflammatory response, interferon-α/γ response, and complement activation, consistent with the heightened immune infiltration observed (Figure 2A). We then mapped the candidate prognostic genes⁴² to clusters (Figure 2D), in which cytokines, transcription regulators, and transmembrane receptors were mainly enriched in Cluster2 samples, and phosphatase CKS2, enzyme TGM2, APOBEC3B in Cluster1. These results indicate distinct molecular signatures and regulatory mechanisms between the two clusters.

Establishment and Validation of a PCD-Related Prognostic Model

To further refine the transcriptional differences between clusters, we identified DEGs (Figure 3A) and intersected them with prognosis-associated genes derived from univariate Cox regression analysis. Using LASSO Cox regression, we constructed a robust seven-gene prognostic signature comprising ELANE, CTSG, MPO, CAMP, TFRC, IL1B and CASP5 (Figure 3B and C). The risk score was calculated using the following formula: risk score = ELANE*0.36917+CTSG*0.22875+MPO*(−0.46844)+CAMP*(−0.31068)+TFRC*(0.28079)+IL1B*(0.10971)+CASP5*(−0.14144). These seven genes showed higher expression levels in the sepsis than healthy samples in GSE65682 cohort (Figure 3D). Stratification based on the median risk score revealed that patients in the high-risk group exhibited significantly worse survival outcomes than those in the low-risk group (Figure 3E), and the 6 PCD patterns were differentially expressed between the two groups (Figure S2). Moreover, the risk scores were significantly elevated in Cluster 1, which was also associated with poor prognosis (Figure 3F). Sex and age were evenly distributed between the high- and low-risk groups, further supporting the prognostic independence of this model (Figure 3G). Subgroup analysis confirmed that the survival probability differences remained consistent across the sex and age groups (Figure 3H).

Figure 3 Construction of a PCD-related prognostic model. (A) Volcano plot of differentially expressed genes (DEGs) between clusters. (B) Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression to identify key prognostic genes. (C) Risk score formula based on seven PCD-related genes. (D) Expression profiles of the seven prognostic genes in healthy and sepsis samples. (E) Kaplan-Meier survival curves showing worse survival in the high-risk group. (F) Distribution of risk scores in Cluster 1 and Cluster 2. (G) Gender and age distribution between high- and low-risk groups. (H) Subgroup survival analysis by gender and age. (I) The correlation of the seven genes with PCD patterns. (J) Time-dependent receiver operating characteristic (ROC) curves demonstrating the prognostic performance of the risk model at 4, 7, 14, and 21 days. ns not significant, *** p < 0.001, **** p < 0.0001.

Functionally, these seven prognostic genes were closely associated with ferroptosis, necroptosis, pyroptosis, mitotic cell death, and lysosome-dependent cell death (Figure 3I), underscoring their relevance in the cell death regulatory network of sepsis. Time-dependent ROC curve analysis demonstrated good predictive performance of the model, with AUC values of 0.768, 0.747, 0.714, and 0.682 at 4, 7, 14, and 21 days, respectively (Figure 3J).

We validated the prognostic robustness of the model across independent external genes (Figure 4A–G). The risk scores were higher in patients with sepsis than in healthy volunteers (Figure 4A), and patients with higher-grade sepsis exhibited higher risk scores than those with lower-grade sepsis (Figure 4C). Survivor patients consistently displayed lower risk scores than non-survivors (Figure 4A, D, and G). Additionally, surgical interventions were associated with reduced risk scores (Figure 4B, E, F), whereas xigris treatment (activated protein C) did not significantly alter the risk scores.

Figure 4 External validation of the PCD-based prognostic model. (A–G) Validation of the risk model across multiple independent GEO datasets. (H) Correlation between the prognostic risk score and immune cell infiltration. (I) Correlation heatmap of individual prognostic genes with various immune cell types. (J) Comparison of immune cell abundance between high- and low-risk groups. ns not significant, * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001.

Further correlation analysis revealed a significant positive association between the risk score and infiltration of activated CD4+ T cells, whereas neutrophils exhibited a strong negative correlation (Figure 4H). At the individual gene level, strong correlations were observed with natural killer T cells, myeloid-derived suppressor cells, monocytes, and neutrophils (Figure 4I). Immune cell abundance across the risk groups further highlighted the immunological divergence associated with the prognostic model (Figure 4J). In summary, the established seven-gene prognostic model provides a reliable survival prediction and offers insights into the immune landscape and mechanisms of cell death in sepsis.

scRNA-Seq and Multi-Organ Temporal Transcriptomics Reveal Cellular Heterogeneity and Dynamic Changes in Sepsis

To further investigate the cellular heterogeneity of sepsis, we integrated scRNA-seq data from human PBMCs in datasets GSE167363 and GSE175453, and obtained 102,689 cells in total, which derived from 7 healthy control (HC, n = 39,204) and 14 sepsis (n = 63,485) samples. Based on the expression of canonical markers, the cells were divided into eight cell types (Figure 5A and B): B cells (CD79A+), dendritic cells (DC, CD1D+CD1C+), monocytes (CD14+VCAN+), neutrophils (CSF3R+MMP9+), plasma (JCHAIN+), platelets (PPBP+), red blood cells (RBC, HBA1+), and T/NK cells (IL7R+). Cellular proportions varied across different samples in the two datasets (Figure S3), in which monocytes and neutrophils mostly increased in abundance, whereas DC were reduced in sepsis samples, suggesting a shift in the immune cell composition during disease progression.

Figure 5 scRNA-seq and multi-organ temporal transcriptomics data analysis. (A) UMAP visualization of 8 cell types. (B) Violin plots show the relative expression levels of canonical markers. (C) Heatmap shows the activities of PCD patterns in each cell type. (D) The expression of seven prognosis genes in the main cell types. (E) The expression of seven prognosis genes in different groups. (F and G) The expression of seven prognosis genes in different tissues (F) and time points (G) after lipopolysaccharide (LPS) injection. (H) The spatial distribution of different organs of ST-seq dataset. (I) The cell type module score from scRNA-seq dataset in different organs of ST-seq dataset. (J) The expression of seven prognosis genes in different organs between control and LPS groups. Abbreviations: PMBC, peripheral blood mononuclear cell; MU, muscle; LU, lung; BF, brown fat; BM, bone marrow; LI, liver; SK, skin; SI, small intestine; BR, brain; CO, colon; KI, kidney; SM, stomach; HE, heart; SP, spleen; PN, pancreas; LN, lymph node; TH, thymus; DC, dendritic cell; RBC, red blood cell.

We then quantified the activity scores of the 14 PCD patterns in the major cell types and found that PCD was predominantly active in plasma, DC, monocytes, and neutrophils, suggesting their central role in driving immune dysfunction and organ damage in sepsis (Figure 5C). Seven prognostic genes showed elevated expression in sepsis samples, particularly in DC, monocytes, and neutrophils (Figure 5D). Notably, these genes exhibited the highest expression levels in neutrophils of sepsis samples (Figure 5E), emphasizing the pivotal role of neutrophil-mediated immune dysregulation and cell death in sepsis.

To extend our findings to tissue-level and temporal dynamics, we analyzed the multi-organ and multi-timepoint transcriptomic profiles of LPS-induced sepsis mice (GSE224127), covering bone marrow (BM), PBMC, brain (BR), heart (HE), kidney (KI), liver(LI), skin (SK), colon (CO), small intestine (SI), spleen (SP), thymus (TH), inguinal lymph node (iLN), lung (LU) at 0, 0.25, 0.5, 1, 2, 3, and 5 days after sublethal LPS injection. The expression of the prognostic genes was predominantly enriched in the BM and PBMC samples (Figure 5F), and their levels gradually increased over time (Figure 5G). Importantly, the transcriptional activity of these genes plateaued around day 5, consistent with previous studies reporting transcriptional stabilization in non-lymphoid tissues after sepsis onset,⁴¹ suggesting that the most active inflammatory and cell death responses occur within the initial 3–5 days.

In addition, we utilized the spatial transcriptome across different organs in control and LPS-treated mice (Figure 5H). The cell type abundance in each organ was estimated using scRNA-seq-derived markers (Figure 5I). For example, TH, LN, and SP showed more T/NK and B cells and SK, LU, BM, Pancreas (PN)-enriched RBC, followed by neutrophils. Similarly, prognostic genes were mainly expressed in the BM and were elevated in the LPS group (Figure 5J). Taken together, single-cell and spatial transcriptomic analyses highlighted the central role of bone marrow and neutrophil-mediated cell death in sepsis progression, providing insights into the spatiotemporal dynamics of immune dysregulation.

Exploration of Potential Therapeutic Targets

To investigate potential therapeutic avenues, we integrated prognostic genes with active components from the SMYAD, a traditional Chinese medicine widely used for its anti-inflammatory and detoxifying properties. Based on the Traditional Chinese Medicine Systems Pharmacology Database, six active compounds with high bioavailability (quercetin, oleic acid, glycyrrhizin, ursolic acid, rutin, and thymol) were identified as candidate molecules that target prognostic genes (Figure 6A and B). Specifically, MPO and IL1B shared targets across multiple herbal components of Lonicerae japonicae flos, Figwort root, and Licorice, while ELANE was targeted by thymol from Lonicerae japonicae flos.

Figure 6 Molecular docking of active compounds from Simiao Yongan Decoction with prognostic target proteins. (A) The network of predicted interactions between Simiao Yongan Decoction (SMYAD) active components and prognostic genes. (B) The structures of active components. (C) Docking results showing the binding energies of each compound-protein complex.

To further explore the binding affinities and potential mechanisms of SMYAD, we performed a molecular docking analysis of the target proteins (MPO, IL1B, ELANE) with the six hub compounds in Figure 6B. We prioritized the docked conformations with favorable hydrogen-bonding interactions, selecting those with the lowest binding energies, ranging from −8.47 kcal/mol (IL1B-Ursolic acid) to −3.36 kcal/mol (MPO-Oleic acid) (Figure 6C). Key amino acid residues were identified as potential binding sites for the small molecules. For example, ARG323 of MPO interacts with quercetin, ARG525 of MPO interacts with oleic acid, and ARG129 of ELANE interacts with thymol. In IL1B, MET20 binds to ursolic acid, and LEU134 interacts with rutin. More importantly, the docking energies of the latter four (Figure 6C) were less than −5 kcal/mol. These residues are likely to contribute to ligand-binding specificity and complex stability.

Discussions

Sepsis is a life-threatening disease that is characterized by a dysregulated host response to infection.⁴³ Currently, clinical assessment of sepsis prognosis mainly relies on scoring systems such as the Acute Physiology and Chronic Health Evaluation (APACHE) II⁴⁴ and Sequential Organ Failure Assessment (SOFA)⁴⁵ but fails to fully reflect clinical complexity or accurately predict disease severity. Previous studies have demonstrated that multiple forms of PCD, including ferroptosis,⁴⁶ cuproptosis,²³ and pyroptosis,²⁸ are involved in sepsis, exhibiting extensive interactions and crosstalk with one another. However, comprehensive investigation of the overall PCD landscape in sepsis is lacking.

In this study, we systematically characterized 14 distinct PCD patterns in sepsis, identified two distinct molecular subtypes, and developed a robust seven-gene prognostic model that was validated using external datasets. Using single-cell and multi-organ transcriptomic data, we emphasized the cellular heterogeneity and dynamic molecular changes associated with sepsis progression. Finally, we explored the potential therapeutic value of SMYAD through a molecular docking analysis.

The seven-gene signatures ELANE, CTSG, MPO, CAMP, TFRC, IL1B, and CASP5 were significantly associated with poor survival outcomes in patients with sepsis. These genes are not only prognostically relevant but also functionally involved in various PCD pathways that contribute to the pathophysiology of sepsis. ELANE, which encodes neutrophil elastase, is a serine protease released by activated neutrophils that plays a critical role in pyroptosis⁴⁷ and netotic cell death.⁴⁸ Its overexpression can lead to excessive neutrophil extracellular trap (NET) formation, resulting in uncontrolled inflammation, tissue damage, and microvascular thrombosis, which are hallmarks of sepsis-induced organ dysfunction.⁴⁹ CTSG is another neutrophil-derived protease implicated in lysosome-dependent cell death, and has been reported to exacerbate inflammation by promoting proteolytic tissue injury and cytokine release in sepsis.⁵⁰ MPO, encoding myeloperoxidase, is a key enzyme in neutrophils responsible for ROS production. MPO is a central regulator of mitotic cell death,⁵¹ and elevated MPO activity has been linked to endothelial damage and increased disease severity in septic patients.⁵² CAMP, which encodes the cathelicidin antimicrobial peptide, not only exhibits potent antimicrobial activity, but also contributes to netotic cell death by modulating neutrophil activation and degranulation,⁵³ further amplifying inflammatory responses during sepsis. TFRC is an essential regulator of cellular iron uptake, and its dysregulation contributes to ferroptosis,⁵⁴ an iron-dependent form of PCD driven by lipid peroxidation. Increased TFRC expression during sepsis may exacerbate oxidative stress and ferroptosis in parenchymal and immune cells, thereby promoting organ injury.⁵⁵ IL1B, a classical pro-inflammatory cytokine and canonical marker of pyroptosis and necroptosis,⁵⁶ is activated and released in response to inflammasome signaling. It directly drives systemic inflammation and vascular leakage in septic shock.⁵⁷ CASP5, a member of the inflammatory caspase family, functions as a critical effector of noncanonical pyroptosis. Upon sensing intracellular LPS, CASP5 activates gasdermin D and initiates pyroptotic membrane rupture, promoting inflammatory cell death and immune dysregulation.⁵⁸ Together, these gene sets embody the convergence of inflammatory amplification and programmed cell death, leading to poor sepsis outcomes.

scRNA-seq data revealed that DCs, monocytes, and neutrophils exhibited high susceptibility to the prognostic gene signature, with neutrophils showing the strongest response to sepsis. In patients with sepsis, the abundance of DCs is reduced via apoptosis and pyroptosis,^59,60 accompanied by impaired maturation and diminished antigen-presenting capacity, ultimately leading to suppression of innate immune responses. Monocyte counts and percentages are significantly decreased in the early stage of sepsis, accompanied by the upregulation of apoptosis- and pyroptosis-related genes.^61,62 Continuous stimulation impairs the ability of monocytes to produce proinflammatory cytokines while promoting anti-inflammatory cytokine expression,⁶³ leading to sepsis-associated immunosuppression. Neutrophil dysfunction during sepsis is particularly notable These cells exhibit impaired chemotaxis,⁶⁴ delayed apoptosis,⁵⁹ and excessive formation of neutrophil extracellular traps (NETs),⁶⁵ resulting in a hyperinflammatory and cytotoxic phenotype. This persistent activation and dysregulated cell death contribute to immunopathology, endothelial dysfunction, and ultimately multi-organ failure.⁶⁶ NETs are web-like extracellular structures decorated with antimicrobial components such as neutrophil elastase (NE, encoded by ELANE), myeloperoxidase (MPO) and histones.⁶⁷ Notably, both ELANE and MPO were present in our prognostic gene signature, underscoring their potential involvement in sepsis-associated immune dysregulation.

Furthermore, the multi-organ, multi-time point transcriptomic dataset revealed that the expression of prognostic genes was predominantly elevated in bone marrow and PBMCs. As a dynamic immune organ, the bone marrow undergoes profound dysregulation during sepsis, leading to significant alterations in hematopoiesis and immune cell function.⁶⁸ Notably, sepsis shifts myeloid progenitors toward emergency myelopoiesis,⁶⁹ resulting in aberrant production of neutrophils and monocytes. PBMCs, as key mediators of the systemic inflammatory response,⁷⁰ show high expression of these prognostic genes, suggesting their critical involvement in peripheral immune activation and immune imbalance during sepsis. Importantly, the gene expression patterns exhibited a progressive increase, peaking around day 3 following LPS-induced sepsis, aligning with transcriptional homeostasis in non-lymphoid tissues after sepsis onset,⁴¹ suggesting that the most active inflammatory and cell death responses occur within the initial days, highlighting a critical therapeutic window during the early stages of sepsis.⁷¹ Morever, this temporal expression pattern was consistent with the kinetics of neutrophil activation and NETs formation,⁷² further supporting the contribution of neutrophil-driven pathology in early sepsis.

In parallel, we explored the pharmacological relevance of SMYAD, a classical TCM formula used for conditions resembling heat, toxin, and blood stasis syndromes in TCM theory.^73,74 Linear-trap-LC/MS^N identified 113 chemical compounds in the extract and medicated serum of SMYAD,³⁰ in which vitexin could inhibit the release of pro-inflammatory cytokines (TNF-α, IL-1β) from macrophages to effectively reduce leukocyte migration in septic mice;⁷⁵ Ferulic acid modulates inflammatory responses, autophagy and apoptosis during sepsis via NF-κB and Nrf2 pathways, thereby alleviating sepsis-related organ injury.^76,77 Furthermore, compounds such as Lonicerae japonicae flos and Figwort root are rich in flavonoids and organic acids, which can markedly decrease TNF-α and IL-6 levels, to reduce cytokine storm and disease exacerbation.^78,79 Clinically, SMYAD combined with antibiotics has been reported to further reduce TNF-α and IL-8 levels compared with antibiotics alone;⁸⁰ and combination with Imipenem/Cilastatin restore cytokine balance (eg, IL-2, TNF-α, IL-4) and attenuate systemic inflammation in septic rats.³⁰ In our study, molecular docking revealed that key SMYAD-derived compounds exhibited strong binding affinities to MPO, IL1B, and ELANE, which are central regulators of PCD pathways in sepsis. These interactions suggest that SMYAD may exert anti-inflammatory and immunomodulatory effects by directly targeting PCD-related proteins,⁸¹ supporting its potential as an adjunctive therapy for sepsis.

Previous studies typically focused on individual PCD pathways in sepsis. For example, Liang et al developed a six-gene (GZMB, CHMP7, NLRP1, MYD88, ELANE, and AIM2) prognostic model related with pyroptosis;²⁸ ferroptosis-related siguratures included TFRC, DPP4 and PRKCA;⁴⁶ hub genes in apoptosis contained BCL2, FASLG, IRF9 and JAK3;⁸² cuproptosis identified a five-gene signature (PTPRJ, FCN1, SVIP, LILRB1 and LST1);⁸³ necroptosis highlighted BACH2, GATA3, LEF1, and BCL2;⁸⁴ and PANoptosis, referring to the interplay of pyroptosis (P), apoptosis (A), and necroptosis (N), involved CD14, GSDMD, IL1B, and FAS.⁸⁵ These studies provided valuable insights, whereas leaving the broader landscape of PCD crosstalk in sepsis unexplored. Notably, accumulating evidence suggests that distinct cell death programs are not isolated but can interact and influence one another,⁸⁶ and the comprehensive analysis promote to uncover the potential therapeutic targets.⁸⁷ Our work systematically integrates 14 distinct PCD patterns into a unified analytical framework, emphasizing the pivotal roles of pyroptosis, ferroptosis, necroptosis, lysosome-dependent cell death, and netotic cell death. Combining previous findings with the relationship between prognostic genes and PCD patterns, we propose that in sepsis, distinct PCD programs form an interconnected network. Activated CASP5 cleaves gasdermin D (GSDMD) to generate a functional N-terminal fragment capable of pore formation on the plasma membrane,⁸⁸ enabling IL1B release into the extracellular space. This process triggers pyroptosis,¹¹ as well as necroptotsis,⁸⁹ thereby amplifying inflammatory responses. The cleavage of GSDMD can also be mediated by ELANE.⁹⁰ During sepsis, excessive inflammatory cytokines upregulate transferrin receptor (TFRC) expression on the cell surface,⁹¹ promoting Fe uptake and triggering ferroptosis.⁵⁴ Loss of lysosomal membrane integrity leads to the release of enzymes such as CTSG, causing lysosome-dependent cell death, meanwhile lysosomal dysfunction further disrupts Fe homeostasis and modulates ferroptosis.⁹² Moreover, pro-inflammatory cytokines activate neutrophils,⁹³ inducing nuclear translocation of MPO and ELANE, which compromise the nuclear envelope and facilitate the formation of NETs (composed of MPO, ELANE, CAMP, etc.), resulting in netotic cell death.⁹⁴ Collectively, these interconnected programmed cell death pathways exacerbate systemic inflammation and tissue damage, driving the progression of sepsis. Although these mechanisms are not directly tested in our study, they provide a plausible biological framework linking PCD heterogeneity with sepsis outcomes, suggesting that coordinated regulation of multiple PCD programs may be critical for therapeutic intervention.

However, there are several limitations should be acknowledged. First, our model was developed using retrospective datasets, potential confounding factors, including patient age, comorbidities, and sample heterogeneity, could influence the observed associations between PCD patterns and clinical outcomes. Second, while our molecular docking results offer promising targets, in vitro and in vivo validation of these compound-protein interactions is necessary. Finally, as our analyses are primarily observational and computational, causality cannot be directly inferred, and model-based predictions may be affected by inherent limitations of the datasets.

Conclusion

This study comprehensively delineated the landscape of 14 PCD patterns in sepsis. We identified two distinct subtypes, and established a robust seven-gene prognostic signature (including ELANE, CTSG, MPO, CAMP, TFRC, IL1B, and CASP5) with poor outcome, which was significantly elevated in sepsis and decreased after treatment in the validation cohorts. Moreover, SMYAD might interact with key PCD proteins according to the molecular docking analysis, supporting its potential as a multi-target adjunctive therapy. Collectively, these findings provide a framework for early risk stratification and biomarker-guided interventions in sepsis. Future research should validate these results in larger and prospective cohorts, and assess the therapeutic efficacy of SMYAD through preclinical and clinical trials.