The target cohort in this study was obtained from the Medical Information Mart for Intensive Care IV (MIMIC-IV) version 2.0 (24), which is a freely accessible critical care database that contains detailed demographics, laboratory tests, various notes and reports, and real-time high-resolution hourly data (vital signs, intravenous medications infusion rates, etc.) acquired from the hardware devices (monitor, infusion pumps, etc.) of over 76,000 patients during ICU admission at Beth Israel Deaconess Medical Center. This study aimed to determine whether VRR is significantly associated with ICU mortality and in-hospital mortality for patients with septic shock. Briefly, we utilized the following inclusion criteria in the identification of the study cohort: (1) patients with sepsis defined by Sepsis-3 criteria; (2) length of stay in ICU over 24 h; (3) age 18 years or older; (4) weight data (if needed) to calculate particular subsections of the VIS; (5) treatment with at least one vasoactive–inotropic agent; (6) therapy with vasoactive–inotropic agents within 24 h after ICU admission; and (7) subsequent ICU stay greater than 24 h after initiation of vasoactive–inotropic therapy. Exclusion criteria: discharge from ICU against medical advice.

The real-time VIS was queried and calculated from the immediate infusion rates of dopamine, dobutamine, epinephrine, milrinone, vasopressin, and norepinephrine. For multiple infusion rates recorded in a given hour for the same vasoactive–inotropic agent, we took the maximum rate as the input value in the specific hour for the dynamic VIS data. The dynamic VIS data with 1 h temporal resolution of patients contained 7 features: 6 subsections of VIS (VIS of dopamine, dobutamine, epinephrine, milrinone, vasopressin, norepinephrine) and the total VIS.

Then, we took the maximum value of each VIS item within 1–24 and 25–48 h and calculated VRR, thus forming a data matrix with 15 features (maximum VIS of dopamine, dobutamine, epinephrine, milrinone, vasopressin, norepinephrine, total VIS within 1–24 and 25–48 h, respectively, and VRR). We called the data matrix with 15 features “static VIS data” by comparison with the high-resolution dynamic VIS data.

Only the data during the first 24 h for each patient after ICU admission were used as covariates in this study for further analysis. We merged the fields of demographic and admission information, severity scores, comorbidities, vital signs, interventions, laboratory tests, VRR, and outcome to form the final cohort data matrix.

The formula of VIS is defined as follows (13): VIS = dopamine dose (μg·kg⁻¹·min⁻¹) + dobutamine dose (μg·kg⁻¹·min⁻¹) + 100 × epinephrine dose (μg·kg⁻¹·min⁻¹) + 10 × milrinone dose (μg·kg⁻¹·min⁻¹) + 100 × norepinephrine dose (μg·kg⁻¹·min⁻¹) + 10,000 × vasopressin dose (U·kg⁻¹·min⁻¹).

VRR is a totally new concept that we proposed to quantify dosage reduction for vasoactive–inotropic agents. The formula of VRR is defined as follows: (VIS_{1-24h max} – VIS_{25-48h max})/VIS_{1-24h max}.

The primary outcome in this study was ICU mortality, and the secondary outcome was in-hospital mortality.

The data were randomly split into a training dataset (80%) and a test dataset (20%) according to the labels of outcome events by the caTools package in R version 4.0.2. Constant “stay_id” in the training data and validation datasets were used in the deep learning model with dynamic data and the machine learning model with static data. A schematic illustration of this study is shown in Figure 1A.

The LSTM is a special modified recurrent neural network (RNN) that is widely used for the prediction of time series data in ICU due to the property of handling previous information for a relatively long period of time. Based on the 1 h high-resolution characteristics of ICU data, the prediction performance of LSTM has been fully developed over time, making it one of the most state-of-the-art deep learning algorithms for dynamic data in ICU. In this study, we processed the dynamic VIS data with the LSTM model with PyTorch version 1.8.0 accelerated by NVIDIA GeForce RTX 2060 SUPER graphics processing units (ZOTAC, Hong Kong). We set 200, 64, 0.0001, 0.25 as parameter “epochs,” “batch size,” “learning rate,” “dropout” for LSTM model, respectively. All the subjects of time series data needed the same length to train the model in batches; however, enrolled patients may have had inconsistent ICU stays within 48 h after initiating vasoactive–inotropic therapy. We referred to Thorsen-Meyer et al.’s method for processing dynamic time series data (21), and patients who did not match the given time point (1–48 h) were removed to generate training data with the same length.

