Both activated and less‐activated regions identified by functional MRI reconfigure to support task executions

Abstract Introduction Functional magnetic resonance imaging (fMRI) has become very important for noninvasively characterizing BOLD signal fluctuations, which reflect the changes in neuronal firings in the brain. Unlike the activation detection strategy utilized with fMRI, which only emphasizes the synchronicity between the functional nodes (activated regions) and the task design, brain connectivity and network theory are able to decipher the interactive structure across the entire brain. However, little is known about whether and how the activated/less‐activated interactions are associated with the functional changes that occur when the brain changes from the resting state to a task state. What are the key networks that play important roles in the brain state changes? Methods We used the fMRI data from the Human Connectome Project S500 release to examine the changes of network efficiency, interaction strength, and fractional modularity contributions of both the local and global networks, when the subjects change from the resting state to seven different task states. Results We found that, from the resting state to each of seven task states, both the activated and less‐activated regions had significantly changed to be in line with, and comparably contributed to, a global network reconfiguration. We also found that three networks, the default mode network, frontoparietal network, and salience network, dominated the flexible reconfiguration of the brain. Conclusions This study shows quantitatively that contributions from both activated and less‐activated regions enable the global functional network to respond when the brain switches from the resting state to a task state and suggests the necessity of considering large‐scale networks (rather than only activated regions) when investigating brain functions in imaging cognitive neuroscience.

