Evaluating phase synchronization methods in fMRI: A comparison study and new approaches.

In recent years there has been growing interest in measuring time-varying functional connectivity between different brain regions using resting-state functional magnetic resonance imaging (rs-fMRI) data. One way to assess the relationship between signals from different brain regions is to measure their phase synchronization (PS) across time. There are several ways to perform such analyses, and we compare methods that utilize a PS metric together with a sliding window, referred to here as windowed phase synchronization (WPS), with those that directly measure the instantaneous phase synchronization (IPS). In particular, IPS has recently gained popularity as it offers single time-point resolution of time-resolved fMRI connectivity. In this paper, we discuss the underlying assumptions required for performing PS analyses and emphasize the importance of band-pass filtering the data to obtain valid results. Further, we contrast this approach with the use of Empirical Mode Decomposition (EMD) to achieve similar goals. We review various methods for evaluating PS and introduce a new approach within the IPS framework denoted the cosine of the relative phase (CRP). We contrast methods through a series of simulations and application to rs-fMRI data. Our results indicate that CRP outperforms other tested methods and overcomes issues related to undetected temporal transitions from positive to negative associations common in IPS analysis. Further, in contrast to phase coherence, CRP unfolds the distribution of PS measures, which benefits subsequent clustering of PS matrices into recurring brain states.

Introduction

It was previously assumed that functional connectivity (FC) in the brain was static during the course of a single resting-state functional magnetic resonance imaging (rs-fMRI) run. Recently, however, several studies (Allen et al., 2014; Chang and Glover, 2010; Hutchison et al., 2013a; Lurie et al., 2019; Preti et al., 2016; Tagliazucchi et al., 2012; Thompson et al., 2013) have pointed to dynamic changes in FC taking place in a considerably shorter time window than previously thought (i.e., on the order of seconds and minutes). Several methods have been proposed to investigate such time-varying connectivity (TVC). These include the widely-used sliding-window approach (Chang and Glover, 2010; Hutchison et al., 2013a; 2013b; Tagliazucchi et al., 2012), change point analysis (Cribben et al., 2012; 2013; Xu and Lindquist, 2015), point process analysis (Tagliazucchi et al., 2011), co-activation patterns (CAPs) (Liu and Duyn, 2013), transient-based CAPs (Karahanoğlu and Van De Ville, 2015), time series models (Lindquist et al., 2014), time-frequency analysis (Chang and Glover, 2010), and variants of hidden Markov models (HMMs) (Bolton et al., 2018; Eavani et al., 2013; Shappell et al., 2019; Vidaurre et al., 2017). Despite development of these promising approaches, estimation of TVC remains a challenging endeavor due to the low signal-to-noise ratio (SNR) of the blood oxygen level dependent (BOLD) signal and the presence of image artifacts and nuisance confounds (Hutchison et al., 2013a; Laumann et al., 2016; Lindquist et al., 2014).

The term synchronization refers to the coordination in the state of two or more systems that can be attributed to their interaction (or coupling) (Rosenblum et al., 1996). Recently, phase synchronization (PS) methods were proposed as a means of measuring the level of synchrony between time series from different regions of interest (ROIs) in the brain (Glerean et al., 2012; Pedersen et al., 2017; 2018). Typically, the phase of a particular time series is computed at each time point through the application of the Hilbert transform, and used to evaluate the phase difference between various time series. Two time series in synchronization will maintain a constant phase difference. In this study, we differentiate between methods that combine a PS metric with a sliding window approach, referred to as windowed phase synchronization (WPS), with those that directly measure PS at each time point, referred to as instantaneous phase synchronization (IPS).

The first class of methods (WPS methods) uses metrics that provide a single omnibus measure of the phase synchronization between two time series obtained using the Hilbert transform. This approach is similar to how correlations provide an omnibus measure of the linear relationship between time series (analogous to the static correlation used in FC). In this approach, a sliding window technique is used to compute the metric locally within a specific time window. As the window is shifted across time, one can obtain a time-varying value of the measure of interest (i.e., the dynamic synchronization between two time series). The use of Phase Locking Value (PLV) (Boccaletti et al., 2018; Glerean et al., 2012; Pauen and Ivanova, 2013) to capture time-varying relationship between a pair of signals has recently been used in this context (Rebollo et al., 2018). In this paper, we propose two other measures that can capture the time-varying relationship between a pair of signals: circular-circular correlation (Jammalamadaka and Sengupta, 2001; Pauen and Ivanova, 2013), and toroidal-circular correlation (Zhan et al., 2017). Importantly, this class of methods suffers from similar issues as sliding-window correlations, such as the need to select an a priori window length for analysis.

The second class of methods (IPS methods) directly analyzes the instantaneous phases extracted using the Hilbert transform. In recent years there has been growing interest in using IPS methods in neuroimaging, with the bulk of the work applied to MEG and EEG data. However, several studies have also applied IPS methods to fMRI data. For instance, Laird et al. (2002) used IPS methods to analyze task-activated fMRI data. However, the lack of narrow band-pass filtering in the study’s analysis pipeline brings into question the validity of the results. Niazy et al. (2011) studied the spectral characteristics of resting state network (RSN) and suggested that it is important to consider the IPS between various RSNs at different frequencies. Glerean et al. (2012) proposed using IPS as a measure of TVC. Finally, Pedersen et al. (2018) examined the relationship between IPS and Correlation-based Sliding Window (CSW) techniques and observed a strong association between the two methods when using absolute values of CSW. Benefits of using an IPS approach is that it offers single time-point resolution of time-resolved fMRI connectivity, and does not require choosing an arbitrary window length.

