Uncovering the heterogeneity and temporal complexity of neurodegenerative diseases with Subtype and Stage Inference — TEI-tagged view

Alexandra L YoungpersName, Razvan V MarinescupersName, Neil P OxtobypersName, Martina BocchettapersName, Keir YongpersName, Nicholas C FirthpersName, David M CashpersName, David L ThomaspersName, Katrina M DickpersName, Jorge CardosopersName, John van SwietenpersName, Barbara BorronipersName, Daniela GalimbertipersName, Mario MasellispersName, Maria Carmela TartagliapersName, James B RowepersName, Caroline GraffpersName, Fabrizio TagliavinipersName, Giovanni B FrisonipersName, Robert LaforcepersName, Elizabeth FingerpersName, Alexandre de MendonçapersName, Sandro SorbipersName, Jason D WarrenpersName, Sebastian CrutchpersName, Nick C FoxpersName, Sebastien OurselinpersName, Jonathan M SchottpersName, Jonathan D RohrerpersName, Daniel C AlexanderpersName, Christin AnderssonpersName, Silvana ArchettipersName, Andrea ArighipersName, Luisa BenussipersName, Giuliano BinettipersName, Sandra BlackpersName, Maura CosseddupersName, Marie FallströmpersName, Carlos FerreirapersName, Chiara FenogliopersName, Morris FreedmanpersName, Giorgio G FumagallipersName, Stefano GazzinapersName, Roberta GhidonipersName, Marina GrisolipersName, Vesna JelicpersName, Lize JiskootpersName, Ron KerenpersName, Gemma LombardipersName, Carolina MarutapersName, Lieke MeeterpersName, Simon MeadpersName, Rick van MinkelenpersName, Benedetta NacmiaspersName, Linn ÖijerstedtpersName, Alessandro PadovanipersName, Jessica PanmanpersName, Michela PievanipersName, Cristina PolitopersName, Enrico PremipersName, Sara PrionipersName, Rosa RademakerspersName, Veronica RedaellipersName, Ekaterina RogaevapersName, Giacomina RossipersName, Martin RossorpersName, Elio ScarpinipersName, David Tang-WaipersName, Hakan ThonbergpersName, Pietro TiraboschipersName, Ana VerdelhopersName, Michael W WeinerpersName, Paul AisenpersName, Ronald PetersenpersName, Clifford R JackpersName, William JagustpersName, John Q TrojanowkipersName, Arthur W TogapersName, Laurel BeckettpersName, Robert C GreenpersName, Andrew J SaykinpersName, John MorrispersName, Leslie M ShawpersName, Zaven KhachaturianpersName, Greg SorensenpersName, Lew KullerpersName, Marc RaichlepersName, Steven PaulpersName, Peter DaviespersName, Howard FillitpersName, Franz HeftipersName, Davie HoltzmanpersName, M Marcel MesulampersName, William PotterpersName, Peter SnyderpersName, Adam SchwartzpersName, Tom MontinepersName, Ronald G ThomaspersName, Michael DonohuepersName, Sarah WalterpersName, Devon GessertpersName, Tamie SatherpersName, Gus JiminezpersName, Danielle HarveypersName, Matthew BernsteinpersName, Paul ThompsonpersName, Norbert SchuffpersName, Bret BorowskipersName, Jeff GunterpersName, Matt SenjempersName, Prashanthi VemuripersName, David JonespersName, Kejal KantarcipersName, Chad WardpersName, Robert A KoeppepersName, Norm FosterpersName, Eric M ReimanpersName, Kewei ChenpersName, Chet MathispersName, Susan LandaupersName, Nigel J CairnspersName, Erin HouseholderpersName, Lisa Taylor-ReinwaldpersName, Virginia LeepersName, Magdalena KoreckapersName, Michal FigurskipersName, Karen CrawfordpersName, Scott NeupersName, Tatiana M ForoudpersName, Steven PotkinpersName, Li ShenpersName, Kelley FaberpersName, Sungeun KimpersName, Kwangsik NhopersName, Leon ThalpersName, Neil BuckholtzpersName, Marylyn AlbertpersName, Richard FrankpersName, John HsiaopersName, Jeffrey KayepersName, Joseph QuinnpersName, Betty LindpersName, Raina CarterpersName, Sara DolenpersName, Lon S SchneiderpersName, Sonia PawluczykpersName, Mauricio BeccerapersName, Liberty TeodoropersName, Bryan M SpannpersName, James BrewerpersName, Helen VanderswagpersName, Adam FleisherpersName, Judith L HeidebrinkpersName, Joanne L LordpersName, Sara S MasonpersName, Colleen S AlberspersName, David KnopmanpersName, Kris JohnsonpersName, Rachelle S DoodypersName, Javier Villanueva-MeyerpersName, Munir ChowdhurypersName, Susan RountreepersName, Mimi DangpersName, Yaakov SternpersName, Lawrence S HonigpersName, Karen L BellpersName, Beau AncespersName, Maria CarrollpersName, Sue LeonpersName, Mark A MintunpersName, Stacy SchneiderpersName, Angela OliverpersName, Daniel MarsonpersName, Randall GriffithpersName, David ClarkpersName, David GeldmacherpersName, John BrockingtonpersName, Erik RobersonpersName, Hillel GrossmanpersName, Effie MitsispersName, Leyla de Toledo-MorrellpersName, Raj C ShahpersName, Ranjan DuarapersName, Daniel VaronpersName, Maria T GreigpersName, Peggy RobertspersName, Marilyn AlbertpersName, Chiadi OnyikepersName, Daniel D’AgostinopersName, Stephanie KielbpersName, James E GalvinpersName, Brittany CerbonepersName, Christina A MichelpersName, Henry RusinekpersName, Mony J de LeonpersName, Lidia GlodzikpersName, Susan De SantipersName, P Murali DoraiswamypersName, Jeffrey R PetrellapersName, Terence Z WongpersName, Steven E ArnoldpersName, Jason H KarlawishpersName, David WolkpersName, Charles D SmithpersName, Greg JichapersName, Peter HardypersName, Partha SinhapersName, Elizabeth OatespersName, Gary ConradpersName, Oscar L LopezpersName, MaryAnn OakleypersName, Donna M SimpsonpersName, Anton P PorsteinssonpersName, Bonnie S GoldsteinpersName, Kim MartinpersName, Kelly M MakinopersName, M Saleem IsmailpersName, Connie BrandpersName, Ruth A MulnardpersName, Gaby ThaipersName, Catherine Mc-Adams-OrtizpersName, Kyle WomackpersName, Dana MathewspersName, Mary QuicenopersName, Ramon Diaz-ArrastiapersName, Richard KingpersName, Myron WeinerpersName, Kristen Martin-CookpersName, Michael DeVouspersName, Allan I LeveypersName, James J LahpersName, Janet S CellarpersName, Jeffrey M BurnspersName, Heather S AndersonpersName, Russell H SwerdlowpersName, Liana ApostolovapersName, Kathleen TinguspersName, Ellen WoopersName, Daniel Hs SilvermanpersName, Po H LupersName, George BartzokispersName, Neill R Graff-RadfordpersName, Francine ParfittpersName, Tracy KendallpersName, Heather JohnsonpersName, Martin R FarlowpersName, Ann Marie HakepersName, Brandy R MatthewspersName, Scott HerringpersName, Cynthia HuntpersName, Christopher H van DyckpersName, Richard E CarsonpersName, Martha G MacAvoypersName, Howard ChertkowpersName, Howard BergmanpersName, Chris HoseinpersName, Bojana StefanovicpersName, Curtis CaldwellpersName, Ging-Yuek Robin HsiungpersName, Howard FeldmanpersName, Benita MudgepersName, Michele AssalypersName, Andrew KerteszpersName, John RogerspersName, Charles BernickpersName, Donna MunicpersName, Diana KerwinpersName, Marek-Marsel MesulampersName, Kristine LipowskipersName, Chuang-Kuo WupersName, Nancy JohnsonpersName, Carl SadowskypersName, Walter MartinezpersName, Teresa VillenapersName, Raymond Scott TurnerpersName, Kathleen JohnsonpersName, Brigid ReynoldspersName, Reisa A SperlingpersName, Keith A JohnsonpersName, Gad MarshallpersName, Meghan FreypersName, Barton LanepersName, Allyson RosenpersName, Jared TinklenbergpersName, Marwan N SabbaghpersName, Christine M BeldenpersName, Sandra A JacobsonpersName, Sherye A SirrelpersName, Neil KowallpersName, Ronald KillianypersName, Andrew E BudsonpersName, Alexander NorbashpersName, Patricia Lynn JohnsonpersName, Joanne AllardpersName, Alan LernerpersName, Paula OgrockipersName, Leon HudsonpersName, Evan FletcherpersName, Owen CarmichaelpersName, John OlichneypersName, Charles DeCarlipersName, Smita KitturpersName, Michael BorriepersName, T-Y LeepersName, Rob BarthapersName, Sterling JohnsonpersName, Sanjay AsthanapersName, Cynthia M CarlssonpersName, Steven G PotkinpersName, Adrian PredapersName, Dana NguyenpersName, Pierre TariotpersName, Stephanie ReederpersName, Vernice BatespersName, Horacio CapotepersName, Michelle RainkapersName, Douglas W ScharrepersName, Maria KatakipersName, Anahita AdelipersName, Earl A ZimmermanpersName, Dzintra CelminspersName, Alice D BrownpersName, Godfrey D PearlsonpersName, Karen BlankpersName, Karen AndersonpersName, Robert B SantullipersName, Tamar J KitzmillerpersName, Eben S SchwartzpersName, Kaycee M SinkpersName, Jeff D WilliamsonpersName, Pradeep GargpersName, Franklin WatkinspersName, Brian R OttpersName, Henry QuerfurthpersName, Geoffrey TremontpersName, Stephen SallowaypersName, Paul MalloypersName, Stephen CorreiapersName, Howard J RosenpersName, Bruce L MillerpersName, Jacobo MintzerpersName, Kenneth SpicerpersName, David BachmanpersName, Stephen PasternakpersName, Irina RachinskypersName, Dick DrostpersName, Nunzio PomarapersName, Raymundo HernandopersName, Antero SarraelpersName, Susan K SchultzpersName, Laura L Boles PontopersName, Hyungsub ShimpersName, Karen Elizabeth SmithpersName, Norman RelkinpersName, Gloria ChaingpersName, Lisa RaudinpersName, Amanda SmithpersName, Kristin FargherpersName, Balebail Ashok RajpersName, Thomas NeylanpersName, Jordan GrafmanpersName, Melissa DavispersName, Rosemary MorrisonpersName, Jacqueline HayespersName, Shannon FinleypersName, Karl FriedlpersName, Debra FleischmanpersName, Konstantinos ArfanakispersName, Olga JamespersName, Dino MassogliapersName, J Jay FruehlingpersName, Sandra HardingpersName, Elaine R PeskindpersName, Eric C PetriepersName, Gail LipersName, Jerome A YesavagepersName, Joy L TaylorpersName, Ansgar J FurstpersName (2018date). “Uncovering the heterogeneity and temporal complexity of neurodegenerative diseases with Subtype and Stage Inference.” Nature Communicationstitle·j.

