key: cord-0504962-w9asj1w7
authors: Honari, Hamed; Choe, Ann S.; Electrical, Martin A. Lindquist Department of; Engineering, Computer; University, Johns Hopkins; Imaging, USA F. M. Kirby Research Center for Functional Brain; Institute, Kennedy Krieger; Injury, USA International Center for Spinal Cord; Radiology, USA Russell H. Morgan Department of; Science, Radiological; Medicine, Johns Hopkins School of; Biostatistics, USA Department of; USA,
title: Evaluating phase synchronization methods in fMRI: a comparison study and new approaches
date: 2020-09-21
journal: nan
DOI: nan
sha: 4b68a9ee250379d4e666664ab055970272d6998a
doc_id: 504962
cord_uid: w9asj1w7

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 here 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 necessity of band-pass filtering the data to obtain valid results. 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.

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 here 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.

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 (Chang and Glover, 2010; Hutchison et al., 2013a; Preti et al., 2016; Tagliazucchi et al., 2012; Thompson et al., 2013; Allen et al., 2014; Lurie et al., 2019) 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 (Tagliazucchi et al., 2012; Chang and Glover, 2010; Hutchison et al., 2013a,b) , change point analysis (Cribben et al., 2012 (Cribben et al., , 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) (Eavani et al., 2013; Vidaurre et al., 2017; Shappell et al., 2019; Bolton et al., 2018) . 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; Lindquist et al., 2014; Laumann et al., 2016) .

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 Pedersen et al., , 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 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) (Glerean et al., 2012; Boccaletti et al., 2018; 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: circularcircular 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) 

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 necessity of band-pass filtering the data to accurately estimate the instantaneous phases. 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 the IPS. We contrast these methods through a series of simulations and application to rs-fMRI data.

To obtain the instantaneous phase (Boccaletti et al., 2018) of an arbitrary real signal x(t) one must first construct an analytic signal:

where j = √ −1 and H represents the Hilbert Transform. This signal can in turn be re-expressed as follows:

where A(t) represents the envelope and φ(t) 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 monocomponent signal and the Hilbert transform will produce an analytic signal with meaningful envelope and phase. 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 Figure 1 . Consider that a pair of time series x(t) and y(t), Figure 1 : A schematic of the approach to calculate the instantaneous phase (IP) framework t = 1, . . . T , from two different ROIs are filtered using a narrow band-pass zero-phase filter, h bp (t), and denote the filtered data by x n (t) and y n (t) 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 φ x (t) and φ y (t) 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.

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.

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.

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 (Rosenblum et al., 2000; Halliday et al., 1998) .

This notion of synchronization can be expressed as follows:

Here the integers m and n are the synchronization indices and ∆Φ m,n (t) the generalized phase difference time series.

In this paper, we assume m = n = 1 and drop the indices and let ∆Φ m,n (t) = ∆Φ(t).

Using the instantaneous phase difference of the signals at each time point, the PLV can be computed as follows:

where the operator · t denotes averaging over time. If the pair of signals are unsynchronized, then P LV = 0 and ∆φ(t) follows a uniform distribution;

otherwise, if the pair are synchronized, P LV 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 Equation (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.

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, Φ x = (φ x (1), . . . φ x (T )) and Φ y = (φ y (1), . . . φ y (T )), while µ and ν represent the mean directions of Φ x and Φ y , respectively. Thus, the terms sin(Φ x − µ) and sin(Φ y − ν) can be interpreted as the deviations of Φ x and Φ y 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μ t andν t 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, µ t is computed as follows (Jammalamadaka and Sengupta, 2001; Bishop, 2006) :μ

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 cos φ x (i) = x i and sin φ y (i) = y i , for i = 1, ..., n. Thus,x = n −1 cos φ x (i),ȳ = n −1 sin φ y (i), and the mean direction thus can be written as expressed in Eq. 11.

Hereν t is computed analogously.

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 φ x (t 1 ) (or similarly φ y (t 1 )) and φ x (t 2 ) (or φ y (t 2 )) be two circular data points and set

The same definition holds for φ y (t 1 ) and φ y (t 2 ).

, then the order function can be expressed as follows:

Now, let us assume that (φ

The circular correlation is then defined as follows:

Based on this definition, the two variables φ x and φ y move on the circumfer-

Similarly, ρ tor < 0 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.

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 (Pedersen et al., , 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.

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 | sin(−∆Φ)| = | sin(∆Φ)|. 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.

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., |∆Φ(t)| ≈ 0), CRP yields a value close to 1. When the phases are dissimilar but in the same direction, their relative phase difference is bounded between [−π/2,π/2], 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.

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 rs-fMRI data 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 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] Hz. Throughout we used a 5 th 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.

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 ϕ x (t). Further, let y(t) have the same angular frequency but with phase ϕ y (t). The signals can be expressed as follows:

Without loss of generality, let ϕ x (t) = 0 and ϕ y (t) be a ramp function,

The time series can then be expressed as follows:

Throughout, ω 0 = 2πf rad/s with f = 0.05Hz, and the transition is set to occur at t 0 = 170 s. The noise terms ε x and ε y 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 t 0 = 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 ϕ x (t) = 0 and ϕ y (t) be a sigmoid function, i.e.

Hence, the time series can be expressed as follows: 

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 Hz. The rs-fMRI data was spatially smoothed using a 6 mm full-width-at-half-maximum 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 

Data from a single run of the Kirby21 was used. Thus, the data consisted of 21

ROIs measured over 210 time points for 20 subjects. The framework described in Section 2.1 (see Fig. 1 ), using a band-pass filter with range [0.03, 0.07] Hz, was applied to the data 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.

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) . 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 as the window size increases, the estimates converge towards their true values (i.e., 0 for circular-circular correlation and toroidal circular correlation). 

