High-order brain functional network for electroencephalography-based diagnosis of major depressive disorder

Brain functional network (BFN) based on electroencephalography (EEG) has been widely used to diagnose brain diseases, such as major depressive disorder (MDD). However, most existing BFNs only consider the correlation between two channels, ignoring the high-level interaction among multiple channels that contain more rich information for diagnosing brain diseases. In such a sense, the BFN is called low-order BFN (LO-BFN). In order to fully explore the high-level interactive information among multiple channels of the EEG signals, a scheme for constructing a high-order BFN (HO-BFN) based on the “correlation’s correlation” strategy is proposed in this paper. Specifically, the entire EEG time series is firstly divided into multiple epochs by sliding window. For each epoch, the short-term correlation between channels is calculated to construct a LO-BFN. The correlation time series of all channel pairs are formulated by these LO-BFNs obtained from all epochs to describe the dynamic change of short-term correlation along the time. To construct HO-BFN, we cluster all correlation time series to avoid the problems caused by high dimensionality, and the correlation of the average correlation time series from different clusters is calculated to reflect the high-order correlation among multiple channels. Experimental results demonstrate the efficiency of the proposed HO-BFN in MDD identification, and its integration with the LO-BFN can further improve the recognition rate.

Introduction

Major depressive disorder (MDD) is a kind of common brain disease that is characterized by persistent and significant low mood, slow thinking, and cognitive function impairment (Holma et al., 2014; LeMoult and Gotlib, 2019; Liu et al., 2020). In the statistics of the World Health Organization (WHO), MDD has become the second largest serious disease in the world (Melchior et al., 2013) and has brought a heavy burden on patients and their families (Zhang et al., 2019). According to medical researches, the accurate early identification of MDD is important, and it can not only effectively relieve the pain of the patients, but also directly reduce the tragedy of suicide (Cipriani et al., 2018; Grunebaum et al., 2018). However, early neuroimaging-based MDD diagnosis is very challenging, because the changes of brain functional connectivity (FC) are considerably complicated. Electroencephalography (EEG) of high temporal resolution (Babajani and Soltanian-Zadeh, 2006) can well describe the temporal evolution of complex FC during brain activity and thereby becomes the best choice for MDD research.

Brain functional network (BFN) constructed based on EEG has been widely used in the diagnosis of brain diseases (Wang et al., 2015; Li et al., 2020). Since brain activity is dynamic in nature, some studies have shown that the dynamic change of FC over the whole scanning time may be the intrinsic feature of brain function (Damaraju et al., 2014; Cohen and D’Esposito, 2016; Kudela et al., 2017). Many studies try to describe the dynamic changes of FC between channels by using sliding windows to construct BFN and the relationship between these dynamic changes and brain diseases (Wee et al., 2016; Guo et al., 2017; Sun et al., 2019; Zhang et al., 2021). Sun et al. (2019) constructed BFN based on EEG signals by sliding windows, confirming that MDD had abnormal cognitive processing. Zhang et al. (2021) used sliding windows to construct BFN based on EEG signals, and the results showed that the brain regions of MDD patients were significantly altered.

Although the aforementioned EEG-based BFN helps us to understand the brain activities of the MDD patients, most of them only reflect the low-order FC (LO-FC) between two channels (as shown in Figure 1A), ignoring the fact that high-order FC (HO-FC) among channels could also be changed for MDD patients (as shown in Figure 1B). For the ease of description, we call BFN based on conventional LO-FC as low-order BFN (LO-BFN). In essence, brain activity is complex, and the HO-FC usually contains more abstract information than the LO-FC, and it helps to reveal high-level and more complex interaction information (Chen et al., 2016). Therefore, it is of clinical significance to investigate effective methods of constructing a high-order BFN (HO-BFN) that better reflects the complex interaction among multiple channels and simulates the mechanisms of the deep brain, providing rich discriminative information for the diagnosis of mental disorders (Plis et al., 2014; Chen et al., 2016; Zhang et al., 2016).

FIGURE 1

Figure 1. The intuitive explanation of brain network research based on EEG signals. (A) The low-order FC (LO-FC) between two channels. (B) High-order FC (HO-FC) among channels. (1) EEG signal acquisition(1) EEG signal acquisition, (2) calculation of correlation between two channels, and (3) calculation of correlation among multiple channels.

In this work, we propose a novel method to construct HO-BFN for MDD classification. Specifically, the entire EEG signals of a subject are divided into multiple overlapping time series by sliding windows, and the correlation between two channels within one window is computed as the LO-FC. The network constructed by the LO-FC is called the LO-BFN, reflecting the dynamic change of FC throughout the whole scanning time. Significantly different from the LO-BFN, each vertex of the HO-BFN represents one pair of channels, and each edge represents the correlation between the channel pairs. In this way, the HO-BFN involves more channels and can reveal high-level and more interaction among brain regions.