Object of study — phenotypic and temporal heterogeneityterm·dom in neurodegenerative diseases

Research question — how to disentangle disease subtypes from disease stages using cross-sectional data

Contribution — SuStaInterm·mth, a technique recovering subtypes with distinct progression patterns

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

Uncovering the heterogeneity and temporal complexity of neurodegenerative diseases with Subtype and Stage Inferenceterm·mth Abstract The heterogeneity of neurodegenerative diseases is a key confound to disease understanding and treatment development, as study cohorts typically include multiple phenotypes on distinct disease trajectories. Here we introduce a machine-learning technique—Subtype and Stage Inferenceterm·mth (SuStaInterm·mth)—able to uncover data-driven disease phenotypes with distinct temporal progression patterns, from widely available cross-sectional patient studies. Results from imaging studies in two neurodegenerative diseases reveal subgroups and their distinct trajectories of regional neurodegeneration. In genetic frontotemporal dementiaterm·key, SuStaInterm·mth identifies genotypes from imaging alone, validating its ability to identify subtypes; further the technique reveals within-genotype heterogeneity. In Alzheimer’s diseaseterm, SuStaInterm·mth uncovers three subtypes, uniquely characterising their temporal complexity. SuStaInterm·mth provides fine-grained patient stratificationterm·key, which substantially enhances the ability to predict conversion between diagnostic categories over standard models that ignore subtype ( p = 7.18 × 10 −4 ) or temporal stage ( p = 3.96 × 10 −5 ). SuStaInterm·mth offers new promise for enabling disease subtype discovery and precision medicine.

Introduction Neurodegenerative disorders, such as frontotemporal dementiaterm·key (FTD) and Alzheimer’s diseaseterm (AD), are biologically heterogeneous, producing high variance in in vivo disease biomarkers, such as volumetric measurements from imaging, protein measurements from lumbar puncture or behavioural measurements from psychometrics, which reduces their utility in disease studies and management. Key contributors to this heterogeneity are that individuals belong to a range of disease subtypes (giving rise to phenotypic heterogeneityterm·dom) and are at different stages of a dynamic disease process (producing temporal heterogeneityterm·dom). Previous studies aiming to explain biomarker variance typically focus on a single aspect of this heterogeneity: phenotypic heterogeneityterm·dom at a coarse, typically late, disease stage or temporal heterogeneityterm·dom in a broad population. However, the inability to disentangle the range of subtypes from the development and progression of each over time limits the biological insight these techniques can provide, as well as their utility for patient stratificationterm·key. Constructing a comprehensive picture separating phenotypic and temporal heterogeneityterm·dom, i.e. identifying distinct subtypes and characterising the development and progression of each, remains a major current challenge. However, such a picture would provide insights into underlying disease mechanisms, and enable accurate fine-grained patient stratificationterm·key and prognostication, facilitating precision medicine in clinical trials and healthcare.

Both FTD and AD exhibit substantial pathologic, genetic and clinical heterogeneity. In FTD a large proportion of cases (around a third) are inherited on an autosomal-dominant basis, with mutations in progranulin ( GRN ), microtubule-associated protein tau ( MAPT ) and chromosome 9 open reading frame 72 ( C9orf72 ) being the most common causes. Of the major genetic groups, GRN mutations are associated with TDP-43 type A pathology, MAPT mutations with tau inclusions, and expansions in C9orf72 with type A or type B TDP-43 pathology 1 . AD instead has a single pathological characterisation: the presence of both amyloid plaques and neurofibrillary tangles, and the proportion of autosomal-dominant cases is much smaller, accounting for between 1 and 6% of cases 2 . The pathological heterogeneity observed in AD consists of variation in the distribution of neurofibrillary tangles, with 25% of patients having an atypical distribution of neurofibrillary tangles (described as hippocampal-sparing or limbic-predominant) on autopsy at the time of death 3 . Both FTD and AD exhibit a diverse range of clinical syndromes. FTD has both behavioural and language presentations, and in genetic FTD the clinical syndromes can further include atypical parkinsonism and amyotrophic lateral sclerosis. In AD, the major clinical syndrome is broadly divided into amnestic and rarer non-amnestic variants, with non-amnestic variants including language variant AD, logopenic progressive aphasia, visuoperceptive variant AD, posterior cortical atrophy and frontal variant AD 4 .

Previous studies of neurodegenerative disease heterogeneity have focussed on either temporal heterogeneityterm·dom (i.e. subjects appear different at different disease stages) or phenotypic heterogeneityterm·dom (i.e. distinct groups of subjects appear different even at the same disease stage), but rarely both. We refer to these two approaches as stages-only modelsterm·mth, which account for temporal heterogeneityterm·dom but not phenotypic heterogeneityterm·dom, and subtypes-only modelsterm·mth, which account for phenotypic heterogeneityterm·dom but not temporal heterogeneityterm·dom. Stages-only modelsterm·mth arise for example from regression against disease stage 5 , 6 , and data-driven disease progression modellingterm·mth 7 – 15 . Although such models have enabled deeper understanding of the temporal progression of a range of conditions, the inherent assumption that all individuals have a single phenotype, i.e. follow approximately the same trajectory, is a key limitation. At best, this limits the biological insight and the accuracy of stratification they can provide, but potentially could also lead to erroneous conclusions. Subtypes-only modelsterm·mth use, for example, clusteringterm·mth (e.g. refs. 16 – 23 ) to identify distinct groups, or group individuals using information independent of the model, such as genetics (e.g. ref. 24 ) or post-mortem examination (e.g. refs. 25 – 28 ) for models based on in vivo imaging. With typical subtypes-only modelsterm·mth, the limitation is the inherent assumption that all subjects are at a common disease stage so that the cohort has no temporal heterogeneityterm·dom. This requires a priori staging and selection of individuals, which is typically crude in practice leaving models that are not specific to subtype differences. Models of both disease subtype and stage heterogeneity have been constructed previously for the small proportion of neurodegenerative diseases that are inherited on an autosomal-dominant basis. For example, RohrerpersName et al. 29 investigate temporal heterogeneityterm·dom within genetic groups by regressing imaging markers against an estimated age of onset (from family history). However, such studies lack the ability to identify within-genotype phenotypes, and the temporal resolution of the recovered genotype progression patterns is limited by inaccuracy of the a priori staging.

This paper presents Subtype and Stage Inferenceterm·mth (SuStaInterm·mth): a computational technique that disentangles temporal and phenotypic heterogeneityterm·dom to identify population subgroups with common patterns of disease progression. We demonstrate SuStaInterm·mth using structural magnetic resonance imaging (MRI) data sets from cohorts of genetic FTD and AD patients. In each case, SuStaInterm·mth provides a data-driven taxonomy (set of subtypes and stages), as well as detailed pictures of the progression of neurodegeneration within each of the data-driven subgroups. From the genetic FTD data set, SuStaInterm·mth identifies subtypes from imaging alone that map closely onto the genotypes and reconstructs patterns of neurodegeneration that reflect analysis of the individual genetic groups. This provides a validation of SuStaInterm·mth’s ability to identify subgroups with distinct temporal progression patterns, as the different genotypes are known to have distinct patterns of neurodegeneration visible as brain atrophy in MRI 29 . However, SuStaInterm·mth further uncovers two distinct within-genotype phenotypes for carriers of a mutation in the C9orf72 gene, while finding the MAPT and GRN mutation groups are more homogeneous. In AD, SuStaInterm·mth identifies three distinct subtypes and reconstructs their previously unseen temporal progression. In both neurodegenerative diseases, we demonstrate strong assignment of individuals to the SuStaInterm·mth subtypes, which is in contrast to subtypes-only modelsterm·mth in the literature (e.g. ref. 23 ). Even at very early stages, at least a proportion of individuals show strong alignment with particular subtypes, which highlights the potential utility in precision medicine. In AD, we show that SuStaInterm·mth subtype and stage enhance the ability to predict conversion between diagnostic categories substantially beyond subtypes-only or stages-only modelsterm·mth.

