Front. Neurosci.Frontiers in NeuroscienceFront. Neurosci.1662-453XFrontiers Media S.A.10.3389/fnins.2022.965937NeuroscienceOriginal ResearchTemporal-spatial dynamic functional connectivity analysis in schizophrenia classificationPanCong1†YuHaifei2†FeiXuan3ZhengXingjuan4YuRenping1*1Henan Key Laboratory of Brain Science and Brain-Computer Interface Technology, School of Electrical Engineering, Zhengzhou University, Zhengzhou, China2Aviation Maintenance NCO Academy, Air Force Engineering University, Xinyang, China3School of Artificial Intelligence and Big Data, Henan University of Technology, Zhengzhou, China4Gaoyou Hospital Affiliated to Soochow University, Gaoyou People’s Hospital, Gaoyou, China
Edited by: Zhengxia Wang, Hainan University, China
Reviewed by: Weikai Li, Chongqing Jiaotong University, China; Sijie Niu, University of Jinan, China
*Correspondence: Renping Yu, yurenping@zzu.edu.cn
†These authors have contributed equally to this work
This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
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.
functional magnetic resonance imagingschizophrenia classificationdynamic functional connectivity networkshigh-order functional connectivity networktemporal variability6190617162006072National Natural Science Foundation of China10.13039/501100001809Introduction
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 methodsMaterialsData acquisition
The rs-fMRI data comes from the Centers of Biomedical Research Excellence (COBRE1), 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.
Demographic information of the participants from COBRE dataset.
SZ
HC
P-value
Numbers
57
64
Gender (male/female)
48/9
45/19
0.072
Ages
36.684 ± 13.620
35.313 ± 11.804
0.554
Handedness (L/R/both)
8/48/1
1/62/1
0.022
The ages are denoted as mean ± SD. L/R/both: left/right/(both left and right). P-value < 0.05 is marked in bold.
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 networkConstruction 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
where xi = [x1i,x2i,⋯,xMi] ∈ ℝM and xj = [x1j,x2j,⋯, xMj] ∈ ℝ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 x¯i and x¯j are the mean value of xi and xj respectively, and Wij reflects the correlation between the i-th and j-th ROIs. The whole brain PC matrix is W = XTX ∈ ℝN×N, where X = [x1,x2,⋯,xN] ∈ ℝ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=⌈M-LΔ⌉+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 xi(k) = [x1i(k), x2i(k),⋯,xLi(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 Wij(k). Therefore, the sub-network under the k-th sub-time window is defined as
where CDL(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 hijs=[Wijs(1),Wijs(2),⋯,Wijs(K)]∈ℝK between i-th and j-th ROIs for the s-th subject can be obtained from CDL(k)(see Eq. 3). hijs 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 Hij,pqs=corr(hijs,hpqs)for each pair of FC series hijs and hpqs, 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 N2, 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
{CDL1,CDL2,⋯,CDLS}∈ℝN×N×KS
where CDLs={CDLs(1),CDLs(2),⋯,CDLs(K)}∈ℝN×N×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 hij=[hij1,hij2,⋯,hijS]∈ℝ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
h¯u=∑(i,j)∈Ωuhij|Ωu|
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
h¯us=∑(i,j)∈Ωuhijs|Ωu|
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 Huv=corr(h¯u,h¯v). The HFCN can be defined as
where CDH 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 selectionTemporal 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
Vi=1-corr(wi(k),wi(l))¯,k,l=1,2,3,…,K,k≠l
where wi(k) = [W11(k),W21(k),⋯,WN1(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 corr(wi(k),wi(l))¯ denotes the average of PC coefficient in different time windows for i-th ROI. The corr(wi(k),wi(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
Di=∑j=1nWij
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 fi for vertex i can be defined as follows
where j and h are the neighbor nodes of node i, ri is the number of neighbor nodes of node i, and ti 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 ℓ1-norm sparsity regularization is used to control the dimension of feature subsets. Lasso can be defined as
minα12||Y-Fα||22+λFS||α||1
where Y = [y1,y2,….,yS] ∈ ℝ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 resultsExperiment 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 resultsEvaluation 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.
The comparison of performance in the classification of SZs and HCs by different methods.
Method
ACC (%)
SEN (%)
SPE (%)
AUC (%)
Static FC
68.60
77.19
60.94
74.42
Time variability
71.90
68.42
75.00
75.27
High-order features
73.55
75.44
71.88
78.84
Time variability + high-order features
81.82
82.46
81.25
89.26
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.
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.
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.
DiscussionModel 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 <italic>L</italic> and <italic>U</italic> 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.
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.
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.
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.
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).
Brain regions corresponding to the 10 time variability features with large weights.
ROI index
ROI name
Anatomical zonation
106
Cerebelum7bR
Cerebellum
23
Gyrus rectus
Frontal lobe
7
Inferior frontal gyrus, opercular part
Frontal lobe
26
Medial orbital gyrus
Frontal lobe
8
Inferior frontal gyrus, opercular part
Frontal lobe
73
Paracentral lobule
Frontal lobe
28
Anterior orbital gyrus
Frontal lobe
46
Amygdala
Limbic system
16
Supplementary motor area
Frontal lobe
117
Vermis7
Cerebellum
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.
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.
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.
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.
ReferencesAndreasenN.C.O’LearyD.S.CizadloT.ArndtS.RezaiK.PontoL. (1996). Schizophrenia and cognitive dysmetria: a positron-emission tomography study of dysfunctional prefrontal-thalamic-cerebellar circuitry.Proc. Natl. Acad. Sci. U.S.A939985–9990. 10.1073/pnas.93.18.99858790444AndreasenN.C.ParadisoSO’learyD. S. (1998). “Cognitive dysmetria” as an integrative theory of schizophrenia: a dysfunction in cortical-subcortical-cerebellar circuitry?Schizophr. Bull.24203–218. 10.1093/oxfordjournals.schbul.a0333219613621BarkhofF.HallerS.RomboutsS. A. (2014). Resting-state functional MR imaging: a new window to the brain.Radiology27229–49. 10.1148/radiol.1413238824956047CaoHCannonT. D. (2019). Cerebellar dysfunction and schizophrenia: from “cognitive dysmetria” to a potential therapeutic target.Am. J. Psychiatry.176498–50010.1176/appi.ajp.2019.1905048031256620ChenX.ZhangH.GaoY.WeeC.Y.LiG.ShenD. (2016). Initiative, High-order resting-state functional connectivity network for MCI classification.Hum. Brain Mapp.373282–3296. 10.1002/hbm.2324027144538CocchiL.ZaleskyA.FornitoA.MattingleyJ. B. (2013). Dynamic cooperation and competition between brain systems during cognitive control.Trends Cogn. Sci.17493–501. 10.1016/j.tics.2013.08.00624021711DamarajuE.AllenE. A.BelgerA.FordJ. M.McEwenS.MathalonD. (2014). Dynamic functional connectivity analysis reveals transient states of dysconnectivity in schizophrenia.NeuroImage5298–308. 10.1016/j.nicl.2014.07.00325161896GalderisiS.MucciA.BuchananR. W.ArangoC. (2018). Negative symptoms of schizophrenia: new developments and unanswered research questions.Lancet Psychiatry5664–677. 10.1016/S2215-0366(18)30050-629602739GaoX.XuX.HuaX.WangP.LiW.LiR. (2020). Group similarity constraint functional brain network estimation for Mild Cognitive Impairment classification.Front. Neurosci.14:165. 10.3389/fnins.2020.0016532210747GraceA. A. (2000). Gating of information flow within the limbic system and the pathophysiology of schizophrenia.Brain Res. Rev.31330–341. 10.1016/S0165-0173(99)00049-1HuangS.LiJ.SunL.YeJ.FleisherA.WuT. (2010). Initiative, Learning brain connectivity of Alzheimer’s disease by sparse inverse covariance estimation.NeuroImage50935–949. 10.1016/j.neuroimage.2009.12.12020079441HutchisonR. M.WomelsdorfT.AllenE. A.BandettiniP. A.CalhounV. D.CorbettaM. (2013). Dynamic functional connectivity: promise, issues, and interpretations.Neuroimage80360–378. 10.1016/j.neuroimage.2013.05.07923707587IchimiyaT.OkuboY.SuharaTSudoY. (2001). Reduced volume of the cerebellar vermis in neuroleptic-naive schizophrenia.Biol. Psychiatry4920–27. 10.1016/s0006-3223(00)01081-711163776InselT. R. (2010). Rethinking schizophrenia.Nature468187–193. 10.1038/nature0955221068826JablenskyA. (2022). The diagnostic concept of schizophrenia: its history, evolution, and future prospects.Dialogues Clin. Neurosci.12271–287. 10.31887/DCNS.2010.12.3/ajablensky20954425KegelesL.S.Abi-DarghamA.FrankleW. G.GilR.CooperT.B.SlifsteinM. (2010). Increased synaptic dopamine function in associative regions of the striatum in schizophrenia.Arch. Gen. Psychiatry67231–239. 10.1001/archgenpsychiatry.2010.1020194823KellyE.EscamillaC. OTsaiP. T. (2021). Cerebellar dysfunction in autism spectrum disorders: Deriving mechanistic insights from an internal model framework.Neuroscience462274–287. 10.1016/j.neuroscience.2020.11.01233253824KupisL.GoodmanZ. T.KornfeldS.HoangS.RomeroC.DirksB. (2021). Brain dynamics underlying cognitive flexibility across the lifespan.Cereb. Cortex315263–5274. 10.1093/cercor/bhab15634145442LiA.ZaleskyA.YueW.HowesO.YanH.LiuY (2020). A neuroimaging biomarker for striatal dysfunction in schizophrenia.Nat. Med.26558–565. 10.1038/s41591-020-0793-832251404LiS.HuN.ZhangW.TaoB.DaiJ.GongY (2019). Dysconnectivity of multiple brain networks in schizophrenia: a meta-analysis of resting-state functional connectivity.Front. Psychiatry10:482. 10.3389/fpsyt.2019.0048231354545LiW.WangZ.ZhangL.QiaoL.ShenD. (2017). Remodeling Pearson’s correlation for functional brain network estimation and autism spectrum disorder identification.Front. Neuroinform.11:55. 10.3389/fninf.2017.0005528912708LiY.LiuJ.PengZ.ShengC.KimM.YapP.-T.WeeC.-Y.ShenD. (2020). Fusion of ULS group constrained high-and low-order sparse functional connectivity networks for MCI classification.Neuroinformatics18124. 10.1007/s12021-019-09418-x30982183LinX.LiW.DongG.WangQ.SunH.ShiJ. (2021). Characteristics of multimodal brain connectomics in patients with schizophrenia and the unaffected first-degree relatives.Front. Cell Dev. Biol.9:631864. 10.3389/fcell.2021.63186433718367LynallM.-E.BassettD. S.KerwinR.McKennaP. J.KitzbichlerM.MullerU. (2010). Functional connectivity and brain networks in schizophrenia.J. Neurosci.309477–9487.QiaoL.ZhangL.ChenSShenD. (2018). Data-driven graph construction and graph learning: A review.Neurocomputing312336–351. 10.1016/j.neucom.2018.05.084RubinovMSpornsO. (2010). Complex network measures of brain connectivity: uses and interpretations.Neuroimage521059–1069. 10.1016/j.neuroimage.2009.10.00319819337SalvadorR.SucklingJ.ColemanM.R.PickardJ. D.MenonD.BullmoreE. (2005). Neurophysiological architecture of functional magnetic resonance images of human brain.Cereb. Cortex151332–1342. 10.1093/cercor/bhi01615635061SmithS. M.MillerK. L.Salimi-KhorshidiG.WebsterM.BeckmannC. F.NicholsT. E. (2011). Network modelling methods for FMRI.Neuroimage54875–891.SpornsO. (2011). The human connectome: a complex network.Ann. N.Y. Acad. Sci.1224109–125. 10.1111/j.1749-6632.2010.05888.x21251014SpornsO. (2022). Graph theory methods: applications in brain networks.Dialogues Clin. Neurosci.20111–121. 10.31887/DCNS.2018.20.2/osporns30250388SuL.WangL.ShenH.FengG.HuD. (2013). Discriminative analysis of non-linear brain connectivity in schizophrenia: an fMRI Study.Front. Hum. Neurosci.7:702. 10.3389/fnhum.2013.0070224155713TibshiraniR. (1996). Regression shrinkage and selection via the lasso.J. R. Stat. Soc. Series B58267–288.Tzourio-MazoyerN.LandeauB.PapathanassiouD.CrivelloF.EtardO. (2002). Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain.Neuroimage15273–289.van den HeuvelM. P.MandlR.C.StamC.J.KahnR.SPolH.E.H. (2010). Aberrant frontal and temporal complex network structure in schizophrenia: a graph theoretical analysis.J. Neurosci.3015915–15926. 10.1523/JNEUROSCI.2874-10.201021106830VecchioF.MiragliaF.RossiniP. M. (2017). Connectome: Graph theory application in functional brain network architecture.Clin. Neurophysiol. Prac.2206–213.WangL.ZouF.ShaoY.YeE.JinX.TanS (2014). Disruptive changes of cerebellar functional connectivity with the default mode network in schizophrenia.Schizophr. Res.16067–72.WardJ. H.Jr (1963). Hierarchical grouping to optimize an objective function.J. Am. Stat. Assoc.58236–244. 10.1080/01621459.1963.10500845WeeC.-Y.YangS.YapP.-T.ShenD. (2016). Sparse temporally dynamic resting-state functional connectivity networks for early MCI identification.Brain Imaging Behav.10342–356. 10.1007/s11682-015-9408-226123390WeinbergerD.R.AloiaM.S.GoldbergT.EBermanK.F. (1994). The frontal lobes and schizophrenia.J. Neuropsychiatry Clin. Neurosci.6419–2710.1176/jnp.6.4.4197841813YanC.-G.WangX.-D.ZuoX.-NZangY.-F. (2016). DPABI: data processing & analysis for (resting-state) brain imaging.Neuroinformatics14339–351. 10.1007/s12021-016-9299-427075850YuR.QiaoL.ChenM.LeeS.-W.FeiX.ShenD. (2019). Weighted graph regularized sparse brain network construction for MCI identification.Pattern Recognition90220–231. 10.1016/j.patcog.2019.01.01531579345YuR.YuH.WanH. (2020). Research on brain network for schizophrenia classification based on resting-state functional magnetic resonance imaging.Shengwu Yixue Gongchengxue Zazhi37661–669. 10.7507/1001-5515.20190800732840083YuR.ZhangH.AnL.ChenX.WeiZ.ShenD. (2017). Connectivity strength-weighted sparse group representation-based brain network construction for M CI classification.Hum. Brain Mapp.382370–2383. 10.1002/hbm.2352428150897ZhangH.ChenX.ShiF.LiG.KimM.GiannakopoulosP (2016). Topographical information-based high-order functional connectivity and its application in abnormality detection for mild cognitive impairment.J. Alzheimer Dis.541095–1112. 10.3233/JAD-16009227567817ZhangJ.ChengW.LiuZ.ZhangK.LeiX.YaoY (2016). Neural, electrophysiological and anatomical basis of brain-network variability and its characteristic changes in mental disorders.Brain1392307–2321. 10.1093/brain/aww14327421791ZhangT.ZhaoZ.ZhangC.ZhangJ.JinZ.LiL. (2019). Classification of early and late mild cognitive impairment using functional brain network of resting-state fMRI.Front. Psychiatry10:572. 10.3389/fpsyt.2019.0057231555157ZhaoF.ZhangH.I.RekikZ. AShenD. (2018). Diagnosis of autism spectrum disorders using multi-level high-order functional networks derived from resting-state functional mri.Front. Hum. Neurosci.12:184. 10.3389/fnhum.2018.0018429867410