Development and Temporal Validation of the SAKARYA Model for Early Cardiorenal Risk Stratification in Critically Ill Adults with Acute Heart Failure: A Credentialed Database Analysis

1. Introduction

Acute heart failure remains a frequent cause of urgent hospitalization, short-term deterioration, and early death across cardiovascular care settings [1,2]. Contemporary European and North American guidance emphasizes rapid syndrome recognition, natriuretic peptide testing, congestion-directed treatment, and early identification of organ dysfunction during the first treatment hours [3,4]. Risk estimation during that interval remains uneven after intensive care unit triage.

Available acute heart failure risk tools were largely derived in emergency department or hospital-wide cohorts, and their outcome definitions, predictor windows, and intended use differ substantially [5,6]. Prospective validation has improved confidence in some models, including EHMRG, but those instruments were not built in critically ill populations and do not directly address early cardiorenal deterioration after intensive care unit admission [7,8].

Cardiorenal interaction is central to acute heart failure. Acute kidney injury is common, worsens short-term prognosis, and reflects the combined effects of venous congestion, impaired perfusion, neurohormonal activation, and treatment intensity [9,10]. An early model that integrates mortality and kidney injury may therefore be more informative than a mortality-only model in the intensive care environment.

Large deidentified electronic health record resources now allow transparent model development in critical care cohorts. The broader PhysioNet framework supports reproducible secondary analyses of physiologic and electronic health record data. The present study developed and temporally validated a parsimonious early model for critically ill adults with an operational acute heart failure phenotype. The objective was to derive the SAKARYA model (representing Sex/Sodium/Saturations, Age, Kidney markers [Creatinine/BUN], Arterial pressure [MAP], Rates [Heart/Respiratory], Yielding early Acid-base/Anemia markers) from variables available within the first 6 hours of intensive care unit admission and to evaluate discrimination, calibration, robustness, and strict comparator performance.

2. Materials and Methods

2.1. Study Design and Data Source

This retrospective cohort study used MIMIC-IV version 3.1, a deidentified critical care database derived from Beth Israel Deaconess Medical Center [11,12]. The source platform is maintained through PhysioNet [13]. Adults were eligible if they were 18 years of age or older, had a first intensive care unit admission, and had a hospital diagnosis code consistent with heart failure. Heart failure coding included ICD-9 code 428* and ICD-10 codes I50*, I110, I130, and I132. Data extraction and cohort assembly were performed programmatically using Structured Query Language (SQL) by a credentialed investigator (H.B.I.) who has completed the required training in human subjects research on the Collaborative Institutional Training Initiative (CITI) platform.

2.2. Cohort Definition

The broad primary cohort required at least one early acute heart failure enrichment feature in addition to heart failure coding. Enrichment was defined by intravenous loop diuretic treatment or BNP/NT-proBNP testing from 6 hours before to 24 hours after intensive care unit admission. This prespecified broad phenotype yielded 8477 patients. After application of the clean core complete-case rules for the 11-variable primary model, 7065 patients remained for modeling, including 5399 in derivation and 1666 in temporal validation. The temporal split used approximate actual year, with 2009-2018 assigned to derivation and 2019-2023 assigned to temporal validation. Figure 1 shows cohort assembly.

A strict sensitivity cohort was defined to address major design concerns. This cohort narrowed phenotype enrichment to the interval from 6 hours before to 6 hours after intensive care unit admission, restricted baseline creatinine to the interval from 24 hours before to 6 hours before intensive care unit admission, and restricted predictor creatinine to the interval from admission to 6 hours after admission. The strict clean core complete-case cohort included 737 derivation patients and 296 temporal validation patients.

2.3. Predictor Definition

Baseline predictors were defined from the window extending from 6 hours before to 6 hours after intensive care unit admission. Candidate variables were selected before model fitting on the basis of clinical availability, physiologic interpretability, and bedside transportability. The final core model included age, sex, heart rate, mean arterial pressure, respiratory rate, oxygen saturation, creatinine, blood urea nitrogen, sodium, bicarbonate, and hemoglobin. Physiologically implausible values were excluded with prespecified cleaning ranges: heart rate 20-250 beats/min, mean arterial pressure 20-200 mmHg, respiratory rate 4-80 breaths/min, oxygen saturation 50-100%, creatinine 0.1-25 mg/dL, blood urea nitrogen 1-250 mg/dL, sodium 95-180 mmol/L, bicarbonate 3-60 mmol/L, and hemoglobin 2-25 g/dL.