The XGBoost is an optimized version of gradient boosted decision trees to implement classification and regression predictive modeling problems (25). The XGBoost algorithm made the following improvements under the gradient boosting framework: regularized boosting function to avoid overfitting, handling of missing values automatically, cross-validation at each iteration, incremental training, and tree pruning. Based on the performance and scalability mentioned above, the XGBoost algorithm is becoming increasingly popular in big data for critical care. We utilized Python version 3.8.8 to train an XGBoost machine learning model to predict ICU mortality with a total of 15 features, including the VRR and static maximum data mentioned above. We set ‘logloss’ as parameter “eval_metric” for XGBoost model. We applied the SHapley Additive exPlanations (SHAP) algorithm (26) in our XGBoost model to obtain explanations of the features that dominate ICU mortality predictions of patients with septic shock.

We divided patients with septic shock into 4 groups according to the VRR value: 50% ≤ VRR, 0 ≤ VRR < 50%, −50% ≤ VRR < 0, and VRR < −50%. Univariate Cox analysis was first performed to identify candidate variables that were considered to be clinically relevant for further multivariate Cox analysis, with a value of p < 0.05 as the cut-off value. We performed survival analysis with the Kaplan–Meier estimator and multivariate Cox proportional hazards regression analysis among these 4 groups. To further control the influence of confounding factors, we also conducted pairwise propensity score matching (PSM) (27, 28) among these 4 groups to handle the missing data and avoid bias. We set 500 and “pmm” as parameter “seed” and “method” in mice function for multiple imputation with mice package in R before PSM for pairwise primary cohorts. PSM was performed with the MatchIt package, and a total of 6 matched cohorts were generated.

We performed the Anderson–Darling normality test to determine whether the data were normally distributed and Bartlett’s test (original cohort, 4 groups) or F test (pairwise cohort, 2 groups) to assess the equality of variances. If the data among groups followed normal distributions and the variances to be compared were homogeneous, we utilized one-way ANOVA (original cohort, 4 groups) or t-test (pairwise cohort, 2 groups) to test the differences for continuous covariates; otherwise, Kruskal–Wallis test (original cohort, 4 groups) or Wilcoxon test (pairwise cohort, 2 groups) was applied as appropriate. The Chi-square test was used to test the differences for categorical covariates. All the statistical approaches mentioned above were executed in R version 4.0.2.

To obtain a robust conclusion, for these 6 matched cohorts, we applied a series of sensitivity analyses including the Kaplan–Meier (K-M) estimator, multivariate Cox regression analysis, and two doubly robust estimation models: the survey-weighted generalized linear model with all covariates using the inverse probability-weighting (IPW) technique calculated from PSM and the survey-weighted Cox model with all covariates using IPW for ICU mortality and in-hospital mortality separately.

We referred to a previous well-designed research (29) on sepsis to define the covariates in our study: (1) demographic information(age, gender, weight); (2) severity scores including simplified acute physiology score (SAPS) II, sequential organ failure assessment (SOFA) score, first nonzero VIS; (3) comorbidities including atrial fibrillation (AFIB), coronary artery disease (CAD); Congestive heart failure (CHF), chronic obstructive pulmonary disease (COPD), liver disease, malignant tumor and chronic renal disease; (4) vital signs including heart rate, mean arterial pressure (MAP), temperature; (5) interventions including mechanical ventilation, sedative therapy; and (6) laboratory tests including white blood cell (WBC) count, hemoglobin, platelet count, sodium, potassium, bicarbonate, partial pressure of oxygen (PO2), partial pressure of carbon dioxide(PCO2); central venous pressure (CVP), etc.

To clarify whether early VIS data could predict patient prognosis, we extracted the first 48 h VIS data after initiating vasoactive–inotropic therapy into dynamic and static VIS data. As the schematic illustration shows (Figure 1A), dynamic VIS data are a real-time, 1 h-resolution data matrix containing 7 features: the VIS of dopamine, dobutamine, epinephrine, milrinone, vasopressin, and norepinephrine, and the total VIS. Static VIS data were characterized with 15 features, including VRR and VIS features at 1–24 h and 25–48 h, respectively, which were extracted based on the dynamic VIS data.

