Natural viewing experiment

We reanalyzed data from our recent natural viewing experiment (J??skel?inen et al., ²⁰⁰⁸). fMRI was measured from 15 healthy right-handed volunteers (age range 19?44, seven females, data from three subjects was discarded previously due to bad quality) when they watched the movie ?Crash? (directed by Paul Haggis, Lions Gate Films, 2005). Ethical approval was granted by the Coordinating Ethical Committee of Hospital District of Helsinki and Uusimaa, Finland. The study was run in accordance with the guidelines of the Helsinki declaration, and a written voluntary consent was obtained prior to participation in the study from each subject. A General Electric Signa 3 Tesla MRI scanner with Excite upgrade located in the Advanced Magnetic Imaging Center at the Helsinki University of Technology was used in recordings. The movie was selected because of its high quality (reflected in three Academy Awards in 2006: best picture, film edit, writing/screenplay), and because it had been identified in the initial screening process as one eliciting strong emotions in viewers. Only those subjects who had not seen the movie previously and who were fluent in English were accepted to the study since the language of the movie was English. In the beginning of the experiment, each subject watched the first 72?min of the movie on a LCD screen in a dimly lit room with a five-speaker surround-sound setting. The movie was then discontinued and the subject was prepared for the fMRI scans, which took about 30?min. In the beginning of the scanning process, initial localizer scan and 45 axial T1-weighted anatomical MRI slices (slice thickness 5?mm, in-plane 2?mm???2?mm) were acquired. During functional imaging, the movie was projected onto a translucent screen, which the subjects viewed through an angled mirror, and the audio track of the movie was played to the subjects via a pneumatic audio system with earplugs. The last 36?min of the movie were shown for a subject in two parts lasting 14 and 22?min, respectively, with a pause of about 30?s in between. The break was due to limitations in the number of slices that could be acquired during a single run in the MRI scanner that was used (J??skel?inen et al., ²⁰⁰⁸). Altogether 630 whole-brain T2* weighted gradient echo EPI volumes (TR?=?3400?ms, TE?=?32?ms, flip angle 90?, 52 axial slices, 3?mm slice thickness without gaps, 3?mm???3?mm in-plane resolution) were acquired during the two movie clips. T2-weighted fluid attenuated inversion recovery and high resolution 3D T1-weighted structural images were obtained immediately after the presentation. The subjects were instructed not to move during the scanning session in order to avoid motion artifacts.

Image preprocessing

We preprocessed the original fMRI experiment data of 12 subjects using the FSL software library (FMRIB, Oxford University; Smith et al., ²⁰⁰⁴). Since we wanted to avoid any discontinuities in the data, we processed and analyzed the two scanning sessions separately. First, motion correction was performed by using MCFLIRT (Jenkinson et al., ²⁰⁰²). Time-series were then detrended using a nonlinear Gaussian weighted least-squares fitting with a cutoff of 100?s, and spatial smoothing was applied with a Gaussian kernel with FWHM of 5?mm. Note that the detrending removes only a part of the low frequency signal components (we observed approximately 10% signal energy loss below the cutoff frequency) and it is reasonable to study the frequency bands below the cutoff value after the detrending. After preprocessing, images were spatially normalized to the ICBM-152 template using the FLIRT registration tool (Jenkinsson and Smith, ²⁰⁰¹). Skull-stripped functional images were first aligned (six-parameter rigid transformation) to the skull-stripped high-resolution T1-weighted image of the same subject, and then the T1-weighted image to the standard (brain only) ICBM-152 template with a 12-parameter affine transformation. The skull-stripping of functional images was performed with BET (Smith, ²⁰⁰²). T1-weighted images were skull-stripped using a customized procedure involving genetic algorithm based voxel classification (Tohka et al., ²⁰⁰⁷) and brain extraction from the voxel classified images using Brain Surface Extractor (Shattuck et al., ²⁰⁰¹).

Filter bank construction with stationary wavelet transform

