Population scale proteomics enables adaptive digital twin modelling in sepsis

2.1 Predicting sepsis patients from non-sepsis mimics

In the cohort, around 70% of the patients had sepsis or septic shock while around 10% of the patients had organ dysfunction unrelated to sepsis (Figure 1d). A uniform manifold approximation and projection (UMAP) of all quantified proteins in the 1364 proteome maps showed little separation between non-sepsis and sepsis patients (Figure 1h). Statistical comparison between the two groups revealed several differentially abundant proteins that were used to train a classifier to predict sepsis with a ROC-AUC of 0.71 based on the test set (Supplementary Figure S1 and Figure 1i). Repeating the same process for septic shock, generated a protein panel that performed 33.8% better (0.95 ROC-AUC) than the original sepsis classifier on a filtered test set (non-sepsis vs. septic shock), but worse on the complete test set (Supplementary Figure S1 and Figure 1i). A possible explanation for the relatively low performance of the classifiers is the heterogeneity within the cohort groups (1-5). To test this, we rescored all sepsis patients using the septic shock protein panel from the training set and separated them into groups of high probability severe sepsis (>0.5 probability) and low probability severe sepsis (≤0.5 probability) (Figure 1i). The more severe sepsis patients correspond to lower survival than the less severe sepsis group (Figure 1k), indicating that this strategy could be useful in identifying patients with potentially worse survival using proteomic markers from time-of-admission samples. Collectively, these results demonstrate that the use of all quantified proteins generates relatively poor separation between the cohort groups. One possible explanation for this is that the predominant signal picked up from the high-dimensional proteomics data is the unspecific inflammatory response that does not separate the patients into discrete subgroups.

2.2 Organ dysfunction panels to construct a low-dimensional latent space

To stratify the cohort more effectively, we speculated that reducing the high-dimensional proteomics data into a low-dimensional interpretable embedding could reveal actionable patterns within the data. To that end, we mirrored the diagnosis of sepsis computationally by investigating the proteomic signatures derived from different SOFA categories (i.e. respiratory, coagulation, liver, renal, cardiovascular, and central nervous system (CNS)), the presence of a verified infection, and sepsis from non-sepsis. For these categories, we identified protein panels using differential analysis and XAI to select for proteins specific to the clinical outcome they were predicting (Figure 2a). The selected panels of proteins could predict future clinical diagnosis from time-of-admission samples with areas under the curve (AUC) ranging from 0.57 to 0.90 (Figure 2b). Aggregating all proteins from the individual protein panels (Figure 2c) generated a combined protein-organ dysfunction network comprised of 65 unique proteins (Figure 2d). Most proteins were singularly identified for a panel, indicating a strong association with their specific category. Two notable examples are the highly lung-specific pulmonary surfactant-associated protein B (PSPB) that was found most important for predicting respiratory dysfunction, and cystatin-c (CYTC), important in renal dysfunction, which has previously been identified as a blood-based biomarker for acute kidney injury [31]. Other examples, such as C-reactive protein (CRP), a common biomarker for infection and inflammation, and lipopolysaccharide-binding protein (LBP), a biomarker for infection [32], were both identified as important in predicting infection. Functional analysis with STRING-DB [33] revealed several significant interactions and biological processes statistically associated with the selected protein panels (Supplementary Figure 2). For example, pathways involved in lipid metabolism were found enriched for the liver panel, coagulation and platelet regulatory pathways were found enriched for the coagulation panel, and opsonization and the response to bacteria were found enriched for the infection panel. The biological relevance between the selected proteins and their corresponding categories emphasizes the robustness of our feature selection method, as the selected proteins are tied to biological functions of the different categories rather than just common inflammatory pathways. In the final step, we used these individual protein panels to predict probabilities for every patient in the training data to create an interpretable probabilistic latent space (ILS) that incorporates probabilities of organ dysfunction, infection, and sepsis to be used for further investigation into the molecular mechanisms of the syndrome (Figure 2e).

