Home | About | Journals | Submit | Contact Us | Français |

**|**Scientific Reports**|**PMC5341062

Formats

Article sections

- Abstract
- Mathematical description
- Results
- Discussion
- Experimental Method
- Additional Information
- Supplementary Material
- References

Authors

Related links

Sci Rep. 2017; 7: 43743.

Published online 2017 March 8. doi: 10.1038/srep43743

PMCID: PMC5341062

Gadi Goelman,^{a,}^{1} Rotem Dan,^{1,}^{2} Filip Růžička,^{3} Ondrej Bezdicek,^{3} Evžen Růžička,^{3} Jan Roth,^{3} Josef Vymazal,^{4} and Robert Jech^{3}

Received 2016 October 10; Accepted 2017 January 30.

Copyright © 2017, The Author(s)

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

This article has been cited by other articles in PMC.

We describe an analysis method that characterizes the correlation between coupled time-series functions by their frequencies and phases. It provides a unified framework for simultaneous assessment of frequency and latency of a coupled time-series. The analysis is demonstrated on resting-state functional MRI data of 34 healthy subjects. Interactions between fMRI time-series are represented by cross-correlation (with time-lag) functions. A general linear model is used on the cross-correlation functions to obtain the frequencies and phase-differences of the original time-series. We define symmetric, antisymmetric and asymmetric cross-correlation functions that correspond respectively to in-phase, 90° out-of-phase and any phase difference between a pair of time-series, where the last two were never introduced before. Seed maps of the motor system were calculated to demonstrate the strength and capabilities of the analysis. Unique types of functional connections, their dominant frequencies and phase-differences have been identified. The relation between phase-differences and time-delays is shown. The phase-differences are speculated to inform transfer-time and/or to reflect a difference in the hemodynamic response between regions that are modulated by neurotransmitters concentration. The analysis can be used with any coupled functions in many disciplines including electrophysiology, EEG or MEG in neuroscience.

The last two decades have demonstrated that brain functionality and architecture can be better understood not only by identifying localized neural activity, but also, and perhaps primarily, by recognizing its connectivity. The distinction between functional segregation and integration and the use of network measures gave new means to understanding brain functionality. Within functional integration, two main classes of connectivity have emerged – functional (directed and undirected) and effective connectivity^{1}. Functional connectivity refers to statistical dependencies amongst measured time-series, while effective connectivity rests on an explicit model of how those dependencies were caused (e.g., dynamic causal modelling^{2} and structural equation modelling^{3}). In fMRI, coherent low frequency fluctuations of the blood-oxygenation-level dependent (BOLD) signal during a resting-state (rs-fMRI) were shown to contain functional neural network information^{4}^{,5}. This information was derived from the correlations between temporal fluctuations of BOLD signals in various brain regions in the absence of external stimuli^{6}^{,7}^{,8}. Multiple resting-state networks were defined^{6}^{,9}^{,10} on the basis of such temporal correlations, and their reliability and robustness were shown at individual subjects and group levels^{11}^{,12}. These networks were shown to be in correlation with individual differences in behavioral performance^{13} and altered in neurological and psychiatric disorders^{14}. Several computational models were proposed to link BOLD signal fluctuations to neuronal communication^{15}^{,16}^{,17}.

Pearson’s correlation or independent component analyses (ICA) are the most commonly used methods to obtain the level of synchrony between time-series functions. However, these methods lack the ability to distinguish between different types of synchrony, such as those that are associated with different oscillation frequencies or those that are associated with a phase difference between the time-series. In here, we propose a new analysis method which enables to characterize functional connectivity by its frequencies and phases. The analysis is general and can be used in various data types and disciplines. In MRI, several approaches were introduced previously to obtain frequency-dependent functional connectivity. Those include coherences and partial coherences in frequency space^{18}^{,19}, undirected frequency dependent graphs^{20}, spectral coherence matrix of pairwise interactions and cluster analysis^{21}, mutual information in frequency space^{22}^{,23} and recently nonlinear coherence between multiple time-series^{24}. Other studies have included dynamic information by sliding-window analysis^{25}^{,26}^{,27}, time-frequency analysis^{28}, instantaneous phase synchronization^{29} or spontaneous coactivation patterns analysis^{30}. Several studies have shown frequency-dependency of the BOLD signal and that this dependency is spatially dependent^{23}^{,31}. The observation that functional connectivity MRI is frequency dependent is supported by several findings. For example, in Parkinson’s disease patients the resting-state functional connectivity patterns of regions in the sensorimotor system were shown to differ between “OFF” and “ON” medication states^{32}. Such differences are in line with known alterations of the frequencies in neuronal firing rates in those regions^{33}. Furthermore, dependency of network topology on frequency^{34} was recently reported. The dependency of functional connectivity on phases (besides a phase of π) however, was not introduced before. In here, a method to observe functional connectivity with arbitrary phase-difference is introduced.

We present a unified framework for simultaneous assessment of frequency and latency of coupled time-series functions. We refer to latency as the phase difference between time-series functions and describe latency and functional connectivity by the same framework and as a function of frequency. Consequently, the latency is coupled in our analysis to the frequency and is represented by the phase-differences between a pair of time-series. The proposed analysis method transfers temporal 4D data (space and time) into a connectivity space in which each ‘functional connection’ (the relation between a pair of time-series) is represented by a cross-correlation with time-lag function. Using the general linear model, the weights of the time-series functions at specific frequencies and phase-differences are estimated. The strength of the proposed analysis is demonstrated on resting-state fMRI data of 34 healthy subjects. The main uniqueness of this method relative to other approaches considering frequency information in resting-state fMRI data^{34}^{,35}^{,36}^{,37} is by providing means to characterize functional connections by their phase-differences. This allows identifying types of connections, anti-symmetric and asymmetric, that were not obtained by any other method. The paper focuses on the methodological description and not on neurobiological findings. We note however that the new functional connections contain important information, such as the functional connectivity of the cerebellum that has biological relevance.

