Graph coarse-graining reveals differences in the module-level structure of functional brain networks

TL;DR

This paper introduces a network coarse-graining method to analyze differences in module-level structure of functional brain networks derived from fMRI data, revealing condition-specific connectivity changes.

Contribution

The study develops a novel coarse-graining approach for comparing brain network modules across conditions, improving interpretability and statistical assessment.

Findings

Increased intra-module links in visual cortex during movie viewing

Method reveals differences not seen with independent partitioning

Provides statistically significant module-level connectivity insights

Abstract

Network analysis is rapidly becoming a standard tool for studying functional magnetic resonance imaging (fMRI) data. In this framework, different brain areas are mapped to the nodes of a network, whose links depict functional dependencies between the areas. The sizes of the areas that the nodes portray vary between studies. Recently, it has been recommended that the original volume elements, voxels, of the imaging experiment should be used as the network nodes to avoid artefacts and biases. However, this results in a large numbers of nodes and links, and the sheer amount of detail may obscure important network features that are manifested on larger scales. One fruitful approach to detecting such features is to partition networks into modules, i.e. groups of nodes that are densely connected internally but have few connections between them. However, attempting to understand how functional networks differ by simply comparing their individual modular structures can be a daunting task, and results may be hard to interpret. We show that instead of comparing different partitions, it is beneficial to analyze differences in the connectivity between and within the very same modules in networks obtained under different conditions. We develop a network coarse-graining methodology that provides easily interpretable results and allows assessing the statistical significance of observed differences. The feasibility of the method is demonstrated by analyzing fMRI data recorded from 13 healthy subjects during rest and movie viewing. While independent partitioning of the networks corresponding to the the two conditions yields few insights on their differences, network coarse-graining allows us to pinpoint e.g. the increased number of intra-module links within the visual cortex during movie viewing.

Tables

Table 1. Table 1: Differences in link numbers within and between rest consensus modules that survive the 0.05 Benjamini-Hochberg FDR correction.

Module pair $i,j$	$p$-value	$\Delta W_{i,j}$	relative increase
VIS,VIS	0.00122	15184.7	46%
VIS,PCUN	0.00024	2624.6	60%
VIS,VISx	0.00171	7865.6	47%
VIS,CRBL/VT	0.00146	3312.1	170%
AUD,DM	0.00293	3646.3	100%

Table 2. Table S1: Spatial correlation of the rest consensus modules with reported modules in the literature.

Module ID	Yeo et al 2011	corr.	Power et al 2011	corr.	Name given	Abbreviation
rest 1	Limbic	0.066	Subcortical	0.84	Limbic	LIM
rest 2	Limbic	0.21	Default mode	0.1	Cerebellum / ventro temporal	CRBL/VT
rest 3	Default mode	0.54	Default mode	0.66	Default mode	DM
rest 4	Visual	0.064	Memory retrieval	0.29	Precuneus	PCUN
rest 5	Visual	0.5	Visual	0.6	Visual	VIS
rest 6	Somatomotor	0.28	Auditory	0.53	Auditory	AUD
rest 7	Ventral Attention	0.27	Cingulo-opercular task control	0.46	Salience	SAL
rest 8	Frontoparietal	0.41	Fronto-parietal task control	0.68	Fronto-parietal	FP
rest 9	Visual	0.31	Visual	0.53	Visual-extrastriate	VISx
rest 10	Somatomotor	0.39	Sensory/ somatomotor hand	0.78	Sensorimotor	SM
rest 11	Default mode	0.079	Fronto-parietal task control	0.13	Language	LAN

Table 3. Table S2: Spatial correlation of the movie consensus modules with reported modules in the literature.

Module ID	Yeo et al 2011	corr.	Power et al 2011	corr.	Name given	Abbreviation
movie 1	Limbic	0.24	Subcortical	0.83	Ventro-temporal limbic	VTL
movie 2	-	-	Default mode	0.02	Cerebellum	CRBL
movie 3	Default	0.39	Default mode	0.46	Default mode	DM
movie 4	Visual	0.14	Visual	0.21	Cuneus	CUN
movie 5	Visual	0.62	Visual	0.80	Visual	VIS
movie 6	Somatomotor	0.19	Auditory	0.38	Auditory	AUD
movie 7	Frontoparietal	0.16	Salience	0.43	Salience	SAL
movie 8	Frontoparietal	0.39	Fronto-parietal task control	0.59	Fronto-parietal	FP
movie 9	Dorsal attention	0.4	Dorsal attention	0.42	Dorsal attention	DA
movie 10	Somatomotor	0.42	Sensory/ somatomotor hand	0.71	Sensorimotor	SM

Table 4. Table S3: Abbreviations of node labels

