Abstract

Resting state fMRI (R-fMRI) is a powerful in-vivo tool for examining the functional architecture of the human brain. Recent studies have demonstrated the ability to characterize transitions between functionally distinct cortical areas through the mapping of gradients in intrinsic functional connectivity (iFC) profiles. To date, this novel approach has primarily been applied to iFC profiles averaged across groups of individuals, or in one case, a single individual scanned multiple times. Here, we used a publically available R-fMRI dataset, in which 30 healthy participants were scanned 10 times (10 min per session), to investigate differences in full-brain transition profiles (i.e., gradient maps, edge maps) across individuals, and their reliability. 10-min R-fMRI scans were sufficient to achieve high accuracies in efforts to “fingerprint” individuals based upon full-brain transition profiles. Regarding test–retest reliability, the image-wise intraclass correlation coefficient (ICC) was moderate, and vertex-level ICC varied depending on region; larger durations of data yielded higher reliability scores universally. Initial application of gradient-based methodologies to a recently published dataset obtained from twins suggested inter-individual variation in areal profiles might have genetic and familial origins. Overall, these results illustrate the utility of gradient-based iFC approaches for studying inter-individual variation in brain function.

Introduction

The delineation of functional brain units, commonly referred to as cortical areas (Van Essen and Glasser 2014), is a fundamental challenge in neuroscience. Central to such efforts, is the determination of areas of transition from one cortical area to the next. Pioneering efforts have used a combination of histological, cytoarchitectural, and myeloarchitectural examinations (Brodmann 1909; Vogt and Vogt 1919; Felleman and Van Essen 1991; Zilles and Amunts 2010), along with lesion studies, to differentiate neighboring cortical territories with respect to their architectonic, connectivity, functional and topographic properties (Felleman and Van Essen 1991; Cohen et al. 2008). Converging evidence supports the notion that boundaries between areas are generally sharp rather than gradual (Kaas 1987; Amunts and Zilles 2012). A well-defined example is the primary visual area (V1), which is clearly different and sharply separated in cell bodies from the second visual cortex (V2) (Sincich and Horton 2005; Buckner and Yeo 2014). Given that these measurements are mainly based on non-human primates, postmortem examinations, or invasive recording experiments, the development of methods capable of mapping cortical areas in vivo remains an active area of research.

Cohen et al. (2008) demonstrated the ability to map transitions between nearby cortical areas through the detection of spatial variation (i.e., gradients) in intrinsic functional connectivity (iFC) profiles estimated from resting-state fMRI (R-fMRI). A number of successful applications have since emerged. At the local level, studies have delineated transition zones between cortical areas within key regions of interest (e.g., cingulate cortex, left lateral parietal cortex, and inferior frontal cortex) (Cohen et al. 2008; Nelson et al. 2010; Hirose et al., 2012, 2013). At the whole brain level, Wig et al. (2014) demonstrated the feasibility of delineating the broader range of boundaries between cortical areas in a single analysis based upon their gradient properties. Overall, these various studies have consistently recapitulated fine-grained cortical boundaries (e.g. between V1 and V2) previously established using histological and cytoarchitectural methodologies (Buckner and Yeo 2014). Additionally, recent studies demonstrated group-level transition-zone patterns to be highly reproducible across independent samples and studies (Gordon et al. 2014; Wig et al. 2014), thereby increasing enthusiasm towards the approach.

However, with few exceptions (e.g., Cohen et al. 2008; Wig et al. 2013; Laumann et al. 2015), studies focusing on detecting transition zones and cortical area boundaries have primarily relied on data averaged across dozens of individuals, or pooled individual-level maps to achieve more robust findings (e.g., Yeo et al. 2011; Kelly et al. 2012; Wig et al. 2014). Although effective in reducing noise, this is problematic, as a growing literature has suggested the presence of meaningful inter-individual variation in iFC patterns (Mueller et al. 2013), which appears to be stable over time and can be related to differences in behavior (e.g., Finn et al. 2015). Individual differences in areal organization are thought to be particularly prominent in higher-order association areas – a finding that is consistent with models suggesting the evolutionary value of cortical expansion (Mueller et al. 2013). Previous studies have suggested the potential importance of inter-individual variation in the transition zones or boundaries between cortical areas detected using iFC, by relating it to differences in task-evoked activation (Mennes et al. 2010) as well as variation in behaviorally quantifiable traits, such as social reciprocity and personality (Di Martino et al. 2009; Adelstein et al. 2011).

Recent work has underscored the feasibility of using iFC-based approaches to capture areal topographies at the individual level. First, Gordon et al. (2015) characterized the areal topological features at the individual-level, demonstrating that the areal sizes and the position of transition zones varied across individuals. Second, Laumann et al. (2015) mapped functional areal boundaries throughout the brain in a single human subject using a highly sampled dataset (>900 min). The individual-level architecture was broadly reflective of the universal architecture revealed in group-level analyses (Biswal et al. 2010; Wig et al. 2014), though with distinct variations that were fine-grained. Importantly, the areal organization was highly repeatable with subsets of data in this rich single subject dataset, suggesting that the generation of such highly sampled datasets is not a prerequisite for appreciating inter-individual variation.

Here, we provide a multifaceted evaluation of test–retest reliability for transition-zones properties revealed by iFC-based mapping approaches (i.e., the gradients measured, and the boundaries drawn between functional areas delineated using them). To accomplish this, we estimated functional gradients and boundaries (i.e., edge maps) for each of 300 publicly available R-fMRI datasets (duration = 10 min; standard EPI sequence; TR = 2000 ms), which were obtained from 30 participants who were scanned 10 times within one month. First we examined the “fingerprinting” potential of transition profiles (i.e., the ability to distinguish individuals from one another based upon their transition profiles). This was accomplished by testing whether the spatial correlation between any two sessions from the same individual was consistently greater than that observed between participants. Second, evaluated test–retest reliability for indices of areal organization at the full-brain and vertex levels, using image intraclass correlation coefficient (I2C2) and intraclass correlation coefficient (ICC) respectively. For each of these perspectives, we repeated our analyses using 10-, 20-, 30-, 40- and 50 min to establish data needs to achieve moderate to high test–retest reliability. To ensure the generalizability of our findings, we replicate our analyses in the publicly available eNKI-TRT dataset, which made use of a multiband R-fMRI sequence (TR  = 645 ms, duration = 10 min). To provide a demonstration of the potential biological determinants of inter-individual variation in R-fMRI derived indices of areal organization, we demonstrated potential familial or genetic associations using a recently published R-fMRI dataset (Yang et al. 2016; Sinclair et al. 2015; TR = 2100 ms, duration = 5 min) obtained from 136 pairs of twins.

Of note, while initial work has based the determination of gradients and boundaries based upon a summary index for iFC, a multitude of alternative indices could be used. As such, additional analyses in the present work explore the similarities and distinctions among findings obtained using alternative iFC indices for the definition of gradients. Finally, we validated that the transitions and boundaries are not derived by the underlying structural architecture by comparing the results from real data to those from random data.

Methods

Data Acquisition