The mathematical description below refers to resting-state functional connectivity MRI but can be modified to fit other types of coupled time-series functions in varies disciplines and particularly stimulus-driven fMRI, electrophysiology, Magnetoencephalography (MEG) or Electroencephalogram (EEG).

Assuming that the BOLD time-series signal can be approximated by a finite sum of weighted cosine and sine functions whose frequencies depend on the repetition time (TR) and the number of collected time points (N), it can be expressed as:

where *i* is a point in space (ROI or voxel), L is the number of functions used in the finite sum; are the normalized weights of the cosine and sine functions respectively and LN due to the filter used in the preprocessing step. Equation 1 is simply the Fourier transform of the temporal signal where are the ‘real’ coefficients and the ‘imaginary’ coefficients. The cross-correlation with time-lag function between two BOLD signals equals:

where *l* is the time-lag between the two BOLD signals; i,j are space indexes and denotes complex conjugate. Note that due to the normalization used in Equation 1, the Pearson’s correlation coefficient equals *CC*^{i, j}(0). The frequency spectrum of Equation 2 equals:

where FT denotes Fourier Transform. Combining Equations 1 and 3 provides an expression for the frequency spectrum of the cross-correlation function. Representing this function in the time domain (the cross-frequency terms are zero for orthogonal basis set) results with a new expression for the cross-correlation function:

This expression implies that symmetric *CC*^{i, j}(*l*) results from ‘in-phase’ weights, i.e., when both BOLD signals are either cosine or sine functions at the same frequency, antisymmetric *CC*^{i, j}(*l*) results from ‘out-of-phase’ weights, i.e., when one BOLD signal is a cosine while the other is a sine function at the same frequency and asymmetric *CC*^{i, j}(*l*) results from any other possible phase difference between the BOLD signals, i.e., when both symmetric and antisymmetric weights are significant. Note that the cross-wavelet transform^{38} and the wavelet transform coherence^{39} as well as our recent derivation^{24}, give a similar expression to Equation 4 while using the wavelet space instead of Fourier space. Using Equation 4, the Pearson’s correlation coefficient is approximated by:

and therefore equals to the in-phase contributions of the cross-correlation function at time-lag zero.

To reduce the complexity of the analysis by assigning a small number of parameters (2L), a general linear model (GLM) is applied to the *CC*^{i, j}(*l*) functions (Equation 4). We chose frequencies for the GLM analysis which cover the entire frequency spectrum conventionally used for resting-state fMRI analysis (0.01–0.1Hz). Discrete Fourier Transform theory indicates that the possible lowest frequency is 2*π*/(*N*·*TR*)=0.0104Hz and the maximum number of terms needed to cover the frequency range is 8. We tested different numbers of basis set functions (from 4 to 8) and compared their agreement (goodness of the GLM fit) with Equation 4. A basis set of 8 functions: 4 cosines and 4 sines, was found to be appropriate. The following frequencies were selected: 0.02, 0.04, 0.06 and 0.08Hz. Basis-set functions were multiplied by a window function. Three different window functions were tested and the Bartlett window was found to be best (see Supplementary Information Figure 1).

The GLM approach on the cross-correlation functions can therefore be written as:

where *B(l*) is the Bartlett window function. Consequently, the new analysis termed hereinafter “Frequency-Phase Analysis” (FPA), results with eight scalars () for each functional connection.

Equation 6 can be written to explicitly express the phase such that a functional connection with any possible phase can be obtained:

where

is the phase-difference between two BOLD signals, and are given by Equation 6 and is given by Equation 1. Note that by using the GLM (Equation 6) to obtain functional connections with their phases, we avoid the need to apply statistics on complex numbers or directly on the phases which simplifies the analysis.

Figure 1 illustrates the proposed analysis method. For each pair of BOLD signals, a cross-correlation function is calculated (*CC*^{i, j}(*l*)). This process transfers the 4D data (time and space) into “interaction space” that contains all pairwise cross-correlation functions. At the next stage, a GLM analysis is performed using a basis set of 4 cosine and 4 sine functions, each with a different frequency. The GLM analysis results with eight values that are the cosine and sine weights of the pairwise cross-correlation function. These values, for a pre-defined seed, are used to construct seed statistical parametric maps (SPMs) for each of the eight GLM weight (β/γ GLM-SPMs).

The *Frequency-Phase Analysis* (FPA) is applied here on resting-state fMRI data of 34 healthy subjects to demonstrate its ability to identify unique functional connections. Statistical parametric maps (SPMs) were calculated for each of the 8 GLM-weights. These maps, referred hereinafter as ‘GLM-SPMs’, allow to characterize the functional connectivity of seed regions according to their phases (β GLM-SPMs vs. γ GLM-SPMs) and frequencies ( GLM-SPM vs. GLM-SPM). Figure 2 shows GLM-SPMs of the left thalamus for β_{1} and β_{4} GLM-weights (SPMs for all GLM-weights are shown in Supplementary Information Figure 2). In the figure, positive and negative t-values are indicated by red and blue colors, respectively, and the left thalamus seed is shown in white. The average (across all voxels and all subjects) F-value corresponding to the goodness of the GLM fit was 13.8±0.006 (mean±standard error), indicating high agreement between the GLM weights and cross-correlation functions in most of the voxels and demonstrating that 8 parameters are suitable to fit the cross-correlation functions. For the thalamus seed, almost no significant volumes were observed for the γ (antisymmetric) GLM-SPMs. GLM-SPMs of β_{1–3} (frequencies of 0.02, 0.04 and 0.06Hz) were similar to each other but very different from β_{4} GLM-SPM (frequency of 0.08Hz, Supplementary Information Figure 2). Note that in this figure (and other SPM figures below) negative GLM-weights are seen in the CSF, in white matter, around the ventricles, in large veins bordering CSF and in brain edges. These negative GLM-weights are thought to result from non-neuronal sources as was previously suggested^{40} and are ignored.

