Hybrid Functional Brain Network With First-Order and Second-Order Information for Computer-Aided Diagnosis of Schizophrenia

Brain functional connectivity network (BFCN) analysis has been widely used in the diagnosis of mental disorders, such as schizophrenia. In BFCN methods, brain network construction is one of the core tasks due to its great influence on the diagnosis result. Most of the existing BFCN construction methods only consider the first-order relationship existing in each pair of brain regions and ignore the useful high-order information, including multi-region correlation in the whole brain. Some early schizophrenia patients have subtle changes in brain function networks, which cannot be detected in conventional BFCN construction methods. It is well-known that the high-order method is usually more sensitive to the subtle changes in signal than the low-order method. To exploit high-order information among brain regions, we define the triplet correlation among three brain regions, and derive the second-order brain network based on the connectivity difference and ordinal information in each triplet. For making full use of the complementary information in different brain networks, we proposed a hybrid approach to fuse the first- and second-order brain networks. The proposed method is applied to identify the biomarkers of schizophrenia. The experimental results on six schizophrenia datasets (totally including 439 patients and 426 controls) show that the proposed method outperforms the existing brain network methods in the diagnosis of schizophrenia.

Brain functional connectivity network (BFCN) analysis has been widely used in the diagnosis of mental disorders, such as schizophrenia. In BFCN methods, brain network construction is one of the core tasks due to its great influence on the diagnosis result. Most of the existing BFCN construction methods only consider the first-order relationship existing in each pair of brain regions and ignore the useful high-order information, including multi-region correlation in the whole brain. Some early schizophrenia patients have subtle changes in brain function networks, which cannot be detected in conventional BFCN construction methods. It is well-known that the high-order method is usually more sensitive to the subtle changes in signal than the low-order method. To exploit high-order information among brain regions, we define the triplet correlation among three brain regions, and derive the second-order brain network based on the connectivity difference and ordinal information in each triplet. For making full use of the complementary information in different brain networks, we proposed a hybrid approach to fuse the first-and second-order brain networks. The proposed method is applied to identify the biomarkers of schizophrenia. The experimental results on six schizophrenia datasets (totally including 439 patients and 426 controls) show that the proposed method outperforms the existing brain network methods in the diagnosis of schizophrenia.