Anatomical name	AAL (cerebral cortex) & HO (sub-cortex and cerebellum) labels	Major region label	Acronym (as per doi:10.1371/ journal.pone.0005226)
Precentral gyrus (L)	Precentral_L	Frontal	PreCG(L)
Precentral gyrus (R)	Precentral_R	Frontal	PreCG(R)
Superior frontal gyrus (L)	Frontal_Sup_L	Frontal	SFGdor(L)
Superior frontal gyrus (R)	Frontal_Sup_R	Frontal	SFGdor(R)
Orbital superior frontal gyrus (L)	Frontal_Sup_Orb_L	Frontal	ORBsup(L)
Orbital superior frontal gyrus (R)	Frontal_Sup_Orb_R	Frontal	ORBsup(R)
Middle frontal gyrus (L)	Frontal_Mid_L	Frontal	MFG(L)
Middle frontal gyrus (R)	Frontal_Mid_R	Frontal	MFG(R)
Orbital middle frontal gyrus (L)	Frontal_Mid_Orb_L	Frontal	ORBmid(L)
Orbital middle frontal gyrus (R)	Frontal_Mid_Orb_R	Frontal	ORBmid(R)
Opercular inferior frontal gyrus (L)	Frontal_Inf_Oper_L	Frontal	IFGoperc(L)
Opercular inferior frontal gyrus (R)	Frontal_Inf_Oper_R	Frontal	IFGoperc(R)
Triangular inferior frontal gyrus (L)	Frontal_Inf_Tri_L	Frontal	IFGtriang(L)
Triangular inferior frontal gyrus (R)	Frontal_Inf_Tri_R	Frontal	IFGtriang(R)
Orbital inferior frontal gyrus (L)	Frontal_Inf_Orb_L	Frontal	ORBinf(L)
Orbital inferior frontal gyrus (R)	Frontal_Inf_Orb_R	Frontal	ORBinf(R)
Rolandic operculum (L)	Rolandic_Oper_L	Frontal	ROL(L)
Rolandic operculum (R)	Rolandic_Oper_R	Frontal	ROL(R)
Supplementary motor area (L)	Supp_Motor_Area_L	Frontal	SMA(L)
Supplementary motor area (R)	Supp_Motor_Area_R	Frontal	SMA(R)
Olfactory cortex (L)	Olfactory_L	Frontal	OLF(L)
Olfactory cortex (R)	Olfactory_R	Frontal	OLF(R)
Medial superior frontal gyrus (L)	Frontal_Sup_Medial_L	Frontal	SFGmed(L)
Medial superior frontal gyrus (R)	Frontal_Sup_Medial_R	Frontal	SFGmed(R)
Orbital medial frontal gyrus (L)	Frontal_Mid_Orb_L	Frontal	ORBsupmed(L)
Orbital medial frontal gyrus (R)	Frontal_Mid_Orb_R	Frontal	ORBsupmed(R)
Gyrus rectus (L)	Rectus_L	Frontal	REC(L)
Gyrus rectus (R)	Rectus_R	Frontal	REC(R)
Insula (L)	Insula_L	Insula	INS(L)
Insula (R)	Insula_R	Insula	INS(R)
Anterior cingulum (L)	Cingulum_Ant_L	Cingulate	ACG(L)
Anterior cingulum (R)	Cingulum_Ant_R	Cingulate	ACG(R)
Middle cingulum (L)	Cingulum_Mid_L	Cingulate	DCG(L)
Middle cingulum (R)	Cingulum_Mid_R	Cingulate	DCG(R)
Posterior cingulum (L)	Cingulum_Post_L	Cingulate	PCG(L)
Posterior cingulum (R)	Cingulum_Post_R	Cingulate	PCG(R)
Parahippocampal gyrus (L)	ParaHippocampal_L	Occipital	PHG(L)
Parahippocampal gyrus (R)	ParaHippocampal_R	Occipital	PHG(R)
Calcarine gyrus (L)	Calcarine_L	Occipital	CAL(L)
Calcarine gyrus (R)	Calcarine_R	Occipital	CAL(R)
Cuneus (L)	Cuneus_L	Occipital	CUN(L)
Cuneus (R)	Cuneus_R	Occipital	CUN(R)
Lingual gyrus (L)	Lingual_L	Occipital	LING(L)
Lingual gyrus (R)	Lingual_R	Occipital	LING(R)
Superior occipital gyrus (L)	Occipital_Sup_L	Occipital	SOG(L)
Superior occipital gyrus (R)	Occipital_Sup_R	Occipital	SOG(R)
Middle occipital gyrus (L)	Occipital_Mid_L	Occipital	MOG(L)
Middle occipital gyrus (R)	Occipital_Mid_R	Occipital	MOG(R)
Inferior occipital gyrus (L)	Occipital_Inf_L	Occipital	IOG(L)
Inferior occipital gyrus (R)	Occipital_Inf_R	Occipital	IOG(R)
Fusiform gyrus (L)	Fusiform_L	Occipital	FFG(L)
Fusiform gyrus (R)	Fusiform_R	Occipital	FFG(R)
Postcentral gyrus (L)	Postcentral_L	Parietal	PoCG(L)
Postcentral gyrus (R)	Postcentral_R	Parietal	PoCG(R)
Superior parietal lobule (L)	Parietal_Sup_L	Parietal	SPG(L)
Superior parietal lobule (R)	Parietal_Sup_R	Parietal	SPG(R)
Inferior parietal lobule (L)	Parietal_Inf_L	Parietal	IPL(L)
Inferior parietal lobule (R)	Parietal_Inf_R	Parietal	IPL(R)
Supramarginal gyrus (L)	SupraMarginal_L	Parietal	SMG(L)
Supramarginal gyrus (R)	SupraMarginal_R	Parietal	SMG(R)
Angular gyrus (L)	Angular_L	Parietal	ANG(L)
Angular gyrus (R)	Angular_R	Parietal	ANG(R)
Precuneus (L)	Precuneus_L	Parietal	PCUN(L)
Precuneus (R)	Precuneus_R	Parietal	PCUN(R)
Paracentral lobule (L)	Paracentral_Lobule_L	Parietal	PCL(L)
Paracentra lobule (R)	Paracentral_Lobule_R	Parietal	PCL(R)
Heschl gyrus (L)	Heschl_L	Temporal	HES(L)
Heschl gyrus (R)	Heschl_R	Temporal	HES(R)
Superior temporal gyrus (L)	Temporal_Sup_L	Temporal	STG(L)
Superior temporal gyrus (R)	Temporal_Sup_R	Temporal	STG(R)
Temporal pole (superior) (L)	Temporal_Pole_Sup_L	Temporal	TPOsup(L)
Temporal pole (superior) (R)	Temporal_Pole_Sup_R	Temporal	TPOsup(R)
Middle temporal gyrus (L)	Temporal_Mid_L	Temporal	MTG(L)
Middle temporal gyrus (R)	Temporal_Mid_R	Temporal	MTG(R)
Temporal pole (middle) (L)	Temporal_Pole_Mid_L	Temporal	TPOmid(L)
Temporal pole (middle) (R)	Temporal_Pole_Mid_R	Temporal	TPOmid(R)
Inferior temporal gyrus (L)	Temporal_Inf_L	Temporal	ITG(L)
Inferior temporal gyrus (R)	Temporal_Inf_R	Temporal	ITG(R)
Thalamus (L)	Left_Thalamus	Subcortex	THA(L)
Caudate (L)	Left_Caudate	Subcortex	CAU(L)
Putamen (L)	Left_Putamen	Subcortex	PUT(L)
Pallidum (L)	Left_Pallidum	Subcortex	PAL(L)
Brainstem	Brain-Stem	Subcortex	BST
Hippocampus (L)	Left_Hippocampus	Subcortex	HIP(L)
Amygdala (L)	Left_Amygdala	Subcortex	AMYG(L)
Nucleus accumbens (L)	Left_Accumbens	Subcortex	NAcc(L)
Thalamus (R)	Right_Thalamus	Subcortex	THA(R)
Caudate (R)	Right_Caudate	Subcortex	CAU(R)
Putamen (R)	Right_Putamen	Subcortex	PUT(R)
Pallidum (R)	Right_Pallidum	Subcortex	PAL(R)
Hippocampus (R)	Right_Hippocampus	Subcortex	HIP(R)
Amygdala (R)	Right_Amygdala	Subcortex	AMYG(R)
Nucleus accumbens (R)	Right_Accumbens	Subcortex	NAcc(R)
Cerebellar lobule I-IV (L)	Left_I-IV	Cerebellum	I-IV(L)
Cerebellar lobule I-IV (R)	Right_I-IV	Cerebellum	I-IV(R)
Cerebellar lobule V (L)	Left_V	Cerebellum	V(L)
Cerebellar lobule V (R)	Right_V	Cerebellum	V(R)
Cerebellar lobule VI (L)	Left_VI	Cerebellum	VI(L)
Cerebellar vermis VI	Vermis_VI	Cerebellum	VI-vermis
Cerebellar lobule VI (R)	Right_VI	Cerebellum	VI(R)
Cerebellar crus I (L)	Left_Crus_I	Cerebellum	XI(L)
Cerebellar crus I (R)	Right_Crus_I	Cerebellum	XI(R)
Cerebellar crus II (L)	Left_Crus_II	Cerebellum	XII(L)
Cerebellar vermis crus II	Vermis_Crus_II	Cerebellum	XII-vermis
Cerebellar crus II (R)	Right_Crus_II	Cerebellum	XII(R)
Cerebellar lobule VIIb (L)	Left_VIIb	Cerebellum	VIIb(L)
Cerebellar lobule VIIb (R)	Right_VIIb	Cerebellum	VIIb(R)
Cerebellar lobule VIIIa (L)	Left_VIIIa	Cerebellum	VIIIa(L)
Cerebellar vermis VIIIa	Vermis_VIIIa	Cerebellum	VIIIa-vermis
Cerebellar lobule VIIIa (R)	Right_VIIIa	Cerebellum	VIIIa(R)
Cerebellar lobule VIIIb (L)	Left_VIIIb	Cerebellum	VIIIb(L)
Cerebellar vermis VIIIb	Vermis_VIIIb	Cerebellum	VIIIb-vermis
Cerebellar lobule VIIIb (R)	Right_VIIIb	Cerebellum	VIIIb(R)
Cerebellar lobule IX (L)	Left_IX	Cerebellum	IX(L)
Cerebellar vermis IX	Vermis_IX	Cerebellum	IX-vermis
Cerebellar lobule IX (R)	Right_IX	Cerebellum	IX(R)
Cerebellar lobule X (L)	Left_X	Cerebellum	X(L)
Cerebellar vermis X	Vermis_X	Cerebellum	X-vermis
Cerebellar lobule X (R)	Right_X	Cerebellum	X(R)