Multiresolution analysis using discrete wavelet transform (DWT) is a widely used technique to analyze time-scale characteristics of discrete-time signals and to construct octave-filter banks (Mallat, ¹⁹⁹⁹). In fMRI signal analysis, wavelets have been advocated for several reasons. For example, wavelet analysis is a prominent approach when analyzing fractal signals or signals with 1/f -like frequency characteristics, both properties that have been associated with cortical resting-state fMRI time-series (Bullmore et al., ²⁰⁰⁴; Maxim et al., ²⁰⁰⁵). The standard implementation of DWT uses a sub-band coding scheme consisting of successive stages of filtering and down-sampling (for a detailed description, see Vetterli and Kova?evi?, ²⁰⁰⁷). However, the down-sampling makes DWT time variant (i.e., DWT is not time-invariant). This means that even a minor shift in the input signal can cause large changes in the transformed signal (Bradley, ²⁰⁰³). The lack of time-invariance is problematic for the ISC analysis because even a small difference in hemodynamic responses between two subjects can lead to inconsistencies in the estimation of correlation between the subjects? time courses.

We constructed a wavelet filter bank using a closely related technique called SWT, also known as undecimated wavelet transform, redundant wavelet transform, or maximum-overlap DWT (MODWT) in the literature, to overcome the problems associated with DWT. Previously, SWT has been used to investigate resting-state fMRI data by Achard et al. (²⁰⁰⁶). SWT is a time-invariant transformation and therefore a good choice for frequency-specific ISC analysis. In SWT, a signal of interest is not downsampled but instead the analysis filter coefficients are upsampled after successive filtering steps (Nason and Silverman, ¹⁹⁹⁵). Because the number of samples in successive decomposition levels remains constant, SWT leads to an overdetermined representation of the original time-series ? hence the name redundant wavelet transform. Next, we briefly describe the algorithm to perform SWT (known as algorithme ? trous; Vetterli and Kova?evi?, ²⁰⁰⁷) to highlight its good properties for the frequency-specific ISC analysis. Let H denote a unit length low-pass filter whose coefficients {h[n]} satisfy the following internal orthogonality relation (Nason and Silverman, ¹⁹⁹⁵):

and let G denote a high-pass filter having coefficients {g[n]} satisfying:

Also the filter G satisfies (1) and together the filters H and G are called a quadrature mirror filter (QMF) pair. To derive SWT, two operations are required. First, we define a filtering operation applied to a sequence x as a convolution sum:

where boundary conditions are treated using periodic signal extension. Second, we define the operator Z that alternates zeros in a discrete sequence x in such a way that for all integers j, (Zx)[2j]?=?x[j] and (Zx)[2j?+?1]?=?0. Now, assume that a low-pass filter H^(r) has the coefficients Z^rh, i.e., h^r[k]?=?h[j] if k?=?2^rj, and h^r[k]?=?0 if k?? 2^rj. Clearly, the filter H^(r) can be obtained by inserting zeros between every adjacent pair of elements of the filter H^(r?1). Similarly to H^(r), we define a high-pass filter G^(r) having coefficients Z^rg. Then, denoting the original signal by c⁰[n], SWT is computed by applying recursively the following operations for j?=?1,2,?,J:

Here, c^j[n] and d^j[n] stand for approximation and detail coefficients at scale j, respectively, and J is the number of decomposition scales.

We wrote customized MATLAB code to implement 4-scale (J?=?4) SWT decomposition using Daubechies 2 QMF pair to filter the fMRI time-series of each subject. With the 4-scale filter bank, five frequency sub-bands were obtained, consisting of four sets of detail (d¹[n], d²[n],d³[n], d⁴[n]) and one set of approximation coefficients (c⁴[n]). The decomposition scheme is illustrated in Figure 2. For notational convenience, we denote decompositions by s^k[n]?=?d^k[n], for k?=?1,2,3,4, and s^k[n]?=?c^k[n] for k?=?0,5.

Table 1 displays approximate frequency characteristics and the corresponding time-scales of the fMRI time-series after the wavelet decomposition. The successive splitting of the low pass-filtered coefficients into new low- and high-pass coefficients leads to a logarithmic scale of the frequency sub-bands.

