Static and dynamic functional network connectivity (FNC) are typically studied separately, which makes us unable to see the full spectrum of connectivity in each analysis. Here, we propose an approach called filter-banked connectivity (FBC) to estimate connectivity while preserving its full frequency range and subsequently examine both static and dynamic connectivity in one unified approach.
First, we demonstrate that FBC can estimate connectivity across multiple frequencies missed by a sliding-window approach. Next, we use FBC to estimate FNC in a resting-state fMRI dataset including schizophrenia patients (SZ) and typical controls (TC). The FBC results are clustered into different network states. Some states showed weak low-frequency strength and as such were not captured in the window-based approach. Additionally, we found that SZs tend to spend more time in states exhibiting higher frequencies compared with TCs who spent more time in lower frequency states. Finally, we show that FBC enables us to analyze static and dynamic connectivity in a unified way. In summary, FBC offers a novel way to unify static and dynamic connectivity analyses and can provide additional information about the frequency profile of connectivity patterns.
Functional connectivity and its cross-network analog, functional network connectivity (FNC), have been the focus of many neuroimaging studies over the past few decades. As the methods used to estimate FNC can be used to estimate functional connectivity (abbreviated as FC in some papers) and vice versa in most cases, in this article we will use the term FNC when talking about connectivity that includes FC too.
Methods developed to estimate FNC can be grouped into two major categories: those that assume connectivity among different networks of the brain is constant through time (static FNC; sFNC) and those that assume temporal variation in connectivity (dynamic FNC; dFNC). The sFNC and dFNC approaches have proven to be extremely informative about both healthy (Allen et al., 2014; Liegeois et al., 2019; Vidaurre, Smith, & Woolrich, 2017) and disordered brain function (Damaraju et al., 2014; de Lacy, Doherty, King, Rachakonda, & Calhoun, 2017; Jin et al., 2017; Kaiser et al., 2016).
While static connectivity has resulted in many interesting findings (Di Martino et al., 2008; van den Heuvel & Hulshoff Pol, 2010), this view of connectivity is limited to the average connectivity patterns over the entire experiment. Approaches designed based on dFNC relax the assumption of static connectivity. The most common way to estimate dFNC uses a sliding window to estimate time-varying connectivity. Typically, a sliding window is paired with Pearson correlation (SWPC) to estimate time-varying connectivity (Allen et al., 2014; Faghiri, Stephen, Wang, Wilson, & Calhoun, 2018; Hutchison, Womelsdorf, Gati, Everling, & Menon, 2013; Kucyi & Davis, 2014). Using a sliding window to estimate dFNC has the benefit of being straightforward but has two major shortcomings. First, we need to choose a window size for any sliding-window approach. We want to choose a window size that is large enough so that the standard deviation is as small as possible. At the same time, the window size should be small enough to allow us to detect faster changes in dFNC (Hutchison, Womelsdorf, Allen, et al., 2013). The second shortcoming of the sliding-window approach is its low-pass nature, which has been reported previously (Leonardi & Van De Ville, 2015; Sakoglu et al., 2010; Thompson & Fransson, 2015). This tells us that, regardless of the chosen window size, the estimated dFNC is subjected to a low-pass filter and therefore the full frequency range of dFNC is not captured. This may be the reason that Shakil, Billings, Keilholz, and Lee (2018) found that using a constant window size for SWPC is not a reliable solution to study dynamic connectivity. Note that sliding window can also be used with other estimators such as multiplication of temporal derivatives (Shine et al., 2015) and weighted average of shared trajectory (Faghiri et al., 2020; Faghiri, Stephen, Wang, Wilson, & Calhoun, 2019). Another dFNC/FC estimator that is getting attraction recently is instantaneous phase synchrony, which defines connectivity as synchrony between two time series (Kaboodvand, Iravani, & Fransson, 2020; Kaboodvand, van den Heuvel, & Fransson, 2019; Pedersen, Omidvarnia, Zalesky, & Jackson, 2018). A sliding window will also act as a low-pass filter if used with these and any other estimators.
Another category of methods that aim to explore the frequency profile of connectivity uses time-frequency analysis ideas. The most well-known methods in this category utilize wavelets (Mallat, 1999). While these methods have resulted in many interesting findings in functional magnetic resonance imaging (fMRI; Chang & Glover, 2010; Yaesoubi, Allen, Miller, & Calhoun, 2015; Yaesoubi et al., 2017), these studies have several limitations. First, the interpretation of the results in these studies can be challenging. The Chang and Glover implementation of wavelet resulted in a large amount of information without a way to succinctly summarize the results (Chang & Glover, 2010). To remedy this, Yaesoubi et al. (2015) proposed an approach using wavelets that can be used to study group differences, but their results are all in the wavelet domain. This presents a difficulty in comparing the results with other dFNC studies, as most dFNC studies work in the time domain (the time domain results are typically considered easier to interpret). In addition, and perhaps more importantly, both wavelet approaches perform frequency tiling in the activity domain instead of in the connectivity domain. We believe that to discuss the frequency properties of dFNC, it is important to implement all time-frequency tiling steps directly in the connectivity domain. The reason behind this statement is that the relationship between the activity and connectivity domains is unknown (and possibly nonlinear); therefore, the frequency information is distorted when transforming from the activity to the connectivity domain.
Apart from the two categories of approaches mentioned above, other methods have been proposed that do not estimate dFNC directly but rather explore different aspects of connectivity dynamics. These methods include hidden Markov models (Ou et al., 2013), hidden semi-Markov models (Shappell, Caffo, Pekar, & Lindquist, 2019), coactivation patterns using clustering approaches (Liu & Duyn, 2013), and window-less dictionary learning approaches (Yaesoubi, Adali, & Calhoun, 2018). There are also methods that use Bayesian inference (Andersen, Winther, Hansen, Poldrack, & Koyejo, 2018; Warnick et al., 2018) and graphical lasso (Cai et al., 2019) for connectivity estimation. Choe et al. (2017) explored the test-retest reliability of several connectivity estimators. For a detailed comparison of various dFNC estimators in the fMRI field, see Xie et al. (2019). A comprehensive review of connectivity metrics in the electrophysiological field can be found in O’Neill et al. (2018).
In this work, we emphasize the importance of differentiating the connectivity domain frequency profile from the activity domain frequency profile (more on this in the Discussion section). As mentioned previously, some studies have implemented frequency tiling in the activity domain and use theses tiles to make inferences about the connectivity frequency profile (which we believe is inaccurate). Here we proposed a new method to estimate connectivity that aims to implement frequency tiling directly in the connectivity domain. As a central part of this approach, we use filter banks for frequency tiling in the connectivity domain (unlike wavelet methods, where frequency tiling is implemented in the activity domain). We do this without making any assumption about the frequency profile of connectivity (unlike SWPC, where only low-pass connectivity values are estimated). Using filter banks, we are able to estimate any desired connectivity frequency bands without being limited to only low-pass bands. In addition, our proposed approach enables us to examine specific frequency bands of connectivity that include both static connectivity (connectivity at zero frequency) and dynamic connectivity (connectivity at nonzero frequencies) in one unifying approach. And, the results are in the time domain, allowing us to more easily compare our results with other studies (unlike wavelet approaches, where the final results are in the wavelet domain).
Another advantage of our approach (filter-banked connectivity; FBC) is that this method removes the impact that incorrect window selection has on the estimated connectivity values. The reason behind this statement is that here we are looking at the full spectrum of connectivity instead of limiting ourselves to any given band without strong prior knowledge. To drive this point home, we did a very simple simulation with two different situations. Assume we have two different pairs of time series, each extracted from different locations in the brain. The first pair has low-frequency connectivity, while the second pair has high-frequency connectivity. Figure 1 illustrates this simulation. For the first scenario, both SWPC and low-pass FBC have estimated the ground truth very nicely (both have high correlation values with the ground truth). But for the second scenario, only high-pass FBC has a high correlation with the ground truth. Therefore, in this situation, the connectivity information is lost if we use SWPC.
To explain the benefits of our proposed approach, we designed a toy example. Furthermore, to showcase the utilization of the proposed approach, we implemented it on a dataset used to study sFNC and dFNC previously (Damaraju et al., 2014; Yaesoubi et al., 2017). This dataset includes both typical controls and individuals diagnosed with schizophrenia.
Schizophrenia is a mental disorder that is associated with functional connectivity abnormalities (Damaraju et al., 2014; Pettersson-Yeo, Allen, Benetti, McGuire, & Mechelli, 2011). There have been a number of studies of functional connectivity in schizophrenia individuals (SZ) using resting fMRI. For example, Camchong, MacDonald, Bell, Mueller, and Lim (2011) reported hyperconnectivity between default mode network (DMN) and the rest of the brain. In another work, Jafri, Pearlson, Stevens, and Calhoun (2008) used a whole-brain approach to study the differences between typical controls (TC) and SZ. They reported that SZ showed increased connectivity between DMN and visual and frontal functional domains compared with TC. Damaraju et al. (2014) reported that SZ compared with TC shows increased connectivity between thalamus and sensory functional domains. Damaraju et al. also reported decreased static connectivity in sensory domains when comparing SZ with TC. Decreased connectivity in SZ compared with TC has been reported in other studies, as well (Dong, Wang, Chang, Luo, & Yao, 2018; Erdeniz, Serin, Ibadi, & Tas, 2017; Friston & Frith, 1995; Lynall et al., 2010; Skudlarski et al., 2010).
A few studies have evaluated the dynamic aspect of connectivity in SZ population. Using SWPC, Damaraju et al. (2014) found that SZ compared with TC tend to stay less time in states that show strong overall connectivity while they tend to spend more time in states showing weak connectivity between different domains. Other studies also reported transient reductions in both functional connectivity and network activities (Iraji, Deramus, et al., 2019; Iraji, Fu, et al., 2019). Miller et al. (2016) reported less dynamism (e.g., less change in transient connectivity patterns) in SZ compared with TC. For a recent review on connectivity-related findings (both static and dynamic findings) in SZ population, see Mennigen, Rashid, and Calhoun (2019).
A different set of studies have explored the spectral properties of fMRI time series between SZ and TC. An earlier study found that the frequency profile of default mode is altered in SZs compared with TCs (Garrity et al., 2007). In addition, Fryer et al. (2016) used voxelwise amplitude of low-frequency fluctuations (ALFF) between frequencies 0.01–0.08 to find that SZs have lower ALFF compared with TCs especially in posterior cortex, occipital, and cerebellar lobes. This observation has been reported in other studies as well (Alonso-Solis et al., 2017; Calhoun et al., 2012; Hare et al., 2017; Hoptman et al., 2010). On the other hand, Alonso-Solis et al. (2017) report that ALFF for SZ is higher than that of TCs in insula. So it seems reasonable to think that the relationship between ALFF differences in SZs and TCs are different for different regions. In addition, Yu et al. (2014) found that these differences are dependent on the frequency; future studies should consider studying different frequency bands.
All these studies point to differences in the frequency profile of activity-level time series in individuals with SZ. A natural evolution of these studies is to explore the frequency profile of connectivity-level information. An earlier attempt at this utilized wavelet coherence methods (Yaesoubi et al., 2017), but as mentioned earlier, this method implemented frequency tiling in activity space; therefore, the relationship between frequency and connectivity patterns is not direct.
In the Methods section, we first introduce the proposed approach and its formulation. Next, we attempt to provide intuition into how our method performs using a toy example. After this, the dataset used in the paper is introduced briefly. We mention findings in the Results section and explore these in more detail in the Discussion section. Finally, we discuss limitations and end the paper with some concluding remarks.
MATERIALS AND METHODS
From the gx,y(t) equation we can see that it is quite similar to the equation for rx,y(t) in Equation 1. Their difference is in how the windowed mean and standard deviation is calculated. In the SWPC equation (Equation 1), for each window we have one mean and standard deviation (the index of μx is t, not i). In contrast, in gx,y(t) (Equation 6), windowed mean and standard deviation are calculated using a window around each sample (the index of μx is i here). Therefore, we can interpret the convolution between h(t) and w(t) as an approximation for SWPC. Based on the system and signal theorem (Oppenheim, 1999), we know that the output of a system with an impulse response function h(t) and input of w(t) is h(t) × w(t). So gx,y(t) (and SWPC that it approximates) is the output of a system with impulse response h(t) and input of w(t) (Figure 2). The defined ht for the SWPC system is a rectangular window that can be viewed as a low-pass filter. In other words, the output of the system (SWPC estimation) is a low-pass signal and we lose the high-frequency information of the connectivity.
A filter bank is an approach that is used frequently in the electrical engineering field (Boashash, 2015). The basic idea behind a filter bank is to design an array of systems to filter one time series into its different frequency sub-bands (usually nonoverlapping bands that cover the entire frequency spectrum).
We have designed a simple toy example to provide insight into how our method works. Note that the aim for designing this toy example was not to comprehensively evaluate FBC for use with fMRI data, but rather to convey intuition. In addition, we believe this approach can be used to analyze many types of datasets, and therefore we designed the toy example to be a general illustration.
Two scenarios with different values for fcorr and analysis parameters were designed. In both scenarios, we use two filter banks (one low-pass filter and one band-pass filter). In the first scenario, we choose SWPC window size such that it covers the same frequency band as the two designed filters. This scenario represents a case where the window size is chosen correctly for SWPC (it covers all the connectivity frequencies where the information resides). In the second scenario, the SWPC only covers the low-pass filter in the filter bank. This scenario represents a case where the window size is chosen larger than what should be used (i.e., we are filtering out some of the relevant connectivity frequencies). This scenario can happen either because of the researcher’s mistake (note that we generally do not know the ground truth about real data), or the technical limitation of the connectivity estimator paired with sliding window. For example, if sample Pearson correlation is used with sliding window (as is the case here), we know that if the number of samples is below a specific number, sample Pearson correlation estimator fails (it gives very skewed results and in the worst case estimates only 1 and −1). Therefore, there is a lower bound on window size if we use this estimator. As the frequency of a rectangular window is tied directly to its window size, the lower bound on window size causes a higher bound on frequency, that is, the connectivity is low-passed.
To compare FBC with SWPC methods, we analyzed the toy example data using both methods. To provide a more direct comparison between FBC and SWPC, the filters designed for toy examples did not cover the whole spectrum of the simulation (we cannot have a window that filters the whole frequency and reach a good estimation of correlation for SWPC). In addition, for simplicity purposes, we used the same window size for both FBC (window used to estimate sample mean and standard deviation) and SWPC. The results of both FBC and SWPC were then clustered into four clusters using the k-means approach. The reason for choosing four as the number of clusters was because we had two original states. Each state has two extreme correlation values (+At and −At because of the sinusoid nature of ρi,j), so we essentially have two pairs of states.
The fcorr value for each state and the designed filter frequency response will be shown when discussing the results.
Real Dataset and Preprocessing
To demonstrate the utilization of FBC, we used it to analyze a resting-state fMRI dataset including SZ and TC individuals. The data were obtained as a part of the Functional Imaging Biomedical Informatics Research Network (FBIRN) project (Potkin & Ford, 2009). The dataset used in this paper includes 163 TC and 151 SZ. The data acquisition and preprocessing steps are explained in our previous work (Damaraju et al., 2014). To summarize, echo planar imaging was used to acquire 162 volumes of bold data at seven sites all using 3T MRI scanners. All scans were acquired using 2 s as TR. Subjects’ eyes were closed during the scanning session.
Preprocessing was started with motion correction, slice-timing correction, and despiking. Next, data were registered to a common Montreal Neurological Institute (MNI) template and smoothed to 6 mm full width at half maximum. For the last step of preprocessing, each voxel time series was variance normalized.
To decompose the data into 100 spatially independent time series and their associated spatial maps, the pipeline proposed by Allen et al. (2014) was used. In the proposed approach group spatial independent component analysis (GICA) implemented in the GIFT (http://trendscenter.org/software/gift) software was used (Calhoun, Adali, Pearlson, & Pekar, 2001b; Erhardt, Allen, Damaraju, & Calhoun, 2011). The 162 time points for each subject were first reduced into 120 dimensions using principal component analysis (PCA). All subjects’ reduced data were then concatenated and another PCA was used to reduce the dimension to 100. Finally, independent components were estimated using the infomax algorithm (Bell & Sejnowski, 1995). ICA was repeated 20 times in ICASSO algorithm (Himberg & Hyvarinen, 2003) and the most central solution was selected for stability purposes (Du, Ma, Fu, Calhoun, & Adali, 2014). Subject-specific time series and their associated spatial maps were calculated using a back reconstruction approach (Calhoun, Adali, Pearlson, & Pekar, 2001a; Erhardt et al., 2011). The spatial maps of these 100 components were visually inspected, and 47 components were chosen as components of interest and were grouped into seven functional domains. These time series where then band-pass filtered between 0.01 Hz and 0.15 Hz using Butterworth filter (fifth order). The data used for the current project are the same 47 components used by Damaraju et al. (2014). The seven functional domains are auditory (AUD), attention/cognitive control (CC), subcortical (SC), cerebellar (CB), default mode (DM), sensorimotor (SM), and visual (VIS). The spatial maps of all the components included in each functional domain can be viewed in Figure S2 in the Supporting Information.
FBC and SWPC Analysis (Real Data)
To analyze the FBIRN dataset using FBC pipeline, each component pair (from the pool of 47 components) was used to calculate a specific wt (Equation 5) for that pair. For calculating wt, a window with size equal to 10 TR (22 s) was used. This step resulted in 1,081 wt time series (47 × (47 − 1)/2). All wt were then filtered in a forward-backward filtering system of the designed filters. For this paper, 10 IIR filters were designed to filter all the values. To obtain the optimal order for filters, we used cheb2ord as implemented in MATLAB to achieve at most 30 dB attenuation in the stopband and 3 dB in the passband (Rabiner & Gold, 1975).
The choice of the number of filters in this analysis is similar to the choice of Fourier transform length in a frequency analysis. As long as the designed filters are stable, the investigators can choose their desired number for any given analysis. We ran the analysis using different filter numbers to ensure the results are consistent. To check for filter effects, we also evaluated the performance of the pipeline using Butterworth and elliptic filters (see the Supporting Information).
The filter bands frequencies are the following:
Band 1: 0.000–0.025 Hz
Band 2: 0.025–0.050 Hz
Band 3: 0.050–0.075 Hz
Band 4: 0.075–0.100 Hz
Band 5: 0.100–0.125 Hz
Band 6: 0.125–0.150 Hz
Band 7: 0.150–0.175 Hz
Band 8: 0.175–0.200 Hz
Band 9: 0.200–0.225 Hz
Band 10: 0.225–0.250 Hz
The filtered values that resulted from all filter banks were then clustered using the k-means approach. For k-means clustering, we used the k-means++ algorithm (Arthur & Vassilvitskii, 2006) with squared Euclidean distance. The clustering was repeated 30 times with different initial cluster centroids, and the one with the lowest within-cluster distance was selected as the best result. To find the best cluster number, we used the elbow criteria paired with within-cluster distance (see Figure S3 in the Supporting Information).
Next, we calculated how many time points each subject has spent in each specific state (fraction rate). This value can be calculated for all 10 frequency sub-bands separately or for all of them combined. In summary, for each state, we produce a plot that is quite similar conceptually to the frequency response of the corresponding state. In addition, we have compared the fraction rate between TC and SZ for each state. Figure 3 illustrates the pipeline used in this paper.
Within this pipeline, there are two analysis parameters that must be selected. First, we have to choose a window size that will be used to calculate wt (i.e., calculate mean and standard deviation of component pairs). Here we have chosen a window size of 11 time points (22 s) for this step. Results from other window sizes are provided in the Supporting Information, Figures S4 through S12. As seen in these figures, the results are similar across all window sizes. Second, the number of clusters (i.e., k) needs to be selected for k-means. We chose k = 8 as the desired number of clusters in all our analyses explained in this manuscript. Our selection was based on within-cluster distance as the metric (for more details, refer to the Supporting Information; Figure S3). In addition, we also performed our analysis using different cluster numbers and provide this information in the Supporting Information (Figures S4 through S12). All statistical tests were corrected for multiple comparison. We used a method that controls false discovery rate (FDR) on the results (Benjamini & Hochberg, 1995).
As mentioned in the Methods section, we designed several toy example scenarios using a multivariate Gaussian probability density function to demonstrate the benefits of FBC compared with SWPC. In addition, we show the use of FBC on real data including SZ and TC.
Toy Example Results
All the toy example results are summarized in Figure 4. Each subfigure (boxes A through F) has seven rows. Row numbers are brought in the left-hand side of boxes A and D. The first and fifth rows demonstrate the normalized frequency response of the covariance matrix for FBC (first row) and SWPC (fifth row). The second and third rows demonstrate the mean and standard deviation of the centroids calculated from the toy example data, while the sixth and seventh rows demonstrate the same measures for SWPC. Conceptually speaking, the mean shows that the estimated values are close to the true values on average, while standard deviation represents how much variation exists in the estimated values. The fourth row shows the total fraction rate of each state in relation to the two filters. This information is exclusive to FBC and is not available for SWPC (SWPC is essentially one filter). Clusters 1 and 2 are the maximum and minimum values for state 1 (connectivity in state 1 oscillates between these values). Clusters 3 and 4 are the maximum and minimum values for state 2. True amplitude matrices, A(t), for the two states are shown in Figure S1 in the Supporting Information.
In the first scenario (Figure 4A through 4C), the SWPC window size (10 time points) covers the same frequency band as the two filters used in FBC. This scenario includes three specific situations. In the first situation, each state fcorr is located in one separate filter (Figure 4A), while in the second and third situations fcorrs are either in the first filter (Figure 4B) or in the second filter (Figure 4C). Looking at the top three rows of Figure 4, we can make several observations. The means of estimated clusters are estimated very distinctly using FBC. In other words, clusters 1 and 2 show only state 1 At while clusters 3 and 4 show only state 2 At. In contrast, if we look at SWPC results, the means of estimated clusters are estimated distinctly only when both states have low frequencies (Figure 4B). When the frequencies of the states are higher (Figure 4A and C), cluster means are not distinct. For example, cluster means of Figure 4C for SWPC show both states 1 and 2 At patterns (Figure S1) in all clusters.
In the second scenario, the window size is longer compared with the first scenario (30 time points) and only covers the frequency band of the first filter of FBC (Figures 4D through F). This scenario represents the case where SWPC window size is not chosen correctly. That is, the passband frequency of the SWPC window does not cover all frequencies where state connectivity frequencies are located. This can happen either because of a user’s incorrect choice or because of the technical limitation of the estimator that is paired with a sliding window (Pearson correlation here). This is shown in Figures 4D and 4F, where at least one of the states has connectivity frequency outside the passband of the SWPC window. Looking at SWPC results in Figures 4D, we can see that clusters 3 and 4 show very weak versions of state 2 At. In a more severe case where all states’ connectivity frequencies are outside the passband of SWPC (Figures 4F), all SWPC cluster mean values are very weak. Contrary to SWPC results, FBC mean values for cluster centroids are quite strong and similar to states 1 and 2 At. Note that the scales of the images are the same across all scenarios.
In addition to the means, we can also interpret the standard deviation values. Based on Figures 4A through C, we can see that SWPC produces higher standard deviation values compared with FBC standard deviation values for null connectivity elements (matrix entries where true connectivity is 0 in each state). In contrast, for the second scenario results, we can see that using longer window sizes results in lower standard deviation values where the true connectivity frequency is very low. This is the case for Figure 4D (first state) and Figure 4E (both states). In addition, when the specific state has a higher frequency than the SWPC frequency band, the standard deviation values are noticeably higher for all matrices’ entries. This can be seen for the second state in Figure 4D and for both state results in Figure 4F.
Another interesting observation can be made for the case where the two states’ fcorrs are located in separate filter bank frequency bands (Figures 4A and 4D). The fraction rate for FBC (row 4 in each of the boxes 4A and 4D) shows that clusters 1 and 2 spend more time in the first filter bank, while clusters 3 and 4 spend more time in the second filter. This is an accurate reflection of the ground truth and shows that FBC can provide correct frequency specificity in the case where connectivity frequencies are distributed in different filter frequency bands. This is also an exclusive and important feature of the FBC approach. This point will be expanded upon more in the Discussion section.
Figure 5 shows the correlation between the estimated clusters and the true centroids (both positive and negative matrices shown in Figure S1). Based on this figure, we can see that in most cases FBC performs better (has a higher correlation) except the case where the connectivity frequencies are well inside the SWPC frequency window (cases B and E). Thus, if even one of the states has a higher frequency than what is covered by the SWPC window, FBC performs better (cases A, C, D, and F). In addition, compared with SWPC, FBC shows more robust performance, that is, FBC performance is similar across different scenarios. This means that FBC performance is not impacted greatly by the true connectivity frequency, which is unknown to us. Finally, FBC results seem to have less variation (spread of correlation values) compared with SWPC results, which have only low variation in cases B and D. Note that these two cases are essentially the best case scenarios for SWPC.
Group Differences Between TC and SZ in the FBIRN Data
As mentioned in the Methods section, FBC was utilized to analyze the FBIRN dataset. After calculating FBC values for all component pairs, we compute eight clusters using the k-means method. The results can be seen in Figure 6. Based on the fraction rate values across all 10 bands (second row in each separate box) we can group the clusters into three groups: low-pass, band-pass, and high-pass clusters.
Figure 7 shows all the SWPC results. Figure 7A depicts the static FNC (FNC calculated over the entire time series using Pearson correlation) for the TC and SZ groups separately. Looking at Figure 7A, we can see that the TC group has a stronger positive connectivity block in AUD, VIS, and SM functional domains compared with SZ static FNC. Comparing these results with the two low-pass clusters in Figure 6 (clusters 1 and 2), we can see a resemblance between static FNCs and clusters 1 and 2. This can also be verified using Figure 7C, where correlation is used to assess the similarity of the cluster centroids. Based on this figure, TC static FNC highest correlation is with FBC cluster 2, while SZ static FNC has a higher correlation with FBC cluster 1. Another observation that supports this comparison is that TCs spend significantly more time in cluster 2 compared with SZs (Figure 6, last row in each box). In addition, SZs tend to stay more in FBC cluster 1 compared with TC.
Figure 7B shows the eight clusters that resulted from SWPC. Comparing Figure 7 with Figure 6, we can see that all the clusters resulting from SWPC are repeated in FBC results. This statement can be verified using Figure 7C. As seen in the aforementioned figure, all SWPC clusters have a high correlation with at least one of the FBC clusters. In contrast, four of the FBC clusters are not visible in SWPC results, namely, clusters 5 through 8. All of these clusters are from either the band-pass or the high-pass group and some show opposite patterns to each other, pointing to an oscillating effect in connectivity patterns.
We used FDR to correct the p values for eight comparisons (number of states) for this analysis. Looking at the comparison between TC and SZ fraction rate of FBC results (Figure 6, last row), we can see that from the eight comparisons, seven are significant after correcting for multiple comparisons. TCs tend to stay more in clusters 2, 3, 4, and 6, while SZs tend to stay more in clusters 1, 7, and 8. Lastly, if we look at the frequency profile of cluster pairs 7–8 (SZ > TC) in Figure 6 and compare them with other clusters, we see that 7–8 clusters have relatively higher frequencies compared with the other pair. This is especially more visible for cluster 8, where the higher frequencies have higher fraction rate values.
In this preliminary work, we developed a new method (FBC) to estimate connectivity dynamics that are not limited to low-frequency connectivity, or even to a specific choice of frequency. We first designed a toy example and showed that FBC enables us to estimate high-frequency changes in functional connectivity, while typical SWPC might miss these changes. We then used FBC to analyze the FBIRN dataset and found eight distinct connectivity states, each with their own unique frequency profile. We showed that FBC is able to estimate the states resulting from SWPC in addition to some other states that go undetected using SWPC because of their higher frequency profile. In addition, FBC enables us to explore the frequency profile of connectivity patterns in the whole frequency range. That is, using FBC we are able to comment on the frequency profile of each state. Applying this approach to real data reveals results in SZ that are consistent with, but extend, previous work and adds to our understanding of functional brain differences in this disorder.
Activation Frequency Versus Connectivity Frequency
As we mentioned in the Introduction section, many studies have tried to explore the frequency profile of connectivity using different methods. The majority of these methods perform frequency tiling in the activity domain and then calculate connectivity. Then they proceed to make indirect inferences about the connectivity frequency profile. This is problematic, as the relationship between the activity domain and the connectivity domain is unclear at best and depends heavily on the specific estimation method used. Many methods estimate connectivity (i.e., transform the activity domain into the connectivity domain) using highly nonlinear systems, and therefore the frequency information is distorted in this transformation. For example, looking at the SWPC formula (Equation 1) we see that μx(t) is subtracted from x(i) and the resulting value is divided by σx(t). This part in itself distorts the frequency profile of x. In addition, [x(i) − μx(t)]/σx(t) is multiplied by [y(i) − μy(t)]/σy(t) to calculate correlation (Equation 1). This step will further distort the frequency information. Therefore, using frequency information within the activity domain to infer frequency-related information specific to connectivity is not straightforward. Some studies overlook this detail when studying connectivity frequency profile. For example, Li, Bentley, and Snyder (2015, p. E2528) used frequency tiling in the activity domain to conclude that “oxygen correlation is band limited.” A similar issue can be found with a more recent paper where they talk about dynamic functional connectivity at specific bands (Luo et al., 2019). In another type of study, Yaesoubi et al. (2017) used the wavelet transform to decompose the activity time series into different time-frequency bands and then calculate coherence. Several observations about connectivity frequencies are then made. In our view, these kinds of statements can be misleading as the frequency profile is not directly studied in the connectivity space and therefore it does not enable us to make claims about the connectivity frequency profile. Rather they show that connectivity is caused by signals (i.e., activity) from specific frequency bands.
We believe that to make correct claims about the connectivity frequency profile, it is important to implement frequency tiling directly in the connectivity space. In addition, we should differentiate between the frequency profile of activity (estimated from time series themselves) and the frequency profile of the connectivity. The relationship between these two is not clear, so using knowledge about time series frequencies to make inferences regarding the connectivity frequency profile is not as straightforward as some studies suggest (Leonardi & Van De Ville, 2015).
FBC Performs Frequency Tiling in the Connectivity Domain
In this work, we have proposed an approach called FBC to estimate dFNC that does not make any assumption about the frequency profile of connectivity. This is in contrast with SWPC, which applies a low-pass filter when calculating dFNC. Note that we have filtered the BOLD time series using a band-pass filter (between 0.01 and 0.15) based on the existing literature (Niazy, Xie, Miller, Beckmann, & Smith, 2011). But to the best of our knowledge, no previous work has made the distinction between activity and connectivity frequency response, and our work is the first one to do so. Because of this, we decided not to assume any prior knowledge regarding the frequency profile of connectivity, which resulted in considering 10 bands covering the entire sampled spectrum. Additionally, we explore the spectrum of activity time series and connectivity time series before filtering, that is, w(t), and presented the results in the Supporting Information (Figure S13). As can be seen in this figure, the activity spectrum is clearly bounded between 0.01 and 0.15 Hz while the w(t) spectrum does not seem to be bounded between any two frequencies. This point should be further explored in a future work.
Using k-means clustering to summarize FBC results, we found that in addition to states that resulted from SWPC, we can estimate some other states exclusive to FBC (see Figure 6). In addition, because of the frequency tiling nature of our method, we are able to discuss the frequency profile of the estimated connectivity patterns. In many dFNC studies, it is reported that one dFNC state tends to be quite similar to the static connectivity (Damaraju et al., 2014; Faghiri et al., 2018). In this study, we see that there are actually three states that have strong low-frequency profiles (states 1, 2, and 3; see Figure 6). State 1 is quite similar to sFNC that resulted from SZ individuals, while the other state, state 2, is similar to sFNC that resulted from TCs (see Figure 7C for correlation between different states). One possible conclusion from this observation is that FBC enables us to distinguish between states that show more static-like behavior (states that show low-pass frequency profile) and others that are more dynamic (states that show a high- or band-pass frequency profile). This feature of FBC can be utilized so that we have a full picture of FNC in all its frequencies, that is, study both static and dynamic FNC patterns simultaneously. As mentioned in the last paragraph, we did not limit our connectivity spectrum, but this can easily be done. To do so, one can design filters to cover the desired bands.
sFNC Repeats in Higher Frequencies With Lower Cognitive Control Connectivity
Similar to the states reported in our earlier work (Damaraju et al., 2014), several of the states that resulted from the SWPC approach (Figure 7B) are quite similar to each other (states that show similar patterns to overall sFNC). We can see two states similar to these in FBC results, too (namely states 2 and 4). Unlike the SWPC results, here we can see that while state 2 shows a very strong low-pass frequency profile, state 4 shows a more broadband frequency distribution. Again because of the frequency tiling of FBC, we have some added information in regard to the frequency profile of connectivity patterns that we can investigate. Cluster 4 has lower CC connectivity in general. One possible observation we can make is that cluster 2 (which is quite similar to TC sFNC and has a very low-frequency profile) occurs in higher frequencies in the form of cluster 4 with smaller CC connectivity values. We speculate that CC (especially intraconnectivity within CC components) is lower in higher frequencies. This point should be examined in more depth in future work.
SZ Connectivity Patterns Have Higher Frequencies Compared With TC Connectivity Patterns
SZ subjects have significantly higher fraction rate for cluster pairs 7 and 8, while TC subjects have significantly higher fraction rate for cluster pairs 4 and 6 (Figure 6). Looking at the frequency profile of these cluster pairs, we can see that cluster pairs 7 and 8 are in the high-pass group while pairs 4 and 6 are in the band-pass group. This seems in line with some previous studies where it was reported that in activity space SZs have higher power at higher frequencies compared with TCs (Garrity et al., 2007; Alonso-Solis et al., 2017). Unlike these studies, we have found alterations in the frequency profile of connectivity instead of activity. This observation has been made possible because of the unique feature of our proposed approach that allows us to study the frequency properties of connectivity directly. We believe these points need further investigation.
Weak Connection Between Somatomotor and Visual/Auditory Networks in SZ
FBC identified two unique states, 7 and 8, which are not visible via the SWPC approach. SZ spend significantly more time in these two states compared with TC (p < 0.01, FDR corrected). Our approach enables us to intuitively capture the frequency profile and connectivity patterns of each state simultaneously, which provide information that is not available to the SWPC approach. If we compare these states with states 4 and 6 (higher fraction rate for TCs compared with SZs) we see a difference in the connection between SM and sensory areas (AUD and VIS functional domains). These connections are strong (regardless of their sign) in states 4 and 6, while they are weak in clusters 7 and 8. Several studies have reported the presence of connectivity between motor and sensory regions in typical controls (D’Ausilio et al., 2009; Londei et al., 2010). A reduction in connectivity between these regions in SZ has been reported in the past years (Berman et al., 2016; Kaufmann et al., 2015).
Dynamic Connectivity Highlights Oscillations Between Two Opposite Patterns
Another observation we can make based on the FBC results (Figure 6) is that there are pairs of states (undetected by SWPC) that show opposite connectivity patterns. In our case, states 4 and 7 show opposite connectivity patterns of states 6 and 8, respectively. One interesting insight about these is that they exhibit a more high-pass frequency profile compared with the other states. This is possibly the reason that these states are not estimated using SWPC. As mentioned in the Methods section, we used forward-backward filtering, therefore it is unlikely that these opposite patterns are due to an alteration in phase caused by the analysis. Another explanation can be that these states are spurious estimations caused by the high-frequency nature of the filter banks, as discussed by Leonardi and Van De Ville (2015). We cannot completely rule out this possibility, but because these states show the opposite patterns of some other low-pass states, this is unlikely. To further investigate, we removed the low-frequency information of SWPC with different high-pass filters (forward-backward filtering) and then performed k-means clustering. Figure 8 shows the resulting clusters using different high-pass filters. When removing the low-frequency information, we obtain these opposite clusters (e.g., C3 in Figure 8 is quite similar to FBC cluster 6). Another possible explanation for these opposite patterns is that, similar to how the toy examples were designed, the dynamic connectivity patterns oscillate between two opposite patterns. Unlike the toy examples, the negative patterns (states 2, 4, and 8) show a higher frequency profile compared with their more low-pass counterparts, but this could reflect asymmetric oscillatory behavior (e.g., more rapid return from the negative patterns). This is an interesting observation; however, future work should be designed to further explore the possibility that connectivity indeed oscillates between two opposite connectivity patterns. In addition, the biological/clinical meaning of this view should be explored further.
A Systematic View of Estimations
We believe that any analysis steps (including statistical estimators) can be viewed as a system and benefit from the extensive work done in the system design field (Bentley, Dittman, & Whitten, 2000). In this work, we examined SWPC and proposed an approximate system diagram for it (Figure 2A). This way of thinking about SWPC facilitated our conclusion that SWPC has a low-pass filter inherent to its formula and led to our proposed method. The use of filter banks as proposed in the current work is only one possibility: for future work we can use custom h(t) functions to extract specific information. For example, one possible choice is the wavelet transform in the place of subsystem 2 (wavelet can be viewed as a system itself). Another possibility is to improve SWPC by replacing its h(t) with an IIR low-pass filter with very sharp transition and flat band-pass (in contrast with SWPC with a rectangular window where the frequency response of the filter is sinc like). This is in line with previous work that has suggested that SWPC window shape can be modulated to achieve better results (Mokhtari, Akhlaghi, Simpson, Wu, & Laurienti, 2019).
Robustness of the Results in Regard to Different Parameter Choices
There are several choices we made that can impact the results. First, there is the choice of filter. We used Chebyshev type 2 filter for this analysis, but we also repeated all analyses with two other filter types (Butterworth and elliptic) with matching characteristics. Figure S14 in the Supporting Information depicts the results for all three filter types. As can be seen in this figure, almost all the cluster repeats in all three filter types. The only difference is in cluster 7, where the fraction rate for Chebyshev filter is high-pass, while the other two filters show a more band-pass fraction rate. This difference is probably caused by the difference in the frequency response of the filters. It is important to note that even though cluster 7 shows band-pass fraction rate for the last two filters, it has higher frequency profile compared with other band-pass clusters (clusters 4, 5, and 6). That is, it started going down in higher frequencies.
Another parameter was window size used for estimating w(t). In Figures S4 through S12 in the Supporting Information, we demonstrated that different window sizes with different cluster numbers show similar results. But to make sure the differences found between SZ and TC hold for different window sizes, we did all the analysis with different window sizes (from 2 TR to 60 TR) and eight clusters and performed all our statistical tests. As can be seen in Figure S15 in the Supporting Information, almost all the results for different clusters hold, with the only distinction being cluster 7. This cluster is not estimated using higher window sizes, which leads us to believe that this cluster is showing very fast changes in variance and/or mean. Both Figures S14 and S15 show that we should be careful when talking about cluster 7 of our main results. On the other hand, cluster 8, which shows an opposite pattern and forms a pair with cluster 7, is repeated in both Figures S14 and S15 with the same statistical test results (higher fraction rate in SZ). This observation reassures us that cluster 7 is a valid connectivity pattern.
Finally, we repeated the analysis with different numbers of filters, while keeping the filter type and cluster number the same as the main results of this paper. Figure S16 in the Supporting Information illustrates these results for different numbers of filters, from 5 to 12. Apart from clusters 4 and 8 (of the main results), all other six clusters are repeated in other rows with similar frequency profile. Cluster 4 is a little different in frequency profiles for filter number 5 and 8 compared with other filter numbers results. But the frequency profile, while categorized as low-pass, has a bump in middle-frequency bands (unlike the first three clusters, which are high only in very low frequencies). In addition, clusters in both columns 5 and 6 in Figure S16 are quite similar visually. Other different clusters are clusters in columns 3 and 10. The reason behind this is possibly how the tiling has been done. When using filter numbers 5–7, we do not have good representation for the higher frequency (near normalized frequency 1); therefore, the high-pass cluster is not estimated with this number of clusters. But for higher filter numbers this cluster has always been estimated (for filter numbers 8 through 12). Therefore, we think that this cluster is strong when we represent the highest frequencies more fairly. The number of filters to be used is a choice that the investigators should make and is akin to selecting the number of Fourier transform points. There are two important properties that should be considered when designing the filters. First, the filters should be stable, that is, all their poles should be outside unit circle in z-plane Oppenheim (1999). Second, the sum of the filtered time series (i.e., rn,x,y) should be very close to the unfiltered time series (i.e., w(t)), ideally equal.
Limitations and Future Directions
The first limitation of the FBC approach is the subsystem used to transform activity space to connectivity space (Figure 2). In this system, we have used a window to calculate the mean and standard deviation. The problem arises when mean and standard deviation move faster than what the window can track. In this case, our estimations will be suboptimal. To remedy this, we can use other instantaneous connectivity estimators suggested (Faghiri et al., 2020; Shine et al., 2015). Another limitation of this method is the possibility of noise contamination in the higher frequencies. This is certainly a valid issue, though our goal in this preliminary work was to refrain from making any strong assumptions about connectivity frequencies. The reason behind this decision was that we have no solid prior knowledge about connectivity frequency. In fact, to the best of our knowledge, this work is the first one that differentiates between activity and connectivity frequency profiles and previous work makes assumptions about connectivity frequency using activity space information. Future investigators can utilize our proposed approach to study any specific frequency bands more specifically.
In this study we have chosen to utilize our novel method to analyze a dataset that has been used extensively in many other studies (Damaraju et al., 2014; Hare et al., 2017; Sui et al., 2018). This can be viewed as a limiting aspect of this study, but doing so allowed us to compare our results with the previously published work by using the same data and the same function for transforming from activity to connectivity (i.e., w(t)) as the previous work (Damaraju et al., 2014). What we changed here is essentially the averaging step (Figure 2 subsystem 2) so that it does not remove frequency band information. We show that this provides a much richer source of information and additional insights into resting brain function in individuals with schizophrenia. One final limitation of this study is that, as shown in Figure 1, the estimated values at some bands might be nonsignificant. However, we have used k-means to cluster the results, which does not distinguish between significant and nonsignificant connections. To remedy this, in future work we can employ some methods that have a built-in statistical test. One such method can be to change point-detection algorithms that have been increasingly applied to fMRI data (Jeong, Pae, & Park, 2016; Xu & Lindquist, 2015). We can pair these methods with our estimation pipeline to find meaningful changes in the connectivity at different bands.
In this paper, we have used k-means paired with Euclidean distance for the summarization step. We selected this algorithm to be consistent with previous work (Damaraju et al., 2014; Yaesoubi et al., 2015). In addition, Abrol et al. (2017) found evidence for reproducibility of k-means results (paired with Euclidean distance). Exploring different clustering approaches (e.g., ensemble clustering) would be an interesting direction for future investigation.
In this work we compare FBC directly with one other estimator (i.e., SWPC). This comparison is useful because SWPC is currently the most widely used estimator, and additionally, FBC was inspired by SWPC. A comprehensive comparison of FBC with other novel connectivity estimation approaches (Pedersen et al., 2018; Shine et al., 2015) is beyond the scope of this preliminary work, which has the primary goal of introducing a new estimator. We leave a more comprehensive comparison between FBC and other estimators for future work. We also suggest that such a comparison should be done with great care to explore the estimators with regard to how they behave in the frequency domain. To the best of our knowledge this has not been done yet.
In this work, we proposed a new approach to estimate dFNC called FBC. Our proposed approach does not make any strong assumptions about connectivity frequency (unlike SWPC) and performs the frequency tiling in the connectivity domain. This is in contrast to previous work where frequency tiling was implemented in the activity domain. FBC aims to estimate connectivity in all frequencies, and it enables us to investigate connectivity pattern frequency profiles. Using toy examples, we showed that FBC is able to even estimate high-frequency connectivities in addition to providing information about the estimated connectivity frequencies. Utilizing FBC, we analyzed an fMRI dataset including TC and SZ. Using FBC we found evidence of both static connectivity and time-varying states (typically identified with SWPC) in addition to some new connectivity states undetected by SWPC (possibly because of their high-pass nature). Finally, FBC points to a possible view of connectivity in which data oscillate between two opposite connectivity patterns. This view should be further explored in future works.
DATA AND CODE AVAILABILITY
Because of limitations imposed by the IRB we are unable to share the raw data, but it is possible to share the derived results. In addition, all the code used in this study will be included in the GIFT software and also shared upon direct request.
Supporting information for this article is available at https://doi.org/10.1162/netn_a_00155.
Ashkan Faghiri: Conceptualization; Methodology; Software; Visualization; Writing - Original Draft. Armin Iraji: Writing - Review & Editing. Eswar Damaraju: Data curation. Jessica Turner: Investigation; Writing - Review & Editing. Vince Calhoun: Conceptualization; Investigation; Supervision; Writing - Review & Editing.
Vince Calhoun, Foundation for the National Institutes of Health (http://dx.doi.org/10.13039/100000009), Award ID: NIH 1 U24 RR021992. Ashkan Faghiri, Foundation for the National Institutes of Health (http://dx.doi.org/10.13039/100000009), Award ID: NIH 1 U24 RR025736-01. Ashkan Faghiri, National Institutes of Health (http://dx.doi.org/10.13039/100000009), Award ID: R01EB020407. Ashkan Faghiri, Foundation for the National Institutes of Health (http://dx.doi.org/10.13039/100000009), Award ID: P20GM103472. Ashkan Faghiri, Foundation for the National Institutes of Health (http://dx.doi.org/10.13039/100000009), Award ID: P30GM122734. Ashkan Faghiri, Office of Integrative Activities (http://dx.doi.org/10.13039/100000106), Award ID: 1539067.
Data collection was supported by the National Center for Research Resources at the National Institutes of Health (grant numbers: NIH 1 U24 RR021992, NIH 1 U24 RR025736-01, R01EB020407, P20GM103472, P30GM122734) and the National Science Foundation (1539067).
- Functional connectivity:
Statistical relationship between two or more time series, each belonging to a different brain region.
- Static connectivity:
Can be called averaged connectivity. A statistical relationship that is calculated using the whole length of the time series; that is, we have one single value for each connectivity pair.
- Dynamic connectivity:
Or time-varying connectivity. A statistical relationship that is calculated using a portion of the whole length of the time series. We have several values for each connectivity pair that changes with time.
- Filter banks:
A combination of several filters that separate a time series into several time series, each belonging to a single frequency sub-band.
- Finite impulse response filter:
A digital filter that depends linearly only on the input (signal being filtered) from previous samples.
- Infinite impulse response filter:
A digital filter that depends linearly on both the input (signal being filtered) and the output (filtered signal) from previous samples.
- Wavelet transform:
A time-frequency analysis method that allows us to explore different frequencies of a time series using different temporal/spatial resolutions. In this article, we are looking at temporal resolution.
Competing Interests: The authors have declared that no competing interests exist.
Handling Editor: Olaf Sporns