Figure 2 demonstrates the differences between β_{1} GLM-SPM and β_{4} GLM-SPM: while the functional connectivity of the left thalamus with voxels within the bilateral thalamus was manifested by all frequencies (i.e., significant for all betas), the left thalamus→occipital, cingulate, temporal and sensorimotor cortex functional connections were characterized only by the highest frequency (i.e., significant only for β_{4}). To explore these differences, we calculated the functional connectivity (cross-correlation with time lag) between the left thalamus seed and two preselected regions. These regions were selected based on the GLM-SPMs according to their significant connectivity with the seed. Specifically, we calculated the following: (i) average cross-correlation functions, (ii) average GLM-weights (betas and gammas) and (iii) average Pearson’s correlation coefficients. All these were done across all voxels in the chosen regions, for each subject separately and presented as the group mean±standard error (N=34). The fitted GLM functions for the group mean cross-correlation functions (using Equation 6) are also shown. Figure 3 shows these calculations for the functional connectivity between the left thalamus and two regions in the occipital-temporal cortex and inferior frontal gyrus, indicated by white circles in the figure. Figure 3A shows the cross-correlation function between the left thalamus and a region in the occipital-temporal cortex (centered at MNI=42, −74, 16). The cross-correlation function exhibits a significant symmetric high frequency component in addition to a tendency for significance of a symmetric and an antisymmetric low frequency components. Note, that in this case β_{4} is significant and β_{1} and γ_{1} show a tendency for significance while the Pearson’s correlation coefficient is not significant. Figure 3B shows the cross-correlation function between the left thalamus and a region in the inferior frontal gyrus (centered at MNI=−50, 14, 16). This functional connection was significant both by Pearson’s correlation coefficient and by the FPA. As seen, all γ (antisymmetric) GLM-weights are close to zero and Pearson’s correlation coefficient can be approximated by the sum of the β (symmetric) values, in line with Equation 5. Note the good agreement between the cross-correlation functions and their GLM-fits.

Figure 4 and Supplementary Information Figure 3 show GLM-SPMs of the left supplementary motor area (SMA) seed as an example for cortical GLM-SPMs. Average (across all voxels and all subjects) F-value was 16.01±0.009 (mean±SE). For this cortical seed, most of significant connections correspond to symmetric cross-correlation functions with positive weights. Antisymmetric connections are mainly with the cerebellum and at low frequencies. In the GLM-SPMs of the left SMA seed, SMA functional connectivity is mainly characterized by the lowest (β_{1,}
Fig. 4A) and the highest (β_{4}, Fig. 4B) frequencies which both show strong connectivity with the sensorimotor cortex. However, some differences are observed between β_{1} and β_{4} GLM-SPMs. For example, SMA-occipital and SMA-thalamic functional connectivity were significant only for the highest frequency (β_{4}, Fig. 4B). Note that most of SMA-cerebellum connectivity is antisymmetric and at the lowest frequency (Supplementary Information Figure 3). Figure 5 shows the average (across voxels and subjects) of the GLM-weights, Pearson’s correlation coefficient, cross-correlation functions and GLM fits for the functional connectivity of the SMA with three different regions: a precentral region (centered at MNI=48, −16, 36) (Fig. 5A), a cerebellar posterior region (centered at MNI=32, −76, −24) (Fig. 5B) and a region in the precuneus (centered at MNI=22, −74, 36) (Fig. 5C). These regions were selected from the GLM-SPMs based on their significant β_{1}, γ_{1} and β_{4} weights respectively. For the SMA-precentral connectivity (Fig. 5A), β_{1}, β_{4} and Pearson’s correlation coefficient values were significant. For the SMA-cerebellar connectivity (Fig. 5B), only the antisymmetric γ_{1} was significant. Note that the SMA-cerebellar connectivity corresponds to a time delay of ~12.5sec (π/2 of a cycle with a frequency 0.02Hz) in agreement with the maximum value of cross-correlation that is ~10sec. For the SMA-precuneus connectivity (Fig. 5C), β_{4} and Pearson’s correlation coefficient values were significant. These three connectivity patterns represent the three different types of functional connections whose characterization is possible by the proposed analysis. Note again the good agreement between the cross-correlation functions and their GLM-fits, indicating the accuracy of the beta and gamma coefficients.

Figure 6 and Supplementary Information Figure 4 show the GLM-SPMs of the left cerebellum crus 1 seed. The average (across all voxels and all subjects) F-value was 16.6±0.008 (mean±SE). For this seed, almost the entire functional connectivity was associated with the lowest frequency. In this low frequency range, both symmetric and antisymmetric GLM-SPMs had large, mostly not overlapping, connectivity volumes. The symmetric functional connectivity of the cerebellum crus 1 (β_{1} GLM-SPMs, Fig. 6A) was mainly with motor cortical areas (precentral gyrus) with negative weights, that corresponds to a time delay of ~25sec since it is associated with a phase of π in a cycle of 0.02Hz. The antisymmetric functional connectivity of the cerebellum crus 1 (γ_{1} GLM-maps, Fig. 6B) was with frontal, precentral, cingulate, SMA and occipital volumes with positive weights, that corresponds to a time delay of ~12.5sec. Interestingly, almost all of the cerebellum crus 1 functional connectivity has a non-zero phase suggesting a significant time delay between cerebellar and cerebral activity. Figure 7 shows the average (across voxels and subjects) of GLM-weights, Pearson’s correlation coefficient, cross-correlation functions and their GLM-fits for the functional connectivity of the cerebellum crus 1 with two adjacent cortical regions: a precentral region (MNI=38, −24, 56) and a frontal region (MNI=24, −2, 56). For the cerebellar-precentral connectivity (Fig. 7A), β_{1} GLM-weight was significant and also the Pearson’s correlation coefficient, as expected since all antisymmetric weights were about zero. The cerebellar-frontal functional connection (Fig. 7B) demonstrates the unique type of antisymmetric connectivity with significant γ_{1} GLM-weight. Note that the sign of symmetric functional connections (such as SMA↔cerebellum) is independent of the order of calculation, while for antisymmetric functional connections the sign of the GLM-weights depends on the order of calculation. For example, connectivity of SMA→cerebellum (SMA is the seed) is negative while connectivity of cerebellum→SMA (cerebellum is the seed) is positive. This is due to the antisymmetric cross-correlation functions.

