Temporally intermittent communities in brain fMRI correlation networks

It is widely agreed that the human brain is organized as a system of segregated modules that reside in separate regions and, through coordinated integration, support different cognitive functions. Through recent breakthroughs in modeling the activity of the brain, it has been demonstrated that each such module can participate in multiple so-called functional networks – networks of brain regions that activate in synchrony during specific types of cognition. If we model the brain as a temporal network, by representing brain regions as nodes and correlations within activity-windows at different times as links, we can formulate the task of finding functional networks as a community detection problem. In spite of this, however, relatively little attention has been given to solving this problem using recently developed techniques for temporal community detection. In this paper, as a proof-of-concept, we apply a novel technique for community detection in temporal networks to a dataset of fMRI measurements from 100 healthy subjects undertaking a working memory task with intermittent fixation (or resting-state) periods. We show that this method recovers two distinct communities that are shared between subjects: one that activates during the fixation period, and another that activates during a period associated with high cognitive load.


Introduction
Global brain activity exhibits a surprising level of organization in both space and time.The consilience of evidence from studies using various imaging technologies and computational methods establishes the existence of so-called "functional networks": sets of brain regions that activate in synchrony and support a vast repertoire of brain functions (Gusnard and Raichle 2001;Fox et al. 2005;Park and Friston 2013).There are many well studied methods for detecting these functional networks, involving non-invasive imaging techniques such as functional magnetic resonance imaging (fMRI), electro-and magnetoencephalography (EEG and MEG).Central to these methods is the inference of synchrony between different brain regions, commonly gauged by measuring mutual information or correlation strength between time-series representing activity in separate regions.The simplest approach is to query the similarity between regions-of-interest (ROI) (Biswal et al. 1995;Lowe et al. 1998;Greicius et al. 2003;Fox et al. 2005), which are typically parcels of a predetermined segmentation of the brain into distinct regions, called a parcellation or an atlas (Jenkinson et al. 2012;Brodmann 1909;Fan et al. 2016).
A key feature of ROI-based methods is that they require a predefined parcellation of the brain as input.There are multiple limitations due to this strategy: (1) Predefined parcellations may not fit individual subject data well due to processing characteristics and population mismatch.(2) There currently exists many different, but mutually inconsistent, atlases (Bohland et al. 2009).(3) A researcher must make a design decision about which parcellation to use, a choice which can severely impact the results.To remedy these limitations, various unsupervised parcellation methods exist, based on -but not limited to -mixture models, k-means clustering, hierarchical clustering, spectral clustering, principal component analysis (PCA) and independent component analysis (ICA) (Thirion et al. 2014).
We employ ICA to obtain an unsupervised representation of the brains functional activity across a set of components.ICA is often used as a dimensionality reduction method before dynamic functional connectivity analysis (Allen et al. 2014;Vidaurre et al. 2017), investigating repeating patterns of time-resolved connectivity between independent components.Furthermore, ICA has previously been used to recover spatial functional networks and is used in many applications (Kiviniemi et al. 2003;McKeown et al. 1998;Beckmann et al. 2005), where some attention is given to identifying temporally fluctuating potentially overlapping functional networks (Smith et al. 2012).The domain of network neuroscience has exploited techniques from network science to uncover structure in connectomics data (Bassett and Sporns 2017).In the context of dynamical functional connectivity analysis of fMRI data, the multiplex modularity framework of Ref. (Mucha et al. 2010) has previously been used to characterize modular structure in ROI based time-resolved dynamic functional connectivity (Bassett et al. 2011;Bassett et al. 2013;Khambhati et al. 2018).However, an important limitation of this method -and to our knowledge all other methods for community detection in temporal networks -is the assumption that dependencies between individual nodes across time can be modeled as the dependence between entire layers.Real systems with multiple asynchronous concurrent events must have varying dependencies between the same layers, and by ignoring these node-level dependencies important temporal dynamics are washed out.
In this paper, we apply the novel neighborhood flow coupling (NFC) method (Aslak et al. 2018) for temporal community detection to a temporal correlation network of independent components (ICs) inferred with ICA on fMRI data (Barch et al. 2013).The temporal network is represented as a sequence of static networks, or layers, that share nodes, but have varying links depending on the inter-IC correlations that the links represent.NFC, in contrast to other methods, inter-couples layers in this temporal network based on link structure at the node-level.Specifically, the time-states of each node are coupled proportional to the similarity of their neighborhoods.This creates a large connected network which we can partition into modules, or communities, using any community detection algorithm.Here, we use Infomap (Edler and Rosvall 2014;Rosvall et al. 2009).This framework enables us to capture both the communities, which may correspond to functional networks in the brain, and their very fine-grained temporal evolution.