Table 5. Table S4: Differences in link numbers within and between movie consensus modules that survive the 0.05 Benjamini-Hochberg FDR correction.

Module pair $i,j$	$p$-value	$\Delta W_{i,j}$	relative increase (+) / decrease (-)
VIS,VIS	0.0010	26204.4	+56%
VIS,CUN	0.0029	2228.9	+59%
DA,SAL	0.0007	-1120.0	-54%
DA,FP	0.0017	-1981.3	-35%

Equations

$$f(G,G^{\prime })=\frac{\text{number of common links in G and G′}}{\frac{1}{2}N(N-1)\rho }\in [0,1].$$

$$⟨\mathbf{W}{⟩}^{\text{rest}}=\frac{1}{13}\sum _{i=1}^{13}{\mathbf{W}}^{\text{rest},i},$$

$$⟨\mathbf{W}{⟩}^{\text{movie}}=\frac{1}{13}\sum _{i=1}^{13}{\mathbf{W}}^{\text{movie},i}.$$

$$S{I}_i=\frac{1}{13}\sum _{J\in \{1,...,13\}}\frac{|C_i^{\text{REF}}\cap C_i^{\text{J}}|}{|C_i^{\text{REF}}|}\frac{|C_i^{\text{REF}}\cap C_i^{\text{J}}|}{|C_i^{\text{J}}|}\in [0,1],$$

