This study included 41 HS patients enrolled from the Department of critical care medicine of the 908th Hospital of PLA joint logistic support force, from June 2022 to February 2024. The patients were divided into two groups: 11 cases in the HSIC group and 30 cases in the NHSIC group. The study was conducted in accordance with the ethical guidelines outlined in the Declaration of Helsinki, and informed consent was obtained from all participants or their family members. The study was approved by the ethics committee of our hospital.

Patients were diagnosed with HS or HSIC based on the Bouchama’s HS criteria (Bouchama and Knochel, 2002) and a newly proposed HSIC score (Lin et al., 2024). HS was defined by meeting the following criteria: (1) medical history of exposure to high temperature, high humidity or high-intensity exercise; (2) core temperature over 40°C; (3) CNS-related changes, including coma, convulsions, delirium, or abnormal behavior. HSIC scores were calculated according to as follows: Core temperature <40°C was scored as 0 points, ≥40°C but <42°C as 1 point, ≥42°C as 2 points; D-dimer of 2-fold increase over baseline was scored as 0 point, ≥ 2-fold but < 5-fold increase over baseline as 1 points, ≥ 5-fold increase over baseline as 2 points; prothrombin time (PT) prolongation of <2 s was scored as 0 point, ≥2 s but <4 s as 1 point, ≥4 s as 2 points. HSIC was diagnosed as a total score at least 3. Patients were excluded from our study if they were younger than 18 years old, they had a congenital coagulation disorder or chronic disease of the liver, or they were using anticoagulant drugs at the time of enrollment.

Venous blood samples were drawn in the early morning after a minimum of 8 h of fasting and placed into a dry blood collection tube. The samples were then centrifuged at 4244 g for 10 min at 4°C. The serum obtained was separated and stored at −80°C.

In our study, we analyzed serum samples from 41 participants using DIA quantitative proteomics. The total protein from each sample was extracted, with a portion used for protein concentration determination and SDS-PAGE, while the remainder underwent trypsin digestion. After desalting, LC-MS/MS was employed to identify the peptides present. LC-MS/MS analysis was conducted on a timsTOF Pro 2 system (Bruker) coupled with Evosep One LC, employing both DDA (scan range m/z 100–1,700, ion mobility 0.75–1.35 Vs/cm², 8 PASEF MS/MS cycles) for spectral library generation and 4-window DIA for quantitative profiling. Spectronaut pulsar software was used to comprehensively search the raw data and compare it against the known UniProtKB human proteome database. The false discovery rate (FDR) for both DDA and DIA data was set at 0.01. DEPs were defined with a p-value <0.05 and |log2(FC)| >1.

GO and KEGG pathway analyses were used to assess protein functions. In the GO annotation process, the protein sequences of the differentially expressed proteins were first searched locally using the NCBI BLAST + client software (ncbi-blast-2.2.28+-win32.exe) and InterProScan to identify homologous sequences. Gene ontology (GO) terms were then mapped, and the sequences were annotated using Blast2GO software, with the results visualized using R scripts. For KEGG annotation, the annotated proteins were blasted against the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://geneontology.org/) to obtain their KEGG orthology identifications, which were subsequently mapped to relevant pathways in KEGG. A p-value less than 0.05 in the pathway enrichment test, along with a protein count greater than 5, were used as criteria for determining significance.

Serum samples from 41 subjects were analyzed using metabolomics. Plasma samples were thawed at 4°C, and 100 µL of each sample was mixed with 400 µL of pre-cooled methanol/acetonitrile (1:1, v/v) to precipitate proteins. The mixture was centrifuged at 14,000 g and 4°C for 20 min. The supernatant was dried under vacuum centrifugation and re-dissolved in 100 µL of acetonitrile/water (1:1, v/v). After centrifugation at 14,000 g and 4°C for 15 min, the supernatant was injected for analysis. An UHPLC system (1,290 Infinity LC, Agilent Technologies) coupled to a quadrupole time-of-flight mass spectrometer (AB Sciex TripleTOF 6600) was used. Hydrophilic interaction chromatography (HILIC) was performed using a 2.1 mm × 100 mm ACQUIY UPLC BEH Amide 1.7 µm column (Waters, Ireland). In both ESI positive and negative modes, mobile phases consisted of A = 25 mM ammonium acetate and 25 mM ammonium hydroxide in water, and B = acetonitrile. The gradient started at 95% B for 0.5 min, decreased to 65% in 6.5 min, to 40% in 1 min (held for 1 min), then returned to 95% in 0.1 min with a 3 min re-equilibration. Source parameters: Gas1 = 60, Gas2 = 60, CUR = 30, temperature = 600°C, ISVF = ±5500 V. MS range: 60–1,000 Da (0.20 s/spectra); MS/MS: 25–1,000 Da (0.05 s/spectra), CE = 35 V (±15 eV), DP = ±60 V, IDA with 10 candidate ions per cycle. An UHPLC (Vanquish, Thermo) coupled to an Orbitrap (Q Exactive HF-X/Q Exactive HF) was also employed. The same HILIC column and mobile phases were used. Gradient: 98% B for 1.5 min, decreased to 2% in 10.5 min, held for 2 min, then returned to 98% in 0.1 min with a 3 min re-equilibration. MS range: 80–1,200 Da, resolution 60,000 (100 ms); MS/MS: 70–1,200 Da, resolution 30,000 (50 ms), exclude time: 4 s. Raw MS data were converted to MzXML files using ProteoWizard MSConvert, imported into XCMS for peak picking (centWave m/z = 10 ppm, peakwidth = c (10, 60)), and annotated with CAMERA. Compound identification was based on m/z and MS/MS comparison with an in-house database. The positive and negative data were merged to create a combined dataset, which was then imported into R software. After sum-normalization, the data were analyzed using the R package (ropls), performing multivariate analysis including PCA and OPLS-DA with Pareto scaling. 7-fold cross-validation and 200 response permutation testing assessed model robustness. VIP values were calculated to determine each variable’s contribution to classification. Student’s t-test was used to evaluate differences between two groups, with VIP >1 and p < 0.05 for screening significant metabolites.

