This study draws data from two primary sources: the eICU Collaborative Research Database (eICU-CRD) and The Medical Information Mart for Intensive Care IV database (MIMIC-IV). The eICU-CRD database encompasses various ICU units across the United States, offering a comprehensive array of clinical data, physiological parameters, and medical events. Spanning from 2014 to 2015, it meticulously documents information for over 200,000 patients, facilitating medical research endeavors and data-informed clinical decision-making (22). On the other hand, MIMIC-IV (version 2.2) represents an extensive repository of intensive care data, featuring detailed records of more than 190,000 ICU patients from 2008 to 2019 (23). This database is characterized by its exhaustive collection of clinical details, including demographic profiles, laboratory findings, and medication histories, serving as invaluable resources for rigorous clinical investigations.

We identified patients with heart failure complicated by sepsis from both the eICU-CRD and MIMIC-IV databases using ICD-9 and ICD-10 codes. The exclusion criteria for the study population are: (1) age under 18 years, (2) ICU stay duration less than 24 h, and (3) clinical information missing rate exceeding 30% at data collection. For patients with multiple hospital admissions or ICU visits, only the first ICU experience during the initial hospital admission is considered. Heart failure was defined as a syndrome resulting from structural or functional cardiac abnormalities that lead to inadequate cardiac output and congestion in the systemic or pulmonary circulation, encompassing all types of heart failure with different ejection fractions. Sepsis was diagnosed based on the Sepsis-3.0 guidelines, which define it as life-threatening organ dysfunction caused by a dysregulated host response to infection. A SOFA score ≥ 2 (or a qSOFA score ≥ 2 for suspected infection in non-ICU settings) was used to diagnose sepsis.

In this study, we included ICU patients diagnosed with heart failure and sepsis, and extracted the following data: (1) Demographics: age, gender, height, and weight; (2) Vital Signs: temperature (T), heart rate (HR), respiratory rate (R), systolic blood pressure (SBP), diastolic blood pressure (DBP), mean blood pressure (MBP), and peripheral oxygen saturation (SpO2); (3) Laboratory parameters: complete blood count, liver and kidney function tests, electrolytes, lipid profile, blood gas analysis, coagulation function, cardiac enzymes, and BNP; (4) Comorbidities: hypertension, diabetes mellitus, hyperlipidemia, chronic obstructive pulmonary disease (COPD), pneumonia, chronic kidney disease (CKD), and atrial fibrillation (AF); (5) Medication data: angiotensin-converting enzyme inhibitors/angiotensin II receptor blockers (ACEI/ARB), beta Blockers, furosemide, spironolactone, dobutamine, dopamine, epinephrine, milrinone, norepinephrine, and phenylephrine; and (6) Other Indicators: Ventilation and 24-h fluid balance. The primary outcome is the 28-day all-cause mortality rate.

Initially, we transformed certain indicators, such as computing BMI from height and weight and determining 28-day in-hospital mortality using hospitalization duration and survival status. Variables with over 30% missing data were removed, and missing values in the remaining features were imputed using KNN. Outliers were identified using the 1.5 times interquartile range method, particularly focusing on BMI, mechanical ventilation time, and 24-h fluid balance, and were subsequently removed. Additionally, Spearman correlation coefficients were calculated to evaluate variable relationships, while VIF values assessed multicollinearity. Variables with high correlation or VIF exceeding 5 underwent pre-screening. Continuous variables were standardized for model stability, and categorical variables were transformed into dummy variables via one-hot encoding. Despite minor sample imbalances in the outcome variable, we chose not to employ sample balancing techniques.