Results
We use a dataset from the Human Connectome Project (HCP) (Barch et al. 2013), with fMRI measurements from 100 subjects undertaking an n-back working memory task.
The task consists of three experimental blocks, 2-back, 0-back and fixation, with varying cognitive load.In the 2-back condition, subjects are presented with images and are asked to press a button if it matches the image presented two instances back.In the (easier) 0back condition, they are asked to press a button when a target image is presented.During the fixation periods, no instructions are given and subjects lie passively with eyes open in a resting-state.For each subject and experimental block, we can measure a performance accuracy as the fraction of correct button presses.The duration of an experiment is just under 5 min, and measurements are taken in 405 ticks with frequency 0.72 s.
The dataset describes activity across hundreds of thousands of voxels in each subject.To reduce noise as well as the network size defined in terms of pairwise correlations we rely on techniques for describing the brain signal using fewer components.We use independent component analysis (ICA) as a method for dimensionality reduction.Using ICA we produce a set of 50 independent components (ICs) each described as a spatial map consistent over subjects (Additional file 1: Figure S1), and an accompanying time-course of component activity.ICs are under few spatial constraints and have a whole-brain spatial distribution.In practice, however, ICs are often peaked around a smaller localized region and components rarely overlap.
We represent component activity as a temporal network, by making the abstraction that ICs of the brain are physical nodes i, and correlations between their time-courses within a window are links at times α such that all intra-layer links are represented by adjacency matrices W α ij .We use a sliding window of size w = 25 (each time-step is 0.72 s).Thus, the 405 original time-ticks in the data results in a temporal networks of 381 layers.Figure 1 gives a schematic overview of the experimental design and the processing procedures.Note that we use a fairly short window for correlating time-courses compared to what is recommended in the resting-state dynamic functional connectivity literature (Leonardi and Van De Ville 2015;Zalesky and Breakspear 2015).This is necessary, as the window-duration must remain shorter than the time-scale of the two tasks and fixation periods, since, if the window exceeds this time-scale, we cannot compare the temporal dynamics of communities to the experiment design.There is, however, evidence from task-based fMRI studies showing that shorter window-lengths have the necessary stability to detect changes in ongoing cognition (Gonzalez-Castillo et al. 2015).The inter-layer adjacency matrix D αβ i that connects states at times α Fig. 1 Experimental design.a Schematic illustration of the order in which subjects undertake memory tasks during the experiment.We apply the same smoothing procedure to this illustration and the data.(1) The initial order of tasks and fixation (405 ticks).( 2) The fMRI method measures blood oxidation rather than the actual neural activity, which offsets the measurement.To offset the sequence by the same amount we convolve it with the hemodynamic response function (Glover 1999) (see (b)). ( 3) We apply a mean filter of width w = 25 across the sequence.The final sequence is of length 381 ticks and β of each physical node i is inferred with neighborhood flow coupling (NFC) (Aslak et al. 2018).In brief, NFC couples temporal states of each physical node by a strength that is proportional to the neighborhood similarities of each state-pair.If at two separate times a node has similar neighborhoods, there will be a strong link between those two time-states of the node.Consequently, time-intermittent communities (communities that repeatedly disappear and re-appear) will have many strong links between them because, at each recurrence, the member nodes have similar neighborhoods.Thus, inferring the inter-layer adjacency matrix D αβ i using NFC creates a temporal network that, when used as input to a community detection algorithm, can reveal intermittent communities.
We use the Infomap algorithm for community detection (Rosvall et al. 2009) which, in brief, partitions a network into modules that have maximal flow within them.Infomap minimizes a cost function called the Map Equation, which is an analytical lower bound on the per-step description length of a random walker on the network given a partition of nodes into modules.On an abstract level, it considers the network as substrate for information flow, which by analogy suits our current application well.We use the software implementation of Infomap which has built-in support for NFC (Edler and Rosvall 2014).