We applied dynamic and static VIS data to train the LSTM deep learning model and XGBoost machine learning model to predict ICU mortality, respectively. The LSTM model was trained with temporal data, which was capable of updating the prediction of ICU mortality based on learning hourly VIS characteristics accumulated over time. Unsurprisingly, the prediction accuracy of mortality by the LSTM model trained with 7 VIS features was increasingly improved to some extent with the extension of the time dimension of the VIS data (Figure 2A). In addition, the improvement of prediction accuracy has reached a bottleneck to a certain extent over time, which may be related to too few features being included. Compared with the LSTM model, which requires extensive and high-resolution temporal data, the XGBoost model trained with less static VIS data had a similar model performance for the prediction of ICU mortality (Figure 2B). VRR simplifies the dosing trends of vasoactive and inotropic agents represented by dynamic VIS data while using less data. Although we converted dynamic data to static data, VRR represented a significant relationship with patient mortality in terms of model predictive performance. Since the VRR is calculated by the VIS data in two consecutive time periods, it still represents the dynamic trend of the dose of vasoactive drugs to some extent. Through the prediction model established by our work above, we were able to predict ICU mortality to some extent by relying solely on VIS data, which led us to wonder which VIS feature made an essential contribution to the predictive output of the model. Next, we applied the SHAP algorithm quantitatively and visually explained the XGBoost model (Figure 2C).

The SHAP summary plot displayed how the top features of VIS data impact our XGBoost model output for ICU mortality. Order on the y-axis represented the magnitude of the average impact of VIS features on our model (Figure 2C; Supplementary Figure 1). The colorful dots from blue to red indicated the value of continuous VIS data from low to high (Figure 2C), and the position on the left or right side of the x-axis implied whether the prediction toward ICU survival (right side: favor to ICU survival, left side: favor to ICU non-survival). The top 6 impactful VIS features of the XGBoost model for ICU mortality were VIS_{25-48h max}, VIS_{25-48h max} of vasopressin, VRR, VIS_{1-24h max} of norepinephrine, and VIS_{25-48h max} of norepinephrine.

Unsurprisingly, VIS_{25-48h max} had the most significant role in our trained model, implying that the higher the VIS during 25–48 h is, the higher the tendency toward non-survival. However, the contribution of VIS_{1-24h max} to the model output was slightly more complex but also made sense: there are some individuals on the left side with higher VIS_{1-24h max} values, but the prediction still favors survival, and a similar paradox exists on the right side of VIS_{1-24h max}. The seemingly contradictory predictions can be explained by some clinically common phenomena, such as patients in ICU with higher VIS at the beginning of ICU admission but whose hemodynamics improve the next day as the treatment progresses so that the dose rate of their vasoactive agents is markedly reduced. Therefore, the dynamic change in VIS highlights the predictive role in the model; it is hardly accidental that the higher the VRR is, the greater the tendency toward survival in ICU will be. The important role of VIS (norepinephrine) at 1–24 h and 25–48 h in the model was also confirmed: norepinephrine is the first choice in hemodynamic support therapy for patients with septic shock, and the effects of the VIS (norepinephrine) were similar to those of VIS_{25-48h max} and VIS_{1-24h max}. In addition to norepinephrine, the max dose rate of vasopressin during 25–48 h also had impactful features of model output, which conformed to the recommendations for vasoactive agents in the latest Surviving Sepsis Campaign Guidelines from 2021 on the one hand, but also indicates that a large dose of vasopressin in late stage is not conducive to the prognosis of ICU outcome.

In general, we successfully illustrated a significant association between VIS data and ICU mortality with a machine learning algorithm with small amounts of static VIS data and a deep learning algorithm with large amounts of dynamic data. We quantified and visualized the dynamic change in VIS data represented by VRR, which is an important indicator of patient prognosis in ICU.

In view of VIS unified 6 widely used vasoactive agents in ICU under one metric, and its spatiotemporal dynamic change was significantly associated with the prognosis of ICU patients, we divided the patients of the cohort into four groups according to the VRR value (Figure 1B): 50% ≤ VRR (n = 5,743), 0 ≤ VRR < 50% (n = 2084), −50% ≤ VRR < 0 (n = 593), and VRR < −50% (n = 476). The basic demographic characteristics of the original cohort are shown in Supplementary Table 1. As shown, the baseline values of some pivotal variables, such as the SAPS II, SOFA score, and first nonzero VIS were significantly imbalanced. Then, we combined the four cohorts for each pairwise comparison, resulting in six new sets of cohorts (cohorts 1–6). To avoid baseline imbalances between groups and the bias caused by missing data, we performed multiple imputation for missing data and matched the imputed cohorts with 1:1 pairwise PSM. Table 1 shows the baseline characteristics of the six cohorts after 1:1 pairwise PSM: cohort 1 (50% ≤ VRR versus 0 ≤ VRR < 50%, 2046 pairs), cohort 2 (50% ≤ VRR versus −50% ≤ VRR < 0, 592 pairs), cohort 3 (50% ≤ VRR versus VRR < − 50%, 475 pairs), cohort 4 (0 ≤ VRR < 50% versus −50% ≤ VRR < 0, 592 pairs), cohort 5 (0 ≤ VRR < 50% versus VRR < −50%, 473 pairs), and cohort 6 (−50% ≤ VRR < 0 versus VRR < −50%, 452 pairs). The covariate characteristics of each cohort before and after PSM are shown in Supplementary Tables 2–7. As shown there, the covariates with standardized mean differences (SMDs) less than 0.1 in the most matched cohorts increased significantly compared with those before matching, indicating that the baseline of the matched cohort was better balanced.