To demonstrate that the phase of a functional connection can be different from 0°, 90° or 180°, we chose post-hoc the caudate which demonstrates this type of functional connectivity, as a seed region. Figure 8 and Supplementary Information Figure 5 show the seed-SPMs of the left caudate. The average (across all voxels and all subjects) F-value was 10.03±0.018 (mean±SE). As shown, most of the caudate’s significant connections are at the lowest frequency with positive and negative symmetric GLM-weights in addition to many antisymmetric GLM-weights. Careful examination of Fig. 8 shows that at the lowest frequency certain voxels have both significant symmetric and antisymmetric GLM-weights. To test this further, we searched for voxels within the precuneus (MNI=6, −54, 56) and paracentral lobule (MNI=6, −40, 56) which were significant for both β_{1} and γ_{1} (symmetric and antisymmetric GLM-weights for the lowest frequency). Figure 9 shows the average±SE (across voxels and subjects) of the GLM-weights, Pearson’s correlation coefficient, cross-correlation function and its GLM-fit for the connectivity of the caudate with these voxels. As seen, both β_{1} and γ_{1} are significant and the cross-correlation function shows a phase difference between the BOLD signals that according to Equation 8 equals to |34°| which corresponds to a time delay of 4.7sec. This demonstrates that the proposed analysis can detect any possible phase-difference between a pair of BOLD functions.

In this paper we introduce a new resting-state functional connectivity analysis method which we name ‘Frequency-Phase Analysis’ (FPA). This analysis enables to characterize the functional connectivity between coupled time-series functions with respect to their frequencies and phases. It presents a simple unified framework for simultaneous assessment of the frequency and phase of functional connectivity under the assumption that the time-series can be expressed by Fourier series. Pairwise cross-correlation with time-lags of the time-series functions enables to infer coherent frequencies and phase-differences between the time-series. The mathematical framework allows to distinguish between symmetric, antisymmetric and asymmetric functional connectivity types. The general linear model (GLM) was used to obtain the weights of symmetric (cosine) and antisymmetric (sine) basis-set functions. These weights were used to calculate the phase-difference between time-series at specific frequencies (Equation 8). The frequency range and resolution of the cosine and sine basis-set functions was selected to cover the entire frequency spectrum and in accordance with the discrete Fourier transform theory. The excellent agreement between the GLM results and cross-correlation functions (the high F-values of the GLM-fits) suggest that the choice of the basis-set is adequate.

The Frequency-Phase Analysis is closely related to previous studies that use time-lags^{41}^{,42}^{,43}^{,44}^{,45}. It was shown that the time-lag propagates in space within conventionally known resting-state networks, depends on neuronal state (e.g., eyes closed or open) and can be used to infer directionality of neural information flow^{43}^{,44}. On the other hand, clinical studies have identified vascular time-lags in patients suffering from hypoperfusion or ischemia^{46}^{,47}^{,48}, suggesting a dependency of the time-lag on cerebrovascular perfusion. It therefore emerges that temporal delays can reveal time-lags of neural origin or reflect regional differences in vascular dynamics. Indeed, it was suggested that extracting neural propagating information from fMRI data is possible only when the variability of the hemodynamic response function is less than the time scale of information flow^{49}. Efforts were made to characterize the differences between perfusion and neuronal contributions to the BOLD signal^{45} and to remove any confounds that can bias the time-lag analysis^{50}^{,51}. It is generally however accepted that for healthy subjects time-lags are related to neuronal transfer times with some neuro-vascular mapping, not entirely known. We claim that focusing on the time-lags could bias data interpretation and that instead the phases and phase-delays provide better measures of functional connectivity and information transfer. Since neuronal communication can be transmitted at certain frequencies and these neuronal frequencies are encoded in the fMRI signal^{23}^{,31}^{,32}, times of information transfer do not necessarily coincide with the volume or efficiency of information transfer. Consequently, equal time-lags could correspond to different volumes of transferred information when they are transferred at different frequencies. For example, it was hypothesized that cortex-hippocampus functional connectivity during memory consolidation occurs at two different frequencies^{44}. Consequently, a certain time-lag in this reciprocal process may correspond to different volumes of transferred information depending on the frequency in which it occurred. Thus, this suggests that phases or phase-differences provide better measures for information flow or for transfer efficiency.

Our analysis is also related to the cross-wavelet transform^{38}, wavelet transform coherence^{39} and dynamic coherence analysis of Yaesoubi *et al*.^{28}, with the first two using coherences module and the latter calculating phases as well. However, the aims and implications of Yaesoubi *et al*. and our proposed analysis are different. Yaesoubi *et al*. focuses on time-frequency analysis for dynamic evaluation of resting-state data, whereas here we focus on frequency-phase analysis to characterize functional connections by their phase differences. We also aim to introduce an analysis method that is similar to conventional correlation analyses. For this reason, the GLM approach was selected as it is a common functional MRI analysis approach. Besides its simplicity, the GLM approach used here reduces the number of parameters to 8 real numbers (the weights of the basis-set), which were shown to be sufficient to obtain good fits for the cross-correlation functions. In addition, since we aim to characterize functional connections by their type of phase differences (symmetric, antisymmetric and asymmetric) and not to define their phases, t-statistics can be used on each of the 8 GLM weights separately, with no need for complex-numbers statistics (e.g., circular Gaussian distributions or Gaussian hypergeometric distributions) which are more difficult and less familiar.