Results Subtype and stage inferenceterm·mth Figure 1 provides a conceptual overview of the SuStaInterm·mth modelling technique. SuStaInterm·mth is an unsupervised machine-learning technique that identifies population subgroups with common patterns of disease progression. SuStaInterm·mth builds on and combines ideas from clusteringterm·mth (e.g. refs. 16 – 23 ) and data-driven disease progression modellingterm·mth (e.g. refs. 7 – 10 , 12 ). The combination uniquely enables SuStaInterm·mth to group individuals with common phenotypes across the range of disease stages. It determines the number of subtypes that the available data can support, reconstructs the trajectory of stages within each subtype, and assigns a probability of each subtype and stage to each subject. These features provide insights into the underlying disease biology and a mechanism for in vivo fine-grained stratification at early disease stages. Fig. 1 Conceptual overview of SuStaInterm·mth. The Underlying model panel ( a ) considers a patient cohort to consist of an unknown set of disease subtypes. The input data (Input data panel, b ), which can be entirely cross-sectional, contains snapshots of biomarker measurements from each subject with unknown subtype and unknown temporal stage. SuStaInterm·mth recovers the set of disease subtypes and their temporal progression (as shown in the Output panel, c ) via simultaneous clusteringterm·mth and disease progression modellingterm·mth. Given a new snapshot, SuStaInterm·mth can estimate the probability the subject belongs to each subtype and stage, by comparing the snapshot with the reconstruction (as shown in the Application panel, d ). This figure depicts two hypothetical disease subtypes, labelled I and II, and the biomarkers are regional brain volumes, but SuStaInterm·mth is readily applicable to any scalar disease biomarker and any number of subtypes. The colour of each region indicates the amount of pathology in that region, ranging from white (no pathology) to red to magenta to blue (maximum pathology) Synthetic data A simulation study (see Supplementary Methods, Supplementary Results, Supplementary Discussion and Supplementary Figures 1 – 12 ) verifies the ability of the SuStaInterm·mth algorithm to recover predefined subtypes and their progression patterns from heterogeneous data sets with comparable numbers of subjects, biomarkers and clusters (subtypes) to those used in this study.

Subtype progression patterns We demonstrate SuStaInterm·mth in two neurodegenerative diseases, genetic FTD and sporadic AD, using cross-sectional regional brain volumes from MRI data in the GENetic Frontotemporal dementiaterm·key Initiative (GENFI) and the Alzheimer’s Diseaseterm Neuroimaging Initiative (ADNI). GENFI investigates biomarker changes in carriers of mutations in GRN , MAPT and C9orf72 genes, which cause FTD. GRN and MAPT mutations are known to be associated with distinct phenotypes, whereas C9orf72 is a heterogeneous group 30 . Here, GENFI serves as a test data set with a partially known ground truth for validation, as we expect SuStaInterm·mth to identify genetic groups as distinct phenotypic subtypes. However, it further supports investigation of the phenotypic and temporal heterogeneityterm·dom within genotypes. Specifically, we ran SuStaInterm·mth on the combined data set from all 172 mutation carriers in GENFI (Fig. 2a ), without genotypes, and compared the resulting subtype assignments and progression patterns with (a) participant’s genotype labels (Fig. 2b ), and (b) subtype progression patterns obtained from each genotype separately (Supplementary Figure 13 ; 76 GRN carriers, 63 C9orf72 carriers, 33 MAPT carriers). Next, we used SuStaInterm·mth to identify sporadic AD subtypes from ADNI (793 subjects, including 524 with mild cognitive impairment (MCI) or AD) and characterise their progression from early to late disease stages (Fig. 3 ). We tested consistency of the SuStaInterm·mth subtypes in a largely independent data set—ADNI 1.5T MRI (576 subjects, including 396 with MCI or AD) scans (Fig. 4 ) rather than the main 3T data set used for Fig. 3 . In each disease, cross-validation tests the reproducibility of the subtypes and estimated progression patterns (Supplementary Figure 14 ). Fig. 2 SuStaInterm·mth modelling of genetic frontotemporal dementiaterm·key using GENFI data. a The progression pattern of each of four subtypes that SuStaInterm·mth identifies. Each progression pattern consists of a sequence of stages in which regional brain volumes in mutation carriers (symptomatic and presymptomatic) reach different z -scores relative to non-carriers. Intuitively (for a more precise description see Methods: Uncertainty Estimation), at each stage the colour in each region indicates the level of severity of volume loss: white is unaffected; red is mildly affected ( z -score of 1); magenta is moderately affected ( z -score of 2); and blue is severely affected ( z -score of 3 or more). The circle labelled A indicates the asymmetry of the atrophy pattern (absolute value of the difference in volume between the left and right hemispheres divided by the total volume of the left and right hemispheres) at each stage for each subtype. CVS is the model cross-validation similarity (see Methods: Similarity between two progression patterns): the average similarity of the subtype progression patterns across cross-validation folds, which ranges from 0 (no similarity) to 1 (maximum similarity). f is the proportion of participants estimated to belong to each subtype. b The contribution of each genotype to each of the SuStaInterm·mth subtypes. This is calculated as the probability an individual has a particular genotype given that they belong to a particular subtype Fig. 3 SuStaInterm·mth modelling of sporadic Alzheimer’s diseaseterm using ADNI data. The rows show the progression pattern of the three subtypes identified by SuStaInterm·mth. Diagrams as in Fig. 2 , but the z -scores are measured relative to amyloid-negative (cerebrospinal fluid Aβ1–42 > 192 pg per ml) cognitively normal subjects, i.e. cognitively normal subjects with no evidence of amyloid pathology on cerebrospinal fluid. The cerebellum was not included as a region in the Alzheimer’s diseaseterm analysis and so is shaded in dark grey Fig. 4 Reproducibility of the SuStaInterm·mth subtypes in Fig. 3 . A largely independent Alzheimer’s diseaseterm data set (only 59 subjects are in both the 576 subject data set used to generate this figure and the 793 subject data set used in Fig. 3 ) consisting of those with regional brain volume measurements from 1.5T MRI scans, rather than 3T MRI scans, are shown. Diagrams are as in Fig. 3 , with the rows showing the progression pattern of one of the subtypes identified by SuStaInterm·mth. SuStaInterm·mth modelling identifies three major subtypes: a typical, a cortical and a subcortical subtype, which are in good agreement with the three subtypes in Fig. 3 , as well as an additional very small outlier group (only 4%) with a subtype we term parietal. This small subgroup may represent outliers with a posterior cortical atrophy phenotype SuStaInterm·mth reveals within-genotype phenotypes in FTD Figure 2 shows that SuStaInterm·mth successfully identifies the progression patterns of the different genetic groups in GENFI, without prior knowledge of genotype, and further suggests that phenotypic heterogeneityterm·dom of the C9orf72 group results from two neuroanatomical subtypes. Figure 2a shows the four subtypes that SuStaInterm·mth finds from the full set of all mutation carriers in GENFI. We refer to them as the asymmetric frontal lobe subtype, temporal lobe subtype, frontotemporal lobe subtype and subcortical subtype. Figure 2b reveals that GRN mutation carriers are the main contributors to the asymmetric frontal lobe subtype, MAPT mutation carriers are the main contributors to the temporal lobe subtype, and C9orf72 mutation carriers are the main contributors to both the frontotemporal lobe subtype and the subcortical subtype. This suggests that there are two distinct subtypes in the C9orf72 group. Application of SuStaInterm·mth to each genetic group separately supports this finding by demonstrating that the GRN mutation carriers are best described as a single asymmetric frontal lobe subtype, the MAPT mutation carriers are best described as a temporal lobe subtype and the C9orf72 mutation carriers are best described as two distinct disease subtypes: a frontotemporal lobe subtype and a subcortical subtype. SuStaInterm·mth additionally finds a subsidiary cluster in the MAPT group for which the progression pattern has high uncertainty. This high uncertainty likely prevents the cluster from being detected when applying SuStaInterm·mth to all mutation carriers in Fig. 2 as this small number of subjects can be sufficiently modelled by the three alternative subtype progression patterns. Supplementary Figure 13 shows that the subtype progression patterns for each genetic group are in good agreement with those found in the full set of all mutation carriers (Fig. 2a ). Supplementary Figure 14A shows that the four subtypes estimated in Fig. 2a are reproducible under cross-validation, with a high average similarity between cross-validation folds of >93% for each subtype. Altogether these results provide strong validation of SuStaInterm·mth’s ability to recover distinct subtypes and their progression patterns from a heterogeneous data set, while simultaneously disentangling the heterogeneity of the C9orf72 group into two distinct subtypes.

