Derivation and Validation of Novel Phenotypes of Multiple Organ Dysfunction Syndrome in Critically Ill Children

Key Points Question Does data-driven phenotyping based on the trajectories of organ dysfunction in the acute phase of critical illness among children with multiple organ dysfunction syndrome uncover phenotypes with prognostic and therapeutic relevance? Findings In this 2-center cohort study of 20 827 pediatric intensive care encounters, a data-driven approach to phenotyping patients with multiple organ dysfunction syndrome using the trajectories of 6 organ dysfunctions uncovered 4 reproducible and distinct phenotypes with prognostic and potential therapeutic relevance. Meaning In this study, data-driven phenotyping based on the type, severity, and trajectory of 6 organ dysfunctions showed promising results in critically ill children with multiple organ dysfunction syndrome.


Introduction
The development of multiple organ dysfunction syndrome (MODS) is a common final pathway for death among critically ill children. 1 Even among survivors, children with more severe organ dysfunction are at higher risk of developing long-term morbidity after critical illness. 2 Approximately one-third of children admitted to a pediatric intensive care unit (PICU) have MODS on presentation or develop it during their illnesses, so moderating the burden of MODS could significantly affect the outcomes of critically ill children. 1,[3][4][5] MODS can develop after many types of injury, but sepsis, trauma, and major surgery are the most common etiologies in the PICU. 1,6,7 The most frequently described shared mechanism of MODS pathophysiology is the development of dysregulated inflammation; however, decades of drug trials targeting mediators of inflammation have failed to show effectiveness in the reduction of the MODS burden associated with sepsis and trauma. [7][8][9][10][11] One of the most likely explanations for this is that MODS is a complex, dynamic, and heterogeneous process with many different phenotypes, and no single strategy will be effective in all phenotypes of MODS.
In recent years, considerable emphasis has been placed on characterizing the phenotypes of heterogenous syndromes in critically ill patients. [12][13][14][15][16][17] These phenotypes may have prognostic and therapeutic implications and form the basis for precision medicine in the critical care setting. 18 The derivation of these different phenotypes has included data-driven approaches using gene expression data 12,17 and clinical data [14][15][16] as well as expert-based approaches using biomarkers and clinical data. 13 In general, these approaches have used a single time (eg, initial presentation) to define the different phenotypes and have not investigated the dynamic patterns of illness, despite evidence that this may be important in phenotyping. [19][20][21][22] Furthermore, a single snapshot approach does not account for the fact that patients may present to a critical care setting at different points in their illness and that organ dysfunctions tend to peak between days 1 and 3 of admission. 23 In this study we aimed to derive, validate, and characterize novel phenotypes of MODS in critically ill children using a data-driven approach based on the type, severity, and trajectory of organ dysfunctions in the acute phase of critical illness. Furthermore, we aimed to determine whether these phenotypes had prognostic and therapeutic relevance.

Study Design and Population
institutions' data warehouses using structured queries and underwent quality checks for conformity, completeness, and plausibility. Data analysis was conducted from March to October 2019. The institutional review boards at the Ann and Robert H. Lurie Children's Hospital of Chicago and The University of Chicago approved this study with a waiver of informed consent because of the retrospective nature of the analysis with minimal risk to patients. The reporting of this cohort study was performed using the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) reporting guideline. 37 Organ dysfunction severity was measured using the subscores of the pediatric Sequential Organ Failure Assessment (pSOFA) score, which measures dysfunction for the respiratory, cardiovascular, coagulation, hepatic, neurologic, and renal systems on a scale from 0 to 4 for each system. 24 Patients had the 6 pSOFA subscores calculated for each 24-hour period between PICU admission and day 3, which was considered the acute phase of critical illness. 23 Individual pSOFA subscores were carried forward for as long as 24 hours if they were not remeasured; otherwise, missing variables were assumed to be normal and the corresponding subscore was assigned a 0.
Patients with MODS were defined as those with a pSOFA subscore of at least 2 in at least 2 organ systems within the first 3 days of admission. Severity of illness on admission was determined using the Pediatric Risk of Mortality (PRISM) III score using variables from the first 24 hours. 25 Chronic comorbidities were based on the classification system developed by Feudtner et al. 26 Patients with immunocompromised status were defined as those with an oncologic disease or transplant recipients. Patients who received antibiotics and microbiological cultures during the first 3 days were considered to have a confirmed or suspected infection. Each patient encounter was treated independently, but only the first PICU admission in a given hospitalization was included in the analysis.

Primary and Secondary Outcomes
The primary outcome was in-hospital mortality. The secondary outcomes included the presence of persistent MODS on day 7 after PICU admission (which included patients who died in the first week), and vasoactive-free, ventilator-free, and hospital-free days at 28 days after PICU admission.