Full text

**This is the pre-peer reviewed version of the following article: Kujala, R., Glerean, E., Pan, R. K., Jääskel inen, I. P., Sams, M., Saramäki, J. (2016), Graph coarse-graining reveals differences in the module-level structure of functional brain networks. European Journal of Neuroscience, 44: 2673 2684, which has been published in final form at http://dx.doi.org/10.1111/ejn.13392. This article may be used for non-commercial purposes in accordance with Wiley Terms and Conditions for Self-Archiving. **

1 Introduction

Methods of network science are increasingly used for analyzing functional magnetic resonance imaging (fMRI) data [1, 2, 3, 4]. In this framework, fMRI time series are mapped to a network, where the nodes correspond to different brain areas and the links between nodes indicate dependencies between the blood-oxygen-level dependent (BOLD) time series of these areas. Whenever the dependency between the BOLD signals corresponding to two brain areas is strong enough, they are thought to be functionally related.

Network analysis has revealed many insights on the functional structure of the brain. For instance, on the scale of network nodes, the use of various network centrality measures has allowed consistent detection of certain hub regions [5]. At the network level, the entire structure of functional brain networks has been found to be of the ‘small-world’ type [6, 7, 8], meaning that the average number of steps required to reach any other node in a network is low while the network is locally clustered. In addition to revealing such general properties, the network framework has also been used for investigating how brain dynamics depend on different stimuli [9], mental health [10, 11, 12], and age [13].

Despite being widely applied, the network approach to brain functional networks can still be considered as somewhat immature, and methodological variations persist in the literature [14, 15]. In particular, there is no standard set of brain areas to be used as network nodes; rather, different studies use different definitions of what constitutes a node [14]. Approaches for defining network nodes range from the use of anatomical atlas-based Region-Of-Interest (ROI) parcellations to the use of the original fMRI imaging voxels of a few mm3 in size [14]. Typically, the number of nodes in atlas-based definition schemes is of the order of $10^{2}$, while the number of nodes in voxel-based node definition schemes is of the order of $10^{3}-10^{4}$. These differences between node definitions make comparison of results across studies challenging.

The use of small number of nodes is understandable as it makes network analysis easier, and a low number of nodes even allows meaningfully visualizing the networks. However, such low resolution also means that the nodes of a network may cover multiple functionally specific areas. As the representative BOLD signal for each node is typically computed as the average over the node’s spatial extent, this can lead to significant loss of information and node-level signals that are not representative of true function. The recommended remedy is to instead use individual voxels as the network nodes [14]. Given that with high-field (e.g. 7 Tesla) fMRI it has been possible to reach a spatial resolution smaller than 1 mm3  [16], the volumes that voxel-based nodes represent can be expected to decrease in the future, while their number increases.

Constructing and analyzing networks that consist of voxel-level nodes is challenging because of the large number of nodes and links. While computing node-wise centrality measures or global network characteristics may be possible for voxel-level networks, understanding the the overall organization of a network’s links becomes difficult. Because the number of potential links is of the order of millions, mere visualization of such networks is extremely challenging, and e.g. comparing two networks on a link-to-link basis becomes meaningless.

One approach for going around this problem is to look at network structure on an intermediate level between the level of nodes and the level of the entire network. This intermediate level is often approached by splitting the network into modules, i.e. groups of nodes that are densely connected internally but have few connections between them [17]. Modules are typically discovered using stochastic algorithms that partition the network into non-overlapping node groups. For examples of applications of module detection to functional brain network analysis, see [18, 19, 20, 3, 21, 22].

Any observed modular structure of functional brain networks may arise from experimental conditions (different stimuli), or reflect more persistent, underlying features of the subject brain. Some recent studies have investigated how network modules differ between healthy controls and patients suffering from schizophrenia [11], patients with autism spectrum disorders [12], or patients who are comatose [10]. However, due to various intricacies in the detection of network modules and the difficulty of comparing different module partitions, results have been difficult to interpret and their statistical significance has remained elusive. Although one can statistically approach the stability of modules [23, 12] or the similarity of partitions between groups of networks [11], verifying the significance of specific differences in modules remains an open problem.

Therefore, there is a need for appropriate methods for analysing and comparing functional network structure at the intermediate level of modules. We argue that instead of focusing on how the network modules themselves differ between groups of networks, it is more fruitful and methodologically sound to assess the differences in the numbers of connections between and within a fixed set of modules that is used as a frame of reference for all groups. The boundaries of these modules are defined using a specific network or group of networks, and applied as such to the rest of available networks. This coarse-graining approach allows transparent, statistically verifiable investigation of module-level differences between networks.

Below, we first demonstrate the difficulties of comparing network modules with the help of a toy example, and proceed to show how the coarse-graining approach overcomes these problems. Then, to show the applicability of our approach in practice, we apply it to networks constructed from fMRI data recorded for 13 subjects during rest and movie viewing.

2 Comparing functional networks: coarse-graining and alternative approaches

