Resting State Network Estimation in Individual Subjects

2.1 Participants

Perceptron training, optimization and validation used data sets previously acquired at the Neuroimaging Laboratories (NIL) at the Washington University School of Medicine. A second, large validation data set was obtained from the Harvard-MGH Brain Genomics Superstruct Project (Yeo et al., 2011). All patients were young adults screened to exclude neurological impairment and psychotropic medications. Demographic information and acquisition parameters are given in Table 1.

2.2 Neuroimaging methods

2.3 Meta-analysis of task fMRI and generation of training data

ROIs representing distinct RSNs were isolated by meta-analysis of task-fMRI responses. We initially targeted 10 functional systems, each represented by a variable number of response foci derived from the literature (Table S1). The initial set of response foci was refined to ensure that all ROIs assigned to the same RSN generated maximally similar correlation maps and that ROIs assigned to different RSNs generated distinct correlation maps. These criteria yielded 169 ROIs representing 7 RSNs with high intra- and low inter-network correlation (Figures 1 and 2 and Table 2; see Supplemental text section 2.3 for algorithmic details). Thus, 3 of the original 10 networks were subsumed into the remaining networks. To these were added a nuisance category consisting of 6 ROIs in cerebrospinal fluid (CSF) spaces. The latter enabled the MLP to separate correlation patterns representing CSF vs. true RSNs. Computing correlation maps for each of the 175 seed regions in all 21 Training subjects produced 3,675 images used as training data. Each image in the Training set was masked to include only grey matter voxels, producing a 3,675 × 30,981 matrix. Similarly, 175 regions 17 subjects = 2,975 images, 175 regions · 10 subjects = 1,750 images and 175 regions · 692 subjects = 121,100 images, were computed in the Optimization, Validation 1, and Validation 2 data sets, respectively.

2.4 MLP design

The core of a perceptron is an artificial neural network consisting of an input, hidden, and output layer, each consisting of nodes fully connecting to the next layer (all-to-all feed-forward). Training samples are passed into this feed-forward network and the output is compared to the a priori assigned RSN label. The error in this comparison is used to update the connection weights, between layers to increase the performance of the MLP (Rumelhart et al., 1986).

2.5 Application of method to individuals

To map RSNs in individual subjects, a correlation map was generated for every voxel in the brain and then propagated through the optimized perceptron. An overall schematic of this process is depicted in Figure 2.

2.6 Comparison of the MLP to linear discriminant analysis and dual regression

Linear discriminant analysis (LDA) is a classification algorithm that operates by projecting data onto vectors that maximize the ratio of between-class scatter to within-class scatter (Appendix B). PCA preprocessing is essential in LDA because the dimensionality of the input space (originally, tens of thousands of voxels) must be substantially less than the number of training examples (Belhumeur et al., 1997). By systematic exploration of the feasible PCA dimensionality range, it was determined that the lowest classification error in the Optimization set was obtained with 20 principal components. All present LDA results were obtained with this LDA design.

Dual regression (DR) is a technique commonly used subsequent to group ICA to extend spatial ICs discovered at the group level to individuals. Group ICA is a technique for discovering RSNs whereas the present work is about mapping known RSNs, not about RSN discovery (see Introduction). In fact, the present LAN (see section 2.3) typically is not obtained by group ICA, e.g., (Damoiseaux et al., 2006). Nevertheless, DR is algebraically well defined regardless of the origin of the group-level RSNs (Appendix B). Here, we used RSNs derived by meta-analysis of task-fMRI. Seed-based correlation maps were averaged across all participants and across all seed regions representing each class, generating 7 spatial RSN components and one nuisance component. Dual regression was then performed as described in Appendix B (cf. (Zuo et al., 2010)).