It is noted that the scale of the HO-BFN is very huge, and it may suffer from the dilemmas of high dimensionality and small size samples. The main reason lies in the fact that the scale of HO-BFN will increase exponentially when the number of EEG channels increases. To address the issue, we introduce hierarchical clustering (Yu et al., 2015) to construct HO-BFN. In other words, similar LO-FC time series are clustered into one group, and the average LO-FC time series are computed for each group. After that, the HO-BFN is constructed based on the correlation between the two groups. As a result, the HO-BFN constructed by hierarchical clustering can not only reduce computing time and memory requirements, but also reflect the HO-FC among multiple channel pairs (more than four channels) and capture more useful and complex information.

In summary, the main contributions of the paper line are twofold: (1) A HO-FC representation strategy is proposed to capture high-order interactions among multiple channels of EEG signals. In fact, the HO-BFN is used to characterize the complex interactions among brain regions, and it has been applied in fMRI and achieved good results (Zhao et al., 2018). However, to the best of our knowledge, few studies have used EEG-based HO-BFN to reveal the complex interactions among EEG channels. (2) The HO-BFN is constructed in both time and frequency domains based on the “correlation’s correlation” strategy. Specifically, we first compute the correlation between two channels to obtain the LO-BFN, and then the HO-BFN is subsequently derived by computing the correlation between each pair of channels from the LO-BFN. Then, we further apply HO-BFN to computer-aided diagnosis for MDD. The experimental results show that HO-BFN provides complementary identification information to the LO-BFN and that combining HO-BFN and LO-BFN can further improve the accuracy of MDD diagnosis.

Materials and data preprocessing

The EEG data used in this study came from the publicly available Multi-modal Open Dataset for Mental-disorder Analysis (MODMA) dataset (Cai et al., 2020). It included 24 patients with MDD (12 male and 12 female) and 29 normal controls (NC) (20 male and 9 female). The MDD group was 16–52 years old, and the NC group was 19–51 years old; all the subjects were right-handed, and their education level was primary school or above. In addition, MDD patients had a health Questionnaire 9-item (PHQ-9) (Spitzer et al., 1999) score greater than or equal to 5 and had not received psychotropic medication for 2 weeks. Table 1 shows the statistics of the subjects.

TABLE 1

Table 1. Demographic information of the subjects.

In data acquisition, 128-channel HydroCel Geodesic Sensor Net and Net Station acquisition software were used to record EEG signals for five minutes. Taking Cz as the reference, the sampling rate was 250 Hz. In order to reduce the interference of EEG data, the subjects were required to close their eyes and keep awake to avoid any unnecessary eye movement, saccade, and blink. The collected EEG data were filtered by 0.1–40 Hz and inhibited by 48–52 Hz to eliminate the data interference caused by baseline drift and electrical interference. The processed data was then re-referenced against REST (Yao, 2001). Finally, after the above steps, some high-power content was contained in the remaining data points and some EEG epochs were removed by the Artifact Subspace Reconstruction (ASR) plugin (Chang et al., 2018; Pion-Tonachini et al., 2018). In this study, theta (4–8 Hz), alpha (8–13 Hz), and beta (13–40 Hz) bands calculated by fast Fourier transform were selected in the frequency domain, which had been proved to be much distinct in the identification of depression (Nyström et al., 1986; Knott et al., 2001; Jaworska et al., 2012).

Methods

Figure 2 shows the overall pipeline of our method, which includes six steps: (1) constructing LO-BFN by sliding window; (2) clustering all low-order correlation time series; (3) constructing HO-BFN by calculating the correlation between clusters; (4) selecting and extracting features from each constructed BFN; (5) constructing support vector machines (SVMs) based on the selected features in both LO-BFN and HO-BFN; and (6) fusing the decision scores of multiple SVMs to predict whether each subject is MDD or NC.

FIGURE 2

Figure 2. The flowchart of the proposed BFN classification framework, including six main steps: (1) constructing LO-BFN; (2) clustering the time series; (3) constructing HO-BFN; (4) feature selection; (5) constructing SVM model; and (6) classification fusion.

Construction of low-order BFN

In order to construct LO-BFN, we first use the sliding window to divide the entire EEG signals into H = [(M−W)/s] + 1 overlapping windows, where M is the image volume during the entire scan period, and W and s are the window width and step size of the sliding window, respectively. Then, we calculate the correlation between $x_i^l\left(h\right)$ and $x_j^l\left(h\right)$, where $x_i^l\left(h\right)$ and $x_j^l\left(h\right)$ denote the i-th and j-th channels under the h-th window of the l-th subject, respectively.

In the frequency domain, we use the phase lag index (PLI) (Stam et al., 2007; Peraza et al., 2012) to calculate the channel correlation, which is robust to volume conduction artifacts. The PLI is denoted as:

$C_{ij}^l=\left|\frac{1}{N}\sum _{n=1}^{N}sign\left({\mathrm{\phi }}_i^l\left(t_n\right)-{\mathrm{\phi }}_j^l\left(t_n\right)\right)\right|(1)$

where N is the sample number, sign is the sign function, and ${\mathrm{\phi }}_i^l\left(t_n\right)-{\mathrm{\phi }}_j^l\left(t_n\right)$ is the phase synchronization of channels $x_i^l$ and $x_j^l$ at time t_n. Among them, φ^l(t_n) can be obtained by analyzing the signal based on Hilbert transform (Bruns, 2004).

