Deep Learning Unlocks the True Potential of Organ Donation after Circulatory Death with Accurate Prediction of Time-to-Death

Modeling Longitudinal Clinical Variables Acquired at Irregular Time Intervals

Our model uses the last 24 hours of clinical observations of a patient prior to terminal extubation, which contains both static and longitudinal variables. The latter pose challenges for statistical analysis and machine learning methods due to their irregular measurements over time and the presence of missing values. We therefore used an Ordinary Differential Equation Recurrent Neural Network (ODE-RNN)¹⁸, a recent state-of-the-art deep learning architecture that addresses both issues. ODE-RNN integrates a recurrent neural network (RNN)¹⁰ component tailored for sequential modeling, with a Neural ODE²¹ component that interpolates between irregularly-sampled time points.

The model operates by accumulating longitudinal variables with static variables to create a summary of the clinical history of each patient, that we call the latent phenotype²². The latent phenotype is then used by a classifier to predict patient outcomes (i.e. TTD). This procedure makes the ODE-RNN particularly effective at processing EHRs containing both static and longitudinal variables²³. A graphical depiction of the architecture is presented in Figure 1. Further details of the model are described in the Methods section.

• Download figure
• Open in new tab

Figure 1. Description of the problem setup and model architecture.

A. Problem setup. Based on the static variables of a specific patient (e.g. age, sex) and the last 24 hours of clinical follow-up prior to extubation (e.g. SpO2, MAP), our model predicts 4 probabilities: the probability that the time-to-death (TTD) is shorter than 30 minutes, between 30 and 60 minutes, between 60 and 120 minutes, and longer than 120 minutes. The sum of these probabilities equals to 1 by design. BMI stands for body mass index, SpO2 for oxygen saturation, MAP for mean arterial blood pressure, Hgb for hemoglobin, and NE for norepinephrine. Note that we consider 5 static variables and 25 longitudinal variables, and only some are shown for illustration purposes. B. Architecture of our ODE-RNN. The set of variables fed to the model consists of a concatenation of the longitudinal variables available at that observation time and a mask specifying which longitudinal variables are observed. Each clinical observation is sequentially processed by a gated recurrent unit (GRU) that incorporates the observation into the hidden state representation from the previous samples in the time series. Between observations, an ordinary differential equation (ODE) models the evolution of the patient’s hidden state continuously over time, which enables processing of variable temporal intervals between subsequent observations. The hidden state obtained after the whole time series is then complemented with the static variables to form the latent phenotype, a vector representation that summarizes the whole available information about the patient. The end classification is performed by using a multi-layer perceptron classifier (MLP) that predicts the TTD probabilities from the latent phenotype.

Predictive Performance Evaluation

We compared our method with the UNOS criteria, the most widely used clinical score for identifying DCD candidates⁸, and existing machine learning models that have been used for the prediction of clinical outcomes, including RNN¹⁰, LSTM¹¹, GRU¹², GRU-D¹², and XGBoost⁹.

All machine learning models were trained on the Yale New Haven Hospital (YNHH) cohort using a temporal data split. Data from patients before 2021 was used for training the models. Patients after 2021 were used for evaluation only to ensure robustness to distribution shifts over time.

The models were trained to predict time-to-death as a categorical variable, i.e., whether TTD fell within a given time frame (0-30 min, 30-60 min, 60-120 min, or >120 min). We evaluated the different models according to the overall categorical accuracy as well as pairwise binary classification for different grouped time frames (e.g., <30 min vs. >30 min). For these binary groupings, we also computed the area under the positive and negative predicted values (PPV, NPV), the area under the receiver operating characteristic curve (AUC-ROC) and the area under the precision-recall curve (AUC-PR). To assess the calibration of the models, that is, how well the predicted values represent the true likelihood, we computed the expected calibration error (ECE)²⁴.

Table 1 displays the comparative performance of various models on the YNHH patient cohort after 2021 and the external validation cohort. We found that ODE-RNN based model consistently outperformed other methods on all metrics, for TTD prediction at 30, 60 and 120 minutes. We note the high performance of the model despite the stringent experimental setup (temporal split and external validation), highlighting the robustness of the method. ODE-RNN also shows the best calibration, suggesting the probability outputs of the model are very reliable.

