Evaluation of the spatial variability in the major resting-state networks across human brain functional atlases.

The human brain is intrinsically organized into resting-state networks (RSNs). Currently, several human brain functional atlases are used to define the spatial constituents of these RSNs. However, there are significant concerns about interatlas variability. In response, we undertook a quantitative comparison of the five major RSNs (default mode [DMN], salience, central executive, sensorimotor, and visual networks) across currently available brain functional atlases (n = 6) in which we demonstrated that (a) similarity between atlases was modest and positively linked to the size of the sample used to construct them; (b) across atlases, spatial overlap among major RSNs ranged between 17 and 76% (mean = 39%), which resulted in variability in their functional connectivity; (c) lower order RSNs were generally spatially conserved across atlases; (d) among higher order RSNs, the DMN was the most conserved across atlases; and (e) voxel-wise flexibility (i.e., the likelihood of a voxel to change network assignment across atlases) was high for subcortical regions and low for the sensory, motor and medial prefrontal cortices, and the precuneus. In order to facilitate RSN reproducibility in future studies, we provide a new freely available Consensual Atlas of REsting-state Networks, based on the most reliable atlases.

INTRODUCTION

Resting‐state functional magnetic resonance imaging (rs‐fMRI) has emerged as the main neuroimaging modality for the examination of the spontaneous functional architecture of the human brain. A substantial body of literature has established that spontaneous brain activity is organized into resting‐state networks (RSNs) defined by their spatiotemporal configuration and functional roles (Damoiseaux et al., 2006; Doucet et al., 2011; Fox et al., 2005; Smith et al., 2009). The spatiotemporal configuration of RSNs is based on their functional connectivity (FC) which represents the temporal correlation of the blood oxygen level‐dependent signals between their constituent brain regions (Biswal, Yetkin, Haughton, & Hyde, 1995). Functionally, RSNs can be divided into those involved in internally guided, higher order mental functions (default‐mode [DMN], central executive [CEN], and salience [SAL] networks) and those supporting externally driven, specialized sensory and motor processing (visual [VIS] and sensorimotor [SMN] networks) (Damoiseaux et al., 2006; Doucet et al., 2011; Smith et al., 2009). Examination of RSNs in healthy populations has been instrumental in identifying processes involved in brain development (Gu et al., 2015) and aging (Ferreira et al., 2016; Shaw, Schultz, Sperling, & Hedden, 2015; Siman‐Tov et al., 2016) while investigation of RSNs in clinical samples has yielded new insights in disease mechanisms (Dong, Wang, Chang, Luo, & Yao, 2018; Doucet, Moser, Luber, Leibu, & Frangou, 2018; Lee, Doucet, Leibu, & Frangou, 2018; Repovs, Csernansky, & Barch, 2011).

Despite progress, reproducibility remains one of the major concerns in neuroimaging (Poldrack et al., 2017; Turner, Paul, Miller, & Barbey, 2018). With regard to RSNs, their reproducibility can be influenced by parameters used to acquire the rs‐fMRI data (Gordon et al., 2016), as well as interindividual variability (Braga & Buckner, 2017; Gordon et al., 2017). Here, we focus on a source of interstudy variability that concerns the method used to identify RSNs. Initially, RSNs were defined as components derived from independent component analyses (ICAs) of the rs‐fMRI data series. However, variations in implementing ICA influence the number and spatial distribution of the RSNs extracted (Beckmann & Smith, 2004; Damoiseaux et al., 2006; Kiviniemi et al., 2009; Smith et al., 2009). Alternate algorithms for partitioning the rs‐fMRI data series have also resulted in variability in the interstudy spatial composition of RSNs (Franco, Mannell, Calhoun, & Mayer, 2013). It has also been suggested that the spatial variability of high‐order networks may reflect inherent functional overlap between at least some of their constituent regions (Yeo, Krienen, Chee, & Buckner, 2014). Finally, the attribution of the same label to RSNs that differ in their spatial composition represents an additional challenge. For example, the DMN defined by Shirer, Ryali, Rykhlevskaia, Menon, and Greicius (2012) includes clusters located in the medial prefrontal cortex/ventral anterior cingulate cortex (ACC), precuneus/posterior cingulate cortex, angular gyri, and medial thalamus while the DMN defined by Yeo et al. (2011) does not include any thalamic regions but instead extensively covers middle temporal and inferior frontal areas.