On the other hand, in the time domain, we calculate the channel correlation by using the Pearson correlation coefficient (PCC) (Eslami and Saeed, 2018) as follows:

$C_{ij}^l=corr(x_i^l,x_j^l)(2)$

Therefore, for the l-th subject, the h-th subnetwork of LO-BFN is constructed as $G_L^l\left(h\right)=(\left\{x_i^l\left(h\right)\right\},\left\{C_{ij}^l\left(h\right)\right\})$ (as shown in Figure 3A), where $\left\{x_i^l(h)\right\}$ is vertices and $\left\{C_{ij}^l(h)\right\}$ is the weights of the edges connecting the i-th and j-th nodes. Then, we construct H subnetworks to form LO-BFN $G_L^l=[G_L^l\left(1\right),G_L^l\left(2\right),\mathrm{\dots },G_L^l\left(H\right)]$ for each subject, which describes the change in FC strength of all channel pairs over time.

FIGURE 3

Figure 3. Schematic diagram of different BFNs. (A) Schematic diagram of the brain region of LO-BFN. (B) Schematic diagram of brain regions of HO-BFN before clustering. (C) Schematic diagram of brain regions of HO-BFN after clustering.

Construction of high-order BFN

In order to capture high-level FCs, we adopt the strategy of “correlation’s correlation”. That is, based on the LO-BFN, the PCC is used to calculate the correlation between the LO-FC time series of the l-th subject, which is called HO-FC, denoted as:

$H_{ij,pq}^l=corr(C_{ij}^l,C_{pq}^l)(3)$

where $C_{ij}^l$ is the LO-FC time series between the i-th and the j-th channels of the l-th subject and $C_{pq}^l$ is the LO-FC time series between the p-th and the q-th channels of the l-th subject. Therefore, $H_{ij,pq}^l$ can represent the HO-FC among the four channels of the l-th subject at most, that is, the correlation between the FC between the i-th and j-th channels and the FC between the p-th and q-th channels. Physiologically, difference in FC among different channels in MDD patients and healthy individual subjects can be used to identify MDD.

If we have 128 channels in our study, the dimensionality of the constructed LO-BFN is 128 × 128. Thus, a large-scale HO-BFN will be constructed (as shown in Figure 3B) in Eq. 3; that is, the dimensionality is (128 × 128)², and the constructed network contains at least thousands of nodes and millions of edges. It is a critical problem that the dimensionality is too large, and it will introduce high computational complexity for the subsequent feature extraction and selection procedures. Besides, the generalization performance of the HO-BFN learning system may also degrade.

To this end, we will reduce the network dimensionality by clustering the LO-FC time series. Specifically, the LO-FC time series of the subjects are clustered into different clusters to find the potential interaction patterns. Then, the HO-FC between the respective average LO-FC time series in clusters is calculated. Compared with the previous large-scale HO-BFN, the dimensionality of the HO-BFN constructed in this way is greatly reduced (as shown in Figure 3C). This method not only preserves important interactive information, but also avoids the problems of high computational complexity and low generalization performance.

In order to ensure that the clustering results in different subjects are consistent, the FC matrices of all subjects are first accumulated together, so that to connect the time series of LO-FC in the same channel pairs in all subjects into a long vector. That is, the long vector connected by the time series ${\{C_{ij}^l\}}_{1\le l\le R}$ of the LO-FC between the i-th and j-th channels of all subjects is $C_{ij}=[C_{ij}^1,C_{ij}^2,\mathrm{\cdots },C_{ij}^R]$, where R is the number of subjects.

After the C_ij is obtained, we divide it into k clusters by hierarchical clustering. The cluster centers are calculated by averaging all the FC long vectors in the k-th cluster of the l-th subject as the following Eq. 4:

${\tilde{C}}_k^l=\frac{{\sum }_{C_{ij}\in {\mathrm{\omega }}_k}{C}_{ij}^l}{\left|{\mathrm{\omega }}_k\right|}(4)$

where |ω_k| denotes the total number of FC long vectors in cluster k. Finally, the PCC between the two clusters ${\tilde{C}}_{k1}^l$ and ${\tilde{C}}_{k2}^l$ of the l-th subject is calculated as follows:

${\tilde{H}}_{k1,k2}^l=corr({\tilde{C}}_{k1}^l,{\tilde{C}}_{k2}^l)(5)$

Finally, we obtain a small-scale HO-BFN $G_{\tilde{H}}^l=(\left\{{\tilde{C}}_{k1}^l\right\},\left\{{\tilde{H}}_{k1,k2}^l\right\})$, taking $\left\{{\tilde{C}}_{k1}^l\right\}$ as vertices and $\left\{{\tilde{H}}_{k1,k2}^l\right\}$ as the weights of edges.

Feature extraction, selection, classification, and fusion

