Unsupervised Clustering Analysis Based on MODS Severity Identifies Four Distinct Organ Dysfunction Patterns in Severely Injured Blunt Trauma Patients

Purpose: We sought to identify a MODS score parameter that highly correlates with adverse outcomes and then use this parameter to test the hypothesis that multiple severity-based MODS clusters could be identified after blunt trauma. Methods: MOD score across days (D) 2–5 was subjected to Fuzzy C-means Clustering Analysis (FCM) followed by eight Clustering Validity Indices (CVI) to derive organ dysfunction patterns among 376 blunt trauma patients admitted to the intensive care unit (ICU) who survived to discharge. Thirty-one inflammation biomarkers were assayed (Luminex™) in serial blood samples (3 samples within the first 24 h and then daily up to D 5) and were analyzed using Two-Way ANOVA and Dynamic Network analysis (DyNA). Results: The FCM followed by CVI suggested four distinct clusters based on MOD score magnitude between D2 and D5. Distinct patterns of organ dysfunction emerged in each of the four clusters and exhibited statistically significant differences with regards to in-hospital outcomes. Interleukin (IL)-6, MCP-1, IL-10, IL-8, IP-10, sST2, and MIG were elevated differentially over time across the four clusters. DyNA identified remarkable differences in inflammatory network interconnectivity. Conclusion: These results suggest the existence of four distinct organ failure patterns based on MOD score magnitude in blunt trauma patients admitted to the ICU who survive to discharge.


INTRODUCTION
Trauma remains the leading cause of mortality and morbidity for individuals under 55 years and accounts for 30% of all life-years lost, with over 190,000 lives lost annually in the USA (1,2). With advances in prehospital transport and resuscitation strategies, patterns of traumatic death have significantly changed over the past three decades (3)(4)(5). Trauma-related deaths have now assumed a largely bimodal distribution, with a vast majority of deaths occurring at the scene or within the first day after injury as a consequence of massive head injury or uncontrolled bleeding (5,6). However, patients who survive beyond the initial traumatic insult are prone to develop a state of persistent critical illness manifested by prolonged intensive care unit (ICU) and hospital length of stays (LOS), and a persistent risk of late in-hospital and post-discharge complications (7)(8)(9)(10).
A common central factor contributing to outcomes following injury is the accompanying immuno-inflammatory response. If this response is appropriate in magnitude and duration it can aid in re-establishing host homeostasis. However, a dysregulated response is associated with multiple organ dysfunction syndrome (MODS) which can evolve to a state of persistent critical illness and a continued increased risk for complications and death after discharge (11)(12)(13). Trauma-induced MODS is widely believed to be the leading cause of death among ICU patients being responsible for 50-80% of ICU mortality (14,15). Typically, MODS peaks within 5 days of injury and is associated with a complicated clinical course; however, the type and number of distinct organ failure patterns that occur after injury are not known. Hence, the current challenge in the early management of severely injured blunt trauma patients is to predict and then prevent MODS (16)(17)(18). However, to do this, there is a need to define the patterns of organ dysfunction in trauma patients that survive the early mortality window and whether distinct MODS patterns are associated with identifiable differences in the early systemic inflammatory response.
In the current study, we utilized an unsupervised clustering strategy to identify the number of MODS-based phenotypes that followed severe blunt injury in a cohort of 493 trauma patients admitted to the ICU. To optimize our study cohort, patients discharged prior to day 5 from the ICU and 21 patients who died before discharge were excluded, which yielded a total of 376 patients used in the current study. The MOD score magnitude over days 2-5 was found to correlate well with adverse in-hospital outcomes and was used to identify four distinct severity-based MODS clusters. The four clusters exhibited differential early inflammation biomarker profiles and correlated with subsequent in-hospital adverse outcomes. These findings provide evidence for the emergence of multiple definable organ dysfunction patterns after severe blunt injury. This information can be useful for the identification of prognostic variables to predict organ dysfunction severity and patterns following blunt traumatic injury.

