Frequency-specific task modulation of human brain functional networks: A fast fMRI study

How coherent neural oscillations are involved in task execution is a fundamental question in neuroscience. Although several electrophysiological studies have tackled this issue, the brain-wide task modulation of neural coherence remains uncharacterized. Here, with a fast fMRI technique, we studied shifts of brain-wide neural coherence across different task states in the ultraslow frequency range (0.01-0.7Hz). First, we examined whether the shifts of the brain-wide neural coherence occur in a frequency-dependent manner. We quantified the shift of a region's average neural coherence by the inter-state variance of the mean coherence between the region and the rest of the brain. A clustering analysis based on the variance's spatial correlation between frequency components revealed four frequency bands (0.01-0.15Hz, 0.15-0.37Hz, 0.37-0.53Hz, and 0.53-0.7Hz) showing band-specific shifts of the brain-wide neural coherence. Next, we investigated the similarity of the inter-state variance's spectra between all pairs of regions. We found that regions showing similar spectra correspond to those forming functional modules of the brain network. Then, we investigated the relationship between identified frequency bands and modules' inter-state variances. We found that modules showing the highest variance are those made up of parieto-occipital regions at 0.01-0.15Hz, while it is replaced with another consisting of frontal regions above 0.15Hz. Furthermore, these modules showed specific shifting patterns of the mean coherence across states at 0.01-0.15Hz and above 0.15Hz, suggesting that identified frequency bands differentially contribute to neural interactions during task execution. Our results highlight that usage of the fast fMRI enables brain-wide investigation of neural coherence up to 0.7Hz, which opens a promising track for assessment of the large-scale neural interactions in the ultraslow frequency range.


Introduction
A prominent attribute of brain dynamics is coherent oscillatory activity across separate areas. The frequency of the coherent oscillation is determined by the anatomical and biophysical properties of the underlying neural circuits. Such properties are hypothesized to specify circuits' functionality in the network ( Siegel et al., 2008 ). Recent studies have found that the task modulation of neural interactions occurs in a frequency-dependent manner ( Bastos et al., 2015;Gregoriou et al., 2009;Haegens et al., 2011 ;Hipp et al., 2011 ;Pesaran et al., 2008 ;Siegel et al., 2008 ), supporting the idea that the frequency reflects the functionality. Therefore, examining large-scale neural interactions based on frequency-specific signals would promote better understanding about the functional network organization of the brain than has the corresponding examination based on the broadband signals. However, because most of the studies have focused on limited regions of interest, the way the whole-brain frequency-specific interactions change with ing state fMRI measurements have better characterized the whole-brain functional network and successfully revealed its attributes -such as the existence of subnetworks consisting of functionally relevant regions with dense interconnections, which are thought to work as functional modules ( Bertolero et al., 2015 ;Crossley et al., 2013 ;Smith et al., 2009 ). Furthermore, a recent study demonstrated that the multiband fMRI technique can detect relatively faster neural oscillations (~0.75 Hz) than previously thought ( Lewis et al., 2016 ). Moreover, neural interactions observed through fMRI have the frequency-specificity in their strength and topology when at rest ( Chen and Glover, 2015 ;De Domenico et al., 2016 ;Liao et al., 2013 ;Sasai et al., 2011Sasai et al., , 2014Thompson and Fransson, 2015 ). These pieces of evidence lead to the idea that fMRI potentially enables the whole brain investigation of neural interactions' frequency-specific task modulation in the ultraslow frequency range that may subserve cognitive functions ( He and Raichle, 2009 ). However, the functional modulation has been studied so far only with broadband fMRI signals.
To tackle this issue, we investigated whether the modulation of fMRIbased whole-brain neural interactions occurs in the frequency dependent manner. This was accomplished by comparing the relative regional variance of frequency-specific neural interactions across multiple cognitive states among different frequency components. The frequencyspecific neural interactions were defined based on the coherence of signals obtained by the multiband fMRI, which we call the frequencyspecific functional connectivity, among 264 regions of interest covering whole brain ( Power et al., 2011 ). To enhance generalizability, we selected five different states including one resting and four task states that require subjects to maintain different foci: (1) detection of changes in luminance of a fixation point on the display (visual attention state), (2) detection of changes in the tone of a sound stimuli (auditory attention state), (3) heartbeat counting (interoceptive attention state), and (4) self-referential memory recollection (self-referential state). Furthermore, we quantified the similarity of the frequency-dependency of the inter-state variance between all pairs of regions, and we examined whether regions in the same module have spectral attributes that are similar to the neural interactions.

Participants
A total of 37 healthy adults (15 men and 13 women; age range, 22-44 years) participated in this study. All participants were right-handed according to the Edinburgh handedness inventory ( Oldfield, 1971 ). No participant had a history of neurological or psychiatric illness. The protocol of this experiment was approved by the institutional review board of the National Institute for Physiological Sciences, Okazaki, Japan, Japan, and the experiments were undertaken in compliance with national legislation and the Code of Ethical Principles for Medical Research Involving Human Subjects of the World Medical Association (the Declaration of Helsinki). To participate in the study, all participants gave their written informed consent.