We find communities that are shared across subjects
One way that we can search for communities which-if they exist-are guaranteed to be universal across subjects, is to time-concatenate the networks of multiple subjects to build a larger temporal network and then find a single community detection solution.This strategy enables a functional network which activates in multiple different subjects, to be detected as a single community.
We start by performing this analysis for the ten subjects with the highest 2-back task performance accuracy, and find two highly distinct communities (Fig. 2).One community (here, community 1) appears regularly during fixation for nearly all subjects.Another, larger, community (community 0) seems to activate during task performance and turn off completely during fixation for most subjects.Communities 0 and 1 are, while not entirely disjoint, distinctly different, both in their component distributions and in their temporal activity series, where they exclude each other for many subjects (Fig. 2c).The community detection method does find more communities, but none that are very large and regular across subjects or otherwise distinguishable from noise.In Additional file 2: Figure S2 we display communities 2 and 3 to support this insight.Furthermore, since Infomap indexes communities by descending size, it is highly unlikely that any of the remaining smaller communities contain signal.
Next, we perform the same analysis on the ten subjects with the lowest 2-back task accuracy, and find that both communities 0 and 1 are recovered with largely the same component profiles (Additional file 3: Figure S3).For community 0 it appears that the time activity series is less modulated by the course of the experiment.Similarly, we perform the analysis on all 100 subjects and recover again communities 0 and 1 with almost identical component distributions (Additional file 2: Figure S2).This repeated detection of communities 0 and 1 strongly suggests that, for the chosen parameter set (see "Materials and methods" section for discussion on parameters), these communities are non-random and prominent in our population.

Task-synchrony reveals two prominent functional networks
From visual inspection of Fig. 2a, it is clear that the activity patterns of communities 0 and 1 are modulated by the experimental tasks that subjects undertake.We now quantify this correlation.Using the community detection solution obtained from all 100 subjects together (Additional file 2: Figure S2), we gauge the level of synchrony between communities and tasks.For a given subject, we measure synchrony between a community and a task as the Pearson correlation between the community's activity level over time and a vector representing when a task is being undertaken (specifically, the proportion of a given task inside each time-window, see Fig. 1a3).We plot the distribution of communitytask synchrony within the population for each combination of task (fixation, 2-back and 0-back) and community (0 and 1).We also scatter-plot synchrony values against the task performance accuracies to investigate if they co-vary.Synchronies that are insignificant after multiple comparisons correction (100 hypotheses) are rendered in gray.
We observe that, for almost all subjects, community 1 is in strong positive synchrony with the fixation period, while community 0 is in strong negative synchrony with the fixation period (Fig. 3a-b).The exact opposite is observed for the 2-back task, where we find that most subjects activate community 0 and shut down community 1 (Fig. 3c-d).Interestingly, there is no such trend in the synchrony distributions for the 0-back task, Fig. 3 Task performance versus community-task synchrony.Gray regions represent insignificant synchrony (measured as Pearson correlation) between tasks and communities after multiple comparisons correction (100 hypotheses).First-axis values in (c-f) are task performance accuracies of the respective task, while the average of 0-back and 2-back is used (a-b).In no case does task performance correlate significantly with community-task synchrony (6 hypotheses) revealing that community 0 is not a general task-related community, but strictly a 2-backtask related community.
We plot the community-task synchrony values for each subject against their task performance accuracy score (recall: fraction of correct presses) and measure the Pearson correlation between these -since it is reasonable to hypothesize that if a community represents brain activity related to a specific task, then the performance accuracy of that task should modulate the dynamics of the community.However, we measure no significant correlations between task performance accuracy and community-task synchrony.Finally, as a way to correct for individual subject ability, we test the relation between performance gain (2-back performance divided by 0-back performance) and synchrony, yet find that it is not significant.
In reviewing the core components of communities 0 and 1, we found a good correspondence when matching them against brain domains that are known to activate during working memory and default mode function, respectively.Both community 0 (working memory) and community 1 (default mode) have components that together cover the regions known to activate in their associated systems.At the same time they also have a number of core components that are not associated with these, especially ones located in the occipital lobe.We perform the matching using the Neurosynth software (Yarkoni et al. 2011;Yarkoni 2011).In Additional file 4: Figure S4, we visualize the core components that best match active regions in their associated activation maps, next to the relevant z-axis slices of those maps.