INTRODUCTION
Resting-state functional magnetic resonance imaging (rs-fMRI) studies indicate that there exists a disorder-related alteration in BFCN (Bluhm et al., 2007;Jafri et al., 2008;Fornito and Bullmore, 2010;Shafiei et al., 2018;Wang et al., 2018). As a severely debilitating mental illness, schizophrenia is usually thought of as connectivity disorders between brain regions (Liang et al., 2006;Salvador et al., 2010). In recent years, many BFCN analysis methods Liang et al., 2007;Honey et al., 2009;Tsuang et al., 2009) have been proposed to explore the biomarkers for schizophrenia. Most of the existing methods first construct the brain functional network by measuring the correlation of brain regions or voxels, and then perform feature extraction on the produced large-scale brain networks for selecting the significant features. These connectivities or sub-networks having significant alterations in some indicators, e.g., topological metrics (Fei et al., 2014) and the alteration degree (Guo et al., 2014;Zhu et al., 2018), are chosen as the biomarker for the disease.
The construction of the brain network plays as an important role in the diagnosis system. The previous BFCN construction methods mainly focus on revealing the low-order information among brain regions or voxels. In other words, the functional connectivities in these methods are usually denoted as the association between each pair of brain regions or voxels; e.g., the Pearson correlation (PC) based method uses the pair-wise linear correlation of brain regions as the connectivity strength (Richiardi et al., 2013). In addition, the Kendall correlation (KC) (Dong et al., 2014) and the Spearman correlation (SC)  are also introduced into the brain network construction process. Considering the correlation between the two brain regions may be affected by other regions, Guo et al. (2014) proposed to eliminate this kind influence from other connectivities through the partial correlation (ParC) test, and applied it to the classification of schizophrenia patients and healthy controls. Li et al. (2019) and Qiao et al. (2016) constructed the low rank and self-weighted brain network by introducing prior knowledge to the model. Zhou et al. (2018) embedded the second-order information among brain regions into brain network and applied it into identifying mild cognitive impairment.
The low-order method has the following natural defects in characterizing the relationship among brain regions. On one hand, there exists the subtle alteration of brain network structure between controls and patients, especially the mild or early patients, which is hard to reveal through the first-order method because, compared to the high-order method, the low-order method is not sensitive to small signal change. On the other hand, the first-order or pair-wise model is not capable of describing the complex multivariable correlation and ignores the ordinal pattern among connectivities (Wu et al., 2018). The topological structure of the brain network is intricate (Bullmore and Sporns, 2009); e.g., one brain region usually interacts with multiple brain regions physiologically, whereas the low-order methods often ignore the relationship among multi-regions. Therefore, it is necessary to investigate effective high-order brain network construction methods to reveal the correlation of multiple brain regions and ordinal patterns from rs-fMRI data.
In this paper, we derived the second-order brain networks by introducing the triplet of brain regions. We calculated the connectivity difference to describe the ordinal information in each triplet. The derived second-order relationship in BFCN construction methods is more capable of capturing the changes in the brain network, making it easier to identify mild or early schizophrenia patients. Considering that the brain network has the property of high regional agglomeration (Alexander-Bloch et al., 2013), we extract the second-order relationship among one brain region and its k nearest neighbors to improve the robustness of the proposed method. Then, we fuse our secondorder BFCN with the most widely used first-order BFCN, i.e., PC based brain network, to construct the hybrid functional brain network. The mixed model not only reveals the subtle differences in brain networks between patients and healthy controls, but also makes full use of complementary brain networks defined in different orders. Our proposed method has the following advantages: 1) We embed high-order information of connectivity groups in the brain network, which may have a beneficial effect on the analysis and diagnosis of mild or early mental illness due to the revelation of a more subtle alteration. 2) Considering the regional aggregation of brain networks, we utilize the local strategy to improve the robustness of the proposed second-order brain network. 3) A hybrid strategy fusing the proposed second-order brain network and the conventional first-order brain network is proposed, which can make full use of the complementary information in these two types of networks. 4) The experiment is performed on six rs-fMRI datasets with schizophrenia patients and controls, including a total of 865 subjects. To the best of our knowledge, it is with the largest population reported in the literature of rs-fMRI analysis for identifying schizophrenia. The results show that the proposed hybrid network is superior to the state-of-the-art methods.
The remainder of the paper is organized as follows: first, we introduce our materials and methods, i.e., the six schizophrenia datasets and the proposed brain network construction method. Then, the experimental results and discussion on six schizophrenia datasets are presented. Finally, we summarize our work and provide the conclusion.

Subjects
Six schizophrenia rs-fMRI datasets we used are Nanjing Brain Hospital (NBH) dataset, The Center for Biomedical Research Excellence (COBRE) dataset, Huaxi dataset, Nottingham dataset, Taiwan dataset, and Xiangya dataset. The subjects have the following requirements: (1) no other Diagnostic and Statistical Manual of Mental Disorders (DSM-IV) disease exists; (2) no history of drug abuse; (3) no clinically significant head trauma. Table 1 summarizes the demographic and clinical characteristics of participants of these datasets (Cheng et al., 2015). The Positive and Negative Symptom Scale (PANSS) (Kay et al., 1987) is used by experts to score participants to obtain the severity of schizophrenia.

