We used the Aortic Trauma Foundation (ATF) Blunt Thoracic Aortic Injury (BTAI) registry, which includes detailed demographic, injury, physiologic, imaging, management, and outcome variables for patients with confirmed BTAI. The source dataset contained 1,375 patients. For the primary analysis cohort, we excluded patients missing more than 50% of the prespecified physiologic variables (defined below) to limit over-imputation from sparse records; 1,083 patients remained for analysis.

Phenotyping targeted eleven admission variables selected a priori to capture early hemodynamic, neurologic, metabolic, and injury-burden signals: systolic blood pressure (SBP), heart rate (HR), Glasgow Coma Score (GCS), temperature, lactate, base deficit value, pH, creatinine, hemoglobin, platelet count, and Injury Severity Score (ISS). Data were inspected for plausibility and harmonized. To preserve true extremes seen during resuscitation, we allowed SBP and HR values down to 0 (e.g., transient loss of pulses) and used wide physiologic bounds for other labs (e.g., lactate up to 50 mmol/L). Non-physiologic entries and unit inconsistencies (e.g., platelet counts recorded in thousands) were corrected when resolvable or set to missing if not. After excluding patients with >50% missing across physiologic variables, remaining missing values were imputed by variable-wise medians and all features were standardized (z-scores) for modeling.

We first performed unsupervised clustering on the full 1,375-patient cohort using K-means clustering, an unsupervised machine learning algorithm that partitions observations into a prespecified number of clusters by minimizing within-cluster variance. This approach groups patients with similar multivariable physiologic profiles without using outcome information and has been widely applied in clinical phenotyping studies across trauma and critical care populations (12–14).

To determine the appropriate number of clusters, we evaluated solutions ranging from k = 2 to k = 6 and selected k = 3 based on silhouette score optimization (silhouette = 0.23). The silhouette score quantifies how similar each observation is to its assigned cluster relative to other clusters, with values ranging from −1 to 1; higher values indicate better separation, whereas values near zero suggest overlap. In complex physiologic datasets with overlapping clinical features, modest silhouette values are common and do not preclude clinically meaningful subgroup identification.

The three clusters were labeled post hoc according to their dominant physiologic characteristics. One cluster demonstrated preserved vital signs and low lactate burden and was labeled “Stable”. A second cluster was characterized by hypotension, tachycardia, and elevated lactate and was labeled “Shock”. The third cluster showed markedly depressed Glasgow Coma Score (GCS) with intermediate hemodynamics and was labeled “Neurologically Compromised.” Principal component analysis (PCA) of standardized physiologic variables was used for visualization only to illustrate separation between clusters and was not used for model derivation.

Because physiologic clustering is sensitive to missing data patterns and sample composition, we derived phenotype centroids using the full 1,375-patient cohort after imputation, then applied these fixed centroids to a stricter analysis cohort. For the primary analysis, we excluded patients with more than 50% missing physiologic variables, yielding 1,083 patients.

Each patient in this analysis cohort was assigned to the nearest phenotype centroid from the k = 3 derivation model using Euclidean distance in standardized feature space. This centroid-based nearest-neighbor classification approach preserves the original cluster definitions while allowing consistent phenotype assignment in a dataset with reduced missingness. Importantly, this avoids re-running clustering on a smaller or altered cohort, which could result in label drift and compromise interpretability and reproducibility.

This two-step approach, first an unsupervised derivation followed by centroid-based assignment, has been successfully used in other biomedical domains, including vaginal microbiota classification and population health phenotyping, to enable reproducible application of data-driven phenotypes across cohorts (18–20).

The primary outcome was in-hospital mortality. Secondary outcomes included treatment modality [medical management, thoracic endovascular aortic repair (TEVAR), open surgical repair], hospital length of stay (LOS), intensive care unit (ICU) LOS, ventilator days, and time from admission to definitive repair (open or endovascular). When multiple treatment fields were present, mutually exclusive “initial modality” was assigned with a prespecified priority of open repair > TEVAR > medical management; otherwise the record was classified as none/unknown.

Categorical variables were compared across physiologic phenotypes using χ² tests or Fisher's exact tests when expected cell counts were <5. Continuous variables were summarized as medians with interquartile ranges and compared across phenotypes using one-way analysis of variance (ANOVA) or nonparametric alternatives, as appropriate.