In this paper, we discuss the concept of phase synchronization in the context of fMRI, with a particular focus on TVC. We begin by reviewing the framework for computing the phase from time series data using the Hilbert transform, and discuss the importance of band-pass filtering the data to accurately estimate the instantaneous phases. In the EEG literature, an alternative to band-pass filtering is to use empirical mode decomposition (EMD) (Huang, 2014) as a preprocessing step. Here we contrast the two approaches in the context of analyzing PS in the fMRI setting. We continue by introducing a number of different methods for evaluating phase synchronization. We focus both on methods already in common use, such as the phase locking value and phase coherence, as well as methods new to the fMRI literature, such as circular-circular correlation and toroidal-circular correlation (Zhan et al., 2017). Finally, we propose a new variant of phase coherence, denoted the Cosine of the Relative Phase (CRP), that can be used to compute IPS. We contrast these methods through a series of simulations and application to two rs-fMRI data sets. The first consists of test-retest data with a temporal resolution of 2 s, while the second is data from the Human Connectome Project with a temporal resolution of 0.72 s.

Methods

A framework for computing instantaneous phase

To obtain the instantaneous phase (Boccaletti et al., 2018) of an arbitrary real signal x(t) one must first construct an analytic signal: where and represents the Hilbert transform. This signal can in turn be re-expressed as follows: where represents the envelope and the instantaneous phase.

Here x(t) is assumed to satisfy Bedrosian’s Product Theorem, which states that a band-limited signal can be decomposed into the product of envelope and phase when their spectra are disjoint. This holds if the signal of interest is first narrow-banded by applying a band-pass filter.

There are two important considerations when choosing the appropriate filter to apply. First, it should not corrupt the phase information in the signal. Thus, it is important to use a filter that does not shift the phase. One class of filters that accomplishes this goal are zero-phase filters. Second, the width of the frequency band must be sufficiently narrow. The narrower the band, the closer the signal will be to a mono-component signal and the Hilbert transform will produce an analytic signal with meaningful envelope and phase. The choice of appropriate band widths in this context have been investigated in previous studies of fMRI data. For example, Ponce-Alvarez et al. (2015) examined band-pass filtering of fMRI data using various frequency bands in the range of 0.01 – 0.13 and reported consistent results for phase statistics at each frequency band. Pedersen et al. (2018) compared using a narrow band-pass filter (0.03 – 0.07 ) with a wider band-pass filter (0.01 – 0.1 ), and found that the narrow-band data yielded stronger associations between the results of CSW and IPS analyses.

A schematic framework for obtaining the instantaneous phase synchronization is shown in Fig. 1.

Fig. 1.

A schematic of the approach to calculate the instantaneous phase (IP) framework.

Consider that a pair of time series x(t) and , from two different ROIs are filtered using a narrow band-pass zero-phase filter, . and denote the filtered data by and respectively, i.e., Here * represents the convolution operator.

If Bedrosian’s theorem holds, the analytical signals of the narrow-banded time series can be expressed as the product of instantaneous envelope and instantaneous phase: Here the subscript a refers to analytical signal. Throughout, we assume that and are the phase time series extracted from a pair of time series x(t), and y(t). Using the instantaneous phases, synchronization can be assessed by studying their differences.

An alternative to band-pass filtering is to use EMD. Here one decomposes the time series into a sum of oscillatory modes referred to as intrinsic mode functions (IMF) that correspond to different frequency contents in the time series. EMD provides a data-driven signal decomposition and does not require an a priori defined basis system. The first IMF consists of the largest frequency oscillation present in the signal, and each subsequent IMF consists of increasingly smaller frequency oscillations than the those previously extracted. A more detailed description of bivariate EMD (BEMD), which is used in this paper, can be found in the Supplementary Materials. In the continuation, we refer to this as the BEMD-based PS framework to differentiate it from the framework described above.

Methods for assessing phase synchronization

Next, we describe how to measure PS based on the extracted phase time series. We discriminate between methods that utilize a PS metric together with a sliding window approach (WPS) from those that directly measure IPS.

Windowed phase synchronization

The first class of methods, place a measure of PS across two time series within a sliding window framework. Here we describe this approach using PLV, circular-circular correlation, and toroidal-circular correlation.

Phase Locking Value

The PLV is a classic metric for assessing phase synchronization based on quantifying to what extent the two signals are phase locked. PLV has found widespread use in the analysis of MEG/EEG data (Halliday et al., 1998; Rosenblum et al., 2000). This notion of synchronization can be expressed as follows: Here the integers m and n are the synchronization indices and the generalized phase difference time series. In this paper, we assume m = n = 1 and drop the indices and let .

Using the instantaneous phase difference of the signals at each time point, the PLV can be computed as follows: where the operator denotes averaging over time. If the pair of signals are unsynchronized, then and follows a uniform distribution; otherwise, if the pair are synchronized, is constant and equal to 1 (Lachaux et al., 1999).

To compute PLV within a sliding window framework, for each time window of length ℓ, the PLV between the pair of the signals can be obtained using Eq. (8). This approach has been previously used for assessing the episodes of elevated gastric-BOLD synchronization by Rebollo et al. (2018) in the study of stomach-brain synchrony.

Next, we introduce two other measures in WPS approach that to our best of knowledge have not been used to assessing the time varying phase synchornization in a sliding window fashion.

Circular-Circular Correlation

The instantaneous phase obtained from each time series are directional data and follow a circular distribution. In this context, the use of the standard Pearson correlation coefficient is no longer appropriate. Instead, a more suitable measure is circular-circular correlation (Jammalamadaka and Sarma, 1988; Jammalamadaka and Sengupta, 2001), defined as follows:

In the equation above, and , while μ and ν represent the mean directions of and , respectively. Thus, the terms and can be interpreted as the deviations of and from their corresponding mean directions.

The circular-circular correlation provides a measure of the static interdependence between the two phase time series. It can also be used within the sliding windows framework to investigate the time-varying PS. This can be expressed as: where and represent the estimated time-varying mean of the two phase time series over the sliding window. It is important to note in the context of directional statistics, is computed as follows (Bishop, 2006; Jammalamadaka and Sengupta, 2001):

This formulation can be understood by representing each directional variable on a unit circle (r = 1) in the polar coordinate system (r, ϕ). Re-expressing to the Cartesian coordinate system (x, y), we can write and , for . Thus, , , and the mean direction thus can be written as expressed in Eq. (11). Here is computed analogously.