Image Acquisitions, Data Preprocessing and Anatomical Parcellation
For NBH, COBRE, Taiwan, and Xiangya datasets (Guo et al., 2014), the rs-fMRI images of each participant were acquired by using a 3-Tesla Siemens Tim-Trio scanner with an eight or 12 channel head coil. For Huaxi dataset, the rs-fMRI images were acquired by using a 3-T General Electric MRI scanner. For the Nottingham dataset, the images were acquired using a 3-T Philips Achieva MRI scanner. During the scanning of images in the NBH, COBRE, and Nottingham datasets, all participants must keep their eyes open and stare intently at the fixed cross in the middle of the screen for 5 min. During scanning of images in the Huaxi, Taiwan, and Xiangya datasets, the participants were instructed to keep their eyes closed but not fall asleep. For all datasets to be preprocessed, we discarded the first 10 volumes to make sure scanner stabilization and the subjects' adaptability to the environment. Rs-fMRI data preprocessing was then performed by statistical parametric mapping (SPM8) software (http://www.fil.ion.ucl.ac.uk/spm) and a Data Processing Assistant for RS-fMRI (DPARSF). The remaining functional scans were first corrected to the difference in acquisition time between on-chip scans, and then readjusted to the intermediate volume to correct for head movement between scans. The functional scan was then spatially normalized to the Montreal Neurological Institute template and resampled to 3 × 3 × 3 mm 3 voxels. After linear detrending, data was filtered using typical temporal bandpass (0.01-0.08 Hz) to reduce low-frequency drift and high-frequency physiological noise. Next, four covariates of no interest (i.e., six rigid-body motion parameters, the global mean signal, white matter, cerebrospinal fluid) were regressed from the data. No subject was excluded under a head motion criterion of 3 mm and 3 • .
After data preprocessing, according to Lynall et al. (2010), the rs-fMRI images were divided into 90 brain regions (excluding the cerebellum) or region-of-interests (ROIs) using the Anatomical Automatic Labeling (AAL) template (Tzourio-Mazoyer et al., 2002). Finally, for each ROI, the average rs-fMRI time series over all voxels was taken as the time series for that ROI. After data preprocessing and anatomical parcellation, for each of the brain regions of each subject, its feature information is represented in the form of the time series. The processed information of all datasets is shown in Table 2.

Triplet-Based Second-Order Functional Brain Network
The traditional functional brain network methods often only consider the pairwise relationship of different brain regions, namely the first-order relationship. However, the alteration of connectivity networks between patients and healthy controls may be subtle, and the first-order brain network is not capable of finding the subtle difference. In addition, the topology of the brain is extremely complex and one brain region usually interacts with multiple brain regions. For revealing the highorder information among brain regions, we construct the secondorder brain network by introducing triplet correlation among three brain regions. Suppose that X = {x 1 , x 2 , . . . , x n } ∈ R n×m denotes the rs-fMRI data from a subject, where n is the number of brain regions and m is the number of time series. x i is the time series for the i-th brain region where i ∈ {1, 2, · · · , n}. Our triplet-based second-order functional brain network is constructed as follows: consists of x i and its neighbors x u and x v . S i uv defines the distance between x i and x v relative to x u : where dist (·, ·) denotes the squared Euclidean distance. It should be pointed out that dist (·, ·) has nothing to do with the anatomical distance, because x i denotes the rs-fMRI time series data of the i-th brain region. Obviously, S i is an antisymmetric matrix. Since one brain region usually interacts with its neighbors rather than distant brain regions, we only consider the k nearest neighbors of x i . Let N i be a set of sequence numbers indicating the k nearest neighbors of x i , and then, relative to all k nearest neighbors of x i , the distance between x i and x v can be expressed as: Based on the above distance, the triplet correlation coefficient (Guo et al., 2017) between x i and x v can be expressed as: where norm (·) denotes normalizing the data to the interval [0, 1]. Thus, the second-order correlation coefficient matrix can be denoted as: Then, we symmetrically handle C as follows: The schematic diagram of triplet-based local brain network is shown in Figure 1.
After all correlation coefficients are calculated, a second-order functional brain network is constructed. It preserves secondorder information among brain regions, which makes it more sensitive to the changes in the brain network and can capture the ordinal information among connectivities. It can also be seen from the above formula that our second-order functional brain network will eventually be converted into a two-dimensional form, which also helps to reduce the space complexity of highorder method.