The poor performance of XGBoost and UNOS can be explained by (1) the inability of UNOS criteria and XGBoost to capture the temporality of the patient’s data, i.e. they only use the last observation at the time of extubation; (2) the limited number of clinical variables used in UNOS (14 variables) compared to our ODE-RNN model (5 static and 25 longitudinal variables).

In Figure 2 panel B, we report the calibration plot of the ODE-RNN model for the binary prediction (< 30 min vs. > 30 min) on the external validation cohort. Calibration plots for other binary tasks are available in the Supplementary Materials. The model tended to give a reliable but conservative estimate of the probability of death within 30 minutes. Importantly, the model appeared well calibrated for low and high predicted probabilities, highlighting its reliability.

• Download figure
• Open in new tab

Figure 2. Model performances and analyses.

A. Graphical representation of performance assessment for the UNOS score, XGBoost, LSTM, and ODE-RNN across different binary tasks. Left: TTD<30 min vs. TTD>30 min; Middle: TTD<60 min vs. TTD>60 min. Right: TTD<120 min vs TTD>120 min. We report the binary accuracy, the area under the receiver operating characteristic curve (ROC AUC), and the precision recall curve (PR AUC). ODE-RNN outperforms all other models for all tasks and all evaluation metrics. B. Calibration plot for the binary classification task TTD <30 min vs. TTD >30 min on the external validation cohort, computed with the R package val.prob.ci.2. The predicted probabilities come from the output of our model and plotted against the fraction of positives observed in the data. The histogram shows the prevalence of patients for different ranges of predicted probabilities. C. Variable importance for the predictions of the ODE-RNN model in the binary classification (left: <30 vs >30 minutes, middle: <60 vs >60 minutes, right: <120 vs >120 minutes). Variable importance was computed using permutation importance testing. The input variables were split into static variables (that do not change over time) and longitudinal variables. ROC-AUC stands for area under the receiver operating characteristic curves.

Variable Importance Assessment

We assessed the importance and impact of the different clinical variables on the prediction of the models with permutation importance testing (Figure 2 panel C). For longitudinal variables, we found that heart rate was the most important variable in the prediction of the ODE-RNN, followed by respiratory rate, mean arterial blood pressure (MAP), oxygen saturation (SpO2), and the Glasgow Coma Scale (GCS) score. Corneal reflex and gag reflex which are frequently used by transplant surgeons and OPOs, appear to have the least impact on the prediction. Notably, static variables were found significantly less important than longitudinal variables, by an order of magnitude. We also found a strong consistency in the variable importance across different binary tasks.

Patient Phenoscape Analysis

The patient phenotypes learned by our model enable accurate TTD predictions because they faithfully represent the patients’ condition and past clinical history. As such, these phenotypes provide richer information about the patients, compared to the single numerical value of TTD prediction. These phenotypes form a continuous landscape of the patient cohort, which we refer to as the phenoscape. Within the phenoscape, we can observe clusters of patients with similar conditions and continuous transitions from one condition to another, thereby uncovering the underlying dynamics of circulatory death. We used PHATE²⁵, a dimensionality reduction method that preserves the underlyding data geometry^26–29, to visualize the phenoscape and provide examples of new insights drawn from such analysis.

Our patient phenoscape visualizations in Figure 3 revealed that patients reside on a continuous spectrum of phenotype that goes beyond the TTD categorization. Patients were organized along an axis that corresponded with TTD but also with heart rate, SpO2, or GCS, among others. Finer investigation allowed us to identify the dynamical patterns most correlated with TTD, complementing the variable importance analysis above. For instance, in Figure 3 panel A, patients were colored according to their average heart rate and to their range of heart rate measurements (defined as the difference between highest and lowest values). While the range correlated with TTD, the average value did not, suggesting the variation in the heart rate is more important than the average.

• Download figure
• Open in new tab

Figure 3. Visualization of the patient phenoscape.