SuStaInterm·mth identifies three subtype progression patterns in AD Figure 3 shows the temporal progression of the three neuroanatomical subtypes that SuStaInterm·mth identifies from ADNI, which we term typical, cortical and subcortical. SuStaInterm·mth reveals that for the typical subtype, atrophy starts in the hippocampus and amygdala; for the cortical subtype in the nucleus accumbens, insula and cingulate; and for the subcortical subtype in the pallidum, putamen, nucleus accumbens and caudate. Supplementary Figure 14B shows that these three subtypes are reproducible under cross-validation, giving an average similarity between cross-validation folds of >92% for each subtype.

AD subtypes are reproducible in an independent data set Figure 4 shows that the three subtypes in Fig. 3 are reproducible in a largely independent data set (<5% subjects in common) consisting of regional brain volumes derived from 1.5T rather than 3T MRI scans. From the 1.5T data, SuStaInterm·mth broadly replicates the three major clusters found in the 3T data, again finding a typical, cortical and subcortical subtype. The origin of atrophy for each subtype is in general agreement with the 3T data: atrophy begins in the hippocampus and amygdala for the typical subtype, in the insula and cingulate for the cortical subtype; and in the pallidum, putamen and caudate for the subcortical subtype. The main difference compared to the 3T data is that the nucleus accumbens is not indicated as an early region to atrophy in the 1.5T data for the cortical and subcortical subtypes. SuStaInterm·mth additionally identifies a small proportion (4%) of outliers with a parietal subtype in the 1.5T data.

Disease subtyping and staging We investigated SuStaInterm·mth’s capability for reliable stratification in each neurodegenerative disease (Fig. 5 ) to determine the potential for homogeneous cohort identification. First, we assessed how reliably SuStaInterm·mth assigns patients to subtypes (Fig. 5a, b ). Specifically, in genetic FTD, we tested the consistency of SuStaInterm·mth subtypes with the different genotypes in symptomatic mutation carriers (Table 1 ), and compared this consistency against models that do not account for temporal heterogeneityterm·dom (Table 2 ). Second, we assessed the reliability of the SuStaInterm·mth stages in each disease (Fig. 5c, d ) by comparison with clinical diagnostic categories. In ADNI, where clinical follow-up information is available, we further examined the ability of SuStaInterm·mth subtypes and stages to predict relevant outcomes, by determining whether SuStaInterm·mth subtype and/or stage modify the risk of conversion between diagnostic categories (Table 3 ). Fig. 5 SuStaInterm·mth subtyping and staging of genetic frontotemporal dementiaterm·key and Alzheimer’s diseaseterm. a , b The assignability of the disease subtypes estimated by SuStaInterm·mth for genetic frontotemporal dementiaterm·key, and Alzheimer’s diseaseterm. Each scatter plot visualises the probability that each individual belongs to each of the SuStaInterm·mth subtypes estimated for a genetic frontotemporal dementiaterm·key (as shown in Fig. 2a ), and b . Alzheimer’s diseaseterm (as shown in Fig. 3 ). In the triangle scatter plots, each of the corners corresponds to a probability of 1 of belonging to that subtype, and 0 for the other subtypes; the centre point of the triangle corresponds to a probability of 1/3 of belonging to each subtype. c , d The probability subjects from each of the diagnostic groups belong to each of the SuStaInterm·mth stages for c genetic frontotemporal dementiaterm·key and d Alzheimer’s diseaseterm. CN = cognitively normal; MCI = mild cognitive impairment; AD = Alzheimer’s diseaseterm Table 1 Ability of subtypes to distinguish between different genotypes in symptomatic mutation carriers in GENFI using the SuStaInterm·mth subtypes in Fig. 2a GRN MAPT C9orf72 Asymmetric frontal (threshold p > 0.65) 93% (13) 9% (1) 4% (1) Temporal (threshold p > 0.35) 0% (0) 91% (10) 21% (5) Frontotemporal 0% (0) 0% (0) 42% (10) Subcortical 7% (1) 0% (0) 33% (8) Accuracy 93% (13/14) 91% (10/11) 75% (18/24) Each entry is the percentage (number) of participants of a particular genotype assigned to that subtype. The final row indicates the percentage (fraction) of participants assigned to the correct subtype from each genotype. The results show that SuStaInterm·mth can accurately discriminate genotypes, validating the ability of SuStaInterm·mth to identify distinct phenotypes that align with known genetic groups Table 2 As Table 1, but for subtypes obtained from a subtypes-only modelterm·mth that accounts for heterogeneity in disease subtype but not disease stage, i.e. the subtypes in Fig. 6 GRN MAPT C9orf72 Severe frontal (threshold p > 0.99) 57% (8) 9% (1) 4% (1) Severe temporal (threshold p > 0.99) 0% (0) 64% (7) 8% (2) Mild frontotemporal 43% (6) 27% (3) 88% (21) Accuracy 57% (8/14) 64% (7/11) 88% (21/24) The results show that SuStaInterm·mth (Table 1 ) provides much better discrimination of the different genotypes than the subtypes-only modelterm·mth shown here, demonstrating the added utility of a model that accounts for heterogeneity in disease stage Table 3 Utility of SuStaInterm·mth subtype and stage for predicting the risk of conversion from mild cognitive impairment to Alzheimer’s diseaseterm SuStaInterm·mth subtype SuStaInterm·mth stage Age Sex Education APOE4 S–C–T 1.57** 1.13 † 0.98 0.98 0.93 ~ 1.82 † S–C 1.76 ~ 1.16 † 0.95* 1.03 0.92 1.53* C–T 1.48 ~ 1.11 † 0.98 0.87 0.97 1.84 † S–T 2.11* 1.13 † 1.02 1.13 0.90* 2.13 † Each row shows Hazards ratiosterm·mth for a different Cox Proportional Hazards model that estimates the risk of conversion from mild cognitive impairment to Alzheimer’s diseaseterm using ADNI data. Each column shows the estimated hazard ratio for each variable. Each hazards ratioterm·mth tells you how the risk of conversion changes for each unit increase of a particular variable: a ratio of 1 means no modification of the risk, a ratio >1 means there is an increase of the risk, and a ratio less than 1 means there is a reduction of the risk. For the first model (S–C–T) it is assumed that the hazard ratio increases multiplicatively from the Subcortical subtype (S) to the Cortical subtype (C) to the Typical subtype (T), i.e. the S–C–T model estimates that each SuStaInterm·mth subtype has a hazards ratioterm·mth 1.57 times that of the previous subtype (i.e. the cortical group have a 1.57 times greater risk of conversion than the subcortical group, the typical group have a 1.57 times greater risk of conversion than the cortical group, and the typical group have a 2.46 (1.57 2 ) times greater risk of conversion than the subcortical group). In the remaining models only two groups are compared at a time to demonstrate that the results are similar without this assumption, although the statistical power is reduced. This result demonstrates the added utility of both disease subtypes and stages obtained from SuStaInterm·mth for predicting conversion between mild cognitive impairment and Alzheimer’s diseaseterm, with both subtype and stage modifying the risk of conversion. Statistical significance is indicated as: ~ p < 0.1, * p < 0.05, ** p < 0.01, † p < 1 × 10 −3 SuStaInterm·mth provides utility for patient stratificationterm·key Figure 5 illustrates the ability of SuStaInterm·mth to provide disease subtyping and staging information for each neurodegenerative disease. Figure 5a shows that the strength of assignment (see Methods: Strength of assignment to subtype) to the SuStaInterm·mth subtypes in genetic FTD increases as the diseases progress, with 88% of the symptomatic mutation carriers in GENFI being strongly assigned (i.e. >50% likelihood of a particular subtype). Figure 5b shows that the strength of assignment to the SuStaInterm·mth subtypes in AD also increases with disease progression, with a strong assignment of individuals to subtypes in 78% of ADNI participants with an AD diagnosis. The strong assignment of the AD subtypes that SuStaInterm·mth achieves by accounting for temporal heterogeneityterm·dom is in contrast to previous studies 23 that model phenotypic but not temporal heterogeneityterm·dom. Moreover, the strong assignment is seen even at early disease stages (MCI), where many subjects cluster around the vertices of the triangles: 37% of MCI subjects are strongly assigned to a subtype. Figures 5c, d show that the distribution of SuStaInterm·mth stages differs between diagnostic groups in both GENFI and ADNI, and provides a good separation of presymptomatic and symptomatic mutation carriers, and cognitively normal (CN) and AD.