Voxel-wise RSN estimates were generated using each technique (MLP, DR, LDA) and converted to a percentile scale as described in section 2.5.1 for visualization. To quantitatively compare performance across methods, RSN estimates were computed for maps generated using the a priori ROIs. Seed-based covariance maps were used to evaluate DR performance (Eq. (B.5) in Appendix B); seed-based correlation maps were used to evaluate LDA and MLP performance (Eqs. (B.6) and (B.7)). These evaluations yielded RSN estimates for each seed for each participant for each method. From these values and the associated a priori RSN labels, ROC curves and AUC scores were computed for each network (see section 2.4.3).

3.1 MLP training

3.2 MLP performance in the Validation datasets

3.3 Comparison of the MLP to alternative RSN estimation schemes

Dual regression (DR) and linear discriminant analysis (LDA) produced RSN topographies that were broadly similar to MLP results at the group level (Fig. 9A). However, DR and LDA topographies generally had greater spatial extent and more overlap across networks than the MLP-derived topographies. For example, the DAN network exhibited strong (greater than 90th percentile) values in dorsal SMN regions for both DR and LDA, but not for the MLP method. Similarly, VIS topography extended to SMN and DAN regions in DR results, less so in LDA results, and minimally in MLP results.

A cross-methodological comparison of classification performance in terms of AUC is given in Table 4. RSN scores were computed as described in section 2.6. These RSN scores were converted to hard classifications by identifying the group-derived map yielding the greatest score for each seed ROI. Linear Projection (LP), the simplest possible method, is included in Table 4 for comparison with the other methods (see Linear Projection in Appendix B). LP performance was evaluated using seed-based covariance maps (Eq. (B.1) in Appendix B). This is essentially nearest neighbor classification. That is, LP AUC is approximately equivalent to the probability that a map in, e.g., in PCA space is closest to the group average RSN template (represented by cluster centers in Figure 3C). Table 4 reports RMS error only for LDA and MLP results because the outputs generated by Linear Projection and DR are not confined to the interval [0,1] and cannot be directly compared to a priori labels.

To illustrate the statistical covariance of RSN estimates, ROIs were plotted in planes defined by pairs of RSN scores (Figure 9B, top row). SMN ROIs (cyan cluster) achieved the highest score along the SMN axis in the SMN vs. VIS plane, and likewise for VIS ROIs. In plots for DR and LDA, DAN ROIs achieved relatively high scores in both dimensions (dark blue cluster in quadrant 1). This observation corresponds to the moderate spatial overlap between DAN, VIS, and SMN in DR results in Fig. 9A. SMN (cyan) and VIS (green) vectors extend from the median score to the center of mass of the ROI clusters. These vectors form an acute angle in the DR results, corresponding to generally correlated SMN vs. VIS estimates of network identity. LDA separates these vectors by a larger angle, which corresponds to less VIS-SMN spatial overlap. However, the DAN ROIs still achieve a high score on the VIS and SMN axes; the significant overlap of the DAN and SMN clusters corresponds to the spatial overlap of the DAN and SMN topographies in Fig. 9A.

In contrast, MLP showed greater within-network scores and lower across-network scores. For example, SMN and VIS ROIs achieved values near unity along their respective axes, while ROIs from other networks, notably DAN, achieved values closer to the median score. MLP output scores across anti-correlated networks also were closer to the median score. As might be expected, the contrast between MLP vs. either LDA or DR was especially marked in the case of difficult to separate RSNs, specifically, LAN vs. VAN and LAN vs. DMN.

Estimates of RSN identity were generally more correlated for DR and LDA methods than MLP. The Pearson correlation coefficient was computed (across all ROIs) for every pair of RSN classes, to generate inter-class correlation matrices in Figure 9C. The DR-derived inter-class correlation (Fig. 9C) retained more of the structure evident in Fig. 3B than the MLP-derived result. The magnitude of inter-class correlation was smaller for the MLP method than for either DR or LDA (p<0.05 and p<10⁻⁵, respectively, Mann-Whitney U-test), indicating that the MLP produced more orthogonal estimates of RSN membership than the other methods. Overall, classification performance was higher for the MLP (0.982 AUC) than for DR (0.954) and LDA (0.970).