In response, researchers have begun to use brain functional atlases to identify normative RSNs, as this approach has the potential to harmonize RSN definition across studies. In principle, the use of a single “gold standard” functional brain atlas would greatly increase the reproducibility of results across studies, as RSNs would be defined using a standardized approach and would be independent of study sample. However, such an atlas is not currently available. Instead, multiple atlases are available which differ in the number and spatial constitution of the RSNs they define. For example, the atlas provided by Yeo et al. (2011) was based on 1,000 participants and used a clustering approach to extract seven RSNs, while the atlas provided by Shirer et al. (2012) is based on 15 individuals and defined 14 RSNs using an ICA. Critically, differences between currently available atlases have yet to be systematically investigated since variations in RSN composition are likely to influence inferences made about the functional properties of these networks and about their role in health and disease.

In this context, we undertook the present study which aims to: (a) conduct a quantitative comparison of the spatial overlap of currently used functional brain atlases and the major RSNs they define; (b) compute the regional flexibility coefficient for each anatomical region, which quantifies the number of times a particular region changes network allegiance across atlases (Bassett et al., 2011); (c) demonstrate that variation in the anatomical constitution of RSNs defined by different atlases results in variation in their FC, in terms of within‐ and between‐network features; and (d) generate a new consensual atlas based on the most reliable atlases, for use in future studies.

METHODS

Selection of brain functional atlases

We chose to focus on functional atlases that are freely available, are embedded in the Montreal Neurological Institute (MNI152), and are easy to use as they are provided as volumetric 3D maps. Consequently, we did not consider atlases that only provide the coordinates of the nodes of each RSN (e.g., Power et al., 2011) or surface‐based atlases that require conversion to a volumetric MNI space (Glasser et al., 2016; Ji et al., 2019).

Based on these considerations, we examined six peer‐reviewed resting‐state functional atlases. In detail, the atlases by Smith et al. (2009) (hereafter referred to as the Smith‐Atlas), by Yeo et al. (2011) (hereafter referred to as the Yeo‐Atlas), and by Shirer et al. (2012) (hereafter referred to as the Shirer‐Atlas) were chosen because they have each been cited more than 500 times by January 31, 2019 according to the Altmetric information available on each article's website. The atlas by Gordon et al. (2016) (hereafter referred to as the Gordon‐Atlas) was chosen because it is used in the Adolescent Brain Cognitive Development Study (https://abcdstudy.org/) which aims to provide national normative standards for brain FC (Casey et al., 2018). We included the atlases by Doucet et al. (2011) and Doucet, Rasgon, McEwen, Micali, and Frangou (2018) (hereafter referred to as Doucet2011‐Atlas and Doucet2018‐Atlas) that were constructed by our group using resting‐state fMRI data from two independent large samples of healthy participants (Doucet et al., 2011; Doucet, Rasgon, et al., 2018). Details of the size of the sample used in the construction of each atlas, the method of RSN extraction and the RSNs defined are provided in Table 1, Figure 1 and Table S1, Supporting Information.

Table 1

Description of the brain functional atlases considered

Reference	Abbreviated atlas name in the manuscript	Sample size (N)	Sample age (min–max)	Sample mean age	Number of networks in the atlas	Analytical approach used to construct the atlas
Smith et al. (2009)	Smith‐Atlas	36	20–35	28.5	10	ICA
Yeo et al. (2011)	Yeo‐Atlas	1,000	18–35	21.3	7	Surface vertex‐based clusteringa
Doucet et al. (2011)	Doucet2011‐Atlas	180	18–57	26	23	ICA
Shirer et al. (2012)	Shirer‐Atlas	15	18–30	Not provided	14	ICA
Gordon et al. (2016)	Gordon‐Atlas	108	18–33	21	12	Surface vertex‐based clusteringa
Doucet, Rasgon, et al. (2018)	Doucet2018‐Atlas	496	22–37	29	13	Region‐based clustering

Note. All studies used data from healthy participants only. Age in years.

Abbreviation: ICA, Independent component analysis.

While the atlas is based on a surface‐based approach, we used its volumetric version provided by the original authors. More details are in Supporting Information.

