This was a retrospective cohort study using the MIMIC-IV v2.2 database (ICU admissions at Beth Israel Deaconess Medical Center, 2008–2019). MIMIC contains de-identified health records maintained by the Laboratory for Computational Physiology at MIT. Institutional review board approvals were in place at MIT and BIDMC, and informed consent was waived owing to de-identification. One author completed the required training and performed all data extractions. The reporting follows the STROBE guidelines.

We initially identified 65,366 unique patients at their first ICU admission (≥18 years). Eligibility required availability of HbA1c and an admission glucose (earliest within 12h of ICU entry), with both recorded before any AF ascertainment. We then excluded patients with (1) an ICU length of stay <2 days (n = 34,142) or >28 days (n = 467); (2) a documented history of atrial fibrillation (AF) or atrial flutter (AFL) prior to ICU admission (n = 263); (3) missing either admission glucose or HbA1c (n = 26,470); (4) missing admission glucose or HbA1c data obtained after NOAF onset (n = 150); (5) absence of any heart rhythm records during the ICU stay (n = 10); and (6) missing key demographics (age, sex, or race). After these exclusions, 3,882 patients remained for analysis. A detailed selection flowchart is shown in Figure 1 .

Data were extracted via Structured Query Language (SQL) with PostgreSQL (version 16.0). Unless otherwise specified, all baseline variables were captured within the first 24 hours of ICU admission and prior to any ascertainment of NOAF.

The primary exposure variables used were the hemoglobin glycation index (HGI) and the stress/hyperglycemia ratio (SHR). Admission glucose was defined as the first plasma glucose measurement within 12 hours of ICU admission. The HGI was calculated as follows: observed HbA1c − predicted HbA1c, where predicted HbA1c = (0.009 × admission glucose [mg/dL]) + 4.940. The SHR was calculated via consistent units as admission glucose [mg/dL]/([28.7 × HbA1c %] − 46.7). Extreme or biologically implausible SHR values (<0.1 or >15) were reviewed and excluded.

Outcome Variable: The primary outcome was incident NOAF within 7 days of ICU admission. The event time was defined as the time of the first documented episode of AF, ascertained from bedside rhythm charting or telemetry records. Patients who did not develop NOAF were censored at the time of ICU discharge, in-hospital death, or on day 7, whichever occurred first.

Covariates: We collected baseline data on covariates that were included in the fully adjusted model (Model 3): (1) demographics (age, sex, race, BMI); (2) major comorbidities (myocardial infarction, congestive heart failure, diabetes, renal disease, liver disease, prior cardiac surgery); (3) illness severity scores (SOFA, SAPS II) and acute conditions (acute kidney injury [AKI], delirium); (4) vital signs (heart rate, mean arterial pressure [MAP], respiratory rate, SpO₂); (5) key laboratory values (white blood cell count, red blood cell count, platelet count, electrolytes, creatinine); and (6) ICU interventions (mechanical ventilation, vasopressor use, CRRT). Continuous variables were assessed for distribution and log-transformed when appropriate.

Time window and temporal ordering. All severity scores were computed from the worst recorded values within the first 24 hours after ICU admission and prior to any ascertainment of new-onset atrial fibrillation (NOAF) to minimize immortal-time and reverse-causation bias.

The Sequential Organ Failure Assessment (SOFA) score was derived across six organ systems—respiratory (PaO₂/FiO₂ with ventilatory support considered), coagulation (platelet count), liver (total bilirubin), cardiovascular (mean arterial pressure and vasopressor dose thresholds per the original SOFA definition), central nervous system (Glasgow Coma Scale, GCS), and renal (serum creatinine or urine output) systems. When preexisting organ dysfunction was unknown, the baseline value was assumed to be zero. For the CNS domain, we preferentially used presedation GCS when sedation information allowed a reliable estimate; otherwise, we used the recorded GCS without sedation adjustment and evaluated the impact in sensitivity analyses (see Supplementary Methods ). Vasopressor categories followed the original SOFA thresholds (dopamine, dobutamine, epinephrine, and norepinephrine) and used weight-normalized doses.

The simplified acute physiology score II (SAPS II) was computed from 17 variables captured in the first 24 hours (the worst value used for each physiological variable): 12 physiologic measurements (temperature, heart rate, systolic blood pressure, PaO₂ or alveolar–arterial oxygen gradient depending on FiO₂, 24-hour urine output, serum bicarbonate, total bilirubin, serum sodium, serum potassium, blood urea nitrogen, white blood cell count, and GCS), age, admission type (medical/scheduled surgical/unscheduled surgical), and three chronic conditions (AIDS, metastatic cancer, hematologic malignancy). Per the original specification, the oxygenation component used PaO₂ when FiO₂ < 0.5 and the A–a gradient when FiO₂ ≥ 0.5; we followed this rule in our implementation. As with the SOFA, we used presedation GCSs where feasible; otherwise, the recorded value was used, and robustness was assessed via sensitivity analyses.