Hybrid Functional Brain Network With First-Order and Second-Order Information
Studies have shown that brain network methods based on correlation coefficients can characterize the interactions between brain regions or neurons (Smith et al., 2011). The conventional first-order method and the proposed second-order method can construct connectivities from different levels. There is some complementary information in these two kinds of brain networks for the diagnosis task. In addition, the low-order method is more robust to noise transmission. Inspired by the above issue, we fuse the proposed brain network with a first-order brain network to construct a hybrid brain network. The PC method is a classical first-order method for constructing brain networks based on correlation coefficients. PC can be calculated using the following formula: where Var (·) means a function of the variance calculated and Cov (·, ·) is a function that calculates the covariance. For convenience, we denote the proposed second-order functional brain network C ′ C ′ mentioned in the previous section FIGURE 1 | The schematic diagram of triplet-based local brain network. The black point x i is the center point we selected. We calculated the distance between x i and its neighbors in connectivity strength. There are three triples based on x i and its neighbors in calculation. The last sub-graph is our triplet-based second-order functional network in the whole brain.
Frontiers in Neuroscience | www.frontiersin.org FIGURE 2 | Overall framework of the hybrid brain network. The first-order brain network based on the Pearson correlation coefficient and the second-order brain network based on the triplet are, respectively, constructed, and then the two are fused with a certain weight to obtain the hybrid functional brain network with firstand second-order information.
Frontiers in Neuroscience | www.frontiersin.org TABLE 3 | Classification results (ACC/SEN/SPE/PPV/NPV/F1/BAC/AUC±STD%, and p-value) by k-fold cross-validation of several Functional fMRI networks algorithms on schizophrenia datasets. Frontiers in Neuroscience | www.frontiersin.org as C 2 . In this case, hybrid functional brain network with firstorder and second-order information can be expressed as: where µ is the weighted coefficient for mixing C 1 and C 2 , and its ranges from 0 to 1. The overall framework of the hybrid brain network is shown in Figure 2.
In the hybrid brain network, C 1 preserves the original connectivity strength information between the two brain regions, and C 2 captures the connectivity strength information of one brain region and its neighbors. The hybrid brain network C synthesizes C 1 and C 2 , and fuses the complementary information in these two kinds of networks.
By introducing the triplet correlation, the ordinal relationship among brain regions is preserved. Through second-order brain network, i.e., the triplet-based brain network, not only highorder information is introduced to extract neighbor ranking information in each brain region, but also connectivity noise is transmitted. From the point of view of signal analysis, the first-order brain network is robust to noise, and the secondorder brain network is sensitive to the subtle differences between the brain network of patients and healthy controls. Thus, the proposed hybrid model can balance robustness and classification performance.

Feature Selection Strategy
There are mainly three stages of mental illness classification based on BFCN: brain network construction, brain network feature selection, and classification. Besides the construction of functional brain networks, feature selection is also a key step in the brain network analysis process. Since this work focuses on brain network construction, the discussion of feature selection algorithm is beyond the research content of this paper. By comparing the classification performance of multiple feature selection algorithms, we use the non-negative elasticnet  in our previous work to choose the significant connectivities.
The functional brain network data are high dimensional, and many subjects are linearly inseparable. To achieve good classification performance, we utilize kernel discriminant analysis (KDA) algorithm to project each sample to the feature space, in which the classification performance of subjects can be improved. Finally, the nearest neighbor classifier is exploited for making the diagnosis decision.
In addition, the topological properties of the brain network can also be used as important features. Similarly, since this work focuses on the construction of the brain network, we use the methods proposed by Narula et al. (2017) and the public code they provided to calculate the topological properties of the brain network.

Classification Performance of Different fMRI Functional Brain Networks
In this experiment, the brain networks were constructed using the proposed method and a series of comparison methods mentioned below. Since this experiment is designed Frontiers in Neuroscience | www.frontiersin.org for comparing the performance of different functional brain networks, the same feature selection and classification algorithm in our method are all performed on these networks to ensure the comparability of the functional brain networks.

Experimental Setup
We adopt the 10-fold cross-validation strategy in this experiment. Specifically, each dataset is equally divided into 10 subsets. One subset is selected as the test set successively, and all remaining subsets are used for training. To avoid possible deviations during sample segmentation, we repeat this process 20 times. Accuracy (ACC), sensitivity (SEN), specificity (SPE), positive TABLE 4 | Classification results (ACC/SEN/SPE/PPV/NPV/F1/BAC/AUC±STD%, and p-value) by k-fold cross-validation of several feature selection algorithms on schizophrenia datasets. predictive value (PPV), negative predictive value (NPV), F1score (F1), class balanced accuracy (BAC), the area under the receiver operating characteristic curve (AUC) and their respective standard deviations (STD) are employed to measure the performance in classification. They can be calculated as follows: where TP, TN, FP, and FN are the number of correctly predicted patients, correctly predicted healthy controls, healthy controls predicted as patients, and patients predicted as healthy controls, respectively. N pat and N nor are the number of patient samples and the number of healthy control samples, respectively. rank i represents the serial number of the i-th sample. In addition, we calculated the p-values (DeLong et al., 1988) of the AUC to measure statistically significant differences between proposed method and comparison methods. If the p-value is <0.05, it indicates that the increase in performance of our method compared to the comparative method is statistically significant . For fair comparison, for our method, we adopt grid search to select the parameter of the number of neighbors from The asterisk (*) denotes the statistically significant differences (p < 0.05).