Perhaps the simplest approach for studying differences in the link structure of functional brain networks is to investigate how the existence or weight of individual links differs between groups of networks. In this case, assessing the statistical significance of the observed differences is straightforward. However, with voxel-based functional brain networks that have thousands of nodes and millions of links this approach is impractical for several reasons. First, the number of statistical tests that need to be performed reaches millions, requiring a large amount of computational resources. Second, the spatial locations of functionally specific parts of the brain differ across people even after imaging data have been transformed into standard coordinates. Consequently, there is no one-to-one correspondence between specific nodes and links across subjects. Third, visualizing all differences between groups of networks becomes impossible.

Because of the above, some studies have focused on quantifying differences between networks at the level of modules that have been discovered using stochastic network partitioning algorithms [11, 12]. It is worth noting that there is no universally agreed definition for network modules [17] – rather, each algorithm introduces its own definition of a network module that dictates how it partitions a network and how modules should look like. Subsequently, different algorithms give different results. Therefore, proper interpretation of module partitions requires a profound understanding of the underlying mathematical ideas and of the actual implementation of the algorithm.

If interpreting what module partitions signify can be difficult, it is even far more difficult to interpret differences between partitions. Often, small differences in a network’s link structure give rise to large differences in the partitioning of the network into modules. Yet, at times, major structural differences result in little or no differences in the discovered modules. We illustrate this with Fig. 1. In this toy example, we consider two networks, A and B, that have a clearly modular structure. Both networks share the same set of nodes and their link structures are almost the same, with 4 links being different wired differently in B. We then apply the Louvain algorithm [24], one of the most popular methods for detecting network modules, to the two networks separately. While the removal of three links between the yellow and brown modules of Network A results in no differences in the modules detected in network B, the addition of only one extra link between the two smaller green and violet modules in Network A forces them to merge into one large blue module in Network B. This simple example illustrates that it is not straightforward to infer the underlying differences in the overall organization of links by comparing only the network partitions: the cause and the apparent effect size can be disproportionate.

When dealing with experimental, noisy fMRI data the challenges are even greater. The randomness in module partitions that arises from the stochasticity of the algorithms is further amplified by the noise in the data, which obscures real structural differences between networks. In addition, there are no statistical frameworks for directly assessing the statistical significance of differences between partitions – e.g., is the merging of the two network A modules into the the larger blue network B module statistically significant or just due to chance? Thus, even when the structure of network modules may shed some light on the functioning of the brain under different conditions, drawing definite conclusions from differences between network partitions is very difficult.

We have seen that directly comparing module partitions is difficult, as is making comparisons at the fine-grained level of individual links. However, as we will show below, combining both points of view yields fruitful results. In particular, we argue that first producing a fixed set of modules, whose boundaries will be the same in all networks, and then comparing inter- and intra-module connections across networks reveals their differences more clearly. Conceptually, this corresponds to network coarse-graining, where each module corresponds to a node of the coarse-grained network, and the number of links between two modules in the original network corresponds to the weight of the link in the coarse-grained network. The number of internal links within each module is taken into account as the weight of the self-link (connection from the node to itself) of each module-node of the coarse-grained network. The usefulness of the coarse-graining approach is demonstrated in Fig. 1 with the same toy networks as earlier. Notably, the comparison of coarse-grained networks reveals differences in a more transparent way than attempting to compare modules. In addition, when groups of networks are to be compared, the coarse-graining approach allows statistical testing of the differences in the mean number of links within a module or between two modules.

3 Materials and Methods

Participants

The participants were 13 healthy native Finnish speakers (ages 22-43 years, 2 females, 2 left-handed, no neurological or psychiatric history, no hearing impairments, normal vision). The ethical committee of the Hospital district of Helsinki and Uusimaa granted permission for this study which was conducted in accordance with the guidelines of the declaration of Helsinki. Each subject gave written informed consent prior to participation.

Stimulus paradigm

The stimulus used in this paper has also been used previously [27, 28] and consisted of an edited version of the Finnish movie ‘The Match Factory Girl’ (Aki Kaurismäki 1990). The film was projected on a semi-transparent screen behind the subject’s head and the audio track was delivered via plastic tubes through porous earplugs. Each subject went through three sessions with the following order: resting state (15 minutes, 450 volumes), free viewing of the film (22 min 58s, 689 volumes), resting state (15 minutes, 450 volumes). In this study we analyze data recorded during rest before viewing the movie, and during movie viewing. After preprocessing, the movie session was truncated to match the length of the rest session to avoid any biases from different scan durations.

Data acquisition

MR imaging was conducted on a 3.0T GE Signa Excite MRI scanner, with a quadrature 8-channel head coil. A total of 29 functional gradient-echo planar axial slices (thickness 4 mm, 1 mm gap between slices, in-plane resolution 3.4 mm $\times$ 3.4 mm, imaging matrix 64 $\times$ 64, TE 32 ms, TR 2000 ms, flip angle 90∘,). T1-weighted images were also acquired (TE 1.9 ms, TR 9 ms, flip angle 15∘, SPGR pulse sequence) with in-plane resolution of 1 mm $\times$ 1 mm, matrix size 256 $\times$ 256 and slice thickness 1 mm with no gap.

Preprocessing