Patient Enrollment, Sampling, and Data Collection
Blunt trauma patients deemed eligible for enrollment were at least 18 years of age at time of the trauma, admitted to the ICU as part of the post-trauma management, and were expected Abbreviations: MODS, multiple organ dysfunction syndrome; ISS, injury severity score; ICU, intensive care unit; IL, interleukin; DyNA, dynamic network analysis. to survive beyond the initial 24 h post-injury as per the oncall trauma surgeon. Reasons for ineligibility were isolated head injury or brain death criteria, or pregnancy. Three plasma samples were collected within the first 24 h following injury, as follows: (1) the initial blood draw upon arrival within 4 h from time of injury; (2) within 4-12 h of admission to the emergency department (ED); and (3) at 24 h of ED admission. Subsequent samples were obtained from day (D) 1 to D5 postinjury. Demographic and clinical data were collected from the inpatient electronic medical record and the trauma registry database. The clinical database and biobank were maintained prospectively from 2004 to 2012, for a total of 8 years period. Totally 493 patients were enrolled in the observational study. The Marshall Multiple Organ Dysfunction (MOD) score (19) was calculated daily during the patients' ICU stay.

Optimization of Study Cohort
Organ failure is known to peak within the first 5 days after injury (20). Therefore, out of the 493 patients enrolled in the observational study (who were admitted to the ICU of the UPMC Presbyterian University Hospital, a Level 1 trauma center), we identified a subset of patients (n = 376) with complete sequential MOD score data who remained in the ICU for at least 5 days postinjury. Patients were excluded because they were discharged from the ICU prior to 5 days and therefore had incomplete MOD score data (n = 96) or because they died prior to hospital discharge (n = 21) and therefore the incidence of in hospital complication rates could not be assessed. The patients that died in-hospital have been previously described in detail (21) and are referenced in this study as a separate group.