The ability of the FPA method to characterize functional connectivity by its phase-differences is demonstrated on resting-state fMRI data. We show that the FPA method identifies unique functional connections that were not defined before. Several examples are given for this attribute. Figure 3A shows functional connectivity of the left thalamus with a region in the occipital-temporal cortex, demonstrating an asymmetric type of functional connectivity that is not observed by the Pearson’s correlation analysis. This functional connection has a significant weight of a cosine function at 0.08Hz along with close to significant weights of cosine and sine functions at 0.02Hz. Their combined effect results with a phase-difference of ~79° corresponding to a time delay of 10.9sec. Since antisymmetric sine functions shift the cross-correlation function in time, the correlation value at time-lag zero (i.e., Pearson’s correlation) is reduced and insignificant. Figures 5B and and7B7B demonstrate examples in which Pearson’s correlation coefficients were not significant but the antisymmetric GLM-weights were significant. These connections were associated with cerebellar functional connectivity. An example for the ability of the FPA method to characterize connections by their dominant frequencies is given in Fig. 5A and C. These figures show functional connections which are significant by Pearson’s correlation analysis and by the FPA. However, the FPA method enables a more subtle characterization of these connections. Whereas the SMA-precentral functional connectivity (Fig. 5A) is mainly apparent in the lowest frequency, the SMA-occipital functional connectivity (Fig. 5C) is mainly apparent in the highest frequency. Figure 6 shows the GLM-SPMs of the left cerebellum crus 1 seed and demonstrates that antisymmetric functional connectivity can be associated with large gray matter volumes. Antisymmetric and symmetric functional connections in this figure correspond to anatomically different brain structures: the cerebellum crus 1 is functionally connected with the motor system by negative symmetric cross-correlation functions, and connected with frontal regions by antisymmetric cross-correlation functions.

Our interpretation of the biological meaning of fMRI time-series frequencies and phases is based on the assumption of a coupling between neuronal activity and BOLD signal and that the communication between different regions can be differently weighted by different frequencies^{52} and phases. Two principle mechanisms, applicable for healthy subjects, are presented below to explain the biological meaning of the measured phase differences in fMRI time-series. One is that phase-differences correspond to neuronal delays and reflect time of information transfer. This assumption was used by others to explain directionality in functional and effective connectivity^{1} and was used by us to show directionality of information flow and temporal organization using nonlinear coherence in wavelet space^{24}. The second mechanism is related to the hemodynamic effect as we previously suggested^{53}. We proposed that the hemodynamic responses of BOLD signals at different brain locations can be modulated by different weights. Specifically, if two remote regions with synchronous neuronal activity have different hemodynamic response functions, such that one is mainly affected by regional cerebral blood flow (rCBF) while the other is mainly affected by regional cerebral blood volume (rCBV), a time-delay between the two BOLD signals will be generated regardless of their neuronal coupling, since changes in rCBV are known to be delayed compared to changes in rCBF^{54}^{,55}^{,56}. Consequently, one of the BOLD signals can, for example, be described by a cosine function with a zero phase, while the other by a cosine function at the same frequency with a non-zero phase. Such scenario may also occur when two regions are affected by different neurotransmitters or by different concentrations of neurotransmitters, since neurotransmitters are known to affect blood circuitry^{53}.

We introduce a new analysis method which characterizes functional connectivity between coupled time-series functions based on their frequencies and phase differences. We speculate that these phase differences in fMRI signals correspond to time of information transfer and/or to regional dependent hemodynamic effects with are modulated by neurotransmitters concentrations. Clearly, much more studies are needed to better understand the capabilities and physiological correlates of the new proposed analysis.

The study was approved by the Ethics Committee of the General University Hospital in Prague, Czech Republic. All subjects provided written informed consent for participation in the study and all methods were performed in accordance with the relevant guidelines and regulations. Magnetic resonance images were acquired with a 3T MR scanner (Magnetom Skyra, Siemens, Germany). Each participant underwent 10-minute resting-state functional MRI during which they were instructed to fixate on a visual crosshair, remain still and awake. Wakefulness was monitored during the whole scan using an MRI compatible camera. Functional images were acquired using T2*-weighted gradient EPI sequence with TR/TE 2sec/30ms, 300 repetitions and voxel size of 3×3×3mm. 34 healthy subjects (age 64.9±8.2, 18 men and 16 women) were included in the study. High resolution anatomical images were acquired using a sagittal T1-weighted magnetization-prepared rapid acquisition gradient echo (MP-RAGE) sequence for coregistration and normalization of the functional images to MNI space.

Resting-state fMRI data was first preprocessed using Statistical Parametric Mapping (SPM8, Welcome Trust Centre for Neuroimaging, London, United Kingdom, http://www.fil.ion.ucl.ac.uk/spm/software/spm8). Standard preprocessing steps included: realignment, coregistration, normalization to MNI-space and spatial smoothing with an 8mm Gaussian kernel. Voxels of 2×2×2mm^{3} were obtained after these steps. An a-priori standard inclusion criterion of maximal head motion <3mm or 3° rotations was chosen although the average displacement for subjects was <1.5mm and <1.5° rotation. Functional connectivity analysis of resting-state data was done using CONN toolbox^{57}. Further preprocessing steps included: removal of confounds by regression (motion parameters, first principal components of CSF and white matter signals), linear detrending and band-pass filtering (0.01–0.1Hz). Removal of the first principle components of the CSF and white matter signals by regression was used rather than removal of global signal in order to minimize biases that might be introduced by global signal regression^{40}^{,58}^{,59}^{,60}^{,61}. Scrubbing was further done using DPARSF toolbox^{62}. Regions of interest (ROIs) for seed analysis were selected using the Automated Anatomical Labeling (AAL)^{63}. Although some AAL regions can be further sub-divided to more specific functional regions (e.g. the thalamus), standard AAL regions were chosen to demonstrate different functional connectivity between cortical, subcortical and cerebellar regions.

All further calculations were performed with IDL version 6.1 using custom-developed software. Cross-correlation functions were generated using 41 time-lag points which correspond to a range of −40 to 40seconds. This time range was chosen to include all relevant basis-set frequencies (0.02–0.08Hz) and to cover the relevant length of the post-stimulus response expected to modulate the hemodynamic response function of the BOLD signal by regional cerebral blood volume (rCBV)^{64}. To study possible effects of lower frequencies (<0.02Hz), a larger range of time-lag points should be selected. Beta and gamma values, i.e., the coefficients of the cosine and sine basis-set functions of Equation 6, were calculated using standard IDL routine for linear regression. Consequently, the FPA method generated 8×N scalars for each seed: 4 beta and 4 gamma values for each of the N functional connections.

The level of significance in the FPA GLM-SPM calculations was defined by one sample t-test on each GLM weight separately. Significance was defined by a combination of a voxel-level p-value of p<0.0001 and a cluster-level threshold of 100 voxels. The cluster size was chosen according to Monte-Carlo simulation to set the threshold at p<0.01 corrected for multiple comparisons^{65}. For the ROI statistics (Figures 3, ,5,5, ,77 and and9),9), one sample t-test on each GLM weight was done separately using a threshold of p<0.0001.