Results and Analysis
We report the comparison results of classification based on different functional brain networks on different datasets in Table 3. In particular, we also report the p-values and mark statistically significant differences (p < 0.05) with the asterisk ( * ). The receiver operating characteristic (ROC) curves are shown in Figure 3. It can be seen that our BFCN construction method has the best diagnostic results on six schizophrenia datasets among all brain network construction methods. Discussion of experimental results is given in section Discussion.

Classification Performance of Different Dimension Reduction and Classification Methods
This experiment is designed to test the performance of our brain network under different dimensionality reduction and classification methods. In this experiment, the brain network is constructed using the proposed method. Some widely used feature extraction and classification methods in functional brain network analysis are compared.

Experimental Setup
Similar to the previous experiment, we adopt 10-fold crossvalidation strategy and repeated this process 20 times. ACC, SEN, SPE, PPV, NPV, BAC, AUC, and their STD are employed to measure the performance in classification. P-value of the AUC is used to measure statistically significant differences between KDA and other feature extraction methods. For fair comparison, for the generation of functional brain networks we adopt a grid search to select the parameter of the number of neighbors from {5, 10, 15, · · · , 45}, and the weighted coefficient of the hybrid brain network is selected from {0, 0.1, 0.2, · · · , 1}. We set threshold from {0, 0.05, 0.1, · · · , 0.4} for hybrid brain network on all datasets. For the feature selection algorithm, we use Gaussian kernel function in KDA. For the similarity-based graph methods, including ISOMAP and NPE, we build an edge between two subjects if and only if they belong to the same class. For SVM, we adopt linear kernel.
FIGURE 4 | SAC between schizophrenia and healthy controls (Top 10). Brain network connectivities with the top 10 discrimination ability. In (A), the position in the i-th row and j-th column indicates connectivity between the i-th brain region and the j-th brain region, and the weight can be judged by the corresponding color.

Results and Analysis
We report the comparison results of classification based on different feature selection and classification methods on different datasets in Table 4. In particular, we mark statistically significant differences (p < 0.05) with the asterisk ( * ). It can be seen that our method has the best ACC and AUC performance on six datasets compared to other feature selection and classification algorithms. Discussion of experimental results is given in section Discussion.

