Single-scale time-dependent window-sizes in sliding-window dynamic functional connectivity analysis: A validation study.

During the past ten years, dynamic functional connectivity (FC) has been extensively studied using the sliding-window method. A fixed window-size is usually selected heuristically, since no consensus exists yet on choice of the optimal window-size. Furthermore, without a known ground-truth, the validity of the computed dynamic FC remains unclear and questionable. In this study, we computed single-scale time-dependent (SSTD) window-sizes for the sliding-window method. SSTD window-sizes were based on the frequency content at every time point of a time series and were computed without any prior information. Therefore, they were time-dependent and data-driven. Using simulated sinusoidal time series with frequency shifts, we demonstrated that SSTD window-sizes captured the time-dependent period (inverse of frequency) information at every time point. We further validated the dynamic FC values computed with SSTD window-sizes with both a classification analysis using fMRI data with a low sampling rate and a regression analysis using fMRI data with a high sampling rate. Specifically, we achieved both a higher classification accuracy in predicting cognitive impairment status in fighters and a larger explained behavioral variance in healthy young adults when using dynamic FC matrices computed with SSTD window-sizes as features, as compared to using dynamic FC matrices computed with the conventional fixed window-sizes. Overall, our study computed and validated SSTD window-sizes in the sliding-window method for dynamic FC analysis. Our results demonstrate that dynamic FC matrices computed with SSTD window-sizes can capture more temporal dynamic information related to behavior and cognitive function.

Introduction

Brain functional connectivity (FC) has been widely studied using resting-state functional magnetic resonance imaging (fMRI), and is most commonly assessed using Pearson’s correlation coefficient among distributed brain regions (Biswal et al., 1995). Many studies have identified a set of statistically interdependent regions where the blood-oxygenated-level-dependent (BOLD) signals are temporally correlated in the absence of an explicit task (Biswal et al., 2010; De Luca et al., 2005; Greicius, 2008; Greicius et al., 2003). Investigating resting-state FC has provided fundamental insight into basic neural function and disease conditions (Bai et al., 2009; Damoiseaux et al., 2006; Fox et al., 2005; Greicius et al., 2004; Smith et al., 2009).

More recently, studies have demonstrated that instead of being temporally static, FC can change periodically over time during an fMRI scan (Allen et al., 2014; Chang and Glover, 2010; Jones et al., 2012; Liu and Duyn, 2013; Majeed et al., 2011; Sakoğlu et al., 2010; Smith et al., 2012). Dynamic FC analysis has been proposed to investigate and understand these periodic temporal changes (see Calhoun et al., 2014; Hutchison et al., 2013; Preti et al., 2017 for reviews). Altered dynamic FC have been reported in various neurological disorders such as schizophrenia (Damaraju et al., 2014), Alzheimer’s disease (Jones et al., 2012), Parkinson’s disease (Cordes et al., 2018), major depression disorders (Holtzheimer and Mayberg, 2011), and autism (Price et al., 2014). All of these findings demonstrate the pathophysiologic relevance of dynamic FC across multiple diseased populations, and thus investigating dynamic FC will lead to better understanding of these conditions.

Many methods have been proposed for dynamic FC analysis, such as the sliding-window method (Allen et al., 2014; Chang and Glover, 2010), multivariate volatility models (Lindquist et al., 2014), temporal independent component analysis (ICA) (Calhoun et al., 2001; Smith et al., 2012), the quasi-periodic pattern method (Majeed et al., 2011; Thompson et al., 2014), the hidden Markov model (Eavani et al., 2013), and co-activation pattern (CAP) analysis (Liu and Duyn, 2013). The sliding-window method captures the FC in a certain time window by calculating pairwise correlation or covariance, then shifts to the next time window and repeats the same procedure until the last fMRI volume. Dynamic FC is then characterized by the windowed covariance or correlation matrices. Multivariate volatility models, such as exponentially weighted moving average and dynamic conditional correlation methods, are parametric models of the conditional correlation between time courses which refine the concept of sliding-window. Temporal ICA decomposes each subject’s fMRI time series into temporally independent components, defined as distinct functional modes. The quasi-periodic pattern method evaluates the spatiotemporal dynamics of the BOLD signals to identify a common whole-brain pattern that occurs periodically during an fMRI scan. The Hidden Markov model decodes connectivity dynamics into a temporal sequence of hidden network “states” for each subject. The CAP analysis examines the dynamics of a specific resting-state network and tracks the variations of FC within each individual time frame.

Among all these methods, the sliding-window method is one of the most widely applied techniques in dynamic FC analysis due to its relative simplicity, and it has been successfully applied to reveal macro-scale spatiotemporal fluctuation of the brain in multiple diseased populations (Damaraju et al., 2014; Gonzalez-Castillo et al., 2015; Kiviniemi et al., 2011; Leonardi and Van De Ville, 2015; Yaesoubi et al., 2017; Yu et al., 2015). However, a few concerns remain unaddressed.

The first concern is the choice of window-size. Currently, there is no consensus on the optimum length of the window used in this method. Ideally, the window-size should be both small enough to capture temporal transients and large enough to produce stable and statistically powerful results (Hutchison et al., 2013). Leonardi and Van De Ville (2015) recommend choosing a window-size that exceeds the longest wavelength in the BOLD signal to avoid spurious connectivity fluctuations (~100s), which is termed the rule-of-thumb in window-size selection. Zalesky and Breakspear (2015) later validate this rule-of-thumb choice but also demonstrate that, theoretically, the dynamic FC can be stably detected with a much shorter window (~40s). The fundamental nature of choices for these fixed window lengths limits the dynamic FC analysis to the fluctuations in the frequency range below the window period, irrespective of the true frequency content of the data (Preti et al., 2017). In this case, a time-dependent window-size, based on the true frequency content of the fMRI time series, could more accurately capture the time-varying FC changes in the data and is necessary to the sliding-window method. Furthermore, adaptive window-sizes have already been successfully applied in other fields, including frequent item sets mining in data stream (Deypir et al., 2012; Li and Wang, 2017) and estimation of data center source utilization (Baig et al., 2020).