Toroidal-Circular Correlation

In a critique of circular-circular correlation, Zhan et al. (2017) argued that the sine of an angle contains less information than the angle itself, as multiple angles can take the same sine value. Furthermore, since the sine function is not monotone within an interval of π, this may lead to unreasonable results. To circumvent these issues, they introduced a circular correlation coefficient for bivariate directional data on a torus, which is the equivalent to the product of two circles (Sojakova, 2016; Zhan et al., 2017).

To elaborate, let (or similarly ) and (or ) be two circular data points and set , , so . When OR , by convention the direction from to is considered to be clockwise. When OR , the direction is considered counter-clockwise. The same definition holds for and .

Let , then the order function can be expressed as follows:

Now, let us assume that and are independent. The toroidal-circular correlation is then defined as follows:

Note that can take values between −1 to +1. Based on this definition, the two variables ϕ and ϕ move on the circumference of the torus in the same direction if , making . Similarly, indicates that the two variables are moving in opposite directions.

An estimator can be obtained as follows:

The main advantage of using toroidal-circular correlation is that no information about the angles are lost. Thus, the estimator circumvents problems due to the non-monotonicity of the sine function that could result in irregular estimation when using the circular-circular correlation. It can be calculated within a sliding window framework similar to the previous sections.

Instantaneous phase synchronization

The second class of methods are based on directly working with the instantaneous phases of two time series. Here we focus on phase coherence (Pedersen et al., 2017; 2018), which has already found wide usage in the field, and the cosine of the relative phase, introduced for the first time in this work.

Phase Coherence

The phase coherence at each time point is defined as follows: Here the absolute value of the sine of the relative phase differences is included to account for phase wrapping and resolve issues with phase ambiguity over time. Note that the range of the values obtained using this metric will take values between 0 and 1, where 0 implies no phase coherence and 1 corresponds to maximal phase coherence.

A shortcoming of this approach is that it discards information about the direction of the relationship as . As these values vary between 0 and 1, Ψ(t) as defined in Eq. (15) does not capture negative association (i.e., when signals are in anti-phase). This may help explain why Pedersen et al. (2018) found that the association between IPS and CSW analysis was strongly dependent on negative correlations obtained from the CSW analysis, and that the association increased when comparing the absolute values of the correlations. In addition, it explains why their analysis was unable to capture temporal transitions from positive to negative associations, and vice versa, that appeared in the CSW analysis.

Cosine of the Relative Phase

To circumvent the issues outlined above, we propose a modification of phase coherence that takes temporal transitions into account and preserves the correlational structure in the data. This can be achieved by not taking the absolute value of the phase difference and using a cosine function instead of a sine function. We refer to this measure as the cosine of the relative phase (CRP), defined as follows: Notably, the range of the values obtained using this metric take values between −1 and 1, and is therefore directly comparable to standard correlation values.

The CRP approach avoids phase unwrapping and takes phase ambiguity into consideration. When the instantaneous phase of two signals are similar to one another (i.e., ), CRP yields a value close to 1. When the phases are dissimilar but in the same direction, their relative phase difference is bounded between , which is the range where the cosine function is positive. As the phases become orthogonal to one another, CRP approaches 0 indicating a lack of coherence. Similarly, the CRP captures negative associations between phases. If the phase difference is greater than ±π/2, this results in negative values of the cosine function. Thus, using CRP as a measure of phase synchrony helps overcome the issue of detecting temporal transitions from positive to negative associations (or vice-versa), and preserves the positive and negative dependence in the data.

Simulations

In this section we introduce three simulations designed to compare the methods presented in Section 2.2 for the two different classes of PS analysis. The first investigates their performance in a null setting, while the second and third investigate PS measures when two sinusoidal signals have the same frequency but differing phase shifts. For all three simulations data was generated with a sampling frequency of 1/T R, where T R represents the repetition time of an fMRI experiment. To be comparable with the main rs-fMRI data set used in this paper we chose T R = 2 seconds.

For each simulation we computed WPS values using the PLV, circular-circular correlation, and toroidal-circular correlation, and IPS values using phase coherence and the CRP method. All simulations were repeated 1000 times, and the mean and variance of the PS measured at each time point was used to construct a 95% confidence interval. Furthermore, the effect of different window lengths in the WPS analysis was evaluated using three different window lengths (30, 60, and 120 TRs).

To illustrate the necessity of band-pass filtering the data, PS analysis was performed on the simulated data both before and after band-pass filtering it in the range [0.03, 0.07] . Throughout we used a 5th order Butterworth filter. The zero-phase version of this filter is implemented in MATLAB by filtering backward in time using MATLAB’s filtfilt function to cancel out the phase delay introduced by this filter.

For comparison purposes we repeated the simulations using BEMD, described in Section S1.1, in place of the band-pass filtering. Using the BEMD-based PS framework, the pairwise intrinsic mode functions (IMF)s are obtained and PS measures are computed on the IMF pair whose frequency is closest to central frequency of 0.05 Hz which corresponds to the peak of the power spectrum.

Simulation 1:

To simulate time series with independent phase dynamics, we generated two independent random signals from a Gaussian distribution with mean 0 and standard deviation 1. Using the logic of surrogate data testing, we generated surrogate data under the assumption of no relationship between the phase from the two signals. To achieve this goal we used cyclic phase permutation (CPP) surrogates (Lancaster et al., 2018), constructed by reorganizing the cycles within the extracted phase of the signals. This destroys any phase dependence between the pair, whilst preserving the general form of the phase dynamics of each time series. For this simulation, the 1000 realizations of signal pairs were generated using CPP surrogates.

Simulation 2:

Here we generated two sinusoidal signals with the same frequency, but with a time-varying phase shift corresponding to a ramp function. To elaborate, consider two sinusoidal signals x(t) and y(t). Let x(t) be the reference signal with an angular frequency of ω0 and phase φ(t). Further, let y(t) have the same angular frequency but with phase φ(t). The signals can be expressed as follows:

Without loss of generality, let and be a ramp function, The time series can then be expressed as follows: Throughout, rad / s with f = 0.05 , and the transition is set to occur at t0 = 170 s. The noise terms ε and ε are Gaussian white noise with mean 0 and standard deviation 1.

To summarize, the two signals start out phase synchronized and remain in this state up to t0 = 170 s. After which the phase difference starts linearly increasing and transition into a non-synchronized state.

Simulation 3:

Here we generated two sinusoidal signals with the same frequency, but with a time-varying phase shift corresponding to a sigmoid function. As in the previous simulation, data was generated according to Eq. (17). Here we let φ(t) = 0 and φ(t) be a sigmoid function, i.e. Hence, the time series can be expressed as follows: Throughout, we set a = 2π, b = −0.01. t0 = 170, and rad / s with f = 0.05 . The noise terms ε and ε are Gaussian white noise with mean 0 and standard deviation 1.

To summarize, the signals are initially in phase, after which the amount of phase shift gradually increases. This continues until t = 170 when the pairs are in anti-phase synchronization. Thereafter, the signals gradually return to being in phase. The transition between the phase of the signals from 0 to 2π occurs smoothly and monotonically increasing.

Application to Kirby21 dataset

Image acquisition

We used the Multi-Modal MRI Reproducibility Resource from the F.M. Kirby Research Center, commonly referred to as Kirby21. It includes data from 21 healthy adults scanned on a 3T Philips Achieva scanner designed to achieve 80 mT/m maximum gradient strength with body coil excitation and an eight channel phased array SENSitivity Encoding (SENSE) (Pruessmann et al., 1999) head-coil for reception. All subjects completed two scanning sessions on the same day, between which they briefly exited the scan room. A T1-weighted (T1w) Magnetization-Prepared Rapid Acquisition Gradient Echo (MPRAGE) structural run was acquired during both sessions (acquisition time = 6 min, TR/TE/TI = 6.7/3.1/842 ms, resolution = 1 × 1 × 1.2 mm3, SENSE factor = 2, flip angle = 8° ). A multi-slice SENSE-EPI pulse sequence (Pruessmann et al., 1999; Stehling et al., 1991) was used to acquire two rs-fMRI runs during each session. Each run consisted of 215 volumes sampled every 2 s at 3 mm isotropic spatial resolution (acquisition time: 7 min. TE = 30 ms, SENSE acceleration factor = 2. flip angle = 75°, 37 axial slices collected sequentially with a 1 mm gap). Subjects were instructed to rest comfortably while remaining still. One subject was excluded from further analyses due to excessive head motion. For a more detailed description of the acquisition protocol see Landman et al. (2011).

Image processing

The data was preprocessed using SPM8 (Wellcome Trust Centre for Neuroimaging, London, United Kingdom) (Friston et al., 1994) and custom MATLAB (The Mathworks, Inc., Natick, MA) scripts. Five initial volumes were discarded to allow for the stabilization of magnetization. Slice-time correction was performed using as a reference the slice acquired at the middle of the TR. Rigid body realignment transformation was performed to adjust for head motion. Structural runs were registered to the first functional frame and normalized to Montreal Neurological Institute (MNI) space using SPM8’s unified segmentation-normalization algorithm (Ashburner and Friston, 2005). The estimated nonlinear spatial transformations were applied to the rs-fMRI data, which were high-pass filtered using a cutoff frequency of 0.01 . The rs-fMRI data was spatially smoothed using a 6 mm full-width-at-half-maximum (FWHM) Gaussian kernel, which is twice the nominal size of the rs-fMRI acquisition voxel.

