Staphylococcus aureus (S. aureus) is a major human pathogen responsible for a wide spectrum of diseases, from superficial skin infections to life-threatening sepsis (Maurer et al., 2015). It is the leading cause of fatal bloodstream infections (BSIs), with mortality rates ranging from 20% to 40%, posing a significant public health burden worldwide (Claes et al., 2014; Cooper et al., 2020; Krasaewes et al., 2022). In intensive care unit (ICU) patients, S. aureus bloodstream infection (SA-BSI) is particularly challenging due to host immunosuppression, frequent invasive procedures, and multiple comorbidities. These factors contribute to prolonged length of stay (LOS) in hospital and ICU, elevated resource utilization, and markedly higher mortality compared to general ward populations (Vincent et al., 2009). While substantial research has centered on Gram-negative bacteremia, especially in regions like China, SA-BSI remains comparatively under-investigated, particularly regarding its dynamic clinical evolution and risk stratification in ICU settings. As the global population ages and antimicrobial resistance escalates, early recognition and individualized management of SA-BSI are increasingly vital to improving outcomes in this vulnerable population.

The Sequential Organ Failure Assessment (SOFA) score quantifies dysfunction across six organ systems—respiratory, coagulation, hepatic, renal, cardiovascular, and central nervous system (CNS)—on a 0–4 scale, with higher scores indicating more severe organ failure (Vincent et al., 1996). It has long been used as a prognostic tool in critically ill patients (Vincent et al., 1998). However, conventional single-time-point SOFA measurements offer only static snapshots of a patient’s condition and may fail to reflect the dynamic and often nonlinear progression of organ dysfunction (OD) in the ICU (Yang et al., 2022). Given the rapid physiological shifts common in critical illness, a longitudinal assessment of SOFA trajectories could provide deeper insight into disease evolution, capturing both subclinical deterioration and unexpected recovery. These temporal trends are particularly valuable for early risk stratification and timely therapeutic decision-making. In the context of SA-BSI, where delays in diagnosis and intervention can be fatal, tracking OD progression over time may offer critical prognostic information and inform personalized management strategies.

Group-based trajectory modeling (GBTM), originally proposed by Nagin (2005), is a statistical method for identifying latent subgroups within a population that share similar longitudinal trends. While first introduced in criminology, GBTM has gained traction in clinical research, where time-series data, such as serial SOFA scores, are increasingly available (Nagin and Odgers, 2010). By enabling data-driven phenotyping based on the evolution of clinical parameters, GBTM offers a promising framework for uncovering disease heterogeneity in the ICU (Wang et al., 2024).

Previous studies have applied GBTM or similar approaches to analyze SOFA trajectories in patients with sepsis (Xu et al., 2022), COVID-19 (Yu et al., 2024), and other critical illnesses, and to investigate temporal trends of inflammatory markers (e.g., C-reactive protein [CRP]) and vital signs in BSI (Gu et al., 2024). However, no prior research has specifically focused on longitudinal SOFA score trajectories in ICU patients with SA-BSI. This represents a critical knowledge gap, particularly given the dynamic and multisystem nature of OD in this population. Notably, the early phase of BSI, which is before microbiological confirmation, is a high-stakes window in which empirical treatment decisions must be made. During this period, early signs of multi-organ dysfunction (MOD) may be subtle or misinterpreted, and prognostically relevant trends may be overlooked. Capturing the temporal evolution of SOFA scores through GBTM may help identify distinct clinical courses that are not evident from initial assessments alone.

In addition, the early identification of such trajectory patterns at ICU admission remains challenging. Integrating machine learning (ML) methods with routinely available clinical data may offer a practical approach to support early trajectory stratification. Furthermore, interpretable ML techniques such as SHapley Additive exPlanations (SHAP) can facilitate clinical understanding by highlighting key baseline features associated with different disease courses.

Accordingly, the objectives of this study were to (1) identify distinct SOFA score trajectory groups in ICU patients with SA-BSI using GBTM, and (2) develop and externally validate an interpretable ML model to predict trajectory group membership using only ICU admission data. Through this two-pronged approach, we aimed to deepen understanding of the dynamic progression of OD in SA-BSI and explore its clinical implications for early phenotyping and individualized ICU management.