All statistical analyses were performed using R software (version 4.1.2). Continuous variables were presented as means ± SD or medians with IQR, depending on distribution. Differences between groups were assessed using t-tests for data following a normal distribution or Mann-Whitney U tests for data following a non-normal distribution. Categorical variables were compared using chi-square or Fisher’s exact tests. A p-value <0.05 was considered significant. Feature selection was conducted with LASSO and Boruta algorithms. Predictive models were developed using LR, XGBoost, and SVM. Model performance was evaluated with metrics such as accuracy, sensitivity, specificity, AUC, kappa, and F1 score. ROC and PR curves were used to assess discriminative power, while DCA validated clinical utility. DALEX was used for model interpretability.

A total of 41 patients with HS were enrolled in this study after screened for eligibility. All of patients with HS were divided into HSIC group (n = 11, 9 males and 2 females), and NHSIC group (n = 30, 25 males and 5 females). Inter-group comparisons revealed statistically significant differences in temperature, heart rate, white blood cells (WBC), platelets (PLT), aspartate transaminase (AST), ALT (alanine transaminase), total bilirubin (TBIL), creatinine (Cr), lactate (Lac), creatine kinase (CK), MB(myoglobin), PT, INR (international normalized ratio), APTT (activated partial thrombin time), TT (thrombin time), d-dimer, APACHE II(Acute Physiology and Chronic Health Evaluation II) score, GCS(Glasgow Coma Scale) score and MODS (multiple organ dysfunction syndrome). While age, gender, risk of coronary disease and diabetes, probability of hypertension, Hemoglobin, HCT (hematocrit), and FIB (fibrinogen) showed no statistical differences between the two groups. The comparisons between HSIC group and NHSIC group are depicted in Table 1.

Compared with NHSIC group, the APACHE II score in HSIC group is significantly increased (Table 1; Figure 2A). The cumulative risk curve of adverse outcome according to the study groups is shown in Figure 2B. The result demonstrated that individuals in NHSIC group were at the higher risk of death compared with HSIC group.

We identified a total of 1907 proteins. For differential expression analysis, p value <0.05 were considered statistically significant and |log2(FC)| >1 was considered up- or downregulated, respectively. A total of 125 DEPs were screened in the NHSIC group vs. HSIC group. A comparison of NHSIC and HSIC patients found that there were 6 and 119 upregulated and downregulated proteins, respectively (Figure 3A). Hierarchical clustering algorithm was used to perform cluster analysis on the DEPs of NHSIC and HSIC groups, and the data was displayed in the form of heat map (Figure 3B).

We performed GO functional annotation on the differentially expressed proteins (DEPs) (Figure 4A). In the NHSIC vs. HSIC group, the difference in protein expression was greatest in Biological Process (BP), Molecular Function (MF), and Cellular Component (CC). Concentrated results were in negative regulation of coagulation, lipoprotein particle receptor binding, and lipoprotein particle. We used KEGG to analyze the signaling pathways of 125 DEPs in the NHSIC vs. HSIC group that were mainly concentrated in complement and coagulation cascades, cholesterol metabolism, and neuroactive ligand−receptor interaction (Figure 4B). From the 125 DEPs in the NHSIC vs. HSIC group, 36 proteins were involved in protein interactions. Furthermore, PPIs analysis showed that APOH and PLG had the highest degree scores in the network (Figure 4C), suggesting that they play a role in the pathological process of HS.

