Aberrant global and local dynamic properties in schizophrenia with instantaneous phase method based on Hilbert transform

Abstract

Background

Emerging functional imaging studies suggest that schizophrenia is associated with aberrant spatiotemporal interaction which may result in aberrant global and local dynamic properties.

Methods

We investigated the dynamic functional connectivity (FC) by using instantaneous phase method based on Hilbert transform to detect abnormal spatiotemporal interaction in schizophrenia. Based on resting-state functional magnetic resonance imaging, two independent datasets were included, with 114 subjects from COBRE [51 schizophrenia patients (SZ) and 63 healthy controls (HCs)] and 96 from OpenfMRI (36 SZ and 60 HCs). Phase differences and instantaneous coupling matrices were firstly calculated at all time points by extracting instantaneous parameters. Global [global synchrony and intertemporal closeness (ITC)] and local dynamic features [strength of FC (sFC) and variability of FC (vFC)] were compared between two groups. Support vector machine (SVM) was used to estimate the ability to discriminate two groups by using all aberrant features.

Results

We found SZ had lower global synchrony and ITC than HCs on both datasets. Furthermore, SZ had a significant decrease in sFC but an increase in vFC, which were mainly located at prefrontal cortex, anterior cingulate cortex, temporal cortex and visual cortex or temporal cortex and hippocampus, forming significant dynamic subnetworks. SVM analysis revealed a high degree of balanced accuracy (85.75%) on the basis of all aberrant dynamic features.

Conclusions

SZ has worse overall spatiotemporal stability and extensive FC subnetwork lesions compared to HCs, which to some extent elucidates the pathophysiological mechanism of schizophrenia, providing insight into time-variation properties of patients with other mental illnesses.

Introduction

Although researchers are still far from understanding the underlying cause of schizophrenia, increasing researches on brain imaging consistently suggest that schizophrenia is a disease closely associated with alterations of brain structure and function (Fitzsimmons, Kubicki, & Shenton, 2013; Friston, Brown, Siemerkus, & Stephan, 2016; Skudlarski et al., 2010). During the last two decades, a hypothesis of functional disconnectivity for schizophrenia has attracted intense attention (Friston, 1998; Friston & Frith, 1995). Researches have revealed abnormalities of functional connectivity (FC) between brain regions or functional integration across brain networks in schizophrenia (Baker et al., 2014; Calhoun, Eichele, & Pearlson, 2009; Du et al., 2015; Lynall et al., 2010; Repovs, Csernansky, & Barch, 2011; Rotarska-Jagiela et al., 2010). However, most of these studies plausibly assume that inter-regional interactions are temporally stationary, and the brain FC is estimated by computing an average of the full-time series and generate a static value to reflect the connection strength (Fornito, Zalesky, Pantelis, & Bullmore, 2012; van den Heuvel & Fornito, 2014).

In recent years, there has been much interest in computing time-resolved connectivity measures and successful applications in identifying biomarkers from dynamic connectivity for mental disorders (Allen et al., 2014; Calhoun, Miller, Pearlson, & Adalı, 2014; Zalesky, Fornito, Cocchi, Gollo, & Breakspear, 2014). In such analysis, brain FC can vary within a short period rather than be considered as static over time. Such results tend to further expand the available information, and avoid the strong assumption that brain activity is static over time. Importantly, the brain dynamic functional architecture has been closely related with a wide range of cognitive and affective processes such as learning (Bassett et al., 2011), executive cognition (Braun et al., 2015), psychological resilience (Long et al., 2019), and emotion (Betzel, Satterthwaite, Gold, & Bassett, 2017), as well as multiple common psychiatric and neurological disorders such as schizophrenia (Long et al., 2020b), autism (Harlalka, Bapi, Vinod, & Roy, 2019), Alzheimer's disease (Schumacher et al., 2019), and major depressive disorder (MDD) (Long et al., 2020a; Wise et al., 2017).

There are numerous methodologies used to estimate dynamic functional connectivity (dFC). Among them, the sliding-window technique is the most widely used approach (Hindriks et al., 2016; Sakoğlu et al., 2010; Shakil, Lee, & Keilholz, 2016). For functional magnetic resonance imaging (fMRI), by assessing FC in different time-windows, one can easily expand existing static connectivity strategies to be time-resolved, and the dFC can be evaluated by measuring FC among regions of interest (ROIs) or voxels in a sliding-window yielding multiple connectivity matrices. Furthermore, the idea of dynamic connectivity based on the segmentation of windows has already been practiced for a while in the electroencephalogram (EEG) community. For instance, Lehmann used a spatial approach of adaptively segmenting EEG map series into time segments of variable length and stationary spatial characteristics to investigate the functional states of brains in detail (Lehmann, Ozaki, & Pal, 1987). Moreover, Khanna et al. also summarized a method of EEG microstate analysis which has been applied suitably to study resting-state EEG and assess global functional states of the brain in health and disease in their review (Khanna, Pascual-Leone, Michel, & Farzan, 2015). Although these approaches have been extensively applied to estimate the dFC, a considerable controversy has arisen over its lacking standards for setting the window length. If the window length is too short, the time points in each window could be too few to generate a robust estimation of connectivity strengths. In contrast, long window length might decrease the temporal variations of FC, consequently hindering from effectively detecting connectivity states (Du et al., 2018b).

Several windowless methods have been proposed to avoid the problem in selecting the window length. Bayesian approach (Robinson, Atlas, & Wager, 2015; Taghia et al., 2017) has been employed to investigate dynamic connectivity, which regards extracting time-varying functional networks as selecting dynamic models in the Bayesian setting. Variable parameter regression model combined with the Kalman filtering method is used to detect the dynamic interactions between different functional networks at all time points (Kang et al., 2011). The recently proposed time-frequency analysis (Yaesoubi, Allen, Miller, & Calhoun, 2015) explored the connectivity by using multiple frequencies, which can be conceptually seen as adapting the observation window to the frequency content of the original time courses.

In this paper, the instantaneous phase method based on Hilbert transform is used to investigate the aberrant global and local dynamic properties in schizophrenia. In the communication system, using Hilbert transform to describe the envelope of amplitude modulation, phase modulation, instantaneous frequency, and instantaneous phase of the signal will make the signal analysis easier and more effective. Moreover, instantaneous phase synchronization of resting-state fMRI based on Hilbert transform, a windowless method with a simple calculation, can be used as a measure for the dFC (Glerean, Salmi, Lahnakoski, Jääskeläinen, & Sams, 2012) and is successfully used to reveal altered variability of FC (vFC) among patients with MDD (Demirtaş et al., 2016). The aims of this research are as follows: (1) explore the overall spatiotemporal stability in schizophrenia patients (SZ) from the perspective of the whole-brain; (2) locate specific subnetworks or brain regions with aberrant dFC; (3) evaluate the ability of global and local dynamic features to differentiate SZ from healthy individuals.

Materials and methods

Participants

This study consists of two independent datasets. Dataset 1 comes from a publicly available dataset (the Center for Biomedical Research Excellence, COBRE), including 72 SZ diagnosed by the Diagnostic and Statistical Manual of Mental Disorders (Fourth Edition) (DSM-IV) and 75 healthy controls (HCs). Dataset 2 comes from OpenfMRI (UCLA Consortium for Neuropsychiatric Phenomics LA5c Study), in which a total of 50 SZ and 75 HCs were initially recruited by community advertisements from the Los Angeles area and completed extensive neuropsychological testing. The detailed inclusion and exclusion criteria of both datasets were described in online Supplementary material.

Data acquisition and preprocessing

Dataset 1 (COBRE dataset) was acquired with a 3-Tesla Siemens Trio scanner (Siemens, Germany) (5 min). All resting-state fMRI data were collected using a single-shot full k-space echo-planar imaging with ramp sampling correction (TR/TE = 2000/29 ms; flip angle = 75 degree; matrix = 64 × 64; slice number = 33; voxel size = 3 × 3 × 4 mm³). Of note, the head motion can introduce substantial alters in the time courses of resting-state fMRI data (Van Dijk, Sabuncu, & Buckner, 2012; Yan et al., 2013; Zeng et al., 2014); therefore, in addition to the head motion correction via Friston 24-parameter model (Friston, Williams, Howard, Frackowiak, & Turner, 1996), totally 20 SZ and 12 HCs were further discarded with a significant head motion from the dataset [excluding criteria: displacement >2 mm, rotation >2 degree, or mean frame-wise displacement (mFD) > 0.5 mm] (Power, Barnes, Snyder, Schlaggar, & Petersen, 2012), and one SZ was excluded from the analysis due to the wrong scanning time points. A total of 114 subjects were finally included (51 SZ and 63 HCs).

Dataset 2 (OpenfMRI dataset) was acquired on one of two 3 T Siemens Trio scanners (304 s). The scan parameters associated with resting-state fMRI are as follows: TR/TE = 2000/30 ms; flip angle = 90 degree; FOV = 192 mm; matrix = 64 × 64; slice thickness = 4 mm; slice number = 34, oblique slice orientation. Using the same head motion excluding criteria as dataset 1 (14 SZ and 15 HCs were excluded), totally 96 subjects were included in the final analysis (36 SZ and 60 HCs). There was no significant difference between the two groups regarding gender, age, and handedness for two datasets. The detailed demographics of datasets were list in Table 1.

Table 1. Subject demographics

Note: For MRI studies, subjects with left-handedness were excluded from the OpenfMRI dataset. M, denotes male; F, denotes female; L, denotes left-handedness; R, denotes right-handedness; and A, denotes ambidextrous.