This study identified several novel and clinically relevant findings. First, by analyzing longitudinal SOFA score trajectories over the early course of SA-BSI, we uncovered three distinct OD patterns that were consistently reproducible across two independent ICU cohorts and strongly associated with clinical outcomes. Higher-risk trajectory groups showed significantly increased mortality, underscoring the prognostic value of dynamic OD assessment. Second, biological markers and SOFA subscores demonstrated progressive deterioration from Group 1 to Group 3, with cardiovascular dysfunction emerging as the predominant contributor to adverse outcomes. Third, in the SAHZU cohort, the availability of inflammatory and coagulation biomarkers (e.g., CRP, PCT, D-dimer) provided additional insight into systemic responses not captured by SOFA alone, highlighting the potential of multidimensional risk profiling. Finally, we demonstrated that ML using early ICU data can accurately predict trajectory group membership. The incorporation of SHAP-based interpretability further revealed physiologically coherent predictors, such as preserved renal function for favorable courses and elevated lactate and acidosis for severe deterioration, supporting both the transparency and clinical relevance of this approach. Together, these findings suggest that combining trajectory modeling with interpretable ML offers a feasible and informative strategy for early risk stratification in SA-BSI.

While prior studies have demonstrated the prognostic value of short-term (e.g., 72-hour) SOFA score changes in sepsis populations (Minne et al., 2008; Xu et al., 2022), our findings extend this understanding by emphasizing the importance of sustained longitudinal trends. By capturing dynamic patterns of OD over a five-day window, our trajectory-based approach identified heterogeneous clinical courses that may be overlooked by static or early SOFA measurements. Rather than modeling the entire ICU stay, this focused early window was intentionally selected to reflect the clinically actionable phase surrounding SA-BSI onset, during which OD evolves rapidly and key management decisions are often required before definitive microbiological confirmation. This is consistent with previous research in other ICU populations (Zhou et al., 2020; Yu et al., 2024), including our own study in critically ill COVID-19 patients, where four distinct SOFA trajectories were identified and the most severe group exhibited markedly higher 28-day mortality (Yu et al., 2024). In SA-BSI, we observed a similar pattern: Group 3, defined by persistently elevated SOFA scores, was strongly predictive of adverse outcomes. An additional insight from our trajectory analysis was the divergent SOFA progression patterns between Groups 2 and 3. Group 2 exhibited a steady increase in SOFA scores, whereas Group 3 peaked early and then declined slightly while remaining critically elevated. Despite this apparent improvement, Group 3 had significantly worse outcomes, suggesting that short-term decreases in SOFA scores do not necessarily indicate recovery and may instead mask persistent OD or ongoing systemic inflammation—key drivers of mortality in SA-BSI (Minne et al., 2008; Zorowitz, 2016). These findings highlight the limitations of relying solely on static or short-interval SOFA changes for prognostication. It should also be noted that the intermediate trajectory group likely encompasses patients with heterogeneous clinical courses, including individuals who subsequently stabilize as well as those who progress toward more severe OD, which may partly explain its intermediate prognostic profile. Notably, these trajectories were not only descriptive but also predictable. Using an XGBoost classifier trained on ICU admission data, we accurately assigned patients to their eventual trajectory groups. SHAP interpretation further enhanced model utility by providing biologically coherent explanations. For example, preserved renal and hepatic function favoring stable courses, while elevated lactate and low pH predicted deterioration. By identifying high-risk phenotypes shortly after blood culture collection (before definitive microbiological results become available), trajectory-informed ML models may support timely triage, guide early intervention, and optimize empiric antimicrobial therapy during this critical diagnostic window. In this context, the SOFA trajectory observation window (Days −1 to +3) was designed to align with the early diagnostic phase of SA-BSI, when clinical decisions frequently need to be made in the absence of definitive microbiological confirmation.