• Download figure
• Open in new tab

Figure 2: Explainable machine learning uncovers specific molecular panels in sepsis:

a All time-of-admission plasma samples from the training data were analyzed using LC-MS/MS to obtain proteome maps for each patient sample. From the global list of identified proteins, 8 different clinical markers used in the diagnosis of sepsis (i.e. organ dysfunction, infection, sepsis) were analyzed and proteins were selected using explainable machine learning that were the most descriptive for each proteome state. b ROC curves displaying the performance of the selected protein panels on the held-out test data set. c Bar plots showing the mean importance, calculated using SHAP, of the selected proteins after 100 bootstrap iterations. d All proteins shown in (c) displayed as a network graph. Edges are unweighted and based on the inclusion of a protein in a panel. e An overview of the workflow used to predict organ dysfunction, infection, and sepsis for all patients in the training cohort to create the interpretable latent space. Probabilities for each of the 8 panels are predicted for each patient in the training data to assemble an 8-dimension latent space where each dimension represents a probability for one of the 8 panels.

2.3 Identifying distinct subphenotypes using the interpretable probabilistic latent space

To determine if the interpretable latent space could reveal distinct subphenotypes, we visualized a UMAP of the ILS using only sepsis and septic shock patients (groups 4-5) from the cohort. In contrast to the full proteome, this UMAP revealed marked substructures of patients. We applied the k-means algorithm directly to the ILS to identify the 5 distinct subphenotypes within the data (Figure 3a). Each cluster was associated with distinct differences in predicted organ dysfunction (Figure 3b) and mortality rates (Figure 3c). Cluster 4 exhibited the best survival and was composed of high probability infection and sepsis but low probability severe organ dysfunction. Cluster 2 had high survival and was composed of high probability sepsis, infection, and respiratory dysfunction, but low probability severe organ dysfunction in the other categories. In contrast, clusters 0, 1, and 3 had worse survival and were all composed of high probability sepsis, infection, and multiple severe organ dysfunction (Figure 3b).

• Download figure
• Open in new tab

Figure 3: Identifying heterogeneous subgroups within sepsis:

a A UMAP of the 8-dimensional latent space, subdivided into 5 clusters using the k-means algorithm on the 8-dimensional latent space and colored by the assigned cluster. The same colors are used throughout the figure to represent the 5 clusters (0-4). b The percentage organ dysfunction, infection, and sepsis was determined for each clusters and plotted in a radar plot and colored according to clusters. c The survival profiles of the 5 clusters across 365 days. The dotted lines highlight survival at 30 and 90 days respectively. d The top-5 most important proteins and clinical parameters selected using explainable machine learning in a one-vs-all analysis for each cluster. The gray box plots represent the values for all other clusters combined. The boxes represent the quartiles of the data, while the whiskers represent the remaining distribution. Outliers are not shown. e Cluster maps for each of the 5 clusters and 8 clinical parameters displaying the raw predicted probabilities of the latent space shown using a green color scale. f A zoomed in UMAP of cluster 1 from (a) colored by 4 subclusters identified using hierarchical clustering. These colors are used for the same subclusters in g-h g A radar plot of the percentage organ dysfunction, infection, and sepsis for the 4 unique subclusters found within cluster 1. h Survival profiles for the 4 subclusters within cluster 1. Similar graphs for the other 4 main clusters shown in (a) are available in the supplementary material.