Conclusion & discussion
By modeling the dynamic dependencies between brain regions in fMRI data, using a temporal network, and inferring inter-temporal dependencies in that network using the NFC principle, we are able to detect two salient and distinctly characteristic communities.The communities are shared across subjects, and for most subjects the communities are strongly modulated by the experimental conditions of the task.The first community primarily activates during the 2-back condition -and importantly not the 0-back conditionwhile the second activates during the fixation periods.
We do not find any significant correlations between the subjects' task performance accuracies and the extent to which they activate certain communities in synchrony with specific tasks.This is something we might have expected to be the case as it is reasonable to hypothesize that when a task is being performed poorly, the associated community activates less in synchrony with that task.A possible explanation for this is that the variation in task performance accuracy is small (almost everyone gets at least 60% right), and that even performing above some of the lowest scores requires significant engagement.
The overall method we present entails two steps, first constructing a temporal network and second clustering it.Together this yields five important hyper-parameters which impact the results: In the first step we (1) fix the edge density of the temporal network, and increasing this yields denser networks which in turn encode more of the weak component correlations as links.(2) We choose a threshold for component activation.Lowering this allows links between nodes representing less activated components.(3) We choose the length of the window inside which we measure correlations between components over time.Shortening this window increases the time-resolution but destabilizes the link dynamics of the temporal network.Note also, that we could have chosen a different measure of similarity rather than the Pearson correlation.In the second step, when performing community detection using NFC with Infomap, we must (4) choose a relax rate which encodes how freely information can flow between layers of time.Lowering the relax rate yields communities that are more constrained in time.(5) Finally, we set a neighborhood similarity threshold, below which inter-temporal links connecting states of each node are removed.Choosing a high threshold restricts information to only flow between temporal states that are highly similar, which yields communities with very well-defined component profiles.
In this proof-of-concept work, we do not set model parameters randomly, but instead sample a large number of solutions for different parameters and choose the set that, based on our qualitative assessment, gives the best results.However, it is clear that future work should include a more rigorous hyper-parameter tuning protocol.There are a two challenges in this regard, the most severe being to choose a cost function since community detection is an unsupervised method and thus difficult to quantitatively evaluate.Second, each evaluation step is slow (minutes) and with five hyper-parameters the search space is large, thus it may be necessary to employ specialized gradient-free techniques for sampling parameters.Our approach, with the chosen parameters, only recovers two communities that are repeatedly detected across different runs and population sizes, i.e. 10 and 100 subjects.It is reasonable to suppose that a better approach could render the modular structure of the temporal network in finer detail.In the current paper, it was, however, not our intention to give the best solution, but merely to demonstrate that a solution with this approach is feasible.
Finally, although computationally impractical, we would ideally want to perform the current analysis entirely without the initial dimensionality reduction step, as it puts a hard bound on the granularity of our results.Future work building on this paper should, therefore, consider using more fine-grained techniques for segmenting the brain into parcels or components.

Dataset
We obtain the dataset from the Human Connectome Project (HCP) (Barch et al. 2013).It originally comprises fMRI measurements for 1200 subjects with various tasks and restingstate data that all are repeated twice (in two separate phase encodings).In this paper, we only use the second run of the first 100 subjects.We analyze the n-back task designed originally to study working memory.Subjects are presented with images and are asked to press a button if it matches the image presented two instances back (2-back condition).Furthermore, they are asked in another experiment block to press when a target image is presented (0-back).Finally, there are periods of rest with no explicit task (fixation).The duration of an experiment is just under 5 min.For each subject, measurements are taken in 405 ticks with frequency 0.72 s.

