Regions of Interest as nodes of dynamic functional brain networks

In 1909, Korbinian Brodmann published the results of his seminal work: maps of brain areas with different cytoarchitectures. His results were among the first to suggest that the brain does not process information as an undivided entity. Instead, cognitive tasks are distributed among specialized brain areas. Since Brodmann’s time, the neuroscientific community has reached consensus on the distributed nature of brain function (see Wig, Schlaggar, & Petersen (2011) for a review). Information processing in the brain is based on the balance between segregation and integration: there are clusters with strong internal connections and weak long-range connectivity between them (Friston, 1994; Sporns, 2013b; Tononi, Sporns, & Edelman, 1994).

Because of the crucial role of connectivity in the brain function, it is natural to model the brain as a network. In a network model of the brain, the nodes represent brain areas and the links represent the anatomical or functional connections between the nodes (Bassett & Sporns, 2017; Sporns, 2013a, 2013b; Wig et al., 2011). Network neuroscience has unveiled several important features of the structure and function of the human brain. For reviews, see, for example, Bassett and Sporns (2017), Betzel and Bassett (2017), Sporns (2013a), and Wig et al. (2011).

Networks of the brain vary across people and in time. Structural and functional brain networks have been reported to differ between people, in particular between diseased subjects and healthy controls, and to change across the lifespan (Bassett & Bullmore, 2009; Chan, Alhazmi, Park, Savalia, & Wig, 2017; Papo, Zanin, Pineda-Pardo, Boccaletti, & Buldú, 2014; Sporns, 2013b). Functional brain networks vary on shorter timescales too, for example, with different cognitive tasks (Bassett et al., 2011; Braun et al., 2015; Göttlich, Ye, Rodriguez-Fornells, Münte, & Krämer, 2017; Honey, Kötter, Breakspear, & Sporns, 2007). However, the traditional tools of connectivity analysis cannot capture this time variation: there is still a lack of appropriate methods for understanding the dynamics of brain networks.

There are two questions of fundamental importance for functional brain networks: what do the nodes represent, and how are their links defined? The common approach is to use Regions of Interest (ROIs) as the nodes. ROIs are collections of fMRI measurement voxels defined on the basis of anatomy, connectivity profiles, or function (for a review, see de Reus & Van den Heuvel, 2013). The BOLD time series of each voxel follows the changes in the voxel’s level of activity. To arrive at a time series that represents an entire ROI, its voxel time series are typically averaged. Then, the weights of the links between ROIs are quantified with some similarity measure of their time series, such as the commonly used Pearson correlation coefficient.

The ROI time series are typically taken as accurate representations of the dynamics of the voxels within the ROI. Consequently, a minimum requirement for an ROI to be reasonably defined is its functional homogeneity: each of the voxels should have similar dynamics. In our previous work (Korhonen, Saarimäki, Glerean, Sams, & Saramäki, 2017), we have used the concept of spatial consistency for quantifying this functional homogeneity. We found that spatial consistency varies widely across ROIs in the commonly used parcellations, indicating that the assumption of functional homogeneity does not hold for all ROIs in functional brain networks.

There are two possible reasons for low spatial consistency. First, it is possible that it indicates technical problems in the investigated parcellations: although functionally homogeneous regions may exist in the brain, the parcellations are not able to capture these regions. Second, spatial consistency may vary in time: averaging over periods of extremely low and moderately high consistency would yield low values of average consistency. In Korhonen et al. (2017), we speculated that the variation of spatial consistency between ROIs may not be just a technical issue that can be overcome by some sophisticated parcellation scheme. Instead, it may carry cognitive meaning and be related to changes in the ROIs’ activation, for example.

In the present work, we generalize the investigation of spatial consistency into dynamic brain networks. We explore how spatial consistency varies in time, and ask how its variation relates to changes in the local network structure around ROIs. To this end, we use two measures: spatiotemporal consistency quantifies temporal changes in spatial consistency, and network turnover measures the amount of turnover in a node’s network neighborhood across time. We use in-house data collected from 13 healthy subjects during free music listening and resting-state data of 28 healthy subjects from the Autism Brain Imaging Data Exchange (ABIDE) initiative (Di Martino et al., 2014). The in-house dataset is a subset of a larger dataset that has been earlier partially described in Alluri et al. (2015, 2017) and Burunat et al. (2015).

With these data, we show that the ROIs exhibit varying levels of spatiotemporal consistency, which indicates that their spatial consistency indeed changes in time. Furthermore, significant turnover takes place in the neighborhoods of many ROIs. Network turnover is high especially for ROIs with low spatiotemporal consistency. Looking at the constituent voxels of ROIs in detail, we see that ROIs often have rich internal correlation structure that varies in time.

These results indicate that the topology of functional brain networks changes continuously on short timescales, which should be taken into account in brain network studies. Furthermore, the significant temporal variation of functional homogeneity may suggest that new, dynamical ways of defining nodes are required for creating an accurate network model of the brain. Importantly, the variation of functional homogeneity should not be seen as a technical issue that should be eliminated with some parcellation approach, but a phenomenon that carries cognitive meaning and that should be taken into account in the analysis of dynamic functional connectivity.