The comorbidity burden was quantified via the Charlson Comorbidity Index mapped from ICD-9/ICD-10 diagnosis codes recorded during the index hospitalization via the Quan et al. coding algorithm; weights were summed to obtain the CCI. Because using only index-stay diagnoses can underascertain preexisting conditions, we flag this as a limitation and probed robustness in sensitivity analyses (e.g., restricting diagnoses flagged as present-on-admission [POA], where available).

Missing covariate data were handled via multiple imputation by chained equations (MICE; m=5; 10 iterations) under a missing−at−random assumption, including all variables in the fully adjusted model (Model 3). The exposure and outcome variables were not imputed. The estimates were pooled via Rubin’s rules.

For descriptive analyses, baseline characteristics were compared between patients with and without NOAF and across the HGI and SHR quartiles. Continuous variables, presented as medians with interquartile ranges (IQRs), were compared via the Mann–Whitney U test for two-group comparisons and the Kruskal–Wallis test for comparisons across quartiles. Categorical variables are presented as counts and percentages and were compared via the chi-square test or Fisher’s exact test, as appropriate.

The cumulative incidence of NOAF across quartiles was visualized with Kaplan–Meier curves and compared via the log-rank test. We then used Cox proportional hazards models to estimate hazard ratios (HRs) and 95% confidence intervals (CIs) for the associations between HGI/SHR quartiles and NOAF risk. A prespecified, progressive adjustment strategy was used: Model 1: Unadjusted. Model 2: Adjusted for age, sex, race, and BMI. Model 3 (Fully Adjusted): Further adjusted for SOFA and SAPS II scores, major comorbidities, key ICU interventions, vital signs, and laboratory parameters, as listed in the definitions section.

To assess nonlinear relationships, we modeled HGI and SHR as continuous variables via restricted cubic splines with three knots (at the 10th, 50th, and 90th percentiles). A likelihood ratio test was used to calculate a P value for nonlinearity. The proportional hazards assumption was checked via Schoenfeld residuals, and multicollinearity was assessed via variance inflation factors (VIFs). All the statistical analyses were prespecified and conducted via R software (version 4.4.0). A two-sided P value < 0.05 was considered statistically significant.

Prespecified subgroup analyses assessed the consistency of associations across age (≥65 vs <65 years), sex, BMI (≥30 vs <30 kg/m²), race, and history of myocardial infarction, congestive heart failure, diabetes, and prior cardiac surgery. Subgroup estimates were obtained from the fully adjusted Model 3, using the same exposure parameterizations (quartiles and, where applicable, continuous RCS). Effect modification was evaluated on the multiplicative scale by adding exposure × subgroup cross-product terms to Model 3; P for interaction was derived from likelihood-ratio tests comparing models with and without the interaction. The results are displayed with forest plots. The primary sensitivity analysis prespecified stratification by diabetes status, repeating the full modeling framework within each stratum.

A total of 3882 patients were enrolled in the analysis. Baseline characteristics are shown in Table 1 . The median age was 68 years, and 2255 (58.1%) patients were male. Compared with those without NOAF, patients with NOAF were older and had a greater comorbidity burden (higher Charlson Comorbidity Index scores and a greater prevalence of congestive heart failure, myocardial infarction, and renal disease (all p < 0.01)). Greater severity (higher SAPS II and SOFA scores) led to higher requirements for intensive interventions (vasopressor use, mechanical ventilation, and continuous renal replacement therapy) in patients with NOAF (all p < 0.001). Vital signs were more unstable in patients with NOAF (significantly lower blood pressure). Key laboratory data revealed lower red blood cell counts, lower platelet counts, and metabolic alterations (lower glucose, calcium, and SHR levels (all p < 0.01)) ( Table 1 , Supplementary Table S1 ).

To investigate the relationship between glycemic variability and clinical outcomes, patients were stratified by HGI and SHR quartiles, revealing distinct association patterns. In the HGI group, Q4 was associated with greater chronic disease burden (higher CCI, diabetes incidence), whereas Q3 presented the highest acute illness severity (SAPS II) ( Supplementary Table S2 ). Interestingly, in the SHR group, Q1 had the most severe acute conditions, with the highest SAPS II and SOFA scores, greater use of CRRT and vasopressors, and more cardiovascular complications ( Supplementary Table S3 ).