Data-Driven Phenotyping
We used subgraph-augmented nonnegative matrix factorization to group patients with MODS into phenotypes based on the type, severity, and trajectory of each of the 6 pSOFA subscores in the first 3 days of PICU stay (Figure 1; eAppendix in the Supplement). 21 Briefly, subgraph mining was used to extract the subgraphs representing the trajectory of individual organ dysfunctions. 27 We then split the cohort of 5054 encounters with MODS and adequate subgraph counts into a derivation set (4044 encounters [80.0%]) and a validation set (1010 encounters [20.0%]). Nonnegative matrix factorization (NMF) was implemented on the patient-subgraph count matrix in the derivation set, which resulted in a patient-phenotype distribution matrix and a phenotype-subgraph mixture coefficient matrix. 21,28 Patients in the derivation and validation set were then assigned phenotype membership based on the highest probability group in the mixture coefficient matrix. Patients assigned to each phenotype were then compared across the derivation and validation sets.
As a sensitivity analysis to determine external validity of the phenotypes, the same procedures were conducted using patients with MODS in the second hospital (PICU B) as the external validation set. Additional details can be found in the eAppendix in the Supplement. laboratory tests and summary measures included were based on those frequently used in prior pediatric organ dysfunction and severity of illness scores. 3,25,29,30

Prognostic Relevance of the MODS Phenotypes
The prognostic relevance of the MODS phenotypes was tested by comparing the association of each phenotype with the outcomes of interest after adjusting for common confounders of PICU-related morbidity and mortality, including age, immunocompromised state, severity of illness on admission using the PRISM-III score, and study site. 4,31 The phenotype with the lowest overall in-hospital mortality was used as reference. Additionally, the association of noncharacteristic organ dysfunctions and in-hospital mortality was analyzed in each of the phenotypes. Noncharacteristic organ dysfunctions were defined as those organ systems that were not used to characterize a specific phenotype and had a pSOFA subscore of at least 2 by day 3 of PICU stay.

Therapeutic Relevance of the MODS Phenotypes in Patients With Vasoactive-Dependent Shock
There are few organ dysfunction-specific treatments routinely used in patients with MODS. One exception is the use of intravenous (IV) hydrocortisone in patients with vasoactive-dependent shock.
Hydrocortisone may contribute to hemodynamic stability through pleiotropic effects; however, mixed results in large randomized clinical trials among adults with septic shock have led to variable use among practitioners. [32][33][34][35][36] Given that shock is a common feature of MODS, we tested the hypothesis that the heterogeneity of treatment effect of IV hydrocortisone in patients with vasoactive-dependent shock is in part explained by their phenotype. Vasoactive-dependent shock was defined as the use of dopamine greater than 5 μg/kg/min or any dose of epinephrine or norepinephrine infusions in patients with MODS. To balance the confounders across treatment groups, we first performed propensity score matching of patients based on the propensity to receive IV hydrocortisone, using age, immunocompromised state, PRISM-III score, study site, and phenotype membership as covariates. We then performed an interaction effect analysis for in-hospital mortality and vasoactive-free days, using phenotype membership, treatment group (ie, receiving Ն2 mg/kg/d of IV hydrocortisone or not), and their interaction as covariates to determine whether there was a differential treatment effect associated with phenotype membership. Additional details can be found in the eAppendix in the Supplement.

M F C ≈ ×
Subgraph mining extracts representative subgraphs from the trajectories of the 6 pediatric Sequential Organ Failure Assessment subscores in patients with MODS, which results in a matrix, M, of patient-subgraph counts. Nonnegative matrix factorization is used to derive a matrix, F, of hidden features (in this case the phenotypes), and a matrix of mixture coefficient C with the coefficients that compose the subgraph-based phenotypes. The factorization is done by iteratively updating F and C using the sparse nonnegative matrix factorization with sparseness in the left factor algorithm to gradually reduce the error between M and F × C. Once the nonnegative matrix factorization is completed, the final phenotypes can be characterized given the highly interpretable nature of the algorithm. Of note, the size of the matrices, number of groups, and the phenotype characterization examples here are used for illustration purposes only. CV indicates cardiovascular; resp, respiratory.

Statistical Analysis
Data were analyzed using R version 3.6.1 (R Project for Statistical Computing). Categorical variables were compared using the χ 2 test and continuous variables using the Kruskal-Wallis test. Regression analysis was used to adjust for confounders when comparing outcomes. Logistic regression was used to compare binary outcomes, and Poisson regression was used to compare count outcomes. Survival analysis to 28 days was performed using Kaplan-Meier curves and adjusted Cox regression analysis.
Statistical significance was set at P < .05, and all tests were 2-tailed.

Study Population
There were 20 827 patient encounters among 14

Data-Driven Phenotyping
After

Characterization of the MODS Phenotypes
The most representative subgraphs for each phenotype were chosen from the mixture coefficient matrix of the NMF model, such that the selected subgraphs accounted for at least 75% of the phenotype group coefficients. The 10 most representative subgraphs for each of the 4 MODS phenotypes and the relative weight of each organ dysfunction subscore in the first 3 days, stratified by phenotype, appear in Figure 2 and eFigure 1 in the Supplement. These representative subgraphs were used to characterize the phenotypes based on the most distinctive characteristics, using the organ dysfunction type, severity, and trajectory in the first 3 days of PICU stay. The clinical characteristics and associated outcomes of patients in each of the 4 MODS phenotypes are presented in Table 2. When compared across age groups, the distribution and outcomes associated with the 4 MODS phenotypes remained consistent (eTable 5 in the Supplement).