Using predefined ROIs as nodes of functional brain networks is based on the assumption of functional homogeneity: all voxels within an ROI are assumed to have similar dynamics that can be accurately represented by the ROI time series. To test this assumption, we calculated the distribution of spatial consistency for five commonly used parcellations of the brain: connectivity-based Brainnetome atlas and Craddock 200/400 parcellations as well as two anatomical atlases: HarvardOxford (HO) and Automated Anatomical Labeling (AAL). Spatial consistency is defined as the average Pearson correlation coefficient between the voxel time series in a ROI (see Equation 2). The results are in concordance with our earlier observations (Korhonen et al., 2017): although the maximum spatial consistency is moderately high (Brainnetome: φ_spatial = 0.53, HO: φ_spatial = 0.53, AAL: φ_spatial = 0.34, Craddock 200: φ_spatial = 0.55, Craddock 400: φ_spatial = 0.65), the distribution of spatial consistency is broad and peaks at low values (Brainnetome: φ_spatial = 0.12, HO: φ_spatial = 0.083, AAL: φ_spatial = 0.083, Craddock 200: φ_spatial = 0.12, Craddock 400: φ_spatial = 0.15) (Figure 1A). For Brainnetome or Craddock 200/400, there is no significant correlation between ROI size in voxels and spatial consistency (Brainnetome: Pearson correlation coefficient r = 0.10, p = 0.12, Figure 4A; Craddock 200: r = −8.31 × 10⁻⁴, p = 0.991; Craddock 400: r = 0.031, p = 0.538). For AAL and HO, there is a weak but significant correlation between ROI size and spatial consistency (AAL: r = −0.32, p = 4.13 × 10⁻⁴; HO: r = −0.33, p = 8.62 × 10⁻⁵). The spatial consistency investigated here was calculated over the whole measurement time series; we will from here on refer to it as static spatial consistency.

At least two different scenarios can explain the low values of static spatial consistency. On one hand, the voxels in an ROI may just have uncorrelated dynamics across the whole measurement time series. On the other hand, a moderately low level of correlation between the voxel time series may result from changes in the overall pattern, for example, there may be periods of highly correlated activity and periods of no correlations at all. In the latter scenario, one would obtain time-dependent changes in spatial consistency by dividing the measurement time series into shorter time windows. Therefore, we divided the measurement time series into five sliding windows of 80 samples each, with 50% overlap between consecutive windows, and investigated the spatial consistency separately for each time window.

We found no visible difference between distributions of spatial consistency calculated in different time windows (Figure 1B). One could assume that decreasing the time window length and the overlap between consecutive windows would cause the spatial consistency distributions to differ more between windows. However, decreasing the window length to 50 samples and the overlap to 25% did not induce more variation between time windows (data not shown). On the other hand, increasing the overlap to the largest possible value, window length - 1, did not affect the distributions of spatial consistency either (see Supporting Information Figure S29, Ryyppö, Glerean, Brattico, Saramäki, & Korhonen, 2018).

At the level of single ROIs, however, the situation is different. The spatial consistency of most ROIs changes between time windows, and the largest relative changes in spatial consistency are around 30% (Figure 2A). These changes have a nonrandom spatial distribution and seem to occur in clusters larger than simple ROIs. This is visible, for example, as the drop in spatial consistency of the frontal regions between the time windows 2 and 3. A possible reason for the similar behavior of spatially close ROIs is their assumed functional similarity; in particular, ROIs belonging to the same functional subsystem of the brain may be expected to behave similarly in terms of spatial consistency.

Similar time-dependent changes also take place in the spatial consistency ranks of ROIs (data not shown), demonstrating that the observations cannot be explained by changes in the overall level of spatial consistency. Note that in ROIs with small spatial consistency, even large relative changes may be rather small in the absolute sense; therefore the distributions of spatial consistency obtained in different time windows are almost identical at lower consistency values (see Figure 1B), although we see relative changes in consistency of almost every ROI. At higher consistency values, large relative changes are larger in the absolute sense as well. Therefore, distributions obtained in different time windows differ more from each other.

In order to quantify the amount of temporal variation in spatial consistency, we defined the spatiotemporal consistency as the inverse of the averaged relative change of spatial consistency across time windows (see Equation 3). In other words, spatiotemporal consistency is a static measure that tells how stable the ROI’s spatial consistency is over time windows on average.

Spatiotemporal consistency is not anatomically uniform (Figure 2B), and identity of the ROIs with the highest and lowest spatiotemporal consistency varies largely across subjects. This intersubject difference is partly explained by technical issues. The investigated parcellations have been defined at the group level, and they match differently with the individual anatomy and brain function of different subjects. On the other hand, differences in spatiotemporal consistency may also reflect the different cognitive responses of different subjects during free music listening.

The Brainnetome ROIs with the highest average spatiotemporal consistency include the right cuneus (5_3), left inferior frontal gyrus (6_4), right occipital gyrus (4_3), right superior occipital gyrus (2_1), and right inferior parietal lobule (6_2). In AAL, among the ROIs with the highest average spatiotemporal consistency are the left medial orbitofrontal cortex, right cerebellar area 6, left middle temporal gyrus, right insula, and left gyrus rectus. In HO, the ROIs with the highest average spatiotemporal consistency include the right supracalcarine cortex, left frontal pole, anterior division of left superior temporal gyrus, right angular gyrus, and posterior division of right middle temporal gyrus. In Craddock 200/400, ROI boundaries do not respect anatomical landmarks, and the ROIs are referred to only by numbers. For the location of the Craddock 200/400 ROIs with the highest and lowest average spatiotemporal consistency, the reader is referred to Supporting Information Figures S18 and S19 (Ryyppö et al., 2018).

In Brainnetome, the ROIs with the lowest average spatiotemporal consistency are the right parahippocampal gyrus (6_5 and 6_2), right paracentral lobule (2_1), and left postcentral gyrus (4_4 and 4_2). In AAL, the ROIs with the lowest average spatiotemporal consistency include the left cerebellar area 4_5, right supplementary motor area, left paracentral lobule, right parahippocampal gyrus, and right thalamus. In HO, the ROIs with the lowest average spatiotemporal consistency include the right hippocampus, posterior division of right parahippocampal gyrus, right precentral gyrus, right thalamus, and brain stem. In all investigated atlases, many subcortical areas are among the ROIs with low spatiotemporal consistency. We will discuss possible reasons for this later in this article.