SuStaInterm·mth subtypes discriminate FTD genotype Table 1 shows the classification accuracy obtained using the SuStaInterm·mth subtypes in Fig. 2 to discriminate the genotype of affected mutation carriers in GENFI. While the use of MRI to identify genotype is not necessary in these subjects, so not clinically relevant, this experiment demonstrates the ability of SuStaInterm·mth to identify subtypes in a data set with a known ground truth. The SuStaInterm·mth subtypes give a balanced accuracy of 95% for the two-way classification task of distinguishing the homogeneous GRN and MAPT carrier groups. For the more challenging three-way classification task of distinguishing all genotypes in the presence of heterogeneity, the SuStaInterm·mth subtypes provide a maximum balanced accuracy of 86%. A high proportion of the homogeneous GRN and MAPT carrier groups are correctly assigned to the asymmetric frontal lobe (93% of affected GRN carriers) and temporal lobe subtype progression patterns (91% of affected MAPT carriers). The heterogeneous C9orf72 carrier group are much more difficult to classify, with a total of 75% of affected C9orf72 carriers being assigned to the frontotemporal lobe and subcortical subtypes. Apart from heterogeneity, the C9orf72 carriers are also more difficult to classify because the frontotemporal lobe and subcortical subtype progression patterns are more similar to the other subtypes; by evaluating the similarity of each pair of subtype progression patterns (see Methods: Similarity between two subtype progression patterns) we find that the asymmetric frontal lobe and temporal lobe subtypes have the most distinct progression patterns of any pair of subtypes; the asymmetric frontal lobe and frontotemporal lobe subtypes have the most similar progression patterns of any pair of subtypes. The precise strategy of assigning subjects to subtype can alter the classification rates somewhat and Supplementary Table 1 examines this effect.

Genotype discrimination out-performs subtypes-only modelsterm·mth Table 2 shows the classification accuracy obtained using a subtypes-only modelterm·mth (Fig. 6 ), which does not account for temporal heterogeneityterm·dom, to discriminate the genotype of affected mutation carriers in GENFI. The SuStaInterm·mth subtypes out-perform the subtypes-only modelterm·mth. The subtypes-only modelterm·mth gives a balanced accuracy of 92% compared to 95% using SuStaInterm·mth for the two-way classification task of distinguishing GRN and MAPT carrier groups; the subtypes-only modelterm·mth gives a maximum balanced accuracy of 69% compared to 86% using SuStaInterm·mth for the three-way classification task of distinguishing all genotypes. In the subtypes-only modelterm·mth the majority of misclassifications arise from the earlier stage affected GRN and MAPT carriers being assigned to the mild frontotemporal subtype associated with C9orf72 carriers. See also Supplementary Table 2 . Fig. 6 Subtypes-only modelterm·mth for GENFI; not accounting for disease stage heterogeneity. Brain diagrams as in Fig. 2 , but here each diagram represents a different subtype, which we refer to as severe frontal, severe temporal and mild frontotemporal. There is no notion of disease stage in the subtypes-only modelterm·mth SuStaInterm·mth subtypes and stages have predictive utility in AD Table 3 shows that the SuStaInterm·mth subtypes and stages have predictive utility for the risk of conversion between diagnostic categories in ADNI. By fitting a Cox Proportional Hazards model, we found significant effects ( t -test) of baseline SuStaInterm·mth subtype ( p = 2.44 × 10 −3 ) and stage ( p = 8.76 × 10 −11 ) on an individual’s risk of conversion from MCI to AD. Of the SuStaInterm·mth subtypes, the subcortical subtype is associated with the lowest risk of conversion, while the typical subtype is associated with the highest risk of conversion. Supplementary Table 3 shows that SuStaInterm·mth out-performs subtypes-only and stages-only modelsterm·mth at estimating the risk of conversion between diagnostic categories in ADNI. By performing likelihood ratio tests comparing SuStaInterm·mth to subtypes-only and stages-only we find that SuStaInterm·mth provides a significantly better fit (likelihood ratio test) than both subtypes-only ( p = 3.96 × 10 −5 ) and stages-only ( p = 7.18 × 10 −4 ) models. This shows that both the subtypes and stages estimated by SuStaInterm·mth provide additional information for predicting the risk of conversion from MCI to AD.

Discussion In this study we introduce SuStaInterm·mth—a powerful tool for data-driven disease phenotype discovery, providing insights into disease aetiology, and enhanced power for patient stratificationterm·key in clinical trials and healthcare. Results from the GENFI data set first validate that SuStaInterm·mth can successfully recover known distinct progression patterns in genetic FTD corresponding to different genotypes. Moreover, SuStaInterm·mth identifies and characterises within-group heterogeneity for carriers of a mutation in the C9orf72 gene as distinct temporal progression patterns in two subtypes. The results demonstrate the utility of SuStaInterm·mth for data-driven disease phenotype discovery, and provide biological insight into the C9orf72 mutation. Application of SuStaInterm·mth to the 3T ADNI data set recovers three distinct AD subtypes with final stages that reflect post-mortem neuropathological findings. Results from a largely independent AD data set (ADNI 1.5T) corroborate these three subtypes. The disease subtype characterisation SuStaInterm·mth provides goes much further than post-mortem neuropathological studies 3 , 28 , or other machine-learning techniques 23 , by characterising the temporal trajectory of each subtype, enabling in vivo stratification of subjects by disease stage as well as disease subtype. We demonstrate the ability of SuStaInterm·mth to stratify in vivo by both subtype and stage in genetic FTD and AD. In genetic FTD, we show that the SuStaInterm·mth neuroimaging subtypes can distinguish affected carriers belonging to different genetic groups with high classification accuracy. In AD, we demonstrate strong assignment of subjects to the SuStaInterm·mth subtypes even at early disease stages (MCI), and that the SuStaInterm·mth subtypes and stages have added utility for predicting conversion between clinical diagnoses, beyond more traditional stages-only or subtypes-only modelsterm·mth.

Previous studies in genetic FTD have found asymmetric frontotemporoparietal lobe volume loss in GRN carriers, temporal lobe volume loss in MAPT carriers, and widespread symmetric grey matter atrophy and volume loss in the cerebellum in C9orf72 carriers 24 . The asymmetric frontal lobe subtype and temporal lobe subtype in Fig. 2 show clear similarities with previous studies of regional volume loss in GRN and MAPT mutation carriers respectively. However, SuStaInterm·mth provides much greater detail and accuracy by avoiding reliance on crude a priori staging, e.g. via mean familial age of onset. The frontotemporal lobe subtype and subcortical subtype in Fig. 2 both have features previously associated with C9orf72 mutation carriers, but SuStaInterm·mth assigns these features to two distinct disease subtypes, and further reveals the temporal progression of each subtype.

Several biological factors may produce the two subtypes observed in C9orf72 mutation carriers, either individually or in combination. Clinically, while there is significant overlap, patients typically present with either a behavioural variant FTD or amyotrophic lateral sclerosis as their main phenotype 31 , and they can progress at various rates; genetically, the expansion length is variable and there are additional genetic modifiers (e.g. TMEM106B and ATXN2 ) that alter phenotype 32 – 34 ; and pathologically, most cases have either type A or type B TDP-43 pathology 31 . While further study is required to determine the biological factors that influence neuroanatomical phenotype, these findings demonstrate the power of SuStaInterm·mth in identifying hitherto unrecognised disease subtypes using clinical data, and thus open up the potential to link variations in genetics, pathology and neuroanatomy.

We also find evidence for the presence of a subsidiary group in the MAPT mutation carriers, but numbers are too small to determine whether this group have a distinct progression pattern. Among individuals with significant evidence of MRI atrophy (SuStaInterm·mth stage of ≥5), four individuals (two pairs of individuals from the same families) of 13 were identified as belonging to the subsidiary group. Although MAPT mutations have been commonly thought to have a very specific pattern of atrophy affecting the anterior and medial temporal lobes predominantly, one previous paper has shown that there can be a second pattern of atrophy in specific mutations, where the lateral temporal lobes are affected more than the medial regions 35 . Interestingly, the two pairs of individuals who constitute the subsidiary group in our analysis all have P301L mutations, a mutation that falls into this second alternate atrophy pattern group in ref. 35 . None of the nine individuals assigned to the predominant progression pattern in our analysis have P301L mutations, or V337M mutations, the other mutation identified in ref. 35 as having an alternate atrophy pattern. This suggests that SuStaInterm·mth may be able to identify particular MAPT mutations that fall into this alternate group, but larger studies will be required to confirm this.

In AD, post-mortem histology 3 and retrospectively analysed MRI scans close to the time of death 28 observe three distinct patterns of atrophy in late-stage AD patients: one focussed on the temporal lobe that is similar to the late stages of the typical SuStaInterm·mth subtype; one affecting predominantly cortical regions cf. late stages of the cortical SuStaInterm·mth subtype; and one with stronger subcortical involvement cf. late stages of the subcortical SuStaInterm·mth subtype. This gives confidence in the SuStaInterm·mth subtypes, which provide much greater information by revealing the progression of each subtype over time, including the earliest sites of regional volume loss. Moreover, and importantly for practical utility, the SuStaInterm·mth subtypes can be assigned in vivo using MRI, enabling linkage of late-stage pathological observations with early-stage neurodegeneration.

The three AD subtypes found in the 3T MRI data set are corroborated by the largely independent 1.5T MRI data set. However, some small differences arise between the subtype progression patterns of the three subtypes recovered in each data set. These differences lie predominantly in how early the nucleus accumbens begins to atrophy in the different subtypes: across all three subtypes the nucleus accumbens shows atrophy earlier in the 3T subtypes than the 1.5T subtypes. A possible explanation for this is that the volume of the nucleus accumbens, which is relatively small, can be estimated more accurately using the higher field strength 3T MRI scans than the 1.5T MRI scans, and thus atrophy in the nucleus accumbens can be identified from an earlier stage in the 3T data set compared to the 1.5T data set.