To understand the molecular and clinical differences between the 5 subphenotypes, we used the previously described XAI methods to identify the five most prominent proteomic and clinical features for each cluster using the data available on time-of-admission (Figure 3d). Interestingly, the patients in cluster 1, associated with high probability liver dysfunction, had higher levels of bilirubin, which is a biomarker for liver dysfunction, than the other clusters as well as higher levels of lactate. The proteins found important for cluster 1 were also found in the liver dysfunction panel (Figure 2c). Cluster 0, the cluster with the highest mortality and high probability of renal dysfunction, was associated with high creatinine levels and lower systolic blood pressure. Although there were distinct differences between the identified subphenotypes, it is clear in Figure 3d that there is considerable overlap of the clinical parameters and sharedness of selected proteins between the clusters. These results indicate that the subphenotypes are heterogeneous and their associated features are not sufficient by themselves to accurately stratify sepsis patients.

To further investigate the degree of heterogeneity within the clusters, we performed hierarchical clustering of the ILS for each patient in each subphenotype. The analysis revealed that each cluster was comprised of distinct subclusters of patients revealing an additional 30 subphenotypes from the original 5 subphenotypes (Figure 3e and Supplementary Figure 3). For example, hierarchical clustering revealed 4 subclusters within cluster 1 that form separate structures in the UMAP (Figure 3f). All 4 subclusters consist of predicted liver dysfunction but have different frequencies of other organ dysfunctions (Figure 3g) and mortality rates (Figure 3h). These results show that the ILS can be used to identify potential subgroups within the sepsis cohort. However, these subgroups are heterogeneous themselves and the number of distinctly unique subgroups could be set arbitrarily high. From these findings, we conclude that the use of subphenotypes for the precision treatment of sepsis patients is not optimal and a more dynamic method that can adapt as new data is acquired could lead to more actionable translational results.

2.4 Modelling new patient outcome with adaptive digital family analysis

As digital twin modelling has previously been shown to enable predictions of future clinical trajectories, we first leverage the population scale proteomics data to predict the ILS for the entire training cohort, building a database of 1091 virtual patient representations. In the second step, probabilities for the 8 panels used to create the ILS were predicted for the 273 previously unseen patients from the held-out test set to provide an unbiased evaluation of the approach. Using these probabilities, we selected 273 digital families each consisting of 2.5% of the nearest neighbors in the database (Figure 4a). These digital families enable the dynamic modelling of all 370 clinical parameters extracted from the cohort using time-of-admission proteomics data.

• Download figure
• Open in new tab

Figure 4: Digital family analysis facilitates accurately modelling of unknown disease trajectories.

a Using the proteome from time-of-admission samples of the training data and the detailed clinical data, the interpretable latent space (ILS) is predicted (Figure 2). This ILS is used as a database to determine groups of patients that are most similar to a new patient. Probabilities for new patients from the held-out test set were predicted for each of the 8 panels from Figure 2 and the nearest neighbors consisting of 2.5% of the database were selected as the digital family for each new patient. b A histogram of the error in predicted SOFA score, determined by taking the absolute value of the median SOFA score increase of the digital family subtracted from the actual SOFA score increase of the patient from the held-out test set. The dotted line represents the median SOFA error at 1.0. c A line plot showing the relationship between the predicted 30-mortality of the digital family, calculated as the percentage of how many in the digital family have survival less than 30-days, against the bootstrapped (n=1000) mortality rates in the test data. d Digital families were created for each patient in the held-out test set, and a UMAP was plotted of that test set using the same projection from Figure 3a. The colors indicate the predicted 30-day mortality of the digital family for each patient. The shaded colored outlines represent the 5 clusters from Figure 3. e Scatter plot of the 7-day survival ratio against the predicted 30-day mortality for each digital family with jitter added so all points are visible. Histograms for each axis are shown above and to the right of the plot. The coloring is based on the predicted cluster from Figure 3 of each patient in the held-out test set. f A selected survival profile for a patient with relatively high predicted 30-day mortality, but a higher 5-day survival gradient. This selected patient died within 5 days after admission to the hospital. Dotted circles show where this particular patient is placed in d and e.