As one possible explanation for the variation of spatiotemporal consistency across ROIs, we found a weak but significant correlation between the ROI size and the spatiotemporal consistency in the Brainnetome atlas (r = 0.24, p = 1.37 × 10⁻⁴; see Figure 4B). In AAL, HO, or Craddock atlases, there is no significant correlation between the spatiotemporal consistency and ROI size (AAL: r = 0.16, p = 0.0963; HO: r = −0.025, p = 0.770; Craddock 200: r = 0.028, p = 0.698, Craddock 400: r = 0.014, p = 0.780). The lack of correlation in the Craddock parcellations is not surprising: these parcellations aim at minimizing the variation of ROI sizes, and they have clearly smaller SDs of ROI size than the other parcellations investigated. There are a few possible reasons for why we observe a correlation in the Brainnettome but not in AAL or HO. First, Brainnetome contains more ROIs than AAL or HO. Second, unlike Brainnetome, AAL and HO contain the cerebellum where ROIs are small because of anatomical reasons, but they do not have systematically higher values of spatiotemporal consistency than the ROIs of the cerebral cortex. Finally, in AAL and HO the ROI size is negatively correlated with static spatial consistency; this is not the case for Brainnetome. On the other hand, in all atlases, ROIs with high static spatial consistency tend to have high spatiotemporal consistency as well (see below). Because of this, one would expect to see a negative correlation between ROI size and spatiotemporal consistency in HO and AAL. This negative correlation may have masked the positive correlation obtained for the Brainnetome atlas.

A sliding window with a one time frame shift is commonly used for studying dynamic functional connectivity (Keilholz, Caballero-Gaudes, Bandettini, Deco, & Calhoun, 2017). In this approach, the overlap between consecutive time windows is as large as possible: window length −1. In our case, however, this large an overlap would hide the changes in local network structure. It would also lead to extremely low values of network turnover (see below). However, we investigated how the 1 TR shift would affect the observed values of spatiotemporal consistency in the Brainnetome parcellation. As expected, using the 1 TR shift sliding window moved the distribution of spatiotemporal consistency slightly to the right (distribution peaking at φ_st = 3.1 vs φ_st = 4.9; see Supporting Information Figure S29B, Ryyppö et al., 2018): as the overlap between consecutive time series increases, there is less room for changes in spatial consistency. However, the 1 TR shift approach did not affect the overall shape of the distribution of spatiotemporal consistency, and low values of spatiotemporal consistency that indicate large relative changes in spatial consistency are observed with this approach too.

Subject motion is known to possibly affect the structure of functional brain networks (Power, Barnes, Snyder, Schlaggar, & Petersen, 2012). Therefore, one may ask if the temporal variation in spatial consistency is of genuine neurophysiological origin or if it could be explained by motion artifacts. To answer this, we investigated the temporal correlation between the mean framewise displacement (FD) and the spatial consistency concatenated across subjects. However, they did not correlate significantly for any ROI in any of the investigated atlases. The correlation between the static spatial consistency and the mean FD over subjects was not significant neither.

Temporal fluctuations in functional connectivity of the brain may underlie changes in cognitive processing (Cocchi et al., 2017). We found a significant correlation between time-resolved functional connectivity (Cocchi et al., 2017; Zalesky, Fornito, Cocchi, Gollo, & Breakspear, 2014) and spatial consistency for some ROIs of Brainnetome and HO. For further details, the reader is referred to Supporting Information Results (Ryyppö et al., 2018).

In the present article, we investigate five atlases: Brainnetome, Craddock 200/400, AAL, and HO. Despite the differences between these atlases, we obtained highly similar results for all of them. In the main article, we concentrate on the results obtained with the Brainnetome atlas; for detailed results and visualizations for the Craddock 200/400, AAL, and HO atlases, the reader is referred to Supporting Information Results (Ryyppö et al., 2018).

To verify that the results generalize, we repeated all analyses for a second, independent dataset from the ABIDE I initiative (Di Martino et al., 2014). The results obtained using the ABIDE data were very similar to those reported here; full details can be found in the Supporting Information Results (Ryyppö et al., 2018).

The structure of functional brain networks is known to change in time. For an individual node, this means that the local structure around the node, that is the identity of its neighbors, may change. This change can be quantified in terms of the Jaccard index between the node’s sets of neighbors in consecutive time windows. We defined an ROI’s closest neighborhood as its 35 most strongly linked neighbors and investigated the Jaccard index. Indeed, we found significant changes in ROIs’ neighborhoods in time (Figure 3A). ROIs with the highest neighborhood turnover may change up to 75% of their closest neighbors between two time windows, corresponding to a Jaccard index of 0.25. Even the ROIs with the most stable neighborhoods only reach a Jaccard index value of 0.55, meaning that half of their closest neighborhood changes between consecutive time windows. For comparison, shuffling the weights of 5% of randomly chosen links in the full network for 1000 times yields an average Jaccard index of 0.89 ± 0.062 (mean ± SD).

The spatial distribution of the Jaccard index over the ROIs appears stable in time. In particular, subcortical ROIs have clearly lower Jaccard index values than cortical ROIs independent of the time window pair investigated. We will discuss possible reasons for this later.

Network turnover, defined as the complement of the mean Jaccard index of the ROI’s neighborhood over time (see Equation 4), quantifies the overall tendency of the neighborhood to change in time. Like spatiotemporal consistency, network turnover varies between ROIs (Figure 3B). As expected because of their low Jaccard index values, subcortical ROIs have higher network turnover than cortical ROIs. Network turnover is spatially strongly correlated: anatomically adjacent ROIs tend to have similar values of network turnover.