Frequency-specific intersubject correlation analysis

We used Pearson's correlation coefficient to derive our multi-subject similarity measure. First, we calculated voxel-wise temporal correlation between every pair of subjects within each frequency band as:

where N is a total number of samples in time-series, r_ij is the sample correlation coefficient between the time-series s_i and s_j obtained from the ith and jth subject, respectively, and si? and sj? denote the mean values of s_i and s_j. Note that superscripts denoting frequency sub-bands are omitted from Eq. 2 for simplicity. Then, we combined r_ij values from all subject pairs into a single ISC statistic by averaging: where m is the number of subjects. Thus, as m was 12 in our study, we averaged correlation coefficients of (12² ??12)/2?=?66 subject pairs. The statistical inference is complicated by the dependency of the 66 correlation coefficients. Wilson et al. (²⁰⁰⁸) suggested a t-test based approach for testing if the sum of Z-transformed intersubject correlation coefficients is different from 0. The parametric form of the null distribution was verified using a goodness-of-fit evaluation based on permutation testing. Because we had adequate computational resources available and the statistic has an intuitive meaning when examining intersubject similarities, we performed a fully nonparametric voxel-wise permutation test for r? statistic similar to the goodness of fit evaluation of Wilson et al. (²⁰⁰⁸). The test was performed against the null hypothesis that r? statistic is the same as for unstructured data, which would be expected if there was no ISC present. To generate the permutation distribution, we circularly shifted each subject's time-series by random amount so that they were no longer aligned in time across the subjects, and then calculated r? statistic. This way we could account for temporal autocorrelations present in the fMRI data. Because calculation of all the possible time shift combinations would be computationally prohibitive, we approximated the full permutation distribution with A?=?100,000,000 realizations. Sampling was randomized over every brain voxel, frequency sub-band and shifting point without any restrictions. We corrected the resulting p-values using false discovery rate (FDR) based multiple comparisons correction with independence or positive dependence assumption (Benjamini and Hochberg, ¹⁹⁹⁵; Nichols and Hayasaka, ²⁰⁰³). FDR corrected thresholds were calculated separately for each frequency sub-band. We calculated two permutation distributions to threshold frequency-dependent ISC maps of the both scanning sessions (separate distributions are needed because of different numbers of data points available to compute the correlation coefficients).

In addition to r?, we calculated indexed ISC maps to visualize the frequency bands showing the highest ISC (according to Eq. 3). For these maps, the frequency band index was calculated for each voxel as follows:

where a?{1,2,3,4,5} indexes a frequency band. The maps were thresholded according to the lowest FDR corrected threshold among the frequency bands, since FDR correction leads to a separate threshold for each frequency band. Investigation of the frequency-specific permutation distributions showed that standard deviations in lower frequency bands were greater compared to those of higher frequency bands. This was expected, because time resolution decreases when moving from higher frequencies towards the lower ones. Difference in standard deviations across the frequency bands is not problematic when investigating individual ISC maps because FDR correction is performed separately for each frequency band. When investigating results simultaneously across several frequency bands, it must be noticed that it is more likely to see spurious ISC values in lower frequency bands compared to higher ones. However, the investigation of the mean values of the frequency-specific permutation distributions showed that the ISC values were not inflated due to filtering (the range of the mean distribution values in different frequency bands was from ?4.06???10^?5 to ?1.96???10^?6 for the first scanning session, and from ?2.67???10^?5 to 2.34???10^?5 for the second scanning session), indicating that it is reasonable to investigate multi-frequency mappings as proposed in Eq. 4.

We also investigated variability of the frequency-specific ISC results across the two scanning sessions. We quantified across-session reproducibility within each frequency band by calculating a mean absolute error between the two sessions across every brain voxel as follows:

where r1? and r2? stand for r? statistic (Eq. 3) across the first and the second scanning session, respectively, and B denotes the set of all brain voxels. There are many sources of variability contributing to the absolute error, including different stimuli during two scanning sessions. Therefore, the interpretation of the variability measure is not straight-forward. However, we wished to provide a quantitative measure of the scanning session differences in addition to qualitative inspection of the statistical parameter maps.