4.1 Inter-individual variability of MLP outputs

The present results (Figures 5–7) exhibit a high degree of face validity with respect to the training data and previously reported RSN results. Thus, for example, components of the DMN used as seeds to generate the training data were classified as DMN in all participants. This was true not only for easily classified networks (e.g., the DMN) but also for networks (e.g., VAN and LAN) that are inconsistently found by unsupervised procedures. The results shown in Figure 5 illustrate that the perceptron reliably classified RSNs in each individual in the Validation 1 group (0.984 worst case AUC), even in cases in which the RMS error was relatively high (> 0.175).

Analysis of voxelwise RSN membership estimates across a large cohort (100 subjects from Validation dataset 2) revealed RSN-specific zones of high as well as low inter-individual variability (Fig. 7, middle row). Voxels with high RSN scores generally showed the least inter-subject variability. Such regions, e.g., the posterior parietal component of the DMN (Fig. 7, right column), were surrounded by zones of high variability (e.g., a ring around the angular gyrus). The pre-and post-central gyri consistently showed high SMN RSN scores but were bordered by regions of high inter-subject SMN variability. Interestingly, inter-subject variability was low also in areas with RSN scores near 0, particularly in areas typically anticorrelated with other networks (e.g., low DAN variance in the angular gyrus, a component of the DMN; low DMN variance in MT+, a component of the DAN) (Fox et al., 2006).

At least four factors potentially contribute to observed inter-subject MLP output variability: (i) limited or compromised fMRI data; (ii) limitations intrinsic to the MLP; (iii) differences in RSN topography attributable to anatomic variability; (iv) true differences in RSN topography independent of variable gyral anatomy. We consider each of the possibilities in turn.

(i) The fMRI data used in the present work were obtained in healthy, cooperative young adults. Hence, the fraction of frames excluded because of head motion (Power et al., 2012) was low (about 4%). The total quantity of fMRI data acquired in each individual was generous by current standards (Van Dijk et al., 2010). However, fMRI data quantity clearly affects MLP performance (see section 4.5.2 below and Figure 10A). Current results suggest that more data generally improves MLP performance. The requirements for fMRI data quality and quantity for acceptable MLP performance in clinical applications remains to be determined. (ii) High variability in RSN boundary regions (e.g., Fig. 7, DMN, middle row) may reflect uncertainty attributable to multiple RSN membership, i.e., voxels with high “participation coefficients” (Guimera and Amaral, 2005; Power et al., 2011). Voxels with multiple RSN membership may be more difficult to classify because the training data included only maps derived from seeds assigned to single networks.

(iii) Figure 6 indicates that the MLP captures a substantial portion of anatomic variability. However, some part of MLP mapping imprecision may be explained by uncorrected anatomical variability. To investigate this possibility, we compared the overall RSN standard deviation map (Fig. S3A) to sulcal depth variability (Fig. S3B) and found a weak spatial correlation (r = 0.2). By inspection, these maps were concordant only at a broad spatial scale: both showed low variability in primary motor/auditory/insular cortices and high variability elsewhere. Little correspondence was evident at finer scales (note lack of annular patterns in Fig. S3B). The degree to which anatomical variability contributes to spurious variance in RSN topography estimates may be addressed by measuring the degree to which non-linear or surface-based registration decreases inter-subject variance and increases overall MLP performance (higher AUC, lower RMS error).