The Brainnetome ROIs that have the highest average network turnover include the left and right parahippocampal gyrus (6_5), left thalamus (8_2), and right and left inferior temporal gyrus (7_1). In AAL, these include vermis 9, the left caudate nucleus, left cerebellar area 3, vermis 1_2, and right olfactory cortex. In HO, the ROIs with the highest average network turnover include the right and left pallidum, anterior division of left temporal fusiform cortex, vermis X, and vermis VIIIb. For the Craddock 200/400 ROIs with the highest and lowest network turnover, the reader is referred to Supporting Information Figures S18 and S19 (Ryyppö et al., 2018). There is some variation in the identity of the highest network turnover ROIs across subjects; however, subcortical areas tend to have high network turnover in all subjects.

Identity of ROIs with the lowest network turnover vary a lot across subjects. The ROIs with the lowest average network turnover include in Brainnetome the left occipital gyrus (4_1), left middle temporal gyrus (4_1), right superior occipital gyrus (2_2), and left superior frontal gyrus (7_7 and 7_3). In AAL, they include the right fusiform cortex, right cerebellar area 6, right superior occipital gyrus, left angular gyrus, and right middle occipital gyrus. In HO, the low average network turnover ROIs include the left frontal pole, left middle frontal gyrus, left angular gyrus, left paracingulate gyrus, and left cuneal cortex.

In addition to spatial variation, we found significant negative correlation between an ROI’s size and network turnover (Brainnetome: r = −0.60, p ≪ 10⁻⁵; HO: r = −0.30, p = 4.02 × 10⁻⁴; AAL: r = −0.42, p ≪ 10⁻⁵; Craddock 200: r = −0.21, p = 0.00232; Craddock 400: r = −0.41, p ≪ 10⁻⁵; Figure 4C). This correlation is most probably dominated by the very high network turnover values of the subcortical ROIs that, for anatomical reasons, tend to be smaller than cortical ROIs. In the AAL and HO atlases, the correlation may have been partly shadowed by the lower number of ROIs and the presence of cerebellar ROIs that are small but do not have systematically lower network turnover values than ROIs of the cerebral cortex.

Next, we asked how spatiotemporal consistency and network turnover relate to each other. At the group level, that is, averaging the spatiotemporal consistencies and turnovers over subjects, we found a significant negative correlation between these measures (Brainnetome: r = −0.42, p ≪ 10⁻⁵; HO: r = −0.44, p ≪ 10⁻⁵; AAL: r = −0.38, p = 2.08 × 10⁻⁵; Craddock 200: r = −0.46, p ≪ 10⁻⁵; Craddock 400: r = −0.42, p ≪ 10⁻⁵) (Figure 5A). In other words, ROIs with the highest spatiotemporal consistency have the lowest amount of turnover in their neighborhoods. These ROIs also have the highest static spatial consistency (Figure 5B).

The correlation between spatiotemporal consistency and network turnover was also visible, albeit weaker, at the level of single subjects (Brainnetome: significant [p < 0.05] negative correlation for 6 subjects, negative but nonsignificant correlation for 5 subjects, nonsignificant positive correlation for 1 subject; HO: significant negative correlation for 7 subjects, negative but nonsignificant correlation for 5 subjects, nonsignificant positive correlation for 1 subject; AAL: significant negative correlation for 4 subjects, negative but nonsignificant correlation for 5 subjects, nonsignificant positive correlation for 4 subjects). A plausible reason for the weaker and less significant correlations obtained for AAL is the number of data points: AAL contains fewer ROIs (116) than Brainnetome (246) or HO (138), which may have made it more difficult to obtain a significant correlation.

The connectivity profiles of ROIs with low network turnover change only little between time windows and they resemble the connectivity profiles obtained over the whole time series. Therefore, ROIs with low network turnover should have stronger links in the networks extracted from the whole time series. As low-turnover ROIs tend to have high static spatial consistency, it is not too surprising that we found a correlation between static spatial consistency and degree in networks extracted from the whole time series. At 2.5% density, this correlation was significant for all investigated parcellations (Brainnetome: r = 0.48, p ≪ 10⁻⁵; AAL: r = 0.31, p ≪ 10⁻⁵; HO: r = 0.41, p ≪ 10⁻⁵; Craddock 200: r = 0.50, p ≪ 10⁻⁵; Craddock 400: r = 0.57, p ≪ 10⁻⁵). The correlation remained significant also for higher network densities; the highest density where the correlation was present varied between parcellations (Brainnetome: d = 45.0, r = 0.037, p = 0.0368; AAL: d = 10.0, r = 0.11, p ≪ 10⁻⁵; HO: d = 10.0, r = 0.057, p = 0.0155; Craddock 200: d = 30.0, r = 0.058, p = 0.00309; Craddock 400: d = 40.0, r = 0.0603, p = 1.67 × 10⁻⁵).

In Brainnetome, the ROIs with the highest spatiotemporal consistency and lowest network turnover tend to be larger than ROIs with lower spatiotemporal consistency and higher network turnover (Figure 5C). This is as one may expect based on the correlations between spatiotemporal consistency and ROI size, and network turnover and ROI size (Figure 4B and C). In AAL, HO, and Craddock 200/400, this relationship is less clear (see Supporting Information Figures S12C, S13C, S14C, and S15C, Ryyppö et al., 2018).

The relationship between spatiotemporal consistency and network turnover strongly depends on how we define spatiotemporal consistency. The definition given in Equation 3 measures relative changes in spatial consistency. To get a more complete picture, we also investigated an alternative definition of spatiotemporal consistency that measures absolute changes. For details, the reader is referred to the reader is referred to Supporting Information Results (Ryyppö et al., 2018).