Experimental design
All participants completed fMRI at rest and during 4 cognitive tasks: [1] a visual attention task, in which participants were instructed to keep looking at a fixation point in the middle of a black screen to capture any change of its luminance, [2] an auditory attention task, in which participants were instructed to keep hearing a pure tone to try to capture any change of the scale, [3] a heartbeat counting task, in which participants were instructed to count all of their heartbeats during the scanning, and [4] a self-referential memory recollection task, in which participants were instructed to reflect on their personality. For visual and auditory tasks, we only used stationary stimuli, which are unchanged over time. This allowed us to estimate the functional connectivity that cannot be explained by task designs. For example, we kept the fixation point's luminance identical over the visual task period. We continuously presented the same pure tone across the scan, with no breaks and no change in frequency, during the auditory task period (see Fig. S1). Although we used these stationary stimuli for both conditions, we instructed participants to try to detect changes in stimuli, as if the stimuli would change during the tasks. This is because we wanted to help participants' attention to the stimuli, while we needed to measure brain activities unaffected by the change of the external inputs. The length of the recording was 5.5 min for all sessions. After each fMRI scanning session, participants were asked to report the degree of drowsiness during the task, the number of stimuli changes they captured (only after [1] and [2]), the number of heartbeats they counted (only after [3]), and what they reflected on (only after [4]). Before any task fMRI session, all participants completed resting state fMRI recordings.

Stimuli presentation
All stimulus presentation was conducted using presentation software (Neurobehavioral Systems, Inc., USA) implemented on a personal computer (dc7900; Hewlett-Packard, Ltd., USA). The resolution of visual stimuli was 1024 × 768, and the refresh rate was 60 Hz. The visual stimulus was projected onto a translucent screen by an LCD projector (CP-SX12000; Hitachi, Ltd., Japan) located outside the shield room. Participants watched the projected visual stimulus via a mirror attached to the 32-ch head coil. The distance between the screen and each participants' eyes was approximately 175 cm; the visual angle of stimuli was 13.8 (horizontal) x 10.4 (vertical). Audio stimuli, generated by the presentation software, were presented to participants through MRI-contingent headphones system (Kiyohara Optics, Inc., JAPAN).

Polysomnographic recording
We carried out the polysomnographic recording, i.e. EEG, vertical and horizontal electrooculography (EOG), electrocardiography (ECG), and respiration, simultaneously with all fMRI recordings, to obtain the objective measure of arousal level based on the standard criteria ( Rechtschaffen and Kales, 1968 ) and to conduct retrospective removal of physiological confounds in the fMRI signals ( Anderson et al., 2011 ).
Twenty-seven channels of EEG, vertical and horizontal EOG, and respiration were acquired using MRI-contingent amplifiers (BrainAmp MR and BrainAmp ExG MR, Brain Products GmbH, Germany) and an electrode cap with Ag/AgCl ring electrodes (BrainProducts GmbH, Germany). All EEG electrodes were placed according to the international 10-20 system. Horizontal EOG was recorded from the left eye with an Ag/AgCl ring electrode placed at the external canthi, and vertical EOG was recorded from the left eye with an Ag/AgCl ring electrode placed at a suborbital area. A reference electrode was placed at the middle point between Fz and Cz. Whole raw data were recorded at 5 kHz with a bandpass filter (0.016 and 250 Hz), using Brain Vision Recorder software (Brain Products GmbH, Germany). In addition to polysomnographic recording, participants' faces were monitored using an MRI-contingent infrared video camera to guarantee participants did not fall asleep.