Both LO-BFN and HO-BFN of the l-th subject, i.e., $G_L^l$ and $G_{\tilde{H}}^l$, are used for the subsequent classification. Due to the possible phase mismatch of all FC matrices of LO-BFN in each subject, the dynamic characteristics of different subjects do not completely correspond. Therefore, to avoid this situation, we calculate the average FC matrix of each subject’s LO-BFN as the low-order feature. For the ease of computation, we vectorize the averaged LO-FC matrix of the l-th subject into $f_L^l$, which is called the low-order feature vector. Similarly, we vectorize the FC matrix ${\tilde{H}}^l$ of HO-BFN of the l-th subject to $f_{\tilde{H}}^l$ as the high-order feature vector.

Both low-order feature vectors and high-order feature vectors have a large number of features, introducing irrelevant or redundant information for subsequent MDD classification. Therefore, we use t-test and Least Absolute Shrinkage and Selection Operator (LASSO) (Tibshirani, 2011) methods to select features for high-order and low-order feature vectors, which can effectively remove redundant features. Specifically, we first perform t-test on both the low-order feature vector f_L and the high-order feature vector $f_{\tilde{H}}$ on the training set and select the features that are significantly different from f_L and $f_{\tilde{H}}$ as preliminary features, denoted ${\overline{f}}_L$ and ${\overline{f}}_{\tilde{H}}$, respectively.

Then, we use LASSO to further remove redundant features and select the features most related to MDD. Let I^l be the labels for the l-th subject. Specifically, if the l-th subject is MDD, I^l = −1; if the l-th subject is NC, I^l = + 1; and α is set to be the weight vector for feature selection. The objective of LASSO is defined as:

$\underset{\mathrm{\alpha }}{min}\frac{1}{2}\sum _{l=1}^N{||I^l⟨-{\overline{f}}^l,\mathrm{\alpha }⟩||}_2^2+\mathrm{\lambda }{||\mathrm{\alpha }||}_1(6)$

where ⟨⋅,⋅⟩ is the inner product operator and λ is the regularization parameter. For simplicity, let ${\overline{\overline{f}}}_L$ and ${\overline{\overline{f}}}_{\tilde{H}}$ denote selected the final feature sets from the feature vectors ${\overline{f}}_L$ and ${\overline{f}}_{\tilde{H}}$, respectively.

Finally, the MDD is classified by using SVM with a linear kernel in this paper. Herein, we train two linear SVM classifiers by using both ${\overline{\overline{f}}}_L$ and ${\overline{\overline{f}}}_{\tilde{H}}$ features, respectively. The final results are obtained by fusing the decision scores of two SVM classifiers by linear combination. Among them, we set the weight β ∈ [0.1,0.2,…,0.9] for each SVM. The weights of the F classifiers to be fused are set to be β₁,β₂,…,β_F and satisfy ${\mathrm{\beta }}_1+{\mathrm{\beta }}_2+\mathrm{\dots }+\mathrm{\beta }{}_F=1$.

Experiments results

To evaluate the effectiveness of our proposed method, we analyze the impact of clustering parameters on HO-BFN. The classification ability of LO-BFN, HO-BFN, and their fusion for MDD is evaluated by six different indicators, i.e., accuracy (ACC), sensitivity or true positive rate (TPR), specificity or true negative rate (TNR), precision or positive predictive value (PPV), negative predictive value (NPV), and F1 score. Finally, we conduct feature analysis to understand the role of each channel in MDD diagnosis.

Furthermore, we adopt a nested ten-fold cross-validation (CV) strategy consisting of two nested loops to evaluate the effectiveness of our proposed method. In the outer loop, all data are randomly divided into ten subsets of roughly the same size, where one subset is selected as the testing set, while the other nine subsets are selected as the training set. In the inner loop, the data of the training set are merged and redivided into ten subsets of similar size, nine of which are used to adjust the hyper-parameters and one for model evaluation. We report the average accuracy of classification results across the ten-fold CV. Then, in order to avoid any possible bias in fold selection, this procedure is repeated 10 times, with a different random partitioning of samples each time. Finally, the average accuracy of 10 repetitions is reported. Since the performance of our method also depends on hyper-parameters, such as W and s in sliding window, k in clustering, p in t-test, λ in LASSO, and c in SVM model. The optimal hyper-parameters can be determined when the average classification accuracy reaches its optimum. In our experiment, we fix the size of the sliding window, i.e., W = 10000, s = 1000, and determine the optimal values of other parameters within the following range: k ∈ [100, 200, …, 800], p ∈ [0.01, 0.02, …, 0.05], λ ∈ [0.1, 0.2, …, 0.9], and c ∈ [2⁻⁴, 2⁻³, …,2⁴].

The clustering effectiveness on high-order BFN

To reduce the computational complexity of the HO-BFN, the hierarchical clustering is employed to construct HO-BFNs. The hyper-parameter k in hierarchical clustering indicates the cluster number, and it has a crucial influence on the constructed HO-BFN, and further affects the final classification results. In the experiment, we optimize HO-BFN by adjusting the clustering number k. Figure 4 shows the classification results when k takes different numbers.

