Blood Urine Nitrogen Trajectories of Acute Pancreatitis Patients in Intensive Care Units

Introduction

Acute pancreatitis (AP) is an inflammatory condition commonly caused by bile stones or excessive dose of alcohol and is one of the most prevalent gastrointestinal diseases that needs hospitalization.^1–3 About 20–40% of AP patients will develop into severe AP (SAP) with persistent organ failure that requires intensive care unit (ICU) admission and had about 30% risk of mortality.^4–7

The blood urea nitrogen (BUN) level can provide valuable insights into the physiological condition in AP, including intravascular volume depletion, prerenal azotemia, kidney injury, potential gastrointestinal bleeding, and the response to fluid resuscitation.^8–10 Therefore, monitoring BUN levels can guide clinicians in taking appropriate interventions to manage AP patients effectively, and BUN has been identified as an important prognostic factor for AP. Previous studies suggested that as a dynamic biomarker, both baseline and early trends of BUN change were significantly associated with mortality in AP patients.^9–11 Wu et al found that the odds ratio (OR)s of in-hospital mortality was 2.9 and 2.2 for every 5 mg/dL increment of admission/24-hour change in BUN in patients with AP.¹⁰ In another study, the pooled ORs of mortality in patients with AP were 4.60 and 4.28 for elevated admission BUN (≥20 mg/dL) and any increment of BUN in the first 24 hours, respectively.⁹ Good prognostic performance has also been demonstrated by both direct BUN levels and BUN-based prediction models for AP patients.^9,12 In the ICU, BUN is a repeatedly measured biomarker with dynamic developmental patterns, while studies to date have only analyzed the association between baseline/first 24-hour difference in BUN levels and AP mortality. To our knowledge, no study has analyzed BUN longitudinal data, although current evidence indicates that the dynamic change of BUN is essential for AP prognosis.¹⁰ Therefore, by fitting with novel trajectory models, it is possible to explore meaningful heterogeneous BUN trajectories in AP patients.

To identify clinically meaningful subclasses of AP patients using BUN trajectories in the ICU, we conducted a retrospective observational study based on a large multi-source intensive care database from the West China Hospital (WCH) system and three multicenter public intensive care databases in the United States.

Materials and Methods

Study Design and Setting

The present study was a multicenter, retrospective observational study using data from China and the United States. ICU inpatients from the Healthcare-Associated Infections of ICUs (ICU-HAIs) Registry, a large multi-source intensive care database of the WCH system, were included as the development cohort. As a major tertiary healthcare system in western China, the WCH system consists of three independent healthcare institutions (University Campus, Wenjiang Hospital, and Jinjiang Hospital) with 4300 beds and more than 5.3 million annual visits. The ICU-HAIs registry included patients from six ICUs (general, surgical, neurological, respiratory, thoracic, and pediatric) with over 8000 annual admissions in the WCH system. Through integrating three databases (ICU-HAI system, ICU system, and electronic medical record system), the ICU-HAIs registry contains detailed information regarding clinical care of ICU patients. The registry included 492 variables of those, many were stored as time-dependent variables. For instance, information regarding vital signs and life supports were recorded on an hourly basis. Quality assessment showed that the registry is a valuable database with high level of quality and comprehensiveness. A detailed description of the ICU-HAIs registry was presented in our previous work.^13–15

For external validation, three databases without patients overlapping, including Medical Information Mart for Intensive Care (MIMIC)-III CareVue subset, MIMIC-IV database that contained ICU patients in Beth Israel Deaconess Medical Center, and eICU Collaborative Research Database (eICU-CRD), a multi-center database that included patients admitted to one of 208 hospitals in the United States were integrated as an external validation cohort.¹⁶ This study was approved by the Ethical Committee of West China Hospital (WCH2018-409). The author completed the required training (Collaborative Institutional Training Initiative ID: 37,078,282) and signed related agreements to access restricted public data.

Study Population

In the ICU-HAIs registry (development) cohort, all patients diagnosed with AP (admitted between 2010 and 2021) met the inclusion criteria. Patients with AP who met any of the following criteria were excluded: 1) without multiple BUN measurements in the 21 days after ICU admission, 2) age <18 years, 3) ICU length of stay <24 hours. To ensure accuracy in the analysis, we excluded AP patients with ICU length of stay <24 hours since very short ICU stays may not reflect the true trend of BUN at 21 days. The external validation cohort (admitted between 2001 and 2019) had identical inclusion and exclusion criteria, except that an additional 727 AP patients without ICU admission in the MIMIC-IV database were excluded.