Another concern is how to validate the observed dynamic FC as true temporal variations instead of confounding noises. Multiple null models have been proposed to test the statistical significance of the dynamic FC variations (Hindriks et al., 2016; Zalesky et al., 2014). However, heterogeneous test results could be obtained from choosing different methods to produce surrogates and applying different test-statistics (Zalesky and Breakspear, 2015). Without a known ground truth, this procedure is data dependent and time-consuming. In this case, indirect validation from investigating the useful information carried by dynamic FC could be more appropriate to determine the significance of dynamic FC analysis. Previous studies have utilized dynamic FC information as features in classifying normal subjects from diseased populations (Rashid et al., 2016; Shen et al., 2010) and predict cognitive function at the individual level (Liégeois et al., 2019; Liu et al., 2018). Therefore, a reliable classification result with dynamic FC matrices as features suggests dynamic FC matrices capture true temporal variations that carry disease-relevant information. Furthermore, multiple studies have revealed associations between static FC matrices and behavioral assessments such as age (Allen et al., 2011), gender (Zhang et al., 2018), intelligence (Song et al., 2008), cognition and emotions (Greicius et al., 2003), which demonstrate the behavioral and clinical relevance of static FC matrices. However, the associations between dynamic FC matrices computed from the sliding-window method and behavioral assessments have not been investigated yet and remain a potential tool for validating the observed dynamic FC as a true temporal variation that carries behavioral-relevant information.

Taken together, in this study, we first computed single-scale time-dependent (SSTD) window-sizes at every time point in the sliding-window approach using empirical mode decomposition (EMD). With SSTD window-sizes we refer to exactly one window-size per time point instead of multiple window-sizes per time point (as in multi-scale window-sizes). EMD is a data-adaptive analysis method for studying the naturally occurring frequency bands in time series, which can be used for non-stationary signals and allows the decomposition of time series into nearly orthogonal modes spanning narrow frequency bands (Huang, 2005). Decomposed signals from EMD track local periodic changes of non-stationary time series and time-dependent window-sizes can be determined at each time point. Multi-scale time-dependent window-sizes based on EMD were first proposed in Chen et al. (2010), where multiple window-sizes were computed for a single time point. Multi-scale time--dependent window-sizes have been successfully applied in dynamic analysis of temperature changes and dissolved oxygen time series (Huang and Schmitt, 2014). However, correlation values computed from multiple window-sizes at a single time point between two time series are difficult to interpret in terms of dynamic FC. In this case, our group first used a weighted average to combine these multi-scale window-sizes in Cordes et al. (2018). In this study, we further improved this method to produce more stable results and termed it SSTD window-sizes. We also demonstrated the reliability of SSTD window-sizes in capturing frequency information at every time-point using simulated time series.

We next computed the standard deviation matrices over all slided windows with SSTD window-sizes to quantify the dynamic FC. We validated dynamic FC matrices computed from SSTD window-sizes with both a regression analysis using Human Connectome Project (HCP) (Van Essen et al., 2013) data with a high sampling rate and a classification framework using Professional Fighters Brain Health Study (PFBHS) (Bernick et al., 2013) data with a conventional sampling rate. Dynamic FC matrices computed with SSTD window-sizes were compared with those computed with conventional fixed window-sizes in both analyses.

Methods

Computation of SSTD window-sizes

EMD decomposition

EMD uses a sifting algorithm to decompose the original multicomponent time series into nearly orthogonal modes spanning narrow frequency bands (Patrick and Goncalves, 2004), which are named intrinsic mode functions (IMFs). Let denote the resting-state time series from a single voxel and N is the total number of time points. To avoid unstable EMD decomposition results at the starting and ending point of y(t), a new time series is first created by: where flip(y(t)) reverses the order of the time points in y(t). is then intensity normalized and decomposed into K IMFs using the sifting algorithm (Patrick and Goncalves, 2004): where c(t) is the ith IMF and r(t) is the residual function. The number of IMFs (K) here is determined from the data itself during the sifting algorithm. IMFs obtained from EMD are naturally oscillating functions and with frequency contents in a descending order, i.e. c1(t) contains the highest frequency component and c(t) contains the lowest. The residual r(t) contains only the trend-level frequency component.

Instantaneous period and average energy

The instantaneous period is computed only for the center part of each IMF c(t), which corresponds to the time points in the original signal. For each IMF, the instantaneous frequency f(t) is first computed using the Hilbert Transform by extending c(t) into the complex domain (Huang, 2005; Huang and Shen, 2014; Huang and Wu, 2008). The instantaneous period p(t) is then defined as . Furthermore, since has been intensity normalized, the energy of each IMF is defined as . In this case, the instantaneous periods p(t) capture the local non-stationarity of the original signal, and the average energy densities E summarize the energy contributions of each IMF to the original signal.

SSTD window-sizes of two fMRI time series y(1) (t) and y(2)(t)