So far, we have investigated the relationship between spatiotemporal consistency and network turnover at the population level. Next, we asked which specific ROIs are the ones with the highest and lowest values of spatiotemporal consistency and network turnover. To this end, we obtained two groups of extreme ROIs by applying principal component analysis (PCA) in the space spanned by spatiotemporal consistency and network turnover. The extreme groups contain the five ROIs with the largest and smallest projected coordinates on the first principal component. ROIs of the first group have high spatiotemporal consistency and low network turnover, and ROIs of the second group have lower spatiotemporal consistency and high network turnover (Figure 6). As the PCA has only two degrees of freedom, the extreme groups could in principle have been defined by visual inspection alone; the main reason for applying PCA was to avoid subjectivity and to ensure that the extreme groups are defined similarly in all investigated parcellations.

In Brainnetome, the first group contains the right cuneus (5_3), right superior occipital gyrus (2_1), left inferior frontal gyrus (6_4), right inferior parietal lobule (6_2) and left occipital gyrus (4_1). The AAL ROIs of this group are the right cerebellar area 6, left medial orbitofrontal cortex, right superior occipital gurys, left angular gyrus and right middle occipital gyrus. In HO, this group comprises the left frontal pole, right and left supracalcarine cortex, left middle frontal gyrus and right angular gyrus. For Craddock 200/400 ROIs belonging to the extreme groups, the reader is referred to Supporting Information Figures S18 and S19 (Ryyppö et al., 2018).

The second group, that is, the ROIs with low spatiotemporal consistency and high network turnover, contains in Brainnetome the left and right parahippocampal gyrus (6_5), right inferior temporal gyrus (7_1), left thalamus (8_8) and left striatum (6_5). The AAL ROIs that belong to this group are the right globus pallidum, left paracentral lobule, right olfactory cortex, right cerebellar area 9 and Vermis 1_2. In HO, this group contains the right and left pallidum, brain stem, right hippocampus, and right thalamus. In all five parcellations, most ROIs of this group are relatively small areas located deep in the brain. Because of the location, the signal-to-noise ratio (SNR) of the fMRI measurement tends to be low for these areas. This may at least partially explain their low spatiotemporal consistency and may also limit the accuracy of estimating their network connectivity, leading to noisy closest neighborhoods and high turnover.

The sets of extreme ROIs in different parcellations are not the same, but this is to be expected. First, the ROIs of different parcellations have different shapes, sizes, and locations. Second, there are many ROIs with spatiotemporal consistency and network turnover that are rather close to those of the five most extreme ROIs; this hard threshold is of course arbitrary.

From both groups of extreme ROIs, we selected two ROIs for a more detailed investigation. We chose the most extreme ROIs that were not exceptionally small or too large for the visualization discussed below. In Brainnetome, the selected ROIs were the left inferior frontal gyrus (6_4) and right superior occipital gyrus (2_1) from the high-spatiotemporal-consistency-low-network-turnover group and the left striatum (6_5) and left thalamus (8_8) from the opposite group. We calculated voxel-level correlation matrices to reveal the internal correlation structure inside these ROIs (Figure 7).

The two groups are visibly different in terms of their correlation matrices: the overall correlation level is clearly higher for ROIs with high spatiotemporal consistency than for ROIs with low spatiotemporal consistency. The voxel-level correlations are not, however, uniformly distributed. Instead, a division into several internally highly correlated subareas is visible inside ROIs with high spatiotemporal consistency and ROIs with low spatiotemporal consistency.

This internal structure of ROIs is seen to change in time. In the right superior occipital gyrus that has high spatial and spatiotemporal consistency, the voxels are uniformly correlated across the whole ROI in time windows 1 and 2 but separate into two clusters between time windows 2 and 3. Similarly, the left thalamus that has low spatial and spatiotemporal consistency shows time-dependent internal cluster structure.

The internal structure of ROIs may affect their spatiotemporal consistency and static spatial consistency in several ways. For example, stable internal structure should manifest itself as high spatiotemporal consistency, because the average voxel-level correlation does not change in time. Fewer and larger subareas lead to larger amounts of correlated voxels within the ROI and should therefore be associated with higher static spatial consistency. On the other hand, a large number of small subareas should lead to low static spatial consistency, as should a total lack of internal structure.

The use of ROIs as nodes of fMRI brain networks assumes functional homogeneity: each of the ROI’s voxels is thought to have similar dynamics, and therefore the ROI time series is considered as an accurate representation of the voxel-level dynamics. Earlier (Korhonen et al., 2017), we have shown that this assumption does not hold for the ROIs of commonly used parcellations. To this end, we used spatial consistency, a measure of functional homogeneity defined as the mean Pearson correlation coefficient between voxel time series inside an ROI.

Functional homogeneity is often considered as a static ROI property. However, functional brain networks change in time, even on short timescales (Bassett et al., 2011; Göttlich et al., 2017; Honey et al., 2007). Here, we investigated the temporal behavior of spatial consistency. We divided fMRI data measured during a free music listening task into time windows and calculated the relative change in spatial consistency between them. For quantifying the temporal variation in spatial consistency, we introduced spatiotemporal consistency as the inverse of the mean relative change in spatial consistency over time windows (see Equation 3). We found that spatial consistency changes significantly in time, the largest relative changes being up to 30%, resulting in low spatiotemporal consistency.

The concept of dynamic functional connectivity has been recently debated among the neuroscientific community. Although many studies have reported time-dependent changes in the structure of functional brain networks (Bassett et al., 2011; Cocchi et al., 2017; Honey et al., 2007; Zalesky et al., 2014), the neurophysiological meaning of these changes is not fully understood (Keilholz et al., 2017; Preti, Bolton, & Van De Ville, 2016). An fMRI measurement is always only a single realization of the underlying stochastic process, and it may therefore show connectivity fluctuations even if the underlying process is stationary (Liegeois, Laumann, Snyder, Zhou, & Yeo, 2017). However, obtaining multiple realizations of the exact same process is impossible; the measurements of different subjects as well as the measurements of same subject at different times are different processes (Liegeois et al., 2017). Therefore, it is hard to construct a proper null model for evaluating the statistical significance of dynamic functional connectivity (Liegeois et al., 2017; Miller, Adali, Levin-Schwartz, & Calhoun, 2017). We do not use a stationary null model in the present study, similarly to many other studies.