Fuzzy C-Means Clustering
To identify the number of distinct MOD score severity-based clusters present in the first 5 days after injury, the sequential MOD scores across days 2-5 for the 376 patients admitted to the ICU after trauma were subjected to fuzzy C-means (FCM) clustering. The FCM is a soft partition, unsupervised clustering method that allows each piece of data to belong to more than one cluster (22). The FCM assigns membership values to each of the data points that indicate the degree to which the data points belong to the different clusters. This feature, in some degree, fits the characteristic of heterogeneous clinical data that exhibits no clear boundaries between clusters. The objective function of the FCM algorithm is to minimize the objective function (see below): . , x n } is a collection of n elements. In our case, n = 376, and x i represents a vector of consisting of ith patient's four MOD scores from days 2-5. µ ij ∈ [0, 1] is the membership of x i to jth cluster, m represents the fuzzifier parameter (set up as 2), c is the number of clusters with cluster centers V = {v 1 , . . . , v c }, and x i − v j represents the distance of x i to the center of jth cluster.
As described in the objective function, the number of clusters (c) needs to be preset. To define the optimal number of clusters, we performed FCM based on the Euclidean distance with c set from 2 (lowest) to 6 (highest) clusters. Next, to determine the optimal number of clusters, we utilized eight internal Clustering Validity Indices (CVI), which evaluate the goodness of a clustering structure relying only on the information in the data and allow for the quantification of intra-cluster compactness or inter-cluster separation (23,24). The CVI analysis included: C index, Dunn index, Gamma index, GDI index, G_plus index, PBM index, S_DBW index, and Tau index (further details of the clustering validity indices are available at: https://cran.r-project.org/web/packages/clusterCrit/vignettes/ clusterCrit.pdf). For each index, specific rule of either computing the greatest index value (max) or the smallest index value (min) must be applied in order to determine the best index value for optimal partition. Among the CVIs aforementioned, the best partition is the one corresponding to the "min" value of C index, G-plus index, S-DBW index and the "max" value of indices including Dunn index, Gamma index, GDI index, PBM index, and Tau index. These analyses were carried out using R (The R Project for Statistical Computing, Version 3.2.2).

Statistical Analysis
All data are expressed as mean ± SEM. Statistical analysis between groups was performed by One-way Analysis of Variance (ANOVA) followed by Tukey post-hoc analysis using SigmaPlot TM 11 software (Systat Software, Inc., San Jose, CA). Fisher's exact test was performed for categorical data using Graphpad PRISM (GraphPad Software, Inc., La Jolla, CA). Group-time interaction of plasma inflammatory mediators' levels was determined by Two-Way ANOVA. P < 0.05 was considered statistically significantly different for all analyses. Dynamic Network Analysis (DyNA) was performed to gain insights into the temporal dynamic changes in network connectivity of the post-traumatic inflammatory response [as we have shown previously (8, 21,25)] among the FCM-defined clusters.

Clustering Validity Indices Identify Four Distinct MOD Score Clusters Following Severe Blunt Trauma
We focused on days 2-5 because this incorporates the known peak in MODS post-injury (20) but avoids the impact of inadequate resuscitation on MOD score sometimes observed in the first 24 h. To determine the number of distinct severitybased MODS clusters that appear after blunt injury, the MOD score data across D2-D5 were subjected to Fuzzy Clustering Analysis followed by eight separate clustering validity indices. Six indices indicated that four clusters was the optimal number, while two indices suggested three clusters (Figure 1). Based on these analyses, FCM segregated the patient cohort into the following four clusters: Cluster 1 (n = 199, mean MOD score = 0.28 ± 0.02), Cluster 2 (n = 99, mean MOD score = 1.97 ± 0.07), Cluster 3 (n = 53, mean MOD score = 3.99 ± 0.12), and Cluster 4 (n = 25, mean MOD score = 7.13 ± 0.23). The MOD score values on individual days were statistically different between the four clusters as determined by Two-way ANOVA (Figure 2).

The Four MOD Score Clusters Differed in Injury Patterns and Presentation Characteristics
In terms of overall demographics, there was no statistically significant difference in average age or gender distribution among the four clusters (Table 1). However, Cluster 1 had a statistically significantly lower average injury severity score (ISS) than Clusters 2-4 (P 2vs.1 = 0.041; P 3vs.1 < 0.001). There was no statistical difference in ISS between Clusters 2, 3, and 4. To determine if injury patterns differed between the clusters, the abbreviated injury scale for six body regions were compared ( Table 1). Cluster 2 exhibited greater rates of abdominal and extremity injuries (P = 0.014 and P = 0.003, respectively) when compared to Cluster 1. Patients in Cluster 3 had higher head and neck injury scores when compared to Cluster 1 and 2 (P = 0.011 and 0.02, respectively) and a statistically significantly higher incidence of brain injury than Cluster 1 (P = 0.003) and Cluster 2 (P = 0.010). There was no difference in brain injury rates between Clusters 3 and 4 ( Table 1). Thus, Cluster 1 patients were less severely injured while Clusters 3 and 4 included patients that were more likely to have traumatic brain injury.
Next, we identified the differences in physiologic and biochemical data on presentation among the four MOD score clusters. Clusters 2-4 had lower average systolic blood pressures and hemoglobin levels upon presentation when compared to patients in Cluster 1 ( Table 2). Patients in Cluster 3 had higher blood creatinine levels on admission compared to Clusters 1 and 2. Admission coagulation parameters (Prothrombin time, International Normalized Ratio, and Partial Thromboplastin  Table 2). Therefore, patients in Clusters 2-4 were more likely to be in shock at admission while patients in Clusters 3 and 4 were more likely to present with evidence of renal dysfunction and coagulation abnormalities.

MOD Score Clusters Differ in Clinical Outcomes
There were statistically significant differences among the four clusters with regards to in-hospital outcomes, including ICU and total hospital length of stay (LOS), days on mechanical ventilation as well as the incidence of NI being all greatest in  Table 3). Surgical intervention rates (within the first 24 h-all types of procedures) were lowest in Cluster 1 and were significantly different between Cluster 1 and Clusters 2-4 ( Table 3). Patients in Cluster 4 were more likely to require a vascular intervention. Patients in Cluster 2-4 were more likely to receive a transfusion in the first 24 h than patients in Cluster 1 and patients in Cluster 4 received significantly greater volumes of packed red blood cells (PRBC) and fresh frozen plasma (FFP) when compared to patients in Clusters 1-3 ( Table 3). These findings further establish that patients that fall into Cluster 1 have less severe injuries than patients in the other clusters, while patients in Cluster 4 are distinguished by a greater need for both PRBC and FFP.