(iv) On the other hand, inter-individual differences may reflect “true” individual variability in RSN topography independent of gyral anatomy (Mueller et al., 2013). Previous work has demonstrated that inter-individual differences in task-evoked activity correspond to “transition zones” in resting state networks (e.g., the boundary between parietal DMN and DAN regions) (Mennes et al., 2010). These same regions appear in our inter-subject variance maps for both DMN and DAN (Fig. 7). We also note that areas of high RSN score variability (pre-frontal, parietal, lateral temporal) broadly correspond to regions exhibiting the greatest expansion over the course of human development and evolution (Hill et al., 2010). This correspondence may be coincidental, but it is consistent with the hypothesis that later developing or evolutionarily more recent areas of the brain tend to be more variable across individuals.

In summary, there are many potential contributions to observed RSN topographic variability. Because we did not use non-linear volume registration in this work, Fig. 6 retains individual differences in gyral anatomy. By inspection, the MLP was able to track these differences. Thus, it is reasonable to expect that the MLP will find “true” differences in RSN topography not attributable to variable gyral anatomy. Future studies are needed to compare MLP-derived topographies with, e.g., task-evoked responses, after correcting for anatomical variability.

4.2 External Validity and Generalizability

Two distinct types of external validity, that is, correct classification outside the Training set (areas covered by the seed ROIs), are evident in our results. First, high overall MLP performance was achieved for a priori seed-based correlation maps in the Optimization (98.2% AUC) and Validation 1 datasets (98.8% AUC). Performance was reliable in all participants (97.1% worst-case AUC), which is critical in clinical applications. Second, and perhaps of greater scientific interest, the RSN estimates in areas not covered by seed regions were strongly concordant with previously reported task-based and resting-state fMRI results. For example, while no temporal FPC seed ROI was included in the training set, a posterior temporal gyrus locus was classified as FPC (Lee et al., 2012; Power et al., 2011; Yeo et al., 2011) at the group level (Fig. 7, FPC column). Similarly, the MLP also identified the parahippocampal gyrus as DMN (Kahn et al., 2008).

The MLP identified as LAN components of dorsal and ventral streams previously associated with cortical speech processing (Hickok and Poeppel, 2004) see also (Binder et al., 2011). Whereas the ventral components were bilateral, the dorsal components were left lateralized (Fig. 8, white arrow, LAN column) (Hickok and Poeppel, 2007). This finding represents a strong demonstration of external validity, as the left dorsal region was not included in the Training set. The right inferior cerebellum was first associated with language function by PET studies of semantic association tasks (Petersen et al., 1988). Identification of this region here as part of the LAN network (Fig. 8, WTA, Z = −30 and −47) is doubly significant: First, no cerebellar seeds were used to generate training data and cerebellar voxels were excluded from the gray matter mask; hence, these voxels were not seen by the MLP. Second, lateralized cerebellar RSN components typically are not found by unsupervised seed-based correlation mapping (Buckner et al., 2011).

These findings highlight the capabilities of supervised learning applied to the problem of identifying RSNs in individuals. The cortical representation of language (primarily Broca's and Wernicke's areas) has been extensively studied using task-based fMRI (Binder et al., 2011) and correlation mapping with a priori selected ROIs (Briganti et al., 2012; Hampson et al., 2006; Pravata et al., 2011; Tomasi and Volkow, 2012). However, the language network, as presently defined, typically is not recovered as such by unsupervised methods, e.g., (Power et al., 2011; Yeo et al., 2011). Rather, components of the LAN are generally found only at fine-scale RSN descriptions. Thus, an RSN including Broca's and Wernicke's areas appears as the 11th of 23 components in (Doucet et al., 2011); these same areas were identified as VAN by Power and colleagues (2011) and DMN by Yeo and colleagues (2011). Lee et al., (2012), found a component consistent with the presently defined LAN at a hierarchical level of 11 (but not 7) clusters. Thus, the present work demonstrates the potential of supervised learning to find networks that are subtle features of the BOLD correlation structure. These may be minor sub-components within hierarchically organized RSNs or functional entities of high scientific interest or clinical value that do not fit within a hierarchical organization, i.e., extend over multiple levels or across multiple branches.