If one wants to investigate in detail whether the observed changes in spatial consistency are meaningful, two different paths can be taken. First, testing the results against a statistically rigorous null model would ensure their significance. Second, the neurophysiological meaning of spatiotemporal consistency could be addressed by comparing consistencies obtained for data measured during rest and during different tasks: obtaining similar changes in response to a common stimulus in a group of subjects can be considered as indicative of some real underlying mechanism, even if no formal null model is applied.

The structure of functional brain networks changes with cognitive tasks (Braun et al., 2015; Chan et al., 2017; Göttlich et al., 2017), and increased local connectivity can be associated with increased activity and cognitive demand (Hearne, Cocchi, Zalesky, & Mattingley, 2017; Jiang & Zuo, 2016; Zang, Jiang, Lu, He, & Tian, 2004). If the temporal variation in spatial consistency is related to changes in brain function, one would expect to see simultaneous changes in network structure as well. Indeed, there was clear turnover in the closest neighborhoods of ROIs, “closest” being defined as the 35 most strongly connected neighbors. This turnover was lower for ROIs with high spatial and spatiotemporal consistency; however, even these ROIs changed up to half of their closest neighbors between consecutive time windows. This indicates that the local structure of functional brain networks truly changes on short time scales. Furthermore, the network turnover investigated here only quantifies the changes in the identities of the closest neighbors, but does not take into account changes in connection strengths within the closest neighborhood. Therefore, prominent changes may take place in the ranks of the closest neighbors of even an ROI with moderately low network turnover.

We saw that network turnover varies across ROIs. One may speculate about how this variation may relate to the ROIs’ different functional roles. It is possible that some ROIs need a diverse and varying set of connections for performing their cognitive tasks, whereas others require a stable set of neighbors. However, there may be a more straightforward explanation for the variation in network turnover. The ROIs with the highest network turnover are subcortical and cerebellar areas that also have low spatiotemporal and spatial consistency. The SNR of the signals originating from these ROIs is known to be low in fMRI measurements (Glasser et al., 2016). This may partially explain their low spatiotemporal and spatial consistency, and also suggests that their connectivity may be inaccurately mapped. Therefore, their extreme network turnover may be partially explained by measurement noise.

We obtained temporal changes in both the spatial consistency and the closest neighborhoods of ROIs. However, we did not investigate the exact timing of these changes. If neighborhood turnover and variation in spatial consistency are both caused by changes in brain activity, these changes should be more or less simultaneous. This would result in a temporal correlation between the variation of spatial consistency and the neighborhood turnover. The datasets used in the present study—free music listening and resting state—may not necessarily be optimal for this kind of investigation. Cognitive responses to the music may differ between subjects, and in the resting-state subjects are instructed to let their mind wander uncontrolled. A more detailed investigation of the connection between spatial consistency and turnover would require a dataset with more control on the timing of putative activity changes. This could be achieved with the traditional block design, where stimuli are repeated at set intervals and the timeline is divided into blocks (see, e.g., Tie et al., 2009). However, the shortness of the blocks is problematic from the viewpoint of network studies: reliable estimation of a functional brain network requires time series significantly longer than typical block lengths.

Functional networks are constructed using only the averaged ROI time series, and the only feature that is used in any subsequent analysis is the ROI’s location on the brain surface. At the same time, their size, shape, and in particular internal connectivity are typically ignored. This view of ROIs as featureless entities may, however, be largely oversimplified. We found rich, time-dependent structure of voxel-level correlations inside ROIs. Considering the complexity of the brain and the small number of ROIs and their connections to which this complexity is reduced, this is not surprising at all.

The ROIs that we investigated have very different looking internal structures. These are not necessarily reflected in their consistency measures; in Figure 7 the two uppermost ROIs have high spatiotemporal consistency, but their correlation matrices display different kinds of structures. The same applies to the two low-consistency bottom rows.

Why do ROIs have different kinds of internal structures? A plausible hypothesis is that the correlation structure inside an ROI relates to an ROI’s functional role. Let us consider local and connector hubs (Bullmore & Sporns, 2009; Guimerà & Amaral, 2005) as an example. Local hub nodes are central in their local network modules and have only few connections to nodes outside of their own module, whereas connector hubs act as bridges between different modules. So, could one separate local and connector hubs from each other in terms of their internal structure? Local hubs are connected only to a relatively stable neighborhood; one might expect that the voxel-level correlation distribution inside them is relatively uniform, and periods of high and low voxel-level correlations reflect changes in the activity of the ROI. Connector hubs, on the other hand, need to be able to connect to several different network modules; an internal structure of diverse subareas could help in this.

When ROIs are used as nodes of functional brain networks, the brain is assumed to contain a set of static functional areas. An optimal parcellation of the brain then maps to these areas, resulting in functionally homogeneous ROIs. If the static-area assumption holds, low functional homogeneity of ROIs then only tells about inaccuracies in ROI definitions, which can be corrected by a more accurate parcellation scheme.

Time-dependent changes have been reported in the module structure of functional brain networks at the ROI level (Khambhati, Sizemore, Betzel, & Bassett, in press). Our results suggest that similar changes also occur in the voxel-level correlation structure inside ROIs. The dynamically changing internal connectivity of ROIs appears to challenge the assumption of static functional areas. Because the ROIs of multiple different parcellations have a time-varying structure, it seems plausible that the changing internal connectivity is not a technical issue that can be fixed by an elaborate parcellation scheme. Instead, it may be a genuine feature and related to how the brain works. If so, it may even be impossible to define ROIs in a way that would make them persistently homogeneous.