The Boruta method is used for initial feature screening, determining feature importance by comparing them with randomly generated “shadow features” (24). The XGBoost method is employed for importance ranking of the preliminarily selected features. Model construction and validation are conducted using the EICU dataset, with 10-fold cross-validation to generate training and validation sets, and Logistic Regression (LR), eXtreme Gradient Boosting (XGBoost), Adaptive Boosting (AdaBoost), and Gaussian Naive Bayes (GNB) models are established and validated. Model performance is evaluated on the validation set using metrics such as the area under the curve (AUC) for discrimination, calibration curve for accuracy, and DCA curve for clinical utility, as well as metrics including accuracy, sensitivity, specificity, positive predictive value, negative predictive value, and F1 score. The final predictive model is optimized using hyperparameter tuning and grid search. Additionally, the MIMIC dataset is utilized as external validation data, following the same data processing methods as the EICU dataset, with evaluation based on metrics including AUC, accuracy, sensitivity, and specificity to assess model generalization performance.

SHAP (SHapley Additive exPlanations) is a technique based on game theory’s Shapley values (25). It’s used to interpret machine learning predictions by dissecting the contribution of each feature. This enhances model transparency and ensures fair decision-making. We employ SHAP to analyze the outcomes of our top-performing machine learning model. This method not only identifies crucial features for optimizing model performance but also provides detailed insights through feature contribution charts, summary plots, and explanations for individual predictions. These tools help us understand the extent of each feature’s influence, whether it’s positive or negative, and how they collectively impact model outcomes.

For continuous variables, we display using mean and standard deviation, and comparison is done using t-tests (or Wilcoxon rank-sum tests); for categorical variables, presentation is in percentages, and comparison is conducted using chi-square tests (or Fisher’s exact tests). A p-value <0.05 is deemed statistically significant. All statistical analyses were performed using R version 4.2.3 and Python version 3.11.4.

According to the inclusion and exclusion criteria, our study cohort comprised a total of 6819 patients with heart failure and sepsis. Among these, 3891 cases from the eICU-CRD were used for model construction, while 2928 cases from the MIMIC-IV database were used for external validation. As shown in Figure 1, the screening process is illustrated. During the selection process, patients with ICU stays less than 1 day or under 18 years old were excluded. Subsequently, data processing involved removing outliers and handling missing values. In the eICU-CRD database, 560 patients (14.4%) died within 28 days, compared to 660 patients (22.5%) in the MIMIC-IV database. Differences in baseline characteristics are summarized in Table 1. In the eICU-CRD database, compared to the survival group, patients in the death group exhibited higher age, white blood cell count, neutrophil count, TBIL (total bilirubin), ALT (alanine aminotransferase), BUN (blood urea nitrogen), respiratory rate, fluid balance, and mechanical ventilation time, and lower BMI (body mass index), calcium, blood pressure, and peripheral oxygen pressure. Differences in comorbidities, such as atrial fibrillation, hypertension, and pneumonia, were also observed between the two groups. Additionally, there were differences in medication usage between the two groups, including the use of ACEI/ARB (ACE inhibitors/angiotensin receptor blockers), beta-blockers, furosemide, spironolactone, dobutamine, dopamine, epinephrine, norepinephrine, and phenylephrine.

We eliminated features with a missing rate exceeding 30%, as demonstrated in Appendix Figure 1. Features with notably high missing rates are primarily found in laboratory tests such as cardiac enzymes, blood gas analysis, lipid profile, and coagulation function. Additionally, guided by the correlation heatmap showing features with correlation coefficients greater than 0.5 and features with VIF exceeding 5, as illustrated in Appendix Figure 2, we conducted further screening. Prior to model construction, we excluded features with high correlation and VIF, including hemoglobin, lymphocyte count, chloride, aspartate aminotransferase, blood urea nitrogen, mean blood pressure, and diastolic blood pressure.

The Boruta method, based on random forests, assesses feature importance by comparing original features with randomly generated “shadow features.” We applied Boruta for initial feature selection, as depicted in Figure 2. Green denotes important features included in the model to enhance predictive capability; red represents unimportant features excluded from consideration; yellow indicates features with uncertain importance requiring further investigation. Blue represents shadow features for comparison but not used in model training. Boruta identified 22 initial features, including age, WBC (white blood cell count), NE (neutrophil count), MONO (monocyte count), PLT (platelet count), sodium, calcium, TBIL (total bilirubin), alanine aminotransferase, creatinine, T (temperature), R (respiratory rate), SBP (systolic blood pressure), oxygen saturation, ventilation time, BMI, atrial fibrillation, dopamine, epinephrine, norepinephrine, phenylephrine, and ventilation.