The SOFA score increase is an important measure for poor prognosis but requires repeated measurements over multiple days. Here, we used the digital families to predict the increase in SOFA score after 24 hours for all 273 patients in the held-out test set and evaluated the accuracy. The digital families could be used to accurately predict the increase in SOFA score with a median error of 1 from time-of-admission samples, showcasing the ability to accurately predict future clinical outcomes for a patient (Figure 4b). Similar results were obtained when predicting the 30-day mortality for each patient within the held-out test set. Here, we calculated the 30-day mortality rate of each digital family and compared that rate to the bootstrapped (n=1000) mortality of test patients for each predicted mortality bin (Figure 4c). This analysis revealed a linear relationship between the predicated and observed mortality rate, although the model overestimated mortality rates in high mortality families. Additionally, by calculating the mortality for each patient in the held-out test set based on their digital families, we can effectively create a mortality map of sepsis (Figure 4d). The map shows that the high mortality rates are concentrated in the top of the high-mortality clusters we identified above, but that the mortality landscapes within each cluster are not evenly distributed. If a new patient is placed in an area of high mortality, it is possible to determine the mortality risk of a patient based on their digital family without placing them in a subphenotype first.

Interestingly, the UMAP in Figure 4d shows that many of the digital families have mortality rates below 10%. To investigate if these predictions could be further subdivided, we modelled the survival profiles for each family (Figure 4f) to identify families with low overall mortality rates, but with a high proportion of mortality within the first 7 days. For each family, the ratio of 7-day mortality over 30-day mortality (mortality ratio) was calculated and plotted against the 30-day mortality of the digital family (Figure 4e), revealing four distinct subgroups of digital families. One group, comprised of high 30-day mortality rates with over 60% mortality ratio, was populated mostly by patients from the high mortality cluster 4 subphenotype. Conversely, there were two groups with low 30-day mortality rate but with drastically different mortality ratios. We highlight a patient in Figure 4d-e who died in less than 3 days in the hospital that was associated with a digital family with an overall predicted mortality of 20% (Figure 4f). However, the family had a proportionally high mortality ratio where more than 80% of the mortality occurred within 7-days. Although future evaluation of these patterns is required to determine the most useful parameters for clinical decision making, these results demonstrate that the digital family model generates new opportunities to analyze and stratify patients previously unavailable using conventional approaches. In conclusion, the methods used here showcase how the complexity of sepsis can be modelled with the ILS and how time-of-admission samples can be used to predict future clinical trajectories with high accuracy using the adaptive digital family approach. The results here suggest that digital family modelling represents an alternative strategy to rule-based subphenotyping and can provide a more dynamic approach to patient stratification. This adaptive approach could prove highly effective in the clinic for future implementations of precision medicine strategies in sepsis.

4.1 Settings and patient population

Patients were eligible for inclusion if ≥ 18 years and admitted to the emergency department (ED) at Skåne University Hospital, a tertiary hospital in Lund, Sweden as sepsis alert. Sepsis alert is a modified triage system to detect patients with suspected sepsis in the ED. The criteria for sepsis alert are patients with any of the following: active seizures, unconsciousness, a respiratory rate higher than 30 breaths/min or lower than 8 breaths/min, oxygen saturation below 90%, a regular heart rate over 130 beats/min, or an irregular heart rate over 150 beats/min, a systolic blood pressure (SBP) below 90 mmHg, p-lactate >3.5 mmol in combination with fever, >38°C or a history of fever or chills [29]. Patients meeting the criteria for sepsis alert are cared for immediately by emergency physician, infection consultant, nurse, and assistant nurse and undergo immediate control of vital parameters, physical examination, blood sampling for biochemical analyses and culture and other microbiological sampling. Sepsis alert patients are also considered for other diagnostic procedures, treatment, level of care and surveillance.

Patients were included prospectively, consecutively from 1st September 2016 to 31st of March 2023 and venous blood samples were drawn in citrate tubes at admission. Samples were centrifuged and stored at −70°C within 2 hours of collection.

4.2 Ethical permission