Preprocessing of the fMRI data was carried out with FSL (release 4.1.6 www.fmrib.ox.ac.uk/fsl). The first 10 volumes of each session were discarded from the analysis. Motion correction was performed with McFlirt and the data were spatially smoothed using a Gaussian kernel with 6 mm full-width half maximum, and high-pass filtered with a 100 s cutoff. Functional data were co-registered with FLIRT to the anatomical image allowing 7 degrees of freedom. Furthermore, the data were registered from the anatomical space to the MNI152 2mm standard template (Montreal Neurological Institute), allowing 12 degrees of freedom. The signal was bandpass filtered with a passband of 0.01–0.08 Hz in accordance with standard functional connectivity procedures. To control for motion artefacts, motion parameters were regressed out from the data with linear regression (36 Volterra expansion based signals, see Ref. [29]). As it is known that head motion affects connectivity results, we controlled for motion with framewise displacement [30]: all subjects had at least 95% of time points under the suggested displacement threshold of 0.5 mm. For this reason, we decided not to use the scrubbing technique and utilize all time points. When looking at the individual mean frame-wise displacement, there was no significant difference between conditions ($p$=0.2732). Finally, to further control for artefacts, voxels at the edge between brain and skull where the signal power was less than 2% of the individual subject ’ s mean signal power were excluded from the analysis. This resulted in 5562 6-mm isotropic voxels of brain grey matter covering the whole cerebral cortex, subcortex and cerebellum. After removing the first and last 15 data points due to bandpass filtering artifacts [29], we obtained for each 6-mm voxel a BOLD time-series with 410 time points corresponding to a duration of 13 min 40 s.

Network construction

The functional dependency of two voxels $i$ and $j$ can be measured in many ways, given their BOLD time series $s_{i}(t)$ and $s_{j}(t)$ [31]. As there is no consensus on the best measure, we opt for simplicity and use the Pearson correlation coefficient, which has been shown to capture a major proportion of pairwise dependencies in fMRI data [32]. For each subject and condition, we then compute a correlation matrix R, whose elements $r_{ij}$ are the estimated Pearson correlation coefficients between voxels $i$ and $j$. Given that we have 2 different conditions and 13 subjects, this yields 26 correlation matrices in total.

There are several ways of constructing networks from such matrices by thresholding them so that only chosen elements remain. One common approach is to use a constant threshold value, so that only nodes pairs whose correlation coefficients (link weights) exceed this value are connected by a link. Another typical approach is to include a fixed fraction of the strongest links in the network. Here, we adopt the latter approach as it has been shown to provide more stable estimates of various network measures [15].

We construct networks from correlation matrices as follows: for each matrix, we begin with an empty network, where each node corresponds to a voxel. We then first compute the maximal spanning tree (MST) of the correlation matrix, and insert the corresponding links to the network. As the MST connects all nodes, this guarantees that no nodes or groups of nodes remain isolated in the network; isolated modules would cause technical difficulties in the later stages of our pipeline. Next, we sort all correlation coefficients, and insert links corresponding to the strongest positive coefficients until the network contains $\frac{1}{2}N(N-1)\rho$ links in total, where $\rho\in[0,1]$ is a pre-defined network density. As the end result, we obtain 26 undirected, unweighted networks that all share the same set of nodes, and have the same number of links.

Selection of network density

Choosing the fixed density for the thresholded networks is not a straightforward problem. There are no commonly-accepted criteria for choosing an optimal density. If the density is very low, too much information is discarded and features of interest may remain hidden. On the other hand, the presence of too many links may obscure relevant structures. Sometimes this problem can be overcome by investigating network structure across different network densities [11, 33]. However, carrying out detailed structural analysis of a large number networks of different densities quickly becomes overwhelming. Then, selecting a reasonable specific network density may be a better option. This is also the case in this study, where we compare the module-level structure of groups of networks in detail.

To guide our choice of network density, we investigate the similarity of pairs of networks as a function of the network density. If two networks share the same set of nodes and have the same number of links, the most straightforward approach to measure their similarity between is to count the number of shared links. This is the approach we adopt. Given two networks $G$ and $G^{\prime}$ that both have $\frac{1}{2}N(N-1)\rho$ links, we monitor how the fraction of links $f(G,G^{\prime})$ common to the networks changes with the network density $\rho$. The fraction of common links $f(G,G^{\prime})$ is defined as the number of shared links divided by the total number of links in one network:

[TABLE]

In Fig. 2 we show $f$ as a function of the network density when averaged over pairs of networks, so that each pair represents (i) the same subject in different conditions, (ii) different subjects in the rest condition, (iii) different subjects in the movie condition, and (iv) different subjects in different conditions. For all cases, the fraction of common links first increases until $\rho\approx 0.1\%$; one might envision that at this density, a common “backbone” shared by networks is well captured. Then $f$ decreases until $\rho\approx$ 2%, after which it begins to monotonously increase as the networks become denser and more and more links are necessarily shared.

Based on the observed behaviour of $\rho$, we pick the value $\rho=2\%$ to be used in all subsequent analyses, as it appears to maximize non-trivial variation between network pairs (of course, there is more variation for excessively low $\rho$). This value was also used in Ref. [11] for networks with smaller numbers of nodes. In our data, the 2% network density translates to 309 303 links in each thresholded network; for different networks, this on average corresponds to a correlation coefficient threshold of 0.56$\pm$ 0.04 (std) (no statistically significant difference between conditions).

Fig. 2 also reveals some further insights. First, the network similarity $f(G,G^{\prime})$ is remarkably higher for networks corresponding to the same subject in different conditions than for different subjects in the same condition. This indicates that individual variation dominates over differences caused by different stimuli. Subsequently, it is essential to take the paired nature of the data into account when validating any results statistically. Second, when the similarity of networks of different subjects is assessed, network pairs corresponding to the movie condition are seen to be more similar than network pairs corresponding to rest condition or different conditions. This is expected, as the viewing of a well-directed movie stimulus has been found to synchronize the subjects’ brains [34], which results in increased functional connectivity compared to the similarities of resting state networks arising from shared functional anatomy and connectivity.