Seed-voxel statistical-parametric-maps (SPM) were obtained from 3 predefined ROIs that represent different types of brain regions (cortex, subcortex and cerebellum) and are all part of the motor system. These were: the left supplementary motor area (SMA), left thalamus and left cerebellum crus 1. The motor system was chosen arbitrarily as an example to demonstrate the capabilities of the proposed analysis. At a later stage, we searched for clusters in the GLM-SPM that showed asymmetric functional connectivity. The left caudate was identified as showing asymmetric functional connections and therefore chosen as a 4^{th} seed region. Note that these seeds were not chosen by their goodness of fit to Equation 6.

**How to cite this article:** Goelman, G. *et al*. Frequency-phase analysis of resting-state functional MRI. *Sci. Rep.*
**7**, 43743; doi: 10.1038/srep43743 (2017).

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This work was supported by a Czech-Israel cooperative research grant issued by the Israel Ministry of science (grant number 3–9814) and the Czech Ministry of education (Grant LH13256 VES13-KontaktII). Additional support was given by the Czech Ministry of Health (Grant IGA MZ CR NT12282-5/2011) and the Charles University in Prague (Research Project PRVOUK-P26/LF1/4).

The authors declare no competing financial interests.

**Author Contributions** G.G.: Conception and design of study, analysis and writing of the paper R.D.: Data preprocessing and assistance in writing of the paper F.R., O.B., E.R., J.R., J.V., R.J.: Data acquisition.