Informed consents were collected through an opt-out procedure. Patients received a mail following the ED admission and were given the opportunity to opt-out of inclusion in the study by sending back a form in a pre-franked envelope. Ethical approval for the study was obtained from the regional ethics board (file numbers 2022-01454-01, 2014/741 and 2016/271).

4.3 Clinical data collection

Clinical data was collected retrospectively through clinical chart review. Chart reviews were conducted by medical researchers using a structured protocol and validated by infectious disease clinicians.

Data were collected on demographics, vital signs, laboratory testing, microbiological investigations, medical history, diagnostic and therapeutic procedures, concomitant or new medications and level of care.

There are 5 groups defined in the cohort overall: no organ dysfunction under the sepsis-3 definition [1] (< 2 increase in sequential organ failure (SOFA)) and no infection based on the LMCI score (1), organ dysfunction under the sepsis-3 definition (≥ 2 increase in SOFA) with no infection according to the LMCI score (2), no organ dysfunction under the sepsis-3 definition with infection according to the LMCI score (3), and sepsis (4) and septic shock (5) defined according to Sepsis-3 definitions with organ dysfunction defined as an increase in sequential organ failure assessment (SOFA) ≥2 and verified infection was determined using the LMCI score. SOFA was modified to be compatible with use outside intensive care units (ICU) (Supplementary Table T3).

The following comorbidities were recorded: cardiovascular disease, liver disease, malignancy, respiratory disease, diabetes mellitus, chronic kidney disease, immunodeficiency and Charlson comorbidity index.

Microbiological diagnostics were performed at the clinical microbiology department (Laboratory Medicine Skåne, Lund, Sweden) per standard clinical practice.

4.4 Sample preparation for LC-MS/MS analysis

Plasma samples were diluted 1:10 with 100 mM ammonium bicarbonate (Sigma-Aldrich), and a total of 10 µl of diluted sample (corresponding to 1 µl plasma) was used for protein digestion. Samples were incubated for 60 minutes at 37°C in 4 M urea (Sigma-Aldrich) and 60 mM dithiothreitol (DTT, Sigma-Aldrich) for denaturation and reduction. Addition of DTT to the samples, as well as the following additions, were performed with an Agilent Bravo liquid handler. The samples were then alkylated using 80 mM 2-iodoacetamide (Sigma-Aldrich) for 30 minutes at room temperature in the dark, followed by digestion with 2 µg LysC (Lysyl Endopeptidase, Mass Spectrometry Grade, Wako) for two hours at room temperature. The samples were further digested using 2 µg trypsin (sequence-grade modified porcine trypsin, Promega) for 16 hours at room temperature. The digestion was stopped with 10% trifluoroacetic acid until pH ∼2. The samples were stored at −80 °C until LC-MS/MS analysis. For peptide clean-up, digested samples were diluted 1:12 in water to generate suitable concentrations and subsequently loaded onto disposable Evotip Pure C18 trap columns (Evosep Biosystems, Odense, Denmark), which were prepared according to the manufacturer’s instructions. Briefly, the Evotips were activated in 0.1% formic acid in acetonitrile, conditioned by wetting the tips in 2-propanol and equilibrated in 0.1% formic acid. Twenty microliters of each diluted sample were transferred to tips, followed by washing with 0.1% formic acid. The Evotips were stored in 0.1% formic acid until LC-MS/MS analysis.

4.5 LC-MS/MS analysis

Evosep One LC system (Evosep Biosystems) was used for separation with nanoflow reversedphased chromatography. All samples were run with the 60 samples per day (SPD) method (gradient length of 21 min, vendor standard settings) with a flow rate of 1 µl/min. Samples were analyzed with timsTOF Pro 2 ion mobility mass spectrometer (Bruker Daltonics) using an 8 cm x 150 µm Evosep column (Evosep Biosystems) packed with 1.5 µm ReproSil-Pur C18-AQ particles. MS data was acquired using data-independent acquisition serial fragmentation (diaPASEF) method. The accumulation and ramp times were set to 100 ms. Variable isolation windows were used for the diaPASEF method and specified in Supplementary Table T2, with an estimated cycle time of 2.76 s. The collision energy was ramped linearly as a function of the mobility from 59 eV at 1/K0 = 1.5 Vs cm-2 to 20 eV at 1/K0 = 0.6 Vs cm-2.