FIGURE 4

Figure 4. Recognition accuracy of HO-BFN with different number of clusters.

As can be seen from Figure 4 that for the clustering parameter k, when the time domain takes the value of 500, alpha takes the value of 200, beta takes the value of 100, and theta takes the value of 600, the HO-BFN generates a relatively satisfactory classification result. The ACC is greatly affected by k; i.e., classification performance is very sensitive to clustering parameters. Different HO-BFNs have different performances, indicating that different HO-BFNs contain different levels of MDD diagnosis information. Therefore, we can conclude that it is necessary to select k carefully toward a better understanding of dynamics in brains.

Comparison of major depressive disorder diagnosis using different brain functional networks

In this subsection, we further study how HO-BFN contributes to MDD diagnosis. We construct and analyze HO-BFN from the time domain and alpha, beta, and theta bands in the frequency domain, respectively. In the experiment, we train and test the classifiers of LO-BFN and HO-BFN, respectively, and determine the parameter combination that can produce the best ACC.

Due to the complex connections in the brain, it is difficult to fully capture the relationship between different brain regions through a single type of BFN. In order to further improve the classification performance, we adopt the linear fusion of the SVM integrated decision score to combine LO-BFN and HO-BFN (Zhao et al., 2020) and analyze their fusion performance. In addition, we also believe that different BFNs constructed by the three bands in the frequency domain can reflect FC between channels from different views, which are complementary. Therefore, we further linearly fuse the decision scores of BFNs in three bands. Table 2 shows the classification performance of different BFNs, and the best classification results are highlighted in bold.

TABLE 2

Table 2. MDD classification using different BFNs.

From Table 2, we can draw the following conclusions: (1) BFNs constructed in different ways have different performance, implying that each BFN provides meaningful and various information for MDD identification; (2) fusing LO-BFN and HO-BFN is better than that of a single network, and the ACC is relatively increased in 4%, indicating that LO-BFN and HO-BFN are complementary to each other in classifying MDD; (3) the performance of each brain network in the frequency domain is obviously better than that of the BFNs in the time domain, indicating that it is more effective to extract features in the frequency domain.

The most discriminative features for major depressive disorder diagnosis

To identify the most discriminative features in MDD diagnosis, we select ten frequently selected LO-FC features and two frequently selected HO-FC features based on the t-test and LASSO regression ten-fold cross-validations of ten times. The higher selection frequency of FC indicates stronger reliability and discriminative ability.

Similar to previous studies (Bian et al., 2014), we divide the brain into five regions, i.e., frontal (F), left temporal (LT), central (C), right temporal (RT), and posterior (P). The frontal region is further divided into left frontal (LF) and right frontal (RF), the central region is further divided into left central (LC) and right central (RC), and the posterior is further divided into left posterior (LP) and right posterior (RP). Table 3 shows the details of these brain regions.

TABLE 3

Table 3. Brain regions corresponding to channels of interest.

Figure 5 shows the 10 most discriminative FC feature maps for the different LO-BFNs. The node color indicates the brain region the channel belongs to, the connection line represents the correlation between two channels, and the line width indicates the frequency. The thicker line represents the higher frequency. It can be seen that the frequently selected FC features in LO-BFN often appear in the LF, RC, RT, and LP regions in the alpha band; left central-right central (LC-RC), RT and left posterior-right posterior (LP-RP) regions in the beta band; the left frontal-right frontal (LF-RF), RC, RT, and LP regions in the theta band, and the LF, RC, left temporal-right temporal (LT-RT) and left posterior-right posterior (LP-RP) regions in the time domain.

FIGURE 5

Figure 5. The 10 most frequently selected connection diagrams for different LO-BFNs.

The HO-BFN is constructed based on clustering. Each cluster is used as the vertex of the network, and the HO-FC is used as the edge of the network. For HO-BFN, we select the two most discriminative cluster pairs and describe the FC feature diagram between these two cluster pairs in Figure 6, where each brain map represents a cluster, the connection between brain maps represents the connection between clusters, that is, HO-FC, and the connections in the brain map represent the FC between channels belonging to this cluster, that is, LO-FC. Table 4 shows the brain regions involved in the most discriminative clusters that are selected from different HO-BFNs. From Figure 6 and Table 4, we can observe that the most frequently selected features in HO-BFN are mainly distributed in LF-RF, LT-RT, and LC-RC regions of the theta band, other brain regions except RT in the time domain, as well as the whole brain region of the alpha and beta band.

FIGURE 6

Figure 6. The two most frequently selected connection diagrams for different HO-BFNs. Each brain diagram represents the FC of interest channels across the brain region in a cluster.

TABLE 4

Table 4. The brain regions selected from different HO-BFNs.

Discussion

In this paper, we propose a novel method to construct HO-BFN based on the “correlation’s correlation” strategy, which can simultaneously capture low-order features reflecting FC information between any two channels and high-order features reflecting FC information among multiple channels, so as to better simulate the mechanisms of the deep brain and provide more discriminative information for the diagnosis of mental disorders. We believe that different BFNs can mine disease-disturbed BFN variation information from different aspects, which has better performance in MDD classification experiments. According to the experimental results, we will discuss HO-BFNs in more detail.