Comparator analyses in the strict sensitivity cohort used the same 11 SAKARYA predictors, the first-day SOFA score, and a simple renal-hemodynamic comparator that included blood urea nitrogen and systolic blood pressure.

2.4. Outcome Definition

The primary outcome was a composite cardiorenal endpoint defined as 28-day death or incident creatinine-based acute kidney injury. Acute kidney injury was defined with KDIGO serum creatinine criteria, using a baseline creatinine obtained from 24 hours before to 6 hours after intensive care unit admission and follow-up creatinine values through day 7 [14]. Urine output was not incorporated into the primary endpoint. Secondary outcomes were 28-day death and incident creatinine-based acute kidney injury. In the strict sensitivity cohort, baseline and predictor creatinine windows were separated as described above to reduce circularity.

2.5. Reporting Framework

The study was designed and reported in accordance with STROBE and RECORD recommendations for observational research using routinely collected data [15,16]. Prediction-model development and validation were reported according to TRIPOD [17].

2.6. Statistical Analysis

Continuous variables are reported as median [interquartile range] and categorical variables as number (percentage). Between-group comparisons used the Mann-Whitney U test or the chi-square test, as appropriate. Logistic regression was used for the primary model and for secondary outcome models. Discrimination was summarized with the area under the receiver operating characteristic curve, calibration with calibration intercept and slope, and global prediction accuracy with the Brier score.

Collinearity was assessed with variance inflation factors, with prespecified concern at values of 3.0 or higher. Restricted cubic spline screening with 4 degrees of freedom was used to test each continuous predictor for nonlinearity in the derivation cohort. The primary specification retained prespecified linear terms because the intended model was an interpretable early bedside model. Internal optimism was evaluated with 1000 bootstrap resamples. Additional sensitivity analyses included an extended complete-case model, multiple imputation by chained equations with five imputations, and prespecified subgroup interaction analyses.

Comparator analyses were performed only in the strict sensitivity cohort with complete availability of SAKARYA, SOFA, and BUN+SBP variables. SAKARYA, SOFA, and BUN+SBP models were refit in strict derivation and tested in strict temporal validation. Comparator performance included area under the receiver operating characteristic curve, Brier score, Brier skill score, calibration intercept, and calibration slope. Validation risk quartiles were created from the strict temporal SAKARYA predictions. Decision curve analysis was examined across threshold probabilities from 0.05 to 0.50. Bootstrap differences in validation area under the curve between SAKARYA and each comparator were estimated with 2000 resamples. All statistical analyses were performed using SPSS version 26.0 (IBM Corp., Armonk, NY, USA) and R software version 4.3.1 (R Foundation for Statistical Computing, Vienna, Austria). Specifically, decision curve analysis and bootstrap validation procedures were implemented using the ‘dca’ and ‘rms’ packages in R.

3. Results

Figure 1 summarizes cohort assembly. Among 94,458 adult intensive care unit stays, 16,641 first admissions carried a heart failure diagnosis. The broad acute heart failure enrichment strategy yielded 8477 patients, and the final clean core complete-case modeling cohort included 7065 patients. The derivation cohort comprised 5399 patients and the temporal validation cohort 1666 patients. The primary composite endpoint occurred in 2289 derivation patients and 810 temporal validation patients. A narrow phenotype cohort included 5618 patients, and the downstream strict sensitivity cohort with narrowed phenotype timing and separated creatinine windows included 737 derivation patients and 296 temporal validation patients.

Baseline characteristics of the broad clean core cohort are shown in Table 1. The temporal validation cohort was younger than the derivation cohort and had higher mean arterial pressure, higher respiratory rate, lower bicarbonate, more frequent early BNP/NT-proBNP testing, and less frequent early intravenous loop diuretic treatment. Sex distribution, heart rate, blood urea nitrogen, and hemoglobin were similar across the two time periods.

Primary model performance is shown in Table 2, with receiver operating characteristic curves in Figure 2 and calibration plots in Figure 3. For the primary composite outcome, the area under the receiver operating characteristic curve was 0.636 in derivation and 0.648 in temporal validation. The corresponding Brier scores were 0.231 and 0.238. Temporal validation calibration intercept and slope were 0.26 and 0.99. Bootstrap internal validation with 1000 resamples yielded an optimism-corrected area under the curve of 0.631 and an optimism-corrected calibration slope of 0.96. Secondary outcome performance is summarized in Table 2 and detailed in Supplementary Tables S6A and S6B.