Computation of condition-wise consensus partitions

Given the networks for all subjects S1-S13 and conditions, our next target is to compute representative network partitions for each condition. As outlined in Fig. 3, this is realized in two steps: First, we partition the networks using a popular partitioning algorithm, the Louvain algorithm [24], which is based on modularity optimization [25]. In more detail, we run the stochastic Louvain algorithm 100 times for each network and select the network partition with the highest value of modularity. Then, for each condition, we summarize all 13 selected partitions with the help of the MCLA meta-clustering algorithm [35] which yields one representative network partition as an output. In addition to the network partitions, MCLA requires the user to input a parameter determining the upper bound for the number of modules in the consensus partition. In this study, the value of this parameter was set to the median number of modules in the 13 partitions.

To summarize, our pipeline transforms the subject-specific networks corresponding to one condition to a single representative consensus partition $\mathcal{P}$ consisting of $m$ modules $C_{1},...,C_{m}$ which in turn are sets of network nodes such that each node belongs to exactly one module.

Comparing groups of coarse-grained networks

The network coarse-graining process briefly introduced in Fig. 1 is defined as follows: Given a network $G$ and a partition $\mathcal{P}$ consisting of a set of modules $\{C_{1},C_{2},...,C_{m}\}$, the coarse-graining process yields a matrix $\mathbf{W}$ that has the following properties: The non-diagonal matrix element $\mathbf{W}_{i,j}$ represents the total number of links between the nodes of $C_{i}$ and $C_{j}$. Similarly, each diagonal element $\mathbf{W}_{i,i}$ represents the number of links within module $C_{i}$.

To evaluate differences between experimental conditions, we first coarse-grain each subject’s network into its matrix representation $\mathbf{W}^{\text{condition,i}}$, where $i$ stands for the index of the subject. Then, we average the coarse-grained matrix representations of all 13 subjects over each condition:

[TABLE]

and

[TABLE]

Then, by investigating the elements of the mean difference matrix $\Delta\mathbf{W}=\langle\mathbf{W}\rangle^{\text{movie}}-\langle\mathbf{W}\rangle^{\text{rest}}$, we can quantify the level of differences in the numbers of connections between and within modules. Thus, the value of $\Delta\mathbf{W}_{1,2}$ indicates how many more connections there are on average between modules 1 and 2 in the movie condition than in the rest condition.

To test the statistical significance of our findings, we perform paired permutation tests separately on each of the matrix elements. As we have a limited number of samples, we use the full permutation distribution yielding in total $2^{13}=8192$ different permutations. All $p$-values we report are two-sided, and we correct for multiple comparisons using the original Benjamini-Hochberg (BH) FDR correction [36].

Code

The Python code used in our analysis is freely available at http:github.com/rmkujala/brainnets.

4 Results and Discussion

Consensus modules are similar to previously reported resting-state modules

The consensus modules computed for the rest and movie-viewing conditions are shown in Fig. 4A with different colors on the cortical surface. A browsable display of each module is available at NeuroVault http://neurovault.org/collections/1080/ [37], where the modules are weighted by their consistency as measured by scaled inclusivity [38] (see Supporting Information for more details). In the center of Fig. 4A, the matching of the rest and movie consensus partitions is also visualized as an alluvial diagram [26].

For the rest condition, we identified 11 consensus modules (Fig. 4, left hand side of alluvial diagram, from bottom to top):

1. Limbic (LIM) subcortical midbrain structures;
2. Cerebellum/ventro-temporal (CRBL/VT);
3. Default mode (DM);
4. Precuneus (PCUN);
5. Visual (VIS);
6. Auditory (AUD);
7. Salience (SAL);
8. Fronto-Parietal (FP);
9. Visual-extrastriate (VISx);
10. Sensorimotor (SM);
11. Language (LAN). For the movie condition, we obtained 10 modules:
12. Ventro-temporal limbic (VTL);
13. Cerebellum (CRBL);
14. Default mode (DM);
15. Cuneus (CUN);
16. Visual (VIS);
17. Auditory (AUD);
18. Salience (SAL);
19. Fronto-parietal (FP);
20. Dorsal attention (DA);
21. Sensorimotor (SM). For details on how the modules were assigned a label, please see Supplementary Information.

There is a good agreement in the literature on the module-level structure of resting networks, which have been identified using various methods such as multidimensional clustering [39], Infomap graph clustering [3], and independent component analysis [40]. Overall, the resting state consensus modules we obtained are in line with the previously reported resting state modules (see Table S1).

However, there is no general agreement on the module structure during movie viewing, or more generally, during a task. While one study has found that task and rest are highly similar [41], another study has found remarkable differences in subcortical, limbic regions as well as primary sensory and motor cortices [42]. In our case, overall, the movie consensus modules are similar to the resting state modules, and e.g. the dorsal attention module that has been previously reported in resting state studies [39, 41] is even better identified amongst in our movie consensus partition than in the rest consensus partition (see Tables S1 and S2). Interestingly, we also identified a VTL subnetwork present only during movie viewing in agreement with [12], possibly suggesting stronger functional couplings between brain areas involved in the processing of social and emotional events in the movie.

Differences in condition-wise consensus modules are difficult to directly interpret