In the experiment to explore the influence of the clustering parameter on the classification accuracy of HO-BFN, we find that the classification result changes with the change of the clustering parameter. When the clustering parameter is too small or too large, the accuracy of diagnosis will gradually decrease. This can be understood from two aspects: (1) when the clustering parameter is too small, the LO-FC time series with different dynamic changes may be divided into the same cluster, which will reduce the similarity of each cluster, so that the HO-BFN constructed by the mean sequence of each class as the vertex is unreliable; and (2) when the clustering parameter is too large, the LO-FC time series with similar dynamic changes may be divided into different clusters, which will increase the number of features extracted from HO-BFN, thus resulting in more redundant features and causing the decrease in accuracy and generalization ability. Thus, choosing the suitable number of clusters for HO-BFN is the key to achieving MDD classification and improving the classification accuracy.

According to the experimental results in Table 2, we find that the proposed BFNs in the frequency domain are usually more discriminative than in the time domain, and the alpha and theta bands in the frequency domain are more discriminative than the beta band, indicating that the brain functional structure of MDD patients has undergone significant changes in these two bands. Therefore, we believe that the BFNs of alpha and theta play an important role in the pathogenesis of MDD. Several previous studies have also reported the same or similar conclusion: abnormal brain function in MDD patients occurs at certain frequency bands. For example, Fingelkurts et al. (2010) showed that the FC in the alpha and theta bands of EEG in MDD was impaired. Hosseinifard et al. (2013) found that there were differences in alpha bands between the MDD group and the NC group.

To further prove the effectiveness of BFNs in MDD disease diagnosis, we trace the BFNs to which the features of the classifier used for training belong. The experimental results in Figures 5, 6 show that: (1) the channel pairs selected by LO-BFN hardly overlap with those of HO-BFN, indicating that the FC features extracted from LO-BFN and HO-BFN are complementary to each other. (2) RC and RT are significantly different in all LO-BFNs, and they are related to the regulation of attention, long-term memory, and emotion. In previous studies, Fan et al. (2013) found that abnormal RT superior gyrus activity could be a potential marker of suicidal tendencies in MDD patients. Sun et al. (2019) observed differences in MDD patients in the F, T, and C of the theta band and in the T and C of the alpha band. Zhang et al. (2021) found significant modifications in brain synchrony of LF, T, and RT in MDD patients. Our observations are generally consistent with these studies (Fan et al., 2013; Sun et al., 2019; Zhang et al., 2021).

However, this study has some limitations. Firstly, either the low-order or the high-order FC is based on the correlation instead of the inherent causality. The relatively small sample size and unbalanced data may also affect the result of the analysis. Although most conclusions obtained by our method are generally consistent with the previous relevant studies, the experimental results may have a potential bias due to the heterogeneity of the experiment (Hassler and Thadewald, 2010). In the future work, we will focus on the influence of nonsensical and biased correlation and investigate more advanced models to calculate correlation, such as causality inference, based on larger datasets to obtain more accurate and sufficient information about the brain changes of MDD patients. Secondly, the sliding window algorithm should set the step size and window width, but we fixed the window width and step length in this study. Further studies will investigate the influence of different parameter settings. Finally, we linearly combine the decision scores of the LO-BFN and HO-BFN at the decision-making level, and this linear combination may not fully mine the complementary information, thereby affecting the classification accuracy. Therefore, in our future work, we will further increase the accuracy of MDD diagnosis by using more advanced information fusion strategies.

Conclusion

In this paper, we propose a framework for constructing HO-BFN based on the “correlation’s correlation” strategy and capture the high-order correlations across different channels for MDD diagnosis. We use hierarchical clustering to reduce the computational complexity of the HO-BFN. Experimental results demonstrate that: (1) the proposed HO-BFN can provide discriminative information for the MDD identification. (2) Fusing high-order and low-order BFNs can significantly improve the recognition rate of MDD patients. (3) The most discriminative brain regions are associated with the regulation of attention, and these findings are consistent with the daily behavior of MDD patients.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: http://modma.lzu.edu.cn/data/index/.

Author contributions

FZ: conceptualization, methodology, writing – review, and editing. HP: conceptualization, software, writing – original draft, methodology, formal analysis, investigation, and validation. NL: validation. XC, HZ, NM, and YR: writing – review, and editing. All authors contributed to the article and approved the submitted version.

Funding

FZ was supported in part by the National Natural Science Foundation of China (62176140, 82001775, 61772319, 61873177, 61972235, 61976125, and 61976124) and Doctoral Scientific Research Foundation of Shandong Technology and Business University (No. BS202016).

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.

References

Babajani, A., and Soltanian-Zadeh, H. (2006). Integrated MEG/EEG and fMRI model based on neural masses. IEEE Trans. Biomed. Eng. 53, 1794–1801. doi: 10.1109/TBME.2006.873748

PubMed Abstract | CrossRef Full Text | Google Scholar

