Here, we use TDA to explore the landscape of brain imaging phenotypes of very young boys with FXS. In contrast to the focus of Hoeft et al. [ 2011 ], which used brain images to discriminate between iAUT and FXS, the specific goal of this project was to identify previously unobserved, yet neurobiologically salient subgroups within FXS that exhibit separable neuroanatomical phenotypes. The output of TDA, taking the form of subgroups of children who share similar brain structure patterns, can then be used to identify brain regions that characterize subgroup differences at the anatomical level as well as possible differences in behavioral profiles. Our primary hypothesis was that application of TDA would uncover brain‐imaging phenotypes that were not just anatomically separable, but also different from a clinical viewpoint.

To investigate putative phenotypic subgroups within FXS as captured by MRI data, we employed topological data analysis (TDA) [Carlsson, 2009 ], a recently‐developed approach that is specifically designed to identify structural characteristics of high‐dimensional datasets. One of the strengths of TDA is its ability to reduce such high‐dimensional data to human‐readable representations that capture essential features, similar to the manner in which a topographical map is able to capture the essential features of a landscape. An example of this approach can be found in Nicolau et al. [ 2011 ], where a TDA‐based approach was successful in elucidating a previously unidentified subgroup of breast cancers that exhibit 100% survival and no metastases. To our knowledge, our study represents the first application of TDA to appear in the literature on neuroimaging.