LASSO regression introduces a penalty function to continuously compress the coefficients, streamline the model, and avoid multicollinearity and overfitting, thereby achieving the effect of variable selection. When a standard error of the minimum distance λ was 0.012, 27 feature coefficients were nonzero (Figures 5A–C). The process and results of feature selection using Random Forest are shown in Figures 5D,E, which identified 16 important variables predisposing to HSIC in heatstroke patients. The two methods identified overlapping proteins, two upregulated (LDHA and NGAL), two downregulated (Prothrombin and GBE) (Figure 5F). Then, these four most vital variables were selected for further analysis. As the principal component analysis result in Figure 5G, the four DEPs aforementioned can clearly distinguished NHSIC and HSIC, which indicated that they may play key roles in the diagnosis of HSIC. The correlations of the proteins were also analyzed as shown in Figure 5H. The absence of significant correlations among these proteins indicates they do not exhibit functional similarities.

XGBoost, LR and SVM algorithms were used to construct models based on selected proteins. We used AUC, accuracy, no information rate, balanced accuracy, kappa, sensitivity, specifcity, precision, and F1 scores to comprehensively evaluate the model’s performance. XGBoost had the largest AUC (0.991) and precision (1.0), followed by LR (AUC: 0.979; precision: 0.90) and SVM (AUC: 0.976; precision: 0.750) (Table 2). Figure 6A described the ROC curves for the three models. The accuracy, kappa, sensitivity, Balanced accuracy and F1 scores of LR were higher than those of XGBoost and SVM, as shown in Table 2. Therefore, The LR had better clinical utility compared to XGBoost and SVM.

The ROC curves and PR curves for the three models indicated that the LR was better than XGBoost and SVM in discrimination (Figures 6A,B). The calibration curve fit was good (Figure 6C), indicating that the model had good diagnostic performance. Furthermore, we found the net benefits (NB) of LR model was higher than 0, with the greater NB in DCA clinical evaluation (Figure 6D), indicating the importance of the LR model for HSIC diagnosis. The DALEX package was used for logistic regression analysis to further demonstrate the importance of the four proteins in the model and the descending order of importance of these features in the model was LDHA, NGAL, Prothrombin and GBE (Figure 6E). Based on the confusion matrix, the assay had a sensitivity, specificity, PPV, and NPV of 81.8%, 96.7%, 90.0%, and 93.6% respectively (Figure 6F).

We present one instances in which interpretability analyses of the model predictions. The ML model predicted a 45.7% risk of HSIC based on four critical predictors. We found prothrombin and GBE being the two contributors to the increased risk of HSIC, whereas LDHA and NGAL reduced the model’s diagnosis of HSIC (Figure 6G). We developed the interface to facilitate the use of the model to explore the relative contribution of HSIC probability factors in ICU patients. In the diagnosis view, the system invokes a diagnosis model, and the LR model diagnoses the patient’s HSIC. The analysis results are visualized in a graphic view, which indicates the HIC probability of the patient input important features values (Figure 6H).

In addition, we also described the effect (positive or negative) of four proteins on the model. Figure 7 showed the relationship between LDHA, Prothrombin, NGAL, GBE and predicted HSIC. Higher LDHA and NGAL were associated with an increased risk of HSIC. Lower Prothrombin and GBE were associated with an increased risk of HSIC (Figure 7).

Figure 8A displays the results of principal component analysis (PCA), which clearly separates the samples into distinct groups. Figure 8B presents the OPLS-DA scores, highlighting a significant difference between the two groups. The robustness of the OPLS-DA model was validated through permutation analysis, yielding an R²Y of 0.979 and a Q² value of 0.545. In Figure 8C, the negative Q² value confirms the absence of over-fitting, further validating the model’s reliability and effectiveness. Differential expression analysis, based on a Benjamini–Hochberg adjusted filter of <0.05 and a log2 fold change (FC) > 1.0, revealed 110 significantly DEMs, including 30 upregulated and 80 downregulated metabolites. These 110 DEMs are shown in the volcano plot in Figure 8D. Finally, Figure 8E illustrates the top 25 enriched metabolic pathways for these DEMs, as identified using the KEGG database.

Through DIABLO algorithm, we identified distinct clustering patterns differentiating NHSIC and HSIC samples (Figure 9A). Subsequent integrative modeling of multi-omics features revealed strong cross-correlations between proteomic and metabolomic profiles, highlighting their coordinated biological regulation (Figure 9B). To identify significant protein-metabolite interactions, we extracted the top 20 differentially expressed proteins (DEPs) and metabolites (DEMs) ranked by p-values and performed Pearson correlation analysis to evaluate their associations, followed by the generation of a correlation heatmap (Figure 9D). Multi-omics clustering analysis revealed strong inter-dataset associations between features derived from the proteomic and metabolomic profiles (Figure 9D). Additionally, the differentially expressed proteins and metabolites were mapped to the KEGG database to assess the potential relationship. It was found that the AMPK signaling pathway, cholesterol metabolism, glycolysis/gluconeogenesis, neuroactive ligand-receptor interaction, propanoate metabolism, glucagon signaling pathway, glycerophospholipid metabolism, and pantothenate and CoA biosynthesis were significantly enriched (Figure 9E).