DISCUSSION
In the first experiment, we compared the performance of our functional networks and other functional networks in ACC, SEN, SPE, PPV, NPV, BAC, and AUC on six schizophrenia datasets and reported the p-value of the AUC between proposed method and comparison methods. Experimental results have proved that our functional network is effective and superior to other functional networks in all the measures and the results are significantly different. In addition, we draw ROC curves of our method and all comparison methods on each dataset. It can be seen that the ROC curves of the comparison methods are almost at the bottom right of the ROC curve of our method, and the results of the area under the curve of the comparison methods are also smaller than our method. The reason why our method can show the above performance may be that our network not only considers the relationship between the two brain regions, but also preserves the second-order information among brain regions. Firstorder information has a certain robustness, and second-order evaluation indicators are ACC, SEN, SPE, PPV, NPV, BAC, and AUC. In addition, p-values of the AUC between proposed method and comparison methods are also reported. The experiment proves that our hybrid brain network is robust and performs well under each classifier. In addition, the experimental results show that our framework performs well, and the ACC and AUC indicators have achieved the best results on the six schizophrenia datasets. This may be due to the complexity of the brain network, which is not simply linearly separable, and our framework chose the KDA classifier, which can select the most significant connectivities and improve the classification performance.
In addition to classification performance, there are many other indicators that evaluate the functional brain networks. The topology metric of the network is a widely used one. We compared the performance of our network and PC based method with some topological properties, and calculated the corresponding p-value, including connectivity strength, average FIGURE 5 | Significant alteration of brain regions between schizophrenia and healthy controls. degree, density, clustering coefficient, characteristic path length, global efficiency, local efficiency, closeness centrality, edge betweenness centrality, node betweenness centrality, radiality, assortativity, structural consistency, and fitted exponent of power-law degree distribution. The brain network that calculates the topological properties consists of the brain networks that obtained the best classification performance in Experiment 3.1. The experimental results are shown in Table 5. The asterisk ( * ) denotes the statistically significant differences (p < 0.05).
The experimental results show that, compared to the PC based method, the brain network based on our functional brain network has better separability in these topological properties between the patient group and the healthy group. Specifically, we found that based on our proposed network, the connectivity strength of the health group is 32.3520 and that of the patient group is 31.5238. The functional network connectivity strength of the patient group is reduced by 0.8282, and there are similar findings in the literature (Shen et al., 2010). We also found that the clustering coefficient of the patient group was lower than that of the healthy group, which was also supported by the work of (Bachiller et al., 2014). Also, the global efficiency of patients with schizophrenia has decreased, indicating that the topology of the brain structure network in schizophrenia is less efficient. Griffa et al. (2014) has also mentioned this point.
The topological properties of the network show the differences in the overall brain network between the patient group and the health group. We also show the local difference with 30 connectivities that have the greatest differences in the patient group and the health group by calculating the significant alteration of connectivity (SAC). The results are shown in Table 6. To visualize these connectivities, we show the top 10 SAC drawn in Figure 4.
As shown in Table 6 and Figure 4, the connectivity between Precuneus and other brain regions altered visibly, and (Faget-Agius et al., 2012) supports our finding. The connectivity between Hippocampus and other brain regions and the connectivity between Inferior frontal gyrus also have visible alteration, and there exist similar findings in some literature (Altshuler, 1990;Zhou et al., 2008;Kubicki et al., 2011). From Figure 4 and Table 6, it also can be seen that some brain regions are recurring, indicating that these brain regions play a key role in the analysis of differences between the patient group and the healthy group. Therefore, for each recurring brain region, we accumulate the weights of the SACs, where the brain region is located, as the weight of the brain region. All weights are normalized and sorted, and we list the top 15 important brain regions in Table 7 and visualize them in Figure 5.
As shown in Table 7 and Figure 5, some brain regions such as Precuneus, Hippocampus and Rolandic operculum are selected. This suggests that these brain regions play an important role in the task of identifying patients with schizophrenia and normal controls, and that these brain regions may be important biomarkers for schizophrenia. Some literatures (Bettus et al., 2009;Salvador et al., 2010;Faget-Agius et al., 2012;Qiu et al., 2018) also prove our findings.
Furthermore, we constructed the connectivity subnets with these 15 brain regions in Table 7. The average patient group subnet, average health group subnet, the difference between average patient group subnet and average health group subnet are shown in Figure 6. It can be seen that although the strength of the patient group connectivity on the overall network is lower than that of the healthy group, the average connectivity strength of the patient group in the subnet is slightly higher than that of the healthy group. This change further demonstrates the importance of the brain regions we have discovered in identifying tasks in patients with schizophrenia, that is, these brain regions may play an important role in schizophrenia classification.

CONCLUSION AND FUTURE WORK
In summary, we propose a hybrid functional brain network with first-order and second-order information for identifying schizophrenia. Specifically, we construct a second-order brain network through triplet correlation, and fuse it with conventional first-order brain network. The second-order brain network is more sensitive to the difference in brain networks between patients and healthy controls, and the first-order brain network is more robust to noise. The proposed method not only captures higher-order information among brain regions, but also reveals the ordinal information of connectivity strength between specific brain regions. Experiments on six schizophrenia datasets show that our method is superior to the existing BFCN construction method. In addition, we analyzed the differences in topological properties between schizophrenia patients and normal controls with the proposed brain network.
This study uses the grid search method in the process of constructing the hybrid functional brain network with first-and second-order information. We will investigate more efficient parameter selection methods in future work. In addition, this study is conducted separately on a dataset from a single site. In our future work, we will focus on how to construct and analyze the high-order brain network on the cross-site dataset, which is helpful to improve the robustness of the model.

AUTHOR CONTRIBUTIONS
QZ conceived the experiment and revised the manuscript. HL was responsible for the experiment and writing the first draft of the manuscript. JH and DG were responsible for preprocessing the fMRI data. XX gave clinical guidance for biomarkers. DZ contributed to data analysis. All authors read and endorsed the final draft.