Primary model coefficients rescaled to clinically interpretable units are shown in Table 3. Older age, higher heart rate, higher respiratory rate, higher creatinine, higher blood urea nitrogen, and higher sodium were associated with higher odds of the composite outcome, whereas higher mean arterial pressure, higher bicarbonate, and higher hemoglobin were associated with lower odds. The full one-unit model specification for direct probability calculation is provided in Supplementary Table S6C.

Outcome-stratified baseline characteristics are summarized in Supplementary Table S1. Collinearity was limited, with all variance inflation factors below 2.0 (Supplementary Table S2A). Restricted cubic spline screening showed departures from linearity in seven of ten continuous predictors, with the strongest evidence for creatinine and bicarbonate (Supplementary Table S2B). An RCS-augmented sensitivity model that applied spline terms to creatinine and bicarbonate produced only a small change in validation performance, with area under the curve rising from 0.648 to 0.653 and validation Brier score changing from 0.238 to 0.236 (Supplementary Table S2C).

Extended complete-case sensitivity results are shown in Supplementary Tables S3A and S3B. Compared with the core model, the extended model changed validation discrimination only slightly. Multiple-imputation sensitivity results are provided in Supplementary Tables S4A and S4B. Prespecified subgroup interaction analyses are shown in Supplementary Table S5 and indicated a lower creatinine-associated odds gradient in patients with higher baseline creatinine, whereas other interactions were small.

In the strict sensitivity cohort, derivation discrimination was 0.653 for SAKARYA, 0.690 for SOFA, and 0.616 for BUN+SBP. In strict temporal validation, SAKARYA achieved an area under the curve of 0.630, compared with 0.682 for SOFA and 0.607 for BUN+SBP (Table 4; Supplementary Table S7A). Validation Brier skill score was 0.006 for SAKARYA, 0.073 for SOFA, and -0.018 for BUN+SBP. Bootstrap differences in validation area under the curve were -0.052 (95% CI -0.129 to 0.023) for SAKARYA minus SOFA and 0.024 (95% CI -0.030 to 0.075) for SAKARYA minus BUN+SBP (Supplementary Table S7B). Supplementary Figure S1 shows the strict-cohort validation receiver operating characteristic curves.

In strict temporal validation, SAKARYA quartile event rates were 36.5% in Q1, 66.2% in Q2, 60.8% in Q3, and 71.6% in Q4 (Supplementary Table S7C). The Q2 and Q3 confidence intervals overlapped substantially, at 54.9-76.0% and 49.4-71.1%, respectively. Supplementary Figure S2 shows decision curve analysis. Net benefit curves overlapped at low thresholds, SOFA was numerically higher between thresholds of 0.20 and 0.30, and SAKARYA was numerically higher between thresholds of 0.35 and 0.40 without consistent dominance across the full threshold range. This minor non-monotonicity reflects the restricted sample size of the strict sensitivity validation subset (n=74 per quartile), with overlapping 95% confidence intervals demonstrating statistical equivalence in the intermediate risk range.

The primary composite outcome was 28-day death or incident creatinine-based acute kidney injury. Bootstrap optimism-corrected AUC for the primary model was 0.631 based on 1000 resamples, and the optimism-corrected calibration slope was 0.96.

Odds ratios are scaled for clinical interpretability. The full one-unit regression specification is provided in Supplementary Table S6C.

The strict comparator cohort included 296 temporal validation patients with 174 composite cardiorenal events. Full derivation and validation results are shown in Supplementary Table S7A.

4. Discussion

The present study developed and temporally validated an early cardiorenal risk model for critically ill adults with an operational acute heart failure phenotype. The broad primary model showed modest discrimination and preserved temporal calibration. The strict sensitivity analyses addressed the main design concerns, but performance remained modest and did not exceed SOFA. The data therefore support a transparent bedside description of early risk rather than a stand-alone replacement for established severity assessment.

Most established acute heart failure risk tools were not built in intensive care populations and were generally designed for short-term mortality rather than combined cardiorenal deterioration. External prospective work on EHMRG improved confidence in emergency-department risk estimation, but its transportability to an intensive care unit phenotype remains limited. The present results fit that gap. The model was derived in a cohort with higher physiologic severity, broader organ dysfunction, and a more heterogeneous clinical trajectory than the populations used for conventional acute heart failure scores.