Because of similarities between the behavioral profiles of children with the single gene disorder fragile X syndrome (FXS) and the criteria used to define autism spectrum disorders (ASD), it had been hoped that FXS might serve as a useful genetically‐defined model for studying ASD. A recent study by our group [Hoeft et al., 2011 ]; however, showed that the gray and white matter profiles of young children with FXS are significantly different from those of children with idiopathic autism (iAUT; i.e., who do not have FXS). These differences were of sufficient magnitude that the two populations could be discriminated with a high degree of accuracy (90%) through the use of machine learning approaches. This finding also supports the hypothesis that there is a high level of neurobiological heterogeneity among individuals meeting diagnostic criteria for iASD [Abrahams and Geschwind, 2008 ]. In the study presented here, we concentrate exclusively on children with FXS in order to explore the degree to which neurobiological heterogeneity may be present within the FXS population itself. To date, there has been little work suggesting the existence of separable neurobiological phenotypes within FXS, the notable exception being a study by Jacquemont et al. [ 2011 ] that suggests methylation status may constitute a biomarker for predicting response to AFQ056, a subtype‐selective mGluR5 inhibitor. Our motivation to search for neuro‐phenotypic subgroups within FXS derived from previous studies suggesting the existence of behavioral subgroups based on whether individuals met criteria for autism [Brock and Hatton, 2010 ; Wolff et al., 2012 ]. Thus, we sought to explore the neuro‐phenotypic “landscape” of FXS to attempt to determine whether there may be previously undiscovered biological underpinnings to this behavioral observation. In this sense, our work may be seen as part of current efforts to understand the biological basis of neuropsychiatric pathology (e.g., the RDoC approach articulated by NIMH: http://www.nimh.nih.gov/research‐priorities/rdoc/index. shtml .)

MATERIALS AND METHODS

Participants and Data Collection The data utilized in our study were collected from infant and toddler boys who participated in a study of brain development in FXS [Hoeft et al., 2011]. Participants in this study included boys with FXS, idiopathic autism (iAUT), and developmental delays (DD) as well as typically developing (TD) children, all of whom were recruited by collaborating research teams at the Stanford University School of Medicine and the University of North Carolina, Chapel Hill; the current study focuses only on the subset of children from this earlier study that were diagnosed with FXS. The study protocols were approved by the human subjects committees at both institutions and consent from parents was obtained. Children with FXS (n = 52; mean [SD] age, 2.90 [0.63] years) were recruited through registry databases maintained by the Stanford University School of Medicine and the University of North Carolina, Chapel Hill, through postings to the National Fragile X Foundation Web site and quarterly newsletter, and through mailings to other regional FXS organizations. All participants had the “full mutation” form of the FMR1 gene known to cause FXS [Hoeft et al., 2008]. Participants completed the Autism Diagnostic Observation Schedule–Generic (ADOS) [Gotham et al., 2009; Lord et al., 2000] and their parents were given the Autism Diagnostic Interview (ADI)–Revised [Lord et al., 1994]. The Mullen Scales of Early Learning was administered to measure child IQ. There were no significant differences between sites in any of the cognitive variables for the participants with FXS (all P's > 0.05). [see Table 1 for a summary of demographic information, including the participant's level of Fragile X Mental Retardation Protein (FMRP)]. Anatomical MRI scan acquisition parameters consisted of a coronal T1‐weighted sequence (inversion recovery preparation pulse = 300 ms; repetition time (TR) = 12 ms; echo‐time (TE) = 5 ms; flip angle = 20°; slice thickness = 1.5 mm; number of excitations = 1; field‐of‐view (FOV) = 20 cm; matrix = 256 × 192) (Phantoms were used to ensure matching calibration of MRI scanners at both sites). Table 1. Demographic information from Hoeft et al. [ ] FXS iAUT DD TD Site, SU:UNC Particpants, no. 28:24 17:46 11:8 11:20 Age, yrs. Mean(SD) 2.90 (0.63) 2.77 (0.41) 2.96 (0.50) 2.55 (0.60) MSEL composite standardized score Mean (SD) 54.94 (9.14) 54.10 (9.41) 55.47 (7.53) 109.55 (17.24) FMRP (%) Mean (SD)a 5.83 (3.94) NA NA NA

Voxel‐Based Morphometry Preprocessing Standard voxel‐based morphometry (VBM) preprocessing of MRI data was carried out using the Statistical Parametric Mapping 5 (SPM5) statistical package (http://www.fil.ion.ucl.ac.uk/spm) and VBM5.1 (http://dbm.neuro.uni‐jena.de/vbm). Images were bias‐field corrected and segmented to GM, WM, and CSF. A Hidden Markov Random Field (HMRF; prior probability weight 0.3) was applied in order to incorporate spatial constraints arising from neighboring voxels. The images were normalized with a 12‐parameter affine transformation with a spatial frequency cut‐off of 25 in all three (x, y, z) directions and resampled to 1 × 1 × 1 mm3 voxels. Linear and nonlinear Jacobian modulation were applied. Customized GM, WM, and CSF templates created using all participants in our previous study [Hoeft et al., 2011] (FXS, n = 52; iAUT, n = 63; DD, n = 19; and TD, n = 31) were used for VBM preprocessing. For each participant, segmentation and normalization accuracy were manually inspected.

Multivariate Pattern Analysis of Magnetic Resonance Images Using Topological Data Analysis TDA is one of a general class of approaches to analyzing high‐dimensional data known in the literature as multivariate pattern analysis (MVPA). Multivariate approaches are designed to detect effects that may be discernable within the relationships (patterns) among variables, but which may elude detection when variables are examined in isolation. One type of multivariate approach that has already been used extensively in brain imaging studies is the use of support vector machines (SVMs) [Bray et al., 2009]. TDA follows the initial steps taken in SVM analyses to prepare image data for analysis by building a data matrix from smoothed and vectorized individual images, but thereafter differs from SVM analyses in two critical ways (see below for a brief discussion of vectorization). One difference is that SVMs are used to compare conditions that are known a priori, such as disorder versus control or stimulus versus rest, and are thus said to be supervised approaches to MVPA. In contrast, TDA is unsupervised since, rather than compare predefined groups, it is used to identify coherent, but possibly heretofore unknown groups within the study population. Other examples of unsupervised approaches to MVPA include independent component analysis [Calhoun et al., 2012] and clustering via correlation matrices [Fair et al., 2012]. A second critical difference, and a key benefit of using TDA not found in other MVPA approaches, is that TDA can produce a compressed and easily readable visual representation of the data, called a Reeb graph, which preserves the underlying geometric structure of the data, and thus facilitates identification of its salient features [Carlsson, 2009]. As can be seen below in Figure 1, these features can encode not only information about clusters that may exist within the data, but also information about spectra (i.e., variation in the data due to continuously varying underlying parameters). To allow the reader to gain insight into the process of Reeb graph construction, we give a brief overview of the TDA pipeline as applied to a toy example from Lum and colleagues [Lum et al., 2013], and provide more a detailed description of our use of TDA in the following subsection. Figure 1 Open in figure viewer PowerPoint 2013 TDA pipeline. Visual depiction of TDA pipeline from point cloud to Reeb graph, reprinted from [Lum et al.,] in accordance with the Creative Commons license under which it was published (CC BY‐NC‐ND 3.0). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com .]

Overview of TDA As with SVMs, the starting point for TDA is a data matrix whose m rows and n columns correspond to the variables of the study and their observed values; in the typical TDA data matrix, each column corresponds to a variable (e.g., voxel), and each row corresponds to an observation (e.g., subject). In geometric terms, we may interpret the n values in a row as giving the coordinates of a point in n‐dimensional space, and can thus interpret the m rows as m points (collectively called a point cloud) in space. In this way, we may also interpret the point cloud in Figure 1A in the shape of a hand as being a visual representation of a data matrix with thousands of rows and three columns, where the columns correspond to the x, y, and z coordinates of the points in the point cloud. The first stage of TDA is to assign numerical values to each point in the point cloud, a process we refer to as selection of filter values; in Figure 1B, these values are represented by a color in the red‐blue spectrum. The next stage of TDA is to separate the point cloud into overlapping regions by separating the filter values into overlapping bins. For example, in Figure 1B we could view the filter values as indicating the distance from the tip of the middle finger as measured by a horizontal ruler, in which case we could then view the regions shown in Figure 1C as defined by the bins [0–1.5 in], [1–2.5 in], [2–3.5], [3–4.5], [4–5.5], [5–6.5]. The last stage of TDA is to represent each region from the previous stage by a dot, called a node, and then for each pair of regions that intersect, to connect the corresponding nodes by a line segment. Figure 1D shows the resulting Reeb graph, where the size and color of each node represent the number and “average” color of the points in the corresponding region. Note how the Reeb graph captures the underlying geometric structure of the hand, providing a schematic representation whose structural features correspond to real physical features, such as fingers, which also illustrates TDA's ability to highlight spectra and clusters within data, as described earlier. Note also that this construction does not involve multiple comparisons, and thus no correction for multiple comparisons is warranted. Unlike the case of the hand, where our familiarity with its shape allows us to validate the Reeb graph of the corresponding point cloud, the “shape” of the point cloud corresponding to MRI data for our study population of young boys with FXS is not known to us a priori. To address this issue, we use TDA to identify subgroups purely on the basis of the MRI data, but then apply standard statistical tests to compare clinical data for these subgroups to provide confirmation they are clinically, as well as anatomically distinct.

Details of Reeb Graph Construction As described earlier, we can interpret a subject's numerical data as the coordinates of a point in a high‐dimensional space S, referred to as subject space, and can thus recast the study population as a point cloud in subject space. In our case, we used the combined gray and white matter voxel data to create the subject space S, and then used the Iris software package to construct a Reeb graph of the data, so that we could explore the geometry of the resulting point cloud to identify its topological (i.e., shape) features [Lum et al., 2013] for an overview of Iris. These features would then correspond to anatomically defined subgroups of the study population, whose clinical profiles could then be compared and subjected to further analysis. In this subsection, we give a detailed description of the process used to construct a Reeb graph (Fig. 2 for a flowchart that summarizes this process). Figure 2 Open in figure viewer PowerPoint Reeb graph construction. Flowchart depicting the process used to construct a Reeb graph from raw MRI data. Note that the last four steps correspond to parts B, C, and D of Figure 1 , with the last two steps corresponding to part D. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com .] As noted above, Reeb graph construction begins, as do SVMs, with a data matrix. Since SVMs have already been used extensively in brain imaging studies, we followed the same approach to data matrix construction; namely, we combined the VBM‐preprocessed imaging data for the all subjects into a single matrix, with the data for each subject entered as a row of the data matrix. Specifically, we placed in each column the voxel intensity at a particular spatial location across all subjects. TThe process of converting a multidimensional object, like a 3D image, into an ordered list, like a row in a matrix, is known as vectorization, and is a common requirement of MVPA approaches [Pereira et al., 2009]. In principle, vectorization results in a loss of information about the spatial relationships between voxels, but in practice this loss does not appear to hamper the success of MVPA approaches in general, as evidenced by their use in numerous studies [Bray et al., 2009] for an extensive review of the use of MVPA in neuroimaging. Before the preprocessed data were vectorized, the intensities for each subject's gray matter 1 × 1 × 1 mm3 voxel images were binned and averaged to produce corresponding gray 4 × 4 × 4 mm3 voxel images, and the same was done for each subject's 1 mm white matter images. This step was undertaken to both reduce noise in the 1 mm images by smoothing and for ease of computation. The voxels in each of these new 4 mm 3D images were then linearly ordered, so that each subject's gray and white matter images were converted into gray and white matter intensity vectors, which were then concatenated to produce a single, combined gray and white matter intensity vector. These combined intensity vectors were then used as rows to build a data matrix M of dimension (# subjects) × (# gray voxels + # white voxels). M is obtained, construction of a Reeb graph proceeds in several steps: (1) column selection, (2) metric selection, (3) filter selection, and finally (4) selection of values of the filter parameter called gain and resolution. The first step, column selection, is used to define the subject space S; this is done by selecting a subset of columns of M, whose size then determines the dimension of S. In our case, the columns of M correspond to voxels, and we chose to define S by selecting those voxels whose variance was at least 0.03. (In modulated and segmented images such as those utilized here [Hoeft et al., 2011 S, and thus provide a measure of similarity between subjects. Euclidean distance, a generalization of the usual distance between points P 1 = (x 1 , y 1 ) and P 2 = (x 2 , y 2 ) in the coordinate plane is only one of many choices for a metric on S. Once the data matrixis obtained, construction of a Reeb graph proceeds in several steps: (1) column selection, (2) metric selection, (3) filter selection, and finally (4) selection of values of the filter parameter called gain and resolution. The first step, column selection, is used to define the subject space; this is done by selecting a subset of columns of, whose size then determines the dimension of. In our case, the columns ofcorrespond to voxels, and we chose to defineby selecting those voxels whose variance was at least 0.03. (In modulated and segmented images such as those utilized here [Hoeft et al.,], voxels intensities correspond to proportions of volume and thus range in value from 0 to 1). The next step, metric selection, is used to define distances within the subject space, and thus provide a measure of similarity between subjects. Euclidean distance, a generalization of the usual distance between points= () and= () in the coordinate plane is only one of many choices for a metric on Because we had no a priori basis for giving more weight in the analysis to a particular subset of the high‐variance voxels, we chose to standardize the voxel data before applying the Euclidean metric. The process of standardization, where data for each variable is first demeaned and then scaled to have variance 1, is a common statistical step taken in this context in order to give each variable equal weight in the subsequent analysis. In Iris, this is achieved by selecting the metric referred to as the variance‐normalized Euclidean metric. The third step, filter selection, is in many ways the most critical choice made in the construction of a Reeb graph since it is the filters that transform the similarity information determined in the previous two steps into a visual representation that is easily grasped by the human eye, namely, a 2D rendering of a 3D graph (Fig. 3, below). Filters provide an alternative notion of similarity among subjects that is distinct from the one defined by the metric chosen for the subject space S, and together these two notions of similarity guide the construction of the vertices of the Reeb graph, in two stages. In the first stage, the filters are used to sort the subjects into bins based on similarity, where bin size is controlled by the resolution parameter from step 4. To capture the structure of variability among the subjects, we chose the first two principal component scores (PC1 and PC2) of the variance‐normalized data as our filters, so that two subjects were judged to be similar (i.e., shared the same bin) if both their PC1 and PC2 scores were sufficiently similar. In the second stage, the metric selected in step 2 is used to further cluster the subjects within each bin, so that two subjects within the same bin are merged into the same cluster if the distance between them is sufficiently small; each cluster obtained in this way is then viewed as a vertex of the Reeb graph. Finally, the edges of the graph are constructed as follows: Although bin size is governed by the resolution parameter, bins are allowed to overlap (i.e., share subjects) to an extent governed by the gain parameter, and this overlap may lead clusters from different bins to share subjects. Any two clusters (i.e., vertices) that share a subject are then joined by an edge. We may think of the resulting graph as a view of the data through a microscope, where gain and resolution play roles analogous to level of light and level of magnification. A group of edge‐connected vertices of the resulting Reeb graph may then serve as a candidate subgroup within the data (Fig. 3). In our case, a choice of parameter values that clearly decomposed the data into subgroups consisted of a gain of 4 for both PC1 and PC2, a resolution of 45 for PC1, and a resolution of 30 for PC2. Figure 3 Open in figure viewer PowerPoint Iris rendering of a Reeb graph of the FXS data, with labels and ellipses added to indicate the subgroups LH and LL. Note that the size of each node corresponds to the number of subjects that were clustered to form that node, and that an edge between two nodes indicates that the corresponding clusters have a subject in common. The light blue node between the subgroups LH and LL is not included in either ellipse because it contains two subjects, one from each subgroup; nevertheless, the edges connected to this node indicate that each of the these subjects is also contained in a neighboring node within the appropriate subgroup. [Color figure can be viewed in the online issue, which is available at TDA results. Anrendering of a Reeb graph of the FXS data, with labels and ellipses added to indicate the subgroups LH and LL. Note that the size of each node corresponds to the number of subjects that were clustered to form that node, and that an edge between two nodes indicates that the corresponding clusters have a subject in common. The light blue node between the subgroups LH and LL is not included in either ellipse because it contains two subjects, one from each subgroup; nevertheless, the edges connected to this node indicate that each of the these subjects is also contained in a neighboring node within the appropriate subgroup. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com .]

Analysis of Behavioral Data Once subgroups were identified with TDA based on neuroimaging data, we used standard t‐tests to compare subgroup differences on behavioral measures.