In the 1.5T MRI data set we additionally find a small proportion (4%) of outliers with a parietal subtype. This small subgroup may represent a posterior cortical atrophy phenotype: comparing the Alzheimer’s diseaseterm Assessment Scale-cognitive subscale (ADAS-cog) scores between individuals with an AD diagnosis that are assigned to the parietal subgroup ( N = 6) and the typical AD subgroup ( N = 65), we find that the parietal subgroup have worse performance (Mann–Whitney U test) on certain praxic (Q6. Ideational Praxis, p = 6.1 × 10 −3 , z = 2.7) and spatially-demanding (Q14. Number Cancellation, p = 4.9 × 10 −3 , z = 2.8) subtests, but similar performance (Mann–Whitney U test) in memory domains (Q8. Word Recognition, p = 0.81, z = −0.2; Q1. Word Recall, p = 0.48, z = 0.70). Additionally, the parietal subgroup is on average 10.3 years younger ( p = 2.8 × 10 −3 , z = −3.0, Mann–Whitney U test) than the typical AD subgroup.

The temporal spreading patterns for distinct subtypes estimated by SuStaInterm·mth offer biological insight. For example, the progression pattern of each subtype provides a view of how neurodegeneration spreads from a distinct origin over the rest of the brain that is uncorrupted by phenotypic heterogeneityterm·dom. A key advantage of SuStaInterm·mth is that it provides a purely data-driven, hypothesis-free, reconstruction of the progression of neurodegenerative disease subtypes. However, these observations also have great potential to inform mechanistic models 36 , 37 of neurodegenerative disease, which explain their temporal progression via various hypothetical mechanisms of disease propagation over brain networks. Current mechanistic models implicitly assume a single-disease progression pattern—an assumption often violated in patient data sets, but much more reasonable if focussed on particular SuStaInterm·mth subtypes.

SuStaInterm·mth shows strong capabilities for patient stratificationterm·key in AD, which we are able to validate in genetic FTD where we expect the subtype assignments to correspond to distinct genotypes. SuStaInterm·mth provides high classification accuracy for differentiating the different mutation types in genetic FTD, and the AD subtypes are clearly assignable. In genetic FTD, SuStaInterm·mth out-performs a subtypes-only modelterm·mth, giving a balanced classification accuracy of 86% for distinguishing genotype compared to 69% for the subtypes-only modelterm·mth. This provides compelling evidence that there is substantial heterogeneity in disease stage within different phenotypes, and that modelling this disease stage heterogeneity is important for better patient stratificationterm·key. This is further demonstrated in AD, in which SuStaInterm·mth’s subtypes and stages substantially out-perform subtypes-only and stages-only modelsterm·mth for predicting conversion between diagnostic categories. These early results are highly promising, particularly given that the particular choice of biomarkers used here (coarse regional brain volumes) is not optimised for stratification. In this initial study, we chose to use MRI to maximise the number of subjects with all available measurements, to simplify the interpretation of the results, and to enhance the clinical utility. However, future work will test the added benefit of including a wider range of biomarkers and a more fine-grained set of regional volumes in SuStaInterm·mth for patient stratificationterm·key. For example in AD, incorporation of amyloid and neurofibrillary tangle measures, e.g. from amyloid and tau positron emission tomography (PET) scans, will enable stratification of individuals at the very earliest disease stages.

The previous study of Zhang et al. 23 also looked at the assignment of AD subtypes using a subtypes-only modelterm·mth that does not account for temporal heterogeneityterm·dom in disease stage. In contrast to the study of Zhang et al. 23 , we observe strong assignment of AD patients to the subtypes (Fig. 5b ) emphasising the importance of accounting for heterogeneity in disease stage. This assignability clearly increases with disease progression, with the subtypes being most strongly assigned in clinically diagnosed AD patients. However, even at early stages (MCI), many subjects cluster around the vertices of the triangles showing strong potential for identifying early-stage cohorts representative of each subtype.

The model underlying SuStaInterm·mth makes several assumptions to enable the simultaneous estimation of subtypes and their progression. One assumption is that biomarker variance is independent. In reality biomarkers tend to co-vary due to shared biological processes. However, simulation experiments (Supplementary Figure 7 ) show that the subtype progression patterns recovered by SuStaInterm·mth are robust to biomarker covariance. Nevertheless, refinements might come from modelling covariance among strongly dependent biomarkers, e.g. using model selection criteria to identify a minimal set of necessary covariance parameters, and future work will explore this idea.

To enable the modelling of purely cross-sectional data, here we make an assumption of an arbitrary timescale. This formulation can also work with longitudinal data when available, although here we reserve the longitudinal information to validate the clinical utility of SuStaInterm·mth to make future predictions at an individual’s first visit. However, extensions to SuStaInterm·mth that utilise longitudinal information to further provide a well-defined timescale are an important area for future work.

Here we make a further implicit assumption that the cohort is correctly diagnosed. While the genetic tests in GENFI ensure this, the clinical diagnoses of AD in ADNI are less reliable and the proportion of misdiagnosis, e.g. of depression or other neurological diseases, is non-negligible. Simulation of the effect of misdiagnosis (Supplementary Figure 9 ) demonstrates that SuStaInterm·mth can robustly recover subtype progression patterns under a substantial proportion of outliers, up to 20%. Nevertheless, future adaptions of the SuStaInterm·mth model might include a broad outlier class to capture individuals that do not fit any of the main clusters.

One caveat on the findings is that the underlying data may come from a spectrum of disease progression patterns, rather than a set of distinct trajectories as SuStaInterm·mth is designed to estimate. Simulations (Supplementary Figure 11 ) demonstrate that SuStaInterm·mth may still recover multiple distinct progression patterns from data generated by a spectrum of progression patterns. In this case the distinct progression patterns identified by SuStaInterm·mth still provide useful information about the extrema within the underlying spectrum of progression patterns. Here, however, the alignment with genotypes in genetic FTD and neuropathological observations in AD provide confidence that the distinct subtypes are genuine. Future work will extend SuStaInterm·mth to be able to represent spectra of progression patterns (e.g. using Mallow’s models as in refs. 38 , 39 ).

We introduce SuStaInterm·mth—a tool to disentangle and characterise the temporal and phenotypic heterogeneityterm·dom of neurodegenerative diseases. We use it to elucidate the temporal and phenotypic heterogeneityterm·dom of both genetic FTD and AD subtypes with previously unseen detail. We further demonstrate SuStaInterm·mth’s potential as a patient stratificationterm·key tool in AD by showing strong alignment of subjects with specific subtypes even at early disease stages, as well as added power to predict conversion between clinical diagnoses. SuStaInterm·mth has the potential to make substantial clinical impact as a tool for precision medicine and is readily applicable to any progressive disease, including other neurodegenerative diseases, respiratory diseases and cancers.