Prognostic Relevance of the MODS Phenotypes
The outcomes associated with each of the 4 phenotypes after adjusting for confounders including age, immunocompromised state, PRISM-III score, and study site are presented in eTable 6 in the Supplement. After adjusting for the same confounders in Cox regression analysis and using phenotype 2, which had the lowest mortality as reference, the adjusted hazard ratios (aHRs) for

Therapeutic Relevance of the MODS Phenotypes in Patients With
Vasoactive-Dependent Shock Interaction analysis in the matched cohort revealed that IV hydrocortisone had a differential treatment effect on vasoactive-free days across phenotypes, most notably with phenotype 3 having most benefit associated with treatment (23 vasoactive-free days in treated patients vs 18 vasoactive-free days in treated patients; P for interaction < .001). The interaction was not significant for in-hospital mortality ( Table 3).

Discussion
We derived and validated 4 novel, data-driven phenotypes of MODS in critically ill children based on organ dysfunction trajectories. These phenotypes have distinct clinical characteristics, are independently associated with outcomes, and have different sets of organ dysfunction-based risk factors for death. In a subset of patients with vasoactive-dependent shock who were matched based on the propensity to receive IV hydrocortisone, treatment was associated with a difference in  The most common subgraphs in each phenotype and the proportion of patients who had each of them.  vasoactive-free days across phenotypes, suggesting that these phenotypes may be associated with differences in response to therapy.

JAMA Network Open | Critical Care Medicine
Our findings share similarities with previously described data-driven phenotypes of sepsis and sepsis-induced MODS, particularly phenotypes 3 and 4, which were the MODS phenotypes most frequently associated with infections in our cohort. Phenotype 3 shares similarities with the "shock with hypoxemia and altered mental status" phenotype described by Knox et al 15  particularly in terms of the thrombocytopenia, hepatobiliary dysfunction, and shock as well as the relatively preserved pulmonary function. To our knowledge, our study is the first to use a data-driven approach to derive MODS phenotypes in critically ill children. Further biomarker and molecular endotyping analyses will be needed to compare our phenotypes with other previously described clinical phenotypes and endotypes in pediatric and adult patients. 12,13,17 The association of IV hydrocortisone administration with differential treatment effect on vasoactive-free days in patients with vasoactive-dependent shock, particularly in phenotype 3, could partially be explained by the enrichment with patients suffering severe hypoxemia, severe shock, and infections. Prior randomized clinical trials and secondary analyses of trial data have found that the use of hydrocortisone was associated with improved outcomes in patients with sepsis-associated acute respiratory distress syndrome. 38,39 However, our results should be interpreted with caution because they are part of an interaction effect analysis on a propensity score matched cohort. The standard approach is to perform an interaction effects analysis across subgroups with a treatment assignment that has been randomized to ensure a better balance of all possible confounders. 16 Given the selection bias inherent in observational data, we used propensity score matching to balance  Similarly, renal dysfunction seems to most negatively affect patients in phenotype 3 who suffer severe hypoxemia and shock, indicating that careful fluid management and nephrotoxin avoidance might be most beneficial in this group of patients. In addition, further endotyping and molecular characterization of these broad phenotypes may uncover relevant MODS subphenotypes that may benefit from novel targeted interventions, including immunomodulatory strategies. [41][42][43] Furthermore, the molecular characterization of these phenotypes may help in the task of assigning phenotype membership shortly after patients are admitted, which would make this phenotyping schema more clinically actionable in the critical care setting. In short, data-driven phenotyping may help develop precision medicine strategies that can reduce the mortality and morbidity associated with MODS, but further investigation into the value of this phenotyping approach in research and clinical care is warranted. 18

Limitations
Our findings are subject to several limitations. First, the data for this study were collected retrospectively from 2 academic PICUs in the same city, limiting the generalizability of our findings.
Second, our data-driven phenotyping approach was intentionally sparse, using only organ dysfunction trajectories to define broad phenotypes and determine their prognostic and therapeutic relevance, but performing deeper phenotyping using other clinical features and biomarkers may result in phenotypes with potentially higher resolution and greater clinical utility. Third, the therapeutic relevance of IV hydrocortisone was studied using an interaction effect analysis on a propensity score matched cohort, and the results must be interpreted with caution, as previously discussed. Fourth, although we described the MODS phenotypes as distinct entities, it is likely that many patients share features of more than 1 phenotype, and future analyses to better understand the implications of these overlaps are warranted.

Conclusions
In this study, we derived and validated 4 data-driven phenotypes of MODS in critically ill children that have distinct clinical characteristics, are independently associated with outcomes, and may be therapeutically relevant in the subset of patients with vasoactive-dependent shock. Our data-driven phenotyping approach in MODS shows promising results, but further characterization of these broad phenotypes and validation in different settings is warranted.