As shown in Fig. 4A, the consensus network partitions obtained for the movie and rest conditions are broadly speaking similar: for most rest modules, there is a clear counterpart among the movie modules. At the same time, almost all rest modules overlap with multiple movie modules (and vice versa) – there are no simple relationships such as one module splitting into two, and the varying amount of overlap between modules results in a diagram that is not straightforward to interpret beyond the clear matches.

As there are no statistical frameworks that can be used for measuring the significance of the relationships between multiple modules, the differences in the alluvial diagram lack statistical validation. Some insights into the significance of the transitions in the alluvial diagram could be obtained by investigating the consistency of the modules [23, 12]. However, these methods do not directly assess the significance of individual splits and joinings of modules between partitions. Thus, it remains challenging to draw conclusions about network differences based on the alluvial diagram alone.

Network coarse-graining allows simple, statistically verifiable interpretations

In Fig 4B, we show the matrix representations of the average coarse-grained networks corresponding to both conditions, as well as the coarse-grained difference network, where the resting-state consensus modules have been used as the basis for coarse-graining. As expected, the matrices representing the coarse-grained networks have high values on their diagonals, indicating that the density of links within consensus modules is higher than between them. Off-diagonal element provide an overview of how strongly the modules are connected.

In the coarse-grained difference matrix, we observe multiple elements that survive the 0.05 Benjamini-Hochberg FDR correction. These are also listed in Table 1. All surviving elements are positive, indicating that there are more connections between the modules in the movie networks. In particular, the VIS module displays significantly more external as well as internal connections in the movie condition. The number of connections between the AUD and DM modules is also increased in the movie condition. The coarse-graining method thus succeeds in highlighting task-driven changes at visual areas as well as inferior temporal structures. For similar coarse-graining results where the movie modules are used as the frame of reference, please see Supporting Information.

There are of course some similarities between the coarse-grained difference network and the alluvial diagram presented in Fig. 4. As an example, in both it is seen that the resting-condition VIS module and part of the resting-condition VISx module merge to form the larger VIS module in the movie condition. A simple explanation for this would be that there are more connections between the VIS and VISx rest modules in the movie condition. However, as motivated in Sec. 2, this is not the only possible reason, and further, the statistical significance of the observation cannot be deduced using the module-matching/alluvial-diagram approach. However, the coarse-grained difference matrix clearly indicates that there is a statistically significant increase (47%) in the number of links between the VIS and VISx modules in the movie condition. Thus, while the alluvial diagram can be used for formulating hypotheses on changes in network structure, the coarse-graining process allows verifying that the observed differences are not due to random chance.

In addition to assessing differences in the connectivity between modules, the coarse-graining approach also allows investigating the internal connectivity of modules. This cannot be directly done by comparison of matched modules; further, if the nodes are from the beginning defined as larger entities (e.g. using anatomical atlases), this information is lost.

Methodological considerations

In this work, we have applied the introduced coarse-graining approach only to the analysis of undirected and unweighted networks. However, the approach itself is more general and could be easily extended to the analysis of weighted networks that also take into account the strength of correlation between voxels, or directed networks that depict causal relationships between brain areas. The coarse-graining approach could also be applied in clinical settings by comparing the module-level differences between different groups of individuals, extending the previous attempts of comparing brain modules [11, 12].

There are many choices to be made in deciding on the particulars of the module detection pipeline, especially the choice of the module detection algorithm (and its parameters, if any). Naturally, these choices will also affect the exact outcome. However, we expect that regardless of the exact way of defining modules, the benefits of using coarse-graining over comparing partitions remain. This is because the problems in interpreting the differences between partitions are universal, e.g. for any module detection method, there is nearly always a borderline case where the existence or absence of one single link affects boundaries between modules. Naturally, the choice of a module detection method also affects the resulting coarse-grained networks; however, when comparing across conditions and using the modules of one condition as the frame of reference, the differences should still be straightforward to interpret. Also, which of the conditions is used as the basis for coarse-graining should not be a critical choice: e.g.  the observed differences in the network structure remain relatively similar for rest and movie modules as basis (see Supporting Information).

It is also essential to point out that, most likely, there is no perfect partition of the brain into distinct areas. Thus, even though different module detection algorithms yield different network partitions, these may still capture at least some meaningful aspect of the network’s organization. Instead of considering this multiplicity of possible partitions as a problem, it is rather an opportunity, as coarse-graining networks using different partitions can actually turn out to be very useful. In the same way as one can deduce the shape of a 3D object from its 2D projections, a network’s structure can be better understood by investigating its different coarse-grained representations. Therefore, coarse-graining could be useful even using an anatomical brain atlas as the basis, as the atlas provides a frame of reference that is well known to researchers in the field.

5 Conclusions

We developed a coarse-graining method to analyze differences in the modular structure of functional brain networks during rest and task. The coarse-graining approach focuses on the differences in the connectivity between and within larger brain areas without sacrificing the spatial accuracy of fMRI data already at the network construction stage. The method yields results that summarizes differences in connectivity on the module level, using a set of modules as a frame of reference across groups or conditions. In contrary to some alternative approaches studying differences in the module-level connectivity of the human brain, the results produced by our method are both easy to interpret and verify statistically – they allow to “see the forest for the trees”. Because data on the structural and functional human connectome are becoming more and more detailed, we believe such methods will play an increasingly important role in understanding the module-level structure of functional networks.

6 Acknowledgments

Computing resources provided by the Aalto Science IT project are acknowledged. JS acknowledges financial support from the Academy of Finland, project n:o 260427. EG acknowledges aivoAALTO Project Grant from the Aalto University, doctoral program “Brain & Mind”. We thank Onerva Korhonen for comments on the manuscript.