Bian, Z., Li, Q., Wang, L., Lu, C., Yin, S., and Li, X. (2014). Relative power and coherence of EEG series are related to amnestic mild cognitive impairment in diabetes. Front. Aging Neurosci. 6:11. doi: 10.3389/fnagi.2014.00011

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruns, A. (2004). Fourier-, Hilbert- and wavelet-based signal analysis: Are they really different approaches? J. Neurosci. Methods 137, 321–332. doi: 10.1016/j.jneumeth.2004.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, H., Gao, Y., Sun, S., Li, N., Tian, F., Xiao, H., et al. (2020). MODMA Dataset: A Multi-Modal Open Dataset for Mental-Disorder Analysis. Available online at: http://modma.lzu.edu.cn/data/index/ (accessed February 18, 2021).

Google Scholar

Chang, C. Y., Hsu, S. H., Pion-Tonachini, L., and Jung, T. P. (2018). Evaluation of artifact subspace reconstruction for automatic EEG artifact removal. Annu. Int. Conf. IEEE Eng. Med. Biol. Soc. 2018, 1242–1245.

Google Scholar

Chen, X., Zhang, H., Gao, Y., Wee, C. Y., Li, G., and Shen, D. (2016). 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

Cipriani, A., Furukawa, T. A., Salanti, G., Chaimani, A., Atkinson, L. Z., Ogawa, Y., et al. (2018). Comparative efficacy and acceptability of 21 antidepressant drugs for the acute treatment of adults with major depressive disorder: A systematic review and network meta-analysis. Lancet (London, England) 391, 1357–1366.

Google Scholar

Cohen, J. R., and D’Esposito, M. (2016). The Segregation and integration of distinct brain networks and their relationship to cognition. J. Neurosci. 36, 12083–12094.

Google Scholar

Damaraju, E., Allen, E. A., Belger, A., Ford, J. M., McEwen, S., Mathalon, D. H., et al. (2014). Dynamic functional connectivity analysis reveals transient states of dysconnectivity in schizophrenia. NeuroImage Clin. 5, 298–308.

Google Scholar

Eslami, T., and Saeed, F. (2018). Fast-GPU-PCC: A GPU-based technique to compute pairwise Pearson’s correlation coefficients for time series data-fMRI study. High Throughput 7:11. doi: 10.3390/ht7020011

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, T., Wu, X., Yao, L., and Dong, J. (2013). Abnormal baseline brain activity in suicidal and non-suicidal patients with major depressive disorder. Neurosci. Lett. 534, 35–40. doi: 10.1016/j.neulet.2012.11.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Fingelkurts, A. A., Fingelkurts, A. A., Rytsl, H., Suominen, K., Isomets, E., and Khknen, S. (2010). Impaired functional connectivity at EEG alpha and theta frequency bands in major depression. Hum. Brain Mapp. 28, 247–261.

Google Scholar

Grunebaum, M. F., Galfalvy, H. C., Choo, T. H., Keilp, J. G., Moitra, V. K., Parris, M. S., et al. (2018). Ketamine for rapid reduction of suicidal thoughts in major depression: A midazolam-controlled randomized clinical trial. Am. J. Psychiatry 175, 327–335. doi: 10.1176/appi.ajp.2017.17060647

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, H., Liu, L., Chen, J., Xu, Y., and Jie, X. (2017). Alzheimer classification using a minimum spanning tree of high-order functional network on fMRI dataset. Front. Neurosci. 11:639. doi: 10.3389/fnins.2017.00639

PubMed Abstract | CrossRef Full Text | Google Scholar

Hassler, U., and Thadewald, T. (2010). Nonsensical and biased correlation due to pooling heterogeneous samples. J. R. Stat. Soc. 52, 367–379.

Google Scholar

Holma, K. M., Haukka, J., Suominen, K., Valtonen, H. M., Mantere, O., Melartin, T. K., et al. (2014). Differences in incidence of suicide attempts between bipolar I and II disorders and major depressive disorder. Bipolar Disord. 16, 652–661.

Google Scholar

Hosseinifard, B., Moradi, M. H., and Rostami, R. (2013). Classifying depression patients and normal subjects using machine learning techniques and nonlinear features from EEG signal. Comput. Methods Programs Biomed. 109, 339–345.

Google Scholar

Jaworska, N., Blier, P., Fusee, W., and Knott, V. (2012). α Power, α asymmetry and anterior cingulate cortex activity in depressed males and females. J. Psychiatric Res. 46, 1483–1491. doi: 10.1016/j.jpsychires.2012.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Knott, V., Mahoney, C., Kennedy, S., and Evans, K. (2001). EEG power, frequency, asymmetry and coherence in male depression. Psychiatry Res. 106, 123–140.

Google Scholar

Kudela, M., Harezlak, J., and Lindquist, M. A. (2017). Assessing uncertainty in dynamic functional connectivity. NeuroImage 149, 165–177.

Google Scholar

LeMoult, J., and Gotlib, I. H. (2019). Depression: A cognitive perspective. Clin. Psychol. Rev. 69, 51–66.