Frequency-band comparison

We designed a permutation test also to assess whether ISC in one frequency band differs significantly from ISC in another frequency band. First, we used the modified Pearson?Filon statistic based on Fisher's z transformation (ZPF) (Raghunathan et al., ¹⁹⁹⁶) to compare ISC across frequency bands between two subjects, which is a recommended test statistic for two, non-overlapping dependent correlations (Krishnamoorthy and Xia, ²⁰⁰⁷). A selection of the statistic was based on the assumption that subject-pairwise correlation coefficients may be correlated in distinct frequency sub-bands. To obtain ZPF statistic, all pairwise correlation coefficients (obtained by Eq. 2) were first transformed into a normally distributed form through Fisher's z transformation:

where superscript a refers to the frequency band. The difference between the two z transformed correlation coefficients was then transformed into a ZPF statistic through the following transformation (Raghunathan et al., ¹⁹⁹⁶):

With certain assumptions, the ZPF statistic is distributed according to the standard normal distribution N(0,1) under the null hypothesis of no difference between the correlation coefficients. Here, notation rijba stands for correlation between the time-courses of ith subject from frequency band a and jth subject from frequency band b, and large-scale covariance is given by Raghunathan et al. (¹⁹⁹⁶):

We derived a group-level statistic by combining pairwise statistics of all subject pairs. We define a sum ZPF statistic as follows:

where m?=?12 is the total number of subjects. No exact parametric assumptions can be made about ZPF?ab statistic due to correlations between the individual pairings without the assumptions about the true, underlying correlation coefficient values between the time-series of subject pairs (see Brien et al., ¹⁹⁸⁴ for the formula of the covariance matrix between different components of a Z-transformed correlation matrix). Therefore, we performed a nonparametric permutation test under the null hypothesis that each ZPF value is drawn from a symmetric distribution with zero mean. This occurs exactly when there is no difference in ISC between two frequency bands. The subject pairs are not independent so this cannot be used to ensure the exchangeability under the null-hypothesis. However, in Eq. 6, the denominator remains the same if the roles of the two frequency bands are interchanged and under the null hypothesis the covariance structure of {zija} and {zijb} are the same. Therefore, under the null hypothesis, the requirement of weak exchangeability (see Good, ²⁰⁰⁵, p.269) is fulfilled. We generated the permutation distribution by randomly flipping the sign of 66 pairwise ZPFijab statistic before calculating Eq. 7. We approximated the distribution with a subsample of 25,000 random labelings (out of 2⁶⁶ possible labelings) and saved maximal and minimal statistic over the entire image corresponding to each labeling (due to symmetry of the distribution, we obtained maximal statistic of 50,000 samples for testing both ZPF?ab and ZPF?ba with this procedure). Note that maximal statistic automatically accounts for multiple comparisons problem by controlling the family wise error rate (FWER, Nichols and Holmes, ²⁰⁰¹), and thus we obtained critical values as the 1+??/A?. largest value of the permutation distribution, where A is the number of samples in the maximal statistic distribution and ? is a significance level. All possible ZPF statistical maps were formed to allow systematic comparison of the frequency bands.

Frequency-specific ISC

Figures 3A?C show frequency-specific ISC maps (according to Eq. 3) for the high frequency band, the mid-frequency band, and the low frequency band. Two additional frequency bands, ISC maps across the full frequency spectrum, and the corresponding results of the second scanning session are presented as Supplementary Material. Visual inspection revealed extensive ISCs in both sensory and association areas, thus supporting the previous observations (Hasson et al., ²⁰⁰⁴, ^2008b; J??skel?inen et al., ²⁰⁰⁸). In the highest frequency band (Figure 3A), ISC was concentrated in visual cortex (the maximum value was 0.34), including occipito-temporal fusiform cortex. In the middle-frequency band (Figure 3B), extensive ISC was present in both auditory and visual areas (showing the maximum value of 0.57). Areas in temporal lobe and visual cortex showed extensive ISC also in the low-frequency band (Figure 3C). In this band, ISC was very high, especially in temporal cortical areas (0.79 at maximum). Interestingly, significant ISC in the low-frequency band was found also in higher-order areas, including frontal pole and superior frontal gyrus. Further, several subcortical areas showed ISC in the lowest frequency band.