Variables Measurement and Outcomes

The longitudinal BUN trajectory in the 21 days after ICU admission was included and summarized each day. If there were multiple measurements in one day, the mean value was calculated for analysis.

Patient’s baseline characteristics, including demographic data (age, sex, ethnicity), first measurements of vital sign or laboratory test after ICU admission [(heart rate, repatriation rate (RR), temperature, mean arterial pressure (MAP), oxygen saturation (SpO2), creatinine, hemoglobin, sodium, potassium, calcium, phosphate, white blood cell (WBC)], and comorbidities (diabetes, renal disease, and liver disease) were extracted. The outcome of this study was 30-day mortality and was identified by the difference between mortality time and ICU admission.

Statistical Analysis

Latent Class Trajectory Model Development and External Validation

In the developed cohort, latent class trajectory modelling (LCTM) was conducted to identify subclasses of AP patients based on the 21 days trajectory of BUN measurements after ICU admission. LCTMs assume a homogeneous population and were operated to determine if subclasses exist within the population that follows qualitatively distinct developmental trajectories.¹⁷ The optimal number of Class(K) was detrended by the Akaike information criterion (AIC), Bayesian information criteria (BIC), sample-adjusted BIC (SBIC), and entropy. Lower AIC, BIC, SABIC, and higher entropy indicated better model fitting and discriminatory power. In addition, to substantiate the clinical meaning and robust result for statistical analysis, the minimum size for each subgroup was defined as 2%.^18,19 Finally, the BUN trajectory of AP patients within 21 days of ICU admission in the external validation cohort was included in the final model analysis. The final model will provide posterior classification individual class-membership probabilities for each patient trajectory, thereby dividing into different subclasses and proceeding to the next step of the analysis.

Comparing Clinical Characteristics and Outcome Between Subclasses

Baseline characteristics of subclasses in both the development and validation cohorts were compared. Continuous variables were described using the median and interquartile range (IQR), and categorical variables were presented as counts and percentages. Kruskal–Wallis tests were used to analyze differences between subclasses of AP patients for continuous variables, and chi-square tests were conducted to compare categorical variables. Additionally, as a commonly concurrently tested indicator, we extracted and analyzed creatinine 21-day trajectories across BUN subclasses.

Multivariate Cox regression models were fitted to estimate the hazard ratios (HR)s of 30-day mortality across subclasses in both cohorts.^20,21 A two-sided P-value <0.05 was considered statistically significant.

Construction of Classifier for Identification of Special AP Subclasses

In this study, a subgroup of patients with AP exhibiting a rapid increase in BUN levels, associated with the highest mortality risk, was identified. Since the 21-day BUN levels in AP patients in the ICU have very complex fluctuation patterns, distinguishing these patients based solely on baseline BUN levels is challenging. To address this, a binary classification eXtreme Gradient Boosting (XGBoost) model was developed using patients’ baseline clinical data. Hyperparameter optimization was performed via 5-fold cross-validation. The area under the receiver operating characteristic curve (AUROC) was calculated for model performance evaluation. Additionally, the SHapley Additive exPlanations (SHAP) algorithm was used for model interpretation.

To enhance the clinical utility and significance of the model, an online interactive server was created (https://wzcjerry.shinyapps.io/AP_Prediction/).

Missing Data Management

The proportion of missing data in the development and validation cohorts was shown in Supplementary Table1. Multiple linear imputation based on decision trees algorithm was used to impute missing values. For longitudinal BUN trajectories data, no imputation process was performed for BUN data because the trajectory model does not constrain individuals to have the same length of observation and patients without any BUN observation were excluded.

All data were analyzed using R statistical software (version 4.1.2). Multivariate imputation was performed with the “mice” package, while LCTMs were developed using the “lcmm” package. Cox regression models were conducted using the “survival” package. The XGBoost model was trained with the “xgboost” package, and the SHAP algorithm was applied using “SHAPforxgboost”. The online server was created with the “Shiny” package.

Results

Subjects and Baseline Characteristics

In the ICU-HAIs registry from WCH systems, 3611 patients were diagnosed with AP. After excluding related patients according to the criteria, 2971 patients were finally included as development cohort. A total of 343, 931, and 611 patients were diagnosed with AP in MIMIC-III CareVue, MIMIC-IV, and eICU-CRD database and finally 300, 157 and 473 patients were included after excluding relevant patients according to the criteria, respectively (Figure 1). Eventually, the included AP patients from MIMIC-III/IV and eICU-CRD were integrated into the validation cohort (N = 930). The baseline characteristics were presented in Supplementary Table2.

Figure 1 The flow chart of patients recruitment and acute pancreatitis subclasses development/validation.

In the development cohort, the median age was 48 years, and 65.9% of patients were male. All patients in this cohort were Asian, 185(6.2%) occurred ICU mortality and 197(6.6%) occurred hospital mortality. In the validation cohort, 58.5% of patients were male, with 68.9% of patients being white, 7.5% being black, 1.4% being Asian and 22.2% of patients in other ethnicities. Within the validation cohort, 9.8% of patients occurred ICU mortality, while 13.1% occurred hospital mortality.

Subclasses LCTMs Modeling and Trajectories Analysis

The fit statistics of LCTMs with the number of Class(es) from 1 to 6 were demonstrated (Supplementary Table3). As the number of Class(es) increased, the log-likelihood also increased, while AIC, BIC, and SBIC declined continuously. However, as the number of Class(es) was greater than 5, the size of minimize subgroup was smaller than 2%, thus the 5-subclasses model (LCTM-5) was considered as the final model. The related fixed effect coefficient based on maximum likelihood estimation was represented in Supplementary Table4.

LCTM-5 split the development cohort into Class 1 [361 (12.2%)], Class 2 [2244 (75.4%)], Class 3 [118 (4.0%)], Class 4 [70 (2.4%)] and Class 5 [178 (6.0%)], respectively. Analysis of LCTM-5 trajectories in the development cohort revealed distinct patterns. Class 2 exhibited consistently low BUN levels, with slight fluctuations between 10 and 20 mg/dL and was designated as the “Non-azotemia” Class (Figure 2A). Class 3 started with the highest initial BUN level, which gradually decreased over 21 days, characterizing it as the “Severe-azotemia, slow decreasing” group. Class 1 had a moderate initial BUN level (50–60 mg/dL), with a gradual decrease over the 21 days, labeling it as the “Moderate-azotemia, slow decreasing” group. Class 4 and 5 both began with moderate initial BUN levels (30–40 mg/dL). However, the BUN levels in these groups exhibited different patterns. Class 5 showed a slow increase, while Class 4 displayed a rapid increase, leading to their categorization as Moderate-azotemia, slow increasing and Moderate-azotemia, rapid increasing, respectively. The trajectories of BUN for individuals were depicted in Figure 2B, revealing intricate dynamic patterns that defy simple description based on baseline levels. Figure 2C illustrates the daily distribution of BUN levels across different AP groups. Notably, as the days following admission to the ICU increase, it becomes evident that the distribution of BUN in Group 4 tends to concentrate around higher values compared to Group 3. Creatinine trajectories showed similar trend to BUN trajectories (Supplementary Figure1).

Figure 2 Trajectories and time distribution of blood urine nitrogen among acute pancreatitis subclasses (Figure 2 (A)) represented the group-based BUN trajectories in the development cohort; (Figure 2 (B)) represented individual BUN trajectories in the development cohort. AP patients in different subclasses were distracted by color. The first five plots showed trajectories for specific subclasses from 1–5 and the last plot represented all trajectories together. Each black point represented a single BUN observation and BUN observations were linked by line individually; (Figure 2 (C)) represented the BUN distribution over 21 days across AP subclasses; In all plot, Class 1 were colored dark blue, Class 2 were colored red, Class 3 were colored green, Class 4 were colored sky blue, and Class 5 were colored purple; (Figure 2 (D)) represented the group-based BUN trajectories the validation cohort.

In the validation cohort, LCTM-5 identified distinct subclasses among acute pancreatitis (AP) patients, comprising 125 (13.4%) in group 1, 669 (1.9%) in group 2, 54 (5.8%) in group 3, 21 (2.3%) in group 4, and 61 (6.6%) in group 5. Remarkably, the trajectories observed in each subgroup within the validation cohort exhibited an identical pattern to those observed in the development cohort (Figure 2D).

Patient Characteristics Between Subclasses in Development and Validation Cohorts

The clinical characteristics across AP subclasses in the development cohort were represented (Figure 3) (Supplementary Table5). Significant differences in multiple parameters were observed across AP subgroups, including age, gender, heart rate, RR, temperature, MAP, SpO2, creatinine, hemoglobin, sodium, potassium, calcium, phosphate, pH, WBC, hypertriglyceridemia AP, and renal disease. Compared to the total population (mean value) in the development cohort, Class 2 showed higher levels of MAP, SpO2, hemoglobin, calcium and pH. Class 1 had lower MAP, SpO2, calcium and pH. Class 3 had lower SpO2, hemoglobin and calcium. Class 4 had lower MAP, SpO2, and pH. Class 5 had lower MAP, SpO2, calcium, and pH (Figure 3).

Figure 3 The clinical characteristics of acute pancreatitis subclasses (Figure 3 (A-E)) showed the clinical characteristic of subclasses 1–5 in the development cohort, respectively. Red links represented a higher value/proportion of variable compared with total AP patients while blue links represented a lower value /proportion of variable compared with total AP patients;The width of links was calculated by the mean value/proportion of variables of each subgroup divided by the mean value/ proportion of total patients if the level is above the total patient. Otherwise, the width of links was calculated by the mean value/ proportion of total patients divided by the mean value/proportion of variables.

The clinical characteristics across AP subclasses in the validation cohort were represented in Supplementary Table6. The AP subclasses in the validation cohort showed clinical characteristics essentially consistent with those in the development cohort. There were statistical differences between the groups in age, sex, RR, temperature, MAP, creatinine, hemoglobin, sodium, potassium, calcium, phosphate, pH, WBC, and renal disease. There was a significant difference in the proportion of liver disease in the validation cohort, with the highest proportion in class 3 (35.2%) and the lowest in class 1 (16.8%).

Difference in Prognosis Between Subclasses in Development and Validation Cohorts

In Cox regression analysis, adjusting for demographic variables, creatinine level and renal disease status, with Class 2 as the reference, a significantly elevated 30-day mortality risk was observed in Class 1 (hazard ratio (HR) 4.08, 95% confidence interval (CI) 2.61–6.37, p < 0.001), Class 3 (HR 3.62, 95% CI 1.80–7.28, p < 0.001), Class 4 (HR 5.32, 95% CI 2.62–10.82, p < 0.001), and Group 5 (HR 3.50, 95% CI 2.07–5.91, p < 0.001) in the development cohort (Figure 4A). The validation cohort yielded consistent results with the development cohort, showing a significantly increased risk of mortality in other subclasses when compared to Class 2. Notably, mirroring the results of the development cohort, Class 4 exhibited the highest mortality risk among all subclasses (HR 6.23, 95% CI 2.93–13.22, p < 0.001) (Figure 4B). The complete results of Cox regression were presented in (Supplementary Tables 7,8).

Figure 4 The hazard ratios of 30-day mortality across acute pancreatitis subclasses in development and validation cohorts (Figure 4 (A)) showed the HRs and 95% CI across AP subclasses in the development cohort; (Figure 4 (B)) showed the HRs and 95% CI across AP subclasses in the validation cohort.

Model Performance and Explanation for Classifier to Identify Class-4 Patients

Classifier to identify Class-4 patients among “Moderate-azotemia” (Class 1, Class 4, Class 5) were trained by the combined dataset of the development and validation cohorts based on the XGBoost algorithm. The process included hyperparameter optimization through 5-fold cross-validation, as depicted in Supplementary Figure2. The final model demonstrated strong discriminative performance, with an AUROC of 0.894 when applied to the original dataset (Figure 5A). The significance of variables was assessed using the SHAP algorithm, revealing that age, phosphate, and potassium were the most influential features in identifying Class 4 patients in the XGBoost classifier (Figure 5B). To enhance clinical accessibility, the classifier was integrated into an interactive online server for the convenience of healthcare professionals (Supplementary Figure3).

Figure 5 Model performance and explanation for class-4 classifier among moderate-azotemia acute pancreatitis patients (Figure 5 (A)) represented the ROC curve for the classifier model; (Figure 5 (B)) represented the SHAP model explanation result for final classifier model. In this figure, each point corresponds to a prediction. The higher the SHAP value, the more the prediction result of the model is biased toward a positive outcome (ICU death). The lower the SHAP value, the more negative the prediction result of the model is. The more purple the color of each point is, the higher the actual value of the feature is, and vice versa. The order of the features used from top to bottom represents the ranking of the feature’s contribution to the model from high to low.

Discussion

We applied LCTMs to analyze 21-day BUN trajectories of AP patients in ICU and revealed five distinct subclasses, each characterized by unique azotemia conditions and different associated prognoses. These subclasses include Class 1: “Moderate-azotemia, slow decreasing”, Class 2: “Non-azotemia”, Class 3: “Severe-azotemia, slow decreasing”, Class 4: “Moderate-azotemia, rapid increasing” and Class 5: “Moderate-azotemia, slow increasing”. Significant disparities in clinical characteristics and prognoses were evident across these subclasses in both the development and validation cohorts. In comparison to Class 2, AP patients in the ICU with baseline azotemia condition exhibited a significantly higher risk of 30-day mortality, with HRs ranging from 4.01 to 6.45 in the development cohort and from 2.24 to 5.88 in the validation cohort.

As a routine collected laboratory test, BUN is more accessible in ICU, and existing evidence shows that BUN is the most accurate single predictor of in-hospital death compared with laboratory tests. Both elevated BUN (20–22 mg/dL) at admission and any rise in BUN within the first 24 hours of hospitalization (despite initial BUN level) significantly increased the risk of mortality of AP patients.^9,10 Time-specific or model-depended ROCs also demonstrated that BUN had a good predictive performance for mortality (AUROC: 0.82–0.91).^9,10,12 However, although existing studies claimed to use continuously measured BUN data, the longitudinal trajectory of BUN has not been analyzed in existing studies. Moreover, approximately half of AP mortality occurs in the first two weeks, while previous studies only analyzed BUN levels during the first 48 hours of hospitalization.^20,21 As a time-sensitive disease, AP had dynamically fluctuating BUN levels; therefore, it is necessary to study the trajectories of BUN in longer disease process based on longitudinal data.

Our study aligns with existing evidence, indicating that AP patients with higher baseline BUN levels faced a greater risk of mortality. However, the HR estimates in the development and validation groups highlight important findings that have not been previously reported. Among AP subclasses, Class 4, representing patients with a baseline of moderate-azotemia and a rapid BUN increasing condition, exhibited the highest risk of 30-day mortality, which notably surpassed Class 3 (highest baseline BUN). This underscores our prior hypothesis that relying solely on baseline BUN levels fails to capture the dynamic pattern of disease processing for AP in the ICU, potentially overlooking unique patient profiles deserving of special attention.

Compared to other AP subclasses, Class 4 exhibited a lower MAP level and a higher proportion of renal disease. Low MAP indicates underlying hypoperfusion, which in turn leads to pancreatic microcirculatory disturbances. This factor significantly impacts the early stages of the disease and is considered a major contributor to the development of pancreatic necrosis and organ failures^22–25 with increased BUN.^9–12 The development pattern of Class 4 AP patients, to some extent, reflects the patient’s possible initial state of low perfusion and the subsequent development of organ damage. This observation aligns with our findings of a higher prevalence of renal in Class 4 patients.

This study represents the first investigation of the longitudinal BUN trajectories of ICU-admitted AP patients, identifying five stable AP subclasses based on BUN trajectories and validating these results with external datasets. However, our findings have several limitations. First, the study’s focus on ICU patients inherently limits the range of its clinical applicability. However, SAP patients often require ICU admission, especially when signs of multiorgan failure are present, a condition that defines AP as severe with more careful management needs. Hence, we meticulously examined the dynamic progression of BUN in ICU-admitted AP patients within the largest international cohort known to date, thereby uncovering a subgroup of patients at the highest risk, a revelation hitherto unreported. Second, although our LCTM-5 model can classify BUN trajectories of varying length and frequency, it was developed and validated using 21-day trajectories, limiting its application for early prediction. However, it is important to emphasize that our primary goal was to analyze 21-day BUN trajectories in ICU-admitted patients with AP and to identify distinct subclasses of AP. We also created an online model to identify the special subclass, which is challenging when relying solely on baseline characteristics, enhanced the clinical utility and significance of our finding. Finally, our study did not examine responses to interventions across subclasses. Nevertheless, we aimed to deepen our understanding of the heterogeneity of AP disease and successfully identified robust AP subclasses with unique clinical characteristics and significantly different prognoses. This provides valuable insights into the development of AP and paves the way for future studies focused on early subclass prediction in the ICU and the discovery of differential responses to treatment.

Conclusion

In conclusion, we demonstrated the important heterogeneity of BUN trajectories among AP patients in ICU. Five robust and clinically meaningful AP subgroups were identified with different statuses of azotemia. We have identified an AP subclass with moderate azotemia statue while had highest risk of mortality, which has been overlooked in previous studies that only considered baseline BUN. The present study promoted our understanding of AP longitudinal process and further work is needed to connect AP subclasses and specific treatment.