Figure 1

Spatial definition of the six brain functional atlases. Each color represents one network. Description of each atlas and each network are provided in Table 1 and Table S1, Supporting Information

Prior to any analyses, each RSN from each functional atlas was created as a binary mask and resliced in a standardized fashion (mask dimension = 91 × 109 × 91; voxel resolution = 2 × 2 × 2 mm3).

Interatlas similarity in overall network spatial overlap

We used the function multislice_pair_labeling.m from the Network Community Toolbox (http://commdetect.weebly.com), to assign a number to each RSN and ensured that networks with the highest spatial overlap across atlases were assigned the same numeric label. Then, we used the function zrand.m to compute the Rand coefficient for each pair of atlases. The Rand coefficient is an index of the overall similarity between two atlases and corresponds to the fraction of voxels that are assigned to the same network in both atlases (Traud, Kelsic, Mucha, & Porter, 2011). The Rand coefficient ranges between 0 and 1 but can be easily skewed toward 1 (Traud et al., 2011). Therefore, we used the z‐score of the Rand coefficient for each voxel i, following the formula:where SD refers to standard deviation across all voxels.

Further, for each pair of atlases, we undertook Spearman correlation analyses between the z‐scores of the Rand coefficient and the sum of size of the samples used in the two atlases considered, the interatlas difference in number of RSNs and the interatlas difference in the total number of voxels. We report significant correlations at a statistical significance of p < .05 after applying a Bonferroni correction.

Interatlas similarity in the spatial composition of the major RSNs

For these analyses, we examined the DMN, the CEN, the SAL, the SMN, and the VIS as these are considered as the major RSNs (Figure S1, Supporting Information). We note that in some atlases, RSNs that comprise identical anatomical regions are described using different terms. Specifically, the SAL is referred to as the ventral attention network in the Yeo‐Atlas and as the executive central network in the Smith‐Atlas despite including the same anatomical regions in both these atlases. This partially reflects the historical evolution of the labeling of RSNs as the SAL network had not been fully established as separate from the executive central network when the Smith‐Atlas was published. In this study, the term SAL is used to refer to the network that includes the insula and the dorsal ACC as defined by the authors that established its presence and its function (Menon & Uddin, 2010; Seeley et al., 2007). Other labels seem to reflect authors' preferences that vary between naming a network based on its anatomical composition or based on its presumed functionality. Consequently, the CEN is referred to as the lateral frontoparietal network in the Doucet2011‐Atlas, in the Gordon‐Atlas, in the Smith‐Atlas and in the Yeo‐Atlas. Moreover, some atlases provide separate maps for subdivisions of the same RSN (e.g., separate maps for the left and right‐sided subdivisions of the CEN as in the Shirer‐, the Doucet2011‐, and the Smith‐Atlases). Such subnetworks were combined to obtain a single map for the entire RSN (Table S1, Supporting Information). We were thus able to construct a map for each of the five major RSNs per atlas (Figure S1, Supporting Information).

We used the Dice's coefficient (D) (Dice, 1945) to quantify the pairwise spatial similarity of RSNs between atlases. For each pair of atlases i and j, the D coefficient for two same‐labeled RSNs, X and X , was computed as:where X and X denote the number of voxels in each network. denote the intersection of the voxels of the two same‐labeled networks from the different atlases. For each pair of atlases, a Dice's coefficient of 100% indicates perfect overlap between the two same‐labeled RSNs and a coefficient of 0% denotes no overlap.

Voxel‐wise variability between the RSNs across atlases

We used the flexibility function of the Network Community Toolbox (http://commdetect.weebly.com) to compute a flexibility coefficient of each voxel (Bassett et al., 2011) as the ratio of the number of times a particular voxel changed network allegiance across atlases to the total number of possible changes (i.e., number of atlases‐1). To account for variation in the degree of total brain coverage in each atlas (Figure S2a and Table S2, Supporting Information), the flexibility coefficient of each voxel was weighted by the number of atlases that included that voxel. Low flexibility coefficients indicate that a specific voxel is likely to retain its network allegiance across atlases, while high flexibility coefficients indicate that a specific voxel is likely to change network allegiance across atlases.

We estimated the flexibility coefficient of each voxel based (a) on the RSNs as originally reported in each atlas and (b) for the five major RSNs as defined in the previous section.

Quantification of atlas‐dependent variation in RSN FC

We hypothesized that within‐ and between‐network FC (WNC and BNC, respectively) would vary as a function of the variable anatomical constitution of the RSNs in the different atlases. To demonstrate this we used publicly available rs‐fMRI datasets from a sample of 100 unrelated healthy participants (58 women) of the Human Connectome Project (HCP; http://www.humanconnectomeproject.org/) (Glasser et al., 2013). The mean and SD of the age of this sample was 29.1 (3.5) years. The acquisition parameters and preprocessing procedures are detailed in the Supporting Information. For each HCP participant, we computed the WNC and BNC for the DMN, CEN, SAL, SMN, and VIS in each of the six atlases. The WNC was estimated by averaging all possible pairwise Pearson's correlations between the time series of the network's voxels. The BNC was estimated by averaging the voxel‐wise correlations between each network's time series and the time series of the other networks. All correlation coefficients were transformed into Fisher's z‐score prior to further analyses.

For each pair of atlases i and j, the functional similarity coefficient of the WNC of two RSNs, X and X , was computed as:where cov is the covariance, WNC and WNC represent the WNC for networks X and X, and and are the corresponding SD. The similarity coefficient of the WNC can range between −1 and 1, respectively, indicating that the WNC of network X and X may show total inverse or total positive covariation while a value of 0 would indicate no covariation. After substituting WNC for BNC, we applied the same equation to calculate the functional similarity coefficient for the BNC of each pair X and X networks.

We used Spearman's correlation coefficient to test the association between the Dice coefficient and the SC coefficients for WNC and BNC of each pair of the major RSNs across the six atlases.

Generation of a consensual atlas of resting‐state networks

The Consensual Atlas of REsting‐state Network (CAREN) was based on the original parcellations provided by the Yeo‐Atlas, the Gordon‐Atlas, the Doucet2011‐Atlas, and the Doucet2018‐Atlas. We selected these atlases on the basis on their substantial overlap as discussed further below and considered only voxels (N = 39,261) included in all four atlases (Figure S2b, Supporting Information).

For each pair of voxels, we computed their binary probability of belonging to the same network in each of the original four atlases. We then averaged these probabilities across atlases to generate a 39,261 × 39,261 P‐matrix. We then applied a hierarchical clustering to iteratively link pairs of voxels with the highest probability of belonging to the same network, forming progressively larger clusters in a hierarchical tree.

The optimal clustering solution, which formed the consensual networks in CAREN, was based on two considerations. First, using the P‐matrix, we conservatively identified clusters for which the average within‐cluster probability was at least 20 times greater than the average between‐cluster probability (Drysdale et al., 2017) (Figure S5, Supporting Information). Second, in accordance with Yeo et al. (2011), we computed a voxel‐wise confidence value to quantify the similarity between the CAREN networks (i.e., the networks that emerged through hierarchical clustering) and the same‐labeled networks in the four original atlases. This voxel‐wise index ranged from 0 (the network assignment of a voxel in CAREN differed from its network assignment in all four atlases) to 100% (the network assignment of a voxel in CAREN was identical to its network assignment in all four atlases). In the main manuscript, we present the optimal partition while details of plausible partitions at lower hierarchical levels are provided in Supporting Information (Section 4.3, Figure S7, Supporting Information).

The networks from CAREN and from the other atlases were visualized using the BrainNet Viewer v1.61 software (http://www.nitrc.org/projects/bnv/) (Xia, Wang, & He, 2013).

RESULTS

Interatlas similarity in overall network spatial overlap

The similarity between atlases in overall network spatial overlap, as inferred by the z‐score of the Rand coefficient, showed no significant association with the interatlas difference in the number of networks or total number of voxels (p > .05) but was positively associated with the size of the samples used in their construction (ρ = 0.72, p = 3 × 10−4). Accordingly, the Doucet2018‐Atlas and the Yeo‐Atlas, which are based on the largest samples, showed the highest interatlas overlap while the Shirer‐Atlas and the Smith‐Atlas which are based on the smallest sample sizes showed the lowest overlap (Figure 2). The degree of similarity between atlases could not be attributed to the method of network extraction employed in each original atlas; the Yeo‐Atlas and the Doucet2018‐Atlas, respectively, used surface‐ and volume‐based approaches yet they were very similar. By contrast, the Shirer‐Atlas and the Smith‐Atlas that were the least similar atlases were both volumetric, based on ICA.

Figure 2

Interatlas similarity in overall network spatial overlap. Higher scores indicate higher similarity

Interatlas similarity in the spatial composition of the major RSNs

Among all major RSNs, the interatlas pairwise comparisons revealed an average spatial overlap of 39% (SD: 15%), ranging from 17 to 76% between pairs of same‐labeled networks (Table S3 and Figure S3, Supporting Information). Only 16 interatlas same‐labeled networks (i.e., 21% of all possibilities) showed a spatial overlap of at least 50% in their spatial constituents. Figure 3 displays the spatial overlap of each major RSN across all atlases. Overall, the common overlap across the majority of the atlases (i.e., four or more) was relatively limited, ranging from 13% (for the SAL network, % of voxels assigned to the SAL network in at least four atlases) to 35% (for the VIS network).

Figure 3

Spatial overlap of the five major RSNs across the six atlases. To create this map, each voxel was assigned a value ranging from 1 to 6 reflecting the degree of constancy in the allocation of this voxel to the same RSN across atlases; a value of 1 signifies that a voxel was assigned to a network in 1 atlas only while a value of 6 signifies that a voxel was assigned to the same network in all six atlases. CEN, central executive network; DMN, default mode network; RSN, resting‐state networks; SAL, salience network; SMN, sensorimotor network; VIS, visual network

VIS network

Among the major RSNs, the VIS networks had the highest spatial similarity in all pairwise comparisons across the six atlases (D range: 17–76%, mean [SD]: 50% [20%]) (Table S3, Supporting Information). The largest overlap was noted between the VIS networks defined in the Doucet2018‐Atlas and the Yeo‐Atlas, while the lowest overlap was observed in VIS networks defined by the Shirer‐ and Smith‐Atlases (Figure S3a, Supporting Information). The regions showing the greatest constancy in being allocated to the VIS network across four or more atlases were located in the medial (calcarine, lingual gyrus) and lateral (inferior/middle) occipital lobe (Figure 3a).

Sensorimotor network

The SMN showed moderate similarity in all pairwise comparisons across the six atlases (D range: 19–71%, mean [SD]: 41% [17%], Table S3, Supporting Information). The largest overlap was noted between the SMN defined by the Doucet2018‐ and the Doucet2011‐Atlases, while the lowest overlap was observed between the SMN defined by the Shirer‐ and Smith‐Atlases (Figure S3b, Supporting Information). The regions showing the greatest constancy in being allocated to the SMN in four or more atlases were located in the ventral part of the pre/postcentral gyri and superior temporal gyri (Figure 3b).

Default‐mode network

The networks defined as DMN across all atlases also showed moderate similarity in their spatial overlap (D range: 22–62%, mean [SD]: 38% [10%]) (Table S3, Supporting Information). The DMN as defined by the Doucet2011‐Atlas and the Shirer‐Atlas showed the largest overlap, while the DMN as defined in the Doucet2018‐Atlas and Smith‐Atlas showed the smallest overlap (Figure S3c, Supporting Information). The regions showing the greatest constancy in being allocated to the DMN in four or more atlases were located in the posterior cingulate cortex/precuneus, the ventral ACC/medial prefrontal cortex, and in the angular gyri (Figure 3c).

Salience network

The SAL network also showed substantial pairwise similarity across atlases (D range: 19–46%, mean [SD]: 32% [7%], Table S3, Supporting Information). The SAL networks defined by the Yeo‐ and the Gordon‐Atlases were the most similar, while the SAL networks defined by the Shirer‐ and Smith‐Atlases were the least similar (Figure S3d, Supporting Information). The regions showing the greatest constancy in being allocated to the SAL in four or more atlases were located in the anterior insula, the dorsal ACC, and in the middle frontal cortex and supramarginal gyri (Figure 3d).

Central executive network

The spatial similarity of the networks defined as CEN across atlases was the weakest among all the major RSNs (D range: 19–50%, mean (SD): 32% (7%), Table S3, Supporting Information). The CENs defined by the Yeo‐ and the Gordon‐Atlases were the most similar, while the CENs defined by the Shirer‐ and Smith‐Atlases were the least similar (Figure S3e, Supporting Information). The regions showing the greatest constancy in being allocated to the CEN in four or more atlases were located in the middle frontal, inferior parietal, and posterior middle temporal cortices (Figure 3e).

Voxel‐wise variability between the RSNs across atlases was assessed using voxel‐wise flexibility coefficient, as shown in Figure 4. Lower flexibility coefficients (i.e., high consistency in RSN assignment across atlases) were noted in the medial prefrontal regions, precuneus, dorsal ACC, precentral and postcentral gyri, and the medial occipital cortex. Higher flexibility coefficients in ventral regions reflect their inconsistent assignment to the same RSN across the six atlases, while the high flexibility coefficients in subcortical regions were attributable to the low coverage of these regions across the six atlases (as shown in Figure 2a, Supporting Information). We obtained the same results when the analyses conducted above were restricted to the five major RSNs (Figure S4, Supporting Information).

Figure 4

Spatial distribution of the voxel‐wise flexibility coefficient. Voxel‐wise flexibility coefficient quantified the constancy with which a voxel is allocated to the same resting‐state network across atlases; lower values indicate higher degree of constancy. It is defined as the ratio of the number of times a particular voxel changed network allegiance across atlases to the total number of possible changes (i.e., number of atlases‐1) weighted by the number of atlases that included that voxel

Across atlases, higher spatial overlap between the major RSNs was associated with higher similarity in their FC properties (WNC: ρ = 0.39, p = 5.2 × 10−4; BNC: ρ = 0.40, p = 3.9 × 10−4).

Consensual Atlas of REsting‐state Networks

The optimal hierarchical clustering solution (based on the probability that any two voxels across any pair of atlases were reliably part of the same network) comprised five networks that we labeled DMN, CEN, SAL, VIS, and SMN in accordance with their spatial distribution as shown in Figure 5 and as detailed in Table S4, Supporting Information. The outline of the CAREN networks over each atlas considered is shown in Figure S6, Supporting Information. Figure 5c shows the confidence value for each voxel with respect to the probability of being assigned to the same‐labeled RSN in CAREN and in each of the original four atlases. The majority of voxels (87%) had confidence values over 75%. Less than 1% of the voxels had a confidence value of 25%. The mean voxel‐wise confidence value per RSN ranged between 94% for the VIS network and 71% for the CEN.

Figure 5

CAREN. (a) Result of the hierarchical clustering. Left panel: dendrogram highlighting a five‐network partition and right panel: intervoxel matrix showing the average probability of belonging to the same network. (b) The spatial distribution of the RSNs defined in CAREN. DMN (dark blue), CEN (red), SAL (orange), SMN (cyan), and VIS network (green). Their constituent regions are described in Table S4, Supporting Information. (c) The spatial distribution of the voxel‐wise confidence value in CAREN. This measure quantifies the probability that a voxel in a CAREN network is assigned to the same‐label RSN in each of the original atlases. Black lines outline the boundaries of the CAREN networks; colors reflect the voxel‐wise confidence values (blue, low; green and orange, intermediate; red, high). CAREN, Consensual Atlas of REsting‐state Networks; CEN, central executive network; DMN, default mode network; RSN, resting‐state networks; SAL, salience network; SMN, sensorimotor network; VIS, visual network

In CAREN, the SAL emerged as the most reproducible network (orange network in Figure 5b); the average probability that voxels assigned to this network were identified across all four atlases was 66%. As defined in CAREN, the SAL comprises the anterior insula bilaterally, the dorsal ACC and posterior regions in the supramarginal gyri and anterior precuneus. The SMN and VIS networks were also reproducible; the average probability that voxels assigned to the SMN and VIS networks were identified across all four atlases was 47 and 63%, respectively (Figure S4, Supporting Information). In CAREN, the VIS network was comprised exclusively of occipital regions while the SMN comprised the sensory and motor regions (precentral and postcentral gyrus and supplementary motor area) and the primary auditory cortex (superior temporal cortex). Of note, at a lower hierarchical level (i.e., six‐network partition), the SMN is the first network to split into two with one subdivision mostly including the primary sensorimotor regions and the other comprising the auditory cortex (see more detail in Figure S6, Supporting Information). For the DMN, the average probability that voxels assigned to this network were identified across all four atlases was 48%. In CAREN, the DMN emerged as a single network comprising the medial prefrontal cortex/ventral ACC, the precuneus/PCC, the inferior frontal cortex, the angular gyri, the middle temporal cortex, and the parahippocampal gyri. The CEN emerged as the least reproducible network; the average probability that voxels assigned to the CEN were identified across all four atlases was 44%. In CAREN, the CEN comprises the dorsolateral prefrontal cortex, the superior parietal cortex, and the posterior inferior temporal cortex.

DISCUSSION

We conducted a systematic comparison of six functional brain atlases derived from independent samples of healthy individuals. We demonstrated spatial variability between RSNs defined by these atlases and the FC properties of the RSNs defined by the different atlases were influenced by their interatlas similarity. In order to promote RSN reproducibility in future studies, we constructed a CAREN, based on the most reliable atlases.

A substantial body of literature has identified a set of RSNs, namely the DMN, SAL, CEN, SMN, and VIS networks, that are reliably identified across studies and are associated with diverse mental operations (Doucet et al., 2011; Kiviniemi et al., 2009; Smith et al., 2009; Yeo et al., 2011). However, the spatial composition of these networks shows interstudy variation which also extends to the currently available functional atlases. This interstudy variation is likely related to multiple factors that include a mixture of interindividual variability (Braga & Buckner, 2017; Gordon et al., 2017), differences in parameters and processing of the rs‐fMRI data (Franco et al., 2013) and parcellation approaches (Gordon et al., 2016). It has also been suggested that brain networks, in particular the high‐order networks, may share regions which lead to their increased spatial variability (Yeo et al., 2014). When assessing the interstudy similarity of six widely used atlases, we found that the method used for the RSN extraction or the parcellation size did not seem to influence the interatlas similarity. By contrast, the most significant correlate of interatlas similarity was the number of participants contributing to the atlas construction. Therefore, this study joins recent calls advocating the importance of large sample sizes for reproducible results in the neuroimaging field (Poldrack et al., 2017; Turner et al., 2018). Our results suggest that atlases derived from samples of at least 100 individuals are conducive to greater reproducibility of atlas‐based RSN studies.

Across atlases, same‐labeled RSNs differed in their spatial distribution and consequently in their WNC and BNC. The interatlas spatial similarity in RSNs was unimpressive and averaged 39% for any of the major RSNs examined. The SMN and VIS showed the highest similarity between the six different atlases. This was not surprising as these networks are comprising mostly of sensory and motor regions, which are known to have low interindividual variability in anatomical morphology (White et al., 1997) and in resting‐state FC (Franco et al., 2013; Li et al., 2017; Mueller et al., 2013) and tend to preferentially participate in single networks (Yeo et al., 2014).

The interatlas similarity of the higher order RSNs, namely the DMN, CEN, and SAL, was generally low. This was unexpected in the case of the DMN is arguably the best characterized RSN and the only RSN that is more active at rest (Mazoyer et al., 2001; Raichle et al., 2001). However, this is in line with a recent study by Braga & Buckner (2017) who showed that these networks showed a relatively high intraindividual variability in terms of spatial definition. In the current study, the DMN showed the most substantial interatlas spatial overlap among the three networks, with an average of 38%. This is in the same range as that described by Jann et al. (2015) when comparing the spatial composition of the DMN as identified using either resting‐state fMRI or arterial spin labeling.

The interatlas variability in the spatial composition of higher order RSNs is thought to reflect an inherent feature of associative regions as these tend to participate in multiple networks (Yeo et al., 2014). However, among the associative regions, we found evidence of a medial–lateral dissociation; voxels in the medial prefrontal cortex and precuneus were highly likely to retain their network assignment across the six atlases, while the opposite was the case for voxels in the lateral frontal regions. These results are aligned with reviews (Eriksson, Vogel, Lansner, Bergstrom, & Nyberg, 2015; Nyberg & Eriksson, 2015) which have confirmed that the lateral aspects of the cortex show high interindividual variability, since these regions are crucial in implementing cognitive control (Moser et al., 2018; Niendam et al., 2012) and their variability is thought to underpin interindividual differences in executive functions (Gordon et al., 2016; Li et al., 2017; Mueller et al., 2013).

Networks labeled as CEN and SAL showed equally low interatlas overlap. The major source of interatlas difference for the CEN was variability in the assignment of lateral prefrontal cortical regions, and for the SAL variability in the assignment of the supramarginal gyrus. The inclusion of this region to the SAL is based on recent studies (Prillwitz et al., 2018; Shirer et al., 2012) that have extended the definition of this network that was initially anchored in the anterior insula and dorsal ACC (Menon & Uddin, 2010; Seeley et al., 2007). The inclusion of this posterior region in the SAL network was also confirmed during the construction of CAREN.

In the case of the DMN, the lateral/inferior frontal and temporal cortex showed the greatest variability in terms of their assignment to the DMN across the six atlases despite their recognized contribution of DMN‐related functions. Buckner, Andrews‐Hanna, and Schacter (2008) have stressed the importance of the lateral temporal cortex as a core constituent of the DMN and commented that it is poorly characterized in humans (Buckner et al., 2008). The inferior frontal region showed the lowest reproducibility across atlases with half of the atlases including it as a DMN region. It is particularly interesting since it has been more often associated with language processing (Vigneau et al., 2006) but may be important for spontaneous cognition when it occurs in the form of inner speech (Christoff, Gordon, Smallwood, Smith, & Schooler, 2009; Doucet et al., 2012).

Moving forward, we propose that the CAREN we constructed has the following advantages over existing alternatives: CAREN (a) is based on reliable atlases derived from more than 100 healthy participants each; (b) is independent of variation in neuroimaging parameters (site, acquisition sequence and analytical methods); and (c) is independent of sample composition. CAREN is composed of five networks that have high spatial similarity with RSNs defined across the main available functional atlases that were considered here. While CAREN offers a realistic option for standardizing the definition of RSNs, we acknowledge its limitations. First, CAREN was based on data from voxels that were common in all four atlases used in its construction; this led to the exclusion of cerebellar and subcortical voxels as these were not covered by all the atlases. Second, in constructing CAREN, we did not model differences in the coverage of cortical regions across the contributing original atlases. The cortical regions that were excluded from CAREN because they were not part of all contributing atlases were located in the most ventral parts of the brain. Ventral brain regions are often problematic because they are subject to artifacts associated with signal loss or inhomogeneity (Yeo et al., 2014). Therefore, we do not consider that the inclusion of a larger number of voxels from these ventral cortical regions would have improved or have influenced the parcellation leading to CAREN. Third, two of the atlases that contributed to CAREN were volumetric (Doucet2011‐ and Doucet2018‐) and two used surface‐based RSN extraction (Yeo‐ and Gordon‐). It is unlikely that this difference influenced the parcellation leading to CAREN as we used the volumetric versions of the surface‐based atlases and because our results indicate that interatlas similarly is mainly associated with sample size and not the method of RSN extraction. As more atlases using a surface‐based approach become available (e.g., Glasser et al., 2016; Ji et al., 2019) it will be possible to examine this issue further in future studies. Fourth, CAREN inherently assumes that a voxel can only be part of a single RSN which mirrors the single‐network allocation of the voxels in the atlases used to create CAREN. Although some study designs may focus on finer or dynamic partitions of the rs‐fMRI data, the purpose of CAREN is to provide a reproducible template of the major RSNs for general use.

CONCLUSION

We presented a systematic and quantitative comparison of the spatial composition and FC of the major RSNs based on their definition in six widely used functional atlases. We demonstrated a large spatial variation in RSN definition between atlases, which directly influences network‐based FC measures. The findings of this study will enable researchers to make informed decisions about their choice of functional brain atlases as templates. Furthermore, in order to bolster reproducibility of the major RSNs in future studies, we propose CAREN, a consensual atlas of RSNs that is publicly available (Supporting Information, https://www.researchgate.net/publication/334042115_Consensual_Atlas_of_REsting-state_Networks_CAREN or upon request).

APPENDIX S1 Supporting Information

Click here for additional data file.

APPENDIX S2 Supporting Information

Click here for additional data file.