Cardiovascular dysfunction emerged as the dominant contributor to poor outcomes, particularly in the highest-risk trajectory group. In Group 3, cardiovascular SOFA subscores remained persistently elevated, indicating sustained circulatory failure closely linked to increased mortality. This finding echoes prior evidence that cardiovascular compromise plays a central role in ICU sepsis mortality (Snell and Parrillo, 1991; Mayr et al., 2006), contributing to impaired perfusion, metabolic derangements, and systemic inflammation. In the context of SA-BSI, this pattern may be further driven by pathogen-specific mechanisms, including toxin-mediated endothelial injury, superantigen-induced hyperinflammatory responses, and sepsis-associated myocardial depression, all of which can exacerbate vasoplegia and circulatory collapse (Powers et al., 2012; Touaitia et al., 2025). Interestingly, this contrasts with our previous work in critically ill COVID-19 patients (Yu et al., 2024), where respiratory dysfunction was the predominant driver of poor prognosis. Such divergence highlights the disease-specific nature of OD evolution and underscores the need for tailored risk assessment strategies in different syndromes. By decomposing SOFA trajectories into organ-specific components, our approach provides granular insight into the pathophysiology of SA-BSI and may better guide resource allocation, such as early circulatory support, in patients at highest risk.

The availability of inflammatory and coagulation biomarkers in the SAHZU cohort offered valuable complementary insights into systemic responses underlying OD progression. In Group 3, markers such as CRP, PCT, and D-dimer were markedly elevated, consistent with a hyperinflammatory and hypercoagulable state that paralleled severe OD. These findings reinforce the well-established role of systemic inflammation and coagulopathy in driving poor outcomes in critical illness (Strauss et al., 2001; Brummel et al., 2024), and support a multidimensional risk stratification approach that integrates both physiological and biochemical markers. Notably, although both Groups 2 and 3 had elevated mortality in both cohorts, only the MIMIC cohort showed a statistically significant survival difference between Groups 1 and 2 (Supplementary Table S6). This discrepancy may reflect sample size limitations or unmeasured cohort-specific factors, but also suggests that the additional biomarker data available in SAHZU allowed finer risk discrimination, particularly among patients with similar SOFA trajectories. Together, these observations highlight the added value of combining trajectory-based modeling with dynamic biomarkers to more comprehensively capture the heterogeneity and severity of SA-BSI. Given that these inflammatory and coagulation biomarkers were available only in the SAHZU cohort, their associations with trajectory patterns were explored in a descriptive and hypothesis-generating manner, and should be interpreted cautiously pending validation in larger, multi-cohort studies.

Beyond trajectory characterization, we further examined the relationship between baseline severity and longitudinal risk patterns. In multivariable Cox regression analyses, baseline SOFA score was no longer independently associated with mortality after adjustment for trajectory group membership. This finding likely reflects the fact that trajectory classification captures longitudinal patterns of OD, integrating both baseline severity and subsequent disease evolution over time. Accordingly, trajectory group membership provides prognostic information beyond that of a single time-point SOFA measurement, rather than diminishing the clinical relevance of baseline SOFA itself. Building on this, our ML model further demonstrated that early ICU data could accurately predict patients’ subsequent OD trajectories. The XGBoost classifier, trained solely on features available within the first 24 hours of admission, achieved excellent discrimination across all trajectory groups, including in external validation. This underscores that early physiological and laboratory indicators already carry prognostic signals reflecting the likely clinical course in SA-BSI. Importantly, SHAP-based model interpretation enhanced clinical transparency by identifying trajectory-specific predictors. Preserved renal and hepatic function characterized stable courses, whereas elevated lactate, low arterial pH, and high APACHE II scores were strong predictors of deterioration. Unlike conventional risk models focused on static outcomes (e.g., mortality), this approach anticipates dynamic disease evolution, enabling earlier recognition of patients at risk for progressive OD. As such, interpretable ML offers a promising tool for real-time triage and personalized management during the diagnostic window when culture results remain pending.