The XGBoost algorithm ranks feature importance based on split frequency and gain in decision trees. Appendix Figure 3 shows our feature importance ranking using XGBoost. The top 10 variables, in descending order of importance, are: age, ventilation, norepinephrine, WBC, TBIL, T, phenylephrine, R, NE, SBP.

This study utilized four binary classification machine learning algorithms, Logistic Regression (LR), eXtreme Gradient Boosting (XGBoost), Adaptive Boosting (AdaBoost), and Gaussian Naive Bayes (GNB), to construct predictive models. Employing the eICU database, we implemented a 10-fold cross-validation technique to establish training and validation sets, followed by evaluation on a separate test cohort. Figure 3 and Table 2 illustrate the performance of these models. The ROC curve (Figure 4A) highlights LR’s superior performance, achieving an AUC of 0.746 in the test cohort, compared to XGBoost (0.726), AdaBoost (0.744), and GNB (0.722). Furthermore, LR’s calibration curve (Figure 4B) closely aligns with the ideal line, indicating excellent calibration. Decision curve analysis (DCA) (Figure 4C) indicates LR’s highest net benefit within the 0–80% threshold range. The precision-recall (PR) curve (Figure 4D) illustrates LR’s higher recall at sustained high precision. Additionally, LR demonstrates robust performance across various metrics, including accuracy (0.685), sensitivity (0.666), specificity (0.712), positive predictive value (0.285), negative predictive value (0.914), and F1 score (0.397). Consequently, we selected the LR algorithm for model construction, incorporating 10 variables: age, ventilation, norepinephrine, WBC, TBIL, T, phenylephrine, R, NE, and SBP. Through hyperparameter tuning and grid search optimization, we established the model parameters as follows: tol (convergence measure): 1e-06, penalty (regularization type): l2, max_iter (number of iterations): 100, C (regularization factor): 1.0.

This study employs the SHAP method to interpret model results, presenting both SHAP summary plots and SHAP force plots. In the SHAP summary plot, the Y-axis represents features, while the X-axis indicates the impact of features on outcomes. Each point represents a sample, with red indicating high-risk values and blue indicating low-risk values. As shown in Figure 3A, the LR model’s feature importance from top to bottom is: age, ventilation, norepinephrine, T, R, TBIL, SBP, WBC, NE, phenylephrine. Older age (red points) correlates with higher SHAP estimated values, predicting an increased risk of mortality. Additionally, higher white blood cell count, total bilirubin, and respiratory rate are associated with increased mortality risk. Patients using ventilation, norepinephrine, and phenylephrine also show increased mortality risk. Furthermore, lower temperature and lower systolic blood pressure are associated with increased mortality risk. In the SHAP force plot, each Shapley value is represented by an arrow, indicating whether it positively (increases) or negatively (decreases) affects the prediction. As illustrated in Figure 3B, increases in white blood cell count, decreases in temperature, and increases in total bilirubin push the predicted mortality risk higher, while younger age and lower neutrophil count push the predicted mortality risk lower. It’s important to note that due to standardization of numerical variables to a mean of 0 and a variance of 1, the data in the plots are not in their original scale.

We selected 2928 cases of patients with heart failure and sepsis from the MIMIC-IV database for external validation. In the MIMIC database, these patients had a 28-day in-hospital mortality rate of 22.5%, slightly lower than that of the eICU-CRD database (14.4%). Prior to external validation, we applied the same data processing methods to the MIMIC data as we did to the eICU-CRD data. The validation results revealed an AUC of 0.699 and a Brier score of 0.169. Additionally, the accuracy, sensitivity, specificity, positive predictive value, negative predictive value, and F1 score were 0.699, 0.156, 0.403, 0.673, 0.648, and 0.261, respectively. With the AUC difference between the external validation and validation/test sets being less than 0.1, we conclude that the LR model demonstrates favorable stability.