Prior to ANOVA, distributional assumptions were evaluated through visual inspection of residual plots and assessment of variance homogeneity across groups. When deviations from normality or heteroscedasticity were observed, the Kruskal–Wallis test, a rank-based nonparametric alternative to ANOVA, was used to compare distributions across phenotypes. When overall group differences were significant, post hoc pairwise comparisons were performed using Dunn's test with Bonferroni correction to control for multiple testing.

In-hospital mortality and treatment modality were analyzed using contingency tables. Overall differences across phenotypes were assessed with χ² tests, and when significant, pairwise comparisons were conducted using Fisher's exact tests with Bonferroni adjustment. Fisher's exact test was selected for pairwise analyses because it provides an exact p-value and remains valid when cell counts are small, which occurred in several subgroup comparisons (e.g., open repair and none/unknown treatment categories). This approach ensured robust inference while controlling for multiple comparisons.

Time-to-intervention variables exhibited skewed distributions and were analyzed using the Kruskal–Wallis test, followed by Dunn's post hoc testing when appropriate. Two-sided p values < 0.05 were considered statistically significant after adjustment.

All analyses were performed in Python 3.12.4 (pandas, scikit-learn, scipy, scikit-posthocs), and figures were generated with matplotlib (Python 3 Reference Manual. Scotts Valley, CA: CreateSpace). PCA plots depict cluster separation but were not used for modeling.

To enhance the interpretability of the unsupervised physiologic clustering, we conducted a post hoc explainable artificial intelligence (XAI) analysis using a supervised surrogate model and SHapley Additive exPlanations (SHAP). XAI methods aim to improve transparency of complex machine learning models by quantifying how individual input features contribute to model outputs. The goal was to assess how well a physiologic classifier could reproduce the K-means–derived phenotypes and to quantify the relative contribution of each physiologic variable to phenotype assignment. Importantly, this step was designed to explain the clustering structure rather than to introduce a novel predictive model.

We trained an XGBoost gradient boosting classifier, a tree-based ensemble learning method well-suited for tabular clinical data, on the standardized admission physiologic variables used for clustering: systolic blood pressure, heart rate, Glasgow Coma Score (GCS), temperature, lactate, base deficit, pH, creatinine, hemoglobin, platelet count, and Injury Severity Score (ISS). The outcome was the three-class phenotype label (Stable, Shock, Neurologically Compromised) originally assigned by K-means. The dataset was randomly split 80:20 into training and held-out test sets with stratification by phenotype. Model performance was evaluated on the held-out test set using accuracy, macro–F1 score, Cohen's κ, and the adjusted Rand index (ARI), which together quantify classification performance and agreement with the underlying cluster structure beyond chance, to quantify agreement with the underlying cluster structure. Five-fold cross-validation was performed on the training set to reduce optimism from a single data split and to estimate generalizable model fidelity.

Global and local feature attributions were obtained using the SHAP TreeExplainer. SHAP is a game theory–based approach that assigns each feature a contribution value based on its marginal impact on model predictions (21). Global importance was defined as the mean absolute SHAP value for each feature across all patients. Class-specific SHAP beeswarm plots were generated to visualize both the magnitude and direction of each feature's contribution within each phenotype. To assess stability of the feature importance rankings, we performed 200 bootstrap refits of the XGBoost model, recomputed global mean absolute SHAP values for each refit, and calculated pairwise Spearman rank correlations between bootstrap-derived rankings and the original model.

After exclusions for >50% missing physiologic data, 1,083 patients remained for analysis. The three physiological phenotypes had been defined in the full 1,375-patient derivation cohort using unsupervised clustering. The silhouette score for the three-cluster solution was modest (0.23), reflecting partial overlap between physiologic profiles rather than discrete separation. To classify the 1,083-patient analysis cohort, each patient's physiologic profile was compared with the average profile (“centroid”) of each phenotype from the derivation step, and they were assigned to the closest match. This preserved the original phenotype definitions and allowed consistent classification without re-running the clustering. The resulting groups showed distinct physiologic patterns and in-hospital mortality rates. One group demonstrated preserved vital signs and low lactate (“Stable,” n = 431, 39.8%), another had hypotension, tachycardia, and elevated lactate (“Shock,” n = 398, 36.7%), and the third had markedly depressed GCS with intermediate hemodynamics (“Neurologically Compromised,” n = 254, 23.5%). These descriptive labels were assigned post hoc by the investigators to reflect the defining characteristics of each group. Principal component analysis (PCA) of standardized physiologic variables demonstrated clear visual separation of the three phenotypes (Figure 1).