The detection of the activation regions for the 7 tasks was implemented by FSL (version 5.0.9)/FEAT (version 6.00) (Jenkinson et al., 2012) using parameter settings that were extremely similar to those in Barch et al.'s report (2013), the key steps of which are described here. Both the fMRI time series and the GLM were temporally filtered with a Gaussian-weighted linear highpass filter with a cutoff of 200 s. The volume data was spatially smoothed using a 3D Gaussian kernel with FWHM = 4 mm. Since the preprocessed data from the HCP had already been skull-stripped and realigned to MNI 2 mm space, these operations, including motion correction and brain tissue extraction in the FSL/FEAT GUI, were not repeated. The FILM prewhitening step was chosen to increase the validity and efficiency of the statistics (Woolrich et al., 2001). Three items, the mean white matter signal, the cerebrospinal fluid signal, and the 12-item head movement parameters (from the S500 dataset), were set as confounding variables.
The subsequent statistical analysis was conducted in two steps. At the subject level, a within-subjects fixed-effects analysis in FSL/FEAT was used to estimate the average effects across the runs, with a cluster based threshold Z = 1.96 as the activation threshold and P = 0.05 as the Monte Carlo-based cluster-level correction. Then at the group level, a mixed-effects analysis implemented in FSL/FLAME (FMRIB's Local Analysis of Mixed Effects) (Beckmann et al., 2003) was used to estimate the average effects of interest separately for the 7 task groups, with a cluster-based threshold Z = 2.32 and P = 0.05. In this study only the positive activation regions were used since, to date, the physiological interpretations of negative activations obtained from BOLD fMRI remain controversial (Bianciardi et al., 2011;Hu and Huang, 2015;Shih et al., 2009).

Section S2: Act and nonAct definitions
As described in the Methods section in the manuscript, for the strategy discriminating Act/nonAct regions, we seek for a balance between the specificity and generalizations of functional activations. This means when try to constrain the regions for a specific task, we should also guarantee that the size of regions are not too small and able to cover a sufficient number of ROIs from the Power264 ROI list (Power et al., 2011) for network modularity evaluation.
The task contrasts used for each task were to detect task-specific functional regions and remove any confounding activations. For example, in the WM task, the subjects were cued by pictures on the screen to make inferences and decisions. The activation maps of [2back -0back] and [2backbaseline] were significantly different. The choice of [2back -baseline] resulted in many areas of the primary cortex being activated, such as the visual cortex (see Fig. S5 in the Supplemental Materials). In this study, however, we were interested in higher cognitive functions, including planning and problem-solving (Cohen et al., 1997), so we adopted the contrast of [2back -0back] instead of [2back -baseline]. In addition, we tried to avoid generating activation maps which were too small, because too few ROIs would lead to unstable network construction and analysis. For example, for the gambling task or the motor task, we chose alternative contrasts rather than contrasts that only activate a very few specific areas. The details of the contrast used for each task are listed in Table 1. The activation masks for the three separate datasets (LR, RL and the averaged) are available on our website (https://github.com/nmzuo/Act-L-Act-network) and the masks largely coincide across the three states.
Here in Table S1 we list all the contrasts for each task state and the descriptions explaining the choice in the manuscript. About the principles how these contrast were designed and potential activated regions please refer to the Human Connectome Project (HCP) website, http://www.humanconnectome.org/documentation/Q1/task-fMRI-protocol-details.html and http://www.humanconnectome.org/documentation/phase1-paradigms-fMRI.html.
All the activation maps can be downloaded (see Appendix B). Table S2 listed the contrast used in the manuscript. Table S1. The contrasts for each task used in the manuscript are based on a series of pre-studies. The name "cope" is in line with the routine of the FSL/FEAT, which means contrast of parameter estimates (Jenkinson et al., 2012)  Use cope1. The activation regions from the cope2 are too small (shown in Fig. S1).
WM cope1: 2bk-0bk; cope2: 2bk-baseline; (combine (body, face, place, tool) as a whole) Use cope1. By contrast to the cope2, the results from cope1 denote the core regions which would be called by more difficult cognitive tasks. Actually, Contrasting Fig. S1 to Fig. 1, they share most areas since they perform very similar tasks except for the difficulty level.   Barch's paper (2013), the slices were chosen so that they would be in similar positions, and the color window is also [2.32 5]. (Throughout the manuscript, the results are presented in the same sequence in the related figures.) In Fig. S2, we show the activation map of those contrasts not included in the manuscript (Table  S2). Only 4 tasks are listed here because for the other 3 tasks either it has a single 'cope', or the combination of all "copes" has been used.

Section S3: The validation results for different thresholds and different datasets
We have 463 subjects (194 male, age 29.1 ± 3.5 years) and 465 subjects (191 male, age 29.1 ± 3.5 years) in the LR and RL dataset, respectively. For simplification, in this section, A, B and C separately indicate the thresholds of 5%, 10% and 15%. So, the label 'LR-A' means the result is based on the dataset LR with the threshold of 5%, and so forth. The results 'AVG-C' is presented in the manuscript as the our main findings.

/ 25
The global efficiency comparisons between different mental states (resting and 7 tasks) (ref. Fig. 2

in the manuscript)
The efficiency indices of the whole brain in different mental states, including the resting state and the 7 tasks (gambling, motor, social cognition, emotion processing, language, and relational processing) are shown in Fig. S3. The results from the three datasets and three thresholds generated very similar trend that the global efficiency in the task states are significantly higher than the one in the resting state. The interactions between Act/nonAct regions for the 7 resting-task pairs (ref. Fig. 4

)
Comparisons of the strength of the interaction between the Act and nonAct regions in the resting and task states are shown in Fig. S4. Within each column pair, the left shows the interaction between the two classes of regions in the resting state and the right shows that in the task state. The overall results are well in line with the main findings in the manuscript except that, for the social cognition network with sparse connectivity (5% and 10% at RL dataset), the interaction between Act/nonAct regions did not decrease. However, with the increasing the density, the trend also coincides well with the main findings that the Act/nonAct interactions are increasing for the first 6 tasks but the opposite trend for the WM task. Correlations between the changes in the interactions between the two regions and the changes in the global efficiency are listed in Table S3. The individual correlation strength R and significance level P-value are shown. The "thr" indicates the network density threshold.
These results collectively demonstrate that, for the first 6 resting-task comparison, the less interaction the more efficiency; and for the resting-WM comparison, the more interaction the more efficiency. The efficiency comparisons for the resting and the task (separately for the Act and nonAct) (ref. Fig. 6) Comparisons of efficiency between the task states and resting state for Act (left panel)/nonAct (right panel) regions. The 7 column pairs in each panel are corresponding to the 7 tasks (Fig. S5).
The left bar in of each column pair shows the task state and the right shows the resting state. The similar results for the Act regions do not show consistent increase trends although they have 13 / 25 statistically significant difference between the resting and task states.
These results collectively demonstrate that, for all the resting-task comparisons, the efficiency in Act regions significantly changes but do not show consistent trend. On the other hand, the nonAct regions shows consistently increasing efficiency. Correlations between the changes in the efficiency indices between the two regions and the changes in the global efficiency are listed in Table S4. The individual correlation strength R and significance level P-value are presented. The panels A and B indicate the Act and nonAct regions, respectively.
These results consistently show that the efficiencies in Act/nonAct regions are increasing coherently with the global network efficiency.

Section S4: The fractional modularity comparisons based on different datasets and density thresholds
When we computed the fractional modularity Q f based on the Eq. S2 in the manuscript, the Qf was normalized by the size, namely the number of nodes for each class of regions, for a fair comparison (see Table 1). Another way for the weighting is to use all possible connections based on the nodes and this is the combination of the node, namely choosing 2 from number of nodes. After we repeated the processing in Table 1 using this weight we reproduced our findings (based on the three datasets and three thresholds, see Table S5) that from the perspective of network modularity measure, the contributions from fractional modularity of the nonAct regions is comparable to the ones from the Act regions. Section S5: The parameter settings for the generalized

Louvain-based modularity computing
The Fig. 3 in the manuscript is for an overall illustration to demonstrate that the entire brain network, including both the Active and nonActive parts, reorganized the network architectures during the transition from the baseline (resting state) to each task (7 different tasks in total). In each panel in Fig. 3, it demonstrates the connectivity flow from the resting state to the task state. The modular partition is based on the generalized Louvain method (Mucha et al., 2010) by repetitions of 100 times to induce the random bias. The resolution parameter γ (In Eq. S2) is related to the final partition scheme through maximizing the modularity, and it has been demonstrated that for the weighted and small network the resolution influence tend to be much less significant (Bassett et al., 2011;Good et al., 2010). In our study a default resolution parameter ( γ = 1 ) was used in the Mucha's implementation (Mucha et al., 2010) (http://netwiki.amath.unc.edu/GenLouvain/GenLouvain). See Fig. S6 for the detailed comparisons and the rationale of the final choice. The number of modules can be affected by the resolution parameter γ and the network density. In the revised manuscript, we use the binarized association matrix and our results demonstrated comparable numbers of modules as reported in (Power et al., 2011;Yeo et al., 2011). Our current modular partition derived from the weighted association matrix usually have only 2-4 main results, which largely coincides with the existing reports in Bassett's work (where the Fig. 2 used the similar methods as ours) (Bassett et al., 2011) and in Vatansever's work ( Fig. 3 using Infomap algorithm (Rosvall and Bergstrom, 2008)) (Vatansever et al., 2015). And this result is also well in line with Power's report (right panel in Fig.  S6 using Newman's algorithm (Newman and Girvan, 2004)).
We further investigated the influences of on the resolution parameter γ and the number of modular partitions on the modularity computation in our manuscript. It is well known that the modularity optimization based methods are prone to be limited by the modular resolution (Fortunato and Barthelemy, 2007). In our study, to make the results comparable, a unified resolution parameter is used, γ = 1, which is also the default parameter for generalized Louvain method (http://netwiki.amath.unc.edu/GenLouvain/GenLouvain). A test was performed for the influence of γ on the resting state network for subject 100307. The left panel of Fig. S6 shows the modularity indices when γ changes from 0.2 to 2.0 by step 0.2. It shows that γ=0.8 or 1.0 could achieve the maximized modularity, providing the rationale for γ= 1.0 used in our study. The right panel of Fig. S6 shows the similarity between the partition results using different γ values. The similarity was measured by the z-score of the Rand similarity coefficient (Traud et al., 2011). Section S6: The illustrations for the connectivity transitions between the resting state to the task states The Fig. 3 in the manuscript showed the connectivity flow when the brain changes the state from the resting to the tasks, where the nodes were grouped as the modular assignments. Here we have generated a series of connectivity matrix maps to describe the connectivity flow between the Act and nonAct communities (see Fig. S7 below). The connectivity matrices were based on the averaged database (described above) and each matrix is the averaged connectivity for each mental state (the resting and the 7 tasks) across the 453 objects. The order for the nodes is based on the Power et al's paper (2011), but we grouped the nodes as Act and nonAct communities for each resting-task pair (labeled as the left vertical texts in each panel). Judging from the visual inspection, in line with the Fig. 3 in the manuscript, they consistently showed global connectivity transitions.

Section S7: Methodological considerations
Definitions of the Act and nonAct regions. In this study, the investigation into the interactions between the Act and nonAct regions was based on an activation detection strategy, which depends on several factors, including the regression methods, contrast design, and the threshold at the participant and group levels. We utilized a traditional GLM model implemented in FSL (Smith et al., 2004), which can be used to generate a generally accepted activation map for the 7 task states ( Barch et al., 2013). Two additional factors, setting an appropriate task contrast (see Table 1) and reserving the regions with a considerable size to avoid too few ROIs remaining in the activation regions, were utilized to seek a balance between sensitivity and specification. To this end, the task contrasts were set as in Table 1 with Z = 1.96 and P = 0.05 for the individual participants and Z = 2.32 and P = 0.05 at the group level in the FSL/FEAT pipeline. In Section 1 of the Supplementary Material, for a better comparison, we listed the activation detection results from different task-baseline contrasts. It should be noted that both different task contrasts and different parameter settings in activation detection could result in different activation maps and different patterns of functional connectivity for the Active and nonAct regions. The activation detection strategy in our study is one of the mostly used ones, and the activation regions are well accepted in the literature for the specific tasks ( Barch et al., 2013). Therefore, it is unnecessary, if not impossible, to enumerate all the possible combinations of contrasts and parameter settings for robust determination of the Act and nonAct regions. From this view, our study could shed light on the collaborations between the Active and nonAct regions. The final activation masks adopted in this work are also available for download.

Construction of functional brain networks.
Since different strategies to measure the functional connectivity and then construct a network exist, in this study, for a complete validation of our findings, three connectivity density thresholds, 5%, 10% and 15%, have been considered (Bullmore and Sporns, 2009;Rubinov and Sporns, 2010;van Wijk et al., 2010) to construct the functional network separately in three datasets, LR phase encoding direction, RL phase encoding direction and the averaged connectivity after the Fisher z-transform. Moreover, before calculating the functional connectivity for the task fMRI time series, the mean task activity has been regressed out and the residuals were used in constructing networks (Cole et al., 2013).
Modularity and efficiency computing. Modularity and efficiency were the two main measures in this study that were used to characterize the global and local (within one class of regions) network evolvement. An extended optimization-based Louvain method was used in this study (Mucha et al., 2010). In addition, a fractional modularity measure was adopted to compare the contributions from the Act and the nonAct regions to the global changes in the modularity separately. The modularity index focuses on the functional segregation property of a network, and the efficiency index was used to characterize the information communicating efficiency of two regions (Latora and Marchiori, 2003;Rubinov and Sporns, 2010). Therefore, the regions assignment and the regions between the Act and nonAct regions of the network have been considerably characterized in the current study.