The latent phenotype is visualized in two dimensions using PHATE. In this plot each point represents a patient and the coloring is based on the value of different clinical variables. A. All patients in the Yale cohort are plotted and colored according to TTD label (top-left), TTD (log-transformed, top-right), heart rate (middle), SpO2 (bottom-left), and GCS (bottom-right). Each point represents a patient. These plots uncover the continuous structure of the patients’ latent phenotype and highlight the correlation between different clinical variables and the time-to-death. The longitudinal variables were transformed into scalar variables using different transformations. (range) computes the average of the five highest observations minus the average of the five lowest observations in the clinical history of the patient. (last) computes the average of the last five observations in the clinical history of the patient. (min) computes the average of the five lowest observations. Different transformations extract different patterns from the time series, enabling a finer interpretation of the dynamical patterns for a given clinical variable. For instance, we observed that heart rate (range) correlates with the TTD label, unlike heart rate (mean), suggesting the variation in heart rate is more important than the average value. The visualization of the whole cohort suggests two distincts groups of patients, characterized by high or low TTD. B. Focus visualization of the identified cluster of patients with TTD<120 min. This uncovers a finer grained structure in this specific cohort of patients. We colored the patients by TTD (log-transformed, top-left), range of heart rate (middle-left), minimum SpO2 (middle-right), average GCS (bottom-left), and BMI (bottom-right). C. We clustered the patients in the zoomed-in group of patients according to the similarity of their latent phenotype. We obtained three clusters: A, B, and C. Guided by the visualization of panel B, we examined the specific phenotype of patients from cluster A, which appear to have higher TTD than the rest of the patients. We show boxplots and corresponding independent t-tests p-values for difference of means between clusters, for various clinical variables. This analysis characterizes cluster A as a subgroup of patients with high TTD, high range of heart rate, low minimum SpO2, high GCS, as well as low BMI.

Our phenoscape visualization also showed a clear distinction between patients with TTD<120 minutes and another group with TTD>120 minutes (categories 0,1,2 vs. 3 in top-left of Figure 3 panel A). This separation suggests an obvious clinical difference between patients with TTD<120 min and TTD>120 min. We performed a fine-grained analysis of the former cluster, to identify clinical drivers that make the difference between short-range (≈30min) and medium-range (≈60min) TTD. Panels B and C in Figure 3 show a visualization of the identified group of patients with TTD<120min. From this zoomed-in analysis, we identified three clusters (A, B, and C). Guided by our visualization, we examined the specific phenotype of patients from cluster A, which appear to have higher TTD than rest of the patients. We found that patients from cluster A were characterized by higher TTD, higher range of heart rate, lower minimum SpO2, higher GCS, higher MAP and lower body mass index (BMI). This reveals that, over the last 24 hours before TE, range of heart rate, minimum SpO2, average GCS, and the maximum MAP observed are all predictive of TTD.

Patient cohort and data preparation

We used two separate cohorts of patients to develop and validate the model. The first cohort contained 3,238 patients at Yale New Haven hospital (YNHH) older than 18 years old, with a recorded TE in the ICU between 2014 and 2023. For unbiased validation, using same inclusion criteria, we formed an external cohort from six different hospitals with 1,908 patients. The hospitals from the external cohort are the Bridgeport hospital, Greenwich hospital, Lawrence and Memorial hospital, Saint Raphael hospital, Westerly hospital and the Yale New Haven Children’s hospital. The median time from extubation to death was 7.15 minutes in the first cohort and 8.28 minutes in the external cohort. The two cohorts are summarized in Table 2.

For each patient, we extracted 5 static variables and 25 longitudinal variables, collected up to 24 hours before extubation. Longitudinal variables were observed at a range of frequencies – from a single observation in 24 hours to one/minute. Missing values in the longitudinal records were imputed by performing a combination of forward-fill, backward-fill and mean-fill imputation⁴⁰. In addition to imputation, presence of missing values was fed to the model by creating binary missingness indicators. TTD was defined as the time from terminal extubation to circulatory death. TTD was converted into an ordinal variable with 4 categories (0: 0∼30 minutes, 1: 30∼60 minutes, 2: 60∼120 minutes, and 3: longer than 120 minutes), that were used as target labels in the machine learning model. We chose these time frames as most transplant centers use 0-60 minutes as “standard criteria” for kidney DCD donation and 60∼120 minutes as an “extended spectrum”⁴¹; similarly, 0∼30 minutes as “standard criteria” for liver DCD donation and 30∼60 minutes as an “extended spectrum”⁴².

List of clinical variables used in the models

The following longitudinal and clinical variables were extracted for each patient.

Model description