Disparate Contribution of Individual Organ Failure Components Among the Four Clusters
The four clusters differed not only in the average magnitude of MOD scores, but also in organ failure patterns (Supplemental Figure 1). Clusters 2-4 exhibited respiratory, cardiovascular, hematologic, and neurologic dysfunction scores that increased significantly between each severity-based cluster. Clusters 3 and 4 were distinguished from Clusters 1 and 2 by worse renal function. A notable increase respiratory and cardiovascular dysfunction scores was observed in Cluster 4 (Supplemental Figure 1). Therefore, organ dysfunction becomes progressively worse across the clusters for all systems except for the hepatic component, with a notable increase in average renal dysfunction scores in Cluster 3 and marked worsening of the severity of respiratory and cardiovascular dysfunction in Cluster 4.

Distinct Inflammatory Patterns Emerge Among the Four Clusters
Seven out of the 31 biomarkers assayed exhibited statistically significant differences among the clusters upon admission and over time (Supplemental Figure 2).  . Notably, Cluster 1 exhibited a highly connected network within the first 16 h post-injury that dissipated rapidly thereafter. The networks in Cluster 2 consisted of sparsely connected networks, a pattern that persisted to day 5. In clear distinction, Cluster 3 had an increase in network connectivity over time with the greatest connectivity among mediators observed at D2-D5. Finally, Cluster 4 also exhibited a unique pattern with highly connected but uncoordinated networks throughout the 5 days. This analysis suggests that the inflammation profiles diverge early and in conjunction with the evolution of the severity and patterns of MODS following severe blunt trauma.

Characteristics of Excluded Patient Cohorts
In order to assure availability of complete MOD score data from D2 to D5 and NI rates through discharge, patients with incomplete MOD D2-D5 data (n = 96) or that died in-hospital (n = 21) were not included in the initial clustering analysis. Among the 96 excluded patients, 91 of them were discharged from ICU prior to day 4, and their characteristics were comparable to patients in Cluster 1 (Supplemental Table 1). We have described previously the characteristics of the non-survivor cohort (21). To provide a comparison of the average MOD score values over time among the four clusters identified for survivors and the nonsurvivor cohort, we inserted the MOD score averages D2-D5 from the published non-survivor cohort with the curves for the survivor cohort (Supplemental Figure 3). It is interesting to note that MOD score values start on average lower in non-surviving patients than surviving patients in Clusters 3 and 4 but then rise steadily.