Recent electronic health record and machine-learning studies in critically ill heart failure cohorts reported higher discrimination for mortality prediction than the composite performance observed here [18,19]. Another intensive care unit heart failure model also reported stronger mortality discrimination than the present composite endpoint [20]. Several differences likely explain this contrast. The present target was a mixed cardiorenal endpoint rather than mortality alone, the phenotype was operational rather than adjudicated, and the derivation strategy favored a fixed bedside variable set and transportable linear specification over algorithmic complexity.

The biological direction of the retained predictors was coherent. Older age, higher heart rate, higher blood urea nitrogen, and higher creatinine marked higher risk, whereas higher bicarbonate and hemoglobin marked lower risk. Those patterns agree with prior work linking renal dysfunction, metabolic derangement, and systemic illness severity to adverse acute heart failure outcomes [21,22]. The coefficients were therefore clinically plausible even though discriminative performance remained limited.

Two revision targets were especially informative. First, the strict cohort narrowed the phenotype window and separated baseline from predictor creatinine. That analysis reduced the main look-ahead and circularity concerns, but it also reduced sample size and left wide uncertainty around comparator differences. Second, restricted cubic spline screening confirmed nonlinearity in multiple predictors, especially creatinine and bicarbonate. An RCS-augmented sensitivity model produced only a small gain in validation discrimination and a small change in Brier score, so the linear model was retained for bedside interpretability. This tradeoff differs from work that pursued more flexible algorithms or artificial-intelligence-based prediction strategies for heart failure outcomes [23,24].

The comparator analyses clarified the clinical position of the model. In the strict cohort, SOFA showed higher validation discrimination and better Brier skill than SAKARYA, while BUN+SBP performed slightly worse. Decision curves showed threshold-dependent overlap rather than consistent superiority by any model. The risk quartiles also reflected limited separation in the middle range, because Q2 and Q3 had overlapping confidence intervals and reversed observed event order. These findings are consistent with broader reviews showing that many heart failure prognostic models lose separation when case mix, endpoint definition, and care setting change [25,26]. Importantly, the higher discrimination achieved by the SOFA score (AUC 0.682) must be interpreted in light of its temporal definition. SOFA represents a cumulative 24-hour assessment of multi-organ failure, incorporating parameters such as the Glasgow Coma Scale, platelet counts, and bilirubin levels, which are frequently unavailable or delayed during the initial hours of critical care. In contrast, the SAKARYA model operates entirely on a 6-hour horizon using readily transportable bedside parameters, offering an early triage window that SOFA cannot address. Furthermore, the overlapping event rates between intermediate risk quartiles in the strict temporal validation cohort highlight the model’s limitations in distinguishing mid-range risk within small sample sizes, pointing to the necessity of larger validation registries to verify calibration across all deciles.

The study has limitations. Acute heart failure status was operationalized rather than clinically adjudicated because structured bedside adjudication is not available in MIMIC-IV. The broad primary phenotype still used a 24-hour enrichment window, so the strict cohort is essential for interpretation rather than optional. The strict comparator analysis was small, with only 296 validation patients and 74 patients per quartile. SOFA was derived from first-day data rather than the same 6-hour horizon, which may have favored SOFA in head-to-head comparison. The database came from a single center, and temporal validation remained internal to that center. The strengths were prespecified variable selection, transparent coefficient reporting, temporal validation, explicit sensitivity analyses for look-ahead bias and creatinine circularity, bootstrap correction, and direct comparison with clinically familiar benchmarks. The next step is external multicenter validation in cohorts that define the phenotype at admission and preserve complete separation between baseline renal function and predictor creatinine.

5. Conclusions

The SAKARYA model was derived and temporally validated as a parsimonious bedside tool for early cardiorenal risk stratification in critically ill adults with an acute heart failure phenotype. The model demonstrated modest discrimination with preserved temporal calibration across both broad and strict validation cohorts. In head-to-head validation, its performance was comparable to established clinical benchmarks but did not exceed the cumulative 24-hour SOFA score. The findings suggest that early cardiorenal risk can be estimated using readily transportable clinical variables from the first 6 hours of ICU admission, providing a useful bedside description of early risk rather than a replacement for comprehensive organ failure assessments. Future external multicenter validation and clinical phenotype refinement are necessary to confirm these results and optimize risk stratification before clinical implementation.