Many problems caused by the functional inhomogeneity of ROIs can be overcome by using fMRI measurement voxels as nodes of brain networks (Fornito, Zalesky, & Breakspear, 2013; Hayasaka & Laurienti, 2010). However, there is evidence for existence of functional areas larger than single voxels (Shen, Tokoglu, Papademetris, & Constable, 2013; Wig et al., 2011), which motivates investigating brain networks above the level of voxels as well. For example, Preti & Van De Ville (2017) have suggested an approach for parcellating the brain based on the dynamic connectivity of voxels; these parcels would be an interesting option for defining network nodes. Similarly, the negative correlation we observed between network turnover and spatial consistency could be used for defining ROIs: defining ROIs as clusters of voxels that have minimal network turnover should produce ROIs with reasonable spatial consistency.

An optimal network model of the brain should measure the dynamic connectivity between clusters of voxels and also quantify the changing internal structure of these clusters. In the coarse-graining approach by Kujala et al. (2016), self-links are used to model changes in internal connectivity. As long as static sets of nodes are used to model the time-dependent connectivity of the human brain, outcomes of brain network analysis may be surprisingly inaccurate. Therefore, network neuroscience would greatly benefit from node definition strategies tailored for dynamic networks.

fMRI data of 13 subjects (7 women, 6 men, age 28.70 ± 10.17 years, mean ± SD, 1 left-handed, 12 right-handed) were used in the present study. The data were collected as a part of a study of functional connectivity during music listening, containing both musicians and non-musicians freely listening to music, and have been earlier described in Alluri et al. (2015, 2017) and Burunat et al. (2015). For the present study, we used the parts of the dataset that were readily available at the Department of Neuroscience and Biomedical Engineering, Aalto University. The subjects used in the present study were considered as non-musicians, that is, they had no formal musical training.

All participants signed an informed consent on arrival to the laboratory and received compensation for the use of their time. All experimental procedures for this study, included in the broad research protocol termed Tunteet, were approved by the Coordinating Ethics Committee of the Hospital District of Helsinki and Uusimaa (the approval number 315/13/03/00/11, obtained on March 11, 2012). All procedures were conducted in agreement with the ethical principles of the Declaration of Helsinki.

Functional magnetic resonance imaging (fMRI) data were acquired using a 3T MAGNETOM Skyra scanner (Siemens Healthcare, Erlangen, Germany) and a standard 32-channel head-neck coil in the AMI Centre (Aalto Neuroimaging, Aalto University, Espoo, Finland). A T2*-weighted whole-brain EPI sequence was measured with the following parameters: TR = 2 s, 33 oblique slices, TE = 32 ms, flip angle = 75°, voxel size = 3×3×4 mm³, FOV = 192 × 192 mm², matrix size = 64 × 64. T1-weighted structural magnetic resonance images (MRI) were acquired with the following parameters: 176 slices, FOV = 256 × 256 mm², matrix size = 256 × 256, slice thickness = 1 mm.

During the measurement, subjects were instructed to fix their gaze on the screen and actively listen to a musical stimulus (Adios, Nonino by Astor Piazzolla) via MR-compatible insert earphones. Foam was used to suppress the noise caused by the imaging gradients. Duration of the stimulus, and therefore of the measured time series, was 8.13 min (244 samples).