As vital signs, laboratory tests, severity scores, complications, and other covariates were greatly influenced, adjustments were made for these variables in the multivariate Cox regression. The K-M curve and multivariate Cox regression analysis showed a significant beneficial effect in the original cohort with increased VRR in ICU mortality (Supplementary Tables 8, 9; Supplementary Figure 2). In the 6 pairwise matched cohorts after multiple imputation and PSM, multivariate Cox regression analysis revealed patients in group 0 ≤ VRR < 50% [hazard ratio (HR) = 1.32, 95% confidence interval (CI) 1.17–1.50, p < 0.001], −50% ≤ VRR < 0 (HR = 1.79, 95% CI 1.44–2.22, p < 0.001), and VRR < −50% (HR = 2.07, 95% CI 1.61–2.66, p < 0.001) had a higher risk of mortality compared with 50% ≤ VRR (Figures 3A–C; Supplementary Tables 10, 11); patients in group −50% ≤ VRR < 0 (HR = 1.41, 95% CI 1.17–1.70, p < 0.001), and VRR < −50% (HR = 1.66, 95% CI 1.34–2.06, p < 0.001) had a higher risk of mortality compared with 0 ≤ VRR < 50% (Figures 3D,E; Supplementary Tables 13, 14); patients in group VRR < −50% (HR = 1.23, 95% CI 1.01–1.50, p < 0.05) had a higher risk of mortality compared with −50% ≤ VRR < 0 (Figure 3F; Supplementary Table 15). The K-M estimator of the matched cohorts showed same results (Figures 4A–F). In summary, ICU mortality was increased with decreased VRR values.

The results for matched cohorts of in-hospital mortality were similar to ICU mortality. Survival analysis indicated patients in group 50% ≤ VRR showed the lowest in-hospital mortality in the original cohort (Supplementary Tables 16, 17; Supplementary Figure 3). For the 6 matched cohorts, multivariate Cox regression analysis revealed significant differences (p < 0.001) in pairwise comparisons between any of the groups (Figures 5A–E; Supplementary Tables 18–22), except for VRR < −50% versus −50% ≤ VRR < 0 (Figure 5F; Supplementary Table 23), which had no statistically significant difference for in-hospital mortality (HR 1.16, 95% CI 0.96–1.39, p = 0.12). The K-M estimator of the matched cohorts showed similar trends (Figures 6A–F). Although in-hospital mortality was significant non-difference in matched cohort 6, it should be noted that VRR values of both two groups in cohort 6 were negative, implying that once the intensity of hemodynamic support therapy quantified by VIS decreased during 25–48 h compared with 1–24 h (namely VRR > 0), it was inclined to be in lower in-hospital mortality.

To ensure the robustness of the results, as shown in Table 2, two doubly robust estimation models: the survey-weighted generalized linear model and survey-weighted Cox model with all covariates using IPW, were also conducted for ICU and in-hospital mortality in each matched cohort. These findings indicated that an increase in VRR had a significant beneficial effect on survival in ICU or hospital. In terms of ICU mortality and in-hospital mortality, the results of two doubly robust estimation models were similar to those of multivariate Cox regression for in-hospital mortality. Except for matched cohort 6 (VRR < −50% versus −50% ≤ VRR < 0), there were statistically significant differences (p < 0.05, Table 2) for ICU mortality and in-hospital mortality between any other groups in the survey-weighted generalized linear model [odds ratio (OR) 1.24, 95% CI 0.94–1.65, p = 0.13; OR 1.06, 95% CI 0.80–1.41, p = 0.7, respectively. Only for cohort 6] and survey-weighted Cox model [OR 1.16, 95% CI 0.96–1.42, p = 0.13; OR 1.06, 95% CI 0.89–1.28, p = 0.5, respectively. Only for cohort 6]. This result may be due to the higher mortality of the two groups in cohort 6. Multivariate Cox regression models for the six matched cohorts mentioned above are also presented in Table 2. All the three models for matched cohort 1 to cohort 5 came to the same conclusion: patients with a higher positive VRR value were highly associated with lower ICU mortality and in-hospital mortality.