Primary analyses in the present study were carried out using the Hangzhou Normal University (HNU) test–retest dataset made publicly available via the Consortium for Reliability and Reproducibility (CoRR: http://fcon_1000.projects.nitrc.org/indi/CoRR/html/data_citation.html) (Zuo et al. 2014). The dataset consists of 300 R-fMRI scans, collected from 30 healthy participants (15 males, age = 24 ± 2.41 years) who were each scanned every three days for a month (10 sessions per individual) (Chen et al. 2015). Data were acquired at the Center for Cognition and Brain Disorders at Hangzhou Normal University using a GE MR750 3 Tesla scanner (GE Medical Systems, Waukesha, WI). Each 10-min R-fMRI scan was acquired using a T2*-weighted echo-planar imaging sequence optimized for blood oxygenation level dependent (BOLD) contrast (EPI, TR = 2000 ms, TE = 30 ms, flip angle = 90°, acquisition matrix = 64 × 64, field of view = 220 × 220 mm2, in-plane resolution = 3.4 mm × 3.4 mm, 43 axial 3.4-mm thick slices). A high-resolution structural image was also acquired at each scanning session using a T1-weighted fast spoiled gradient echo sequence (FSPGR, TE = 3.1 ms, TR = 8.1 ms, TI = 450 ms, flip angle = 8°, field of view = 220 × 220 mm, resolution = 1 mm × 1 mm × 1 mm, 176 sagittal slices). Foam padding was used to minimize head motion. Participants were instructed to relax during the scan, remain still with eyes open, fixate on a “+” symbol, stay awake and not think about anything in particular. After the scans, all participants were interviewed to confirm that none of them had fallen asleep. Data were acquired with informed consent and in accordance with ethical committee review.

Two additional collections were included for secondary analyses. First, the Enhanced NKI-Rockland Sample test–retest dataset (eNKI-TRT), which is also available via CoRR (Zuo et al. 2014). Second, a recently published twin dataset obtained as part of the Queensland Twin Imaging (QTIM) study (de Zubicaray et al. 2008). The data acquisition and image preprocessing have been described in previous studies (Nooner et al. 2012; Zuo et al. 2013; Yang et al. 2016) and thus are only briefly outlined here.

The eNKI-TRT sample consists of data obtained from 22 participants (16 males, age = 33 ± 12.24 years) who were each scanned twice (1 week apart) at the Nathan S. Kline Institute for Psychiatric Research (Nooner et al. 2012). Each session includes a 10-min R-fMRI scan (multiband EPI, TE = 30 ms, TR = 645 ms, flip angle = 65°, acquisition matrix = 112 × 112, in-plane resolution = 3 mm × 3 mm, 64 axial 3-mm thick interleaved slices, number of measurements = 404, multiband factor = 4) and a structural MRI (MPRAGE, TE = 2.52 ms, TR = 1900 ms, TI = 900 ms, flip angle = 9°, resolution = 1 mm × 1 mm × 1 mm). Data acquisition was performed on a Siemens 3 T scanner with a 32-channel head coil located at NKI.

The QTIM study consists of data obtained from 200 twin pairs, including 86 monozygotic (MZ) and 114 dizygotic (DZ) twin pairs using a 4 T Bruker Medspec scanner (Bruker Medical). For each participant, the data consisted of 150 volumes of T2*-weighted BOLD contrast (EPI, TR = 2100 ms, TE = 30 ms, flip angle = 90, acquisition matrix = 64 × 64, field of view = 230 mm, 36 axial 3.6-mm thick slices). High-resolution T1-weighted images were also acquired (MPRAGE, TE = 3.83 ms, TR = 2500 ms, TI = 1500 ms, flip angle = 15°, slice thickness = 0.9 mm, acquisition matrix = 256 × 256 × 256).

Image Preprocessing

All data were processed using the Connectome Computation System (CCS) (https://github.com/zuoxinian/CCS), which provides a platform for multimodal image analysis by combining components of AFNI (Cox 1996), FSL (Smith et al. 2004; Jenkinson et al. 2012), FreeSurfer (Dale et al. 1999; Fischl et al. 1999a), and SPM (http://www.fil.ion.ucl.ac.uk/spm). The CCS includes various implementations of quality control, surface-based R-fMRI measures, reliability and reproducibility assessments (Xu et al. 2015). Below summarize the structural and functional steps of the analyses employed.

Structural MRI preprocessing included spatial noise removal by a non-local mean filtering operation (Xing et al. 2011; Zuo and Xing 2011), followed by brain extraction, tissue segmentation, and surface reconstruction, as implemented in the standard recon-all processing pipeline in FreeSurfer 5.1 (http://freesurfer.net/fswiki). Tissue segmentation resulted in gray matter (GM), white matter (WM) and cerebrospinal fluid (CSF) masks for each hemisphere. For each participant, triangular meshes representing white matter and pial surfaces were reconstructed by tessellating the GM-WM and GM-CSF interfaces, and averaging the surfaces (white matter, pial) to create a middle cortical surface (Dale et al. 1999; Fischl et al. 1999a). The resulting surface in native space was then inflated into a sphere for alignment to the fsaverage template by shape-based spherical registration (Fischl et al. 1999b). Following calculation of the transform for registration to the fsaverage, which has 164 000 vertices, we then down-sampled the surface to a representation with 10 242 vertices (fsaverage5) to ensure computational feasibility of the functional analyses to be performed in the present work. Quality control was carried out using a series of screenshots generated by CCS to facilitate visual inspection of the outputs from essential processing procedures. In particular, for the preprocessing of a given structural MR image to be considered valid, the outputs of the following steps had to pass inspection: brain extraction, tissue segmentation, surface reconstruction, co-registration and normalization.

Functional preprocessing included: discarding the first five volumes of the R-fMRI time series and storing all images as double precision floating point numbers (to reduce variability of outcomes across operating systems, see Glatard et al. 2015), detecting and compressing temporal spikes (AFNI 3dDespike), slice timing correction (this was not performed for multiband images), motion correction (3DVolReg), and 4D intensity normalization to 10 000. Nuisance variable regression (Fox et al. 2005; Lund et al. 2006) was performed to remove the mean time series for individual-specific WM and CSF masks (derived from the FreeSurfer segmentation), as well as 24 motion parameters (Friston et al. 1996; Yan et al. 2013); linear and quadratic trends were also removed from the data. The time series residuals were band-passed filtered (0.01–0.1 Hz) to restrict the signals to those frequencies previously implicated in resting state functional connectivity (Biswal et al. 1995; Cordes et al. 2001). Functional images were co-registered to the native high-resolution anatomical images using 6-parameter boundary-based registration (BBR) in FreeSurfer (Greve and Fischl 2009).

Quality Control Procedure

Following functional preprocessing, the frame-wise displacement (FD) and mean FD were calculated to quantify the head micromovements (Power et al. 2012, 2014; Patriat et al. 2013). Any scans with mean FD greater than 0.2 mm were excluded from our analysis. For the test–retest HNU dataset employed in our primary analyses, mean-FD was less than 0.2 mm for all 300 scans (mean = 0.057, SD = 0.019). In addition, the percentage of excessive motion frames (FD > 0.2 mm) for each of the scans was found to be less than 1%.

For the eNKI-TRT, 20 participants (age = 33.40 [SD = 12.58], 14 male) passed our criteria (mean FD < 0.2 mm) and were included in the current studies (mean = 0.107, SD = 0.036). For the QTIM twin dataset, 379 of 400 participants passed the head motion criterion (mean FD < 0.2 mm). Given that previous studies have raised concerns about the potential contributions from gender difference in brain structure (e.g. cortical area, thickness, gyrification, etc.) and head motion (Schmitt et al. 2007; Chen et al. 2013; McKay et al. 2014; Docherty et al. 2015) to findings in twin studies, Yang et al. (2016) included only same sex adult twin pairs in which both twins exhibited low motion (mean FD < 0.2 mm). We used the same criteria and resulting participant list in the present work. This resulted in inclusion of 272 participants (age = 22.12 [SD = 2.35], 68 male), consisting of 78 MZ twin pairs and 58 DZ twin pairs in our analysis. These criteria controlled the potential heritability of the head motion and spatial information; no significant differences between MZ and DZ pairs were found in the mean FD (p = 0.52); only a non-significant trend was noted between MZ and DZ pairs with respect to the level of spatial correlation for anatomical images in MNI space (p = 0.08).

Computation of iFC Metrics on Native Surface

There are several advantages for analyzing fMRI data on the surface rather than in volume space. A cortical sheet provides a more accurate representation of the morphology and topology of brain structure (Dale et al. 1999). It provides a more accurate registration between individual data and the template surface (Van Essen 2004; Ghosh et al. 2010; Klein et al. 2010; Yeo et al. 2010). Also, it enables the visualization of spatial relationships between brain regions in terms of their geodesic distances along the cortical surface, which is more neurobiologically meaningful than 3D Euclidean distance in volume space. In light of these putative benefits, all functional indices in the present study were calculated on the surface. First, the volumetric fMRI data were aligned to anatomical space using the BBR transformation matrix and then projected to their corresponding native middle cortical surface (about 160 k vertices for each hemisphere). To make the resolution comparable between the volumetric and surface data, a coarser fsaverage5 template surface was used (the average triangle edge length between adjacent vertices of fsaverage5 is about 3.5 mm) for computing functional metrics. The transformation calculated during spherical registration was applied to the fsaverage5 template to convert it into a standard mesh of the native surface, for which each node has a direct correspondence with a node on the fsaverage5 surface (Van Essen et al. 2012b). The standard mesh was calculated by interpolating the native surface to a transformed version of the fsaverage5 template in native space.

Consistent with previous studies of connectivity gradients, iFC similarity was measured by the spatial similarity between whole brain functional connectivity maps calculated from a seed vertex and those calculated from every other voxel‘s time courses (Cohen et al. 2008; Wig et al. 2011, 2014). Specifically, for each individual, the time course for each vertex was extracted and used to calculate a whole brain functional connectivity profile that consists of 20 484 vertices (10 242 for each hemisphere) of cerebral gray matter and 9413 voxels in subcortical regions and the cerebellum. The distributions of the resulting correlation values were standardized to the normal distribution using Fisher‘s r-to-z transform. For each hemisphere, an iFC similarity map was measured for each vertex by calculating the spatial correlation between the vertex‘s iFC profile and the iFC profile of every other vertex – resulting in a 10 242 vertices × 10 242 vertices symmetric matrix. Each column (or row) of this matrix represents the vectorized iFC similarity map for each surface vertex. These 10k iFC similarity maps were then smoothed along the native surface with a Gaussian Kernel (FWHM = 8 mm, about twice that of the surface resolution [~3.5 mm], see Supplementary Figure S1–2 to address the smoothness effect). The gradient (i.e., the first spatial derivative) of each smoothed iFC similarity map was computed on the native middle surface to measure transitions in iFC profile across vertices, resulting in 10k gradient maps for each hemisphere. The details of the gradient computation are described in the next section.

Beyond iFC similarity, we adopted three additional functional indices that capture local-, global- and network-scale aspects of intrinsic brain function, respectively, and can each be used to generate gradient maps (independent of iFC similarity). Two-dimensional Regional Homogeneity (2dReHo) was employed as a local-level measure of iFC that characterizes the temporal synchronization of the BOLD signal at each vertex across the cortical mantle (Zang et al. 2004; Zuo et al. 2013; Jiang et al. 2014; Jiang and Zuo 2015). For each vertex, ReHo was defined by the Kendall‘s coefficient of concordance of the time series within the given vertex and its one-step neighboring vertices. The preprocessed R-fMRI data for ReHo calculation was not spatially filtered, but temporally band-pass filtered to avoid artificial increases of ReHo intensity (Zuo et al. 2013).

The network-scale measurement of iFC was captured using a dual regression (DR) procedure to map 10 intrinsic connectivity networks previously defined from a meta-analysis of activations from the BrainAtlas database (Smith et al. 2009). These networks includes Medial Visual network (DR-MedVis), Occipital Visual network (DR-OccVis), Lateral Visual network (DR-LatVis), Default Mode network (DR-DMN), Cerebellar network (DR-Cerebellar), Sensorimotor network (DR-SenMot), Auditory network (DR-Audi), Executive Control network (DR-Control), Left Frontoparietal network (DR-FrontL), and Right Frontoparietal network (DR-FrontR). The spatial network maps were conducted in the first regression on normalized R-fMRI data in MNI152 space to generate the characteristic time series of 10 networks for each individual. The resultant characteristic time series were further entered into the second   regression to yield the spatial map of 10 networks on cortical surface in individual level.

At the full-brain scale, Degree centrality (DC) and Eigenvector centrality (EC) were employed to capture features of iFC as a whole (i.e., full-brain connectivity patterns) (Bullmore and Sporns 2009; Rubinov and Sporns 2010; Zuo et al. 2012). Specifically, a pair-wise iFC matrix was calculated for all possible pairings of vertices on the cortical surface using Pearson‘s correlation. We then thresholded and binarized the full brain iFC matrix at p(r) < 0.001.

Transition Zones, Boundaries and Parcellation

Gradients were calculated on the convoluted (i.e., not flattened) standardized native surface mesh, followed by an edge detection procedure that sharpened the gradients to identify putative boundaries. The initial framework for gradient mapping follows the Canny edge detection procedure (Canny 1986) and involves three steps: (1) each vertex‘s iFC connectivity profile is smoothed with a 2D Gaussian filter to reduce noise, (2) a gradient map is calculated for each vertex‘s smoothed profile to identify regions of rapidly changing IFC profiles, and (3) performing edge detection by finding local maxima in the gradient maps. Cohen et al. (2008) first applied Canny edge detection to iFC similarity profiles calculated from data projected onto flattened surfaces. Flattening the cortical surface requires it to be cut in five places; findings associated with edges near these cuts are prone to artifact and are thus difficult to interpret. To avoid this issue, Wig et al. (2013, 2014) extended the method to un-flattened surfaces by performing all three steps (smoothing, gradient computation, edge detection) on the middle cortical surface (mid-way between white matter and pial surfaces). The major challenge of this approach is that the gradient calculation must account for the curvature of the cortical surface. This was handled by calculating the gradient of each vertex on a plane perpendicular to the surface normal vector extending from that vertex (Caret 5.65 “metric-gradient” function http://brainvis.wustl.edu/wiki/index.php/Caret:About). The gradient is calculated by solving a linear model relating the geodesic distances between the vertex and its nearest neighbors to differences in their image intensity. The magnitude of the resulting gradient is the change of image intensity per unit distance. Here, the intensity image is a functional feature of the cortical surface described in the above section, including iFC similarity, ReHo, DC, EC and DR-Networks. The gradients were calculated for each of these functional indices to detect abrupt changes at the local, global, and network scales.

Finally, the parcels and corresponding edges were identified from the resulting gradient maps using a “watershed by flooding” algorithm (Gordon et al. 2014). This algorithm identifies the local minima on a gradient map as starter seeds for parcel creation, and then gradually grows the seeds outward until they meet another parcel. The barriers built in this procedure for separating parcels are identified as edges. With the exception of iFC similarity, an individual-specific binary edge map was generated for each gradient map (ReHo, DC, EC, and DR-networks).

For the iFC similarity map, a more complex set of procedures is carried out, generating an edge density map for each individual, rather than a simple binary edge map. Specifically, for each vertex, a gradient map was calculated, as well as a corresponding edge map using the watershed algorithm – resulting in 10k gradient/edge maps for each hemisphere. For each individual, the 10k gradient/edge maps were then averaged to generate a final gradient map and an edge density map. While an edge density map can be more preferable depending on the specific application, it is also possible to generate a binary map by applying the same watershed algorithm to the edge density map.

To make the resultant maps from ReHo, DC, EC and DR-networks more comparable to those iFC similarity, we created edge density maps for these R-fMRI measures by bootstrapping. Specifically, we generated alternative observations of the data using circular block bootstrap to ensure the temporal dependencies of R-fMRI data (Bellec et al. 2006, 2010). The procedure concatenated temporal blocks of original data (block length = the square root of the number of total time points) to form a time series with same temporal points as the original. Each temporal block was randomly extracted from the original time series with replacement. We utilized 100 bootstraps to yield 100 edge maps. The final edge density maps were calculated by averaging 100 bootstrapped binary edge maps for ReHo, DC, EC and DR-networks.

For each of the R-fMRI measures (iFC similiarity, ReHo, DC, EC, DR-networks), the resulting discrete parcels were assigned to large-scale system networks driven using group-averaged parcellation generated by Gordon et al. (2015). The assignment procedures are described in the next section.

Network Systems Assignment

Recent studies have demonstrated that the network structures of parcellation in individuals were similar to the group averaged system (Gordon et al. 2015; Laumann et al. 2015). To characterize large-scale systems for each edge map generated in the present work, we assigned system identity to parcels created in the above procedure using the previously established network definitions from Gordon et al. (2014). The matching procedure was modified based on Gordon et al. (2015). Specifically, we first detected the vertex with local minimal edge density for each discrete parcels. Then we averaged that time series with the time series of its one-step neighbors within that parcel. The averaged time series was considered as a representative time series of the chosen parcel, and correlated against all other time series across the cortical surface to obtain a connectivity map. After that, we thresholded and binarized the connectivity map at the top 5% of connectivity strengths. This resulted in a binarized map of regions with high connectivity to that local minimal vertex. Then we measured the overlap of this binarized map to the 12 network maps from Gordon, the best match network with highest Dice similarity was assigned to that vertex and its representative parcels. Consistent with prior work, the neighboring parcels with the same network identifier were merged together to create a “network-patch”.

Estimating Intra-, Inter-Individual Variability and Intra-Class Correlation

Intra- and inter-individual variability were simultaneously estimated using separate Linear Mixed Effects Models (LMEs) for gradient and edge density maps. The gradient or edge density for a given vertex v can be denoted as Yij(v), where i indicates the participant and j represents the measurement (for i = 1, 2,…, 30 and j = 1, 2). 

(1.1)
Yij(v)=μ00(v)+γi0(v)+mYij+meanFDij+mcBBRij+agei+genderi+Jacobiani(v)+ϵij(v)
In this model, μ00(v) is a constant term that represents the intercept or fixed effect of the group average in gradient/edges at vertex v, while γi0(v) is the random effect term for ith participant at vertex v. Due to the potential influence of various confounding factors to our final estimation, we included several covariates in this mixed model at both the session and individual level. At the session level j, the model included three covariates: the mean frame-wise displacement (meanFDij) of head motion, minimum cost of BBR (mcBBRij), and the global mean of gradients/edges across the whole brain within group mask (mYij) for each participant i. At the individual level, agei, genderi and Jacobian determinant of the spherical transformation (Jacobiani(v) at vertex v) were included for each participant i as well. The model estimations were implemented using the lme function from the nlme R package (http://cran.r-project.org/web/packages/nlme). The variance estimation σγ2 of the random effect term γi0(v) is the inter-individual variance across all the participants, while the variance σϵ2 of the residual ϵij is the intra-individual variance for a single participant across all the sessions. Meanwhile, the intra-class correlation (ICC) was calculated by dividing inter-individual variance to the sum of inter- and intra-individual variance.

Surface Geometry Constraints Testing

To determine whether transition zones and detected boundaries could have arisen from the surface geometry, we replicated the boundary mapping analyses using surrogate data. Specifically, for each participant, we replaced the preprocessed R-fMRI data with temporal Gaussian random data (mean = 0, standard deviation = 1, the same number of time points as the original data) in native functional volume space. The surrogate data were then projected onto the participant‘s native middle surface, registered and down resampled to “fsaverage5” version of native surface using the same co-registration and spherical registration transformation calculated during processing of the real data. The functional measures (iFC similarity) were computed, followed by gradient computation and edge detection using the same procedure and registration transformations used for the real data. The gradient and boundaries maps from this surrogate data were compared with the results from real data to investigate whether surface geometry and registration errors would contribute to the functional areal organization detected.

Results

Boundary Map Characteristics

For each of the 300 scans (30 participants, 10 sessions each) from the HNU data collection, a gradient map and a corresponding edge map were generated. This was accomplished using the gradient-based boundary mapping approach and iFC similarity measure specified by Wig et al. (2014), though at the individual level. Specifically, for each individual, the maps generated were averages of those calculated for each of the 20k vertices (10k per hemisphere). Prior to examining differences among participant-specific maps, individual level maps were averaged to generate a group-level gradient map that was strikingly similar to those previously published by Wig et al. (2014) and Gordon et al. (2014) (Fig. 1A). To measure the similarity of our group map versus previous work, we calculated the percentile of the group-level gradient scores for the vertices that were identified as brain systems boundaries in Gordon et al. (2014). The median percentile of gradient within Gordon's system boundaries was 70.62%. The group-level edge density map was also similar to those previously published, though somewhat less comparable (the median percentile of edge density was 63.79% within Gordon's boundaries); this was largely due to the fact that we calculated the edge maps on individual participants prior to averaging, as opposed to group-average maps (e.g., Gordon et al. calculated an edge map for each vertex using the group-level gradient map, and then averaged across edge maps). The latter approach averages out the contributions of individual variation in the gradients prior to the generation of the edge maps (see supplementary materials Figure S3 for a direct comparison of the two approaches from prior methodological testing by our group using the QTIM dataset).

Figure 1.

Individual areal organization is unique to each individual and highly repeatable across scan sessions. (A) Group-level gradient and edge density maps for iFC similarity. The individual maps were calculated on the standardized “fsaverage5” version of native surface for each individual, then averaged and visualized here on the “fsaverage5” template surface. (B) Between-scan spatial correlation matrices for the gradient (top) and edge density (bottom) maps derived from the 300 datasets (30 participants * 10 sessions); the rows/columns of the matrix are ordered by participant (i.e., first ten rows are sessions 1–10 for participant one, second ten rows are sessions 1–10 for participant two and so on; total number of participants = 30). (C) Distribution of between-scan spatial correlations (from B). Black dots are the within-individual (i.e., same participant, different sessions) correlations for gradient (top) and edge density maps (down), while the gray dots are between-individual correlations (i.e., different participants, same session, or different participants, different sessions).

Figure 1.

Individual areal organization is unique to each individual and highly repeatable across scan sessions. (A) Group-level gradient and edge density maps for iFC similarity. The individual maps were calculated on the standardized “fsaverage5” version of native surface for each individual, then averaged and visualized here on the “fsaverage5” template surface. (B) Between-scan spatial correlation matrices for the gradient (top) and edge density (bottom) maps derived from the 300 datasets (30 participants * 10 sessions); the rows/columns of the matrix are ordered by participant (i.e., first ten rows are sessions 1–10 for participant one, second ten rows are sessions 1–10 for participant two and so on; total number of participants = 30). (C) Distribution of between-scan spatial correlations (from B). Black dots are the within-individual (i.e., same participant, different sessions) correlations for gradient (top) and edge density maps (down), while the gray dots are between-individual correlations (i.e., different participants, same session, or different participants, different sessions).

High Reproducibility of Individual Areal Boundaries (Fingerprinting)

For both gradient maps and edge maps, we next investigated the similarity of the maps generated across individuals and time (i.e., scan sessions). In order to accomplish this, we calculated the spatial correlation between the 300 gradient maps generated (30 participants * 10 sessions), as well as between the 300 edge maps. Figure 1B demonstrates the spatial correlation matrix obtained for gradients (upper panel) and edge maps (lower panel). Regardless of participant, the correlations between different scan sessions (gradient: mean r = 0.609 (SD = 0.110); edge: mean r = 0.450 (SD = 0.093)) are notably greater than those between participants (gradient: mean r = 0.318 (SD = 0.078); edge: mean r = 0.123 (SD = 0.037)). To facilitate appreciation of this point, the correlation matrix is ordered by participant (i.e., first 10 rows belong to participant 1, second 10 rows belong to participant two, and so on). As can be seen in Figure 1C, the distribution of between-individual correlation scores for gradient maps overlap minimally with that for within-individual; the correct identification rate (calculated as the percent of within-individual correlations exceeding between-individual correlations) was 90.44% (1221 in total 1350 within-individual correlations of 300 scans exceeded any between-individual correlations); For edge maps, fewer overlaps were present for within- and between-individual correlations, resulting in a 99.26% correct identification rate. These findings suggest that the full-brain areal transition characteristics indexed by gradient and edge measures are unique to each individual, and repeatable across scan sessions.

Figure 2 depicts data from multiple 10-min sessions for three representative participants to provide illustrative examples of distinct patterns associated with each participant, as well as their repeatability. The distribution of within- and between- individual spatial correlations for each individual was plotted in Figure S4. We further examined whether the within-individual spatial correlations were related to age, gender, head motion or the functional-to-structural co-registration. We did not find any significant relationships between the mean or standard deviation of within-individual correlations and the factors listed above.

Figure 2.

Examples of gradient and edge density of iFC similarity in 3 individual participants. Left three columns are gradients from three sessions and right two columns are gradients from two 50-min subsets. To facilitate comparison, for each participant we mark off (using colored boxes) an example of a feature that differs across participants and is relatively consistent within individual.

Figure 2.

Examples of gradient and edge density of iFC similarity in 3 individual participants. Left three columns are gradients from three sessions and right two columns are gradients from two 50-min subsets. To facilitate comparison, for each participant we mark off (using colored boxes) an example of a feature that differs across participants and is relatively consistent within individual.

We replicated the general pattern of findings in another independent R-fMRI dataset (eNKI-TRT data collection: 20 participants * 2 sessions, multiband EPI sequence TR = 645 ms). Specifically, we again observed substantially higher within-individual correlations (gradient: mean r = 0.56 (SD = 0.159); edge: mean r = 0.364 (SD = 0.095)) than between-individual correlations (gradient: mean r = 0.219 (SD = 0.129); edge: r = 0.075 (SD = 0.032)).

High Reproducibility of Individual Networks

Using Dice Similarity, we assessed the within- and between-individual reproducibility of the patches associated with each of the 12 networks defined by Gordon et al. (2014). Five of the 12 networks (visual, auditory, dorsal sensorimotor, ventral sensorimotor, default) exhibited relatively high reproducibility across both sessions and individuals (i.e., average pairwise Dice Similarity > 0.5) (see Fig. 3) – suggesting a relatively conserved architecture for these networks.

Of note, for 11 of the 12 networks, within-individual similarity was found to be greater than between-individual similarity (Fig. 3) – suggesting that despite similarities, consistent inter-individual differences exist. The salience network was the notable exception. This network could not be successfully assigned, largely due to the small size of the patches comprising the salience network and their spatial variation in individuals.

Figure 3.

The reproducibility of system-matched patches across 10 sessions and 30 participants. The color-coded bars indicate Dice similarity of system patches (mean and standard deviation) for within and between individual at each system network.

Figure 3.

The reproducibility of system-matched patches across 10 sessions and 30 participants. The color-coded bars indicate Dice similarity of system patches (mean and standard deviation) for within and between individual at each system network.

Individual Areal Boundaries Have Familial and Genetic Associations

To demonstrate the feasibility of relating measures of similarity for areal organizations to known sources of inter-individual or phenotypic or biological variation, we made use of the R-fMRI dataset (Yang et al. 2016) generated as part of QTIM (de Zubicaray et al. 2008). Specifically, we examined whether inter-individual differences in full-brain areal measures (gradient maps, edge maps) could be related to familial (i.e., sibling vs. non-sibling pair) or genetic (monozygotic twins vs. dizygotic twins) relationships. To accomplish this, for every pairing of individuals in the sample, we calculated the spatial correlation between the iFC similarity-based gradients maps, as well as for the edge maps. Our a priori predictions were that gradient/edge maps would be more similar among siblings than non-siblings, and more similar among monozygotic twin pairs than dizygotic twin pairs. As expected, the correlation between siblings was significantly higher than non-siblings (gradient: p < 10–23; edge map: p < 10–50) (see Fig. 4). The correlation between MZ pairs was significantly higher than DZ pairs for edge maps (p = 0.002), and marginally significant for gradient maps (p = 0.075).

Figure 4.

Spatial correlations of gradients and edge density in iFC similarity between MZ, DZ, and unrelated pairs. Blue lines are the mean correlation across MZ pairs for gradient (A) and edge density (B) maps. Red lines are the mean correlation across DZ pairs for gradient (A) and edge density (B). Green lines are the mean correlation across all unrelated pairs of QTIM dataset for gradient (A) and edge density (B).

Figure 4.

Spatial correlations of gradients and edge density in iFC similarity between MZ, DZ, and unrelated pairs. Blue lines are the mean correlation across MZ pairs for gradient (A) and edge density (B) maps. Red lines are the mean correlation across DZ pairs for gradient (A) and edge density (B). Green lines are the mean correlation across all unrelated pairs of QTIM dataset for gradient (A) and edge density (B).

Given that a trend was found between MZ and DZ pairs with respect to the level of spatial similarity for anatomical images (p = 0.08), we further tested whether the observed heritable functional transitional profiles could have been confounded by structural information. However, we did not find any significant correlations between similarity of structure and gradient/edges in twin pairs (all p > 0.3). The significance comparisons between MZ and DZ pairs showed almost identical p-values after being controlled for structural similarity (gradient: p = 0.072, edge: p = 0.002).

Data Requirements for Mapping Inter-Individual Differences in Areal Organization

Recent work has emphasized the value of using longer data acquisitions (or combinations of shorter ones) to increase the stability of estimates of functional connectivity. In particular, using a correlation matrix based on 380 min of data from a single participant as a reference, Laumann et al. (2015) found that the consistency of parcel-based full-brain correlation matrix estimates (with respect to the reference) increased substantially until 27 min of data were included; at this point, gains were notably more modest, with estimates nearing an asymptote at 90 min. In the present work, we attempted to further inform this issue by examining the stability of gradient and edge maps, as data is added in 10-min increments. More specifically, first, for each participant, we randomly selected 5 of the 10 sessions and used their data to derive reference gradient and edge maps. Then, for each participant, we randomly selected samples of 10, 20, 30, and 40 min of data from the remaining five 10-min sessions, for spatial correlation of the resultant gradient/edge maps with the same participant‘s 50-min reference maps.

Regardless of whether 10, 20, 30, 40 or 50 min of data were used, the resultant gradient or edge map for an individual consistently showed a higher spatial correlation with the reference map calculated for the same participant, than with the reference map for any other participant. In other words, for all participants, one could readily distinguish whether the map derived was properly matched with that from the same participant as opposed to different participants. The distribution of the spatial correlations to the reference gradient or edge map is demonstrated in Figure 5.

Figure 5.

The stability of individual gradient and edge maps as a function of the amount of data used for gradient estimation. (A) Each dotted-line represents a different participant; the spatial correlation between maps generated from each incremental amount of data (10, 20, 30, 40 and 50 min) and those from the remaining 50 min session data (reference) are depicted (top: gradient; bottom: edge density). Red dots are the mean correlation across participants at each data subsets (B) The distribution of within-individual correlations (black dots) and between-individual correlations (gray dots) at selected data subsets to the 50-min references for each individuals.

Figure 5.

The stability of individual gradient and edge maps as a function of the amount of data used for gradient estimation. (A) Each dotted-line represents a different participant; the spatial correlation between maps generated from each incremental amount of data (10, 20, 30, 40 and 50 min) and those from the remaining 50 min session data (reference) are depicted (top: gradient; bottom: edge density). Red dots are the mean correlation across participants at each data subsets (B) The distribution of within-individual correlations (black dots) and between-individual correlations (gray dots) at selected data subsets to the 50-min references for each individuals.

It is important to note that the consistency of findings obtained across the various quantities of data should not be taken to infer that there were no differences. Consistent with predictions based upon prior work, across participants, each 10-min increment in data produced an incremental increase in the mean correlation with the reference images – for both gradient and edge maps (Fig. 5A). Focusing first on gradient maps, we found that the average within-participant correlation between the map generated from a single scan session (10 min, 295 time points) and an individual‘s reference image was r = 0.70 (SD = 0.08); this progressively increased to a high of r = 0.86 (SD = 0.05) when data from five sessions (50 min, 1475 time points) was used (Fig. 5B, upper). Statistical testing found that the observed increase was significant for each 10-min increment (p < 0.001), except for the 20- to 30-min increment (p = 0.108). For edge density, we found that the average within-participant correlation between the map generated from a single scan session and an individuals' reference image was r = 0.44 (SD = 0.06); this progressively increased to a high of r = 0.68 (SD = 0.07) with 50 min (1475 time points) (Fig. 5B, lower). Statistical testing found that increases in correlation for each 10-min increment were significant (all p < 0.001).

To gain greater insights into the influence of scan time in terms of test–retest reliability, we calculated the ICC of gradient and edge density maps for each of the scan durations of 10 min, 20 min, 30 min, 40 min and 50 min. This was accomplished by randomly selecting 1, 2, 3, 4, or 5 10-min sessions of the 10 available sessions for each participant (i.e., 10, 20, 30, 40 or 50 min of data, respectively). As shown in Figure 6, as a measure of functional local transition, the gradient of iFC similarity achieved good test–retest reliability (ICC > 0.5 for over 50% vertices across cortex) with a scan duration of 20 min, which is consistent with prior work (Laumann et al. 2015). A relationship was observed between gradient strength and ICC across vertices (see Supplementary Figure S5).

Figure 6.

The influences of scan time on test–retest reliability of gradient and edge density of iFC similarity. The reliability of gradient and edge density were estimated as intra-class correlation maps for each of the scan durations of 10 min, 20 min, 30 min, 40 min and 50 min (top panel). To access the influences of scan durations on the ICC distribution, the histograms of ICC scores are plotted for each scan of the durations in the bottom.

Figure 6.

The influences of scan time on test–retest reliability of gradient and edge density of iFC similarity. The reliability of gradient and edge density were estimated as intra-class correlation maps for each of the scan durations of 10 min, 20 min, 30 min, 40 min and 50 min (top panel). To access the influences of scan durations on the ICC distribution, the histograms of ICC scores are plotted for each scan of the durations in the bottom.

To provide a single measure of reliability, we complemented the spatial correlation-based measures of repeatability with image intra-class correlation coefficient (I2C2) (Shou et al. 2013). The I2C2 was increased with scan time for gradient (10 min: 0.42, 20-min: 0.54, 30-min: 0.61, 40-min: 0.67, 50-min: 0.73) and edge density (10 min: 0.35, 20-min: 0.48, 30-min: 0.55, 40-min: 0.61, 50-min: 0.65). In addition, we replicated the I2C2 findings in eNKI test–retest dataset (multiband EPI sequence TR = 645 ms, 10 min, 895 time points). The I2C2 was nearly identical to that observed for the 10-min HNU dataset (0.49 for gradient and 0.34 for edge density). However, the distributions of vertex-wise ICC values for gradient and edge maps derived from the eNKI-TRT were superior to those observed for the 10-min HNU scans (gradient maps: mean ICC [eNKI-TRT] = 0.40 [SD = 0.22], mean ICC [HNU] = 0.37 [SD = 0.18]; edge density: mean ICC [eNKI-TRT] = 0.30 [SD = 0.21], mean ICC [HNU] = 0.32 [SD = 0.18]); see supplementary Figure S6.

Within- and Between-Individual Variability

It is generally accepted that anatomical components of higher order association areas tend to have greater anatomical variability among individuals than those comprising lower-order sensory and motor areas (Mueller et al. 2013). To gain similar insights into potential regional variation in areal organization, we examined the within- and between-participant variation separately. Figure 7A and B depicts the percentage of the total variability that is the intra- vs. inter-individual variability (after accounting for that attributable to nuisance signals) for each, gradient maps and edge maps. Intra-individual variability indexes the temporal stability of transition scores within participant while inter-individual variability indexes the stability of transition scores across the participants. Both intra- and inter-individual variability demonstrated a non-uniform distribution across brain regions. For both gradients and edges, between-individual differences were largest in high-order association cortex including the lateral prefrontal lobe, lateral parietal lobe, and around the border of the PCC-precuneus default mode networks, while lower in sensory-motor cortices and medial occipital visual cortices. Within-individual variability depicted the complementary pattern from inter-individual variance, with the maximal variation of gradient and edges being noted within sensory-motor and visual regions. Regions within the default-mode network demonstrated a moderate level of within-individual variability.

Figure 7.

Within- and between-individual variability for gradient and edge density maps based on iFC similarity. (A) The percentage of within individual variation (relative to total variation) in gradients (top) and edge density (down) scores (based on IFC similarity) across sessions. (B) The percentage of between individual variation in gradient (top) and edge density (down) scores (based on iFC similarity). (C) Within- and between- individual variation (standard deviation) at network borders (defined as confidence <0.3 along the network borders) and within network region (defined as confidence > = 0.3 within each network). The limbic network was not included here due to the substantial signal loss in temporal areas.

Figure 7.

Within- and between-individual variability for gradient and edge density maps based on iFC similarity. (A) The percentage of within individual variation (relative to total variation) in gradients (top) and edge density (down) scores (based on IFC similarity) across sessions. (B) The percentage of between individual variation in gradient (top) and edge density (down) scores (based on iFC similarity). (C) Within- and between- individual variation (standard deviation) at network borders (defined as confidence <0.3 along the network borders) and within network region (defined as confidence > = 0.3 within each network). The limbic network was not included here due to the substantial signal loss in temporal areas.

To provide a network-level perspective of intra-individual variations in gradient and edge densities we made use of the confidence maps for the 7 networks in Yeo et al. (2011) (note: we used the Yeo et al. (2011) network definitions rather than those of Smith et al. (2009) networks due to: (a) the lack of confidence maps, and (b) overlaps that occur among the network maps when transformed into surface space). Specifically, the standard deviations were averaged at network borders (defined as confidence < 0.3 along the network borders) and within network region (defined as confidence >=0.3 within each network). Without exceptions, the between-individual differences contributed more than 60% of the variance in transitions at the border region of 7 networks (Fig. 7C). For the within network regions, the between-individual differences still accounted for more variance than within-individual differences, however, the somatomotor and visual networks exhibited increases of within-individual variation to 45% and 46%, respectively (Fig. 7C). This suggested the spatial iFC profile in the somatomotor and visual networks were more state dependent than the other networks. This finding is consistent with those of Mennes et al. (2010), Craddock et al. (2013), and Laumann et al. (2015).

Transition Pattern Based on Local-, Global, and Network-Scale iFC

Previous studies using stimulus-based functional neuroimaging methods have suggested that areal discrimination maps could vary depending on the specific property or properties used in their definition, e.g. angular and eccentricity representation maps for distinct areas of early visual cortex (Wandell and Winawer 2011; Buckner and Yeo 2014). Here, we examined whether transition patterns are dependent upon the specific intrinsic brain measure used for their definition. Specifically, we repeated our analyses focused on the definition of gradient and edge maps using commonly examined measures in the literature (i.e., ReHo, DC, EC, and DR-Networks). We then measured the convergence of transition patterns for both gradients and edge density from iFC similarity, ReHo, DC, EC and DR-Networks. Figure 8 demonstrates the average pair-wise correlations of gradient between all functional metrics. Among all the metrics, the centrality maps (DC and EC) showed the most remarkably similar transition scores (r = 0.94, SD = 0.03). Edge density of iFC similarity exhibited low correlation with ReHo (mean r = 0.11, SD = 0.08) but had notable correlations with DC (mean r = 0.30, SD = 0.19), EC (mean r = 0.45, SD = 0.18), DR-MedVis (mean r = 0.26, SD = 0.15), DR-OccVis (mean r = 0.25, SD = 0.13), DR-LatVis (mean r = 0.43, SD = 0.14), DR-DMN (mean r = 0.60, SD = 0.11), DR-Cerebral (mean r = 0.32, SD = 0.13), DR-SenMot (mean r = 0.64, SD = 0.10), DR-Auditory (r = 0.61, SD = 0.12), DR-Control (r = 0.47, SD = 0.12), Dr-FrontL (r = 0.49, SD = 0.08), and DR-FrontR (mean r = 0.51, SD = 0.10). These findings are more reflective of global- and network-level iFC characteristics, than local features (Fig. 8). The Kendall Coefficient of edge density for all the functional indices was 0.29, SD = 0.05. Similar finding were observed for edge density (Figure S7).

Figure 8.

Group-level gradient maps for each functional indices from two 50-min subsets. For each participant, we calculated the spatial correlation matrix of the gradient maps for various functional indices, and then averaged across individuals to provide the spatial correlation matrix (bottom left corner). The group-level edge density maps were demonstrated in Supplementary material in Figure S7.

Figure 8.

Group-level gradient maps for each functional indices from two 50-min subsets. For each participant, we calculated the spatial correlation matrix of the gradient maps for various functional indices, and then averaged across individuals to provide the spatial correlation matrix (bottom left corner). The group-level edge density maps were demonstrated in Supplementary material in Figure S7.

Within- Versus Between-individual Variance of Transition Pattern based on Local-, Global, and Network-Scale iFC

Reminiscent of the finding for iFC similarity above, for all the functional metrics, gradient (and edge) maps calculated from different sessions in the same individual exhibited a consistently higher amount of spatial correlation with one another than those from differing individuals (Figure S8). However, with the exceptions of ReHo, DR-DMN, and DR-FrontL, the spatial correlation of gradient and edge maps generated in different sessions for an individual were generally lower than iFC similarity. Those measures characterized by a higher between-participant correlation in gradient or edge maps, typically exhibited a higher within-individual correlation (i.e., repeatability) over time. Given that the gradient of iFC similarity was averaged from 20k gradient maps across vertices, this may in part reflect the cleaning of data within an individual by averaging.

Reliability of Transition Properties based on Local-, Global, and Network-Scale iFC

While our analyses primarily focused on full-brain transition patterns, we did take the opportunity to provide insights into the test–retest reliability of vertex-based gradient and edge map scores. Specifically, we calculated vertex-wise intra-class correlation coefficients for each, gradients and edges, using the two 50-min subsets. Gradients in iFC similarity showed considerably high reliability across cortex, especially in PCC-DMN, LP-DMN, and frontoparietal association cortex (Fig. 9). The pattern was more apparent in reliability of edge maps (Figure S9). Intriguingly, highly reliable edges were found around the borders of the PCC-DMN, LP-DMN, and the frontoparietal regions, which were observed as a common transition zones in variant iFC metrics. Among different functional indices, ReHo exhibited the highest ICC in transition scores. The intra- and inter-individual variation analysis confirmed that high ICC of ReHo was due to the considerably small amount of intra-individual variation. For DR networks, the ICC of gradient and edge scores showed non-uniform spatial distribution in different networks. The vertex-wise ICC of gradient/edges scores obtained for a given DR network tended exceptionally high when looking at vertices that were part of the network (Fig. 9). Similarly, in contrast to iFC similarity, vertex-wise ICC was strongly related to gradient strength for the nearly all of the DR networks (the DR Cerebellar network is an exception) (see Supplementary Figure S10).

Figure 9.

The intra-class correlation maps for the gradient of 14 functional indices are depicted here (ICC was calculated between two 50-min subsets of data generated through random selection [without replacement]). The distribution of ICC for edge density maps for each functional metric is shown in the bottom-left corner. The ICC maps for edge density maps were demonstrated in Supplementary material in Figure S9.

Figure 9.

The intra-class correlation maps for the gradient of 14 functional indices are depicted here (ICC was calculated between two 50-min subsets of data generated through random selection [without replacement]). The distribution of ICC for edge density maps for each functional metric is shown in the bottom-left corner. The ICC maps for edge density maps were demonstrated in Supplementary material in Figure S9.

Parcel Homogeneity based on Local-, Global, and Network-Scale iFC

In evaluating and comparing the optimality of parcellation strategies, unit homogeneity is commonly viewed as a desirable property. As such, we opted to calculate the homogeneity of full-brain iFC patterns among voxels within a given parcel as a means of comparing the parcellations produced using the differing R-fMRI measures. More specifically, we first computed the gray matter (cortical and subcortical) connectivity map at each vertex using the reference 50-min data. Then for each parcel within a given brain parcellation from the other 50-min data, the Kendall‘s coefficient of concordance was calculated from the functional pattern of all vertices within a parcel. We then averaged the KCC-homogeneity scores for each approach at the individual-level and test which approach performed most accurate parcellation. Figure 10 demonstrates the KCC-homogeneity scores across all the participants. Locally Weighted Scatter-plot Smoother (LOWSS) was applied to demonstrate the average KCC-homogeneity for each of the functional indices. Parcellations generated using the iFC similarity methodology achieved the highest within-unit homogeneity, followed by the DR-network. Parcellations generated using ReHo and DC exhibited the lowest homogeneity. This suggested that local- and network-scale measures might not be well suited at generating homogenous parcels for depicting the complexity of cortical areal organization as a whole system.

Figure 10.

Homogeneity of individual parcels by parcel size based on 14 iFC indices. Red dots are homogeneity scores for each parcel driven from iFC similarity across individuals. Dark blue dots are homogeneity scores for each parcel based on ReHo across individuals. Each line is the LOWESS fit represents the effect of parcel size on homogeneity across individuals.

Figure 10.

Homogeneity of individual parcels by parcel size based on 14 iFC indices. Red dots are homogeneity scores for each parcel driven from iFC similarity across individuals. Dark blue dots are homogeneity scores for each parcel based on ReHo across individuals. Each line is the LOWESS fit represents the effect of parcel size on homogeneity across individuals.

Effect of Surface Geometry and Surface Registration

A potential concern regarding the methodologies presented is that the underlying surface geometry may influence the presence and locations of estimated gradients and edges. To test for this, we repeated our gradient-based analyses using surrogate fMRI data that consisted exclusively of white Gaussian noise. Figure 11A and B summarizes the spatial correlation of iFC similarity and its gradient between real and random data (ordered by participants and scan sessions, real to random). For the surrogate data, within-participant spatial correlations were again higher than between participant spatial correlations. However, spatial correlations of gradients between maps derived from real and surrogate R-fMRI showed a relatively low degree of correlation – even for the same participant (gradient: mean r = 0.03, SD = 0.06). These findings suggest that the surface geometry explained only less than 1% of variance in the gradient of iFC similarity.

Figure 11.

Effect of surface geometry and surface registration. (A) The spatial correlations matrix of mean iFC similarity from real and random data. (B) The spatial correlations matrix of mean gradient in iFC similarity from real and random data. The correlation matrixes are ordered by participant: top left sub-matrix is from real data (first ten rows for ten sessions of participant 1, second ten rows for ten sessions of participant 2, and so on) and bottom right sub-matrix is from random data in the same order. (C) The 2D histogram depicting the relationship between cortical area and mean gradients of iFC similarity computed from random data.

Figure 11.

Effect of surface geometry and surface registration. (A) The spatial correlations matrix of mean iFC similarity from real and random data. (B) The spatial correlations matrix of mean gradient in iFC similarity from real and random data. The correlation matrixes are ordered by participant: top left sub-matrix is from real data (first ten rows for ten sessions of participant 1, second ten rows for ten sessions of participant 2, and so on) and bottom right sub-matrix is from random data in the same order. (C) The 2D histogram depicting the relationship between cortical area and mean gradients of iFC similarity computed from random data.

Thus, while inter-individual variations in the surface geometry in and of themselves can produce repeatable gradients, they are unrelated to findings observed with functional data. Although beyond the scope of the present work, preliminary analyses suggest that the strength of the findings obtained for surrogate data is largely attributable to inter-individual variations in the distance between adjacent neighbors, which impacts gradient computation (Fig. 11C). However, it is worth noting that the gradient values from noise were generally minute; for nearly all vertices, the real gradient score exceeded the 90th percentile of scores generated from noise. In those cases where this was not the case, the failure to exceed the 90th percentile was due to the presence of a relatively low gradient score in the real data (thus failing to rise above the noise) (see the supplementary Figure S11).

To address any concerns about the possible effects of noise with regard to comparison of within- and between-individual variance, we repeated our analyses focused on establishing the similarity of gradient maps across subjects and sessions; though this time, we regressed out noise-based gradient maps for each before testing the differences of within- versus between-individual variance. The correct identification rate was almost the same as 90.44% for real data and 90.96% after regressed out the noise-based gradient. In sum, the whole transition patterns were predominantly driven brain function, rather than structure.

As a final sanity check, we repeated our edge detection analyses using the similarity of iFC maps rather than gradient maps. Spatial correlations between sessions (for the same participant) were quite low – again, reinforcing the validity of the findings obtained in our primary analyses.

Discussion

Using publically available datasets, the present work assessed the ability of gradient-based iFC boundary mapping approaches (Cohen et al. 2008; Wig et al. 2013, 2014; Gordon et al. 2014, 2015; Laumann et al. 2015) to reliably characterize inter-individual variation in the transitional properties of cortical areas throughout the brain. First, from the perspective of “finger-printing”, across the 10 scan sessions, we were able to demonstrate the ability to accurately identify R-fMRI data obtained from a given individual based upon the similarity of full brain gradient and full-brain edge maps derived from each of the scan sessions. Supplementary analyses demonstrated that these abilities were not driven by structural contributions. Second, using image intraclass correlation coefficient (I2C2), we were able to demonstrate moderate full-brain test–retest reliability for a 10-min scan; I2C2 progressively increased with the amount of data available, reaching relatively high levels when 50 min of data were used. Third, using vertex-wise intraclass correlation coefficient, we demonstrated low to high test–retest reliability depending on region; again, as the amount of data included was increased, ICC values approached impressively high levels. Of note, within-participant variation in transition zone properties appeared to be greater in lower order networks than higher order, while between-individual variation was greater in higher order networks; these findings are consistent with those of prior efforts examining iFC profiles. Intriguingly, using a previously published twin dataset, the present work was able to push this finding one step further, showing associations between the similarity of edge maps and each, familial relations (i.e., sibling vs. non-sibling) and genetic (i.e., monozygotic vs. dizygotic twin pairs). Finally, the present work extended prior work by employing alternative iFC measures to define gradients, revealing similar but distinct gradient profiles for each of the iFC features – all of which exhibited impressive repeatability and reliability over time if sufficient data is used. Importantly, the various findings obtained in the present work did not appear to be related to potential sources of artifactual variation, such as head-motion or surface registration. As such, our findings increase confidence in the potential for gradient measures to be used as a novel feature upon which inter-individual differences in brain function can be mapped and eventually related to phenotypic variation.

Individual Areal Organizations are Unique

Long focused on the comparison of groupings of individuals that differ on one or more features (e.g., diagnostic status), the delineation of inter-individual variation has emerged as a central focus in emerging neuroscientific and clinical agendas (e.g., biomarker identification). Several studies have emphasized the potential utility of taking variations in cortical area transition zone properties into account when attempting to catalog phenotypic variation (e.g., Di Martino et al. 2009; Adelstein et al. 2011). Our analyses revealed full-brain transition patterns on an individual basis, which distinguished participants from one another across scan sessions. Findings that the topological architecture differs between participants are consistent with recent work highlighting inter-individual variation in key features of the topological architecture (e.g. area size and shape) (Gordon et al. 2015). Perhaps more interesting, were findings that full-brain gradient and edge maps can distinguish nearly all individuals using only 10 min of data (Fig. 1C). This finding echoes recent studies, which suggested the potential for full-brain characterizations of the intrinsic brain to “fingerprint” individuals (Poldrack et al. 2013; Miranda-Dominguez et al. 2014; Finn et al. 2015) – a capability that is encouraging for efforts towards biomarker discovery. Although not a focus of the present work, confirmatory analyses replacing actual functional data with random noise were able to demonstrate the fingerprinting abilities of structural brain properties embedded in the cortical surface as well. Importantly, these structural fingerprints were largely unrelated to the various indices derived from true functional data. Taken together with prior work, it appears that future work may find potential value in the development of fingerprinting profiles.

Individual Functional Areal Organizations Are Reliable With Sufficient Scan Time

As suggested by prior work, more is generally better when discussing the test–retest reliability for R-fMRI findings. While fingerprinting could be performed with relatively high accuracy using only 10 min of data, the bar for drawing reliable functional boundaries appears to be higher. Recent work using a highly sampled individual dataset found that the reproducibility of the connectivity matrix improved steeply as the amount of data available increased from 9 to 27 min, beyond which more modest increases were observed (Laumann et al. 2015). Consistent with this report, we found that 20 min or more of data were required for gradient and edge density scores to achieve moderate to high test–retest reliabilities at the majority of vertices.

A key caution is that while studies attempting to determine the minimum sufficient data required to optimize an analysis commonly report their findings in terms of time, most do not allow this construct to be disentangled from the number of time points included. In our findings, although the number of time points for eNKI-TRT dataset (TR = 645 ms, 10 min, multiband sequence) was about three times that of the HNU dataset (TR = 2000 ms, 10 min, standard sequence), only the vertex-wise reliabilities showed an advantage over the 10-min HNU dataset, and they were inferior to those obtained with 20 min of HNU data. This may suggest the amount of time sampled could be more important than the number of time points. However, it is important to note that the eNKI-TRT dataset used an early version of the CMRR multiband sequence, rather than standard EPI as in the HNU dataset. Future work will be required to disambiguate the various factors further. It is also possible that continued refinement of preprocessing strategies may help to further decrease data needs for such analyses. Nonetheless, we emphasize that at the present time, it appears prudent to obtain more data rather than less if the goal of a study is to reliably characterize fine-grained functional boundaries.

Functional Boundaries Have Familial and Genetic Associations

Recent work has begun to explore genetic and environmental influences on inter-individual variation in indices of intrinsic brain function (Glahn et al. 2010; Schutte et al. 2013; Fu et al. 2015; Yang et al. 2016). Already, connectivity within a number of functional networks has been shown to have moderate to high heritability (e.g. default, frontoparietal, somatosensory, visual, and attention networks). Complementing these findings are those of studies suggesting the heritability of some portions of between-network functional connectivity (e.g., Yang et al. 2016), as well as larger global network architecture properties (e.g. modularity, clustering coefficient, global efficiency) (van den Heuvel et al. 2013; Sinclair et al. 2015). The present work helps to build upon this growing body of literature by suggesting possible familial (sibling vs. non-sibling) and genetic (monozygotic vs. dizygotic) associations for measures of areal brain organization. Importantly, these associations did not appear to be driven by the underlying neuroanatomy. Consistent with previous studies suggesting that different genetic mechanisms are responsible for the structure and function of brain areal organization (e.g., Glahn et al. 2010). Of note, familial effects (sibling vs. non-sibling) were found to be relatively robust, possibly suggesting environmental contributions in addition to genetic effects, though more specific constraints would need to be built into the design to optimally establish them –e.g. monozygotic twins raised apart vs. monozygotic twins raised together, or non-sibling raised together vs. non-siblings raised apart.

In light of the relatively limited amount of R-fMRI data available per participant in the QTIM dataset (scan duration = 5 min), the results obtained were taken to be particularly encouraging. We considered pursuing quantification of heritability for the indices of areal organization examined in the present work, but we decided against it as: (1) comprehensive exploration of the topic was beyond the scope of the present work, and (2) the limited data available for each participant is still likely to compromise the quality of gradients and edge maps obtained, thereby leading to an underestimation. Future work with larger datasets focused on the collection of large amounts of R-fMRI data from each participant (e.g., the Human Connectome Project, Van Essen et al. 2013) would be more appropriate for such determinations.

Sources of Within- and Between-individual Variability in Areal Organization

A growing number of studies are appreciating regional and network-level differences in the stability of connectivity patterns across individuals and time, which can be informative to our understanding of brain development and function. Intra-individual variability characterizes the temporal stability of transitional zone properties for cortical areas (i.e., low variation infers high stability), while inter-individual variability captures the conservation of these zones from one person to the next (i.e., low variation infers conservation). Intra- and inter-individual variability of transitions in iFC features were not uniformly distributed across the cortex in the present work. Specifically, multimodal association networks (e.g., default, dorsal attention, and executive control) exhibited greater variability between individuals than within; in contrast, unimodal networks (e.g., visual, sensorimotor) were found to have lower inter-individual variability whilst high intra-individual variability. These results echoed previous findings highlighting greater intra-individual variation in lower-order networks (e.g. Mennes et al. 2010; Craddock et al. 2013), whose functional interactions are more heavily influenced by current task demands (Mennes et al. 2010). Similarly, they echo the findings of studies suggesting greater inter-individual variation in higher order networks (Mueller et al. 2013; Wang and Liu 2014; Langs et al. 2015; Laumann et al. 2015), which appear to be more affected by genetic and environmental factors (Anderson and Finlay 2014; Gao et al. 2014). The greater susceptibility of higher order multimodal association networks to environmental influences is not surprising given their more protracted developmental period relative to unimodal (Mueller et al. 2013; Zilles and Amunts 2013). Additionally, the later evolutionary development and enlargement of association cortices supporting higher-order network function may contribute in part to our findings of increased variation (i.e., decreased conservation) across individuals (Van Essen and Dierker 2007; Brun et al. 2009; Van Essen et al. 2012a; Chan et al. 2014).

Potential confounds that can arise in consideration of inter-individual variation in functional transition zones come from differences in cortical folding patterns (Hill et al. 2010). For example, sulcal depth can exhibit patterns of inter-individual variation that are similar to those observed in iFC (Mueller et al. 2013). Additionally, given that current surface-registration algorithms are based on the anatomical curvature, the cross-individual registration could lead to non-uniform misalignment across cortex, thereby contributing to individual variability into functional transitions (Robinson et al. 2014). In the present study, we explored these potential confounds through the replacement of functional MRI data with random noise. The underlying anatomical architecture had only a slight relationship on functional transition profile of iFC features (Fig. 8B), reaffirming confidence that our findings are driven by iFC.

Different Functional Features Shows Similar but Distinct Areal Transition Profiles

Recent years have witnessed the emergence of a growing number of R-fMRI measures, each capturing a unique aspect of the intrinsic functional architecture. The distinctions among some indices can be readily delineated based upon differences in their definitions (e.g., centrality measures) or the networks being examined (e.g., dual regression of network components); while for others, the exact positioning of one measure relative to another can be more challenging (e.g., amplitude of low frequency fluctuation [ALFF], voxel-mirrored homotopic connectivity [VMHC], ReHo). While the iFC similarity measure employed for initial gradient-based mapping efforts has several desirable features (e.g., ability to overcome noise through averaging across vertices, utility for full-brain parcellation, most homogeneous parcels), it is only one aspect of the intrinsic brain. It can be argued that studies may benefit from selection of their measures based upon the specific question or purpose at hand. For example, if the question is specifically about the transitions of frontoparietal network, then dual-regression may be preferable. As might be expected, for DR-networks, the ICC of a given voxel showed profound dependencies on network membership and gradient strength; such properties should be considered in their application. Alternatively, more exploratory approaches may derive benefit from the calculation of multiple features for each cortical area simultaneously (e.g., fALFF, ReHo, VMHC, DC), providing gradient profiles which can be employed for analyses.

Among the functional indices examined in the present work, the gradients derived from ReHo (an index of local synchrony) were the most reliable, followed by similarity of iFC. Relative to other measures, these two measures are unique in that they involve averaged values, which will provide some degree of noise reduction. While the number of values being averaged into the ReHo score for a given vertex are dramatically fewer than those that go into the same vertex‘s similarity of iFC, the values being averaged are inherently more similar to one another. Additionally, it is worth noting that ReHo is consistently noted to engender a high degree of reliability as a measure, which may in turn lend to its more reliable gradients (Jiang and Zuo 2015; Chen et al. 2015). It may be worth noting that ReHo relies on the non-parametric Kendall‘s Concordance Coefficient rather than Pearson‘s Correlation, which may also be a contributing factor (Zuo et al. 2013). Regarding similarity of iFC, two observations are worth noting. First, that the most closely related function indices were the dual regression components for the default and sensorimotor networks, which are most similar to the task-negative and task-positive networks defined by Fox et al. (2005). Second, those gradients found to be highly reliable for similarity of iFC, tended to have reliability across the other R-fMRI indices as well. These findings suggest summary measures (e.g., similarity of iFC) have clear strengths, though may have limitations with respect to their reliability to findings outside the cardinal large-scale networks (e.g., task-positive, task-negative).

Limitations and Future Directions

Although promising, there are several limitations to the findings of the present work beyond those already discussed, which merit consideration. First, R-fMRI scans did not cover the whole brain for all the participants, particularly in the inferior temporal and orbital frontal cortex; as such, our ability to map transition profiles for these areas was limited. This reduced coverage could have impacted our findings (e.g., DR-auditory network), though we do not believe to a high degree. Second, while several studies have demonstrated that the areal transition profile based on R-fMRI correspond to the tasks activation (Wig et al. 2013; Laumann et al. 2015), the present work is limited in its ability to make generalizations regarding transition zone repeatability beyond R-fMRI; future serial scanning efforts would benefit from the inclusion of task fMRI as well. If the areal transition profile from task and rest are reliable, it may facilitate alignment across participants by offering spatial functional variability information to registration progress (Robinson et al. 2014). Finally, the current findings of reliable transition profile are based on datasets sampling a one-month period; future work may investigate the dynamic changes of areal organization across more extended periods of time, as well as address questions regarding potential age-related differences across the lifespan.

Conclusions

In sum, the present work demonstrates the ability to map repeatable, individual-specific areal transition profiles for an array of iFC features, confirming their potential for fingerprinting and biomarker discovery. Data needs for achieving moderate to high test–retest reliability appeared to be greater than those for fingerprinting, though achievable. Points of convergence were noted in the cortical area transition zones defined using differing iFC indices. However, the distinctiveness of the full-brain transition zones profiles obtained using differing iFC indices suggests the merits of considering multiple indices to provide a more comprehensive characterization of cortical area transition zone properties.

Supplementary Material

Supplementary material can be found at: http://www.cercor.oxfordjournals.org/

Funding

This work was supported by the National Basic Research (973) Program (2015CB351702); Funds of International Exchange for Post-doctoral Scientists of China (2013), grants from the National Institiute of Health (NIMH U01MH099059 to M.P.M.; NIMH BRAINS R01-MH101555 to R.C.C.), the Natural Science Foundation of China (81270023, 81278412, 81171409, 81000583, 81471740, 81220108014), the National Science Foundation for Post-doctoral Scientists of China (No. 2013M530073), the Key Research Program and the Hundred Talents Program of the Chinese Academy of Sciences (KSZD-EW-TZ-002, X.N.Z), as well as gifts from Phyllis Green, Randolph Cowen and Joseph Healey to M.P.M.

Notes

We would also like to thank Ian Hickie, Katie McMahon, Paul Thompson, and Greig de Zubicaray, who are investigators for the QTIM project and generously provided the data (Australian Twin Study; R01HD050735). We also thank Evan Gordon and Timothy Laumann for valuable discussions on surface-based gradient computation. Conflict of Interest: All of the authors declare no conflict of interest.

1
We made use of  Yeo et al. (2011) network definitions rather than those of Smith et al. (2009) networks due to: (a) the lack of confidence maps, and (b) overlaps that occur among the network maps when transformed into surface space.

References

Adelstein
JS
,
Shehzad
Z
,
Mennes
M
,
Deyoung
CG
,
Zuo
X-N
,
Kelly
C
,
Margulies
DS
,
Bloomfield
A
,
Gray
JR
,
Castellanos
FX
, et al
.
2011
.
Personality is reflected in the brain's intrinsic functional architecture
.
PLoS One
 .
6
:
e27633
.
Amunts
K
,
Zilles
K
.
2012
.
Architecture and organizational principles of Broca's region
.
Trends Cogn Sci
 .
16
:
418
426
.
Anderson
ML
,
Finlay
BL
.
2014
.
Allocating structure to function: the strong links between neuroplasticity and natural selection
.
Front Hum Neurosci
 .
7
:
918
.
Bellec
P
,
Perlbarg
V
,
Jbabdi
S
.
2006
.
Identification of large-scale networks in the brain using fMRI
.
Neuroimage
 .
29
:
1231
1243
.
Bellec
P
,
Rosa-Neto
P
,
Lyttelton
OC
,
Benali
H
,
Evans
AC
.
2010
.
Multi-level bootstrap analysis of stable clusters in resting-state fMRI
.
Neuroimage
 .
51
:
1126
1139
.
Biswal
B
,
Yetkin
FZ
,
Haughton
VM
,
Hyde
JS
.
1995
.
Functional connectivity in the motor cortex of resting human brain using echo-planar MRI
.
Magn Reson Med
 .
34
:
537
541
.
Biswal
BB
,
Mennes
M
,
Zuo
X-N
,
Gohel
S
,
Kelly
C
,
Smith
SM
,
Beckmann
CF
,
Adelstein
JS
,
Buckner
RL
,
Colcombe
S
, et al
.
2010
.
Toward discovery science of human brain function
.
Proc Natl Acad Sci USA
 .
107
:
4734
4739
.
Brodmann
K
.
1909
. Brodmann's localisation in the cerebral cortex. Springer. Leipzig.
Brun
CC
,
Leporé
N
,
Pennec
X
,
Lee
AD
,
Barysheva
M
,
Madsen
SK
,
Avedissian
C
,
Chou
Y-Y
,
de Zubicaray
GI
,
McMahon
KL
, et al
.
2009
.
Mapping the regional influence of genetics on brain structure variability – a tensor-based morphometry study
.
Neuroimage
 .
48
:
37
49
.
Buckner
RL
,
Yeo
BTT
.
2014
.
Borders, map clusters, and supra-areal organization in visual cortex
.
Neuroimage
 .
93
:
292
297
.
Bullmore
E
,
Sporns
O
.
2009
.
Complex brain networks: graph theoretical analysis of structural and functional systems
.
Nat Rev Neurosci
 .
10
:
186
198
.
Canny
J
.
1986
.
A computational approach to edge detection
.
IEEE Trans Pattern Anal Mach Intell
 .
8
:
679
698
.
Chan
MY
,
Park
DC
,
Savalia
NK
,
Petersen
SE
,
Wig
GS
.
2014
.
Decreased segregation of brain systems across the healthy adult lifespan
.
Proc Natl Acad Sci USA
 .
111
:
e4997
e5006
.
Chen
B
,
Xu
T
,
Zhou
C
,
Wang
L
,
Yang
N
,
Wang
Z
,
Dong
H-M
,
Yang
Z
,
Zang
Y-F
,
Zuo
X-N
, et al
.
2015
.
Individual variability and test–retest reliability revealed by ten repeated resting-state brain scans over one month
.
PLoS One
 .
10
:
e0144963
.
Chen
C-H
,
Fiecas
M
,
Gutiérrez
ED
,
Panizzon
MS
,
Eyler
LT
,
Vuoksimaa
E
,
Thompson
WK
,
Fennema-Notestine
C
,
Hagler
DJ
,
Jernigan
TL
, et al
.
2013
.
Genetic topography of brain morphology
.
Proc Natl Acad Sci USA
 .
110
:
17089
17094
.
Cohen
AL
,
Fair
DA
,
Dosenbach
NUF
,
Miezin
FM
,
Dierker
D
,
Van Essen
DC
,
Schlaggar
BL
,
Petersen
SE
.
2008
.
Defining functional areas in individual human brains using resting functional connectivity MRI
.
Neuroimage
 .
41
:
45
57
.
Cordes
D
,
Haughton
VM
,
Arfanakis
K
,
Carew
JD
,
Turski
PA
,
Moritz
CH
,
Quigley
MA
,
Meyerand
ME
.
2001
.
Frequencies contributing to functional connectivity in the cerebral cortex in “Resting-state” data
.
AJNR Am J Neuroradiol
 .
22
:
1326
1333
.
Cox
RW
.
1996
.
AFNI: software for analysis and visualization of functional magnetic resonance neuroimages
.
Comput Biomed Res
 .
29
:
162
173
.
Craddock
RC
,
Milham
MP
,
LaConte
SM
.
2013
.
Predicting intrinsic brain activity
.
Neuroimage
 .
82
:
127
136
.
Dale
AM
,
Fischl
B
,
Sereno
MI
.
1999
.
Cortical surface-based analysis. I. Segmentation and surface reconstruction
.
Neuroimage
 .
9
:
179
194
.
de Zubicaray
GI
,
Chiang
M-C
,
McMahon
KL
,
Shattuck
DW
,
Toga
AW
,
Martin
NG
,
Wright
MJ
,
Thompson
PM
.
2008
.
Meeting the challenges of neuroimaging genetics
.
Brain Imaging Behav
 .
2
:
258
263
.
Di Martino
A
,
Shehzad
Z
,
Kelly
C
,
Roy
AK
,
Gee
DG
,
Uddin
LQ
,
Gotimer
K
,
Klein
DF
,
Castellanos
FX
,
Milham
MP
.
2009
.
Relationship between cingulo-insular functional connectivity and autistic traits in neurotypical adults
.
Am J Psychiatry
 .
166
:
891
899
.
Docherty
AR
,
Sawyers
CK
,
Panizzon
MS
,
Neale
MC
,
Eyler
LT
,
Fennema-Notestine
C
,
Franz
CE
,
Chen
C-H
,
McEvoy
LK
,
Verhulst
B
, et al
.
2015
.
Genetic network properties of the human cortex based on regional thickness and surface area measures
.
Front Hum Neurosci
 .
9
:
440
.
Felleman
DJ
,
Van Essen
DC
.
1991
.
Distributed hierarchical processing in the primate cerebral cortex
.
Cereb Cortex
 .
1
:
1
47
.
Finn
ES
,
Shen
X
,
Scheinost
D
,
Rosenberg
MD
,
Huang
J
,
Chun
MM
,
Papademetris
X
,
Constable
RT
.
2015
.
Functional connectome fingerprinting: identifying individuals using patterns of brain connectivity
.
Nat Neurosci
 .
1
11
.
Fischl
B
,
Sereno
MI
,
Dale
AM
.
1999
a.
Cortical surface-based analysis. II: inflation, flattening, and a surface-based coordinate system
.
Neuroimage
 .
9
:
195
207
.
Fischl
B
,
Sereno
MI
,
Tootell
RB
,
Dale
A M
.
1999
b.
High-resolution intersubject averaging and a coordinate system for the cortical surface
.
Hum Brain Mapp
 .
8
:
272
284
.
Fox
MD
,
Snyder
AZ
,
Vincent
JL
,
Corbetta
M
,
Van Essen
DC
,
Raichle
ME
.
2005
.
The human brain is intrinsically organized into dynamic, anticorrelated functional networks
.
Proc Natl Acad Sci USA
 .
102
:
9673
9678
.
Friston
KJ
,
Williams
S
,
Howard
R
,
Frackowiak
RS
,
Turner
R
.
1996
.
Movement-related effects in fMRI time-series
.
Magn Reson Med
 .
35
:
346
355
.
Fu
Y
,
Ma
Z
,
Hamilton
C
,
Liang
Z
,
Hou
X
,
Ma
X
,
Hu
X
,
He
Q
,
Deng
W
,
Wang
Y
, et al
.
2015
.
Genetic influences on resting-state functional networks: a twin study
.
Hum Brain Mapp
 .
36
:
3959
3972
.
Gao
W
,
Elton
A
,
Zhu
H
,
Alcauter
S
,
Smith
JK
,
Gilmore
JH
,
Lin
W
.
2014
.
Intersubject variability of and genetic effects on the brain‘s functional connectivity during infancy
.
J Neurosci
 .
34
:
11288
11296
.
Ghosh
SS
,
Kakunoori
S
,
Augustinack
J
,
Nieto-Castanon
A
,
Kovelman
I
,
Gaab
N
,
Christodoulou
JA
,
Triantafyllou
C
,
Gabrieli
JDE
,
Fischl
B
.
2010
.
Evaluating the validity of volume-based and surface-based brain image registration for developmental cognitive neuroscience studies in children 4 to 11 years of age
.
Neuroimage
 .
53
:
85
93
.
Glahn
DC
,
Winkler
AM
,
Kochunov
P
,
Almasy
L
,
Duggirala
R
,
Carless
MA
,
Curran
JC
,
Olvera
RL
,
Laird
AR
,
Smith
SM
, et al
.
2010
.
Genetic control over the resting brain
.
Proc Natl Acad Sci USA
 .
107
:
1223
1228
.
Glatard
T
,
Lewis
LB
,
Ferreira da Silva
R
,
Adalat
R
,
Beck
N
,
Lepage
C
,
Rioux
P
,
Rousseau
M-E
,
Sherif
T
,
Deelman
E
, et al
.
2015
.
Reproducibility of neuroimaging analyses across operating systems
.
Front Neuroinform
 .
9
:
12
.
Gordon
EM
,
Laumann
TO
,
Adeyemo
B
,
Huckins
JF
,
Kelley
WM
,
Petersen
SE
.
2014
.
Generation and evaluation of a cortical area parcellation from resting-state correlations
.
Cereb Cortex
 .
26
:
288
303
.
Gordon
EM
,
Laumann
TO
,
Adeyemo
B
,
Petersen
SE
.
2015
.
Individual variability of the system-level organization of the human brain
.
Cereb cortex
 . .
Greve
DN
,
Fischl
B
.
2009
.
Accurate and robust brain image alignment using boundary-based registration
.
Neuroimage
 .
48
:
63
72
.
Hill
J
,
Inder
T
,
Neil
J
,
Dierker
D
,
Harwell
J
,
Van Essen
D
.
2010
.
Similar patterns of cortical expansion during human development and evolution
.
Proc Natl Acad Sci
 .
107
:
13135
13140
.
Hirose
S
,
Watanabe
T
,
Jimura
K
.
2012
.
Local signal time-series during rest used for areal boundary mapping in individual human brains
.
PLoS One
 .
7
:
e36496
.
Hirose
S
,
Watanabe
T
,
Wada
H
,
Imai
Y
,
Machida
T
,
Shirouzu
I
,
Miyashita
Y
,
Konishi
S
.
2013
.
Functional relevance of micromodules in the human association cortex delineated with high-resolution fMRI
.
Cereb Cortex
 .
23
:
2863
2871
.
Jenkinson
M
,
Beckmann
CF
,
Behrens
TEJ
,
Woolrich
MW
,
Smith
SM
.
2012
.
FSL
.
Neuroimage
 .
62
:
782
790
.
Jiang
L
,
Xu
T
,
He
Y
,
Hou
X-H
,
Wang
J
,
Cao
X-Y
,
Wei
G-X
,
Yang
Z
,
He
Y
,
Zuo
X-N
.
2014
.
Toward neurobiological characterization of functional homogeneity in the human cortex: regional variation, morphological association and functional covariance network organization
.
Brain Struct Funct
 .
220
:
2485
2507
.
Jiang
L-L
,
Zuo
X-N
.
2015
.
Regional homogeneity: a multimodal, multiscale neuroimaging marker of the human connectome
.
Neurosci
 . .
Kaas
JH
.
1987
.
The organization of neocortex in mammals: implications for theories of brain function
.
Annu Rev Psychol
 .
38
:
129
151
.
Kelly
C
,
Toro
R
,
Di Martino
A
,
Cox
CL
,
Bellec
P
,
Castellanos
FX
,
Milham
MP
.
2012
.
A convergent functional architecture of the insula emerges across imaging modalities
.
Neuroimage
 .
61
:
1129
1142
.
Klein
A
,
Ghosh
SS
,
Avants
B
,
Yeo
BTT
,
Fischl
B
,
Ardekani
B
,
Gee
JC
,
Mann
JJ
,
Parsey
RV
.
2010
.
Evaluation of volume-based and surface-based brain image registration methods
.
Neuroimage
 .
51
:
214
220
.
Langs
G
,
Wang
D
,
Golland
P
,
Mueller
S
,
Pan
R
,
Sabuncu
MR
,
Sun
W
,
Li
K
,
Liu
H
.
2015
.
Identifying shared brain networks in individuals by decoupling functional and anatomical variability
.
Cereb Cortex
 . .
Laumann
TO
,
Gordon
EM
,
Adeyemo
B
,
Snyder
AZ
,
Joo
SJ
,
Chen
M-Y
,
Gilmore
AW
,
McDermott
KB
,
Nelson
SM
,
Dosenbach
NUF
, et al
.
2015
.
Functional system and areal organization of a highly sampled individual human brain highlights
.
Neuron
 . .
Lund
TE
,
Madsen
KH
,
Sidaros
K
,
Luo
W-L
,
Nichols
TE
.
2006
.
Non-white noise in fMRI: does modelling have an impact?
.
Neuroimage
 .
29
:
54
66
.
McKay
DR
,
Knowles
EEM
,
Winkler
AAM
,
Sprooten
E
,
Kochunov
P
,
Olvera
RL
,
Curran
JE
,
Kent
JW
,
Carless
MA
,
Göring
HHH
, et al
.
2014
.
Influence of age, sex and genetic factors on the human brain
.
Brain Imaging Behav
 .
8
:
143
152
.
Mennes
M
,
Kelly
C
,
Zuo
X-N
,
Di Martino
A
,
Biswal
BB
,
Castellanos
FX
,
Milham
MP
.
2010
.
Inter-individual differences in resting-state functional connectivity predict task-induced BOLD activity
.
Neuroimage
 .
50
:
1690
1701
.
Miranda-Dominguez
O
,
Mills
BD
,
Grayson
D
,
Woodall
A
,
Grant
KA
,
Kroenke
CD
,
Fair
DA
.
2014
.
Bridging the gap between the human and macaque connectome: a quantitative comparison of global interspecies structure-function relationships and network topology
.
J Neurosci
 .
34
:
5552
5563
.
Mueller
S
,
Wang
D
,
Fox
MD
,
Yeo
BTT
,
Sepulcre
J
,
Sabuncu
MR
,
Shafee
R
,
Lu
J
,
Liu
H
.
2013
.
Individual variability in functional connectivity architecture of the human brain
.
Neuron
 .
77
:
586
595
.
Nelson
SM
,
Cohen
AL
,
Power
JD
,
Wig
GS
,
Miezin
FM
,
Wheeler
ME
,
Velanova
K
,
Donaldson
DI
,
Phillips
JS
,
Schlaggar
BL
, et al
.
2010
.
A parcellation scheme for human left lateral parietal cortex
.
Neuron
 .
67
:
156
170
.
Nooner
KB
,
Colcombe
SJ
,
Tobe
RH
,
Mennes
M
,
Benedict
MM
,
Moreno
AL
,
Panek
LJ
,
Brown
S
,
Zavitz
ST
,
Li
Q
, et al
.
2012
.
The NKI-Rockland sample: a model for accelerating the pace of discovery science in psychiatry
.
Front Neurosci
 .
6
:
152
.
Patriat
R
,
Molloy
EK
,
Meier
TB
,
Kirk
GR
,
Nair
VA
,
Meyerand
ME
,
Prabhakaran
V
,
Birn
RM
.
2013
.
The effect of resting condition on resting-state fMRI reliability and consistency: a comparison between resting with eyes open, closed, and fixated
.
Neuroimage
 .
78
:
463
473
.
Poldrack
RA
,
Barch
DM
,
Mitchell
JP
,
Wager
TD
,
Wagner
AD
,
Devlin
JT
,
Cumba
C
,
Koyejo
O
,
Milham
MP
.
2013
.
Toward open sharing of task-based fMRI data: the OpenfMRI project
.
Front Neuroinform
 .
7
:
12
.
Power
JD
,
Barnes
KA
,
Snyder
AZ
,
Schlaggar
BL
,
Petersen
SE
.
2012
.
Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion
.
Neuroimage
 .
59
:
2142
2154
.
Power
JD
,
Schlaggar
BL
,
Petersen
SE
.
2014
.
Recent progress and outstanding issues in motion correction in resting state fMRI
.
Neuroimage
 .
105
:
536
551
.
Robinson
EC
,
Jbabdi
S
,
Glasser
MF
,
Andersson
J
,
Burgess
GC
,
Harms
MP
,
Smith
SM
,
Van Essen
DC
,
Jenkinson
M
.
2014
.
MSM: a new flexible framework for multimodal surface matching
.
Neuroimage
 .
100
:
414
426
.
Rubinov
M
,
Sporns
O
.
2010
.
Complex network measures of brain connectivity: uses and interpretations
.
Neuroimage
 .
52
:
1059
1069
.
Schmitt
JE
,
Eyler
LT
,
Giedd
JN
,
Kremen
WS
,
Kendler
KS
,
Neale
MC
.
2007
.
Review of twin and family studies on neuroanatomic phenotypes and typical neurodevelopment
.
Twin Res Hum Genet
 .
10
:
683
694
.
Schutte
NM
,
Hansell
NK
,
de Geus
EJC
,
Martin
NG
,
Wright
MJ
,
Smit
DJA
.
2013
.
Heritability of resting state EEG functional connectivity patterns
.
Twin Res Hum Genet
 .
16
:
962
969
.
Shou
H
,
Eloyan
A
,
Lee
S
,
Zipunnikov
V
,
Crainiceanu
AN
,
Nebel
NB
,
Caffo
B
,
Lindquist
MA
,
Crainiceanu
CM
.
2013
.
Quantifying the reliability of image replication studies: the image intraclass correlation coefficient (I2C2)
.
Cogn Affect Behav Neurosci
 .
13
:
714
724
.
Sincich
LC
,
Horton
JC
.
2005
.
The circuitry of V1 and V2: integration of color, form, and motion
.
Annu Rev Neurosci
 .
28
:
303
326
.
Sinclair
B
,
Hansell
NK
,
Blokland
GAM
,
Martin
NG
,
Thompson
PM
,
Breakspear
M
,
de Zubicaray
GI
,
Wright
MJ
,
McMahon
KL
.
2015
.
Heritability of the network architecture of intrinsic brain functional connectivity
.
Neuroimage
 .
121
:
243
252
.
Smith
SM
,
Fox
PT
,
Miller
KL
,
Glahn
DC
,
Fox
PM
,
Mackay
CE
,
Filippini
N
,
Watkins
KE
,
Toro
R
,
Laird
AR
,
Beckmann
CF
.
2009
.
Correspondence of the brain‘s functional architecture during activation and rest
.
Proc Natl Acad Sci
 .
106
:
13040
13045
.
Smith
SM
,
Jenkinson
M
,
Woolrich
MW
,
Beckmann
CF
,
Behrens
TEJ
,
Johansen-Berg
H
,
Bannister
PR
,
De Luca
M
,
Drobnjak
I
,
Flitney
DE
, et al
.
2004
.
Advances in functional and structural MR image analysis and implementation as FSL
.
Neuroimage
 .
23
(
Suppl 1
):
S208
S219
.
van den Heuvel
MP
,
van Soelen
ILC
,
Stam
CJ
,
Kahn
RS
,
Boomsma
DI
,
Hulshoff Pol
HE
.
2013
.
Genetic control of functional brain network efficiency in children
.
Eur Neuropsychopharmacol
 .
23
:
19
23
.
Van Essen
DC
.
2004
.
Surface-based approaches to spatial localization and registration in primate cerebral cortex
.
Neuroimage
 .
23
(
Suppl 1
):
S97
S107
.
Van Essen
DC
,
Dierker
D
.
2007
.
On navigating the human cerebral cortex: response to “in praise of tedious anatomy”
.
Neuroimage
 .
37
:
1050
1054
.
Van Essen
DC
,
Glasser
MF
.
2014
.
In vivo architectonics: a cortico-centric perspective
.
Neuroimage
 .
93
(
Pt 2
):
157
164
.
Van Essen
DC
,
Glasser
MF
,
Dierker
DL
,
Harwell
J
.
2012
a.
Cortical parcellations of the macaque monkey analyzed on surface-based atlases
.
Cereb Cortex
 .
22
:
2227
2240
.
Van Essen
DC
,
Glasser
MF
,
Dierker
DL
,
Harwell
J
,
Coalson
T
.
2012
b.
Parcellations and hemispheric asymmetries of human cerebral cortex analyzed on surface-based atlases
.
Cereb Cortex
 .
22
:
2241
2262
.
Van Essen
DC
,
Smith
SM
,
Barch
DM
,
Behrens
TEJ
,
Yacoub
E
,
Ugurbil
K
.
2013
.
The WU-Minn human connectome project: an overview
.
Neuroimage
 .
80
:
62
79
.
Vogt
C
,
Vogt
O
.
1919
.
Allgemeinere Ergebnisse unserer Hirnforschung
.
J Psychol Neurol
 .
292
398
.
Wandell
B a
,
Winawer
J
.
2011
.
Imaging retinotopic maps in the human brain
.
Vision Res
 .
51
:
718
737
.
Wang
D
,
Liu
H
.
2014
.
Functional connectivity architecture of the human brain: not all the same
.
Neurosci
 .
20
:
432
438
.
Wig
GS
,
Laumann
TO
,
Cohen
AL
,
Power
JD
,
Nelson
SM
,
Glasser
MF
,
Miezin
FM
,
Snyder
AZ
,
Schlaggar
BL
,
Petersen
SE
.
2013
.
Parcellating an individual subject's cortical and subcortical brain structures using snowball sampling of resting-state correlations
.
Cereb Cortex
 .
2036
2054
.
Wig
GS
,
Laumann
TO
,
Petersen
SE
.
2014
.
An approach for parcellating human cortical areas using resting-state correlations
.
Neuroimage
 .
93
:
276
291
.
Wig
GS
,
Schlaggar
BL
,
Petersen
SE
.
2011
.
Concepts and principles in the analysis of brain networks
.
Ann N Y Acad Sci
 .
1224
:
126
146
.
Xing
X-X
,
Zhou
Y-L
,
Adelstein
JS
,
Zuo
X-N
.
2011
.
PDE-based spatial smoothing: a practical demonstration of impacts on MRI brain extraction, tissue segmentation and registration
.
Magn Reson Imaging
 .
29
:
731
738
.
Xu
T
,
Yang
Z
,
Jiang
L
,
Xing
X-X
,
Zuo
X-N
.
2015
.
A Connectome Computation System for discovery science of brain
.
Sci Bull
 .
60
:
86
95
.
Yan
C-G
,
Cheung
B
,
Kelly
C
,
Colcombe
S
,
Craddock
RC
,
Di Martino
A
,
Li
Q
,
Zuo
X-N
,
Castellanos
FX
,
Milham
MP
.
2013
.
A comprehensive assessment of regional variation in the impact of head micromovements on functional connectomics
.
Neuroimage
 .
76
:
183
201
.
Yang
Z
,
Zuo
X-N
,
McMahon
KL
,
Craddock
RC
,
Kelly
C
,
de Zubicaray
GI
,
Hickie
I
,
Bandettini
PA
,
Castellanos
FX
,
Milham
MP
, et al
.
2016
.
Genetic and environmental contributions to functional connectivity architecture of the human brain
.
Cereb Cortex
 .
26
:
2341
2352
.
Yeo
BTT
,
Krienen
FM
,
Sepulcre
J
,
Sabuncu
MR
,
Lashkari
D
,
Hollinshead
M
,
Roffman
JL
,
Smoller
JW
,
Zöllei
L
,
Polimeni
JR
, et al
.
2011
.
The organization of the human cerebral cortex estimated by intrinsic functional connectivity
.
J Neurophysiol
 .
106
:
1125
1165
.
Yeo
BTT
,
Sabuncu
MR
,
Vercauteren
T
,
Ayache
N
,
Fischl
B
,
Golland
P
.
2010
.
Spherical demons: fast diffeomorphic landmark-free surface registration
.
IEEE Trans Med Imaging
 .
29
:
650
668
.
Zang
Y-F
,
Jiang
T-Z
,
Lu
Y-L
,
He
Y
,
Tian
L-X
.
2004
.
Regional homogeneity approach to fMRI data analysis
.
Neuroimage
 .
22
:
394
400
.
Zilles
K
,
Amunts
K
.
2010
.
Centenary of Brodmann's map — conception and fate
.
Nat Rev Neurosci
 .
11
:
139
145
.
Zilles
K
,
Amunts
K
.
2013
.
Individual variability is not noise
.
Trends Cogn Sci
 .
17
:
153
155
.
Zuo
X-N
,
Anderson
JS
,
Bellec
P
,
Birn
RM
,
Biswal
BB
,
Blautzik
J
,
Breitner
JC
,
Buckner
RL
,
Calhoun
VD
,
Castellanos
FX
, et al
.
2014
.
An open science resource for establishing reliability and reproducibility in functional connectomics
.
Sci Data
 .
1
:
140049
.
Zuo
X-N
,
Ehmke
R
,
Mennes
M
,
Imperati
D
,
Castellanos
FX
,
Sporns
O
,
Milham
MP
.
2012
.
Network centrality in the human functional connectome
.
Cereb Cortex
 .
22
:
1862
1875
.
Zuo
X-N
,
Xing
X-X
.
2011
.
Effects of non-local diffusion on structural MRI preprocessing and default network mapping: statistical comparisons with isotropic/anisotropic diffusion
.
PLoS One
 .
6
:
e26703
.
Zuo
X-N
,
Xu
T
,
Jiang
L
,
Yang
Z
,
Cao
X-Y
,
He
Y
,
Zang
Y-F
,
Castellanos
FX
,
Milham
MP
.
2013
.
Toward reliable characterization of functional homogeneity in the human brain: preprocessing, scan duration, imaging resolution and computational space
.
Neuroimage
 .
65
:
374
386
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com