DISCUSSION
In the current study, we set out to identify organ dysfunction phenotypes in severely injured blunt trauma patients that survive to discharge using an unsupervised clustering strategy. The FCM followed by CVIs defined four distinct MODS patterns. These four clusters were not only different in the magnitude of MODS, but also in organ failure patterns exhibited by unique patterns of six Marshall MODS components. The clusters also differed both in the clinical features and inflammatory profiles upon admission and over time up to day 5 post-injury. This analysis suggests that it may be feasible to stratify critically ill trauma patients early in the clinical course into sub-groups at risk for multiple clinical trajectories defined by specific patterns and magnitude of organ dysfunction, which in turn may be useful in supporting tailored research and clinical therapies for blunt trauma patients. Improvements in medicines have led to improved prognosis after severe injury (3,26). Data from recent studies on large patient cohorts shows that MODS peaks within the first 5 days after injury (16,20). Shepherd et al. (20) describe three broad patterns of organ function scores after injury including patients without MODS, and patients with either resolving MODS or non-resolving MODS. Others have defined a state of persistent critical illness characterized by prolonged ICU LOS associated with isolated organ dysfunction and risk for complications, including nosocomial infection (11,13). These patients probably coincide with the non-resolving MODS in the work from Shepherd et al. (20).
In our study, we utilized a cohort of trauma patients deemed injured severely enough to be admitted to the ICU. We narrowed the cohort further to analyze data only from patients that remained in the ICU for at least 5 days to capture the 5 day "peak" for MODS and to assure that we included patients at risk for some degree of organ dysfunction. Others have shown adding MOD score information beyond day 5 may not add be useful to predicting MODS based phenotypes (27). However, this would require further analysis.
Using an unbiased clustering strategy, we determined that organ failure magnitude values between days 2 and 5 after injury fit into one of four distinct clusters. One of prime determinants of cluster destination is likely to be the injury and early shock characteristics. Cluster 1 was comprised of patients with moderate injury (average ISS = 18.7 ± 0.7) and minimal shock on presentation. ISS significantly increased in Cluster 2, indicating that Cluster 2 patients were more likely to have moderate-severe injury and, based on shock index, have some degree of shock; thus, these are defined as patients more likely to present with moderate-severe injury + mild shock/blood    loss. Although ISS values were not different among Clusters 2-4, Cluster 3 patients had a higher likelihood of head and traumatic brain injury (TBI), and therefore were characterized by moderate-severe injury + mild shock + TBI. Cluster 4 included the patients who progressed to the highest MOD scores and had the same head injury rates as Cluster 3. These patients, despite having the same head injury rates as Cluster 3, had significantly higher transfusion requirements; therefore, this cluster appears to comprise patients presenting with moderate -severe injury + TBI + severe shock/blood loss. While some overlap in the presenting characteristics can be seen across the clusters, it is highly likely that the differences in presenting injury and early physiologic characteristics drive most of the subsequent organ dysfunction dynamics. Based in the prolonged ICU stays and complication rates observed in Clusters 3 and 4, it appears that these patient group would fall into the "non-resolving MODS" or "persistent critical illness" categories described by others (11,13). Clusters 1 and 2 would likely fit into the phenotype of resolving MODS, with many of the patients in Cluster 1 exhibiting little evidence of MODS. However, our studies raise the possibility that even more precise designations than those described previously might be feasible based on MODS score magnitude and duration, or organ-specific responses. Clinical manifestations of organ failure following trauma are diverse. Previous work examining the patterns of organ failure among individual trauma patients suggested a cumulative temporal sequence of single organ failures (28)(29)(30). Since the lungs are highly sensitive to mediator-induced inflammation, they are often the first system to show signs of failure (31)(32)(33)(34). However, acute respiratory failure can precede MODS and may represent a trigger factor of subsequent dysregulated immune events that may contribute to remote organ failures (28,29,34). Although the typical sequence of organ failure in trauma patients is difficult to predict, the cardiovascular system is the second system to fail followed by renal failure (28,35). It is apparent that there is a progressive increase in organ dysfunction with the increase in the MOD scores from cluster to cluster. While this is not surprising, our results show some interesting patterns, including a marked increase in hematologic and neurologic dysfunction between Clusters 1 and 2; a significant increase in the appearance of renal dysfunction in Clusters 3 and 4; and a large increase in the severity of cardiovascular and respiratory dysfunction in Cluster 4. Thus, there appears to be a hierarchy of organ dysfunction that may be driven in part by injury severity or pattern and the concomitant presence of shock and/or need for transfusion. Notably, age does not appear to a major distinguishing characteristic across the four clusters defined in the current study, despite multiple indications that age is a major complicating factor in trauma outcomes (6,36,37). This may be due in part, to the range of injury severity represented within our patient cohort. For example, we have recently shown that young and old patients that had moderate severity injury experienced a similar level of MODS (38). Two age-based sub-cohorts were identified, namely young (age: 18-30 years) and aged (age: 65-90 years) matched for ISS in the moderate range. Analysis of the average MOD score between these two groups showed that there were no statistically significantly differences with regards to the average MOD score between the two groups.
In other work, we have shown that levels of specific inflammatory mediators and DyNA patterns correlate with injury severity (39), hypotension (40), death (21), and complications such as nosocomial infection (8). Certain single nucleotide polymorphisms may also segregate patients at risk for organ dysfunction after trauma (41). Others have provided evidence that gene expression patterns in circulating leukocytes could be used to identify patients as risk for MOF after trauma (9,42). Whether more precise stratification tools can be developed to prognosticate for the MODS clusters identified here is unknown. Our analysis of inflammation biomarker levels show clear differences the levels of certain inflammation biomarkers early and over time among the clusters. In addition, DyNA provides a visual window into the differences-and inferred degrees of coordination or lack thereof-in the dynamic changes in biomarker levels between the clusters over the early time periods following injury. Although significant differences in biomarker levels could be found among the clusters, we expect that levels of biomarkers alone may be insufficient to prognosticate for multiple clusters, and that other parameters as well as assessment of changes in biomarkers over time may be needed. For example, we have demonstrated the utility of patient-specific Principal Component Analysis based on inflammatory mediators assayed in the first 24 h following admission in differentiating later courses of MOD (43). In general, our observations point to stepwise increases in the magnitude and duration of the systemic inflammatory response with the MODS severity clusters.
As part of our inclusion criteria, we excluded patients died within first 24 h, and as expected the in-hospital mortality in the study cohort was low (21 out of 493 patients or 4.4%). We have previously reported the characteristics of the nonsurvivors (21), and found that the MODS trajectory was unique in these patients compared to those seen in the four clusters identified in the survivors. Although the average MOD score of the non-survivors would have segregated these patients into Cluster 2 or 3, the pattern was unique in that the MOD score started low and then steeply ascended. Whether patient destined to die during the initial hospitalization would constitute an additional MODS cluster will require additional work and larger patient numbers.
We have recently proposed a time window-based platform as a hypothetical construct for outcome stratification of trauma patients (44). This platform takes advantage of the fact that the time of the traumatic event is typically known, and facilitates adopting Precision Medicine methodology to quantify individualized injury response indices and thereby better prognosticate for adverse outcomes. We propose that optimization of this quantitative platform using MODS and subsets of MODS, i.e., cluster-based MODS, as a composite endpoint will lead to more informed early decision making, guide early interventions, improve quantifiable short-and longterm outcome indices, and could potentially facilitate tailored treatments or directed research.
There are several limitations to note in our study. First, this is a single institution study and 376 blunt trauma patients may not be sufficient to identify all of the organ failure phenotypes. Second, the study excluded blunt trauma non-survivors, and only patients with MOD score between D2 and D5 who survived to discharge were included. However, comparing the MOD score trajectory of the non-survivors to the four identified clusters in the survivor cohort suggests that non-survivors exhibit a distinct MODS pattern that can be differentiated from the four clusters by day 5. Third, our biomarker panel included only 31 inflammatory mediators that represent only a fraction of all potential circulating biomarkers. Finally, confirmatory studies involving multiple institutions that include contemporary patient data sets that incorporate penetrating trauma patients will be needed to confirm the number of MODS clusters.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by University of Pittsburgh Institutional Review Board. The patients/participants or next of kin provided their written informed consent to participate in this study, as per IRB regulations.

AUTHOR CONTRIBUTIONS
DL: literature search, data analysis, and writing. RN: data collection, data analysis, and writing. YV: study design and critical revision. AP and RS: study design. HY: critical revision. QM: data analysis and data interpretation. TB: study design, data interpretation, and critical revision.  Supplemental Figure 3 | Average MOD scores of the four FCM-derived Clusters compared to the average MOD score of the non-survivor patients (n = 21) from day 2 through day 5 post-injury. Data presented as mean with 95% confidence interval (CI).
Supplemental Table 1 | Comparison of demographic and injury pattern characteristics, admission physiological and biochemical parameters between Cluster 1 patients (n = 199) and the excluded survivor patients (n = 96). Values are expressed as mean ± SEM. Mann-Whitney U-test, One-Way ANOVA, and Fisher exact test were used as appropriate with statistical significance set at P < 0.05.