Figure 4 shows an indexed multi-frequency ISC map (Eq. 4, p?<?0.001, FDR-corrected) across the two scanning sessions. Frequency bands showing strongest ISC are indicated with color coding. The highest frequency band (0.07?0.15?Hz) was omitted because it did not show maximum ISC anywhere in the brain. This visualization scheme provides an overview of the major frequency component showing ISC in different brain areas. This presentation is not suitable for the comparison of the frequency bands, because statistical inference was performed individually within each sub-band.

During the first scanning session (Figure 4A), several occipital areas, including visual areas, occipital fusiform gyrus, and superior parietal lobule, showed maximum ISC in the 0.02?0.04?Hz range (dark yellow). In some areas, such as in lingual gyrus and temporal occipital fusiform cortex, ISC was evident mostly in the 0.01?0.02?Hz range (blue). Some areas, such as occipital pole, temporal occipital fusiform cortex and precuneus cortex showed ISC across a wider frequency range (high ISC in both 0.01?0.02 and 0.02?0.04?Hz bands). ISC in temporal areas was mostly concentrated within the lowest frequency band, i.e., 0?0.01?Hz range (green). Further, several other areas, such as angular gyrus, supramarginal gyrus, and frontal areas (e.g., orbito-frontal cortex) showed maximum ISC mainly in the lowest frequency band. The analysis results across the second scanning session (Figure 4B) appeared to be qualitatively similar to those obtained across the first scanning session.

Table 2 presents reproducibility analysis (see Eq. 5) results across the two scanning sessions. Results across the full spectrum as well as across each individual frequency sub-band are presented. Reasonably small mean absolute errors were obtained especially across the full spectrum (0?0.15?Hz), and the higher frequency bands (0.07?0.15, 0.04?0.07, and 0.02?0.04?Hz), indicating that ISCs were reasonably well reproducible across the sessions in these frequency bands. Note, however, that there was a gradual increase in q values when moving from higher frequencies towards the lower ones [the highest variability q?=?0.049 was obtained in the lowest frequency band (0?0.01?Hz)]. This analysis does not reveal the sources of variability and their respective contribution to the total variability over the two scanning sessions. More advanced methods would be needed to allow conclusions about different sources of variability.

Frequency band comparison

In Figure 5, three selected sum ZPF maps (Eq. 7), p?<?0.05, FWER corrected) are presented to compare ISCs in the lowest frequency band against the ISCs in higher frequency bands. The rest of the results, including the remaining sum ZPF maps of the first scanning session and all the sum ZPF maps of the second scanning session are provided as Supplementary Material. Figure 5A shows the comparison between the lowest frequency band against the highest frequency band. Warm colors indicate where the low-frequency ISC was significantly greater than the high-frequency ISC. Cold colors show where the ISC in the high-frequency band was significantly greater than the low-frequency ISC. Low-frequency ISC was dominant in temporal cortical areas, and also in some visual and medial frontal cortical areas. Significantly stronger ISCs were not found in the highest compared with the lowest frequency band.

Figure 5B shows the comparison between the lowest and the middle frequency bands. Again, the low-frequency band showed significantly higher ISC in temporal and frontal areas. In the visual cortical areas, ISC was more dominant in the mid-frequency band compared to the lowest-frequency band. Comparison of the two lowest frequency bands is presented in Figure 5C. ISCs in the lowest frequency band were not significantly stronger in most brain areas compared to the slightly higher frequency band. This implies, together with the other comparisons, that ISC in frontal and temporal regions was most prominent across the 0?0.02?Hz frequency range. Visual cortex showed significantly higher ISC in the 0.01?0.02?Hz range compared to the lowest frequency band, implying (together with the other comparisons) that visual cortical ISC was concentrated in the 0.01?0.04?Hz range.