Following the above steps, we have computed and to denote the instantaneous period at time t for IMFs corresponding to y(1) (t) and corresponding to y(2) (t), respectively. We also have computed and for the average energy corresponding to each and . We then define a weighted instantaneous period for each y( (t), (m = 1, 2) as

SSTD window-size of y(1)(t) and y(2)(t) is finally defined as:

Furthermore, for noisy time series with limited durations, IMFs with low frequencies may not be sufficiently covered by the length of time series, and EMD may lead to non-orthogonal components at very low frequencies. These IMFs usually contain longer instantaneous period but lower energy. Therefore, in computing the weighted instantaneous period for each y((t), (m = 1, 2) using Eq. (3), IMFs with energy less than of the total energy have been excluded.

Validation with simulation

We demonstrate through simulation that SSTD window-sizes can represent the true period (frequency) information at every time point. To this end, we considered two pure sinusoidal signals. The first signal was a cosine wave with frequency f1 suddenly switching to f2 at a certain time point: where TR was the repetition time (inverse of the sampling rate) and N was the total number of time points. The other signal was a simple cosine wave with a fixed frequency f3:

In our simulation, we chose f1 = 0.02 Hz, f2 = 0.05 Hz and f3 = 0.04 Hz, respectively. The total duration of time series was set to 300s (5mins), and the frequency shift of y(1) happened at t = 200s. We chose the TR to be 0.72s and 2.8s, respectively, to mimic the two real fMRI data used in this study (described in more details in section 2.3 and 2.4). N was set to be in each case. SSTD window-sizes were determined as described in section 2.1. Dynamic FC values between y(1) and y(2) were computed using the sliding-window approach with SSTD window-sizes and fixed window-sizes of 30s, 60s and 100s.

Validation with classification analysis using fMRI data with low sampling rate

Subjects

65 cognitively non-impaired fighters and 68 cognitively impaired fighters from the PFBHS (Bernick et al., 2013) were included in our study. Their cognitive impairment status was predefined based on fighters’ performance of the Finger Tapping test and Symbol Digit Coding test from the CNS Vital Signs tests (Gualtieri and Johnson, 2006). Cognitively non-impaired and impaired fighters were matched for age, gender, years of education, and fighting exposures (number of fights, years of fighting and knock-out counts). Detailed demographics are listed in Table 1(A).

Table 1 (A)

Demographics of PFBHS subjects.

	Nonimpaired fighters	Impaired fighters	Group difference (p-value)
No. of Subjects	65	68	NA
Gender	58 Men	65 Men	0.16
	7 Women	3 Women
Age at Imaging (years)	28.78 ± 5.27	29.78 ± 6.20	0.32
Years of Education (years)	13.28 ± 1.63	13.03 ± 2.12	0.45
Processing Speed score	58.28 ± 7.30	40.85 ± 8.34	<0.001
Psychomotor Speed score	183.12 ± 15.95	153.16 ± 16.11	<0.001
Number of Fights	14.45 ± 12.97	14.47 ± 12.68	0.99
Years of Fighting	6.03 ± 4.02	6.85 ± 4.45	0.27
Knock-outs	0.78 ± 1.14	1.07 ± 1.83	0.27
fMRI motion (mm)	0.23 ± 0.09	0.25 ± 0.11	0.30

MRI data were collected for all subjects on a 3T Siemens Verio scanner with a 32-channel head coil. Resting-state fMRI data were collected with following parameters: TR = 2.8s, TE = 28 ms, in-plane resolution 2mm × 2mm, slice thickness 4 mm, 30 axial slices and 137 time frames (6mins and 24s). In addition, a high resolution T1-weighted structural image was acquired using a standard 3D MPRAGE sequence.

Preprocessing

Each T1-weighted image was input into the FreeSurfer 6.0 processing pipeline (Fischl, 2012) to generate a subject-specific cortical and subcortical parcellation. 66 cortical labels from Desikan-Killiany atlas (Desikan et al., 2006), and 12 sub-cortical labels were defined as regions of interest (ROIs). This anatomical labeling was wrapped to each subject’s fMRI space using a 12-parameters affine transform in Advanced Normalization Tools (ANTs) software (http://stnava.github.io/ANTs/).

The first 4 time frames (~12s) of fMRI data were removed to allow the MR signal to achieve T1 equilibrium. Remaining time frames were slice-timing corrected and realigned to the mean echo-planar image in SPM12 (http://www.fil.ion.ucl.ac.uk/spm/). FMRI data were further spatially smoothed using a 6 mm 3D-Gaussian filter in the subject fMRI space and high-pass filtered with a cosine filter with a cut-off frequency of 0.008 Hz. Nuisance regression was performed to further denoise the fMRI data, using six head motion parameters and CompCor generated white-matter (WM) and cerebrospinal fluid (CSF) signals (Behzadi et al., 2007). Finally, all voxel time courses were variance normalized.

Dynamic FC estimation

Average time series of each ROI were obtained in the subject’s fMRI space. Static FC were estimated by computing the Pearson’s correlation coefficient between time series from each ROI pair. Dynamic FC were estimated using the sliding-window method with both SSTD window-sizes and multiple conventional fixed window-sizes. Pearson’s correlation between time series from all ROI pairs were computed within the defined window and the window was slided over 1 TR (2.8s). The standard deviation of Pearson’s correlations over all windows were defined as the dynamic FC values for all pairs of ROIs and formed the dynamic FC matrices. This dynamic FC matrices then were used to represent whole-brain temporal dynamics.

Classification framework

Both static and dynamic FC were used as features to classify cognitive impairment status of fighters. Classification framework included an automated feature selection step and a radial basis functional classifier step (see Mishra et al., 2017, for details). A ten-fold cross-validation was used to determine the classification accuracy. The classification process was repeated 100 times to rule out the potential division bias in the ten-fold division. Classification accuracy, sensitivity (correctly classified impaired fighters), specificity (correctly classified nonimpaired fighters), and the area under the receiver operating characteristic curve (ROC) for all 100 iterations were used to evaluate the classifier performance.

Validation with regression analysis using fMRI data with a high sampling rate

The same subjects used in our previous study (Yang et al., 2019) were used in this study. Briefly, 88 male subjects from HCP 1200 Subject Release (Van Essen et al., 2013; WU - Minn Consortium Human Connectome Project, 2017), who were 26 to 30 years old and completed both resting-state fMRI and T1-weighted structural scans, were included. Detailed demographic information is listed in Table 1 (B).

Table 1 (B)

Demographics of HCP subjects.

	Normal Subject
No. of Subjects	88
Gender	88 Men
Age at Imaging (years)	27.71 ± 1.28
Years of Education (years)	15.03 ± 1.65
fMRI motion	Passed HCP quality control

Resting-state fMRI data were acquired on 3T Siemens scanners using a Gradient-echo EPI sequence with the following parameters: multiband factor = 8, TR = 0.72s, TE = 33.1 ms, 2.0 mm isotropic voxels, 72 axial slices and 1200 time points (14mins and 24s). A high resolution T1-weighted structural image was acquired using a 3D MPRAGE sequence.

Behavioral measurements

56 behavioral measurements were used to evaluate if the dynamic FC matrices computed with the proposed optimum window-size can better explain the behavioral variance than the dynamic FC matrices estimated with conventional fixed window-sizes. Initially, 64 behavioral measurements covering cognitive, social, emotion, personality traits, tobacco consumption, drug consumption, rule breaking, antisocial behavior, pain intensity, and externalization were included. Out of these 64 behavioral measurements, 8 were excluded from our analysis based on the following two criteria: 1) more than 5% of subjects did not have valid measurements; 2) the maximal value was more than 100 times the mean value, which suggested extreme outliers. Details of these 56 behavioral variables are listed in Supplementary Table S1.

The minimally preprocessed resting-state fMRI data (in the standard MNI space) were downloaded from the HCP website (https://www.humanconnectome.org/study/hcp-young-adult/document/1200-subjects-data-release) and were treated as the raw fMRI data in our analysis. A more detailed description of HCP minimal preprocessing steps can be found in Glasser et al. (2013). Additionally, the first 15 vol (~11s) of fMRI data were discarded to avoid data with an unsaturated T1 signal, and the same nuisance regressors as used for the PFBHS data were used to further denoise the HCP data.

Average time series from 112 regions of the automated anatomical labeling atlas (Tzourio-Mazoyer et al., 2002) were obtained. Static and dynamic FC matrices were estimated the same way as for the PFBHS data. Specifically, both SSTD window-sizes and multiple conventional fixed windows were used in the sliding-window approach. The window was slided over 4 time points (2.88s) due to the fast sampling rate.

Regression analysis

Static and dynamic FC matrices were rearranged as a single vector for each subject, and these vectors for all subjects were arranged in a matrix. Principal component analysis (PCA) was applied on the matrix to reduce data dimension and only those components that explain more than 1% of the variance were retained. Linear regression was applied for each behavioral measurement separately with PCA components as independent variables and the behavioral measure as the dependent variable. The coefficient of determination, namely R-squared, was used to measure the proportion of behavioral variance explained by the FC matrix (Draper and Smith, 1998). Subjects’ age and handedness were used as covariates and regressed out from behavioral measurements prior to PCA.

Results

Simulated signals with TR = 0.72s and 2.8s give similar results, and therefore we include simulation results with TR = 0.72s here in Fig. 1. The simulation results with TR = 2.8s is included in supplementary Fig. S1. Simulated signals y(1) and y(2) are plotted in Fig. 1 (A). Corresponding frequency spectrums are shown in Fig. 1 (B). For time points before t = 200s, y(1) is dominated by frequency f1 = 0.02 Hz and y(2) is dominated by frequency f3 = 0.04 Hz, and therefore the corresponding periods of y(1) is 50s (~70 TR) and of y(2) is 25s (~35 TR) in this time segment. For time points after t = 200s, y(1) is dominated by frequency f2 = 0.05 Hz and y(2) is dominated by frequency f3 = 0.04 Hz, and therefore the corresponding periods of y(1) is 20s (~28 TR) and of y(2) is 25s (~35 TR) in this time segment. As shown by the solid blue (y(1)) and red lines (y(2)) in Fig. 1 (D), the instantaneous periods computed from our method capture exactly these periods at almost all time points for both signals. The SSTD window-sizes are then set to be the larger period of the two signals, as shown by the dashed green lines in Fig. 1 (D). Multiple fixed window-sizes are also shown as dashed yellow, purple and light blue lines in Fig. 1 (D) for references. Fig. 1 (C) plots dynamic FC values computed from the sliding-window approach with SSTD window-sizes (solid green line) and multiple fixed window-sizes (dashed yellow, purple and light blue lines). As shown by the solid green line, using SSTD window-sizes in the sliding-window approach captures both the stable and small correlation between two signals at time points before t = 200s, and large variations of correlations between two signals at time points after t = 200s.

Fig. 1.

Simulation results of TR = 0.72s. (A). Simulated time series y(1) (blue) and y(2) (red). (B). Corresponding frequency spectrums. (C). Static FC values between two time series (dashed pink line). Dynamic FC values computed using the sliding-window method with SSTD window-sizes (solid green line) and multiple fixed window-sizes (dashed yellow, purple and light blue lines). (D). Instantaneous periods of y(1) (solid blue line) and y(2) (solid red line), SSTD window-sizes (dashed green line), and multiple fixed window-sizes (dashed yellow, purple and light blue lines).

Validation with classification analysis

Average SSTD window-sizes

Average SSTD window-sizes computed from our method are 34.66s±2.29s for cognitively nonimpaired fighters and 34.20s± 2.29s for cognitively impaired fighters (Table 2). Histograms of average SSTD window-sizes for each subject are also shown in Fig. S2 (B). Thus, the main fixed window-size for comparison is chosen to be 35s (~12 TR). We further compute the dynamic FC matrices with fixed window-sizes of 60s (~21 TR) and 90s (~32 TR) as suggested in Hutchison et al. (2013) for comprehensive comparisons.

Table 2

Repetition time (TR) and average SSTD window-sizes computed for both HCP data and PFBHS data.

	HCP data	PFBHS data
		Nonimpaired fighters	Impaired fighters
TR (s)	0.72	2.8
Duration	14mins and 24s	6mins and 24s
Average SSTD window- sizes (s)	32.73 ± 3.52	34.66 ± 2.29	34.20 ± 2.29

Classification results

Fig. 2 shows the classifier performance of the PFBHS data, with static FC (sFC, red), static and dynamic FC with fixed window-size of 35s (sFC&dFC-35s, green), and static and dynamic FC with SSTD window-sized (sFC&dFC-SSTD, purple) as features, respectively. Boxplots of the area under the ROC curve (A), accuracy (B), sensitivity (C), and specificity (D) are shown for 100 iterations. As shown in Fig. 2, adding temporal dynamics information significantly improves the classifier performance, as compared to using static FC information alone (green and purple boxes vs. red box). Using dynamic FC matrices computed with SSTD window-sizes further significantly improves the classifier performance, as compared to using dynamic FC matrices computed from the comparable fixed window-size (p = 1.52 × 10−10, purple box vs. green box).

Fig. 2.

Classification results of the PFBHS data. Area under the ROC curve (A), accuracy (B), sensitivity (C) and specificity (D) of classification between cognitively nonimpaired and impaired fighters, when using static FC alone as features (red), using both static FC and dynamic FC computed with fixed window-size (35s) as features (green), and using both static FC and dynamic FC computed from SSTD window-sizes as features (purple). Boxplots show measurements of 100 iterations. Abbreviations: sFC: static functional connectivity; dFC-fixed-35s: dynamic functional connectivity matrix computed with the 35s window-size; dFC-SSTD: dynamic functional connectivity matrix computed with single-scale time-dependent window-sizes.

Fig. 3 further plots areas under the ROC curves for classifications using static FC and dynamic FC with 35s (sFC&dFC-35s), 60s (sFC&dFC-60s), 90s (sFC&dFC-90s), and SSTD windows (sFC&dFC-SSTD) as features separately. Table 3 lists p-values of statistical comparisons (two-sample t-test) between classification results using different features. As shown in Fig. 3 and Table 3, using sFC&dFC-SSTD as features significantly outperforms sFC&dFC-35s (p = 1.52 × 10−10), sFC&dFC-60s (p = 1.33 × 10−11), and sFC&dFC-90s (p = 3.04 × 10−12) in classifying cognitive impairment status in fighters.

Fig. 3.

Classification results of the PFBHS data. Areas under the ROC curves of classification between cognitively nonimpaired and impaired fighters, when using both static FC and dynamic FC computed with SSTD window-sizes as features (purple box) and using both static FC and dynamic FC computed with multiple fixed window-sizes (green boxes). Boxplots show measurements of 100 iterations. * indicates statistically significant differences. Abbreviations: sFC: static functional connectivity matrices; dFC-fixed-35/60/90s: dynamic functional connectivity matrix computed with the 35/60/90s window-size (~12/21/32 TR); dFC-SSTD: dynamic functional connectivity matrix computed with single-scale time-dependent window-sizes.

Table 3

Classification results of the PFBHS data: statistical significances (p-values) of different dynamic FC matrices as features in the classification analysis. Abbreviations: sFC: static functional connectivity matrices; dFC-fixed-35/60/90s: dynamic functional connectivity matrix computed with the 35/60/90s window-size (~12/21/32 TR); dFC-SSTD: dynamic functional connectivity matrix computed with single-scale time-dependent window-sizes.

	sFC&dFC-35s	sFC&dFC-60s	sFC&dFC-90s	sFC&dFC-SSTD
sFC&dFC-35s	NA	0.55	0.43	1.52E-10
sFC&dFC-60s	0.55	NA	0.86	1.33E-11
sFC&dFC-90s	0.43	0.86	NA	3.04E-12
sFC&dFC-SSTD	1.52E-10	1.33E-11	3.04E-12	NA

Validation with regression analysis

As listed in Table 2, the average SSTD window-sizes computed from the proposed method for HCP data is 32.73 ± 3.52s. Histograms of average SSTD window-sizes for each subject are also shown in Fig. S2 (A). In this case, the main fixed window-size for comparison is chosen to be 32.4s (~45 TR). Furthermore, since HCP data are acquired with a fast sampling rate, we can use a smaller window-size but still keep enough time points to compute the windowed FC matrix. Therefore, the dynamic FC matrices computed using the sliding-window approach with fixed window-sizes of 21s (~30 TR) and 65s (~90 TR) are further used for comprehensive comparisons.

Fig. 4 (A) plots the explained behavioral variance of HCP data using sFC, sFC&dFC-32s, and sFC&dFC-SSTD as features. As shown in Fig. 4 (A), significant larger variance is explained by both static and dynamic FC matrices than by static FC matrices alone (green and purple boxes vs. red box). Using dynamic FC matrices computed with SSTD window-sizes further increases the explained behavioral variance by 19.7% (purple box vs. green box) on average. Fig. 4 (B) further compares the performance of SSTD window-sizes with multiple fixed window-sizes in the regression analysis. As shown by the purple box, larger behavioral variances can be explained more by the dynamic FC matrices computed with SSTD window-sizes than by multiple fixed window-sizes (green boxes). Statistical significance levels (two-sample t-test) are listed in Table 4.

Fig. 4.

Regression analysis results of the HCP data. Boxplots of explained behavioral variance, when using static FC matrices as features alone (red box in A), both static and dynamic FC matrices computed with SSTD window-sizes as features (purple boxes in A and B), and both static and dynamic FC computed with multiple fixed window-sizes as features (green boxes in A and B). Abbreviations: sFC: static functional connectivity matrices; dFC-fixed-21/32/60s: dynamic functional connectivity matrix computed with the 21/32/60s window-size (~30/45/90 TR); dFC-SSTD: dynamic functional connectivity matrix computed with single-scale time-dependent window-sizes.

Table 4

Regression results of the HCP data: statistical significances (p-values) of explained behavioral variances when using different dynamic FC matrices as features. Abbreviations: sFC: static functional connectivity matrices; dFC-fixed-21/32/60s: dynamic functional connectivity matrix computed with the 21/32/60s window-size (~30/45/90 TR); dFC-SSTD: dynamic functional connectivity matrix computed with single-scale time-dependent window-sizes.

	sFC&dFC-21s	sFC&dFC-32s	sFC&dFC-60s	sFC&dFC-SSTD
sFC&dFC-21s	NA	0.39	0.01	1.18E-08
sFC&dFC-32s	0.39	NA	0.10	4.93E-07
sFC&dFC-60s	0.01	0.10	NA	5.15E-04
sFC&dFC-SSTD	1.18E-08	4.93E-07	5.15E-04	NA

Discussion

In this study, we computed SSTD window-sizes in the sliding-window method for dynamic FC analysis. SSTD window-sizes of two time series were based on their frequency contents and were computed from the instantaneous period and average energy of each IMF obtained from the EMD method. Using simulated time series, we first demonstrated that SSTD window-sizes can capture the time-dependent frequency information in the time series. The optimum performances of dynamic FC matrices computed with SSTD window-sizes were further validated through both a classification analysis using the PFBHS data with a low sampling rate and a regression analysis using the HCP data with a high sampling rate. Specifically, both a higher classification accuracy in predicting cognitive impairment status in fighters and a larger explained behavioral variance in healthy young adults were achieved using dynamic FC matrices computed with SSTD window-sizes as features, as compared to using dynamic FC matrices computed with several conventional fixed window-sizes.

SSTD window-sizes are data-driven, time-dependent, and can capture the frequency content in simulated time series

SSTD window-sizes are based on the frequency content of the time series itself, and thus are purely data-driven. The frequency content of the time series is represented by the instantaneous period at every time point and therefore SSTD window-sizes are time-dependent. Simulated time series with known frequency contents are used to demonstrate that time-dependent and data-driven window-sizes can accurately capture the frequency information of the time series.

The simulated time-series are basic oscillatory and stationary functions with local means equal to zero, therefore the Hilbert Transform is well defined and produces meaningful non-negative instantaneous frequencies for these two time series (Huang and Shen, 2014). As shown in Fig. 1 (A), in our simulation, a sudden frequency switch from f1 = 0.02 Hz to f2 = 0.05 Hz happens at t = 200s for signal y(1) (Fig. 1 (A)), and the frequency for signal y(2) stays at f3 = 0.025 Hz for the entire time series. The lowest frequency contained in simulated signals is therefore 0.02 Hz. According to the rule-of-thumb criteria (Leonardi and Van De Ville, 2015), the window-size used in the sliding-window approach should be larger than . However, dynamic FC values computed with choices following this rule miss the large periodic temporal variations after t = 200s and lead to a larger transition period (dashed light blue and purple lines in Fig. 1 (C)). Furthermore, dynamic FC values computed with a relatively small window are able to capture the large temporal variations after t = 200s, but create unstable and unreal large temporal variations before t = 200s (dashed yellow line in Fig. 1(C)). In some cases, dynamic FC values computed with a small window-size may even pass the statistical significance test, but the strong variability is still not real and is mainly the reason for the lack of degree of freedom in the calculation (Chen et al., 2010).

SSTD window-sizes, however, automatically capture the frequency switch in our simulation by outputting a window-size of 50s before t = 200s and another window-size of 25s after t = 200s (dashed green line in Fig. 1 (D)). This result demonstrates that SSTD window-sizes are time-dependent and capture the exact frequency information of the time-series. Furthermore, no prior information is required to calculate SSTD window-sizes and thus the procedure is purely data-driven. Dynamic FC values computed from SSTD window-sizes also give optimum performance, as they not only capture both large temporal variations in functional correlation between y(1) and y(2) after t = 200s but also correctly identify small correlation between y(1) and y(2) before t = 200s (solid green line in Fig. 1 (C)).

The simulation with TR = 2.8s (Fig. S1) gives similar results as the simulation with TR = 0.72s (Fig. 1), but the less time points of TR = 2.8s results in less smoothed curves in Fig. S1. Overall, as shown in Fig. 1 and Fig. S1, our simulation results demonstrate that SSTD window-sizes can capture frequency changes in time series.

SSTD window-sizes capture more meaningful temporal dynamics information of fMRI time series

Functional MRI time series are usually non-stationary. Therefore, in our method, we first decompose the original time series into multiple IMFs, which represent different scales of the original time series and form the adaptive and oscillating basis of the data (Chen et al., 2010; Huang, 2005; Huang et al., 1998). An IMF is defined as a basic oscillatory function with an extreme value followed by a zero-crossing (Patrick and Goncalves, 2004). Thus, the Hilbert Transform is well defined for an IMF (Huang and Shen, 2014). Furthermore, since the fMRI time series is first intensity normalized in our study and EMD is essentially a sifting algorithm that estimates an IMF and subtracts it successively (Eq. (2)), the average energy of each IMF represents weight and contribution of each IMF to the original signal in terms of intensity. Therefore, at each time point, we then use the average of instantaneous periods of all IMFs weighted by the corresponding energy (Eq. (3)) to represent the period information of the original fMRI time series. We further demonstrate that the dynamic FC matrices computed using this energy-weighted period measurement can capture more cognitively and behaviorally meaningful information.

Both static and dynamic FC matrices contain disease and cognitive relevant information

Previous studies have demonstrated that FC matrices, both static and dynamic, can be used as feature sets to classify normal subjects from subjects with a cognitive disorder or disease. For example, using static FC matrices alone, Shen et al. (2010) achieved an accuracy of 92.3% in classifying schizophrenia patients from normal subjects. However, the small sample size (~50 in total) in Shen et al. (2010) limited its generalization. Rashid et al. (2016) used a much larger sample size (273) and demonstrated that using both static and dynamic FC as features, an accuracy level of 89% was achieved in classifying schizophrenia patients, bipolar patients, and normal subjects, which was significantly higher than using static FC matrices as features alone (59%). These studies have demonstrated that both static and dynamic FC contain disorder-related information. Furthermore, using static FC derived measurements, multiple studies have shown that high accuracy can be achieved in classifying normal subjects, subjects with mild cognitive impairment (MCI), and subjects with Alzheimer’s disease (Challis et al., 2015; Chen et al., 2016; Khazaee et al., 2017; Wang et al., 2018). These studies further suggest that FC matrices contain cognitively meaningful information, as MCI and AD patients usually show significant loss in cognition (see Jack et al., 2013, for a review).

Both static and dynamic FC matrices contain behavioral meaningful information

Disease is not the only factor driving the structure involved in strength and flexibility of brain FC. Multiple studies have demonstrated that a battery of personal behaviors are encoded in functional brain connections as well, including intelligence (Liu et al., 2018; Song et al., 2008), personality, socialization, cognition and emotions (Finn et al., 2015; Greicius et al., 2003; Liégeois et al., 2019). Furthermore, Liu et al. (2018) have shown that fluid intelligence is significantly correlated with the predicted score based on the strength of dynamic FC matrices, and the observed cognitive flexibility and executive inhibition are significantly correlated with the predicted scores based on the variability of dynamic FC matrices. Therefore, all of these studies suggest that FC matrices contain behaviorally meaningful information.

Dynamic FC matrices computed with SSTD window-sizes capture cognitively and behaviorally meaningful information

As shown by Figs. 2–4 and Tables 3 and 4, better classification performances and higher explained behavioral variances are achieved when using dynamic FC matrices computed with SSTD window-sizes as features in the classification and regression analysis separately, as compared to using dynamic FC matrices computed with conventional fixed window-sizes. The frequency information of a non-stationary fMRI signal remains unstable and, therefore, a fixed window-size that is optimal for one time segment may not be suitable for other time segments. SSTD window-sizes, in contrast, are data-driven and time-dependent, allowing these subtle variations of the temporal dynamics to be captured.

Technical perspectives

Improvements over previous methods

As we introduced above, SSTD window-sizes are based on the multi-scale time-dependent intrinsic correlation method that was first proposed by Chen et al. (2010), where multiple window-sizes were computed for a single time point. Specifically, at a single time point, one window-size was determined for each IMF and a 2D heat plot was generated to represent the multi-scale intrinsic correlations between two time series (Chen et al., 2010). This method captures the time-dependent frequency information in the original time series but is difficult to interpret for BOLD signals in terms of FCs.

In Cordes et al. (2018), our group proposed a weighted average to combine these multiple window-sizes at each time point in the sliding-window dynamic FC analysis. It has been observed that the Hilbert Transform produces unstable results at cutting time points, and therefore extreme and oscillated window-sizes were obtained for the starting and ending time points in Cordes et al. (2018), leading to an average window-sizes of 97.75s ± 41.36s (Cordes et al., 2018). In the current study, we first improved the method to compute SSTD window-sizes by creating longer time series to ensure the continuity of the starting and ending time points and to avoid unstable and extreme results. In addition, from the method perspective, the number of IMFs in Cordes et al. (2018) was a fixed and predefined parameter. In the current method, the number of IMFs is determined from the data itself, and therefore the entire process in computing SSTD window-sizes is data-driven without any predefined parameters. As a result, average SSTD window-sizes computed for both data in the current study is around 34s, which is closer to values reported in Zalesky and Breakspear (2015). Furthermore, from a validation perspective, the combined window-size in Cordes et al. (2018) was only validated for BOLD signals with a TR of 2.4s using the dynamic states routine following Allen et al. (2014), where the number of dynamic states was another tuning parameter which may affect the results (discussed in more detail in the next section). In the current study, we validate SSTD window-sizes using fMRI data with TR = 2.8 and 0.72s respectively, which demonstrates the optimum performance of SSTD window-sizes over fixed window-sizes in the sliding-window analysis for BOLD signals with both low and high sampling rates. Finally, the current study quantifies dynamic FC matrices using the standard deviation of each FC over all slided windows, which is computed without any tuning parameters and is a more direct measure to compare performances between different window-sizes, as compared to the dynamic states routine.

Methods to estimate instantaneous frequency

In computing SSTD window-sizes, we used EMD combined with the Hilbert transform to estimate the instantaneous frequency of BOLD signals, which is the Hilbert Huang Transform (Huang and Shen, 2014). Analytically, the instantaneous frequency of a time signal is defined as the derivative of the phase (Cohen, 1989). For a real signal, the most common method to compute its phase is to build a complex signal with the original signal as the real part and the Hilbert Transform of the original signal as the imaginary part (Zhao et al., 2004). However, this method works only for the signal with a narrow frequency range, which is not the case for multicomponent BOLD signals. The Hilbert Huang transform works for the BOLD signals by applying EMD to decompose the original signal into IMFs, where each IMF represents a simple oscillatory mode within a limited frequency band that meets the pre-request of the analytical method. Instead of EMD, narrow-band pass filters could also be applied to decompose the original BOLD signals into multiple time-series with limited frequency bands that are suitable for the Hilbert transform.

Other signal decomposition methods

As we stated above, EMD decomposes original BOLD signals into multiple IMFs with different frequency components. Other methods also can be applied to represent BOLD signals in terms of multiple frequency components, such as the Fourier transform and the Wavelet transform (Chang and Glover, 2010). Both Fourier and Wavelet transforms require a predefined fixed basis function and give a bounded time-frequency relationship, which work well for time series with known underlying physical processes. However, underlying process of BOLD signals are usually difficult to define. On the other hand, EMD requires no a priori defined basis functions and uses data itself to determine the underlying frequency information. Therefore, in our study, we used EMD, instead of Fourier or wavelet transform, to decompose BOLD signals.

Using standard deviation to quantify temporal dynamics over all slided windows

Following computing windowed FC matrices, most studies perform a clustering analysis to identify different brain-states based on spatial similarities among all windowed FC matrices, and used the clustered brain-states to quantify temporal dynamics (Allen et al., 2014; Damaraju et al., 2014; Sakoğlu et al., 2010). We do not use the classified brain states, mainly because, in this method, the number of brain-states is another parameter to be pre-determined. Functional connections within each brain-state, time points spent in each brain-state and frequency of brain-state transitions significantly depend on this predefined number. Furthermore, in a classification framework with dynamic FC as input features similar to the framework used in our study, the number of features increases linearly with the pre-defined number of brain-states, which could easily exceed the number of subjects and thus make it difficult for classification framework to produce meaningful results. In addition, Choe et al. (2017) have demonstrated that the fluctuation of dynamic FC (i.e. its variance) has a strong potential to provide individual-specific differences. Finally, Liu et al. (2018) compute dynamic FC strength, stability, and variability to characterize temporal dynamics over all slided windows and show that these measurements could successfully identify individuals with high accuracy and reliably predict higher individual cognitive performance. In our study, similar ideas are applied, and the temporal dynamics of FC between an ROI pair is quantified by the standard deviation of this connection over all slided windows.

Using both static FC and dynamic FC matrices as dynamic FC features in classification and regression analysis

In our study, we always kept both static and dynamic FC matrices as dynamic features in both classification and regression analysis. Liu et al. (2018) defined the average of FC matrices over all slided windows as dynamic FC mean strength and considered it as one of the dynamic characteristic measurements. The same measurement was computed in our study for both PFBHS data and HCP data. High spatial similarities between the static FC matrix and dynamic FC mean strength matrix were observed for both data (Supplementary Fig. S3). Therefore, in our study, a static FC matrix was used to represent the dynamic FC mean strength, and was used together with the standard deviation matrix as features for both classification and regression analysis.

Using two datasets to validate SSTD window-sizes

We used two fMRI data to validate SSTD window-sizes in this study, one with a conventional sampling rate (TR = 2.8s, PFBHS data) and a total duration of 6mins and 20s; and the other one with a faster sampling rate (TR = 0.72s, HCP data) and a total duration of 14mins and 20s. For fMRI data with a low sampling rate, temporal changes may be missed due to under-sampling of the fMRI signal, whereas for fMRI data with a high sampling rate, the captured temporal changes may be sensitive to spurious noises. Therefore, we have validated SSTD window-sizes in both cases in the current study. The larger explained behavioral variance and higher classification accuracy when using SSTD window-sizes suggest that SSTD window-sizes were able to capture temporal dynamics in fMRI data with both low and high sampling rates.

Limitations and future directions

One limitation of the SSTD window-sizes method is that the performance may deteriorate as the noise level in BOLD signals increases. The first step of the SSTD window-sizes method is to decompose BOLD signals into IMFs using the EMD method. Despite the capability of dealing with nonlinear and non-stationary data, EMD still cannot resolve the signal from the noise in complicated cases, especially when the noise has the same time scale as the signal (Wu and Huang, 2004). Future performance comparisons of the SSTD window-sizes method on BOLD signals with different de-noising strategies are needed to test and clarify the noise effects of the SSTD window-sizes method.

Furthermore, EMD is currently run on the entire concatenated time series; therefore, the number of IMFs is the same for all time points and remains static. For longer fMRI time series, the number of IMFs may also vary at different time points. One way to achieve a dynamic number of IMFs is to run EMD separately on each time series segment, which will increase the computational cost and is also less stable with shorter segments. After the EMD step, instantaneous period at every time point and average energy across all time points are used to compute the SSTD window-sizes in the current method. In this case, while the period is dynamically calculated at every time point, the energy is still assumed to be static across all time points. Static assumptions on both the energy and the number of IMFs may not be the most appropriate for a dynamic analysis. Future efforts towards computing an accurate time-dependent energy measure and determining dynamic number of IMFs at each time point, if possible, may further improve the performance of SSTD window-sizes. However, whether these extra dynamics reflect the intrinsic temporal variations of human brain function or just random noises also require further evaluations.

Another limitation of the current study is our use of simple sinusoidal signals with frequency shifts in our simulation to demonstrate that SSTD window-sizes capture the exact period (frequency) information at every time point (sections 2.2 and 3.1). Real fMRI signals are much more complicated than sinusoidal signals due to mixed frequency contents, physiological and electronic noise sources. We refrained from using more elaborate simulations with mixed frequency content because the ground-truth of the signal becomes difficult, if not impossible, to define. Nevertheless, future efforts toward generating simulated signals that both have known frequency contents and represent signal distributions mimicking more real fMRI data may provide alternatives to validate SSTD window-sizes. In addition, previous studies have suggested that the variance of dynamic FC values derived from sliding-window approach with fixed window-sizes are less reliable than the dynamic conditional correlation derived measures (Choe et al., 2017). Future efforts towards the reproducibility analysis of sliding-window approach with SSTD window-sizes and the comparisons with other dynamic FC analysis methods are required to further validate our method.

Conclusion

Single-scale time-dependent window-sizes (SSTD) have been computed in the sliding-window method for dynamic FC analysis. SSTD window-sizes are based on the frequency contents of the time series and computed with EMD and a spectral analysis using the Hilbert Huang Transform. Simulation has demonstrated that SSTD window-sizes are data-driven, time-dependent and capture the exact frequency content in the time series. Using dynamic FC matrices computed with SSTD window-sizes as features, both a higher classification accuracy between cognitively impaired and nonimpaired fighters and a larger explained behavioral variance in normal healthy adults are achieved, as compared to using dynamic FC computed with multiple conventional fixed window-sizes as features. Both classification and regression results have demonstrated that SSTD window-sizes are able to capture cognitive-function-related and behaviorally meaningful temporal dynamics information in resting-state fMRI data. MATLAB codes used to compute SSTD window-sizes and generate the simulation in this manuscript are available on GitHub (https://github.com/EEHULULU/dynamicFC_optimum_windowsize.git).