ICA
We apply a group independent component analysis (ICA) (Calhoun et al. 2001;Varoquaux et al. 2010) to extract 50 spatially independent components using the nilearn Python library (Abraham et al. 2014).In group (spatial) ICA, subjects data are concatenated in time such that we end up with a matrix of size #voxels×#timepoints•#subjects•#runs.This is then first reduced using PCA to 50 components, and the spatial maps rotated according to some independence criterion to make them as spatially independent as possible.In our case, we used the cubic-approximation to the negative entropy, which is the default in nilearn.The ICA algorithm is run on the first 100 subjects in the data repository using both runs of the experiment.Preprocessing includes spatial smoothing (FWHM 6 mm) and high-pass filtering at 0.008 Hz.The ICA algorithm (fastICA (Hyvarinen 1999)) is restarted 10 times and the best solution is chosen using the cost-function as a measureof-fit.We finally rank the components by their temporal consistency over subjects.This is done by calculating for each component separately the temporal correlation between all pairs of subjects, then taking the average and ranking components in descending order.One component is removed due to a large overlap with regions associated with cerebrospinal fluid.
Component activations are first standardized such that the distribution of activities in each component over time was centered at 0 with standard deviation 1, then rescaled to [ 0, 1].

Temporal networks
We define a link between two components to be the correlation between their activity within a window of time, w.While comparing at each instant of time is meaningful for creating temporal networks from time-series in other datasets, fMRI measurements are too noisy for this.As discussed in the "Results" section, we choose the window-duration w = 25, and the window then slides one tick at a time, thus for T = 405 timesteps the resulting number of time-windows after this procedure is 381 = 405+1−25.When highpass filtering the fMRI signal at 0.008 Hz this time-window is short enough to induce spurious correlations between components (Leonardi and Van De Ville 2015;Zalesky and Breakspear 2015), and we, therefore, threshold away components that are not activated slightly above the mean (specifically 0.6).Moreover, for each subject we choose the top 1% links with highest positive correlation coefficient out of all possible links in the network, such that link density is consistent across subjects.Since we concatenate subjects in time and perform community detection on this new larger network, a constant density is necessary as subject variation would otherwise be interpreted by the model as temporal variation.

Community detection
We use Infomap with the Neighborhood Flow Coupling (NFC) method (Aslak et al. 2018) for community detection in our temporal networks.NFC infers dependencies across time between states of each node based on their neighborhood similarity and embeds inter-temporal links in the network proportional to these.Infomap minimizes the per-step description length of a random walker on the network, lower-bound by the Map Equation, by finding modules and submodules that capture large amounts of flow (Rosvall et al. 2009).This method takes two important hyperparameters: First, one has to choose the relax rate (we use 0.5), which describes how frequently the walker can transition between layer.A high relax rate leads to large amounts of flow across time, which gives rise to large communities that span many time-layers.With this dataset, we find that all nodes are classified as being in the same community when the relax rate is too high.Second, one has to choose a similarity threshold for connecting sibling-nodes across time, above which interlayer links may be created (we use 0.75).A low threshold allows two states of a node, residing in different temporal layers, to have an interlayer link that connects them even when they only have weakly similar neighborhoods.Reversely, a high threshold sets a high requirement for neighborhood similarity for interlayer links to be created.In this data, a too high threshold yields communities that are only detected when the few nodes that make them up are linked.These communities in turn share many nodes but have different labels.Low thresholds create the opposite problem, which is that only few communities with very broad component profiles are detected.
In this work, we do not search for parameters in any rigorous way, but simply sample a large number of solutions and choose a corresponding parameter set that gives meaningful results.There are two reasons for this.First, we cannot easily define a good objective function for an optimization program, since community detection is an unsupervised problem.Second, the primary purpose of this study is to demonstrate the feasibility of using the described temporal community detection framework on fMRI data.To this end, we do not strictly need the optimal parameter set, but only a set of parameters that serves our point.

Fig. 2
Fig. 2 Temporal communities of 10 "best" performers concatenated.a Two largest communities in the subject-concatenated network, reshaped with one subject per row to enable comparison of temporal profiles.Each row pertains to an individual subject and visualizes the temporal dynamics of the community as a white strip of varying height.The background color illustrates the task order (see Fig. 1 legend).Notice that the task order deviates for subject 96.b Inter-subject correlation of temporal activity for the two largest communities.Mean correlations (not including diagonal of ones) are 0.33 and 0.41 respectively.c Component distributions.See Additional file 1: Figure S1 for component brain maps