- Friston K. J. Functional and effective connectivity: a review. Brain Connect 1, 13–36, doi: (2011).10.1089/brain.2011.0008 [PubMed] [Cross Ref]
- Friston K. J., Harrison L. & Penny W. Dynamic causal modelling. NeuroImage 19, 1273–1302 (2003). [PubMed]
- McIntosh. A. R. & Gonzalez-Lima F. Structural equation modeling and its application to network analysis in functional brain imaging. Human brain mapping 2, 2–22 (1994).
- Buckner R. L., Andrews-Hanna J. R. & Schacter D. L. The brain’s default network: anatomy, function, and relevance to disease. Annals of the New York Academy of Sciences 1124, 1–38, doi: (2008).10.1196/annals.1440.011 [PubMed] [Cross Ref]
- Fox M. D. & Raichle M. E. Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nat Rev Neurosci 8, 700–711 (2007). [PubMed]
- Biswal B., Yetkin F. Z., Haughton V. M. & Hyde J. S.
Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn Reson Med
34
**(4)**, 537–541 (1995). [PubMed] - Fox M. D. et al. . The human brain is intrinsically organized into dynamic, anticorrelated functional networks. Proceedings of the National Academy of Sciences of the United States of America 102, 9673–9678 (2005). [PubMed]
- Greicius M. D., Krasnow B., Reiss A. L. & Menon V. Functional connectivity in the resting brain: a network analysis of the default mode hypothesis. Proceedings of the National Academy of Sciences of the United States of America 100, 253–258 (2003). [PubMed]
- Raichle M. E. et al. . A default mode of brain function. Proceedings of the National Academy of Sciences of the United States of America 98, 676–682 (2001). [PubMed]
- Deco G., Jirsa V. K. & McIntosh A. R. Emerging concepts for the dynamical organization of resting-state activity in the brain. Nat Rev Neurosci 12, 43–56 (2011). [PubMed]
- Shehzad Z. et al. . The resting brain: unconstrained yet reliable. Cerebral cortex 19, 2209–2229, doi: (2009).10.1093/cercor/bhn256 [PMC free article] [PubMed] [Cross Ref]
- Zuo X. N. et al. . Reliable intrinsic connectivity networks: test-retest evaluation using ICA and dual regression approach. NeuroImage 49, 2163–2177, doi: (2010).10.1016/j.neuroimage.2009.10.080 [PMC free article] [PubMed] [Cross Ref]
- Hampson M., Driesen N. R., Skudlarski P., Gore J. C. & Constable R. T. Brain connectivity related to working memory performance. J Neurosci 26, 13338–13343, doi: (2006).10.1523/JNEUROSCI.3408-06.2006 [PMC free article] [PubMed] [Cross Ref]
- Seeley W. W., Crawford R. K., Zhou J., Miller B. L. & Greicius M. D. Neurodegenerative diseases target large-scale human brain networks. Neuron 62, 42–52 (2009). [PMC free article] [PubMed]
- Honey C. J., Kotter R., Breakspear M. & Sporns O. Network structure of cerebral cortex shapes functional connectivity on multiple time scales. Proceedings of the National Academy of Sciences of the United States of America 104, 10240–10245 (2007). [PubMed]
- Ghosh A., Rho Y., McIntosh A. R., Kotter R. & Jirsa V. K. Noise during rest enables the exploration of the brain’s dynamic repertoire. PLoS Comput Biol 4, e1000196 (2008). [PMC free article] [PubMed]
- Deco G., Jirsa V., McIntosh A. R., Sporns O. & Kotter R. Key role of coupling, delay, and noise in resting brain fluctuations. Proceedings of the National Academy of Sciences of the United States of America 106, 10302–10307 (2009). [PubMed]
- Cordes D. et al. . Frequencies Contributing to Functional Connectivity in the Cerebral Cortex in “Resting-state” Data. Am J Neuroradiol 22, 1326–1333 (2001). [PubMed]
- Sun F. T., Miller L. M. & D’Esposito M. Measuring interregional functional connectivity using coherence and partial coherence analyses of fMRI data. NeuroImage 21, 647–658, doi: (2004).10.1016/j.neuroimage.2003.09.056 [PubMed] [Cross Ref]
- Salvador R., Suckling J., Schwarzbauer C. & Bullmore E. Undirected graphs of frequency-dependent functional connectivity in whole brain networks. Philos Trans R Soc Lond B Biol Sci 360, 937–946, doi: (2005).10.1098/rstb.2005.1645 [PMC free article] [PubMed] [Cross Ref]
- Thirion B., Dodel S. & Poline J. B. Detection of signal synchronizations in resting-state fMRI datasets. NeuroImage 29, 321–327, doi: (2006).10.1016/j.neuroimage.2005.06.054 [PubMed] [Cross Ref]
- Salvador R. et al. . Frequency based mutual information measures between clusters of brain regions in functional magnetic resonance imaging. NeuroImage 35, 83–88, doi: (2007).10.1016/j.neuroimage.2006.12.001 [PubMed] [Cross Ref]
- Salvador R. et al. . A simple view of the brain through a frequency-specific functional connectivity measure. NeuroImage 39, 279–289, doi: (2008).10.1016/j.neuroimage.2007.08.018 [PubMed] [Cross Ref]
- Goelman G. & Dan R. Multiple-region directed functional connectivity based on phase delays. Human brain mapping, doi: (2016).10.1002/hbm.23460 [PubMed] [Cross Ref]
- Handwerker D. A., Roopchansingh V., Gonzalez-Castillo J. & Bandettini P. A. Periodic changes in fMRI connectivity. NeuroImage 63, 1712–1719, doi: (2012).10.1016/j.neuroimage.2012.06.078 [PMC free article] [PubMed] [Cross Ref]
- Hutchison R. M., Womelsdorf T., Gati J. S., Everling S. & Menon R. S. Resting-state networks show dynamic functional connectivity in awake humans and anesthetized macaques. Human brain mapping 34, 2154–2177, doi: (2013).10.1002/hbm.22058 [PubMed] [Cross Ref]
- Kiviniemi V. et al. . A sliding time-window ICA reveals spatial variability of the default mode network in time. Brain Connect 1, 339–347, doi: (2011).10.1089/brain.2011.0036 [PubMed] [Cross Ref]
- Yaesoubi M., Allen E. A., Miller R. L. & Calhoun V. D. Dynamic coherence analysis of resting fMRI data to jointly capture state-based phase, frequency, and time-domain information. NeuroImage 120, 133–142, doi: (2015).10.1016/j.neuroimage.2015.07.002 [PMC free article] [PubMed] [Cross Ref]
- Glerean E., Salmi J., Lahnakoski J. M., Jaaskelainen I. P. & Sams M. Functional magnetic resonance imaging phase synchronization as a measure of dynamic functional connectivity. Brain Connect 2, 91–101, doi: (2012).10.1089/brain.2011.0068 [PMC free article] [PubMed] [Cross Ref]
- Liu X. & Duyn J. H. Time-varying functional network information extracted from brief instances of spontaneous brain activity. Proceedings of the National Academy of Sciences of the United States of America 110, 4392–4397, doi: (2013).10.1073/pnas.1216856110 [PubMed] [Cross Ref]
- Baria A. T., Baliki M. N., Parrish T. & Apkarian A. V. Anatomical and functional assemblies of brain BOLD oscillations. J Neurosci 31, 7910–7919, doi: (2011).10.1523/JNEUROSCI.1296-11.2011 [PMC free article] [PubMed] [Cross Ref]
- Esposito F. et al. . Rhythm-specific modulation of the sensorimotor network in drug-naive patients with Parkinson’s disease by levodopa. Brain 136, 710–725, doi: (2013).10.1093/brain/awt007 [PubMed] [Cross Ref]
- Devergnas A., Pittard D., Bliwise D. & Wichmann T. Relationship between oscillatory activity in the cortico-basal ganglia network and parkinsonism in MPTP-treated monkeys. Neurobiol Dis 68, 156–166, doi: (2014).10.1016/j.nbd.2014.04.004 [PMC free article] [PubMed] [Cross Ref]
- Thompson W. H. & Fransson P. The frequency dimension of fMRI dynamic connectivity: Network connectivity, functional hubs and integration in the resting brain. NeuroImage 121, 227–242, doi: (2015).10.1016/j.neuroimage.2015.07.022 [PubMed] [Cross Ref]
- Majeed W. et al. . Spatiotemporal dynamics of low frequency BOLD fluctuations in rats and humans. NeuroImage 54, 1140–1150, doi: (2011).10.1016/j.neuroimage.2010.08.030 [PMC free article] [PubMed] [Cross Ref]
- Thompson G. J., Pan W. J., Magnuson M. E., Jaeger D. & Keilholz S. D. Quasi-periodic patterns (QPP): large-scale dynamics in resting state fMRI that correlate with local infraslow electrical activity. NeuroImage 84, 1018–1031, doi: (2014).10.1016/j.neuroimage.2013.09.029 [PMC free article] [PubMed] [Cross Ref]
- Pan W. J. et al. . Broadband local field potentials correlate with spontaneous fluctuations in functional magnetic resonance imaging signals in the rat somatosensory cortex under isoflurane anesthesia. Brain Connect 1, 119–131, doi: (2011).10.1089/brain.2011.0014 [PMC free article] [PubMed] [Cross Ref]
- Chang C. & Glover G. H. Time-frequency dynamics of resting-state brain connectivity measured with fMRI. NeuroImage 50, 81–98, doi: (2010).10.1016/j.neuroimage.2009.12.011 [PMC free article] [PubMed] [Cross Ref]
- Torrence C. & Webster P. Interdecadal changes in the ENSO-Monsoon system. J Clim 12, 2679–2690 (1999).
- Bianciardi M., Fukunaga M., van Gelderen P., de Zwart J. A. & Duyn J. H. Negative BOLD-fMRI signals in large cerebral veins. J Cereb Blood Flow Metab 31, 401–412 (2011). [PMC free article] [PubMed]
- Mitra A., Snyder A. Z., Hacker C. D. & Raichle M. E. Lag structure in resting-state fMRI. J Neurophysiol 111, 2374–2391, doi: (2014).10.1152/jn.00804.2013 [PubMed] [Cross Ref]
- Mitra A., Snyder A. Z., Constantino J. N. & Raichle M. E. The Lag Structure of Intrinsic Activity is Focally Altered in High Functioning Adults with Autism. Cerebral cortex, doi: (2015).10.1093/cercor/bhv294 [PubMed] [Cross Ref]
- Mitra A., Snyder A. Z., Blazey T. & Raichle M. E. Lag threads organize the brain’s intrinsic activity. Proceedings of the National Academy of Sciences of the United States of America 112, E2235–2244, doi: (2015).10.1073/pnas.1503960112 [PubMed] [Cross Ref]
- Mitra A. et al. . Human cortical-hippocampal dialogue in wake and slow-wave sleep. Proceedings of the National Academy of Sciences of the United States of America 113, E6868–E6876, doi: (2016).10.1073/pnas.1607289113 [PubMed] [Cross Ref]
- Amemiya S., Takao H., Hanaoka S. & Ohtomo K. Global and structured waves of rs-fMRI signal identified as putative propagation of spontaneous neural activity. NeuroImage 133, 331–340, doi: (2016).10.1016/j.neuroimage.2016.03.033 [PubMed] [Cross Ref]
- Lv Y. et al. . Identifying the perfusion deficit in acute stroke with resting-state functional magnetic resonance imaging. Annals of neurology 73, 136–140, doi: (2013).10.1002/ana.23763 [PubMed] [Cross Ref]
- Amemiya S., Kunimatsu A., Saito N. & Ohtomo K. Cerebral hemodynamic impairment: assessment with resting-state functional MR imaging. Radiology 270, 548–555, doi: (2014).10.1148/radiol.13130982 [PubMed] [Cross Ref]
- Christen T. et al. . Noncontrast mapping of arterial delay and functional connectivity using resting-state functional MRI: a study in Moyamoya patients. Journal of magnetic resonance imaging: JMRI 41, 424–430, doi: (2015).10.1002/jmri.24558 [PMC free article] [PubMed] [Cross Ref]
- Garg R., Cecchi G. A. & Rao A. R. Full-brain auto-regressive modeling (FARM) using fMRI. NeuroImage 58, 416–441, doi: (2011).10.1016/j.neuroimage.2011.02.074 [PubMed] [Cross Ref]
- Power J. D., Barnes K. A., Snyder A. Z., Schlaggar B. L. & Petersen S. E. Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion. NeuroImage 59, 2142–2154, doi: (2012).10.1016/j.neuroimage.2011.10.018 [PMC free article] [PubMed] [Cross Ref]
- Power J. D., Plitt M., Laumann T. O. & Martin A. Sources and implications of whole-brain fMRI signals in humans. NeuroImage, doi: (2016).10.1016/j.neuroimage.2016.09.038 [PubMed] [Cross Ref]
- Zuo X. N. et al. . The oscillating brain: complex and reliable. NeuroImage 49, 1432–1445, doi: (2010).10.1016/j.neuroimage.2009.09.037 [PMC free article] [PubMed] [Cross Ref]
- Goelman G., Gordon N. & Bonne O. Maximizing negative correlations in resting-state functional connectivity MRI by time-lag. PLoS One 9, e111554, doi: (2014).10.1371/journal.pone.0111554 [PMC free article] [PubMed] [Cross Ref]
- Buxton R. B., Wong E. C. & Frank L. R. Dynamics of blood flow and oxygenation changes during brain activation: the balloon model. Magn Reson Med 39, 855–864 (1998). [PubMed]
- Buxton R. B., Uludag K., Dubowitz D. J. & Liu T. T. Modeling the hemodynamic response to brain activation. NeuroImage 23 Suppl 1, S220–233 (2004). [PubMed]
- Zong X., Kim T. & Kim S. G. Contributions of dynamic venous blood volume versus oxygenation level changes to BOLD fMRI. NeuroImage 60, 2238–2246, doi: (2012).10.1016/j.neuroimage.2012.02.052 [PMC free article] [PubMed] [Cross Ref]
- Whitfield-Gabrieli S. & Nieto-Castanon A. Conn: a functional connectivity toolbox for correlated and anticorrelated brain networks. Brain Connect 2, 125–141, doi: (2012).10.1089/brain.2012.0073 [PubMed] [Cross Ref]
- Fox M. D., Zhang D., Snyder A. Z. & Raichle M. E. The global signal and observed anticorrelated resting state brain networks. Journal of neurophysiology 101, 3270–3283 (2009). [PubMed]
- Murphy K., Birn R. M., Handwerker D. A., Jones T. B. & Bandettini P. A. The impact of global signal regression on resting state correlations: are anti-correlated networks introduced? NeuroImage 44, 893–905 (2009). [PMC free article] [PubMed]
- Weissenbacher A. et al. . Correlations and anticorrelations in resting-state functional connectivity MRI: a quantitative comparison of preprocessing strategies. NeuroImage 47, 1408–1416 (2009). [PubMed]
- Chai X. J., Castanon A. N., Ongur D. & Whitfield-Gabrieli S. Anticorrelations in resting state networks without global signal regression. NeuroImage 59, 1420–1428 (2011). [PMC free article] [PubMed]
- Chao-Gan Y. & Yu-Feng Z. DPARSF: A MATLAB Toolbox for “Pipeline” Data Analysis of Resting-State fMRI. Front Syst Neurosci 4, 13, doi: (2010).10.3389/fnsys.2010.00013 [PMC free article] [PubMed] [Cross Ref]
- Tzourio-Mazoyer N. et al. . Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. NeuroImage 15, 273–289, doi: (2002).10.1006/nimg.2001.0978 [PubMed] [Cross Ref]
- Arichi T. et al. . Development of BOLD signal hemodynamic responses in the human brain. NeuroImage 63, 663–673 (2012). [PMC free article] [PubMed]
- Forman S. D. et al. . Improved assessment of significant activation in functional magnetic resonance imaging (fMRI): use of a cluster-size threshold. Magn Reson Med 33, 636–647 (1995). [PubMed]

Articles from Scientific Reports are provided here courtesy of **Nature Publishing Group**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |