Temporal-spatial dynamic functional connectivity analysis in schizophrenia classification

With the development of resting-state functional magnetic resonance imaging (rs-fMRI) technology, the functional connectivity network (FCN) which reflects the statistical similarity of temporal activity between brain regions has shown promising results for the identification of neuropsychiatric disorders. Alteration in FCN is believed to have the potential to locate biomarkers for classifying or predicting schizophrenia (SZ) from healthy control. However, the traditional FCN analysis with stationary assumption, i.e., static functional connectivity network (SFCN) at the time only measures the simple functional connectivity among brain regions, ignoring the dynamic changes of functional connectivity and the high-order dynamic interactions. In this article, the dynamic functional connectivity network (DFCN) is constructed to delineate the characteristic of connectivity variation across time. A high-order functional connectivity network (HFCN) designed based on DFCN, could characterize more complex spatial interactions across multiple brain regions with the potential to reflect complex functional segregation and integration. Specifically, the temporal variability and the high-order network topology features, which characterize the brain FCNs from region and connectivity aspects, are extracted from DFCN and HFCN, respectively. Experiment results on SZ identification prove that our method is more effective (i.e., obtaining a significantly higher classification accuracy, 81.82%) than other competing methods. Post hoc inspection of the informative features in the individualized classification task further could serve as the potential biomarkers for identifying associated aberrant connectivity in SZ.

Introduction

Schizophrenia (SZ) is a serious chronic mental disorder (Insel, 2010) that affects the brain, with the symptoms typically including hallucinations, delusions, emotional disturbances and confusion in language and behavior (Galderisi et al., 2018). Due to the complexity of the brain, the heterogeneity of disease, and the overlapping symptoms between different psychiatric disorders, clinical diagnosis is still difficult (Jablensky, 2022). Despite decades of research, the pathological mechanism of SZ is still not well understood. Research suggests that SZ may have differences in the connectivity of brain regions (Lynall et al., 2010; Lin et al., 2021).

Benefiting from the development of neuroimaging and computer technology, resting-state functional magnetic resonance images (rs-fMRI), measuring the low-frequency fluctuations in the blood-oxygen-level-dependent (BOLD) signals, have been proved to capture spontaneous neural activity of the brain that reflects functional organization (Barkhof et al., 2014; Zhang et al., 2019). Graph theory-based brain networks can effectively represent the coordinated interactions of neural activities between different ROIs/nodes and reflect the functional connectivity (FC; edges between nodes) in the complex brain (Vecchio et al., 2017; Sporns, 2022). A series of functional brain network construction methods based on BOLD signals such as correlation-based methods (Sporns, 2011) and sparse-learning-based methods (Yu et al., 2017; Gao et al., 2020) have been widely applied in brain disease diagnosis. For instance, Li et al. (2017) remodeled the functional brain network for autism spectrum disorder diagnosis based on Pearson’s correlation [PC; the most commonly used correlation-based method (Smith et al., 2011)] using sparsity and scale-free prior, respectively. Yu et al. (2019) constructed a weighted graph regularized sparse brain network for mild cognitive impairment diagnosis. However, these methods assumed the temporal stationarity of functional brain network (aka, static functional connectivity network, SFCN) across the duration of the scan, which ignore the rich dynamic information during scan. The dynamic interactions of ROIs may be critical for diagnosing brain disorders.

Researchers suggest that through dynamic reconfiguration of the brain, its various parts can adaptively coordinate and integrate in response to rapidly changing stimuli (Cocchi et al., 2013; Kupis et al., 2021). To capture dynamic interactions between ROIs, dynamic FCN (DFCN), constructed with the sliding window method (Hutchison et al., 2013), could consider the changing characteristics of the brain and simulate these changes by measuring the correlation of brain regions in a short period (Damaraju et al., 2014). For instance, Wee et al. (2016) used the sliding window method and PC to construct a DFCN for the identification of early mild cognitive impairment. Zhang J. et al. (2016) proposed a DFCN-based temporal variability measure to discover differences in the temporal variability of brain networks in subjects with different diseases through the diagnosis of multiple brain diseases. DFCN can simply and effectively represent the dynamic changes in the interaction patterns between ROIs, but cannot represent the complex interaction patterns between multiple brain regions in a deeper way. Neurological findings have demonstrated that a brain region predominantly interacts directly with multiple ROIs in neurological processes (Huang et al., 2010).

High-order functional connectivity network (HFCN) has great potential for mining complex spatial interactions across multiple brain regions. Many studies have proposed different HFCN construction methods. For instance, through extracting the FC series from the DFCN, Chen et al. (2016) proposed to cluster them and constructed the HFCN based on the pairwise correlation of averaged FC series. Zhao et al. (2018) proposed a multi-level HFCN construction method based on SFCN. A multi-level high-order network is obtained through the “correlation’s correlation” strategy, and then features are extracted from each level of the network for autism classification. Li Y. et al. (2020) proposed an ultra-least squares group constrained ultra-orthogonal least squares regression algorithm to construct low-order and high-order brain function networks, and successfully realized the prediction of mild cognitive impairment. High-order networks can provide additional spatial information for disease identification by characterizing complex brain region interactions. However, most traditional functional brain network analysis focus on DFCN or HFCN, which may ignore the complementarity that exists between them.

In this article, we propose a temporal-spatial dynamic functional connectivity method for the diagnosis of SZ. Specifically, the sliding window method is used to construct the DFCN of each subject, and the HFCN is constructed based on DFCN, and then two different features are extracted from these two networks respectively. The features extracted from DFCN reflect the relationship of brain functional network connectivity over time, which called temporal variability. The features extracted from HFCN represent the functional interaction between different modules in the brain network to a certain extent, reflecting the modularity of the brain network, which called spatial variability. Through fusing these two kinds of features that characterize the brain FCNs from region and connectivity aspects, we utilize Least Absolute Shrinkage and Selection Operator (Lasso) (Tibshirani, 1996) for feature selection and build a linear classifier for identification of SZ and HC.

Materials and methods

Materials

Data acquisition

The rs-fMRI data comes from the Centers of Biomedical Research Excellence (COBRE¹), which contains 57 SZ subjects and 64 HC subjects. The detailed information of the dataset is shown in Table 1. According to the fourth edition of Diagnostic and Statistical Manual of Mental Disorders (DSM-IV) criteria, the diagnosis of SZ is identified by psychiatrist and symptom is assessed using the Positive and Negative Syndrome Scale (PANSS). All subjects are screened and excluded if (1) history of DSM-IV disorders; (2) history of significant head trauma; and (3) history of substance abuse. The single-shot full k-space echo-planar imaging (EPI) with ramp sampling correction using the inter-commissural line (AC-PC) as a reference is used to obtain rest data, where repetition time (TR)/time-of-echo (TE) = 2,000/30 ms, in-plane voxel = 64 × 64, 32 slices, voxel size = 3 mm × 3 mm × 4 mm, field of view (FOV) = 256 mm × 256 mm and number of volumes = 150.

TABLE 1

Table 1. Demographic information of the participants from COBRE dataset.

Data preprocessing

The first 10 volumes are discarded to allow for scanner stabilization and the subjects’ adaptation to environment. The remain volumes are preprocessed using DPABI toolbox (Yan et al., 2016). The processing flow includes slice timing correction, realignment, spatial normalization [a standard Montreal Neurological Institute (MNI) template, resampled to 3 mm × 3 mm × 3 mm], and spatial smoothing using an 8 mm FWHM Gaussian kernel. The band-pass filter (0.01–0.08 Hz) is used to reduce low frequency drift and high frequency physiological noise. The automated anatomical labeling atlas (AAL2) (Tzourio-Mazoyer et al., 2002) is used to parcellate the brain into 120 ROIs. After data preprocessing, the mean time series are extracted from each ROI, and the time point signal of each ROI is normalized.

Construction of low-order functional connectivity network

Construction of static functional connectivity network

Pearson’s correlation measures the functional connectivity of ROIs by calculating the correlation coefficient between the average rs-fMRI time series of ROIs, defines as

$W_{ij}=\frac{{(x_i-{\overline{x}}_i)}^T(x_j-{\overline{x}}_j)}{\sqrt{{(x_i-{\overline{x}}_i)}^T(x_i-{\overline{x}}_i)}\sqrt{{(x_j-{\overline{x}}_j)}^T(x_j-{\overline{x}}_j)}}(1)$

where x_i = [x_1i,x_2i,⋯,x_Mi] ∈ ℝ^M and x_j = [x_1j,x_2j,⋯, x_Mj] ∈ ℝ^M denote the average rs-fMRI time series from the i-th and j-th ROIs respectively. M = 140 is the length of the time series. The ${\overline{x}}_i$ and ${\overline{x}}_j$ are the mean value of x_i and x_j respectively, and W_ij reflects the correlation between the i-th and j-th ROIs. The whole brain PC matrix is W = X^TX ∈ ℝ^N×N, where X = [x₁,x₂,⋯,x_N] ∈ ℝ^M×N denotes the whole-brain BOLD signals, and N = 120 denotes the number of ROIs. W characterizes the static correlation of a pair of ROIs throughout the scan time, which ignores dynamic interaction information in the rs-fMRI signal.

Construction of dynamic functional connectivity network

To capture the dynamic nature of neural activity, the entire BLOD time series is partitioned into multiple segments of overlapping subseries to construct the sub-networks using a sliding window method. As shown in Figure 1, the number of subseries K is

$K=\lceil \frac{M-L}{\mathrm{\Delta }}\rceil +1(2)$

FIGURE 1

Figure 1. Framework for construction of DFCN and HFCN. The DFCN is first constructed by computing the PC coefficient under k different sub-windows based on the BOLD signal. Then the DFCNs of all subjects are stacked to extract FC time series and clustered. Finally, the HFCN is obtained by calculating the PC of different FC sequence clusters. The figure in the upper right illustrates the dimension reduction of higher-order networks through clustering.

where ⌈⌉ denotes the ceiling function, and Δ denotes the sliding step. L represents the length of the window. It is worth noting that a smaller L can capture short-term inter-ROIs dynamic interactions, but is also more susceptible to noise. Correspondingly, a larger L is not conducive to the detection of dynamic interaction information between brain regions. L = {20, 30, 40, 50, 60, 70} will be automatically searched in subsequent experiments. The step size is uniformly set to 1, since a smaller Δ can obtain more subsequences to make the results more continuous.

Based on the above, suppose that x_i(k) = [x_1i(k), x_2i(k),⋯,x_Li(k)] ∈ ℝ^L(k = 1,2,⋯,K) denotes the average region time series of i-th brain region under the k-th sub-time window, the PC coefficient of i-th and j-th ROIs are W_ij(k). Therefore, the sub-network under the k-th sub-time window is defined as

$C_{DL}(k)=\left[\begin{array}{c}{W}_{11}(k)W_{12}(k)\mathrm{\cdots }{W}_{1N}(k)\hfill \\ W_{21}(k)W_{22}(k)\mathrm{\cdots }{W}_{2N}(k)\hfill \\ \mathrm{⋮}\mathit{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{⋮}\mathit{\hspace{1em}\hspace{0.5em}\hspace{1em}}\mathrm{\ddots }\mathit{\hspace{1em}\hspace{1em}}\mathrm{⋮}\hfill \\ W_{N1}(k)W_{N2}(k)\mathrm{\cdots }{W}_{NN}(k)\hfill \end{array}\right](3)$

where C_DL(k) represents the low-order DFCN of the subjects under the k-th sub-time window. In essence, DFCN reflect the temporal variability of functional connections between ROIs.

Construction of high-order functional connectivity network

To capture high-order functional interactions across ROIs, the “correlation’s correlation” principle is adopted to obtain HFCN (Chen et al., 2016; Zhang H. et al., 2016). Specifically, the FC series $h_{ij}^s=[W_{ij}^s(1),W_{ij}^s(2),\mathrm{\cdots },W_{ij}^s(K)]\in {ℝ}^K$ between i-th and j-th ROIs for the s-th subject can be obtained from C_DL(k)(see Eq. 3). $h_{ij}^s$ reflects the short time-dependent series of the correlation between i-th and j-th ROIs. Then the high-order connectivity for the s-th subject can be computed as $H_{ij,pq}^s=corr(h_{ij}^s,h_{pq}^s)$for each pair of FC series h${}_{ij}^s$ and $h_{pq}^s$, where corr denotes PC coefficient. Thus, H can extract interaction information from up to four ROIs, whereas the correlation in Eq. 3 involves only two ROIs. This suggests that high-order connectivity can characterize more complex interactions.

However, since the number of FC series is N × (N − 1) / 2, the dimension of the high-order connectivity network grows exponentially with N², which will cause a large amount of computation complexity and the poor generalization performance. Therefore, the mean clustering algorithm (Ward, 1963) is employed to group the FC series into different clusters. Specifically, the DFCN of all S subjects are stacked together in the following

$\{C_{DL}^1,C_{DL}^2,\mathrm{\cdots },C_{DL}^S\}\in {ℝ}^{N\times N\times KS}(4)$

where $C_{DL}^s=\{C_{DL}^s(1),C_{DL}^s(2),\mathrm{\cdots },C_{DL}^s(K)\}\in {ℝ}^{N\times N\times K}$ denotes the stacked sequence of all K sub-networks along the time for the s-th subject. Then FC of i-th and j-th ROIs is denoted by $h_{ij}=[h_{ij}^1,h_{ij}^2,\mathrm{\cdots },h_{ij}^S]\in {ℝ}^{KS}$ for all S subjects. After that, the N × N FC series in the stacked network are clustered into U clusters, and the corresponding average FC series of U clusters are calculated respectively as follows

${\overline{h}}_u=\frac{{\sum }_{(i,j)\in {\mathrm{\Omega }}_u}{h}_{ij}}{\left|{\mathrm{\Omega }}_u\right|}(5)$

where Ω_u denotes the set of the u-th cluster, and |Ω_u| denotes the number of elements in the set. Thus, the average FC series of the s-th subjects in the u-th cluster is obtained as

${\overline{h}}_u^s=\frac{{\sum }_{(i,j)\in {\mathrm{\Omega }}_u}{h}_{ij}^s}{\left|{\mathrm{\Omega }}_u\right|}(6)$

It is important to note that the size of U directly affects the difference of different clusters. In order to simulate the real high-order interaction of ROIs, we set the parameter U as one of {300, 400, 500, 600, 700, 800}. The cluster parameter U will be automatically searched in subsequent experiments, and the effect will be given in the discussion section. Based on the average FC series obtained in Eq. 6, the high-order functional connectivity between the u-th cluster and the v-th cluster can be expressed as $H_{uv}=\text{corr}({\overline{h}}_u,{\overline{h}}_v)$. The HFCN can be defined as

$C_{DH}=\left[\begin{array}{c}{H}_{11}{H}_{12}\mathrm{\cdots }{H}_{1U}\hfill \\ H_{21}{H}_{22}\mathrm{\cdots }{H}_{2U}\hfill \\ \mathrm{⋮}\mathit{\hspace{1em}}\mathrm{⋮}\mathit{\hspace{1em}}\mathrm{\ddots }\mathrm{⋮}\hfill \\ H_{U1}{H}_{U2}\mathrm{\cdots }{H}_{UU}\hfill \end{array}\right](7)$

where C_DH denotes the inter-modulation relationship between the functional connections of multiple ROIs, and reflects the spatial variability of the brain network. Furthermore, the HFCN obtained by “correlation’s correlation” strategy assumes that the functional brain network is fully connected. In order to better represent the real characteristics of the brain, we sparse the HFCN with percentage thresholding (edges with top λ_HD% connection strength are retained) (Qiao et al., 2018). The high-order network thresholding parameter λ_HD range is {0:10:90}, which will be automatically searched in subsequent experiments, and the effect will be given in the discussion section.

Feature extraction and selection

Temporal variability feature extraction

Temporal variability is a measure that describes the degree to which the FC of a particular ROI to all other ROIs changes over time (Zhang J. et al., 2016). The temporal variability of i-th ROI can be expressed as

$V_i=1-\overline{corr(w_i(k),w_i(l))},k,l=1,2,3,\mathrm{\dots },K,k\ne l(8)$

where w_i(k) = [W₁₁(k),W₂₁(k),⋯,W_N1(k)]^T ∈ ℝ^N×1 denotes FC under the k-th window between the i-th ROI and all other ROIs, k and l represent different time windows. The $\overline{corr(w_i(k),w_i(l))}$ denotes the average of PC coefficient in different time windows for i-th ROI. The $\overline{corr(w_i(k),w_i(l))}$ compares the functional architecture, i.e., overall functional connectivity profile associated with ROI k across different time windows. Then, a deduction from 1 indicates temporal variability of a ROI. In this way, it is possible to both target specific ROIs and assess the extent to which the functional architecture has changed over time at the global level.

High-order network topology feature extraction

We use the node degree and weighted-graph local clustering coefficient (Rubinov and Sporns, 2010) to extract features, separately. The node degree is the number of nodes directly connected to i-th node. The degree of the i-th node is defined as

$D_i=\sum _{j=1}^n{W}_{ij}(9)$

The weighted-graph local clustering coefficient quantifies the probability that neighbors of this vertex are also connected to each other, which can better characterize the clique structure of the FCN. The weighted-graph local clustering coefficient f_i for vertex i can be defined as follows

$t_i=\frac{1}{2}\sum _{j,h\in N}{\left(W_{ij}{W}_{ih}{W}_{jh}\right)}^{1/3}(10)$

$f_i=\frac{1}{n}\sum _{i\in N}\frac{2t_i}{r_i\left(r_i-1\right)}(11)$

where j and h are the neighbor nodes of node i, r_i is the number of neighbor nodes of node i, and t_i is the number of edges between all nodes connected to node i. The number of HFCN nodes depends on the number of clusters U, and the dimensions of the weighted-graph local clustering coefficient and node degree are both U.

Feature selection

The temporal variability extracted from DFCN is low-order features (N dimensional vectors), and the weighted-graph local clustering coefficient and node degree extracted from HFCN are high-order features (2U dimensional vectors). In order to remove irrelevant and redundant features for improving the generalization performance, Lasso regression algorithm (Tibshirani, 1996) is employed to find the most discriminative features relevant to SZ.

For considering the interaction among features, Lasso regression combines all features (with the dimension of m = N + 2U) to select feature subsets and the involved ℓ₁-norm sparsity regularization is used to control the dimension of feature subsets. Lasso can be defined as

$\underset{\mathrm{\alpha }}{min}\frac{1}{2}{||Y-F\mathrm{\alpha }||}_2^2+{\mathrm{\lambda }}_{FS}{||\mathrm{\alpha }||}_1(12)$

where Y = [y₁,y₂,….,y_S] ∈ ℝ^S×1 denotes the ground-truth label of the subjects (i.e., 1 for SZ, −1 for HC), and F ∈ ℝ^{S×(m + 1)} denotes a sparse dictionary that consists of the original feature matrix of S subjects and a column vector with all elements 1. α denotes the regression coefficient vector, and the position of its non-zero elements is the index of the selected feature in original feature set. λ_FS represents the penalty coefficient, which is used to control the sparsity of the feature subset. λ_FS is automatically searched in the experiments and its impact on model performance is given in the discussion section.

Experiments and results

Experiment settings

After selecting important features with Lasso, the support vector machine (SVM) with a linear kernel is trained for SZ identification. There are four hyperparameters in our proposed model, including window length L (range of {20:10:70}), clustering number U (range of {300:100:800}), high-order network thresholding λ_HD (range of {0:10:90}), and Lasso feature selection parameter λ_FS (range of {0.1:0.1:0.6}). The nested leave-one-out cross-validation (LOOCV) is used to evaluate classification performance and optimize those hyperparameters. Specifically, for S subjects in the dataset, S-1 subjects are used for training while the left-out one is used for testing. This procedure is repeated S times for evaluating the classification performance. To determine the optimal combination of the above four parameters, LOOCV is executed again on S-1 training subjects in the above process. Then, by applying the combination of optimal parameters on the S-1 different training subsets, we train S-1 classifiers to classify the test subject, and the final classification result is determined via majority voting. After repeating the above process S times, an overall cross-validation classification accuracy is calculated.

Evaluation methodology and results

Evaluation metrics

The five metrics are employed to evaluate the performance, including accuracy (ACC), sensitivity (SEN), specificity (SPE), area under curve (AUC), and receiver operating characteristic (ROC). Those evaluation metrics can be defined as ACC = (TP + TN)/(TP + FN + TN + FP), SEN = TP/(TP + FN), SPE = TN/(TN + FP), where TP, TN, FP, and FN are true positive, true negative, false positive, and false negative, respectively.

Method comparison

We compare three methods with features extracted from SFCN, DFCN, and HFCN, respectively. The feature selection and classification procedures are same with our proposed method. Specifically, the upper triangular in static FC network (i.e., the matrix W) is extracted and converted into a long vector. The vectors of all subjects are standardized and used as input to SVM after feature selection. Similarly, the temporal variability features and high-order features extracted from the low-order DFCN and the HFCN are separately used as the input of SVM after feature selection. The classification results are listed in Table 2 and the ROC curve is shown in Figure 2.

TABLE 2

Table 2. The comparison of performance in the classification of SZs and HCs by different methods.

FIGURE 2

Figure 2. Receiver operating characteristic curves achieved by four different methods in SZ vs. HC classification.

It can be seen from Table 2 that the performance of the static FC method is less than 70%. This may be because the SFCN only delineates the functional connectivity between paired ROIs throughout the scan, ignoring the dynamic functional interactions between multiple brain regions during the scan. The temporal variability method is improved by 3.3% compared to the static FC method. This illustrates the DFCN and temporal variability method of simulating the dynamic interaction of the brain by measuring the short-term correlation of brain regions, which is beneficial to the capture of discriminative features. Through the clustering of FC series and the calculation of high-order correlations, HFCN can further characterize the modulation relationship between FC of multiple ROIs, thereby characterizing the complex abstract interactions of the brain. Therefore, the high-order features (node degree and graph-weighted local clustering coefficients) extracted from HFCN can further improve the classification performance. Finally, the classification performance can be effectively improved by fusing low- and high-order features. This may be because temporal variability and high-order topological features describe properties between brain regions from different perspectives, indicating that the complementary information is critical for disease diagnosis.

Connectivity network analysis

One participant is randomly selected to investigate the constructed connectivity network, shown in Figure 3. It can be seen from Figures 3F,G that there is less difference between the brain network constructed by SFCN using the full time series and the averaged DFCN. This may be because both are essentially correlation measures for BOLD signals, while subtle differences may be the short-term FC properties obtained by DFCN. From Figures 3A–E, it can be seen that measuring the correlation of BOLD signals at a fine-grained level can effectively reflect the dynamic functional interaction of ROIs. It is worth noting that the FC intensity varied significantly throughout the scanning phase. This change is reflected in spatially and temporally, that is, the FC has undergone interactive reorganization in different ROIs, which changes the modular structure of the brain network. Figure 3H shows that network topology in the HFCN is more complex, implying it containing more high-order ROI’s interaction.

FIGURE 3

Figure 3. Visualization of SFCN, DFCN, and HFCN when the window length L = 40. (A) All the rearranged DFCN upper triangular elements along the time window, where the horizontal axis represents the 101 time windows, and the vertical axis represents the 7,140 functional connectivities. (B–E) The DFCN sub-networks under the k-th time windows, where k = 1, k = 40, k = 80, and k = 101, respectively. (F,G) The average of the sub-networks under whole time windows of DFCN and the static FC network based on pair-wise correlation of BOLD signals, respectively. (H) The constructed high-order FC network.

Furthermore, we investigate the connectivity differences in the HFCN of SZ and HC from both individual and group perspectives. From Figures 4A,B, it can be seen that the high-order functional connectivity of HC is overall higher than that of SZ, which is consistent with the existing researches on lower-order brain networks (Su et al., 2013; Yu et al., 2020). This illustrates that there is a wide range of functional disconnections in SZ patients (Lynall et al., 2010). This disconnection is manifested in complex functional interactions across multiple ROIs (Li et al., 2019), leading to abnormalities in the functional separation and integration characteristics of SZ.

FIGURE 4

Figure 4. Visualization of HFCN for SZ and HC (L = 40, U = 600, λ_HD = 0.5). (A) The HFCN and local details of the randomly selected SZ patient and HC. (B) The averaged HFCN and local details of SZ and HC groups.

Discussion

Model parameter sensitivity analysis

To further analyze the influence of different parameters on the classification results, the classification performance under different window L, different cluster U, as well as different parameter combinations are discussed.

The effect of L and U on the performance

We fix L at one of {20, 30, 40, 50, 60, 70}, and use nested LOOCV to optimize several other parameter combinations to analyze the effect of L on model performance. It can be seen from Figure 5A that ACC first increases and then decreases with L, and reaches the highest when L = 40. This indicates that the size of the L is important for the detection of dynamic interactions between ROIs. A smaller window length can capture shorter-term functional changes, but is more susceptible to noise. A larger window length would make the performance more stable, but the dynamic information cannot be effectively detected. This is consistent with previous research (Li Y. et al., 2020). In addition, when the L becomes large, the length of the FC series decreases, which may lead to a decrease in the reliability of HFCN. The choice of cluster number U also has a greater impact on the classification results. As shown in Figure 5B, the impact of U on performance is analyzed by fixing the number of clusters U to be one of {300, 400, 500, 600, 700, 800}. It can be seen that ACC is more sensitive to changes in the number of clusters, and the highest ACC can be obtained when U is 600.

FIGURE 5

Figure 5. The effect of different parameters on classification performance. (A) The different window length L; (B) the different cluster number U.

The clustering algorithm can reduce the dimension of HFCN and effectively decrease the computational complexity as well as the redundant features. The FC series of some ROI pairs are depicted in Figure 6A. For many ROI pairs, their temporal correlations are significantly different, such as the FC strength and change trend. From Figures 6B,C, it can be seen that the functional connectivity order of three clusters has a consistent change trend within the cluster, but has a significantly different change trend among the clusters. The clustering algorithm is able to find the underlying dominant dynamic patterns in all FC sequences.

FIGURE 6

Figure 6. Some selected FC series and their clustering results (L = 40, U = 600). (A) The original FC series. (B) Red, green, and blue indicate three different clusters. (C) The averaged FC series of the three clusters, separately.

The effect of different parameter combinations on the performance

To further analyze the impact of different parameter combinations on the performance, we conduct an additional LOOCV experiment, and the accuracy under different combinations of parameters (L, U, λ_FS) with HFCN threshold parameter λ_HD = 0.5 is shown in Figure 7. It can be seen that when L = 40 and U = 600, the accuracy is generally higher, which is consistent with the previous results when L and U are analyzed separately. In addition, the Lasso feature selection parameter λ_FS has a greater impact on the classification results. When λ_FS is set to a small value (0.1–0.3), high classification performance can be obtained. This is because λ_FS controls the number of features in the feature selection process. When λ_FS is too large, some discriminative features are discarded together with redundant features, which may be detrimental to the classification performance.

FIGURE 7

Figure 7. Statistical classification accuracy is estimated by our proposed method with different values of hyperparameters. For clear visualization, HFCN thresholding parameter λ_HD is set to 0.5. These results are obtained by LOOCV on whole data.

Most discriminative features

To find out the discriminative features or clusters for SZ identification, we fix the number of clusters U at 600 and analyze the importance of each feature in the classification. Specifically, during the feature selection process in each validation step, the selected features for classification might be different for different training datasets. We count the total weight of each feature in all training and classification, and use it to measure the classification contribution of features, as shown in Figure 8. It can be seen that only a few features have a large weight, indicating that they are the discriminative features for classification. Furthermore, most of the features with large weight belong to the time variability features and node degree features of high-order network, which reveals that these two kinds of features play a more important role in SZ diagnosis.

FIGURE 8

Figure 8. The normalized classification weights of all features (U = 600).

To further find the most important classification features, we conduct a separate analysis of time variability features and high-order features. For the time variability features, we take the absolute value of the normalized feature weights and sort them. Table 3 lists the brain regions corresponding to the top ten features with large weights. We can see that the brain regions with highest contribution including the frontal lobes, Cerebellum7bR, Vermis7, and limbic system, which are visualized in Figure 9. These regions have been suggested to be related to SZ by previous studies (Weinberger et al., 1994; Grace, 2000; Ichimiya et al., 2001).

TABLE 3

Table 3. Brain regions corresponding to the 10 time variability features with large weights.

FIGURE 9

Figure 9. The projections of ROIs corresponding to the ten time variability features with the large weights on the cortical surface.

For high-order features, we add the feature weights of the node degree and the weighted-graph local clustering coefficient of HFCN in classification to obtain the importance of a certain cluster. The five most important clusters are filtered out for classification and presented in the form of a chord diagram, as shown in Figure 10. It can be seen that the five clusters with the largest weights involve 58 functional connections and 21 brain regions (20 mirrored brain regions of both hemispheres among them). Five brain regions are consistent with Figure 9 and Table 3, including the insular inferior frontal gyrus, paracentral lobule, retroorbital gyrus, supplementary motor area, and cerebellum7b. Figure 11 shows the related FC within the five clusters. It can be seen that these functional connections mainly involve the frontal lobe, parietal lobe, and cerebellum. Among them, the frontal lobe-cerebellar and the parietal lobe-cerebellar connection are more numerous and weighted. Figure 12 shows the specific top twenty clusters with the highest classification weights. It can be seen that the cerebellum and frontal lobes still have constant importance in these top clusters. In addition, the visual network regions of the occipital lobe, such as the fusiform gyrus and the cortex around the talar fissure are also of greater importance for the diagnosis of SZ.

FIGURE 10

Figure 10. Functional connectivity contained in the five clusters with the largest classification weights. Each node in the figure represents a ROI, and the line between the two nodes represents functional connectivity. The functional connections contained in different clusters are shown in different colors, and the importance increases sequentially from Cluster 5 to Cluster 1.

FIGURE 11

Figure 11. The most discriminative connections and brain regions in the whole brain view. The size of a node, calculated from the total weight of all clusters that containing the node, indicates the importance of the node during the classification. The thickness of the edge represents the classification weight of the cluster that connects the edge.

FIGURE 12

Figure 12. The top 20 clusters with the highest classification weight. The functional subnetworks represented by different colors are the same as in Figure 11.

In this article, cerebellum (especially Cerebelum7bR) and frontal lobe-related features play the most important roles in SZ classification, which confirmed by features extracted from both DFCN and HFCN. In fact, abnormalities in the cerebellum and frontal lobe of patients with SZ have been identified in other studies (van den Heuvel et al., 2010; Wang et al., 2014). For example, Andreasen et al. (1996) found that SZ patients showed dysfunction of the prefrontal-thalamic-cerebellar circuit, and proposed the Cognitive Dysmetria hypothesis (Andreasen et al., 1998), which believed that abnormalities in the neural circuits connecting the thalamus and cerebellum may lead to abnormalities in cognitive control, coordination and affective aspects. This may be the underlying psychopathological cause of the complex symptoms of SZ. Functionally, the frontal lobe is involved in attention regulation, abstract rules, social behavior, etc., and is also linked to the symptoms of SZ. In terms of the cerebellum, it is not only mainly responsible for the function of motor control but also has been reported in the research of SZ that there are abnormalities of neural circuits between it and other cortex (Cao and Cannon, 2019; Kelly et al., 2021). In addition, the striatum and multiple regions of the parietal-premotor cortex (paracentral lobule, postcentral gyrus, precuneus, supramarginal gyrus, inferior parietal gyrus, and superior parietal gyrus) also show important roles in SZ classification (Salvador et al., 2005). The striatum has been a key region in the study of SZ, and has an important relationship with the dopamine hypothesis of the etiology of SZ (Li A. et al., 2020). Many antipsychotic drugs rely on blockade of dopamine receptors in the striatum (Kegeles et al., 2010). In the AAL2 template, the striatum is divided into three brain regions: caudate nucleus, putamen, and pallidum, two of which appear in the five clusters with the largest classification weight and show abnormal connectivity to the orbitofrontal cortex. The paracentral lobules, postcentral gyrus, and other brain regions of the parietal-premotor cortex are mainly involved in sensory-motor, spatial attention, etc. Abnormalities in these brain areas may explain the sensory, thinking and behavioral disorders in SZ.

Conclusion

In this article, we propose a temporal-spatial dynamic functional connectivity analysis method that combines DFCN and HFCN for SZ diagnosis. Specifically, we first construct a DFCN and extract temporal variability features that reflect temporal dynamic information. Based on DFCN, the HFCN is constructed by means of clustering algorithm and “correlation’s correlation” strategy, and the node degree and graph-weighted local clustering coefficient features of the sparse HFCN are extracted. Then, the most discriminatory features are selected with Lasso method for SZ classification. The proposed method is verified on a real SZ dataset and demonstrates promising performance via comparison with the competitive brain network analysis methods. The abnormal brain regions detected in this paper can provide a direction for more detailed research on the pathology of SZ and the search for biological markers in the future.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository can be found in the article.

Ethics statement

The studies involving human participants were reviewed and approved by IORG0000256 – University of New Mexico. Written informed consent was obtained from all participants for their participation in this study.

Author contributions

CP, HY, and RY have developed the modeling algorithm and architecture. XF and XZ have preprocessed the data. All authors have contributed to preparation of the article, figures, and tables, approved it for publication, and listed have made a substantial contribution to the work.

Funding

This work was partly supported by the National Natural Science Foundation of China (61906171 and 62006072), Key Technologies Research and Development Program of Henan Province (222102210108), and Key Project at Central Government Level (2060302).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

1. ^ http://fcon_1000.projects.nitrc.org/indi/retro/cobre.html

References

Andreasen, N.C., O’Leary, D.S., Cizadlo, T., Arndt, S., Rezai, K., Ponto, L., et al. (1996). Schizophrenia and cognitive dysmetria: a positron-emission tomography study of dysfunctional prefrontal-thalamic-cerebellar circuitry. Proc. Natl. Acad. Sci. U.S.A 93, 9985–9990. doi: 10.1073/pnas.93.18.9985

PubMed Abstract | CrossRef Full Text | Google Scholar

Andreasen, N.C., Paradiso, S and O’leary, D. S. (1998). “Cognitive dysmetria” as an integrative theory of schizophrenia: a dysfunction in cortical-subcortical-cerebellar circuitry? Schizophr. Bull. 24, 203–218. doi: 10.1093/oxfordjournals.schbul.a033321

PubMed Abstract | CrossRef Full Text | Google Scholar

Barkhof, F., Haller, S., and Rombouts, S. A. (2014). Resting-state functional MR imaging: a new window to the brain. Radiology 272, 29–49. doi: 10.1148/radiol.14132388

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, H and Cannon, T. D. (2019). Cerebellar dysfunction and schizophrenia: from “cognitive dysmetria” to a potential therapeutic target. Am. J. Psychiatry. 176, 498–500 doi: 10.1176/appi.ajp.2019.19050480

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Zhang, H., Gao, Y., Wee, C.Y., Li, G., Shen, D. et al. (2016). Initiative, High-order resting-state functional connectivity network for MCI classification. Hum. Brain Mapp. 37, 3282–3296. doi: 10.1002/hbm.23240

PubMed Abstract | CrossRef Full Text | Google Scholar

Cocchi, L., Zalesky, A., Fornito, A., and Mattingley, J. B. (2013). Dynamic cooperation and competition between brain systems during cognitive control. Trends Cogn. Sci. 17, 493–501. doi: 10.1016/j.tics.2013.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Damaraju, E., Allen, E. A., Belger, A., Ford, J. M., McEwen, S., Mathalon, D., et al. (2014). Dynamic functional connectivity analysis reveals transient states of dysconnectivity in schizophrenia. NeuroImage 5, 298–308. doi: 10.1016/j.nicl.2014.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Galderisi, S., Mucci, A., Buchanan, R. W., and Arango, C. (2018). Negative symptoms of schizophrenia: new developments and unanswered research questions. Lancet Psychiatry 5, 664–677. doi: 10.1016/S2215-0366(18)30050-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, X., Xu, X., Hua, X., Wang, P., Li, W., and Li, R. (2020). Group similarity constraint functional brain network estimation for Mild Cognitive Impairment classification. Front. Neurosci. 14:165. doi: 10.3389/fnins.2020.00165

PubMed Abstract | CrossRef Full Text | Google Scholar

Grace, A. A. (2000). Gating of information flow within the limbic system and the pathophysiology of schizophrenia. Brain Res. Rev. 31, 330–341. doi: 10.1016/S0165-0173(99)00049-1

CrossRef Full Text | Google Scholar

Huang, S., Li, J., Sun, L., Ye, J., Fleisher, A., Wu, T. et al. (2010). Initiative, Learning brain connectivity of Alzheimer’s disease by sparse inverse covariance estimation. NeuroImage 50, 935–949. doi: 10.1016/j.neuroimage.2009.12.120

PubMed Abstract | CrossRef Full Text | Google Scholar

Hutchison, R. M., Womelsdorf, T., Allen, E. A., Bandettini, P. A., Calhoun, V. D., Corbetta, M., et al. (2013). Dynamic functional connectivity: promise, issues, and interpretations. Neuroimage 80, 360–378. doi: 10.1016/j.neuroimage.2013.05.079

PubMed Abstract | CrossRef Full Text | Google Scholar

Ichimiya, T., Okubo, Y., Suhara, T and Sudo, Y. (2001). Reduced volume of the cerebellar vermis in neuroleptic-naive schizophrenia. Biol. Psychiatry 49, 20–27. doi: 10.1016/s0006-3223(00)01081-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Insel, T. R. (2010). Rethinking schizophrenia. Nature 468, 187–193. doi: 10.1038/nature09552

PubMed Abstract | CrossRef Full Text | Google Scholar

Jablensky, A. (2022). The diagnostic concept of schizophrenia: its history, evolution, and future prospects. Dialogues Clin. Neurosci. 12, 271–287. doi: 10.31887/DCNS.2010.12.3/ajablensky

PubMed Abstract | CrossRef Full Text | Google Scholar

Kegeles, L.S., Abi-Dargham, A., Frankle, W. G., Gil, R., Cooper, T.B., Slifstein, M., et al. (2010). Increased synaptic dopamine function in associative regions of the striatum in schizophrenia. Arch. Gen. Psychiatry 67, 231–239. doi: 10.1001/archgenpsychiatry.2010.10

PubMed Abstract | CrossRef Full Text | Google Scholar

Kelly, E., Escamilla, C. O and Tsai, P. T. (2021). Cerebellar dysfunction in autism spectrum disorders: Deriving mechanistic insights from an internal model framework. Neuroscience 462, 274–287. doi: 10.1016/j.neuroscience.2020.11.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Kupis, L., Goodman, Z. T., Kornfeld, S., Hoang, S., Romero, C., Dirks, B., et al. (2021). Brain dynamics underlying cognitive flexibility across the lifespan. Cereb. Cortex 31, 5263–5274. doi: 10.1093/cercor/bhab156

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, A., Zalesky, A., Yue, W., Howes, O., Yan, H., Liu, Y, et al. (2020). A neuroimaging biomarker for striatal dysfunction in schizophrenia. Nat. Med. 26, 558–565. doi: 10.1038/s41591-020-0793-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, S., Hu, N., Zhang, W., Tao, B., Dai, J., Gong, Y et al. (2019). Dysconnectivity of multiple brain networks in schizophrenia: a meta-analysis of resting-state functional connectivity. Front. Psychiatry 10:482. doi: 10.3389/fpsyt.2019.00482

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, W., Wang, Z., Zhang, L., Qiao, L., and Shen, D. (2017). Remodeling Pearson’s correlation for functional brain network estimation and autism spectrum disorder identification. Front. Neuroinform. 11:55. doi: 10.3389/fninf.2017.00055

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Liu, J., Peng, Z., Sheng, C., Kim, M., Yap, P.-T., Wee, C.-Y., and Shen, D. (2020). Fusion of ULS group constrained high-and low-order sparse functional connectivity networks for MCI classification. Neuroinformatics 18, 124. doi: 10.1007/s12021-019-09418-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, X., Li, W., Dong, G., Wang, Q., Sun, H., Shi, J., et al. (2021). Characteristics of multimodal brain connectomics in patients with schizophrenia and the unaffected first-degree relatives. Front. Cell Dev. Biol. 9:631864. doi: 10.3389/fcell.2021.631864

PubMed Abstract | CrossRef Full Text | Google Scholar

Lynall, M.-E., Bassett, D. S., Kerwin, R., McKenna, P. J., Kitzbichler, M., Muller, U., et al. (2010). Functional connectivity and brain networks in schizophrenia. J. Neurosci. 30, 9477–9487.

Google Scholar

Qiao, L., Zhang, L., Chen, S and Shen, D. (2018). Data-driven graph construction and graph learning: A review. Neurocomputing 312, 336–351. doi: 10.1016/j.neucom.2018.05.084

CrossRef Full Text | Google Scholar

Rubinov, M and Sporns, O. (2010). Complex network measures of brain connectivity: uses and interpretations. Neuroimage 52, 1059–1069. doi: 10.1016/j.neuroimage.2009.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Salvador, R., Suckling, J., Coleman, M.R., Pickard, J. D., D. Menon, and E. Bullmore (2005). Neurophysiological architecture of functional magnetic resonance images of human brain. Cereb. Cortex 15, 1332–1342. doi: 10.1093/cercor/bhi016

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, S. M., Miller, K. L., Salimi-Khorshidi, G., Webster, M., Beckmann, C. F., Nichols, T. E., et al. (2011). Network modelling methods for FMRI. Neuroimage 54, 875–891.

Google Scholar

Sporns, O. (2011). The human connectome: a complex network. Ann. N.Y. Acad. Sci. 1224, 109–125. doi: 10.1111/j.1749-6632.2010.05888.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Sporns, O. (2022). Graph theory methods: applications in brain networks. Dialogues Clin. Neurosci. 20, 111–121. doi: 10.31887/DCNS.2018.20.2/osporns

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, L., Wang, L., Shen, H., Feng, G., and Hu, D. (2013). Discriminative analysis of non-linear brain connectivity in schizophrenia: an fMRI Study. Front. Hum. Neurosci. 7:702. doi: 10.3389/fnhum.2013.00702

PubMed Abstract | CrossRef Full Text | Google Scholar

Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. R. Stat. Soc. Series B 58, 267–288.

Google Scholar

Tzourio-Mazoyer, N., Landeau, B., Papathanassiou, D., Crivello, F., and Etard, O. (2002). Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. Neuroimage 15, 273–289.

Google Scholar

van den Heuvel, M. P., Mandl, R.C., Stam, C.J., Kahn, R.S, and Pol, H.E.H. (2010). Aberrant frontal and temporal complex network structure in schizophrenia: a graph theoretical analysis. J. Neurosci. 30, 15915–15926. doi: 10.1523/JNEUROSCI.2874-10.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Vecchio, F., Miraglia, F., and Rossini, P. M. (2017). Connectome: Graph theory application in functional brain network architecture. Clin. Neurophysiol. Prac. 2, 206–213.

Google Scholar

Wang, L., Zou, F., Shao, Y., Ye, E., Jin, X., Tan, S et al. (2014). Disruptive changes of cerebellar functional connectivity with the default mode network in schizophrenia. Schizophr. Res. 160, 67–72.

Google Scholar

Ward, Jr, J. H. (1963). Hierarchical grouping to optimize an objective function. J. Am. Stat. Assoc. 58, 236–244. doi: 10.1080/01621459.1963.10500845

CrossRef Full Text | Google Scholar

Wee, C.-Y., Yang, S., Yap, P.-T., and Shen, D. (2016). Sparse temporally dynamic resting-state functional connectivity networks for early MCI identification. Brain Imaging Behav. 10, 342–356. doi: 10.1007/s11682-015-9408-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Weinberger, D.R., Aloia, M.S., Goldberg, T.E, and Berman, K.F. (1994). The frontal lobes and schizophrenia. J. Neuropsychiatry Clin. Neurosci. 6, 419–27 doi: 10.1176/jnp.6.4.419

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, C.-G., Wang, X.-D., Zuo, X.-N, and Zang, Y.-F. (2016). DPABI: data processing & analysis for (resting-state) brain imaging. Neuroinformatics 14, 339–351. doi: 10.1007/s12021-016-9299-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, R., Qiao, L., Chen, M., Lee, S.-W., Fei, X., and Shen, D. (2019). Weighted graph regularized sparse brain network construction for MCI identification. Pattern Recognition 90, 220–231. doi: 10.1016/j.patcog.2019.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, R., Yu, H., and Wan, H. (2020). Research on brain network for schizophrenia classification based on resting-state functional magnetic resonance imaging. Shengwu Yixue Gongchengxue Zazhi 37, 661–669. doi: 10.7507/1001-5515.201908007

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, R., Zhang, H., An, L., Chen, X., Wei, Z., and Shen, D. (2017). Connectivity strength-weighted sparse group representation-based brain network construction for M CI classification. Hum. Brain Mapp. 38, 2370–2383. doi: 10.1002/hbm.23524

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Chen, X., Shi, F., Li, G., Kim, M., Giannakopoulos, P et al. (2016). Topographical information-based high-order functional connectivity and its application in abnormality detection for mild cognitive impairment. J. Alzheimer Dis. 54, 1095–1112. doi: 10.3233/JAD-160092

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Cheng, W., Liu, Z., Zhang, K., Lei, X., Yao, Y, et al. (2016). Neural, electrophysiological and anatomical basis of brain-network variability and its characteristic changes in mental disorders. Brain 139, 2307–2321. doi: 10.1093/brain/aww143

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, T., Zhao, Z., Zhang, C., Zhang, J., Jin, Z., and Li, L. (2019). Classification of early and late mild cognitive impairment using functional brain network of resting-state fMRI. Front. Psychiatry 10:572. doi: 10.3389/fpsyt.2019.00572

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, F., Zhang, H. I., Rekik, Z. A, and Shen, D. (2018). Diagnosis of autism spectrum disorders using multi-level high-order functional networks derived from resting-state functional mri. Front. Hum. Neurosci. 12:184. doi: 10.3389/fnhum.2018.00184

PubMed Abstract | CrossRef Full Text | Google Scholar