EEG processing
We used two MATLAB toolboxes for data preprocessing. First, gradient-switching and pulse artifacts induced by simultaneous MR imaging were removed by using FASST toolbox ( Leclercq et al., 2011 ; http://www.montefiore.ulg.ac.be/~phillips/FASST.html/ ). Then the following processes were conducted by using EEGLAB toolbox ( http://sccn.ucsd.edu/eeglab/ ): (1) a band-pass filter was used to keep components between 1 and 30 Hz, (2) data periods showing large amplitude over 6 standard deviations from the mean were excluded, and (3) independent components reflecting eye movement, remaining MR artifacts, and anything generated by single channel were removed. To make sure that subjects were awake during the experiments, manual sleep scoring was performed over 30-s epochs for each state based on standard criteria ( Berry et al., 2015 ).

fMRI preprocessing
Functional MR volumes were motion-corrected using the SPM8 package (Wellcome Department of Imaging Neuroscience, London, UK). After spatially smoothing the volumes with a 5-mm full-width-at-halfmaximum Gaussian blur, we registered them to the standard MNI space. fMRI data are generally contaminated with physiological noises, such as subject motion, respiration, and cardiac activities, resulting in coherent signal fluctuations across the brain (e.g., global signals). To minimize these effects on the functional connectivity estimation, we modified the phase-shifted soft-tissue correction method (PSTCor; Anderson et al., 2011 ). PSTCor is a general linear model-based denoising method that subtracts nuisance regressors after optimally adjusting their time-lags from the global signal. A limitation of PSTCor is that it does not use some regressors significantly correlated with fMRI signals, such as heart rate variability or heart rate time-course convoluted with cardiac response function ( Chang et al., 2009 ). Furthermore, it does not help addressing large amplitude spikes caused by head movements ( Power et al., 2012 ). Due to these limitations, we modified by using new regressors on top of the original set of regressors introduced by Anderson et al. The added regressors include: 1. RETROICOR regressors ( Glover et al., 2000 ), 2. white matter time series, 3. cerebrospinal fluid time series, 4. soft tissue time series, 5. respiration volume per time, 6. respiration volume per time convolved with respiration response function, 7. heart rate per time, 8. heart rate variability per time, 9. heart rate per time convolved with cardiac response function, 10. time series of motion parameters from realignment procedure, 11. temporal derivatives of 1-10, and 12. spike regressors due to Satterthwaite et al. (2013) . We followed the method introduced by Chang et al. (2009) to obtain regressors 5-9. We selected this set of regressors after validating their effectiveness in reducing confounds (see Figures S1-S3). The motion spike regressors were made by identifying significant spike events present in each fMRI time-series. Spikes were identifying by calculating a measurement called framewise displacement (FD), which is the sum of the absolute values of the differentiated realignment estimates (by backward differences) at every time point ( Power et al., 2012 ). We made a regressor for every spike, so the number of regressors equals the number of identified spikes ( Satterthwaite et al., 2013 ). It has been shown that denoising through the general linear model could artifactually transform structured functional patterns of low-frequency signals to the high-frequency band ( Chen et al., 2017 ). This effect can be mitigated by band-limiting the low-frequency components. Thus, we denoised using signals and regressors filtered before the regression ( Hallquist et al., 2013 ). Specifically, we decomposed the full frequency range into the low-frequency band (0.01-0.15 Hz) and the high-frequency band (0.15-0.7 Hz) based on the power spectrum. After processing signals and regressors with bandpass filters to remove components outside of the frequency band of interest, we denoised using the general linear model at each frequency band separately.

ROI
A set of 264 regions of interest (ROI) introduced by Power et al. (2011) was used to define nodes of the functional connectivity network. Signals in voxels within each ROI (5 mm radii) were averaged to obtain representative signals for the ROI.

Frequency-specific functional connectivity
We estimated functional connectivity by the magnitude-squared coherence for all pairs of the ROI. The coherence measures the linear and time-invariant relationship between two signals x and y at frequency and is defined as follows: where C xy ( ) refers to the coherence, P xy ( ) is the cross-spectrum of x and y, P xx ( ) is the power spectrum of x , and P yy ( ) is the power spectrum of y . For each pair of signals, we obtained frequency-specific functional connectivity matrices by averaging the coherence values within 68 narrow, 50%-overlapping frequency bands with bandwidths of 0.02 Hz ( Fig. 1 A). We have used the same definition of functional connectivity in our previous studies ( De Domenico et al., 2016 ;Sasai et al., 2011Sasai et al., , 2014.

Strength and diversity
The strength is defined as the sum of connectivity strengths of a node in the functional connectivity network, quantifying the degree to which the corresponding region coordinates with all other regions in the brain. We define the diversity as the variance of the strength across all states. As such, the diversity measures how much the region differentiates its interactions with the rest of the brain across states on average. We obtained strengths of all nodes at each frequency band as a row vector containing the sum of each column of the frequency-specific functional connectivity matrix, resulting in five 68-by-264 matrices ( Fig. 1 B, strength matrices). Then we calculated the variance of strengths across five states for each node at each frequency band and made the variance matrix ( Fig. 1 C). Finally we standardized each element of the variance matrix so that each row vector has a mean of 0 and a standard deviation of 1. We used this standardized variance to measure the relative diversity of the nodal strength at each frequency. We refer to the standardized variance matrix as the "diversity matrix, " its row vectors as "diversity profiles, " and its column vectors as "diversity spectra " ( Fig. 1 D).

Frequency cluster and brain module detection
We used a network community detection method called the "Louvain method " ( Blondel et al., 2008 ) to find both frequency clusters and brain modules. The Louvain method is one of the greedy methods. It searches for clusters of nodes in a network by optimizing a cost function called Modularity, the density of edges inside clusters to edges outside clusters. We used a MATLAB function, "community_louvain.m, " implemented in the Brain Connectivity Toolbox ( Rubinov and Sporns, 2010 ), to search the optimal cluster partitions based on similarity matrices (see the next section). We applied this method to the similarity matrices we chose to group, ones showing similarity between all pairs of units (frequency or ROI). Two key parameters affect grouping results, (1) network resolution parameter, and (2) threshold applied to the similarity matrix. We optimized these parameters to maximize the stability of our results as follows. First, we used the group-averaged similarity matrix instead of individual similarity matrices separately. Specifically, in the brain module detection analysis, we obtained similarity matrices of diversity spectra for all subjects and averaged them to get the group-mean similarity matrix. We also used the group-average similarity matrix in the frequency cluster detection analysis. Then, we searched for optimal parameters to provide the most stable partitions. The optimal network resolution parameter was chosen by using the method introduced by Traag et al. (2013) . In the paper, they proposed the method to quantify the significance and stability of community structure for a network resolution parameter. So, we used the product of the significance, which is minus the logarithm of the probability that identified clusters appear in a random graph of the same size, and the stability, which is the similarity among 25 trials of community detection, as an evaluation function to find the optimal resolution. As the stability, we used a recently proposed, improved measure for the partition similarity called adjusted mutual information instead of Traag et al. 's original measure ( Jeub et al., 2018 ). We changed the network resolution parameter from 0.1 to 10 in increments of 0.1 and found the optimal resolution parameter maximizing the above measure. For the threshold applied to the similarity matrix, we used the sparsity, the ratio of the number of edges over the possible maximum number of edges. We explored optimal sparsity by evaluating the stability of the detected partitions by resampling the similarity matrix through bootstrapping subjects. We quantified the similarity of detected communities across the resamples using normalized mutual information of the partition. We looked for optimal sparsity by changing it from 0.05 to 1 in increments of 0.05. Although we identified 1 as the optimal sparsity for brain module detection, multiple sparsity values achieved the highest stability for the frequency cluster detection. Thus, we further carried out the clustering using the summary matrix of frequency cluster assignments across all sparsity values achieving the highest stability. Specifically, we first made a binary matrix whose entry shows that the corresponding pair of frequencies is in the same cluster (1) or not (0). We made this matrix for the cluster assignment obtained by using each threshold achieving the highest stability. Then, we averaged these binary matrices together to get one summary matrix. Using this matrix as a new similarity matrix, we conducted frequency detection with the Louvain method.

EEG-fMRI coherence analysis
We conducted a coherence analysis between EEG and fMRI signals. To do so, we first obtained an EEG power time-course using timefrequency analysis. Specifically, we estimated the EEG power spectrum for each time-frame of fMRI acquisition (700 ms), and thus got a timecourse of the power for each EEG frequency component of an EEG channel. We repeated power estimation for all twenty-seven EEG channels. The EEG power time-course was used for two analyses in this study: (1) clustering analysis and (2) examination of the frequency-wise relationship between fMRI and EEG. For the first analysis, we obtained a coherence spectrum between an ROI's fMRI time-course and an EEG channel's power time-course of an EEG frequency component. This estimation of the coherence spectrum was repeated for all combinations of fMRI ROI, EEG channel, and EEG frequency component, which makes a 3-dimensional matrix (dimension 1: fMRI frequency, dimension 2: EEG channel, and dimension 3: EEG frequency) for each fMRI ROI. Using the 3-dimensional correlation of these matrices, we obtained a similarity matrix, which we used for the clustering analysis (see the next section). For the second analysis, we first obtained the mean fMRI timecourses of each brain modules and those across all ROIs. Using the mean time-courses instead of ROI-wise time-courses, we repeated the estimation of the coherence spectrum mentioned above, which creates the 3-dimensional matrices for brain modules. Averaging the matrix over the EEG channels, we obtained a matrix showing coherence between each fMRI frequency component and each EEG frequency component. To summarize the relationship with respect to each fMRI frequency band identified in this study, we obtained the means of the matrix rows corresponding to the identified frequency band.
It is known that there are substantial time-lags between neural events captured by EEG and corresponding BOLD response. Usually, this issue is addressed using the convolution of hemodynamic response function (HRF) with EEG events. However, it remains unclear what type of HRF is best suited here because the convolution of the canonical HRF functions as a low-pass filter and significantly reduces the high-frequency components of fMRI signals. Furthermore, Lewis et al. (2016) have already raised the possibility that high-frequency stimuli would produce a hemodynamic response that reaches its peak in a much shorter duration (~2 s) in comparison to the canonical HRF. Thus, we repeated the coherence analysis without the HRF, while changing the time-lag between EEG and fMRI time-courses increments of 1 second, from 1 second to 5 s. Finally, we selected an optimal time-lag that maximizes the deviation of the global mean coherence value from its null distribution obtained through the phase shuffling of both EEG and fMRI time-courses.

Similarity matrices
The similarity matrix used to the frequency cluster detection was obtained by quantifying the similarity of diversity profiles between all pairs of frequency bands with Spearman's correlation. We only used this matrix for the frequency cluster detection. On the other hand, we used three different similarity matrices for the brain module detection, which are the similarity matrix of diversity spectra, functional connectivity matrix, and the similarity matrix of EEG-fMRI coherence (see the previous section). The matrix of diversity spectra was obtained with Spearman's correlation for all pairs of ROIs. The functional connectivity matrix was obtained, calculating Pearson's correlation coefficients of time-series for all pairs of ROIs. The matrix of EEG-fMRI coherence was made as described in the previous section.

Comparison of diversities among modules
We extracted the mean diversity spectrum for each detected module by averaging the diversity spectra of all assigned ROIs. These values were then compared among modules after being averaged in the following 4 frequency bands identified by the clustering analysis in the frequency band detection (see above): (1) 0.01-0.15 Hz, (2) 0.15-0.37 Hz, (3) 0.37-0.53 Hz, and (4) 0.53-0.70 Hz. We conducted the analysis of variance (ANOVA) to test if variances are the same across modules and frequency ranges. We performed post-hoc tests to test the null hypothesis that the average diversities for all pairs of modules have equal means in each frequency range. The modified sequentially rejective multiple test procedure was used to account for the multiple comparisons ( Shaffer, 1986 ).
We also tried to investigate if any of the identified modules corresponds to the flexible hubs reported by Cole et al. (2013) . In their paper, Cole et al. used Power et al. (2011) 's partition of brain modules and argued that the flexible hub is the frontoparietal network. Thus, we also identified the most diversified module, according to Power's definition of module parcellations. Specifically, we obtained diversity values for each ROI at each band, and we averaged these values for ROIs assigned to the same module. Then, we identified the most diversified module at each band based on the group-mean diversity value. We conducted a paired t -test to statistically examine whether the variance of the most diversified module is the same as those of other modules. We applied the Bonferroni correction to get corrected p -values for these multiple comparisons.

Comparison of inter-state variance in strengths among modules
To test if all modules differentiate connectivity strengths across states in the same way, we standardized the mean functional connectivity strength of a frequency range across states for each node and averaged them in each module. By regarding the standardized strength in each state as the coordinate of a dimension, we projected the strengths of all states of a module onto a point of a five-dimensional Euclidean space. We calculated Euclidean distances among all projected points to get quantities of dissimilarity of inter-state variance of strengths among all modules. The mean dissimilarity across all subjects was statistically compared to those obtained by permuting module labels to test the null hypothesis of no dissimilarity.

Task-specific functional connectivity detection
We examined whether there is any task-specific functional connectivity modulation at each frequency band. To do so, we obtained a functional connectivity matrix for each task at each frequency band by averaging coherence within the frequency band. We applied the network-based statistical (NBS) analysis ( Zalesky et al., 2010 ) to test whether there are any significant task-specific FC patterns at each identified frequency band. To run the NBS, we needed to determine the test statistic threshold. As recommended in the reference manual for the NBS connectome ( https://www.nitrc.org/ frs/download.php/5331/Reference_Manual_NBS_v1.2.pdf ), we changed this threshold value from 2 to 4 in increments of 0.1. After confirming that the results are qualitatively consistent across all the thresholds we tried, we selected 2.5 as the threshold.

Diversity profile and spectrum
We first developed a method to characterize the inter-state variance of neural interactions in terms of their frequency-specific spatial profile and local spectral property. These variance features were derived from a matrix summarizing each ROI's frequency-specific variance of functional connectivity (diversity matrix, Fig. 1 D), which was achieved through the following procedure. We estimated coherence spectra for all pairs of ROIs and obtained frequency-specific functional connectivity matrices in all task states by calculating band-limited coherence coefficients (Fig.1A). We assessed strength of an ROI's neural interactions with all the rest by sum of coherence coefficients (Fig.1B), and we defined a measure of diversity as the standardized variance of the strength across all task states (Fig.1C, D). The strength quantifies the degree to which a ROI coordinates with the rest of the brain, while the diversity measures how much the region differentiates the coordination across task states. We used row vectors of strength and diversity matrices to characterize whole brain profiles (strength or diversity profile), while we used column vectors to extract their spectral features (strength or diversity spectrum).

Diversity profile determines four frequency ranges
We attempted to clarify whether cognitive tasks differentiate regions' neural interactions in a frequency-specific manner. We evaluated the similarity of diversity profiles between different frequency bands through Spearman's correlation ( Fig 2 A). Then, we tested whether there are clusters of frequency bands based on similarity. We applied a cluster detection algorithm called the Louvain method to the correlation matrix at an optimal network resolution scale (see Method) and identified four clusters of contiguous frequency ranges: (1) 0.01-0.15 Hz, (2) 0.15-0.37 Hz, (3) 0.37-0.53 Hz, and (4) 0.53-0.70 Hz. We also confirmed the network has a statistically significant cluster organization ( p =5 . 7 × 10 −172 ).

Diversity spectra detect functional modules
The spectrum of regional diversity characterizes the regional frequency-dependency of the inter-state variance. As such, its interregional correlation enables the clustering of regions based on the frequency-dependency of diversity. By applying the module detection analysis to the correlation matrix, we found three separate modules ( Fig. 3 A, B, p =1 . 0 × 10 −498 ). The first module contains the pre-and postcentral gyri, the superior temporal gyrus, and some cingulate cortical areas (Module 1). The second module includes the posterior parietal cortices, the posterior cortices, and the cerebellum (Module 2). The third module mainly consists of frontal cortical regions (Module 3). The

Fig. 1. Strength and diversity profiles and spectra. (A)
We estimated coherence for all pairs of ROIs in each task state with each individual data and averaged them at each frequency band to obtain a frequencyspecific functional connectivity matrix. (B) We calculated the nodal strength as the sum of each row vector of the frequency-specific functional connectivity matrix and obtained a 68-by-264 matrix for each state. We call row vectors "strength profiles " in frequency bands, and we call column vectors "strength spectra. " (C) We estimated the variance of the nodal strengths across states for each node at each frequency, and we obtained the variance matrix. (D) We standardized the variance matrix such that each row vector has a mean of 0 and a standard deviation of 1, and we made the diversity matrix. We call the row vectors "diversity profiles " and the column vectors "diversity spectra. " fourth module is mostly formed by the temporal, subcortical, cerebellar, and cingulate cortical regions (Module 4). Furthermore, we compared this partition of brain modules with another partition obtained based on the functional connectivity by using NMI ( Fig. 3 C). We found that the NMI value was significantly higher than 0 (0.5623, p < 10 − 5 ), indicating that the modular partition is consistent with those based on functional connectivity.
We also compared the observed modular partition with one previously reported by Power et al. (2011) . We confirmed that both modular partitions are significantly similar (NMI = 0.356, p < 10 − 5 ). By referring to the module labels in the previous study, we found that Module 1 included many areas of the somatosensory and auditory networks, as well as some of the cingulo-opercular and ventral attention networks. Module 2 contained all brain regions in the visual and cerebellar networks, and many of the dorsal and ventral attention networks, on top of some areas in the somatosensory, auditory, frontoparietal, and default networks. Module 3 included many regions forming salience and frontoparietal networks, as well as some of those forming cingulo-opercular, default, and attention networks. Module 4 encompassed all brain regions of subcortical networks; it also included small portions of default, auditory, salience, cingulo-opercular, and somatosensory networks. Table 1 summarizes the correspondence between modules identified here and those reported by Power et al. (2011) .

The most diverse module switches between high and low-frequency ranges
We next sought to recapture the frequency-specificity in the statedependent modulation of functional connectivity in terms of identified modules' diversity. We defined the diversity spectrum of a module as the mean diversity spectra of regions assigned to the module (see Method and Fig. 1 ). Fig. 4 A shows the diversity spectra of modules. We found that the diversity values of Modules 1 and 2 are the highest in the verylow frequency range ( < 0.1 Hz). However, it decreased almost monotonically along frequency and got lower than the diversity of Module 3 above 0.15 Hz. On the other hand, the diversity of Module 3 increased along frequency and became the highest above 0.15 Hz. The diversity of Module 4 was as low as or lower than those in other modules across all frequency. To ensure the statistical significance, we averaged the diversity spectra of modules in the frequency ranges identified based on similarity of diversity profiles (see above) and compared the mean values among modules. Specifically, we conducted ANOVA for the set of frequency ranges to test if there was any significant difference in the diversity across modules and states. As a result, we found significant main effects for modules, states, and their interactions ( p < 0.0001 for all). By conducting post-hoc multiple comparisons among modules, we found that Module 1 and 2 showed significantly higher diversity than other modules at 0.01-0.15 Hz, but values of these modules were significantly lower than that of Module 3 at 0.15-0.7 Hz ( Fig. 4 B). We confirmed that the diversity of Module 4 was significantly lower than other modules in all frequency bands below 0.53 Hz ( p < 0.0001 for all comparisons). Fig. 4 C shows the change of the diversified areas across frequency bands.
We also investigated if any of the identified modules corresponds to the flexible hub reported by Cole et al. (2013) . We calculated Power et al. (2011) 's modules' diversity values and identified the most diversified module at each frequency band. As a result, the most diversified module is the sensorimotor network in 0.01-0.15 Hz; it is replaced with the salience network in the 0.15-0.70 Hz ( Fig. 4 D). In 0.01-0.15 Hz, the visual system, visual attention system, and dorsal attention system also show high diversity values comparable to that of the sensorimotor system, while the frontoparietal system has high diversity values comparable to that of the frontoparietal system in 0.15-0.70 Hz.

The relative connectivity strength across states is specific to the functional module
We investigated whether all modules differentiate connectivity strengths across states in the same way or in module-specific manners. To do so, we standardized the functional connectivity strength across states for each node and averaged them in each module. Fig. 5 A shows the mean standardized strengths in four frequency ranges. We evaluated the dissimilarity of the strengths between all pairs of modules by calculating distances in the five-dimensional Euclidean space; the strength    in each state was regarded as the coordinate of one spatial dimension. The mean dissimilarity across all subjects was statistically compared by those obtained by permuting module labels to test the null hypothesis of no dissimilarity. We found that the inter-state variance (ISV) of Module 3 was significantly dissimilar to those of modules consisting of parieto-occipital regions, Module 1 or 2. For example, the difference of ISV between Module 1 and 3 is significant in the first and the second frequency bands (0.01-0.15 Hz and 0.15-0.37 Hz). We found a difference between Module 2 and 3, in the first and the third frequency bands (0.01-0.15 Hz and 0.37-0.53 Hz). We also observed significant a difference in ISV between Module 1 and 2 in the first and the third frequency bands. ISV of Module 4 also shows a significant difference from those of Modules 1 and 2.

Task-specific functional connectivity analysis
We further examined whether there is any task-specific functional connectivity (FC) in any identified frequency band. Thus, to test for any significant task-specific FC patterns, we conducted a network-based statistical analysis ( Zalesky et al., 2010 ). Fig. 6 shows FCs whose strengths increased or decreased in a task-specific manner. We found significant FCs in visual, auditory, heart, and resting conditions, while we could not find such FCs in the self-task state. Interestingly, we found that the network wiring patterns in the visual condition depend on the frequency band where functional connectivity is estimated.

Module detection using EEG-fMRI coherence spectra
Lastly, we examined brain network organization in relation to EEG-fMRI coherence. Specifically, we tested whether we can find consistent brain modules using the similarity of the coherence among ROIs. As a result, we successfully detected brain modules that are significantly similar to those found by the similarity of diversity spectra ( Fig. 7 A, NMI = 0.5383, p < 10 − 5 ). Next, we asked whether fMRI signals reflect similar electrophysiological signals across all the fMRI frequency bands. To do so, we obtained the mean EEG-fMRI coherence spectrum for each module at each frequency band ( Fig. 7 B). We found that spectra within the same frequency band are significantly correlated among modules ( Fig. 7 C). While we also found significant correlations between different frequency bands; these correlations were significantly lower than those within the same frequency band, indicating that EEG frequency components underlying fMRI signals depend on the fMRI signal's frequency band. To test this, we obtained the mean spectra for all modules at each frequency band and examined the statistical significance of coherence for each EEG frequency. As a result, we found that each fMRI frequency band shows significant coherence with a specific set of EEG frequency components ( Fig. 7 D).

Discussion
We aimed to clarify whether cognitive states differentiate wholebrain functional connectivity in a frequency-dependent manner. We recorded whole brain activities with multiband fMRI at a high sampling frequency of 1.4 Hz during resting and four task states. We quantified the variance of the regional connectivity strengths, diversity, in the range of 0.01-0.7 Hz by dividing the frequency range into 68 narrow frequency bands. By calculating the correlation matrix of the network-wide diversity profiles between all pairs of frequency bands, we found that the matrix is modular so that it can be decomposed to four contiguous clusters of frequency bands (0.01-0.15 Hz, 0.15-0.37 Hz, 0.37-0.53 Hz, and Fig. 6. Task-specific functional connectivity varies across frequency bands We detected increased and decreased functional connectivity in each task at each frequency band through the network-based statistical analysis. Magenta edges indicate functional connectivity significantly strengthened, while cyan ones indicate those significantly weakened. We only found such functional connectivity in (A) visual task state, (B) auditory task state, (C) heart task state, and (D) resting task state. We only show the frequency bands where we found significant task-specific functional connectivity. 0.53-0.70 Hz). Then, we examined whether the frequency-specificity of the diversity relates to functional network organization. We evaluated the similarity of diversity spectra between all pairs of regions and found three modules. The sets of regions forming the modules were substantially similar to those forming functional modules of the whole-brain network. Interestingly, the modules showing the highest diversity were switched from Modules 1 and 2, which largely consist of posterior regions, to Module 3, which is mainly formed by frontal regions, across frequency. Furthermore, we found that Modules 1 and 2 show a dissimilar shifting pattern of connectivity strengths across states compared to those of Module 3 below 0.15 Hz and above 0.15 Hz. Moreover, the EEG-fMRI coherence analysis showed that fMRI signal fluctuations at 0.01-0.15 Hz correlate with the power of EEG waves different from those correlated with the fluctuations above 0.15 Hz.
To quantify the variance of functional connectivity, we used the variance of band-limited coherence strengths across resting and four task states. This measure was selected to characterize the state-general modulation of the functional connectivity with reducing the state-specific effect. It has been shown that there are both task-general and task-specific modulations of the functional connectivity ( Cole et al., 2014 ). The taskgeneral modulation refers to changes in strengths common across tasks, seen in a large number of functional connectivity. On the other hand, the task-specific modulation refers to changes specific to the tasks, occurring locally. Because of the task-specific modulation, the standard pair-wise comparison between a single task state and the resting state can introduce biases in evaluating the state-general modulations. Thus, we selected four conditions requiring subjects to maintaining different foci and used variance across multiple tasks as the measure of the stategeneral modulation.
We selected the steady-state task design rather than the event-related design. The more sophisticated connectivity analysis using event-related designs is better for capturing instantaneous neural interactions in the time domain than are those that use steady-state designs. This is because the connectivity analysis using event-related designs enables detection of both event-dependent and event-independent connectivity by using techniques such as generalized psychophysiological interactions ( McLaren et al., 2012 ). However, it remains unknown in what way the event-related design can be combined with techniques estimating frequency-dependent functional connectivity. Therefore, instead of the event-related design, we used a steady-state design that requires subjects to keep their cognitive states stable in order to assume that neural interactions during the task are as stable as possible.
The fMRI's blood-oxygenation-level-dependent (BOLD) signal can be modeled as time-series of neural events convolved with a canonical hemodynamic response function that basically has the property of a low-pass filter ( Friston et al., 2000( Friston et al., , 1998. Because of the sluggishness of the hemodynamic response, it has been thought that fMRI cannot measure neural activities occurring in higher frequency ranges even if the sampling rate is improved. However, this point of view is challenged by the existence of a frequency cluster consisting of only highfrequency components ( > 0.53 Hz) that are strongly attenuated by the canonical hemodynamic response function. In line with this, a recent paper reported fast oscillatory fMRI signals in response to stimuli oscillating at up to 0.75 Hz and showed by a simulation that BOLD responses Fig. 7. EEG-fMRI correspondence is specific to the fMRI frequency band (A) We detected four brain modules based on the similarity of EEG-fMRI coherence spectra (red, blue, yellow, and orange colors). Although the number of detected modules is different from what we found using the similarity of diversity spectra, the similarity of network partition quantified by normalized mutual information is significant ( p < 0.001). (B) The mean EEG-fMRI coherence spectrum of each module at each frequency band. (C) The mean Spearman's correlation of intra-band/interband pairs of coherence spectra. (D) The significant coherence in the averaged EEG-fMRI coherence spectrum across all modules. become quicker for rapidly changing stimuli ( Lewis et al., 2016 ). Another experimental study has shown that the coupling between cerebral blood flow and volume changes during dynamic stimulation ( Simon and Buxton, 2015 ), indicating that hemodynamic response can be specific to the time-scale of stimuli. Furthermore, additional previous studies demonstrated that fMRI-based functional connectivity at rest persists in the frequency range above 0.1 Hz ( Chen and Glover, 2015 ;Gohel and Biswal, 2015 ;Liao et al., 2013 ). In the current study, we also revealed the existence of task-specific functional connectivity in this frequency range. Remarkably, task-specific functional connectivity observed during both the visual task and the heart task state showed different spatial patterns between frequency range below 0.15 Hz and that above 0.15 Hz. Moreover, we confirmed that fMRI signal components significantly correlate with EEG signals all over the frequency bands ( Figure  S4). Thus, despite fMRI's low signal-to-noise ratio, it is likely that it can capture neural interactions in frequency ranges higher than 0.1 Hz so that we can investigate whole-brain functional organization in the slow-delta frequency range, which has been shown to be involved in attentional selection ( Lakatos et al., 2008 ).
It has been reported that the fMRI signal correlates with fluctuations of EEG power ( He et al., 2008 ;Sadaghiani et al., 2010 ). A recent study has shown that EEG and fMRI dynamic functional connectivity exhibit shared dynamics ( Wirsich et al., 2020 ). However, to the best of our knowledge, previous studies examined the relationship using fMRI signals obtained using conventional time-resolution (TR~2 s). Utilizing the high temporal resolution, we investigated the correlation up to 0.7 Hz and found that the frequency of EEG power correlation with fMRI is dependent on the fMRI's frequency ranges ( Fig. 7 B). It has been shown that fMRI signals in the frequency range below 0.1 Hz have a significant correlation with the EEG power time-courses of frequency components around 10 and 20 Hz ( Sadaghiani et al., 2010 ). Similarly, we found that the fMRI signals in the first frequency band showed significant coherence with the EEG power time-courses of components whose frequencies are around 10, 17.1, and 21.4 Hz. Yet, the fMRI signals in the frequency bands above 0.15 Hz demonstrated a differential relationship with EEG frequency. For example, significant coherence with the fMRI signals in the second and fourth frequency bands occurred more broadly. In the second frequency band, we found significant coherence in the delta (1-4.3 Hz), high alpha (11.4-17.1 Hz), and high beta bands (27.1-30 Hz). The coherence in the fourth frequency band occurred in the theta (4.3-7.1 Hz) and high alpha (11.4-17.1 Hz). In the third frequency band, significant coherence with the fMRI signals occurred selectively around 1.4 Hz, 8.6 Hz, and 30 Hz. Notably, none of these high fMRI frequency bands' signals showed significant coherence around 10 and 21.4 Hz. This indicates that fMRI frequency components in the high-frequency range reflect different frequency components of electrophysiological activity from those in the low-frequency range, which are typically used for functional connectivity analysis. It should be noted that we obtained our results by averaging the fMRI-EEG coherence across all EEG channels to summarize the results. A potential caveat is that the summarized results for different EEG frequency bands may reflect different contributions from different channels. Moreover, our results indicate that task-specific functional connectivity patterns will vary with fMRI frequency components. In sum, our results support the notion that the characterization of functional connectivity in a frequency-specific manner is a promising approach to obtain rich information about brain network organizations that cannot be captured by conventional broadband characterization. We assumed that the change in signal-to-noise ratio (SNR) along frequency is similar across distinct regions. However, since our task design is not suitable for estimating task-modulated activation, we could not validate this assumption. Future studies exploring the network organization in frequency ranges above 0.1 Hz should clarify the SNR changes' inter-regional relationship.
We found four brain modules based on the correlation matrix of the diversity spectra. We confirmed that the modularity, a quantity evaluating the degree to which each module is separate from others, was much higher than the canonical criterion for judging if a network has a modular structure. The regional assignments to the modules were significantly similar to those of the modules identified based on functional connectivity. We also ensured that the regional assignments are significantly similar to those reported in the previous study ( Power et al., 2011 ). These results demonstrate that the frequency-specificity of the inter-state variance is similar among regions forming the same functional modules. Interestingly, we found that the modules consisting of posterior regions show the highest diversity in 0.01-0.15 Hz; whereas, a different module formed by frontal regions shows the highest diversity above 0.15 Hz. This relationship between these modules was also observed by comparing mean frequency-specific functional connectivity strengths. Altogether, the results indicate that a module playing a key role in the network-wide interactions switches between posterior and frontal modules across frequency.
Because each module consists of regions playing specific functional roles in the brain, modules might show different task-specific modulations of functional connectivity. However, as mentioned above, such effects cannot be captured with the diversity. Thus, we further investigated the relative connectivity strengths across states. We found the patterns in Module 3, a module consisting of frontal regions, were significantly dissimilar to those in Modules 1 and 2, mainly including parietooccipital regions, both in the range of 0.01-0.15 Hz and above 0.15 Hz. Taking it into account that in the low-frequency range (0.01-0.15 Hz) the modules with parieto-occipital regions contribute more to the brainwide neural interactions than does Module 3 and that they are switched in the high-frequency range (above 0.15 Hz), it is likely that neural interactions in these frequency ranges play distinct roles in brain function. Interestingly, the spatial layout of these modules also implies a functional distinction between frequency ranges; Module 1 and 2 basically consist of sensor and motor areas, while Module 3 is mainly formed by higher-order association areas. Expecting the inter-frequency functional distinctions to be highlighted by the difference in the attentional modulation, we used 2 exteroceptive and 2 interoceptive tasks in this study. However, we could not find any commonalities amongst the 2 tasks of a given type; rather, each task seems to have a specific modulatory effect on the functional network organization. It will be important for future research to investigate the association of fMRI-based frequency-specific connectivity by using task sets designed to perturb brain states in many directions in order to clarify the canonical functions of neural interactions in the frequency ranges.
Recent advances in the framework for multilayer networks have enabled a new vista for network science. This framework can deal with networks whose nodes are interacted in multiple different ways, and it can provide an improved characterization of networks by taking the multiple interactions into account. The brain is one of the multilayer networks, for it has several ways of interactions among its parts, such as timevarying interactions and frequency-dependent interactions. The multilayer brain modeling has been getting popular recently, while most of the application is that of time-varying interactions. The main shortcoming of constructing a multi-frequency-layer network with fMRI is its poor frequency-resolution. However, the current study demonstrates that fast fMRI techniques can derive rich information about frequency-specific functional networking at least in the ultraslow frequency range. Furthermore, our previous study demonstrated that the multilayer modeling can provide unique information about the spatial distribution of brain hubs, bringing better discriminability between the healthy and schizophrenic populations ( De Domenico et al., 2016 ). Since the frequency domain would contain different information about inter-regional interactions compared to those in the time domain, our results encourage the multifrequency-layer modeling of fMRI-based networks for a better understanding of the large-scale functional brain organization.

Declaration of Competing Interest
None