Because the MLP is a universal function approximator, it is subject to over-fitting. Over-fitting refers to learning features that are particular to the training set but that do not generalize to other sets. To minimize over-fitting, we halt training when Optimization dataset error reaches a minimum (Figure 4A). However, because early stopping is implemented throughout optimization by simulated annealing, over-fitting may still occur with respect to the Optimization dataset. Therefore, we demonstrate the performance of the fully trained MLP in two Validation datasets (Figs. 5–8).

4.3 Comparison of the MLP to dual regression (DR) and linear discriminant analysis (LDA)

The MLP, DR, and LDA all represent different strategies for extending RSN mapping results obtained at the group level to individuals. The MLP showed better spatial specificity of RSN estimates than either DR or LDA (Fig. 9A) and produced more orthogonal estimates of RSN identity in a priori ROIs (Fig. 9B). These comparative results indicate greater statistical independence of RSN estimates obtained by the MLP. However, perfect separation of classes was not achieved by any method, as indicated by residual correlation of RSN estimates (non-zero off-diagonal elements in Fig. 9C). This residual was greatest in the least separable RSNs (compare Figs. 9C and 2B).

Performance differences across the three methods can be related to their underlying algebraic structures, which reveal interesting commonalities as well as differences (Appendix B). All three methods operate by transforming individual voxel-wise covariance or correlation matrices into 7 dimensional RSN membership estimates at each voxel. However, the methods differ in the strategies used to achieve separation of RSN estimates. LDA requires considerable dimensionality reduction by PCA preprocessing. The optimal number of PCs was found to be 20, but this accounted for only 70% of the variance. In contrast, the MLP does not strictly require PCA preprocessing although this step greatly reduces the computational load without sacrificing information. In fact, the MLP performance optimum was obtained with 2,500 PCs (Figure S2), which included 99.97% of the variance. DR, as conventionally implemented (Zuo et al., 2010), does not involve PCA preprocessing or training.

Both LDA and MLP optimize separation of classes using thousands of labeled training samples, which captures variability across brain regions and individuals. Therefore, the superior performance of the MLP is not simply attributable to a large number of training samples. Rather, high-dimensional non-linear classification boundaries allow the MLP to extract arbitrary features from a large (2,500 PCs) input space. Future work will compare the MLP to LDA and DR in patients with focal lesions, i.e., stroke or brain tumors. At issue is performance under circumstances in which RSNs exhibit altered topography.

4.4 Dynamics of perceptron training reflect hierarchical brain organization

Several studies have demonstrated that resting state networks are hierarchically organized (Boly et al., 2012; Cordes et al., 2002; Doucet et al., 2011; Lee et al., 2012; Marrelec et al., 2008). Here, the hierarchical scale of an RSN is reflected in training performance trajectories (Fig. 4B): the DMN was the first to be separated from other RSNs. The DMN arguably is the most robust feature in the correlation structure of intrinsic brain activity. Its topography is very similar across RSN mapping strategies, specifically, spatial ICA (Beckmann et al., 2005) and seed-based correlation mapping (Yeo et al., 2011). Here, the DMN and regions anti-correlated with the DMN (Fox et al., 2005) were well separated along the first principal component of the training data (Figure 3)¹.

After the DMN, the sensorimotor and visual networks were next to achieve separation (Fig. 4B). These networks are often seen at the next level down in the RSN hierarchy as offshoots of the anti-DMN (Lee et al., 2012) or `extrinsic system' (Doucet et al., 2011). The DAN achieved only a small peak in error descent compared to other `extrinsic' networks, although this occurred in close proximity to SMN and VIS. In contrast, the LAN and VAN were last to achieve separation. This corresponds to the observation that LAN and VAN systems are typically found by analyses extending to lower levels of the RSN hierarchy (Lee et al., 2012; Power et al., 2011).

4.5 MLP performance applied to the evaluation of fMRI methodology

4.6 Limitations