4.6 LC-MS/MS data analysis

Data was analyzed using DIA-NN (1.8.1) [36] using a spectral library created from human tissue and plasma samples and match-between-runs enabled. Downstream analysis was performed using DPKS [21]. Samples were randomly split into a training (80%) and test (20%) set keeping the distributions of the included groups (1-5) intact and were processed separately using the following procedure to obtain protein quantities. Results were filtered for 1.0% false discovery rate at the precursor and library level. Precursors were normalized using a retention time sliding window mean normalization method [37] using DPKS. Proteins were quantified from precursors using the relative quantification iq [38] algorithm in the DPKS package.

4.7 Statistical analysis and explainable machine learning

All statistical analysis, explainable machine learning, and feature selection was performed using DPKS [21]. Statistical differential tests were performed between 2 groups using linear regression models for proteomics data, and ANOVA for clinical data. Multiple testing correction was performed using the Storey-Tibshirani method [39]. Features (proteins or clinical features) below a 0.1 corrected p-value were selected for downstream machine learning analysis. Feature importance was estimated using SHAP [40] with 100 bootstrapping iterations (n=100) to account for stochasticity in explainable machine learning methods caused by the background data using an XGBoost classifier [41]. The background data was resampled to correct for any class imbalances between groups during interpretation.

4.8 Sepsis prediction

Comparisons between groups of non-sepsis (1, 2, 3) (n=346) and sepsis (4, 5) (n=745) were performed using the described statistical and explainable machine learning methods. Protein ranks were calculated for 100 bootstrap iterations and the number of times a protein appeared in the top-20 was calculated. Proteins were selected if they were above 2 standard deviations from the mean number of occurrences for a protein in the top-20. An XGBoost classifier was trained to predict sepsis using the minimized protein panel. This classifier was validated using the unbiased held-out test set.

Comparisons between non-sepsis (1, 2, 3) (n=346) and septic shock (5) (n=67) were performed using the above described statistical and explainable machine learning methods. Protein ranks were calculated for 100 bootstrap iterations and the number of times a protein appeared in the top-20 was calculated. Proteins were selected if they were above 2 standard deviations from the mean number of occurrences for a protein in the top-20. An XGBoost classifier was trained to predict sepsis using the minimized protein panel. This classifier was validated using the unbiased held-out test set (n=273) and a filtered test dataset consisting of just the sepsis groups 1, 2, 3, and 5 (non-sepsis and septic shock).

4.9 Septic shock stratification

Probabilities for all patients in the training data were calculated using the septic shock classifier. The non-septic shock sepsis patients (4) were split into 2 groups (≤ 0.5 and > 0.5 probability) and their survival rates over 30-days were plotted.

4.10 Protein panel identification

Eight clinical outcomes were identified and analyzed to investigate the underlying proteins associated with them: 6 Organ dysfunction groups based on SOFA scores (respiratory, renal, liver, cardiovascular, coagulation, and central nervous system), presence of infection, and sepsis-3 (1, 2, 3 vs. 4, 5). This analysis was performed using the training data and validated using the unbiased test set. The following comparisons were made for the 6 organ dysfunction groups:

• Respiratory organ dysfunction: respiratory SOFA 0 (n=92) to respiratory SOFA >2 (n=102)

• Coagulation organ dysfunction: platelets >200 x 109/L (n=502) to platelets below 150 x 109/L (n=109)