Summary

Figure 3 revealed that the frequency range showing significant ISC contracted when moving from lower-order sensory areas towards higher-order cortical areas. In visual cortex, ISC was present across all frequency bands. In temporal areas ISC occurred in all but the highest (0.07?0.15?Hz) frequency band, and frontal cortical ISC was present only in the two lowest frequency bands (0?0.01 and 0.01?0.02?Hz). Further, the combined ISC maps across all frequency bands (Figure 4), and the results from the specific contrasts between the different frequency bands (Figure 5), revealed that the strongest ISC component occurs at progressively lower frequencies as one moves from lower-order sensory visual cortical areas towards functionally higher-order cortical areas such as the anterior frontal lobes. However, it should be noted that even in the visual cortical areas, the dominant ISC component occurred at relatively low frequencies (0.01?0.04?Hz). The mean ISC was relatively high in some sensory areas (visual cortical areas in Figure 3B and temporal areas in Figure 3C) in contrast to higher order areas, which showed lower ISC values (frontal cortical areas in Figure 3C). Although this may suggest that neuronal activation in frontal cortical areas was weaker compared to the sensory areas during the movie, the difference may also be explained by inter-individual variability, which might be greater in higher-order areas. Further investigation is required to assess variability of ISC between subjects across different brain areas.

There are a number of possible explanations for the spatial distribution of the ISCs occurring at different frequency bands we found in this study. A recent study, utilizing short silent movie clips, suggested that there is a hierarchy of TRWs in the human brain, with sensory visual cortical areas showing short TRWs, and the TRWs becoming progressively longer as one ascends to functionally higher-order cortical areas (Hasson et al., ^2008c). While a more limited range of cortical areas was found to be activated than in the present study, perhaps due to the unimodal nature of the stimulus and the fact that the movie clips were relatively short, the frontal and temporal cortical areas that were significantly activated (frontal eye field, posterior lateral sulcus, and temporal?parietal junction) showed the longest TRWs. Thus, it is possible that the present low-frequency ISCs in frontal and temporal areas is explained by longer TRWs in these areas. On the other hand, visual features of the movie that sensory cortical areas process occur at higher frequencies than, for instance, interpersonal communication and/or development of social relationships between movie characters that the temporal and prefrontal cortical areas should be more sensitive to (see Carrington and Bailey, ²⁰⁰⁹). However, further investigations are required to find out whether and how the low frequency frontal cortical ISCs/BOLD fluctuations specifically relate to the stimulus.

Visual cortex showed ISC extensively across all frequency bands, possibly reflecting different phenomena. For example, high frequency ISC component can be related to the TRWs and low frequency component may be partially induced by modulatory and top-down inputs (Monto et al., ²⁰⁰⁸). The frequency components can also reflect other stimuli related phenomena such as temporal correlations within stimuli itself.

Conclusion

We extended ISC analysis to frequency domain by combining wavelet analysis to non-parametric permutation methods for making group-level voxel-wise statistical inferences about frequency-band specific ISC. Computationally intensive non-parametric methods were necessary for statistical inference because of the dependency among subject pair-wise correlation coefficients that were averaged to yield the ISC-statistics. Our results suggest that several regions within the frontal and temporal lobes show ISC predominantly at low frequency bands, whereas visual cortical areas exhibit ISC also at higher frequencies. It is possible that these findings relate to recent observations of a cortical hierarchy of TRWs, or that the types of events that the temporal and prefrontal cortical areas are processing (e.g., social interactions) occur over longer time-scales than the stimulus features processed in the visual areas. Software tools to perform frequency-specific ISC analysis, together with the visualization software, are made available as open source Matlab code at Supplementary Material.

The original time-series is denoted by s⁰[n].

Mean absolute error across the whole brain is presented across the full frequency spectrum (0?0.15?Hz) as well as across every individual frequency sub-band.