Abstract
Objective:
Early identification of mortality risk in critically ill children with suspected infection remains challenging. Whether severity scores (SIRS, pSOFA, Phoenix) add incremental value beyond data-derived physiologic features in machine-learning models is uncertain. We compared the predictive value of physiologic dynamics, treatment features, and continuous laboratory trends expert-derived severity representations for 24-hour PICU mortality prediction, and derived a parsimonious, interpretable model.
Study Design:
We analyzed 3310 PICU encounters (107 deaths, 3.2% mortality) using first-24-hour data. A staged, domain-structured feature selection strategy was used to build sequential XGBoost models evaluated via 5-fold stratified cross-validation. Performance was assessed by AUROC, AUPRC, Brier score, and specificity at 80% sensitivity, with SHAP-guided pruning at each stage.
Results:
Vital sign extrema provided baseline discrimination. Incorporation of variability and slope features improved calibration and specificity. Vasopressor and supplemental oxygen timing and intensity features produced substantial gains in AUROC and AUPRC. Continuous laboratory extrema and temporal trends reflecting evolving organ dysfunction yielded further improvement. Expert-threshold-derived metrics, including SIRS, pSOFA, and Phoenix composite scores, did not improve performance beyond data-derived features.
Conclusion:
Routinely collected hourly PICU data from the first-24-hour are sufficient for XGBoost to extract deterioration patterns strongly predictive of early mortality, without requiring severity thresholds.
Keywords: Physiologic data, Temporal resolution, Infection, Pediatric intensive care unit
Introduction
Early risk stratification for mortality in critically ill children with a suspected infection remains a central challenge in pediatric intensive care. Although overall sepsis-associated mortality in the pediatric intensive care unit (PICU) is relatively low, early identification of children at elevated risk within the first 24 hours of admission is clinically valuable for triage, escalation of care, and resource allocation [1]. Contemporary PICU environments generate dense, time-resolved electronic health record (EHR) data streams that include frequent vital signs, repeated laboratory assessments, and detailed treatment records [2]. These data create an opportunity for machine learning (ML) models to detect patterns in evolving physiologic instability not captured by single “worst-value” assessments.
Early identification and risk stratification in pediatric sepsis management has traditionally relied on expert-derived thresholds and composite severity scores. Established frameworks such as SIRS [3], pSOFA [4–6], and Phoenix [7,8] encode physiologic and laboratory assessments into interpretable summaries of illness severity and organ dysfunction, and are widely used for bedside decision support, quality benchmarking, and research. These scores commonly use age-adjusted reference ranges for vital signs and laboratory values to quantify physiologic derangement and mortality risk [9]. These constructs were developed to standardize assessment and reduce cognitive burden in settings such as Emergency Departments where data are sparse, irregular, or difficult to integrate in real time.
However, the availability of dense, high-frequency physiologic data fundamentally alters the prediction problem. Threshold-based representations discretize continuous measurements and may discard information about magnitude, trajectory, and response to therapy, features that modern ML models can exploit effectively when trained on temporally resolved data. Composite scores also impose fixed weighting schemes and organ system abstractions that may not align with patient history, local patient populations, practice patterns, or evolving clinical trajectories. They may also obscure short-term dynamics critical for early risk prediction.
Although these scoring systems share a reliance on worst-value (extrema) physiologic and laboratory abnormalities relative to age-adjusted norms, they differ in their incorporation of treatment variables. While SIRS is intervention-independent, pSOFA and Phoenix incorporate cardiovascular and respiratory therapies, specifically vasopressor use and mechanical ventilation as indicators of organ dysfunction and septic shock, alongside expert-defined thresholds on physiologic and biochemical variables [5,7].
Recent work has demonstrated that early use of vasopressor or invasive mechanical ventilation is strongly associated with in-hospital mortality in PICU patients with a suspected infection [10]. Another recent ML study predicting Phoenix criteria in PICU patients found pSOFA to be the most important predictive feature [8]. However, it remains unclear whether predictive performance attributed to pSOFA derives primarily from its embedded intervention content, expert threshold defined physiologic abnormalities, or their composite structure.
More broadly, the value of expert knowledge in ML-based clinical prediction remains an active area of investigation. Expert-augmented machine learning approaches suggest that incorporating domain knowledge in the form of iterative expert rankings of machine generated rules may improve predictive model performance, particularly under dataset shift or limited sample sizes [11,12]. However, the benefit of a priori expert knowledge was not evaluated here.
A further, often under-recognized challenge is that commonly used pediatric severity scores—including pSOFA, and Phoenix—incorporate neurologic components that are frequently unobserved or inconsistently documented in infection-driven PICU admissions without trauma or primary neurologic disease. Formal neurologic assessments such as the Glasgow Coma Scale (GCS) are often omitted in routine care when neurologic injury is not suspected, reflecting standard clinical workflow rather than data quality limitations. Although scores such as Phoenix are designed to permit calculation in the absence of documented neurologic assessments, missing neurologic data implicitly contribute no neurologic dysfunction points, potentially underrepresenting neurologic involvement in this population. Neurologic dysfunction may therefore remain latent or incompletely captured in EHR-derived datasets. This motivates evaluation of alternative feature representations—physiologic trajectories, laboratory dynamics, and treatment timing—that can capture early deterioration when formal neurologic scoring is unavailable or incomplete.
Here, we examined whether externally defined expert-knowledge features provide meaningful incremental predictive value for early mortality risk when paired with time-resolved PICU data and modern gradient-boosted tree models (XGBoost) [13]. Expert features considered include threshold-based abnormality flags, burden metrics capturing persistence and duration of vital-sign abnormality [14], and established risk-stratifying severity scores.
To address this question, we implemented a staged, domain-structured modeling strategy that incrementally introduces distinct classes of features while monitoring discrimination, calibration, and feature importance. Using data from the first 24 hours of PICU admission, we compared static extrema measures, dynamic physiological trajectories, treatment-based variables, continuous laboratory features, and expert-derived representations to determine when–and to what extent–expert knowledge augments the predictive performance of contemporary machine learning models for pediatric critical illness.
Methods
Study design and cohort:
We conducted a retrospective observational study of children aged >3 months and <18 years with suspected infection admitted to the PICU at a single tertiary, freestanding pediatric hospital during the period March 2022-February 2025. Suspected infection was defined by the presence of a blood culture order. A total of 3310 PICU encounters met inclusion criteria, of which 107 (3.2%) resulted in in-hospital mortality. Predictor variables were derived exclusively from data available within the first 24 hours of PICU admission. Mortality occurring within the initial 24-hour window was excluded from the prediction target to ensure temporal separation between feature extraction and outcome ascertainment. The primary outcome was in-hospital mortality, identified by either a documented death timestamp or discharge disposition of “expired.” All data were deidentified and collected under a Pearl IRB-approved protocol (22-CTAI-101)
Data Sources:
Time-stamped EHR data were extracted across multiple domains, including demographics, vital signs, laboratory results, microbiology, medications, and respiratory support. Vital signs were recorded at high temporal density (median >20 measurements per variable over 24 hours), and laboratory and treatment data were captured as recorded in routine clinical care.
Feature engineering and domain construction.
Candidate predictors were derived from data available within the first 24 hours of PICU admission and organized into clinically coherent feature domains representing progressively richer representations of patient physiology and clinical context:
Vital sign extrema, defined as minimum and maximum values over the 24-hour window, consistent with traditional severity score representations.
Vital sign dynamics, including variability (standard deviation) and linear slopes, to capture physiologic instability and short-term trajectories.
Physiologic burden metrics, defined as severity–duration summaries based on expert age-adjusted thresholds (e.g., cumulative time or area beyond abnormal ranges).
Treatment features, including timing and intensity of vasopressor therapy and supplemental oxygen, representing clinician response to physiologic instability.
Laboratory features, represented by continuous extrema and short-term trends.
Laboratory abnormality flags defined using expert age-adjusted normality thresholds.
Expert-derived severity constructs, including SIRS, pSOFA and Phoenix scores.
Preliminary statistical screening
Prior to model development, candidate features were compared between survivors and non-survivors using Mann–Whitney U tests for continuous variables and χ2 or Fisher’s exact tests for categorical variables. Features demonstrating minimal univariate separation were deprioritized but not excluded solely on statistical significance, as multivariable interactions were evaluated during model training. These univariate comparisons were used to support feature screening rather than as formal hypothesis tests and were therefore reported descriptively without adjustment for multiple comparisons.
Sequential model development
Model development followed a staged, domain-structured approach. An initial model based on vital sign extrema served as the baseline. At each subsequent stage, candidate features from a new domain were added to the best-performing model from the prior stage, allowing explicit assessment of incremental predictive value.
At each stage, feature retention was guided by changes in discrimination (ΔAUROC), precision–recall performance (AUPRC), calibration (Brier score), and specificity at a prespecified sensitivity of 0.80. This approach ensured that each feature domain was evaluated in the context of previously retained information rather than through independent re-optimization.
Feature pruning and parsimony
Within each stage, feature pruning was performed using a combination of quantitative and clinical criteria. Feature contributions were assessed using SHapley Additive exPlanations (SHAP), and features with consistently low importance or high collinearity were iteratively removed. When multiple features captured overlapping physiologic constructs, preference was given to features that improved interpretability or preserved performance with fewer predictors. Because gradient-boosted trees employ greedy split selection that assigns each split to the single most informative predictor, collinear features typically contribute little marginal information and are surfaced as low-SHAP features within the staged framework; SHAP-based pruning therefore served as the primary mechanism for identifying and removing redundancy, supplemented by clinical judgment for features capturing overlapping physiologic constructs (e.g., hours-above-threshold and area-above-threshold metrics derived from the same signal).
After evaluation of all domains, a final SHAP-guided reduction was performed to derive a parsimonious model minimizing dimensionality while preserving discrimination and calibration.
Model training and evaluation
Mortality prediction models were developed using Extreme Gradient Boosting (XGB) [13]. Class imbalance was addressed using the scale_pos_weight parameter, defined as the ratio of non-events to events in the training data. XGBoost’s native handling of missing values was used, allowing informative missingness patterns common in clinical workflows to be incorporated without explicit imputation.
Models were trained using five-fold stratified cross-validation, generating out-of-fold (OOF) predictions for all encounters. Discrimination was assessed using the area under the receiver operating characteristic curve (AUROC) and the area under the precision–recall curve (AUPRC). Calibration was evaluated using the Brier score. Sensitivity was fixed near 0.80 to reflect a clinically relevant rule-out operating point, with corresponding specificity and negative predictive value reported.
Bootstrap resampling of OOF predictions (500 iterations) was used to estimate 95% confidence intervals for AUROC and ΔAUROC. All models used consistent hyperparameters across experiments (learning rate 0.05, maximum tree depth 4, 300 boosting iterations).
Model interpretability
SHAP values were computed to quantify global feature importance and the directionality of feature effects, enabling transparent comparison between models using continuous physiologic trajectories, threshold-based representations, and expert-derived severity constructs.
Results
Cohort and feature domains
The study cohort included 3310 PICU encounters meeting criteria for suspected infection, of which 107 (3.2%) resulted in in-hospital mortality Table 1. All predictor variables were derived from data available within the first 24 hours of PICU admission.
Table 1.
Study cohort age distribution.
| Age group | Cases (Mortality = 0) | Cases (Mortality = 1) | Total Cases (Age Group) |
|---|---|---|---|
| 3–11 months | 511 (15.4%) | 20 (0.6%) | 531 (16.0%) |
| 12–23 months | 490 (14.8%) | 12 (0.4%) | 502 (15.1%) |
| 24–59 months | 686 (20.7%) | 20 (0.6%) | 706 (21.3%) |
| 60–143 months | 898 (27.1%) | 25 (0.8%) | 923 (27.9%) |
| 144–216 months | 618 (18.7%) | 30 (0.9%) | 648 (19.6%) |
| Total | 3203 96.8%) | 107 (3.2%) | 3310 (100%) |
Vitals-based domains: We examined the independent predictive performance contribution of 3 domains based on vital signs captured and recorded within the 24-hour period: vitals extrema (min/max), vitals dynamics (std/slope) and vitals severity/burden (hours/area of abnormality). Table 2 shows that virtually all patients had high density data recorded in the 24-hr. window and Table S1 in the supplement details the specific vital-signs based features derived and examined in each domain, together with the statistical distribution of medians and IQR against the mortality outcome that, together with SHAP analysis helped guide feature selection for inclusion in predictive model development.
Table 2.
Vitals data.
| Vital sign name | Count of patients with non-null values | Avg # of readings in 24 hr. window |
|---|---|---|
| Pulse | 3310 | 31.7 |
| SpO2 | 3310 | 31.2 |
| Resp | 3309 | 31.4 |
| SIa | 3309 | 19.7 |
| MAP | 3309 | 20.1 |
| SBP | 3309 | 20.2 |
| Temp | 3299 | 12.8 |
Abbreviations: SpO2 pulse oximetry, Resp respiratory rate, SI shock index derived from pulse and SBP, MAP mean arterial pressure, SBP systolic blood pressure, Temp body temperature.
SI derived.
Table S2 in the supplement details the treatment-based features available for predictive modeling, specifically vasopressors and supplemental O2. These features can be viewed as providing “context” to key vitals such as MAP and SpO2. From a statistical significance perspective all aspects of vasopressor treatment (e.g. timing of initiation, number of administrations, number of different vasopressors used, etc.) were strongly associated (P < 0.001) with mortality, whereas this was not the case for supplemental O2. A potential explanation may be the administering supplemental O2 (e.g., low flow cannula) in cases of mild hypoxia. The feature “resp_support_severe_24h” shows statistically significant association when either invasive mechanical ventilation, high flow nasal cannula, or positive pressure ventilation was inferred from the available data.
Table S3 in the appendix summarize the level off missingness of routine laboratory data available for predictive modeling. To minimize imputation-induced bias, laboratory features were restricted to those with <26% missingness in the first 24 hours. For included variables, missing values were handled using XGB native missing-value routing, in which decision trees learn default directions for absent observations at each split, allowing informative missingness patterns to be incorporated without explicit imputation. Table S3 provides the details of derived static and dynamic lab features evaluated for predictive performance. Finally, Tables S4 provide the definition and distribution of the composite SIRS features. Univariate comparisons across Supplementary Tables S1–S4 are descriptive and unadjusted for multiple testing, as described in the Methods.
Table 3 summarizes the stepwise evaluation of feature domains for early prediction of PICU mortality, highlighting changes in discrimination, calibration, and clinical operating characteristics as progressively richer clinical information was incorporated.
Table 3.
Sequential evaluation of feature domains for early PICU mortality risk prediction.
| Stage | Feature Domain included | RepresentativeHigh Value Added Features | AUROC (95% CI) | ΔAUROC vs prior (95% CI) | AUPRC (95% CI) | Brier | Spec @ 80% Sens | PPV | NPV | Retained for next stage (Y/N) |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Vitals Extrema | sbp_min, resp_max, map_min, si_max, spo2_min, pulse_max, temp_max | .66 (.61/.71) | .07 (.05/.1) | .06 | .40 | .04 | .98 | Yes | |
| 1 | + Vitals Dynamics | sbp_min, spo2_std, temp_std, pulse_std, map_std, resp_slope, resp_max, map_slope, temp_max, map_min | .72 (.67/.77) | .06 (.01/.10) | .09 (.07/.12) | .05 | .47 | .05 | .99 | Yes |
| 2 | + Severity/Burden | si_abnormal_hrs rr_hrs_above_uln, rr_area_above_uln, total_time_map_abnormal | .72 (.66/.77) | −.00 (−.03/.03) | .08 (.06/.10) | .04 | .50 | .05 | .99 | Yes |
| 3 | + Treatment | vaso_delay_ post PICU_hrs, total_administrations, min_qsfr, min sfr, q hourly slope, O2 delay hrs, sfr hourly slope, max flow rate | .82 (.77/.86) | .01 (−.01/.03) | .17 (.13/.23) | .03 | .66 | .07 | .99 | Yes |
| 4 | + Labs | hco3_std, ast_max, plt_min, bun_slope, wbc_max, bun_max, | .86 (.82/.90) | .05 (.01/.08) | .21 (.17/.29) | .03 | .76 | .10 | .99 | Yes |
| 5 | + Lab abn flags | hco3_abn, bun_abn, plt_abn, wbc_abn, ast_abn | .85 (.81/.89) | −.01 (−.02/.001) | .21 (.17/.29) | .03 | .74 | .09 | .99 | No |
| 6 | + SIRS | First_SIRS_positive_hrs, SIRS_risk_events_n_24 | .84 (.80/.88) | −.02 (−.04/−.001) | .21 (.16/.29) | .03 | .71 | .09 | .99 | No |
| 7 | + pSOFA | pSOFA approximate total score | .86 (.82/.90) | −.001 (−.01/.01) | .23 (.18/.31) | .03 | .77 | .10 | .99 | Yes |
| 8 | + Phoenix | Phoenix approximate total score | .86 (.82/.90) | −.001 (−.01/.01) | .21 (.17/0.29) | .03 | .79 | .11 | .99 | No |
| 9 | + Phoenix resp score element | min_phoenix_sfr_hourly | .86 (.82/.90) | .002 (−.01/.01) | .20 (.16/.27) | .03 | .79 | .11 | .99 | Yes |
| 10 | Parsimmonious model | Dropped clinically redundant features (29à10) | .83 (.82/.90) | −.03 (−.06/.00 | .17 (.13/.22) | .04 | .73 | .09 | .99 | Yes |
At Stage 0, a baseline model using only vital sign extrema demonstrated modest discrimination (AUROC 0.66, 95% CI: 0.61–0.71), consistent with prior observations that static physiologic snapshots provide limited prognostic resolution early in critical illness. Despite low positive predictive value (PPV 0.04), the model achieved high negative predictive value (NPV 0.98), indicating utility primarily for ruling out near-term mortality risk. In Stage 1, the addition of vital sign dynamics—including variability (standard deviation) and slopes—resulted in a statistically meaningful improvement in discrimination (AUROC 0.72, ΔAUROC +0.06, 95% CI: 0.01–0.10). This finding underscores the importance of short-term physiologic trajectories over isolated extrema, particularly for respiratory and hemodynamic variables. Improvements were also observed in specificity at fixed sensitivity (47% at 80% sensitivity), with preserved calibration (Brier score 0.05).
The inclusion of severity and burden metrics in Stage 2 did not materially improve model performance (ΔAUROC −0.00, 95% CI: −0.03 to 0.03). While these features captured cumulative physiologic derangement, their prognostic signal appeared largely redundant with already incorporated dynamic vital sign features. Nonetheless, these variables were retained for subsequent stages to preserve clinical completeness and assess downstream interactions.
A substantial performance gain emerged in Stage 3 with the addition of early treatment and respiratory support features. Discrimination increased markedly (AUROC 0.82, 95% CI: 0.77–0.86), accompanied by a pronounced rise in AUPRC (0.17 vs 0.08–0.09 in prior stages). This stage yielded the largest absolute improvement in clinical utility, with specificity increasing to 66% at 80% sensitivity and PPV rising to 0.07, while maintaining excellent NPV (0.99). These results suggest that early treatment intensity and response encode latent severity information not fully captured by physiologic data alone.
Further gains were observed in Stage 4 following the addition of selected laboratory features, with AUROC increasing to 0.86 (ΔAUROC +0.05, 95% CI: 0.01–0.08) and AUPRC reaching 0.21. Notably, laboratory trends and extrema (e.g., bicarbonate variability, platelet nadir, renal indices) improved discrimination without degrading calibration (Brier score 0.03), reinforcing their complementary role in early mortality risk stratification.
In contrast, the introduction of binary laboratory abnormality flags (Stage 5) did not confer additional benefits and slightly reduced discrimination (AUROC 0.85). Similarly, incorporation of SIRS-based features in Stage 6 resulted in a statistically significant decrement in AUROC (ΔAUROC −0.02, 95% CI: −0.04 to −0.001), indicating limited incremental prognostic value once granular physiologic, treatment, and laboratory data were available.
Composite organ dysfunction scores were evaluated in Stages 7-9. Addition of approximate pSOFA scores preserved discrimination (AUROC 0.86) and modestly improved AUPRC (0.23), whereas Phoenix composite scores did not yield further gains. However, a single Phoenix respiratory component (minimum S/F ratio) demonstrated comparable performance to the full score, supporting the notion that select subdomains may be more informative than aggregate indices in early prediction tasks.
Finally, Stage 10 presents the parsimonious model, in which clinically redundant features were removed, reducing dimensionality from 29 to 10 predictors. Although discrimination decreased modestly (AUROC 0.83), performance remained robust with preserved calibration and strong negative predictive value (0.99). This trade-off was considered acceptable given substantial gains in interpretability, stability, and feasibility for real-time clinical deployment.
Receiver operating characteristic (ROC) curves (Figure 1A) are shown for the best-performing all-domain model incorporating 29 features across vitals, dynamics, treatment, laboratory values, and expert-derived thresholds (solid blue line) and for the reduced parsimonious model limited to 10 data-driven features (dashed orange line). The all-domain model demonstrated higher overall discrimination (AUROC 0.864 [95% CI: 0.823–0.896]) compared with the parsimonious model (AUROC 0.833 [95% CI: 0.789–0.875]), with a modest but statistically significant reduction in AUROC (ΔAUROC −0.030 [95% CI: −0.064 to −0.002]). Despite feature reduction, the parsimonious model preserved sensitivity at the prespecified operating point (~80%) and maintained excellent negative predictive value (>0.99), supporting its suitability for early risk stratification with improved interpretability and deployment feasibility. SHAP analysis of the full feature model (Stage 9) demonstrated that early mortality risk was primarily driven by markers of hemodynamic instability, metabolic derangement, respiratory failure, and early treatment intensity rather than static thresholds alone. Delayed vasopressor initiation emerged as the strongest contributor to increased predicted risk, followed by markers of metabolic acidosis (bicarbonate variability), systolic blood pressure trajectories, pulse variability, and thrombocytopenia, reflecting both physiologic instability and evolving organ dysfunction. Laboratory features such as AST elevation, worsening renal indices, and markers of infection and inflammation contributed to additional prognostic signal, while respiratory impairment was captured through minimum S/F ratios and their temporal trends. Notably, dynamic features (slopes, variability, and delay metrics) consistently exerted greater influence than corresponding extrema, supporting the importance of short-term physiologic trajectories in early risk stratification. Composite scores (e.g., pSOFA) contributed comparatively modest incremental information once granular physiologic and treatment features were included.
Figure 1.
(A) Comparison of predicted performance between the full all-domain model and the final parsimonious model. (B) SHAP summary beeswarm plot of full feature model. (C) SHAP analysis of parsimonious model.
SHAP analysis of the parsimonious Stage 10 model (Figure 1B–C) demonstrated that a compact set of 10 features preserved the dominant risk signals identified in the full model, with delayed vasopressor initiation remaining the strongest contributor to early mortality risk. Measures of acid-base disturbance, thrombocytopenia, cardiovascular and autonomic variability (pulse variability, systolic blood pressure slope, MAP slope), and evolving renal and hepatic dysfunction (BUN slope, AST maximum) continued to exert substantial influence on model output. Notably, the retained features were predominantly dynamic or trend-based rather than static extrema, reinforcing the importance of short-term physiologic trajectories. Several of these slope-based features (sbp_slope, map_slope, bun_slope, q_hourly_slope) exhibited weak univariate associations with mortality in the supplementary tables but contributed substantially to multivariable model performance, consistent with slope features deriving their prognostic value through interaction with concurrent extrema, variability measures, and treatment context rather than as independent predictors. For example, a declining MAP slope carries different prognostic weight depending on whether vasopressor support has already been initiated—an interaction captured naturally by XGBoost’s tree-based architecture, which can condition feature effects on other variables. Despite the reduction in feature dimensionality, the relative ordering and directionality of SHAP effects were consistent with the full model, supporting the robustness, interpretability, and clinical plausibility of the parsimonious approach.
Discussion
We evaluated the incremental predictive value of clinically derived, knowledge-based features—burden metrics, laboratory abnormality flags, and composite organ dysfunction scores—within a modeling framework using dense, high-frequency physiologic and treatment data from children with suspected infection. While our cohort focused on pediatric sepsis, the findings have implications for other forms of life-threatening critical illness. Our principal finding is that when rich temporal data are available, data-driven feature dynamics dominate predictive performance, and the addition of expert-derived summaries provides little, if any, incremental benefit. The predominance of dynamic features in our models closely mirrors routine clinician reasoning at the bedside, where decisions are driven less by isolated values and more by whether physiologic and laboratory measures are improving, stable, or deteriorating over time.
Knowledge-based features, including burden metrics and composite scores, were developed to address data sparsity and cognitive load in clinical settings where continuous monitoring or repeated laboratory assessments are unavailable. In contrast, when granular temporal data are available, machine learning models can directly exploit raw physiologic trajectories and treatment responses without intermediate expert summarization. Under these conditions, the model learns its own representations of severity and deterioration, rendering handcrafted abstractions largely redundant.
Clinically, the proposed approach is viewed as a continuously updating risk stratification tool, rather than a static severity score. Rather than assigning patients to fixed risk categories based on worst observed values, the model evaluates response to therapy over time, identifying patients whose physiology is failing to improve despite intervention—even when conventional thresholds for hypotension or organ dysfunction have not yet been crossed. This framing aligns naturally with real-time clinical decision-making and supports early recognition of treatment non-response.
An additional finding of clinical relevance was the dominant contribution of vasopressor initiation delay to early mortality risk. Current pediatric resuscitation practices commonly emphasize repeated fluid boluses prior to vasopressor initiation, particularly in the absence of overt hypotension. While intended to optimize preload, this approach may delay vasoactive support in children with evolving, fluid-refractory shock. The strong predictive signal associated with vasopressor delay suggests that treatment timing itself captures latent severity not fully reflected in contemporaneous vital signs or laboratory values. This signal may also reflect delayed recognition of shock physiology rather than treatment intensity alone. Vasopressor initiation timing reflects not only patient physiology and individual clinician recognition but also institutional escalation thresholds and resuscitation practices. The predictive contribution of treatment-timing features should therefore be interpreted as associational rather than causal: SHAP values quantify how strongly these features predict mortality given the observed data, not the effect of altering treatment timing on outcomes.
Traditional scoring systems such as SIRS, pSOFA, and Phoenix remain valuable for standardized assessment and retrospective benchmarking but may be less effective for early decision support in data-dense environments. These scores rely on discretization of continuous variables and fixed weighting schemes that do not explicitly incorporate rich, high-frequency physiologic trajectories or treatment timing. Moreover, composite scores often commingle direct indicators of clinical intervention (e.g., vasopressor use) with expert-defined thresholding of physiologic measurements, making it difficult to disentangle evolving severity from clinician response. In contrast, the present approach models physiologic dynamics and treatment timing directly, allowing risk estimates to evolve continuously as new data become available.
Several respiratory rate–based features contributed less to model performance than might be anticipated clinically. Likely contributors include the inherent measurement variability of pediatric respiratory rate, redundancy with stronger oxygenation features such as minimum S/F ratio and SpO2 burden, and a cohort in which mortality was dominated by hemodynamic rather than primary respiratory pathophysiology—with minimum respiratory rate (often reflecting peri–arrest bradypnea) remaining a notable exception (P < 0.001).
Sepsis provides a clear, illustrative use case. Current clinical workflows emphasize repeated reassessment during early resuscitation, yet existing tools offer limited support for determining when a patient is failing to respond to fluids and may require earlier escalation to vasoactive therapy. The strong predictive signal associated with vasopressor initiation delay in this study suggests that trajectory-based models could serve as a decision support adjunct. Such models could flag patients with persistent metabolic acidosis, worsening hemodynamics, or lack of physiologic improvement despite treatment—signals that may precede overt shock. Importantly, this approach is not limited to sepsis and may generalize to other forms of pediatric critical illness characterized by evolving cardiopulmonary or metabolic instability.
The model relies exclusively on routinely collected PICU data, including high-frequency vital signs, commonly ordered laboratory tests, and time-stamped treatment events. No specialized devices, novel biomarkers, or institution-specific data elements are required, and the approach is agnostic to EHR vendor. Required resources are limited to real-time data extraction, periodic model inference, and integration into existing clinical dashboards or alerting workflows.
These findings support a shift from static, threshold-based scoring toward trajectory-aware decision support that augments clinician judgment rather than replacing it. By formalizing how clinicians already assess response to therapy, data-driven models such as the one presented here may enable earlier recognition of treatment failure and more timely escalation of care, with potential implications for improving outcomes across a broad range of critical illnesses.
Limitations
Several limitations warrant emphasis. First, this is a single-center study conducted at a tertiary, freestanding pediatric hospital. Several of the most influential features in our models—vasopressor initiation timing, respiratory support type and timing, and laboratory ordering frequency—are closely tied to local clinical practice patterns and escalation thresholds. Practice variability across institutions may therefore influence both the absolute importance of individual features and the calibration of model outputs. Although the broader finding that dynamic, treatment-aware representations outperform threshold-based summaries in dense-data settings is likely to be robust, external validation in geographically and operationally diverse multi-center cohorts will be required before clinical deployment, and recalibration may be needed at sites with materially different patient populations or workflows.
Second, our primary outcome mortality was selected because it is objective, clinically meaningful, and comparable across institutions and prior literature. However, mortality is a coarse endpoint that may fail to capture clinically relevant deterioration short of death—such as new-onset organ dysfunction, escalation to extra-corporeal support, or prolonged mechanical ventilation—that may be equally relevant for early triage and decision support. Future work should evaluate composite deterioration outcomes that better reflect the full spectrum of severe illness.
Finally, the 3.2% mortality rate poses challenges of training-data imbalance. We mitigated this through stratified five-fold cross-validation and class-weight adjustment in XGB training, and characterized the residual uncertainty through AUPRC reporting alongside AUROC and bootstrap 95% confidence intervals on all metrics. Even with these measures, performance estimates carry wider uncertainty than at higher event rates, and PPV is bounded by prevalence (0.115 and 0.090 for the full and parsimonious models at sensitivity 0.804)—appropriate for decision support rather than diagnosis. Small between-model differences should be interpreted in light of overlapping intervals, and external validation in higher-event-rate cohorts will further test stability.
Supplementary Material
Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.ijid.2026.108899.
Funding
This study was supported by NIAID (National Institute of Allergy and Infectious Diseases) R41AI167224 and NIGMS (National Institute of General Medical Sciences R35GM157142 (I.K.).
Footnotes
Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Ethical approval
The study was approved by the Pearl institutional review board (IRB Number:22-CTAI-101) with waiver of informed consent.
CRediT authorship contribution statement
Tom Velez: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Software, Visualization, Writing – original draft. Oluwakemi Badaki-Makun: Data curation, Investigation, Writing – review & editing. Danielle Hirsch: Data curation, Investigation, Writing – review & editing. Danielle Claire Mercurio: Data curation, Investigation, Writing – review & editing. Holly Depinet: Data curation, Investigation, Writing – review & editing. Maya Dewan: Data curation, Investigation, Writing – review & editing. Rishikesan Kamaleswaran: Data curation, Investigation, Writing – review & editing. Jocelyn Grunwell: Data curation, Investigation, Writing – review & editing. Maria Triantafyllou: Data curation, Investigation, Writing – review & editing. Nora Wolff: Data curation, Investigation, Writing – review & editing. Fehima Abdelrahman: Data curation, Investigation, Writing – review & editing. Charles Macias: Data curation, Investigation, Writing – review & editing. Ioannis Koutroulis: Conceptualization, Data curation, Funding acquisition, Investigation, Supervision, Validation, Writing – review & editing.
References
- [1].Frankel LR, Hsu BS, Yeh TS, Simone S, Agus MSD, Arca MJ, et al. Criteria for critical care infants and children: PICU admission, discharge, and triage practice statement and levels of care guidance. Pediatr Crit Care Med 2019;20(9):847. doi: 10.1097/PCC.0000000000001963. [DOI] [PubMed] [Google Scholar]
- [2].Sinha A, Stevens LA, Su F, Pageler NM, Tawfik DS. Measuring electronic health record use in the pediatric ICU using audit-logs and screen recordings. Appl Clin Inform 2021;12(4):737–44. doi: 10.1055/s-0041-1733851. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [3].Eisenberg MA, Balamuth F. Pediatric sepsis screening in US hospitals. Pediatr Res 2022;91(2):351–8. doi: 10.1038/s41390-021-01708-y. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [4].Malik A, Taksande A, Meshram R, Malik A, Taksande A, Meshram R. Pediatric sequential organ assessment score: a comprehensive review of the prognostic marker in the pediatric intensive care unit. Cureus 2024;16. doi: 10.7759/cureus.60034. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [5].Matics TJ, Sanchez-Pinto LN. Adaptation and validation of a pediatric sequential organ failure assessment score and evaluation of the sepsis-3 definitions in critically ill children. JAMA Pediatr 2017;171(10):e172352. doi: 10.1001/jamapediatrics.2017.2352. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [6].Balamuth F, Scott HF, Weiss SL, Webb M, Chamberlain JM, Bajaj L, et al. Validation of the pediatric sequential organ failure assessment score and evaluation of third international consensus definitions for sepsis and septic shock definitions in the pediatric emergency department. JAMA Pediatr 2022;176(7):672–8. doi: 10.1001/jamapediatrics.2022.1301. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [7].Sanchez-Pinto LN, Bennett TD, DeWitt PE, Russell S, Rebull MN, Martin B, et al. Development and validation of the Phoenix criteria for pediatric sepsis and septic shock. JAMA 2024;331(8):675–86. doi: 10.1001/jama.2024.0196. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [8].Chanci D, Grunwell JR, Rafiei A, Brown SR, Ripple MJ, Bishop NR, et al. Machine learning model for daily prediction of pediatric sepsis using Phoenix criteria. Pediatr Res 2025:1–7. doi: 10.1038/s41390-025-04221-8. [DOI] [PubMed] [Google Scholar]
- [9].Soghier LM, Fratantoni K, Reyes C. Reference range values for pediatric care. Mullins K, editor. 2nd ed.. American Academy of Pediatrics; 2019. 10.1542/9781610022811. [DOI] [Google Scholar]
- [10].Liang H, Carey KA, Jani P, Gilbert ER, Afshar M, Sanchez-Pinto LN, et al. Association between mortality and critical events within 48 hours of transfer to the pediatric intensive care unit. Front Pediatr 2023;11. doi: 10.3389/fped.2023.1284672. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [11].Gennatas ED, Friedman JH, Ungar LH, Pirracchio R, Eaton E, Reichmann LG, et al. Expert-augmented machine learning. Proc Natl Acad Sci U S A 2020;117(9):4571–7. doi: 10.1073/pnas.1906831117. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [12].Digitale J, Franzon D, Ge J, McCulloch C, Pletcher MJ, Gennatas ED. Expert-augmented machine learning for predicting extubation readiness in the pediatric intensive care unit. BMC Med Inform Decis Mak 2025;25(1):232. doi: 10.1186/s12911-025-03070-z. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [13].Chen T, Guestrin C. XGBoost: a scalable tree boosting system. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, in KDD ’16. New York, NY, USA: Association for Computing Machinery; 2016. p. 785–94. doi: 10.1145/2939672.2939785. [DOI] [Google Scholar]
- [14].Agulnik A, Gossett J, Carrillo AK, Kang G, Morrison RR. Abnormal vital signs predict critical deterioration in hospitalized pediatric hematology-oncology and post-hematopoietic cell transplant patients. Front Oncol 2020;10. doi: 10.3389/fonc.2020.00354. [DOI] [PMC free article] [PubMed] [Google Scholar]
Associated Data
This section collects any data citations, data availability statements, or supplementary materials included in this article.