During the seven-day follow-up period, 750 patients (19.3%) presented with NOAF. A Kaplan–Meier analysis revealed significant disparities in the cumulative incidence of NOAF across the quartiles of both the HGI and the SHR (log-rank test, both P < 0.001; Figure 2 ). For HGI, the highest cumulative incidence was found in the third quartile (Q3), followed by the second quartile (Q2), suggesting a nonlinear, inverted U-shaped relationship. In contrast, SHR demonstrated a consistent linear inverse pattern, such that the lowest quartile (Q1) was associated with the highest risk of NOAF, with a progressive linear decrease in incidence across quartiles.

To further investigate the relationships among HGI, SHR, and the risk of NOAF, we constructed multivariate Cox proportional hazards models. The results of the associations between quartiles of HGI and SHR with NOAF are presented in Table 2 .

For HGI, in the unadjusted (Model 1) and partially adjusted (Model 2) models, both the second (Q2) and third (Q3) quartiles were significantly associated with an increased risk of NOAF compared with the reference group (Q1). In the fully adjusted Model 3, HGI Q2 (HR, 1.26; 95% CI, 1.02–1.57; P = 0.034) and Q3 (HR, 1.36; 95% CI, 1.09–1.70; P = 0.006) remained independent risk factors for NOAF. Notably, the risk associated with the highest quartile (Q4) was not statistically significant after full multivariable adjustment (HR, 1.28; 95% CI, 0.95–1.72; P = 0.102).

For SHRs, a negative association with the risk of NOAF was consistently observed. In both the unadjusted (Model 1) and partially adjusted (Model 2) models, the third (Q3) and fourth (Q4) quartiles of SHR were associated with a significantly lower risk of NOAF than the lowest quartile (Q1). This association remained robust even after full adjustment in Model 3, where Q3 (HR, 0.76; 95% CI, 0.62–0.93; P = 0.007) and Q4 (HR, 0.69; 95% CI, 0.55–0.85; P < 0.001) were still linked to a significantly reduced risk of developing NOAF.

Using restricted cubic spline (RCS) models, we evaluated the continuous associations between HGI, SHR, and the risk of NOAF, adjusting for the predefined covariate set used in Model 3. As illustrated in Figure 3A , HGI showed a significant inverted U−shaped relationship with NOAF (P for nonlinearity < 0.001), with risk peaking at mid−range values and attenuating toward both tails. As shown in Figure 3B , SHR displayed a predominantly linear inverse association (P for nonlinearity = 0.236), indicating a progressively lower NOAF risk at higher SHR levels. These spline patterns are concordant with the quartile analyses (higher risk in HGI Q2–Q3 vs Q1; lower risk across higher SHR quartiles) and avoid reliance on arbitrary cutoff points.

We performed prespecified risk subgroup analyses and multiplicative interaction tests according to routinely reported strata: age (e.g., <65 vs. ≥65 years), sex, BMI (e.g., <30 vs. ≥30 kg/m²), race (white, black, and other), and clinical history (MI, CHF, DM, prior cardiac surgery), using the same covariate set as Model 3. Figure 4 presents the adjusted hazard ratios with 95% CIs for each stratum and the corresponding P values for interaction.

For HGI, the association with NOAF was directionally consistent across strata, with the excess risk concentrated in mid−range categories (Q2–Q3 vs Q1) in multiple groups (e.g., females, individuals aged ≥65 years, and patients without prior MI/CHF/DM). When sample sizes were smaller, CIs widened, and statistical significance varied, but the pattern of a mid−range elevation was preserved.

For SHRs, higher quartiles were generally associated with lower NOAF risk, mirroring the main analysis. The inverse gradient was most evident in clinically common strata (e.g., older adults, White patients, and those without MI/CHF/DM or with prior cardiac surgery), whereas some subgroups presented wider CIs and non−significant estimates.

Crucially, no effect modification was detected—all P values for interaction > 0.05 for both HGI and SHR—indicating that the observed subgroup differences are descriptive rather than confirmatory. Accordingly, the pooled estimates provide the most reliable summary of the associations, and subgroup findings should be interpreted with caution given the varying precision across strata.

We observed distinct patterns according to diabetes status ( Supplementary Figure S1 , Supplementary Table S3 ). Among nondiabetic patients, higher HGI quartiles were associated with increased NOAF risk after full adjustment (Q3 vs Q1: HR 1.37, 95% CI 1.04–1.81, P = 0.024; Q4 vs Q1: HR 1.37, 95% CI 1.04–1.81, P = 0.024). Conversely, a higher SHR was protective (Q3 vs. Q1: HR 0.75, 95% CI 0.58–0.96, P = 0.024; Q4 vs. Q1: HR 0.67, 95% CI 0.51–0.87, P = 0.003). Among diabetic patients, HGI quartile was not significantly associated with NOAF; for SHRs, only Q4 was associated with a lower risk (HR 0.64, 95% CI 0.43–0.94; P = 0.024). Where available, P for the exposure × diabetes interaction is shown in Figure 4 (numerical values in Supplementary Table S3 ).