Google Scholar

Li, G., Jiang, Y., Jiao, W., Xu, W., Huang, S., Gao, Z., et al. (2020). The maximum eigenvalue of the brain functional network adjacency matrix: Meaning and application in mental fatigue evaluation. Brain Sci. 10:92. doi: 10.3390/brainsci10020092

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Q., He, H., Yang, J., Feng, X., Zhao, F., and Lyu, J. (2020). Changes in the global burden of depression from 1990 to 2017: Findings from the Global Burden of Disease study. J. Psychiatric Res. 126, 134–140. doi: 10.1016/j.jpsychires.2019.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Melchior, M., Chastang, J. F., Head, J., Goldberg, M., Zins, M., Nabi, H., et al. (2013). Socioeconomic position predicts long-term depression trajectory: A 13-year follow-up of the GAZEL cohort study. Mol. Psychiatry 18, 112–121. doi: 10.1038/mp.2011.116

PubMed Abstract | CrossRef Full Text | Google Scholar

Nyström, C., Matousek, M., and Hällström, T. (1986). Relationships between EEG and clinical characteristics in major depressive disorder. Acta Psychiatr. Scand. 73, 390–394.

Google Scholar

Peraza, L. R., Asghar, A. U., Green, G., and Halliday, D. M. (2012). Volume conduction effects in brain network inference from electroencephalographic recordings using phase lag index. J. Neurosci. Methods 207, 189–199. doi: 10.1016/j.jneumeth.2012.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Pion-Tonachini, L., Hsu, S. H., Chang, C. Y., Jung, T. P., and Makeig, S. (2018). Online automatic artifact rejection using the real-time EEG source-mapping toolbox (REST). Annu. Int. Conf. IEEE Eng. Med. Biol. Soc. 2018, 106–109. doi: 10.1109/EMBC.2018.8512191

PubMed Abstract | CrossRef Full Text | Google Scholar

Plis, S. M., Sui, J., Lane, T., Roy, S., Clark, V. P., Potluru, V. K., et al. (2014). High-order interactions observed in multi-task intrinsic networks are dominant indicators of aberrant brain function in schizophrenia. NeuroImage 102(Pt 1), 35–48. doi: 10.1016/j.neuroimage.2013.07.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Spitzer, R. L., Kroenke, K., and Williams, J. B. (1999). Validation and utility of a self-report version of PRIME-MD: The PHQ primary care study. Primary Care Evaluation of Mental Disorders. Patient Health Questionnaire. JAMA 282, 1737–1744.

Google Scholar

Stam, C. J., Nolte, G., and Daffertshofer, A. (2007). Phase lag index: Assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources. Hum. Brain Mapp. 28, 1178–1193. doi: 10.1002/hbm.20346

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, S., Li, X., Zhu, J., Wang, Y., La, R., Zhang, X., et al. (2019). Graph theory analysis of functional connectivity in major depression disorder with high-density resting state EEG data. IEEE Trans. Neural Syst. Rehabil. Eng. 27, 429–439. doi: 10.1109/TNSRE.2019.2894423

PubMed Abstract | CrossRef Full Text | Google Scholar

Tibshirani, R. (2011). Regression shrinkage and selection via the lasso: A retrospective. J. R. Stat. Soc. Ser. B (Statistical Methodology) 73, 267–288.

Google Scholar

Wang, N., Zhang, L., and Liu, G. (2015). EEG-based research on brain functional networks in cognition. Biomed. Mater. Eng. 26(Suppl. 1), S1107–S1114.

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

Yao, D. (2001). A method to standardize a reference of scalp EEG recordings to a point at infinity. Physiol. Meas. 22, 693–711.

Google Scholar

Yu, M., Hillebrand, A., Tewarie, P., Meier, J., van Dijk, B., Van Mieghem, P., et al. (2015). Hierarchical clustering in minimum spanning trees. Chaos (Woodbury, N.Y.) 25:023107.

Google Scholar

Zhang, B., Yan, G., Yang, Z., Su, Y., Wang, J., and Lei, T. (2021). Brain functional networks based on resting-state EEG data for major depressive disorder analysis and classification. IEEE Trans. Neural Syst. Rehabil. Eng. 29, 215–229. doi: 10.1109/TNSRE.2020.3043426

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. Alzheimers Dis. 54, 1095–1112. doi: 10.3233/JAD-160092

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Shen, J., Din, Z. U., Liu, J., Wang, G., and Hu, B. (2019). Multimodal depression detection: Fusion of electroencephalography and paralinguistic behaviors using a novel strategy for classifier ensemble. IEEE J. Biomed. Health Inform. 23, 2265–2275. doi: 10.1109/JBHI.2019.2938247

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, F., Chen, Z., Rekik, I., Lee, S. W., and Shen, D. (2020). Diagnosis of autism spectrum disorder using central-moment features from low- and high-order dynamic resting-state functional connectivity networks. Front. Neurosci. 14:258. doi: 10.3389/fnins.2020.00258

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, F., Zhang, H., Rekik, I., An, Z., 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