This study focuses on the early-stage SOFA score trajectories in patients with SA-BSI, aiming to integrate trajectory-based assessment into clinical practice during the period when microbiological results are still pending. This approach offers several potential benefits. First, it provides a framework for continuous and dynamic evaluation of OD, enabling clinicians to capture evolving physiological changes rather than relying solely on static SOFA values. Second, trajectory group classification serves as a practical tool for early ICU risk stratification. During this early clinical period, trajectory-based information should be interpreted as a dynamic indicator of evolving risk rather than as a definitive prognostic endpoint. Trajectories trending toward higher-risk patterns may signal the need for closer monitoring and heightened clinical vigilance before microbiological confirmation, whereas stable trajectories may suggest a lower immediate risk of rapid deterioration. By identifying patients most likely to deteriorate, this method may help optimize resource allocation and guide the timing and intensity of empiric antimicrobial therapy, particularly during the diagnostic window before culture confirmation. For instance, patients classified in Group 3 may benefit from early escalation of cardiovascular support and broad-spectrum antibiotics, whereas Group 2 patients may require closer hemodynamic monitoring and tailored interventions to prevent further deterioration. In contrast, Group 1 patients with stable OD may be managed less intensively, thereby facilitating efficient ICU triage. Third, trajectory-based modeling holds promise for improving clinical trial design. Stratifying patients by trajectory-defined risk can enhance cohort homogeneity, enabling interventions to be evaluated within more phenotypically consistent subgroups and thereby improving both statistical power and translational relevance in infectious disease trials. Finally, integrating automated trajectory calculations into electronic health records (EHRs) could enable real-time decision support for clinicians. By flagging patients who begin to follow high-risk trajectories, such systems may prompt timely interventions and more personalized treatment strategies, including targeted antimicrobial stewardship. Collectively, these applications align with the broader goals of dynamic infection surveillance and precision critical care in the ICU. Importantly, the trajectory-based framework is intended to support early risk stratification and clinical decision-making, rather than to replace clinician judgment or to serve as direct guidance for specific therapeutic interventions. Focusing specifically on SA-BSI provides additional insight beyond sepsis in general, as SA-BSI represents a distinct clinical entity characterized by unique virulence factors, a propensity for persistent bacteremia, and a high burden of cardiovascular complications. By examining SOFA trajectories within this pathogen-specific context, our study highlights patterns of OD evolution that may be diluted or obscured in more heterogeneous sepsis populations.

This study has several limitations. First, as a retrospective analysis, selection bias and missing data inevitably influenced our findings, particularly in the MIMIC-IV cohort, where detailed microbiological and laboratory information, such as inflammatory markers, was unavailable. These gaps reduced the granularity of trajectory characterization and limited direct comparability between cohorts. To mitigate missingness, we leveraged the inherent strengths of GBTM in handling incomplete longitudinal data via MLE and provided detailed data availability summaries (Supplementary Table S7) to enhance transparency and reproducibility. Second, differences between the MIMIC and SAHZU cohorts in demographics, comorbidities, ICU admission thresholds, clinical management, and biomarker availability, as well as differences in laboratory measurement frequency and data collection practices, may have contributed to variations in trajectory distribution and survival patterns. Although multivariable Cox regression adjusted for known confounders, residual confounding cannot be excluded. Nonetheless, the consistent three-group pattern and reproducible associations with mortality across cohorts support the robustness and cross-setting generalizability of our findings. Future studies should incorporate additional biomarkers, such as lactate, CRP, and pathogen-specific data, to refine trajectory-based phenotyping. Moreover, microbiological details, including resistance profiles and toxin gene expression, were unavailable; thus, we could not assess pathogen-driven modulation of OD or treatment response. More recently, SOFA-2 has been proposed to better reflect contemporary ICU practices (Ranzani et al., 2025). The present study was designed and conducted using the original SOFA score due to data availability and consistency across cohorts, which also facilitates comparison with prior trajectory-based studies. Future work will be needed to evaluate whether SOFA-2–based trajectories provide additional or incremental prognostic value in patients with SA-BSI. Despite these limitations, our study provides novel insight into the dynamic patterns of OD in SA-BSI and highlights the potential of trajectory-informed and interpretable ML approaches to enhance individualized risk assessment in critical care.

In conclusion, this study identified three distinct SOFA score trajectory groups in ICU patients with SA-BSI, each representing unique patterns of OD progression and associated with markedly different clinical outcomes. Higher-risk trajectories, particularly Group 3, were characterized by severe and sustained cardiovascular compromise, which emerged as a prominent prognostic contributor across cohorts. Trajectory groups were also associated with meaningful clinical and biological gradients, supporting their phenotypic relevance. Moreover, we demonstrated that ML models based on early ICU admission features can accurately predict OD trajectories, and SHAP-based interpretation provided clinically meaningful insights into key physiological drivers. Together, these findings highlight the potential of trajectory-informed phenotyping for early risk stratification in critically ill patients with SA-BSI, complementing conventional severity scoring systems. Further prospective, multicenter studies are warranted to validate these findings, and future integration with microbiological profiling may enhance infection-specific precision care and antimicrobial stewardship in ICU settings.