Baseline demographic, injury, and physiologic characteristics stratified by phenotype are summarized in Table 1. The three groups differed significantly across multiple admission physiologic variables, including systolic blood pressure, heart rate, lactate, Glasgow Coma Score, and Injury Severity Score (all p < 0.001), consistent with their descriptive phenotype labels. The “Stable” phenotype demonstrated preserved hemodynamics and lower injury burden, whereas the “Shock” phenotype exhibited hypotension, tachycardia, and the highest lactate levels. The “Neurologically Compromised” phenotype showed markedly depressed neurologic status with intermediate hemodynamic derangements and the highest overall injury severity.

Demographic characteristics were largely similar across phenotypes, though modest differences were observed in smoking status, hypertension prevalence, coronary artery disease, and interhospital transfer rates (Table 1). These findings indicate that the observed phenotypic separation was driven primarily by physiologic and injury-related variables rather than baseline demographic differences.

Patients in the “Stable” group had preserved hemodynamics [median SBP 131.0 [118.0–145.5] mmHg, HR 88.0 [78.0–100.0] bpm] and low lactate [2.8 (1.9–4.0) mmol/L]. Neurologic status was generally normal [median GCS 15.0 (12.5–15.0)]. This group had the lowest injury burden [median ISS 29.0 (24.0–38.0)] and the highest prevalence of hypertension (25.5%). Smoking was reported in 18.8%, and coronary artery disease in 5.3%.

“Shock” patients were defined by hypotension [SBP 101.0 (80.0–128.0) mmHg], tachycardia [HR 110.0 (90.0–129.0) bpm], and the highest lactate levels [5.3 (3.7–7.6) mmol/L]. Injury severity was greater than in Stable patients [ISS 34.0 (24.0–43.0)]. Hypertension was less common (16.1%), and smoking prevalence was 15.3%. Coronary artery disease was rare (2.0%).

This group was characterized by markedly depressed neurologic function [median GCS 14.0 (3.0–15.0)] with intermediate hemodynamics [SBP 102.5 [90.0–115.0] mmHg, HR 104.5 [93.0–117.0] bpm] and moderate lactate elevation [3.8 (2.8–5.0) mmol/L]. Injury severity was highest overall [ISS 36.0 (29.0–43.0)]. Hypertension prevalence was 17.7%, smoking was least common (10.6%), and coronary artery disease occurred in 3.1% of patients.

Mortality differed significantly by phenotype (χ² = 27.821, p < 0.001). Observed mortality was 7.2% in “Stable” (30/431), 18.9% in “Shock” (72/398), and 18.1% in “Neurologically Compromised” (46/254) (Figure 2). Pairwise Fisher's exact tests with Bonferroni correction demonstrated significantly higher mortality in both “Shock” vs. “Stable” (adj. p = 0.000004) and “Neurologically Compromised” vs. “Stable” (adj. p = 0.000041), with no difference between “Shock” and “Neurologically Compromised” (adj. p = 1.000).

Treatment modality varied significantly by phenotype (χ² = 37.415, p < 0.001), with TEVAR being the most frequent intervention across all groups but with differing relative distributions. “Stable” patients most often underwent TEVAR (58.2%) or medical management (31.1%), whereas “Shock” and “Neurologically Compromised” patients had lower use of medical management (19.1% and 22.0%, respectively) and higher proportions of “none/unknown” treatment (20.1% and 11.4%, respectively). Pairwise testing confirmed that “Stable” patients were significantly more likely than “Shock” patients to receive medical management (adj. p = 0.0010) and less likely to have no documented intervention (adj. p < 0.0001). “Shock” patients also had a higher rate of no documented intervention compared with “Neurologically Compromised” patients (adj. p = 0.0437) (Figures 3, 4).

