IMPORTANCE: Multiple organ dysfunction syndrome (MODS) is a dynamic and heterogeneous process associated with high morbidity and mortality in critically ill children. OBJECTIVE: To determine whether data-driven phenotypes of MODS based on the trajectories of 6 organ dysfunctions have prognostic and therapeutic relevance in critically ill children. DESIGN, SETTING, AND PARTICIPANTS: This cohort study included 20 827 pediatric intensive care encounters among 14 285 children admitted to 2 large academic pediatric intensive care units (PICUs) between January 2010 and August 2016. Patients were excluded if they were older than 21 years or had undergone cardiac surgery. The 6 subscores of the pediatric Sequential Organ Failure Assessment (pSOFA) score were calculated for the first 3 days, including the subscores for respiratory, cardiovascular, coagulation, hepatic, neurologic, and renal dysfunctions. MODS was defined as a pSOFA subscore of at least 2 in at least 2 organs. Encounters were split in a 80:20 ratio for derivation and validation, respectively. The trajectories of the 6 subscores were used to derive a set of data-driven phenotypes of MODS using subgraph-augmented nonnegative matrix factorization in the derivation set. Data analysis was conducted from March to October 2019. EXPOSURES: The primary exposure was phenotype membership. In the subset of patients with vasoactive-dependent shock, the interaction between hydrocortisone and phenotype membership and its association with outcomes were examined in a matched cohort. MAIN OUTCOMES AND MEASURES: The primary outcome was in-hospital mortality. Secondary outcomes included persistent MODS on day 7, and vasoactive-free, ventilator-free, and hospital-free days. Regression analysis was used to adjust for age, severity of illness, immunocompromised status, and study site. RESULTS: There were 14 285 patients with 20 827 encounters (median [interquartile range] age 5.2 years [1.5-12.7] years; 11 409 [54.8%; 95% CI, 54.1%-55.5%] male patients). Of these, 5297 encounters (25.4%; 95% CI, 24.8%-26.0%) were with patients who had MODS, of which 5054 (95.4%) met the subgraph count threshold and were included in the analysis. Subgraph augmented nonnegative matrix factorization uncovered 4 data-driven phenotypes of MODS, characterized by a combination of neurologic, respiratory, coagulation, and cardiovascular dysfunction, as follows: phenotype 1, severe, persistent encephalopathy (1019 patients [19.2%]); phenotype 2, moderate, resolving hypoxemia (1828 patients [34.5%]); phenotype 3, severe, persistent hypoxemia and shock (1012 patients [19.1%]); and phenotype 4, moderate, persistent thrombocytopenia and shock (1195 patients [22.6%]). These phenotypes were reproducible in a validation set of encounters, had distinct clinical characteristics, and were independently associated with outcomes. For example, using phenotype 2 as reference, the adjusted hazard ratios (aHRs) for death by 28 days were as follows: phenotype 1, aHR of 3.0 (IQR, 2.1-4.3); phenotype 3, aHR of 2.8 (IQR, 2.0-4.1); and phenotype 4, aHR of 1.8 (IQR, 1.2-2.6). Interaction analysis in a matched cohort of patients with vasoactive-dependent shock revealed that hydrocortisone had differential treatment association with vasoactive-free days across phenotypes. For example, patients in phenotype 3 who received hydrocortisone had more vasoactive-free days than those who did not (23 days vs 18 days; P for interaction < .001), whereas patients in other phenotypes who received hydrocortisone either had no difference or had less vasoactive-free days. CONCLUSIONS AND RELEVANCE: In this study, data-driven phenotyping in critically ill children with MODS uncovered 4 distinct and reproducible phenotypes with prognostic relevance and possible therapeutic relevance. Further validation and characterization of these phenotypes is warranted.