To capture the impact of longitudinal variables on the target while accommodating for the irregular sampling of clinical trajectories, we used an Ordinary Differential Equation Recurrent Neural Network (ODE-RNN)¹⁸. ODE-RNNs combine two powerful architectures, recurrent neural networks (RNN) and neural ordinary differential equations (Neural ODE), making them exceptionally apt at processing clinical time series^19,23. RNNs are neural networks specialized for processing sequences. At each time step, they maintain a hidden state which represents the whole previous information in the time series. Upon reading a clinical record, the RNN updates its hidden state by combining the previous hidden state with the new observation, using an update unit (here a gated recurrent unit (GRU)). However, RNNs assume continuous time intervals between observations, an assumption typically not met in clinical time series. To address this limitation, ODE-RNN uses a Neural ODE, that describes dynamics in continuous time, to model the dynamics of the hidden state between observations. This uniquely allows the model to capture the full span of the clinical records of each patient, correctly accounting for the time interval between observations, and improving upon previous methods such as logistic regression, XGBoost, or the UNOS criteria, among others.

The model accumulates the last 24 hours before extubation of the patient’s vitals, medications used (e.g. vasopressors), neurological assessments, and lab results, together with their demographic records, and produces a representation of the patient, which we refer to as the patient’s latent phenotype. This latent phenotype can be understood as a learnt compact clinical summary of a particular patient. This representation is then used as an input to a multi-layer perceptron classifier that predicts the probability of each label category (0, 1, 2, or 3, corresponding to the 4 time ranges). Figure 1 depicts our model pipeline.

Our model contains the following neural networks:

• f_static: a multi-layer perceptron (MLP) that processes the static variables.

• f_ODE: an MLP that predicts the derivative of the hidden state dynamics between observations.

• f_GRU: a gated recurrent unit (GRU) that updates the hidden states at each observation point.

• f_longitudinal: an MLP that processes the final hidden state of the longitudinal variables.

• f_fusion: an MLP that fuses the latent states derived from static and longitudinal variables, producing a latent variable called the latent phenotype.

• f_classifier: an MLP that performs classification using the latent phenotype.

Algorithm 1 describes how the model predicts outcomes for a single patient. The model takes as input static variables s ∈ ℝ^ℓ, longitudinal variables x₁, …, x_n ∈ ℝ^k, observation times t₁, …, t_n ∈ ℝ, and boolean observation masks m₁, …, m_n ∈ ℝ^k, where m_i = (m_i1, …, m_ik). The value m_ij = true if x_ij, the j^th variable at time i, is observed, and m_ij = false otherwise. Including these observation masks enables the model to utilize “informative missingness,” which is correlated with the patient’s condition. For instance, a patient is unlikely to be suspected of heart failure if B-type Natriuretic Peptide (BNP) is not frequently measured. The model iterates through the observation time points, updating the hidden state h of the longitudinal variables using f_ODE and f_GRU. Between observation points, f_ODE is integrated to continuously update h, while at observation times, f_GRU updates h using x_i, m_i, t_i. The final h contains accumulated information from the entire history of the longitudinal variables, which is then processed by f_longitudinal and fused with the processed static variables s (via f_static) using f_fusion. This results in a latent variable representing the patient’s condition, referred to as the latent phenotype. f_classifier uses the latent phenotype to make the final classification prediction.

Algorithm 1

ODE-RNN using GRU cell update (using one patient for illustration)

• Download figure
• Open in new tab

Evaluation of model performance

We compared the performance of our approach with machine learning methods and clinical scores such as the UNOS criteria⁸. Machine learning methods directly learn associations from the available data while clinical scores consist of clinical criteria designed by experts. Within the machine learning methods, we distinguish between static methods and longitudinal methods. Unlike longitudinal methods, static methods cannot process a time series of clinical information. They are therefore trained on the last available clinical observation at the time of extubation only.

Visualizing structures in high-dimensional patient phenoscape with PHATE

To predict TTD, our ODE-RNN model produces a latent phenotype for each patient, which can be intuitively understood as a learnt summary of the clinical history of the patient20. We explored the space of latent phenotypes of all patients in the cohort, the patient phenoscape, showing its potential to provide new clinical insights.

The latent phenotype of each patient is high-dimensional, and thus cannot be directly visualized. Therefore, we first produced a two-dimensional representation of the phenotype using PHATE²⁵, a non-linear dimensionality reduction and visualization method that stays faithful to the geometry of the data and retains the inherent similarity between patients. In this visualization, each patient is represented as a point in a two-dimensional phenotypic space. The patient phenoscape is the set of the representations of all patients in the cohort. We colored each point according to their TTD and the value of certain clinical variables, enabling a fine-grained exploration of the impact of different clinical factors on the TTD.