The data were preprocessed with FSL software (www.fmrib.ox.ac.uk, version 5.0.9) and custom in-house MATLAB code (BRAMILA pipeline v2.0, available at https://version.aalto.fi/gitlab/BML/bramila) following the standard fMRI preprocessing steps. This included EPI slice time correction as well as head motion correction using MCFLIRT. The data were coregistered to the Montreal Neurological Institute (MNI) 152 2 mm standard template using FLIRT two-step procedure where the EPI volumes were first registered to the anatomical image of participants brain (9 degrees of freedom) and the participants anatomical image was then registered to the standard template (12 degrees of freedom). No spatial smoothing was applied, but a 240-sec-long cubic Savitzky-Golay filter (Çukur, Nishimoto, Huth, & Gallant, 2013) was used to remove scanner drift, and the BOLD time series were filtered using a Butterworth bandpass filter at 0.01–0.08 Hz. For increased control of motion and physiological artifacts, 24 motion-related regressors, signal from deep white matter, ventricles and cerebrospinal fluid were regressed out of the BOLD time series (Power et al., 2014).

Voxels with over 70% of their variance explained by motion or signal from tissues other than the grey matter were removed from the analysis.

Some of the parcellations used in this study, in particular AAL and Craddock 200/400, are known to show rather low mean functional homogeneity across ROIs (Gordon et al., 2014). However, the ROIs of these parcellations are commonly used as nodes of functional brain networks, and therefore we have chosen to include them in our study.

The Brainnetome atlas (Fan et al., 2016) is based on combination of structural and functional connectivity measured by multimodal imaging techniques. In the present study, we used 246 Brainnetome ROIs. Of these, 210 ROIs were located in the cerebral cortex, while 36 ROIs covered subcortical gray matter. Note that the Brainnetome atlas does not include cerebellar ROIs.

Size of the Brainnetome ROIs varied between 6 and 1,102 with a median of 414. Mean ROI size was 424.02 ± 222.76 (mean ± SD).

The connectivity-based Craddock parcellations (Craddock, James, Holtzheimer, Hu, & Mayberg, 2012a) have been obtained by applying a two-level normalized cut spectral clustering algorithm on the voxel-level resting-state connectivity matrix. In the present study, we investigate the Craddock 200 and Craddock 400 parcellations that contain 200 and 392 ROIs, respectively, covering the cerebral cortex, subcortical areas, and cerebellum.

The sizes of the Craddock 200 ROIs varied between 202 and 1,239 with a median of 706.5. The mean ROI size was 689.96 ± 168.12 (mean ± SD). The sizes of the Craddock 400 ROIs varied between 56 and 600 with a median of 354.5. The mean ROI size was 352.02 ± 88.787 (mean ± SD).

AAL (Tzourio-Mazoyer et al., 2002) is an anatomical atlas that has been obtained by parcellating a spatially normalized high-resolution single-subject structural volume based on the main sulci. After the parcellation, each ROI has been automatically associated with a label. We used 116 AAL ROIs, 90 of which were located in the cerebral cortex, 8 in the subcortical gray matter, and 18 in the cerebellum.

Size of the AAL ROIs varied between 44 and 4,370 with a median of 1,158.5 and a mean of 1,366.01 ± 929.64.

The HO atlas (http://neuro.debian.net/pkgs/fsl-harvard-oxford-atlases.html; Desikan et al., 2006) is a probabilistic atlas, where the brain is divided into ROIs based on macroanatomical boundaries. We used HO ROIs at the probability level of 30%, meaning that each voxel belonged to the ROI it is associated with in 30% or more of the subjects used to construct the atlas. We used 138 HO ROIs, out of which 96 were located in the cerebral cortex, 15 covered subcortical gray matter, and 27 were located in the cerebellum. Note that one of the cerebellar ROIs of the HO atlas (Vermis Crus I) is not defined at the probability level of 30%. Therefore, this ROI is not included in the present study.

Size of the HO ROIs varied between 28 and 5,578 with a median of 633.5 and a mean of 915.63 ± 921.83 (mean ± SD).

To construct the dynamic functional brain networks, the time series were divided into time windows of 80 samples. This corresponds to a duration of 160 s. The consecutive time windows had an overlap of 50%. This resulted in a total of five time windows along the duration of the scan.

The window length was selected so that we were able to investigate the changes of spatial consistency and local network structure (see below) across as many windows as possible, but the values of spatial consistency were not affected by the short window length. The window length that we used was selected so that further increasing it did not increase the value of spatial consistency obtained in the window (see Supporting Information Methods, Ryyppö et al., 2018, for details on selecting the window length). It has been suggested that the time window length should be equal or larger than 1/f_min where f_min is the lowest signal frequency present in the data (Leonardi & Van De Ville, 2015; Shakil, Billings, Keilholz, & Lee, 2017; Shakil, Keilholz, & Lee, 2015); the selected window length fulfills this condition. Furthermore, time windows of similar length have been used for constructing dynamic brain networks in the literature (Bassett et al., 2011, 2013).

We computed the ROI-level adjacency matrix A separately in each of the time windows. The elements A_IJ of the adjacency matrix quantified the connectivity between ROIs I and J, defined as Pearson correlation coefficient between their ROI time series. The diagonal of the adjacency matrix was set to zero in order to remove self-links. No thresholding of the correlation values was performed at this stage.

The behavior of turnover as a function of the size of the neighborhood varies between ROIs, especially in small neighborhoods (for details, see Supporting Information Methods, Ryyppö et al., 2018). We selected the neighborhood size so that this variation associated with small neighborhoods has stabilized but the trivial decrease of turnover due to large neighborhood size had not yet started.

To ensure that our results are not explained by any feature of our in-house dataset, we repeated all analyses of the present study for a secondary, independent dataset to which we from now on will refer as the ABIDE dataset. This dataset was part of the Autism Brain Imaging Data Exchange I (ABIDE I) initiative (Di Martino et al., 2014) and contained resting-state data of 28 healthy adult subjects. Importantly, data of these subjects were collected with the same TR as our in-house data (TR = 2.0s); differences in TR could have caused unexpected effects in correlation-based measures. Details about the ABIDE data can be found in Supporting Information Methods (Ryyppö et al., 2018).

Ryyppö, E., Glerean, E., Brattico, E., Saramäki, J., & Korhonen, O. (2018). Supporting Infor mation for “Regions of Interest as nodes of dynamic functional brain networks.” Network Neuroscience, 2(4), 513–535. https://doi.org/10.1162/netn_a_00047

Elisa Ryyppö: Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Software; Validation; Visualization; Writing - original draft; Writing - review and editing. Enrico Glerean: Conceptualization; Methodology; Software; Validation; Writing - review and editing. Elvira Brattico: Data curation; Resources; Writing - review and editing. Jari Saramäki: Conceptualization; Supervision; Writing - review and editing. Onerva Korhonen: Conceptualization; Formal analysis; Methodology; Project administration; Software; Supervision; Validation; Visualization; Writing - original draft; Writing - review and editing.

We wish to thank the Danish National Research Foundation (project DNRF 117).

We acknowledge the computational resources provided by the Aalto Science-IT project. Moreover, we are grateful to Brigitte Bogert, Benjamin Gold, Marina Kliuchko, David Ellison, Taru Numminen-Kontti, Mikko Heimälä, Jyrki Mäkelä, Marita Kattelus, and Toni Auranen for their contribution to data collection. We wish also to thank Petri Toiviainen, Vinoo Alluri, and Iballa Burunat for their help in several aspects of this project.

Functional homogeneity:

Every voxel of an ROI is performing some particular function and therefore has (roughly) similar dynamics, yielding strongly correlated voxel time series.

Turnover:

Differences between a node’s sets of neighbors in two consecutive time windows.

Static spatial consistency:

A measure of functional homogeneity; defined as the mean intra-ROI Pearson correlation coefficient calculated over whole voxel time series.

Sliding windows:

A set of time windows defined so that subsequent windows overlap with each other.

Jaccard index:

A measure of similarity between two sets; defined as the proportion of intersection of the sets to their union.

Closest neighborhood:

The set of nodes most strongly connected to the focal node; here, the size of the closest neighborhood is 35.