Within-phenotype mortality analyses demonstrated strong associations between treatment modality and survival. In the “Shock” phenotype, mortality was highest with open repair (40.0%) and “none/unknown” treatment (33.8%), intermediate with TEVAR (13.7%), and lowest with medical management (15.1%) (χ² = 24.745, p < 0.0001). In “Neurologically Compromised” patients, mortality was highest with open repair (77.8%) and “none/unknown” treatment (57.1%), compared with much lower mortality after TEVAR (9.8%) or medical management (14.5%) (χ² = 56.293, p < 0.0001). “Stable” patients had low overall mortality, but deaths were more frequent in the “none/unknown” group (21.6%) compared with TEVAR (7.1%) or medical management (3.7%) (χ² = 14.432, p = 0.0024).

Hospital LOS showed a nonsignificant trend across phenotypes (ANOVA F = 2.235, p = 0.107). Exploratory post hoc testing indicated slightly longer hospital stays in “Neurologically Compromised” vs. “Stable” (Bonferroni-adjusted p = 0.0426), whereas differences between “Shock” and the other groups were not significant (Figure 5). ICU LOS differed significantly overall (ANOVA F = 4.854, p = 0.008); “Neurologically Compromised” required longer ICU care than “Stable” (p = 0.000587), with no significant differences between “Shock” and the other groups after adjustment (Figure 6).

Ventilator days did not differ significantly overall (ANOVA F = 1.571, p = 0.208). Exploratory post hoc comparisons suggested more ventilator days in “Neurologically Compromised” vs. “Stable” (p < 0.0001) and in “Shock” vs. “Stable” (p = 0.000517); these pairwise findings should be interpreted cautiously given the nonsignificant overall test.

Time from admission to open repair did not differ significantly across phenotypes (Kruskal–Wallis H = 2.400, p = 0.301). Median time to open repair was shortest in the “Neurologically Compromised” group [1.3 h (IQR 1.0–4.5)] compared with “Shock” [8.9 h (0.6–134.7)] and “Stable” [13.0 h (4.5–25.7)], though pairwise Dunn's tests with Bonferroni correction found no statistically significant differences between groups (all adjusted p > 0.36).

Similarly, time to TEVAR did not differ significantly (Kruskal–Wallis H = 2.222, p = 0.329). Median time to TEVAR was 4.5 h (3.0–14.0) for “Neurologically Compromised”, 6.0 h (3.0–25.0) for “Shock”, and 5.5 h (3.0–16.0) for “Stable”. Dunn's post hoc testing again revealed no significant pairwise differences (all adjusted p > 0.44).

The XGBoost surrogate classifier demonstrated high fidelity to the original K-means–derived phenotype labels. Because the supervised model was trained on the same physiologic variables used for clustering, its high classification performance was expected and is presented solely to support explainability, not predictive novelty. On the held-out test set, the model achieved an accuracy of 0.917, macro–F1 score of 0.906, Cohen's κ of 0.873, and an adjusted Rand index (ARI) of 0.788, indicating strong agreement with the unsupervised cluster structure. Five-fold cross-validation on the training set yielded similar results (mean accuracy 0.913, macro–F1 0.905, κ 0.867, ARI 0.770), supporting the generalizability of the classifier.

Global SHAP-based feature importance (Figure 7) ranked systolic blood pressure (mean |SHAP| = 1.16), lactate (0.87), and heart rate (0.85) as the dominant physiologic drivers of phenotype assignment, followed by temperature (0.14), Glasgow Coma Score (0.14), creatinine (0.11), base deficit (0.11), hemoglobin (0.10), pH (0.08), Injury Severity Score (0.08), and platelet count (0.07). Stability analysis across 200 bootstrap refits demonstrated high reproducibility of these rankings, with a median Spearman rank correlation of 0.845 between bootstrap-derived importance orders and the original model.

Class-specific SHAP beeswarm plots (Figure 8) illustrated distinct physiologic signatures within each phenotype. The Shock phenotype was driven by low systolic blood pressure and high lactate and heart rate, whereas the Stable phenotype was characterized by higher blood pressure and lower lactate. The Neurologically Compromised phenotype exhibited intermediate hemodynamics with relatively higher creatinine and base deficit and modestly lower GCS, features that primarily differentiated this group from the others rather than contributing uniformly across the full cohort. Together, these XAI findings confirm that the clusters reflect physiologically coherent patterns primarily along axes of perfusion and metabolic derangement, with neurologic status providing additional within-phenotype refinement.