The present MLP results were obtained using a particular set of RSNs selected for their scientific and clinical value in the context of ongoing research. A different training set including other RSNs might be optimal in other circumstances. MLP performance may be expected to vary with the number and separability of RSNs being classified as well as the total duration (Fig. 10A) and quality of acquired resting state fMRI data.

As discussed above (section 4.1), inter-individual differences in computed RSN topographies may reflect multiple factors. Cross-gyral contamination due to the relatively large voxels used in this study (3–4 mm acquisition, 3 mm post-processing analyses) may limit the precision of RSN estimation in our dataset (see (Yeo et al., 2011) for a discussion of this issue). Of greater theoretical importance is that we cannot distinguish between correctly identified geodesic (parallel to the cortical surface) shifts of RSN topography versus morphological differences in gyral/sulcal anatomy. Potential strategies for validating MLP-derived results individuals include comparison with measures of structural connectivity (Damoiseaux and Greicius, 2009) and invasive electrophysiologic recording (He et al., 2008).

As shown in Figure 3B as well as in prior work (Boly et al., 2012; Doucet et al., 2011; Lee et al., 2012; Marrelec et al., 2008), RSNs are hierarchically organized. Thus, a major dichotomy separates the `extrinsic' (DAN, VAN, VIS, SMN) from `intrinsic' networks (FPC, LAN, DMN). At the next level of granularity, RSN blocks are distinguishable within the major dichotomy. However, the cognitive domains corresponding to these RSNs (e.g., `motor', `language', `attention', etc.) are not conventionally regarded as hierarchically organized. Our training labels are categorical because they were generated from task-fMRI responses corresponding to conventionally understood cognitive domains (Table 2). Thus, this work is fundamentally limited in that we impose a `flat' conceptualization of brain function onto an intrinsically hierarchical system. This is reflected in the fact that ROIs are almost perfectly classified (by the AUC measure) in all datasets, yet RMS error always remains substantial (15%–20%) (Table 3, Fig. 10A). The existence of this residual may indicate that resting-state brain networks are inherently non-separable in the sense of classification. Indeed, this is consistent with the notion of “near decomposability” of hierarchical systems formed by multiple, sparsely inter-connected modules (Simon, 1995); this concept has been extended to brain networks (Meunier et al., 2010).

Another potential limitation is that all present analyses assume temporal stationarity, as do the overwhelming majority of extant papers on intrinsic BOLD activity. However, a growing literature (Chang and Glover, 2011; Kiviniemi et al 2011; Hutchinson et al 2012; Allen et al 2012) indicates that resting state correlation patterns fluctuate on a time scale of minutes. Our analysis of the effect on classifier performance of varying the total quantity of fMRI data (Fig. 10A) utilized contiguous epochs to the extent possible given fMRI runs of duration ~5 minutes. Less data reliably yielded worse performance. This result may be attributable either to sampling error (i.e., limited signal to noise ratio) or true temporal non-stationarities or both. Similarly, the performance floor may reflect non-stationarities. Disambiguating these possibilities would require much longer contiguous acquisitions.

4.7 Summary

The MLP estimates RSN membership at the voxel level via computed correlation maps. After training, it reliably identifies RSN topographies in individuals. RSN estimation is rapid (2 minutes using Matlab running on Intel i7 processors) and automated, hence suitable for deployment in clinical environments. After training, operation is independent of any particular seed. Therefore, the trained MLP is expected to be robust to anatomical shifts and distortions, for example, owing to enlarged ventricles and mass effects or even loss of neural tissue (e.g., stroke).

In this work, the MLP was trained to operate in 3D image space for compatibility with clinical imaging formats. However, the MLP concept can be readily adapted to operate on correlation maps represented on the cortical surface. Similarly, an MLP can be trained to ignore anatomical abnormalities (e.g., brain tumors) by altering the domain of the training set, i.e., excluding tumor voxels. This possibility will be explored in future work.