The Data Processing Assistant for Resting-State fMRI toolbox (DPARSF, http://rfmri.org/dpabi) was used to carry out data preprocessing (Yan, Wang, Zuo, & Zang, 2016), which included the removal of first 10 volumes, slice timing correction, head motion correction, normalized to Montreal Neurological Institute (MNI) space and further resampled to 3 × 3 × 3 mm³, spatially smoothed by convolution with an isotropic Gaussian kernel (FWHM = 4 mm), linear detrending, nuisance signal regression (including Friston 24 head motion parameters, white matter, and cerebrospinal fluid signals), and bandpass filtering. The bandpass range was selected as 0.04–0.07 Hz, which could reduce the influence of low-frequency drift and high-frequency noise and make Hilbert transform play a better performance (see online Supplementary material for the detailed study of bandpass). Next, the brain was divided into 90 ROIs by using the Automated Anatomical Labeling atlas (Tzourio-Mazoyer et al., 2002), thus the average BOLD signal in each brain region was obtained respectively.

Dynamic functional connectivity construction

The instantaneous phase method based on Hilbert transform was used to assess the dFC in this study, which allowed us to define and extract dynamic connectivity coefficient at each time point. The Hilbert transform, S(t) = Acos((φ(t))) of the preprocessed BOLD signals broke the signal down to an analytical signal S(t) with an instantaneous phase φ(t) and amplitude A, where t indicates the time point. Phase differences Δφ_ij(t) at each time point were then calculated between ROI i and ROI j (1 ⩽ i, j ⩽ 90). The following formula was used to adjust the phase difference between 0 and π:

$$\Delta \varphi _{ij}( t{\rm \;) } = \left\{{\matrix{ {\vert {\varphi_i( t) -\varphi_j( t) } \vert {\rm \;, \;\;\;\;\;\;\;\;\;\;if\;\;0\;} \le \vert {\varphi_i( t) -\varphi_j( t) } \vert {\rm < \;}\pi } \cr { 2\pi -\vert {\varphi_i( t) -\varphi_j( t) } \vert {\rm \;, \;\;\;if\;\;}\pi {\rm \;} \le \vert {\varphi_i( t) -\varphi_j( t) } \vert {\rm < \;2}\pi } \cr } } \right.$$

Then, instantaneous coupling matrices (ICMs) C(t) were constructed using the phase difference normalized between 0 and 1:

$$C_{ij}{\rm ( t) } = 1-\Delta \varphi _{ij}( t{\rm ) \;/\;}\pi , \;1\,\le \,i, \;j \le 90, \;t = 1-140$$

The closer Δφ_ij(t) is to π (approximately perfect anti-synchronization), the closer the value of C_ij(t) is to 0, whereas the closer the value of C_ij(t) is to 1, the closer Δφ_ij(t) is to 0 (approximately perfect synchronization). Thereby, the 90 × 90 matrix C(t) could be used to denote the dynamic correlation coefficient at time point t for each subject.

Global dynamic features

Global synchrony

ICMs were binarized subsequently for each subject. When the difference between instantaneous phases of BOLD signals is small, signals can be considered highly synchronized at that time point. Therefore, a series of different thresholds were compared for the criterion of high phase synchronization. Finally, π /8 was chosen as the threshold to obtain the binarized matrix $C_{ij}^b ( t )$, where i and j are indexes of each ROI (1 ⩽ i, j ⩽ 90), and t indicates time points. The percentage of non-zero connections at each binary matrix ICM $C_{ij}^b ( t )$ was then defined as global synchrony G(t), and average global synchrony ($\overline {G( t) }$) was finally computed for each subject when considering all scanning time points. In our study, the method of instantaneous phase synchrony between brain regions measured the phase similarity or synchronization between the BOLD time series signals of two brain regions (see online Supplementary material for more details).

Intertemporal closeness

In order to quantify the temporal state stability of the whole-brain dFC, the measure of intertemporal closeness (ITC) was defined. The Pearson's correlation coefficient was used to measure the similarity between ICMs at different time points. A mean matrix of all ICMs was firstly obtained (denoted as av-FC), then we got 140 coefficients between av-FC and ICMs of each time point, respectively. Finally, 140 values were averaged to obtain a reference for each subject. Therefore, we defined ITC as the proportion of observing greater similarity between ICMs than the reference value of each subject. The larger the value of this proportion, the higher the overall similarity between ICMs at different time points. Moreover, due to the high correlation between neighboring ICMs (online Supplementary Figs S1 and S2), ITC with different time-lags (1–20 s) was also considered to eliminate the effect of adjacent time points. In order to illuminate this situation more clearly, the ITC level of 0.05 was taken as an example to characterize the time-lags differences between two groups in our study. Finally, ITC without time-lag and ITC with the given time-lags were extracted, respectively. ITC considers the similarity between each temporal state characterized by ICMs. In such a way, lower ITC indicates that the fluctuations in temporal states of subject's whole-brain dFC are higher (the dFC of whole-brain over time is less stable), while higher ITC indicates that the dFC of whole-brain is more stable over time (see online Supplementary material for more details about the definition of ITC).

Local dynamic features

Strength and variability of functional connectivity

We further investigated the local alterations of the whole-brain dFC. The strength of FC (sFC) between any pair of brain regions was defined by averaging the ICM at all time points, in other words, the sFC was defined as the mean matrix of ICMs at all time points.

The vFC was defined to describe the dispersion degree or fluctuation of the dFC (i.e. sample variance/sample mean) at all time points between any pair of brain regions. Similarly, the 90 × 90 matrix of the vFC could be obtained for each subject (see online Supplementary material for the relevant mathematical formulas about the definition of sFC and vFC).

Primary analysis

The statistical analyses included in this study are as follows: (1) For global features ($\overline {G( t) }$ and ITC), the non-parametric permutation test was performed to test significant differences between SZ and HCs (10 000 permutations, p < 0.05). (2) For local features (sFC and vFC), the Network-Based Statistics (NBS) method (Zalesky, Fornito, & Bullmore, 2010) was performed to explore subnetworks with significant differences between SZ and HCs. In this context, a component in connected graph (i.e. a subnetwork) was a set of supra-threshold connections for which a path could be found between any two nodes. All connected components could be efficiently identified with a breadth or depth search. For each pair, a non-parametric permutation test was performed 5000 times. The corrected p value for the connected subnetwork was established as p < 0.05, and the significant dynamic subnetworks were visualized through BrainNet toolbox (http://www.nitrc.org/projects/bnv/) (Xia, Wang, & He, 2013). (3) In order to further estimate the classification ability of global and local features extracted by our instantaneous phase method between two groups (SZ v. HC), the support vector machine (SVM) analysis was used in our study as an auxiliary analysis (please see online Supplementary material for more details about the SVM modeling).

Repeatability and robustness analysis

To further validate the robustness of the instantaneous phase method in our study, the statistical analysis processes in two datasets were implemented identically. Initial analysis was performed using the COBRE dataset as a discovery dataset, then analysis was repeated in the OpenfMRI dataset as a replication dataset. Moreover, we compared the existing results with the results obtained by the sliding-window approach frequently applied in prior functional dynamics studies. The analysis flow chart of this study was shown in Fig. 1.

Fig. 1. The dynamic analysis flow chart of this study. We performed Hilbert transform on original blood oxygen level-dependent (BOLD) signals and obtained instantaneous phases after narrow band filtering. Instantaneous phase differences between all brain regions were obtained, and the value of differences was finally normalized between 0 and 1. The instantaneous coupling matrix (ICM) was defined to extract global dynamic features (average global synchrony and intertemporal closeness) and local dynamic features [strength of FC (sFC) and variability of FC (vFC)] of all subjects. These features were finally used for support vector machine (SVM) classification.

Results

Global dynamic features

Two global dynamic features were defined to explore and quantify the overall spatiotemporal stability of the dFC. We firstly compared the distributions of average global synchrony between SZ and HC using the Kolmogorov–Smirnov distance between cumulative distribution functions. However, no significant difference was found in the cumulative distribution function of average global synchrony between two groups (Kolmogorov–Smirnov test, D = 0.20635, p = 0.153, Fig. 2a). Then we compared the average global synchrony between two groups, as shown in Fig. 2b, the SZ group showed a significantly lower mean of average global synchrony than did the HC group (p = 0.0107, permutation test).

Fig. 2. The results of the global and local dynamic properties in the COBRE dataset. (a) The comparison of cumulative distribution function between schizophrenia patients (SZ) and healthy controls (HC) (Kolmogorov–Smirnov test: p = 0.153); (b) comparison of mean of average global synchrony (10 000 permutations, p = 0.0107); (c) intertemporal closeness (ITC) with different time-lags in two groups (excluding nearby τ time points, 1–20 s), ITC value below the level of 0.05 was about 5.5196 s for HC group, and about 4.9588 s for SZ group; (d) comparison of mean ITC without time-lag (10 000 permutations, p = 0.0294); (e, f): dynamic brain networks analysis of significantly decreased connections in the strength of functional connectivity (sFC) and significantly increased connections in the variability of functional connectivity (vFC) (SZ v. HC). The results are based on NBS using 5000 permutations, p value < 0.05 and maximum component threshold t > 4.0 (sFC) and t > 3.4 (vFC).

Subsequently, ITC was defined to quantify the stability of the dFC over time. As shown in Fig. 2d, the SZ group showed a significantly lower mean of ITC than did the HC group (p = 0.0294, permutation test). On the other hand, in the SZ group, the ITC value fell below the level of 0.05 at a lag of approximately 4.9588 s, while for the HC group, the time lag was approximately 5.5196 s (Fig. 2c). And significant differences were also found in the mean of time-lag between two groups (p = 0.0219, permutation test). Furthermore, ITC in the SZ group was always lower than the HC group under the same time-lags condition, suggesting the greater stability of the dFC in the HC group.

Local dynamic features

The NBS method was used to investigate the local alterations in the whole-brain dFC. Compared to the HC group, the SZ group showed significantly decreased sFC and increased vFC widespread across the brain dynamic subnetworks (Fig. 2e, f, Table 2), suggesting that the dFC of SZ is characterized by an abnormal pattern of decreased strength and increased variability (5000 permutations). In addition, more than 80% of regions in significant subnetworks of sFC and vFC were consistently located at the anterior cingulate cortex (ACC), temporal cortex (including the temporal pole), frontal cortex, fusiform gyrus, cuneus and the supplementary motor area (SMA), suggesting that for SZ, decreased strength and increased variability of the dynamic FC occur in similar brain regions.

Table 2. Local alteration of brain networks for COBRE dataset (schizophrenia patients v. healthy controls)

Note: L, denotes the left cerebral hemisphere; R, denotes the right cerebral hemisphere.

SVM classification

In order to meticulously estimate the classification ability of global and local dynamic features between SZ and HCs, SVM with different types of dynamic features at global and local levels was constructed. Under the condition of narrow filter range (0.04–0.07 Hz), SVM models revealed 85.75% of balanced accuracy (Brodersen, Ong, Stephan, & Buhmann, 2010) when classifying SZ and HC by using all aberrant global and local dynamic features. The relatively lower degree of balanced accuracy (79.75%) when only using aberrant local features (sFC and vFC) was presented as input (see online Supplementary Table S2). The other SVM results were presented in online Supplementary material.

Repeatability and robustness analysis

For the replication dataset, significant differences were found in the cumulative distribution function of average global synchrony between two groups (Kolmogorov–Smirnov test, D = 0.32779, p = 0.0125, Fig. 3a). Furthermore, the SZ group also showed a significantly lower mean of average global synchrony than did the HC group (p = 0.0005, permutation test, Fig. 3b). Subsequently, SZ group showed a significantly lower mean of ITC than did HC group (permutation test, p = 0.0016, Fig. 3d). In SZ group, the ITC value fell below the level of 0.05 at a lag of approximately 5.1284 s, while for HCs, the time lag was approximately 5.4791 s (Fig. 3c). Furthermore, although significant differences were not found in the mean of time-lag between the two groups (p = 0.1655, permutation test), ITC in SZ group was always lower than the HC group under the same time-lags condition. Therefore, we still suggested that the overall spatiotemporal stability of SZ group was lower than that of HC group, which was consistent with the discovery dataset. Therefore, it should be considered that two independent datasets consistently showed an abnormal pattern with less stable dynamic functional architecture in SZ at the global level.

Fig. 3. The results of the global and local dynamic properties in the OpenfMRI dataset. (a) The comparison of cumulative distribution function between schizophrenia patients (SZ) and healthy controls (HC) (Kolmogorov–Smirnov test: p = 0.0125); (b) comparison of mean of average global synchrony (10 000 permutations, p = 0.0005); (c) intertemporal closeness (ITC) with different time-lags in two groups (excluding nearby τ time points, 1–20 s), ITC value below the level of 0.05 was about 5.4791 s for HC group, and about 5.1284 s for SZ group; (d) comparison of mean ITC without time-lag (10 000 permutations, p = 0.0016); (e, f): dynamic brain networks analysis of significantly decreased connections in the strength of functional connectivity (sFC) and significantly increased connections in the variability of functional connectivity (vFC) (SZ v. HC). The results are based on NBS using 5000 permutations, p value < 0.05 and maximum component threshold t > 3.55 (sFC) and t > 3.3 (vFC).

Subsequently, the NBS analysis on the dFC consistently found significantly decreased sFC and increased vFC in SZ group (Fig. 3e, f, online Supplementary Table S1). Specifically, more than 55% of regions in significant subnetworks of sFC and vFC were consistently located at the hippocampus, parahippocampal gyrus, temporal cortex (including the temporal pole) and the rolandic operculum, while alterations of the prefrontal cortex and ACC were only found in sFC and vFC, respectively. In conclusion, although specific regions with significant differences (sFC and vFC) were different in the replication sample from those in the discovery sample, the abnormal pattern of the dFC with decreased strength and increased variability was consistent across two datasets.

To test the robustness of the results, we further used the sliding-window approach to repeat the analysis of the dFC. Consistently, the sFC was significantly decreased in the SZ group compared to the HC group, although the alterations of vFC were not found possibly due to the low temporal resolution and analytical sensitivity of this approach (see online Supplemental material and Fig. S4, Table S3 for more details about the results derived from the sliding-window approach).

Discussion

In this study, the instantaneous phase method based on Hilbert transform was used to investigate the overall spatiotemporal stability and local alteration of the brain dynamic functional architecture in SZ. The instantaneous phases of BOLD signals were extracted and an ICM was provided for each time point, thus dynamic interactions at global and local levels were finally investigated. Our results firstly demonstrated decreased global synchrony and temporal stability in SZ at the global level; secondly, we further found that the abnormal pattern of local alterations in the dFC was characterized by significantly decreased strength and increased variability; thirdly, our replication dataset reinforced the robustness of the findings, which also showed decreased global synchrony and less stable dFC at the global level, along with the decreased strength and increased variability of dFC at the local level in SZ; finally, machine learning models as auxiliary analysis revealed a high degree of accuracy when classifying SZ and HC by using all dynamic features (average global synchrony and ITC at global levels, significantly decreased strength and increased variability at local levels) as classification features. Our study extended the disconnectivity investigation in schizophrenia from static to dynamic functional architecture, which may help to unravel the pathophysiologic mechanisms of this severe mental disorder.

The dynamic properties were investigated in schizophrenia by using the instantaneous phase method, which could adequately consider the fluctuation or vFC over time. This method has advantages of high temporal resolution, simple calculation, and less time consumption. The dynamic analysis procedures in this approach allow us to explore global and local properties of brains in SZ, which are far different from the mainstream method such as graph theory analysis (Du et al., 2016; Sun, Collinson, Suckling, & Sim, 2019) and provide new insights into the aberrant dynamic interaction in schizophrenia. Although the sliding-window approach also detected the aberrant sFC to some extent in this study, no aberrant subnetwork in vFC was found possibly due to its low temporal resolution, which may support the advantages of the instantaneous phase method on the exploration of the dynamic functional architecture of the brain.

For global features, the average instantaneous phase difference between BOLD signals was greater in SZ, and the overall phase coupling was lower, indicating the lower global synchrony in SZ during the scanning. Similarly, the ITC was significantly lower in SZ, indicating worse state stability or more fluctuations of the dFC over time, which was consistent with the results of Sun in the global properties (higher temporal global efficiency represents more fluctuations of dynamic brain network backbone in schizophrenia) (Sun et al., 2019). Moreover, in the case of ITC with different time-lags, SZ was always lower than HCs under the same condition. Therefore, our study based on two independent datasets suggests that the overall spatiotemporal stability in schizophrenia is significantly decreased.

It is noteworthy that the indicator of global synchrony in our study only considered the phase similarity (synchronization), but not the case of inverse coupling (anti-synchronization), which may lead to the indicator's insensitivity to the synchronization property of inverse coupling. It is also worth mentioning that, in comparison with inter-regional high phase synchrony method used in the current study, although the common cross-site phase coherence in the EEG measure (Reiser, Wascher, Rinkenauer, & Arnau, 2020) that similarly uses time-point specific phase angle data known as inter-site phase coherence (i.e. coherence analysis) is very similar with the correlation analysis in principle, coherence analysis can only define and detect the overall correlation between two signals in the frequency domain instead of exploring any correlations between two signals at different time points. It focuses not on the specific phase angle of the oscillation at a specific time point at each of the two locations, but rather on the phase angle difference between the locations over time, which suggest that the global synchrony used in our study may be more suitable for exploring the overall spatiotemporal stability of dFC over time.

For local features, our findings repeatedly highlighted an abnormal pattern of the dFC with decreased strength and increased variability in SZ. What should be noted is that the local alterations were mainly located at the prefrontal cortex, ACC, temporal cortex, and the visual cortex, no matter what the strength or variability of the connections. Similarly, the OpenfMRI dataset showed that the abnormal pattern of the dFC was consistent with the COBRE dataset, although the main regions with disconnectivity were different. Meanwhile, the aberrant subnetwork of the sFC by using the sliding-window approach was also identified and located at the ACC and temporal cortex. Moreover, regions with disconnectivity found in SZ were largely overlapped with the regional findings of Sun in schizophrenia-related significant increment of temporal regional efficiency (Sun et al., 2019). Interestingly, our finding was in line with the meta-analyses of structural imaging studies in schizophrenia, suggesting that the most robust gray matter changes occurring in the frontal and temporal regions may parallel with those with aberrant dynamic functional properties (Ellison-Wright, Glahn, Laird, Thelen, & Bullmore, 2008; Fornito, Yücel, Patti, Wood, & Pantelis, 2009). The lesions in the prefrontal regions were prone to mental disorders, and investigations related to the prefrontal of the brain were not rare (Duncan & Owen, 2000; Pucak, Levitt, Lund, & Lewis, 1996). FFG was thought to be critical for face recognition, which may be associated with impaired facial recognition and interpretation of facial expression in schizophrenia. Researches showed that schizophrenia was associated with bilateral reduction in FFG gray matter volume (Lee et al., 2002), and neuroanatomic FFG abnormalities underlie at least some of the deficits associated with facial recognition (Onitsuka et al., 2003). Heuristically, the disconnectivity associated with the ACC was predominantly found in SZ, which suggested that this area may be a key region in the functional characterization of schizophrenia. However, significant abnormalities of the posterior cingulate were consistently discovered in most brain function studies related to schizophrenia (Fletcher, McKenna, Friston, Frith, & Dolan, 1999; Garrity et al., 2007). Thus, while the interpretation of case-control differences in resting-state fMRI may not always be straightforward (Fornito & Bullmore, 2010), aberrant dynamic subnetworks identified in frontal, temporal, cingulate, or hippocampus regions using our method in this study were consistent with pathological alterations reported in the literature.

Moreover, models including dynamic analytical approaches or measures have demonstrated a sound relationship with the psychopathology of schizophrenia. Results from Kottaram at the resolution of whole-brain networks showed that the severity of positive symptoms was associated with a longer proportion of time spent in states characterized by inactive default mode and executive networks (Kottaram et al., 2019). Dynamic connectivity analysis of auditory hallucinations found that connection-wise variability in SZ was reduced between the left auditory perception and speech-production brain areas (Zhang et al., 2018). Dynamic graph analysis indicated that reduced stability of lateral occipital cortex connectivity may be an important factor underlying neurocognitive dysfunctions and symptom severity in schizophrenia (Li, Sweeney, & Hu, 2020). Meanwhile, the dFC also shows functional specificities in participants considered to be at risk of developing schizophrenia (Barber, Lindquist, DeRosse, & Karlsgodt, 2018; Bolton et al., 2020; Briend, Armstrong, Kraguljac, Keilhloz, & Lahti, 2020; Du et al., 2018a). For instance, using innovation-driven coactivation patterns (iCAPs) for dynamic large-scale brain network analysis, Zöller et al. uncovered the patterns of brain network activation duration and coupling that were relevant with clinical risk for psychosis in 22q11.2 deletion syndrome, demonstrating that longer durations and couplings of iCAPs were associated with the severity of positive psychotic symptoms and anxiety, thus reinforcing the case for relationships between decreased dFC and schizophrenia (Zöller et al., 2019).

Our investigation was subject to various limitations. Firstly, due to the absence of more detailed clinical and demographic information (many clinical variables were not available in the public dataset), we were unable to further explore the relevance between clinical or behavioral information and our findings in global or local features in our study; however, such information and the corresponding analysis would still be important to evaluate the clinical utility of this approach in future studies. Secondly, Hilbert transform is only applicable to relatively narrow bandpass signals, and the narrower the bandwidth of the target signal, the better does the Hilbert transform produces an analytic signal with meaningful envelope and phase (Bedrosian, 1963); thus, the performance will be relatively reduced for some non-narrow bandpass signals. Finally, our research was subject to the short scan duration (scanning time of two datasets are 300 and 304 s respectively), which might be insufficient for the stabilization of the dFC in the resting-state and bias the results. Therefore, subjects with relatively longer scan duration could be investigated to improve the reliability of the dFC in future experiments.

In conclusion, quantitative assessment of the dFC in terms of the global and local perspectives, as performed in this study, provided new insights into which to better our understanding of the anomalous pattern of brain function during resting-state in schizophrenia. Our study showed that beyond a significant decrease in the overall spatiotemporal stability, there were local aberrations of strength and variability of the dFC in schizophrenia. Further investigation of the relationship between the aberrant dFC networks and physiological signification may provide more clinical value with which to deepen comprehension in the pathophysiology of schizophrenia.