• Liver organ dysfunction: bilirubin <20 µmol/L (n=556) to bilirubin >32 µmol/L (n=70)

• Cardiovascular organ dysfunction: cardiovascular SOFA 0 (n=618) to cardiovascular SOFA >0 (n=469)

• Central nervous system (CNS) organ dysfunction: CNS SOFA 0 (n=766) to CNS SOFA >0 (n=321)

• Renal organ dysfunction: renal SOFA 0 (n=615) to renal SOFA >0 (n=291)

The cutoffs are preferably chosen to provide a zone of uncertainty of 1 SOFA in-between to separate a more severe presence from the absence of organ dysfunction when possible. For respiratory dysfunction the zone of uncertainty between presence and absence of organ dysfunction is expanded since oxygen saturation is a volatile examination. As there is no definition for infection within the Sepsis-3 definition, we used the Linder-Mellhammar Criteria of Infection [30] (no infection = 139, verified infection = 741). The sepsis panel was determined using groups 1, 2, 3 (n=346) vs. 4, 5 (n=745) in the cohort.

For each of the 8 panels, proteins were selected using the statistical analysis and explainable machine learning methods described above. The percentage a protein is found in the top-10 in a bootstrap iteration (n=100) was calculated, and the top 90th percentile were selected as the representative proteins for a particular comparison.

To verify the biological function of the selected panels, we predicted high-confidence (> 0.9) functional interactions using STRING-DB [33] with up to 20 interactions and performed GO Biological Function overrepresentation tests to infer the biological function associated with the selected protein panels.

4.11 Creation of the probabilistic latent space

The probabilistic interpretable latent space is created by predicting probabilities for every patient in the training data for each of the 8 clinical outcome protein panels. This created a 1091 x 8 matrix that contains no missing values where each column represents the probability that a patient will develop one or several organ dysfunctions, infection, or sepsis based on time-of-admission samples. These probabilities are interpretable and can be associated with the molecular panel profiles identified using the explainable machine learning approach described above.

4.12 Identification and analysis of subphenotypes

The 8-dimensional ILS was clustered using the k-means algorithm (k = 5) into 5 different groups. The survival profiles were constructed for each subgroup over 365 days and were compared at 3 days, 30 days, and 365 days. To obtain the qualitative descriptions of the ILS for each of the clinical outcome protein panels, the average rate of panel prediction was calculated by dividing the number of positively predicted instances (> 0.5) by the total number of patients in a particular cluster.

The most important proteins and clinical parameters were identified for each subgroup using the statistical and explainable machine learning methods described above by comparing each cluster to a combination of all other clusters (one-vs-all). The 5 most important clinical parameters were visualized for each cluster against the values from all other clusters.

4.13 Creation of the digital family model and analysis of the unbiased held-out test set

Probabilities were predicted for the held-out test set (n=273) and neighborhoods consisting of 2.5% of the ILS were extracted the nearest neighbor algorithm to build digital families consisting of 27 patients per patient in the held-out test set. The predicted SOFA increase was derived for each family by calculating the median SOFA increase for day 1. The SOFA error was calculated as the difference between the predicted SOFA increase and the actual SOFA increase of the test patient. 30-day survival was calculated for 3 bins of predicted 30-day mortality for each digital family ([0.0-0.1], [0.1-0.3], [0.3-0.6] predicted mortality bins). For each predicted 30-day mortality rate, the actual 30-day mortality rate of the test set was calculated for the current bin. This was bootstrapped (n=1000) to get the estimated 30-mortality rate in the held-out test set.

The mortality map for sepsis was created by calculating the proportion of patients within the digital family that died within 30-days. Each patient in the test data was colored by this rate and reduced using the same UMAP used to visualize the subphenotypes. The mortality ratio was calculated by modelling the 7-day mortality rate for each family and dividing it by the 30-day mortality rate for the family. These were plotted against each other to visualize how much of the 30-day mortality occurred earlier in the hospital visit.