The Group ICA of fMRI toolbox (GIFT) (https://trendscenter.org/software/gift/) was used to estimate the number of independent components (ICs) present in the data, perform data reduction via principal component analysis (PCA) prior to independent component analysis (ICA), and perform group independent component analysis (GICA) (Calhoun et al., 2001) on the PCA-reduced data. The number of ICs was estimated using the minimum description length (MDL) criterion (Li et al., 2007). Across subjects and sessions, 56 was the maximum estimated number of ICs and 39 the median. Prior to GICA, the image mean was removed from each time point for each session, and three steps of PCA were performed. Individual session data were reduced to 112 principal components (PCs), which were concatenated within subjects in the temporal direction and further reduced to 56 PCs. Finally, data were concatenated across subjects and reduced to 39 PCs. The dimensions of the individual session PCA (112) was chosen by doubling the estimated maximum IC number (56), to ensure robust back-reconstruction (Allen et al., 2011; 2012) of subject- and session-specific spatial maps and time courses from the group-level ICs. ICA was repeated on these 39 group-level principal components 10 times, utilizing the Infomax algorithm with random initial conditions (Bell and Sejnowski, 1995). The resulting 390 ICs were clustered across iterations using a group average-link hierarchical strategy, and 39 aggregate spatial maps were defined as the cluster modes. Subject- and session-specific spatial maps and time courses were generated from these aggregate ICs using the GICA3 algorithm.

The spatial distribution of each IC was compared to a publicly available set of 100 unthresholded t-maps of ICs estimated using rs-fMRI data collected from 405 healthy participants (Allen et al., 2014). These maps were pre-classified as resting-state networks (RSNs) or noise by a group of experts, and the 50 components classified as RSNs have been organized into seven functional groups, namely visual (Vis), auditory (Aud), somatomotor (SM), default mode (DMN), cognitive-control (CC), subcortical (SC) and cerebellar (Cb) networks. For each spatial map, we calculated the percent variance explained by each of the seven sets of RSNs. The functional assignment of each Kirby component was determined by the set of components that explained the most variance, and if the top two sets of RSNs explained less than 50% of the variance in a Kirby component, the component was labeled as noise. In total 21 of the 39 components were assigned to a RSN. Subject- and run-specific time series from these components served as input for our analyses.

Analysis

The Kirby21 dataset consists of 21 ROIs measured over 210 time points for 20 subjects in two repeated sessions. The framework described in Section 2.1 (see Fig. 1), using a band-pass filter with range [0.03, 0.07] , was applied to the data for each session separately to compute the region-wise instantaneous phase for each of the 20 subjects. For each pair of subject-specific phase time series, we applied the WPS and IPS methods. We also applied CSW for comparison purposes. To facilitate comparison between the WPS and CSW methods, we used a common window length of 28 time points. We further compared the results with a prewhitened Correlation-based Sliding Window (PW-CSW) assuming an AR(1) model. This comparison was performed as a previous study (Honari et al., 2019) showed that prewhitening the data prior to analysis can lower the variance of the estimated TVC and improve brain state estimation.

For each of the two sessions, application of each method gave rise to a subject-specific 21 × 21 × 210 array of PS measures. Following the approach of Allen and colleagues (Allen et al., 2014), we applied k-means clustering to estimate recurring brain states across subjects. First, we reorganized the lower triangular portion of each subject’s dynamic correlation data into a matrix of dimension 210 × 210. Here the row dimension corresponds to the number of elements in the lower triangular portion of the matrix (i.e., 21(21 − 1)/2), and the column dimension corresponds to the number of time points. Then we concatenated the data from all subjects into a matrix with row dimensions 210 and column dimensions (210 × 20 = 4200). Finally, we applied k-means clustering to the concatenated data, where each of the resulting cluster centroids were assumed to represent a recurring brain state. The k-means clustering was repeated 200 times, using random initialization of centroid positions, in order to increase the chance of escaping local minima. In this study, we set the number of centroids to two, representing two distinct brain states, as determined using the Davies-Bouldin Index (DBI) (Davies and Bouldin, 1979). This is consistent with the number of the clusters for this dataset used in previous studies by Choe et al. (2017) and Honari et al. (2019). In order to validate the reproducibility of the brain states estimated using each of the phase synchronization measures, the results obtained from each of the two sessions were further compared and contrasted.

Application to the Human Connectome Project (HCP)

To validate the results obtained using the Kirby21 dataset, we applied the proposed methodology to a second rs-fMRI dataset from the Human Connectome Project (HCP). This data set has a different repetition time and scan length than the Kirby21 data, thus providing valuable insight into the performance of our method in this alternate setting. We used 50 subjects from the Human Connectome Project 500 Parcellation+Timeseries+Netmats (HCP500-PTN) release (Van Essen et al., 2013) to evaluate the performance of the PS measures. For space reasons, the image acquisition, image processing and analysis are discussed in Section S1.4 of the Supplementary Materials.

Results

Simulation 1

Fig. 2 shows a single realization of Simulation 1 for illustration purposes. Panel (a) shows a randomly generated pair of time series, and (b) the extracted instantaneous phases between the two time series at each time point. Since this is null data, the phase difference should vary uniformly in the interval [0, 2π] as illustrated in Panels (c) and (d).

Fig. 2.

A single realization of Simulation 1. (a) A pair of signals x(t) and y(t) generated from an independent Gaussian process. (b) The difference in the estimated phase between the signals at each time point. (c) The circular distribution of the phase difference time course in a polar coordinate system. (d) Same results in histogram form.

Fig. 3 shows a summary of the results for 1000 surrogate data sets with the analysis performed on the data prior to band-pass filtering. The mean and 95% confidence interval for each measure are shown at each time point. Results for the WPS measures (Panels (a)–(c)) are shown for each window length (30, 60, and 120 time points). The results illustrate that all measures of PS are roughly constant across time.

Fig. 3.

Results of Simulation 1 without band-pass filtering. The bold line indicates the estimated value, while the shaded area represents the 95% confidence interval. Results are shown for: (a) PLV using a sliding window; (b) circular-circular correlation using a sliding window; (c) toroidal-circular correlation using a sliding window; (d) phase coherence; and (e) CRP. The sliding window techniques are evaluated at three different window lengths.

Note that measures such as PLV and phase coherence take values between 0 and 1. The mean value using phase coherence is roughly 0.35 (Panel (d)), which is consistent with the results obtained using PLV for a window size of 60 (Panel (a)). As the window size decrease, the value of PLV tends to be lower. In contrast, circular-circular correlation, toroidal-circular correlation, and CRP all take values between −1 and 1 The mean of the CRP is 0 at each time point (Panel (e)), which is consistent with the results of the WPS obtained using circular-circular and toroidal-circular correlation (Panels (b) and (c)).

It is important to note that phase coherence and CRP preserve the temporal resolution of the phase difference as they are not estimated using a sliding window. However, this appears to come at the cost of increased variability as indicated by the relatively wider 95% confidence intervals.

The effect of the chosen window length on various WPS measures shows that for this simulation as the window size increases, the estimates converge towards their true values (i.e., 0 for circular-circular correlation and toroidal-circular correlation). The average bias and variance of each method is reported in Table S1.

Fig. 4 shows the comparison between various measures of PS used on the surrogate data after band-pass filtering. The results indicate that synchronization measures remain roughly constant across time. However, the WPS measures (Panels (a)–(c)) show a noticeable difference compared to the results without band-pass filtering. Significantly, the results of the WPS measures show inflated values, indicating a higher degree of phase synchronization than would be expected in a null setting. PLV, circular-circular correlation, and toroidal-circular correlation take values around 0.84, 0.53, and 0.62, respectively. These is driven by the fact that both signals have a center frequency of 0.05 after applying band-pass filtering, leading to a situation where the signals are constrained to remain phase locked throughout the time course. The IPS measures are not similarly affected. The mean value using phase coherence is again roughly 0.35 (Panel (d)), and the mean of the CRP is 0 at each time point (Panel (e)). Both values are roughly equivalent to those seen before band-pass filtering. The average bias and variance of each method is reported in Table S2.

Fig. 4.

Results of Simulation 1 with band-pass filtering. The bold line indicates the estimated value, while the shaded area represents the 95% confidence interval. Results are shown for: (a) PLV using a sliding window; (b) circular-circular correlation using a sliding window; (c) toroidal-circular correlation using a sliding window; (d) phase coherence; and (e) CRP. The sliding window techniques are evaluated at three different window lengths.

Figure S1 in the Supplementary Materials show the results of this simulation using the BEMD-based PS framework. The results are similar to those obtained using band-pass filtering, with the exception that the variance of the estimates are slightly higher; see Table S3. This can be explained by the well-known sensitivity of the EMD-based approach to noise (Imaouchen et al., 2017; Wu and Huang, 2009; Zheng et al., 2014), which in turn leads to a higher variance of the phase synchronization measures.

Simulation 2

Figs. 5 and 6 illustrate the results of Simulation 2 performed on the data before and after band-pass filtering, respectively. Recall that in this simulation the two signals are designed to have the same phase up to time t = 170, after which a phase shift is introduced that varies linearly from 0 to 4π (see Fig. 5a). Thus, the signals should gradually move in and out of phase during the second half of the time course. Here the signals will be in-phase when the phase difference is 2π and 4π, and in anti-phase when the difference is π and 3π.

Fig. 5.

Results of Simulation 2 without band-pass filtering. (a) The ground truth phase shift between the two signals as a function of time. Results are shown for: (b) PLV using a sliding window; (c) circular-circular correlation using a sliding window; (d) toroidal-circular correlation using a sliding window; (e) phase coherence; and (f) CRP. The sliding window techniques are evaluated at three different window lengths.

Fig. 6.

Results of Simulation 2 with band-pass filtering. (a) The ground truth phase shift between the two signals as a function of time. Results are shown for: (b) PLV using a sliding window; (c) circular-circular correlation using a sliding window; (d) toroidal-circular correlation using a sliding window; (e) phase coherence; and (f) CRP. The sliding window techniques are evaluated at three different window lengths.

The costs of not band-pass filtering the data are apparent in Fig. 5, as all five of the methods return results consistent with those seen in the null setting. None of the methods does a good job of either detecting the fact that the signals are in phase in the first half of the time course, or that they gradually move in and out of phase in the second half. This can be explained by the fact that the signal is contaminated with noise from all frequencies, which in turn corrupts the estimated instantaneous phase.

Contrast this with the results after band-pass filtering shown in Fig. 6. Here all of the measures of PS correctly predict a value close to 1 in the first half of the signal, indicating that all methods are picking up on the fact that the signals are in phase. In Panels (b)–(d), which represents the WPS measures, the phase shift occurring after t = 170 leads to a decrease in phase synchronization from this time point on. The toroidal-circular correlation appears to perform best, showing more sensitivity in detecting the episodes of phase synchronization compared to PLV and circular-circular correlation. It can also be observed that the circular-circular correlation is more susceptible and sensitive to the noise than the other measures (see the increased wiggles in the estimates values). Interestingly, the PLV results appear to be more sensitive to the window length used than the other two metrics. However, it is important to note that none of the WPS methods are able to detect that the signals are in phase when the phase shift equals 2π and 4π.

In contrast, both the phase coherence and CRP better captures the PS variation than the WPS measures. This is partly due to the fact that using a sliding window deteriorates the resolution of the PS depending on the window size. In particular, note how well CRP detects that the signal is in phase at points when the phase shift equals 2π and 4π. This is in contrast to phase coherence that erroneously assumes that signals are also in phase when the shifts are equal to π and 3π. The latter is due to the fact that phase coherence cannot differentiate between when the signals are in phase from when they are in anti-phase. The benefits of band-pass filtering can be clearly seen in Tables S1 and S2, which show a reduction in bias when the data is properly filtered.

Figure S2 and Table S3 in the Supplementary Materials show the results of this simulation using the BEMD-based PS framework. Again, the results are consistent with those seen using band-pass filtering, though the variance of the estimates are somewhat higher.

Simulation 3

Figs. 7 and 8 illustrate the results of Simulation 3 performed on data before and after band-pass filtering. As illustrated in Fig. 7a, the two signals are designed to initially be in phase, after which they gradually go out of phase. At time t = 170 when the phase difference is 2π, the two signals will be in anti-phase, before returning to being in phase at the end of the time course.

Fig. 7.

Results of Simulation 3 without band-pass filtering. (a) The ground truth phase shift between the two signals as a function of time. Results are shown for: (b) PLV using a sliding window; (c) circular-circular correlation using a sliding window; (d) toroidal-circular correlation using a sliding window; (e) phase coherence; and (f) CRP. The sliding window techniques are evaluated at three different window lengths.

Fig. 8.

Results of Simulation 3 with band-pass filtering. (a) The ground truth phase shift between the two signals as a function of time. Results are shown for: (b) PLV using a sliding window; (c) circular-circular correlation bsing a sliding window; (d) toroidal-circular correlation using a sliding window; (e) phase coherence; and (f) CRP. The sliding window techniques are evaluated at three different window lengths.

The costs of not band-pass filtering the data are again apparent in Fig. 7, as all of the methods show results consistent with the null setting. This is in contrast to the results obtained after band-pass filtering shown in Fig. 8. Here all measures of phase synchronization pick up on the fact that the signals start out in phase, gradually goes out of phase (culminating at time t = 170), before gradually return to being in phase. Tables S1 and S2 show a reduction in bias when the data is properly filtered, except for PLV.

In Panels (b)–(d), which represents the WPS measures, we see that using a longer window length tends to capture phase dynamics better than using a smaller window length. Again, the toroidal-circular correlation performs best, showing increased sensitivity in detecting the episodes of phase synchronization compared to circular-circular correlation and PLV. Increased window lengths provide better results. Both phase coherence and CRP capture the manner in which phase synchonization varies more clearly than the WPS measures. In particular, CRP provides the most reliable measures in this simulation and clearly detects both when the signals are in and out of phase. In comparison, phase coherence cannot separate when the signals are in phase and anti-phase, illustrating one of the shortcomings of the approach.

Figure S3 and Table S3 in the Supplementary Materials shows the results of this simulation using the BEMD-based PS framework. Again, the results are consistent with those seen using band-pass filtering, though the variance of the estimates are somewhat higher.

Analysis of rs-fMRI data

After applying each method to the Kirby21 rs-fMRI data, two brain states were extracted using k-means clustering. Fig. 9 contrasts the estimated brain states obtained using the different methods for assessing phase synchonization, as well as with correlation-based sliding window analysis (both with and without pre-whitening). Brain states are organized into seven functional groups as described earlier in Section 2.4.2.

Fig. 9.

Analysis of the Kirby21 Data. After applying each method the time-varying connectivity measures were clustered into 2 reoccurring brain states. Results are shown for both sessions (top to bottom) for: PLV using a sliding window; circular-circular correlation (CIRC) using a sliding window; toroidal-circular correlation (TORC) using a sliding window; phase coherence (PC); cosine of the relative phase (CRP); correlation-based sliding window (CSW); and prewhitened correlation-based sliding window (PW-CSW). The sliding window techniques are evaluated with window length 30 time points.

Beginning with the WPS methods, PLV (top row), circular-circular correlation (second row) and toroidal-circular correlation (third row) show roughly similar results with regards to the relationship between functional groups in each brain state. However, the PLV derived brain states in general take higher values than those obtained using toroidal-circular correlations, which in turn take higher values than those obtained using circular-circular correlations. These results are largely consistent with those seen in the simulation studies, and the fact that toroidal-circular and circular-circular correlations take a wider range of values (compared to PLV which is constrained between 0 and 1).

Turning to the IPS methods, the brain states obtained using phase coherence (fourth row) tends to provide higher values in general compared to CRP (fifth row). This is not necessarily surprising as the range of potential values are different (phase coherence takes values between 0 and 1, while CRP takes values between −1 and 1). In addition, as seen in the simulation studies phase coherence has problems differentiating between when signals are in phase versus when they are in anti-phase. Together, this provide higher values in the estimated brain states. For example, State 2 shows a hyper-connected state where all PS measures are close to 1.

Interestingly, the results obtained using sliding windows without any prewhitening (sixth row) is very similar to those obtained using the CRP. Across both sessions and estimation methods, State 2 was characterized by stronger correlations (both positive and negative) relative to State 1. Moderate to strong negative correlations between sensory systems (auditory, somatomotor, and visual) components were present in State 2 but were reduced in State 1. Similarly, negative correlations within the DMN (grey) components were present in State 2 and were reduced in State 1. This similarity between methods indicates that CRP may be finding similar brains states as CSW, but using more high-resolution data as it does not use a predefined window. These findings are largely consistent with those of Pedersen et al. (2018) who found that IPS and CSW conveyed comparable information of time-resolved fMRI connectivity, though IPS provided finer temporal resolution. The results obtained using sliding windows without prewhitening (seventh row) show lower estimated values than the results without prewhitening, which is consistent with results found in Honari et al. (2019).

Finally, we found that both States 1 and 2 were highly reliable across sessions regardless of estimation method. The results for CRP and the sliding window approaches (described above) appear more interpretable than those obtained using the other methods. In particular the results for PLV and PC, which takes values between 0 and 1 appear to show extremely hyperconnected states, perhaps reflective of the lower range of values that these measures return.

Figure S8 in the Supplementary Materials shows the results for the analysis of the HCP dataset. Here after applying each method to the rs-fMRI data, four brain states were extracted using k-means clustering. Brain states are organized into the same seven functional groups as before. The results show a great deal of similarity between the HCP and Kirby21 results with respect to the properties of the estimated brain states obtained using each PS measure. Again, PLV derived brain states (top row) take higher values than those obtained using toroidal-circular correlations (third row), which in turn take higher values than those obtained using circular-circular correlations (second row). Similarly, the brain states obtained using phase coherence (fourth row) tend to provide higher values compared to CRP (fifth row). Finally, the results obtained using sliding windows without any prewhitening (sixth row) is very similar to those obtained using the CRP.

Focusing on the CRP results, States 1, 2 and 4 all show moderate to high correlations among signal components representing sensory systems (i.e., visual, somatomotor, and auditory components). In States 1 and 4, a set of components in the cerebellum showed negative correlations with sensory components. The HCP states were largely similar to those obtained from the Kirby21 data, though it is important to note that the number and placement of the components in each HCP RSN do not map directly onto one another.

Discussion

There is growing interest in measuring time-varying functional connectivity between time courses from different brain regions using rs-fMRI data. One such approach is to measure their phase synchronization across time. In this paper, we evaluate a number of methods for measuring PS and contrast them with one another. In discussing methods, we differentiate between two classes of methods: windowed phase synchronization and instantaneous phase synchronization.

WPS methods combine a static PS measure between two different signals with a sliding window to obtain a time-varying measure of PS. In principal, any metric that allows one to calculate an omnibus measure of PS can be used within this framework. Since phase information is circular data, the use of circular-circular correlation and toroidal-circular correlation were deemed natural candidate methods to use as a measure of PS. To the best of our knowledge, neither approach has previously been used to study PS in fMRI. The PLV in WPS method, has in contrast previously been used to assess episodes of elevated gastric-BOLD synchronization (Rebollo et al., 2018).

IPS methods directly use the phase difference time series obtained from applying the Hilbert transform, allowing one to compute an instantaneous measure of PS. This has the benefit of providing a higher temporal resolution, as there is no need to choose an arbitrary window size as for the WPS methods. However, there remains a related somewhat arbitrary choice of filter bandwidth to narrow-band the signals prior to analysis. Here we focused on two measures of IPS, phase coherence, which has already found wide usage in the field, and CRP, a newly developed method.

The three simulations illustrate several important points regarding the performance of these methods. Simulation 1 shows that the WPS methods are highly affected by band-pass filtering. To illustrate, Fig. 4 shows that these methods tend to provide estimates that indicate that signals are consistently in phase, even when the phases are designed to randomly vary. These results hold because the two signals being compared have a center frequency of 0.05 Hz after applying band-pass filter with cutoff frequencies [0.03, 0.07] Hz. This leads to a situation where the signals are constrained to remain relatively phase locked and thus have constant PS throughout the time course. Importantly, the IPS results appear to perform similarly on the data both before and after band-pass filtering (see Figs. 3 and 4), and thus appear to be less sensitive to filtering in the null setting.

While at first glance, the results of Simulation 1 appear to indicate that band-pass filtering is not beneficial, and may in fact be detrimental, Simulations 2 and 3 put this notion to rest. Here, the results performed on the non band-pass filtered data indicate that none of the methods are able to pick up changes in real PS present in the data, and instead appear to erroneously indicate that the data behave in manner consistent with null data. This is largely corrected after band-pass filtering the data (see Figs. 6 and 8). This result holds both for WPS and IPS methods, and indicates that band-pass filtering is a necessary step in the analysis of PS.

This result corresponds to theoretical findings (Bedrosian’s theorem) that suggest using band-pass filters in the study of PS is critical for the signal to have physically meaningful demodulation into its envelope and instantaneous phase components. However, it is important to note that band-pass filtering comes at the cost of introducing further autocorrelation into the phase of the signal. In addition, band-pass filtering increases the risk of spurious detection of phase synchronization (Rosenblum et al., 2000). While a band-pass filter denoises the signal, it can also lead to an increase in the degree of synchronization by narrowing the band width; see Fig. 4.

The results of Simulations 2 and 3 together show that all methods to a certain extent were able to detect changes in PS. Focusing on the WPS measures, toroidal-circular correlation performed best, showing increased sensitivity to detecting episodes of PS compared to PLV and circular-circular correlation. Circular-circular correlation was the most susceptible and sensitive to noise. A previous study comparing PLV and circular-circular correlation (Pauen and Ivanova, 2013) suggested that circular-circular correlation is appropriate for estimating the phase coupling reliably and not restricted to bivariate analyses. It also indicated that using it as a measure of phase coupling could show slightly lower estimates than its counterpart. This result is consistent with what we found in our simulations.

When assessing WPS measures, we also investigated a variety of window lengths. The simulations indicated that shorter windows yielded a higher estimate of phase synchronization and increased risk of detecting spurious PS. However, longer windows made it harder to detect subtle changes. In general, longer window lengths tend to provide more accurate estimates of PS as they lead to a decrease in the variation of the estimates. However, this may be only the case if the pairwise phase coupling between the underlying signals is of semi-stationary nature with no abrupt changes.

It is important to note that if overlapping windows are used, an autocorrelation (beyond that already present due to the smoothness of the BOLD signal) is induced in the estimated WPS values, making changes in connectivity appear artificially smooth. Overlapping windows do provide a finer temporal resolution than non-overlapping windows, as the latter is a subset of the former. However, the exact gain in resolution is hard to quantify as it depends heavily on the autocorrelation between windows. One of the goals of IPS is to gain a time-resolved resolution. However, it is worth noting that the band-pass acts as an implicit window due to Parseval’s therorem.

The simulations clearly point to the benefits of band-pass filtering when computing PS. An alternative approach is to instead use empirical mode decomposition to extract frequencies of interest. This approach has been widely used to analyze EEG data, but also found recent use in application to fMRI data (Niazy et al., 2011; Yuen et al., 2019; Zhou et al., 2020). In the Supplemental Materials the results of the simulations described above were contrasted with results obtained using an EMD-based PS framework. Both approaches yielded similar results for the different PS measures. However, the variance of the estimates obtained using EMD were somewhat higher for the same level of noise compared to band-pass filtering.

IPS measures consistently outperformed WPS measures in the simulations, and were able to better pick up changes in PS across time. While phase coherence offers more accurate and sensitive results than the WPS methods, it still discards information about the direction of the relationship. In contrast, CRP was not only able to detect phase synchronization but also preserve the directional information contained in the relative phase difference of the signals. However, one should note that the variation present in the IPS methods appears larger than WPS methods as evidenced by the wider confidence bounds.

It is interesting to consider the range of values each method returns. PLV and phase coherence take values between 0 and 1. In contrast, circular-circular correlation, toroidal-circular correlation, and CRP all take values between −1 and 1. This has to be taken into consideration when interpreting the results of each method. For example, in Simulation 1 where we analyzed null data, the latter methods returned values that lay symmetrically around 0. This makes it easier to interpret null results compared to PLV and phase coherence whose null values were around 0.35. Understanding what null values look like is a critical component towards understanding the performance of a method, as it is otherwise difficult to differentiate signal from noise.

Application to real rs-fMRI data showed results that were consistent with the simulations. The WPS methods showed roughly equivalent results with respect to the relationship between functional groups in each estimated brain state. As described above, the shorter range of values for PLV and phase coherence made it more difficult to pick up subtle differences between brain states, and they both returned a hyper-connected brain state where all PS measures are close to 1. Interestingly, the brain states estimated using CRP closely resembled those estimated using sliding window correlations. Thus, it appears that CRP is finding similar brains states but using more high-resolution data as it does not require the use of a sliding window. Finally, there was a large degree of similarity between states estimated from the Kirby21 and HCP datasets, even though they differ in important aspects including repetition time, scan length, and sample size.

The focus of this work is on showing the potential of the PS framework and highlighting the necessity of using appropriate PS measures to estimate time-varying functional connectivity. However, we believe the approach can be extended to allow for the analysis of directed connectivity between regions of interest. For example, this could be studied using instantaneous phase time series from fMRI data in conjunction with transfer entropy (phase-based transfer entropy) (Lobier et al., 2014).

In summary, we recommend the use of CRP as a measure of PS as it is able to separate when the signals are in phase from when they are in anti-phase. In addition, it returns a range of values similar to correlation, which makes it possible to interpret results similarly. Of all the methods tested, it showed the greatest concurrence with CSW, with the benefit of not having to predefine a window length.