Figures 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 Figures 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π. The costs of not band-pass filtering the data are apparent in Figure 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 Figure 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 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.

Figures 7 -8 illustrate the results of Simulation 3 performed on data before and after band-pass filtering. As illustrated in Figures 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. The costs of not band-pass filtering the data are again apparent in Figure 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 Figure   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.

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 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.

After applying each method to the rs-fMRI data, two brain states were extracted using k-means clustering. Figure 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: visual (Vis); auditory (Aud); somatomotor (SM); default mode (DMN); cognitive-control (CC); sub-cortical (SC); and cerebellar (Cb) networks.

Beginning with the WPS methods, PLV (top row), circular-circular correlation (second row) and toroidal 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 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 and circular-circular correlations take a wider range of values (compared to PLV which is constrained between 0 and 1). Figure 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 (top to bottom) for: PLV using a sliding window; circular-circular correlation using a sliding window; toroidal correlation using a sliding window; phase coherence; cosine of the relative phase; correlation-based sliding window; and prewhitened correlation-based sliding window. The sliding window techniques are evaluated with window length 30 time points.

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. 

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 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 bandpass 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 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 phase synchronization. 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.

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 preserved 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 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.

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.

Tracking whole-brain connectivity dynamics in the resting state

A baseline for the multivariate comparison of resting-state networks

Capturing inter-subject variability with group independent component analysis of fmri data: a simulation study

Unified segmentation

An information-maximization approach to blind separation and blind deconvolution

Pattern recognition and machine learning

Synchronization: from coupled systems to complex networks

Interactions between large-scale functional brain networks are captured by sparse coupled hmms

A method for making group inferences from functional mri data using independent component analysis

Time-frequency dynamics of resting-state brain connectivity measured with fmri

Comparing test-retest reliability of dynamic functional connectivity methods

Dynamic connectivity regression: determining state-related changes in brain connectivity

Detecting functional connectivity change points for single-subject fmri data

A cluster separation measure

Unsupervised learning of functional network dynamics in resting state fmri

Statistical parametric maps in functional imaging: a general linear approach

Functional magnetic resonance imaging phase synchronization as a measure of dynamic functional connectivity

Using electroencephalography to study functional coupling between cortical activity and electromyograms during voluntary contractions in humans

Investigating the impact of autocorrelation on time-varying connectivity

Dynamic functional connectivity: promise, issues, and interpretations

Resting-state networks show dynamic functional connectivity in awake humans and anesthetized macaques

A correlation coefficient for angular variables. Statistical theory and data analysis II

Topics in circular statistics

Transient brain activity disentangles fmri resting-state dynamics in terms of spatially and temporally overlapping networks

Measuring phase synchrony in brain signals

Characterizing instantaneous phase relationships in whole-brain fmri activation data

Surrogate data for hypothesis testing of physical systems

Multiparametric neuroimaging reproducibility: a 3-t resource study

On the stability of bold fmri correlations

Estimating the number of independent components for functional magnetic resonance imaging data

Evaluating dynamic bivariate correlations in resting-state fmri: a comparison study and a new approach

Time-varying functional network information extracted from brief instances of spontaneous brain activity

Questions and controversies in the study of time-varying functional connectivity in resting fmri

Spectral characteristics of resting state networks

Circular correlation coefficients versus the phase-locking-value

Spontaneous brain network activity: Analysis of its temporal complexity

On the relationship between instantaneous phase synchrony and correlation-based sliding windows for time-resolved fmri connectivity analysis

Resting-state temporal synchronization networks emerge from connectivity topology and heterogeneity

The dynamic functional connectome: State-of-the-art and perspectives

Sense: sensitivity encoding for fast mri

Stomach-brain synchrony reveals a novel, delayed-connectivity restingstate network in humans

Detection of phase locking from noisy data: application to magnetoencephalography

Phase synchronization of chaotic oscillators

Improved state change estimation in dynamic functional connectivity using hidden semi-markov models

The equivalence of the torus and the product of two circles in homotopy type theory

Echo-planar imaging: magnetic resonance imaging in a fraction of a second

Spontaneous bold event triggered averages for estimating functional connectivity at resting state

Dynamic bold functional connectivity in humans and its electrophysiological correlates

Shorttime windows of correlation between large-scale functional brain networks predict vigilance intraindividually and interindividually

Brain network dynamics are hierarchically organized in time

Dynamic connectivity detection: an algorithm for determining functional connectivity change points in fmri data

On circular correlation for data on the torus

The work presented in this paper was supported in part by NIH grants R01 EB016061 and R01 EB026549 from the National Institute of Biomedical Imaging and Bioengineering and R21 NS104644 from the National Institute of Neurological Disorders and Stroke.