Methods GENFI data set We used cross-sectional volumetric MRI data from GENFI ( http://www.genfi.org.uk/ ). Subjects were included from the second data freeze of GENFI, which in total consisted of 365 participants recruited across 13 centres in the United Kingdom, Canada, Italy, The Netherlands, Sweden and Portugal. A total of 313 participants had a usable volumetric T1-weighted MRI scan for analysis (15 participants did not have a scan and the other participants were excluded as the scans were of unsuitable quality due to motion, other imaging artefacts or pathology unlikely to be attributed to FTD). The 313 participants included 141 non-carriers, 123 presymptomatic carriers and 49 symptomatic carriers. Of the 123 presymptomatic mutation carriers there were 62 GRN , 39 C9orf72 and 22 MAPT carriers. Of the 49 symptomatic carriers, there were 14 GRN , 24 C9orf72 and 11 MAPT carriers. The acquisition and post-processing procedures for GENFI have been previously described in ref. 29 . Briefly, cortical and subcortical volumes were generated using a multi-atlas segmentation propagation approach 40 , combining cortical regions of interest to calculate grey matter volumes of the entire cortex, separated into the frontal, temporal, parietal, occipital, cingulate and insula cortices. In addition to regional volumetric measures, we also included a measure of asymmetry, which is calculated as the absolute value of the difference between the volumes of the right and left hemispheres, normalised by the total volume of both hemispheres. This asymmetry measure was log transformed to improve normality. See Supplementary Table 4 for a summary of the biomarkers used in the SuStaInterm·mth modelling.

ADNI data set Data used in the preparation of this article were obtained from the Alzheimer’s Diseaseterm Neuroimaging Initiative (ADNI) database ( http://adni.loni.usc.edu ). The ADNI was launched in 2003 by the National Institute on Aging (NIA), the National Institute of Biomedical Imaging and Bioengineering (NIBIB), the Food and Drug Administration (FDA), private pharmaceutical companies and non-profit organisations, as a $60 million, 5-year public-private partnership. For up-to-date information, see http://www.adni-info.org . Written consent was obtained from all participants, and the study was approved by the Institutional Review Board at each participating institution.

We downloaded data from Laboratory Of Neuro Imaging (LONI; http://adni.loni.usc.edu ) on 11 May 2016 and constructed two cross-sectional volumetric MRI data sets for SuStaInterm·mth model fitting: those with higher (3T) and lower (1.5T) field strength. The inclusion criteria for the 3T and 1.5T data sets were having cross-sectional FreeSurfer volumes available that passed overall quality control from either a 3T (processed using FreeSurfer Version 5.1) or a 1.5T (processed using FreeSurfer Version 4.3) MRI scan. The 3T data set consisted of 793 subjects (183 CN, 86 significant memory concern, 243 early MCI, 164 late MCI, 117 AD), of which 73 were enroled in ADNI-1, 99 were enroled in ADNI-GO and 621 were enroled in ADNI-2. The 1.5T data set consisted of 576 ADNI-1 subjects (180 CN, 274 late MCI, 122 AD). The 1.5T and 3T data sets are largely independent: only 59 subjects (14 CN, 33 late MCI, 12 AD) have both 1.5T and 3T scans. We downloaded processed cross-sectional FreeSurfer volumes for 1.5T and 3T scans, using FreeSurfer Versions 4.3 and 5.1, and quality control ratings. We retained only the volumes that passed overall quality control, and normalised them by regressing against total intracranial volume. We further downloaded demographic information for covariate correction: age, sex, education and APOE genotype from the ADNIMERGE table. We downloaded diagnostic follow-up information to test the association of the SuStaInterm·mth model subtypes and stages with longitudinal outcomes. We also downloaded baseline cerebrospinal fluid (CSF) measurements of Aβ1–42, which we used to identify a control population. Again, see Supplementary Table 4 for a summary of the biomarkers used in the SuStaInterm·mth modelling.

z -scores We expressed each regional volume measurement as a z -score relative to a control population: in GENFI we used data from all non-carriers, in ADNI we used amyloid-negative CN subjects, defined as those with a CSF Aβ1–42 measurement >192 pg per ml 41 . This gave us a control population of 48 amyloid-negative CN subjects for the 3T data set, and 56 amyloid-negative CN subjects for the 1.5T data set. We used these control populations to determine whether the effects of age, sex, education or number of APOE4 alleles (ADNI only) were significant, and if so to regress them out. We then normalised each data set relative to its control population, so that the control population had a mean of 0 and standard deviation of 1. Because regional brain volumes decrease over time the z -scores become negative with disease progression, so for simplicity we took the negative value of the z -scores so that the z -scores would increase as the brain volumes became more abnormal.

SuStaInterm·mth modelling We formulate the model underlying SuStaInterm·mth as groups of subjects with distinct patterns of biomarker evolution (see Mathematical Model). We refer to a group of subjects with a particular biomarker progression pattern as a subtype. The biomarker evolution of each subtype is described as a linear z -score model in which each biomarker follows a piecewise linear trajectory over a common timeframe. The noise level for each biomarker is assumed constant over the timeframe and is derived from a control population (see Mathematical model). This linear z -score model is based on the event-based model in refs. 7 , 8 , 38 , but reformulates the events so that they represent the continuous linear accumulation of a biomarker from one z -score to another, rather than an instantaneous switch from a normal to an abnormal level. A key advantage of this formulation is that it can work with purely cross-sectional data because it requires no information about the timescale of change, but instead uses events as control points of piecewise linear segments with arbitrary duration. The model fitting considers increasing number of subtypes C , for which we estimate the proportion of subjects f that belong to each subtype, and the order S C in which biomarkers reach each z -score for each subtype c = 1 … C . We determine the optimal number of subtypes C for a particular data set through ten-fold cross-validation (see Cross-validation).

Mathematical model The linear z -score model underlying SuStaInterm·mth is a continuous generalisation of the original event-based model 7 , 8 , which we describe first.

The event-based model in refs. 7 , 8 describes disease progression as a series of events, where each event corresponds to a biomarker transitioning from a normal to an abnormal level. The occurrence of an event, E i , for biomarker i = 1 … I , is informed by the measurements x ij of biomarker i in subject j , j = 1 … J . The whole data set X = { x ij | i = 1 … I , j = 1 … J } is the set of measurements of each biomarker in each subject. The most likely ordering of the events is the sequence S that maximises the data likelihood 1 P X ∣ S = ∏ j = 1 J ∑ k = 0 I P ( k ) ∏ i = 1 k P x i j ∣ E i ∏ i = k + 1 I P x i j ∣ ¬ E i , where P ( x | E i ) and P ( x | ¬ E i ) are the likelihoods of measurement x given that biomarker i has or has not become abnormal, respectively. P ( k ) is the prior likelihood of being at stage k , at which the events E 1 , ..., E k have occurred, and the events E k +1 , …, E I have yet to occur. The model uses a uniform prior on the stage, so that P ( k ) = 1/( I + 1), k = 0 … I , i.e. a priori individuals are equally likely to belong to any stage along the progression pattern. The likelihoods P ( x | E i ) and P ( x | ¬ E i ) are modelled as normal distributions.

The linear z -score model we use in this work reformulates the event-based model in (1) by replacing the instantaneous normal to abnormal events with events that represent the (more biologically plausible) linear accumulation of a biomarker from one z -score to another. The linear z -score model consists of a set of N z -score events E iz , which correspond to the linear increase of biomarker i = 1 … I to a z -score z ir = z i 1 … z i R i , i.e. each biomarker is associated with its own set of z -scores, and so N = ∑ i R i . Each biomarker also has an associated maximum z -score, z max , which it accumulates to at the end of stage N . We consider a continuous time axis, t , which we choose to go from t = 0 to t = 1 for simplicity (the scaling is arbitrary). At each disease stage k , which goes from t = k N + 1 to t = k + 1 N + 1 , a z -score event E iz occurs. The biomarkers evolve as time t progresses according to a piecewise linear function g i ( t ), where g t = z 1 t E z 1 t , 0 < t ≤ t E z 1 z 1 + z 2 - z 1 t E z 2 - t E z 1 t - t E z 1 , t E z 1 < t ≤ t E z 2 ⋮ z R - 1 + z R - z R - 1 t E z R - t E z R - 1 t - t E z R - 1 , t E z R - 1 < t ≤ t E z R z R + z m a x - z R 1 - t E z R t - t E z R , t E z R < t ≤ 1 . Thus, the times t E i z are determined by the position of the z -score event E iz in the sequence S , so if event E iz occurs in position k in the sequence then t E i z = k + 1 N + 1 .

To formulate the model likelihood for the linear z -score model we replace Eq. ( 1 ) with 2 P X ∣ S = ∏ j = 1 J ∑ k = 0 N ∫ t = k N + 1 t = k + 1 N + 1 P ( t ) ∏ i = 1 I P x i j ∣ t ∂ t , where, P x i j ∣ t = NormPDF x i j , g i t , σ i . NormPDF( x , μ , σ ) is the normal probability distribution function, with mean μ and standard deviation σ , evaluated at x . We assume the prior on the disease time is uniform, as in the original event-based model.

The SuStaInterm·mth model is a mixture of linear z -score models, hence we have P X ∣ M = ∑ c = 1 C f c P X ∣ S c , where C is the number of clusters (subtypes), f is the proportion of subjects assigned to a particular cluster (subtype), and M is the overall SuStaInterm·mth model.

Model fitting Supplementary Figure 15 provides a flowchart detailing the processes involved in the SuStaInterm·mth model fitting. Model fitting requires simultaneously optimising subtype membership, subtype trajectory and the posterior distributions of both. In particular, the cost function here depends on the sequence ordering, which to our knowledge standard algorithms do not handle. We therefore derive our own algorithm to fit SuStaInterm·mth, based on the well-established methods developed for the event-based model ( 7 , 8 , 42 , 43 ), for which we demonstrate convergence and optimality in simulation (see Supplementary Results: Convergence) and in the data sets used here (see Convergence). As shown in the blackpersName box in Supplementary Figure 15 , the SuStaInterm·mth model is fitted hierarchically, with the number of clusters being estimated via model selection criteria obtained from cross-validation. The hierarchical fitting initialises the fitting of each C -cluster (subtype) model from the previous C-1 -cluster model, i.e. the clusteringterm·mth problem is solved sequentially from C = 1 … C max (where C max is the maximum number of clusters being fitted), initialising each model using the previous model. For the initial cluster ( C = 1), we use the single-cluster expectation maximisation (E-M) procedure shown in the greenpersName box in Supplementary Figure 15 , and described subsequently. We fit subsequent cluster numbers ( C > 1) hierarchically by generating C-1 candidate C -cluster models using the split-cluster E-M procedure shown in the blue box in Supplementary Figure 15 , and described subsequently. From these C-1 candidate C- cluster models, the model with the highest likelihood is chosen.

The split-cluster E-M procedure shown in the blue box in Supplementary Figure 15 is used to generate each of the C-1 candidate C cluster models. For each of the C-1 clusters, the split-cluster E-M procedure first finds the optimal split of cluster c into two clusters. To find the optimal split of cluster c into two clusters, the data points belonging to cluster c are randomly assigned to two separate clusters. The optimal model parameters for these two data subsets are then obtained using the single-cluster E-M procedure (greenpersName box in Supplementary Figure 15 ). These cluster parameters are used to initialise the fitting of a two-cluster model to the subset of the data belonging to cluster c , using E-M. This two-cluster solution is then used together with the other C-2 clusters to initialise the fitting of the C -cluster model. The C -cluster model is then optimised using E-M, alternating between updating the sequences S c for each cluster and the fractions f c . This procedure is repeated for 25 different start points (random cluster assignments) to find the maximum likelihood solution (see Convergence).

The single-cluster E-M procedure shown in the greenpersName box in Supplementary Figure 15 is used to find the optimal model parameters (the sequence S in which the biomarkers reach each z -score) for a single-cluster. In the single-cluster E-M procedure the sequence S is initialised randomly. This sequence is then optimised using E-M by going through each z -score event E in turn and finding its optimal position in the sequence relative to the other z -score events, i.e. by fixing the order of the subsequence T = S /E and maximising the likelihood of the sequence by changing the position of event e in the subsequence T . The sequence S is updated until convergence. Again the single-cluster sequence S is optimised from 25 different random starting sequences to find the maximum likelihood solution (see Convergence).

Convergence At several points in the model fitting we perform a greedy optimisation from a number of different starting points and choose the maximum likelihood sequence or set of sequences. The multiple runs safeguard against local minima. However, in fact, we find that the optimisation displays good convergence: runs from all start points typically converge to a solution with likelihood within a 1 × 10 −3 % of the maximum likelihood, and within the uncertainty estimated by the uncertainty estimation procedure (see Uncertainty estimation). The convergence of the SuStaInterm·mth algorithm and ability to locate the global minimum and correct solution is further demonstrated in simulation using synthetic data in (Supplementary Results: Convergence).

Uncertainty estimation In addition to estimating the most probable sequence S c for each subtype, we can determine the relative likelihood of all sequences for each subtype by evaluating the probability of each possible sequence. This gives us an estimate of the uncertainty in the ordering S c , which we summarise by plotting the probability that each z -score event appears at each position in the sequence for each subtype. We visualise this probability (see Fig. 2 for example) using different colours to indicate the cumulative probability each region has reached a particular z -score: the cumulative probability of a region going from a z -score of 0-sigma to 1-sigma ranges from 0 in white to 1 in red, the cumulative probability of a region going from a z -score of 1-sigma to 2-sigma ranges from 0 in red to 1 in magenta, and the cumulative probability of a region going from a z -score of 2-sigma to 3-sigma ranges from 0 in magenta to 1 in blue. In practise the number of sequences is too large to evaluate all possible sequences so we use Markov Chain Monte Carlo (MCMC) sampling to provide an approximation to this uncertainty, as in 7 , 8 . As in refs. 7 , 8 , we take 1,000,000 MCMC samples initialised from the maximum likelihood solution, checking that the MCMC trace shows good mixing properties.

Cross-validation We use ten-fold cross-validation here for two distinct purposes: (i) to evaluate the optimal number of subtypes and (ii) to evaluate the consistency of the subtype progression patterns. We evaluated the optimal number of subtypes using the Cross-Validation Information Criterion (CVIC) 44 , i.e. by evaluating the likelihood of each c -subtype model from c = 1 … C on the test data for each fold and choosing the model with the highest out-of-sample likelihood P ( X | M ), or equivalently the lowest value of the CVIC, across all folds. The CVIC is defined as CVIC = −2 × log( P ( X | M )), where P ( X | M ) is the probability of the data for a particular SuStaInterm·mth model, M , i.e. P ( X | M ) = ∑ c = 1 C P X ∣ S c P S c . In cases where the evidence for a more complex model was not strong (a difference of less than 6 between the CVIC and the minimum CVIC across models, or equivalently a difference of less than 3 between the out-of-sample log-likelihood and the minimum out-of-sample log-likelihood across models), we favoured the less complex model to avoid over-fitting 45 . To evaluate the consistency of the subtype progression patterns we performed ten-fold cross-validation by dividing the data into ten folds and re-fitting the model to each subset of the data, with one of the folds retained for testing each time. We report the consistency of the models across folds by computing the similarity between the progression patterns of two subtypes (see Similarity between two subtype progression patterns): the model fitted to each fold and the model fitted to the whole data set.

Similarity between two subtype progression patterns To enable the comparison of subtype progression patterns in data subsets (Supplementary Figure 13B ) and across cross-validation folds (CVS in Figs. 2 – 4 , and Supplementary Figures 13A and 14 ), we measure the similarity of pairs of subtype progression patterns using the Bhattacharyya coefficient 46 . We evaluate the Bhattacharyya coefficient between the position of each biomarker event in the two subtype progression patterns, averaged across biomarker events and MCMC samples. The Bhattacharyya coefficient measures the similarity of the distribution of the position of biomarker events in the subtype sequences and ranges from 0 (maximum dissimilarity) to 1 (maximum similarity).

Patient subtyping and staging We assigned subjects to subtypes and stages predicted by the SuStaInterm·mth model (Fig. 5 , Tables 1 and 3 ) by first evaluating the likelihood that they belonged to each subtype (by integrating over disease stage) and choosing the subtype with the highest likelihood, and then evaluating the probability they belonged to each stage of the most probable subtype and choosing the stage with the highest likelihood. When evaluating the likelihood we integrated over the set of MCMC samples to account for the uncertainty in the model parameters, rather than just evaluating the likelihood at the maximum likelihood parameters. This means that a patient’s model stage indicates the average position over the posterior distribution on the sequence given the data.

Strength of assignment to subtype We evaluated the strength of an individual’s assignment to a particular subtype by comparing the probability that they were at stage ≤2 (i.e. they had no major imaging abnormalities and therefore could not be assigned to a particular subtype) with the probability that they belonged to each SuStaInterm·mth subtype (probability for each subtype summed over stages 3+). The strength of the assignment was evaluated as their maximum probability of belonging to one of the subtypes. We considered those with a maximum probability of belonging to a particular subtype of greater than a half as having a strong assignment to a subtype.

Comparison to subtypes-only and stages-only modelsterm·mth We compared our SuStaInterm·mth model to a subtypes-only modelterm·mth and a stages-only modelterm·mth. In the subtypes-only modelterm·mth, individuals are clustered together into groups based on the similarity of their biomarker measurements—without accounting for heterogeneity in disease stage. The stages-only modelterm·mth is a disease progression model where all subjects are assumed to be samples of a single common progression pattern—without accounting for heterogeneity in disease subtype. We formulated the subtypes-only and stages-only modelsterm·mth so that they were as close as possible to the SuStaInterm·mth model, but did not model heterogeneity in disease stage or disease subtype, respectively. This allows us to assess the benefit of accounting for this disease stage or subtype heterogeneity in the SuStaInterm·mth model. The subtypes-only modelterm·mth consists of a mixture of Gaussians with unknown mean and variance. The subtypes-only modelterm·mth is fitted to symptomatic mutation carriers for GENFI, and AD subjects for ADNI, so that the subtypes correspond to a single diagnostic group. As done for the SuStaInterm·mth model, we evaluated the optimal number of clusters (subtypes) using the CVIC 44 . The stages-only modelterm·mth is a special case of the SuStaInterm·mth model outlined in Mathematical Model, where only a single subtype is modelled, i.e. C = 1.

SuStaInterm·mth modelling of GENFI data set Supplementary Tables 4 – 6 provide a summary of the settings of the SuStaInterm·mth algorithm. We applied SuStaInterm·mth modelling to various subgroups of the GENFI data set: all 172 mutation carriers, 76 GRN mutation carriers, 63 C9orf72 mutation carriers, 33 MAPT mutation carriers. For all mutation carriers we fitted SuStaInterm·mth models of up to a maximum of 5 subtypes. For the GRN, C9orf72 and MAPT mutation carriers we fitted SuStaInterm·mth models of up to a maximum of 3 subtypes. We chose the z -score events for the GENFI data set to include z -scores of 1, 2 and 3 for each volume, but excluded z -score events where fewer than 10 mutation carriers had values that were greater than that z -score. The maximum z -score, which is reached at the final stage of the progression, was set to be 2, 3 or 5 depending on whether the maximum z -score event was 1, 2 or 3, respectively. We maintained the same z -score events across each of the GENFI experiments.

SuStaInterm·mth modelling of ADNI data set We applied SuStaInterm·mth modelling to two largely independent (overlap of 59 individuals) subgroups of the ADNI data set: 793 individuals with 3T MRI scans and 576 individuals with 1.5T MRI scans, for which we tested SuStaInterm·mth models of up to a maximum of 5 subtypes. As we did for GENFI, we chose the z -score events to include z -scores of 1, 2 or 3 for each volume, but excluded z -score events where fewer than 10 subjects had values that were greater than that z -score. Again the maximum z -score, which is reached at the final stage of the progression, was set to be 2, 3 or 5 depending on whether the maximum z -score event was 1, 2 or 3 respectively. Full details of the settings of the SuStaInterm·mth algorithm can be found in Supplementary Tables 4 – 6 .

Classification of mutation groups using subtypes We performed two experiments to compare the ability of subtypes obtained from SuStaInterm·mth and the subtypes-only modelterm·mth to classify mutation carriers in GENFI into their different mutation groups. In the first experiment (Tables 1 and 2 ) we optimised the probability required for assignment to each of the subtypes. This accounts for different amounts of heterogeneity within the different subtypes. In the second experiment (Supplementary Tables 1 and 2 ) we simply assigned individuals to their most probable subtype and compared their assigned subtype with their mutation group. In both experiments the classification results are reported as out-of-sample accuracies obtained through 10-fold cross-validation.

Electronic supplementary material Supplementary Information Peer Review File