Multi-graph Fusion for Functional Neuroimaging Biomarker Detection

1 Introduction

Functional magnetic resonance imaging (fMRI) characterizes brain activity by detecting the synchronized time-dependent changes of the blood oxygenation level dependent (BOLD) signals. Recently, fMRI data has been becoming one of popular sources to improve neuro-disease diagnosis because neuroimaging biomarker detection with fMRI data has the potentiality to comprehensively understand neurological disorders at a whole-brain level [Shu et al., 2019a].

Given BOLD signals, a functional connectivity network (FCN) is constructed for each subject. Usually, a FCN is represented by a symmetric matrix, where each element implies the correlation of the BOLD signals between two nodes (i.e., brain regions) and is calculated by either Pearson analysis methods or wavelet correlation methods [Shu et al., 2019b]. After this, two steps are designed for conducting neuro-disease diagnosis with fMRI data, i.e., representation learning and disease diagnosis (i.e., classification). Representation learning is designed to fine-tune the full FCN, where each node connects all nodes and the value of each connectivity represents the correlation between two nodes. Disease diagnosis usually employs existing methods to conduct classification tasks on the representations of all subjects.

� Corresponding author

In the process of representation learning, full FCN methods (e.g., [Karmonik et al., 2019]) are designed to extract the upper triangle of the symmetric matrix (i.e., the full FCN) to represent the subject by a vector. Full FCNs have been verified being vulnerable to false or irrelevant functional connectivity [Kong et al., 2015; Zille et al., 2017]. Therefore, sparse FCN methods [Li et al., 2017; Zhang et al., 2019a] are designed to connect each node to a part of nodes to possibly remove unimportant functional connectivity. For example, [Eavani et al., 2015] and [Zille et al., 2017] proposed to directly transfer the dense matrix representation in the full FCNs to a sparse matrix. Furthermore, a number of studies employ traditional classifiers (e.g., support vector machine (SVM) and logistic regression) to conduct neuro-disease diagnosis. To avoid the issue of the curse of the dimensionality on high-dimensional data, previous methods of disease diagnosis usually conduct dimensionality reduction before the classification tasks [Zhang et al., 2019b; Zhang et al., 2017].

Previous FCN methods have a number of issues to be addressed due to all kinds of reasons, such as inter-subject variability, heterogeneity across subjects, and discriminative ability. First, previous sparse FCN methods (e.g., [Wee et al., 2012]) often make the assumption that every node has the same connectivity number. Actually, human brain is a complex system and human brain contains the inter-subject variability where every subject or every node within one subject has individual characteristics. The inter-subject variability makes the assumption of equivalent connectivity number unreasonable. Moreover, it is di�cult to decide the connectivity number for each node in real applications because we usually have litter prior knowledge about the brain functional connectivity. Second, existing FCN methods (e.g., [Li et al., 2017; Eavani et al., 2015]) ignore the heterogeneity across subjects for representation learning. Specifically, they generate the representation of each subject independent on other subjects without taking the group e�ect into account. In practice, different subjects may be obtained from di�erent places or operated by di�erent doctors, and thus have di�erent distributions. Third, the independent process for representation learning ignores to consider the group e�ect so that the outputted representation has limited discriminative ability.

In this paper, we propose a functional connectivity analysis framework to conduct representation learning and personalized disease diagnosis on fMRI data in a semi-supervised manner. Specifically, we first propose a multi-graph fusion method to generate homogeneous and discriminative representations for all subjects, and then employ '1-SVM to conduct joint brain region selection (i.e., feature selection) and disease diagnosis (i.e., classification). In the multi-graph fusion method, we employ Pearson correlation analysis to output a full FCN as well as an extremely sparse FCN for every subject, denoted two FCNs as multi-graph in this paper. We use the obtained multi-graph to automatically learn a sparse FCN for each subject where di�erent nodes have different connectivity numbers and the subjects within the same class have maximal similarity while the subjects with di�erent class labels have maximal dissimilarity.

2 Method

In this paper, we denote matrices, vectors, and scalars, respectively, as boldface uppercase letters, boldface lowercase letters, and normal italic letters. Given the BOLD signal of the m-th subject among M subjects B m 2 R n�t (m = 1; :::; M) where n and t, respectively, represent the number of brain regions and the length of signals, in this paper, we first obtain multiple graphs (i.e., FCNs) A m;v 2 R n�n (v = 1; :::; V) by Pearson correlation analysis where V is the graph number, and then propose to learn a sparse FCN S m for each subject so that it could automatically learn the connectivity number of every node as well as is homogenous and discriminative to other sparse FCNs S m 0 (m , m 0 ).

2.1 Multi-graph Fusion

Previous studies demonstrated that the sparse FCN is preferred in representation learning of brain function connectivity analysis ([Karmonik et al., 2019]), compared to the full FCN, duo to that 1) the full FCN lacks interpretability; 2) the connectivity between two nodes may contain noisy connectivity (i.e., either irrelevant or spurious connectivity) to a�ect brain functional connectivity analysis [Whitwell and Josephs, 2012]; and 3) neurologically, a brain region predominantly interacts only with a part of brain regions. Existing methods of functional connectivity analysis usually obtain sparse FCNs from the full FCNs. Specifically, previous methods design di�erent techniques to learn sparse FCNs based on the full FCNs, such as sparse learning [Zhang et al., 2019a; Eavani et al., 2015] and clustering [Zhang et al., 2019b]. However, previous methods have limitations in brain functional connectivity analysis.

First, existing methods usually assume that each node connects a fixed number of nodes out of all nodes. To achieve this, the sparse k-nearest neighbor (kNN) graph is constructed so that each node connects with k nodes. Such an assumption obviously ignores the fact that a brain region predominantly interacts only with a part of brain regions. Second, previous methods generate the sparse FCN of a subject independent from other subjects. On one hand, by considering the heterogeneity across subjects, the FCNs obtained from these heterogenous subjects possibly have di�erent distributions. On the other hand, the independent process of representation learning makes it di�cult to consider the group e�ect, e.g., the discriminative ability across classes or subjects.

Given the full FCN connecting each node with all nodes, we obtain an extreme sparse FCN, i.e., 1NN graph (excluding itself). By this way, we could obtain multiple graphs for each subject to solve the first issue of existing functional connectivity analysis. Moreover, in this paper, we only use 2 graphs for every subject, i.e., a full FCN and an extremely sparse FCN. The full FCN contains all connectivity information (i.e., the most complex connectivity) and the extremely sparse FCN contain the least information (i.e., the simplest connectivity). We expect to obtain a flexible connectivity number for every node based on the data distribution in the range [1; n] where n is the node number. To do this, we design the following objective function to automatically learn specific connectivity number for the m-th subject S m by fusing the information from multiple graphs.

where k � kF indicates Frobenius norm. s m i;� and s m i; j , respectively, represent the i-th row of S m and the element in the i-th row and the j-th column of S m . 1 and N(i), respectively, indicate the all-one-element vector and the set of nearest neighbors of the i-th node. The constraint s m T i;� 1 = 1 keeps the shift invariant similarity. After optimizing s m i;� by our proposed optimization method in Section 2.3, we could obtain di�erent non-zero numbers for every row, i.e., s m i;� in S m . This indicates that di�erent nodes have di�erent connectivity numbers for every subject.

Eq. (1) employs multiple graphs to conduct representation learning, aim at selecting an optimal connectivity number between 1 and n. However, the optimization of S m is independent on the optimization of S m 0 (m , m 0 ), which explores the inter-subject variability, but does not touch the issue of the heterogeneity across subjects. To address this issue, we propose the following objective function.

where H and G are two variables, R1(H;G) and R2(S 1 ; :::; S M ) are regularization terms. We use the summation operator in the first term of Eq. (2) to learn the representations of all subjects in a unified framework, and design two regularization terms to achieve the group e�ect, e.g., discriminative ability across subjects.

First, we expect that positive subjects are similar or close to the positive template G while negative subjects are similar to the negative template H. Hence, the subjects within the same class are close. Moreover, the outputted templates could be widely applied in medical imaging analysis, such as guiding parcellations for new subjects and measuring the group di�erence [Reyes et al., 2018]. To achieve this, we design R1(H;G) as follows

where D, E, and U, respectively, represent the set of negative subjects, positive subjects, and unlabeled subjects. Moreover, jDj and jEj, respectively, indicate the cardinality of D and E.

Eq. (3) has at least two advantages: 1) preserving the global structure since all the subjects are close to their template and 2) outputting practical templates. However, Eq. (3) does not take the local structure of the data, which has been regarded as the complementary of the global structure of the data [Wang et al., 2017; Yang et al., 2015]. In this paper, we design R2(S 1 ; :::; S M ) as follows

where G(i) and F (i), respectively, are the set of near-neighbor and the set of distant-neighbor, of the i-th subject. In the proposed framework, i.e., semi-supervised learning, the training subjects include labeled subjects and unlabeled subjects, we denote the set G(i) of the i-th unlabeled subject as its k nearest neighbors including labeled subjects and unlabeled subjects, and the set G(i) of the i-th labeled subject as its k nearest neighbors with the same label to the i-th subject. We further define the set F (i) of the i-th unlabeled subject as its k furthest subjects including labeled subjects and unlabeled subjects, and the F (i) of the i-th labeled subject as its k nearest neighbors with di�erent labels to the i-th subject. It is noteworthy that the value of k is insensitive in our experiments, so we fixed k = 10 for all subjects.

Eq. (4) minimizes the ratio of two terms, similar to linear discriminative analysis [Shen et al., 2015]. Specifically, the subjects have the same label with their nearest neighbors, while the subjects with far similarity have di�erent labels. In this way, the local structure of the subjects is preserved. The optimization of Eq. (4) is very challenging, so we follow [Shen et al., 2020] to convert the minimization of Eq. (4) to minimize the following objective function:

where � m can be updated as � m = P p2G(m) jjS m �S p jj 2 F P q2F (m) jjS m �S q jj 2 F in the implementation based on [Shen et al., 2020].

Compared to previous literature, Eq. (2) outputs the representation of every subject dependent on other subjects as well as taking into account the following constraints, such as multi-graph information and the preservations of the global as well as the local structure.

2.2 Joint Regions Selection and Disease Diagnosis

Our method generates a sparse FCN S m (m = 1; :::; M) from the multi-graph, i.e., a full FCN and a 1-NN graph, for each subject. Moreover, we follow previous methods to transfer the matrix representation to its vector representation, i.e., extracting the upper triangle of the symmetric matrix S m (m = 1; :::; M) to form a row vector xm;� 2 R 1�[n(n�1)=2] . In this way, we have the data matrix X 2 R M�[n(n�1)=2] and the corresponding label vector y 2 f�1; 1g M�1 .

Many existing studies separately conduct feature selection and disease diagnosis (i.e., classification). The goal of feature selection is to remove the redundant features from highdimensional data because the vector representation is a 4005dimensional vector for 90 nodes in our data sets. However, the optimal results of feature selection cannot guarantee the optimal classification in two separated processes. In this paper, we employ '1-SVM to simultaneously conduct feature selection and classification, where the result of feature selection will be iteratively updated by the optimized classifier so that outputting significant classification performance.

2.3 Optimization

In this paper, we employ the alternating optimization strategy [Shen et al., 2020] to optimize S m (m = 1; :::; M), H, and G, as well as list the pseudo of our optimization method in Algorithm 1.

The variables S 1 ; :::; S M include the representations of positive subjects, negative subjects, and unlabeled subjects, so we explain the optimization process one by one.

When m-th subject is a negative subject, we obtain the objective function with respect to S m as follows:

Algorithm 1 The pseudo of our proposed functional connectivity analysis framework.

C;

• 1: Obtain A m;v (v = 1; :::; V) by B m ;
• 2: Initialize S m as the average of A m;v (v = 1; :::; V);
• 3: while Eq:(2) not converges do
• 4: Update � m = P p2G(m) jjS m �S p jj 2 F P q2F (m) jjS m �S q jj 2 F ;
• 5: Update H and G via Eq. (14);
• 6: Update S m (m = 1; :::; M) via Eq. (11);
• 7: end while
• 8: Obtain X by extracting the upper triangle of S m ;
• 9: Run '1-SVM on X and y to output the classifier C;

The optimization of each row s m i;� (i = 1; :::; n) in S m is independent on other rows s m i 0 ;� (i , i 0 ), so the objective function with respect to s m i;� is:

After conducting mathematical transformation, we have

where

Based on the complementary slackness of the KarushKuhn-Tucker (KKT) conditions [Bertsekas, 1995], we have the closed-form solution for s m i; j

where f m � i; j is the j-th element of f m � i;� .

By following the same process from Eq. (6) to Eq. (10), we have

where

�1, �2 and �3 are the Lagrange multipliers. (ii)Update H and G by fixing S 1 ; :::; S M

When S 1 ; :::; S M are fixed, the objective function with respect to H and G are:

According to the KKT conditions, we have:

where ˆ s m � i; j = ( P m2D s m T i; j )=jDj, ˆ s m + i; j = ( P m2E s m T i; j )=jEj, �4 and �5

are Lagrange multipliers.

The values of the Lagrange multipliers �1, �2, �3, �4, and �5, can be obtained based on [Duchi et al., 2008]. For simplicity, we list the details of �3 as follows and the values of �1, �2, �4, and �5 can be obtained by similar principles.

2.4 Convergence, Initialization, and Complexity

The optimizations of the variables, such as S 1 ; :::; S M , H, and G, in Eq. (2), have close-form solutions. Moreover, Eq. (2) iteratively updates Eq. (11) and Eq. (14) based on the alternating optimization strategy [Shen et al., 2020], which has been proved to achieve convergence. Hence, the proposed multi-graph fusion model converges and '1-SVM achieves convergence based on [Yang et al., 2015].

In Algorithm 1, we initialize S m (m = 1; :::; M) as the average of A m;v (v = 1; :::; V), which makes the optimization of Eq. (2) converge within tens of iterations. Moreover, the result of Eq. (2) is insensitive to the initialization of S m (m = 1; :::; M).

The generation of multi-graph can be finished o�ine. Hence, we ignore to calculate the time complexity and the space complexity. The multi-graph fusion framework takes a closed-form solution for the optimization of S m (m = 1; :::; M), H and G. The time complexity of S m is O(Mn 2 ) and the time complexity of either H or G is O(n 2 ), where M and n, respectively, represent the number of the subjects and the number of brain regions. Hence, the time complexity of our multi-graph fusion method is O(lMn 2 ), i.e., linear to the subject size, where l is the iteration number and is less then 50 in our experiments. Moreover, the proposed multi-graph fusion method needs to store S m (m = 1; :::; M), H, and G in the memory with the space complexity O(Mn 2 ). The time complexity of '1-SVM is linear to the subject size, while its space time complexity is O(Mn 2 ) [Yang et al., 2015].

3 Experiments

We experimentally evaluated our proposed method, compared to four state-of-the-art methods, on two real neurodisease data sets with fMRI data in terms of binary classification performance.

3.1 Experimental Setting

Data Sets

The data set fronto-temporal dementia (FTD) contains 95 FTD subjects and 86 age-matched healthy control (HC) subject-

1. FTD was derived from the NIFD database managed by

the frontotemporal lobar degeneration neuroimaging initiative. The data set obsessive-compulsive disorder (OCD) has 20 HC subjects and 62 OCD subjects.

For all imaging data, we followed the automated anatomical labeling (AAL) template [Tzourio-Mazoyer et al., 2002] to construct the functional connectivity network for each subject with 90 nodes. The region-to-region correlation was measured by Pearson correlation coe�cient.

Comparison Methods

The comparison methods include the baseline method '1SVM embedded (L1SVM) in Liblinear toolbox [Fan et al., 2008], two popular methods in neuro-disease diagnosis, i.e., high-order functional connectivity (HOFC) [Zhang et al., 2017] and sparse connectivity pattern (SCP) [Eavani et al., 2015], and a deep learning method, i.e., simplify graph convolutional networks (SGC) [Wu et al., 2019].

L1SVM and SGC extract the upper triangle of the full FCN for each subject as the vector representation of the classifier. The methods (e.g., HOFC, SCP, and our proposed method) designed di�erent methods to transfer full FCNs to sparse FCNs, followed by extracting the vector representation. It is noteworthy that all methods can be directly applied for supervised learning and only two methods (e.g., SGC and our method) can be used for personalized classification.

Setting-up

In our experiments, we repeated the 10-fold cross-validation scheme 10 times for all methods to report the average results as the final results. In the model selection, we set �; � 2 f10 �3 ; 10 �2 ; :::; 10 3 g in Eq. (2), and fixed k = 10 since the value of k is insensitive to the result of Eq. (2). We further set C 2 f2 �10 ; 2 �9 ; :::; 2 10 g for '1-SVM. We followed the literature to set the parameters of the comparison methods so that they outputted the best results.

We designed 4 experiments to evaluate all methods, i.e., classification performance of supervised learning, classification performance of personalized classification, e�ectiveness of multi-graph fusion and e�ectiveness of brain region selection of our method. The evaluation metrics include ACCuracy (ACC), SENsitivity (SEN), SPEcificity (SPE), and Area Under the ROC Curve (AUC).

3.2 Result Analysis

Supervised Learning

In the experiments of supervised learning, we used all labeled subjects as the training set. We report the results of all methods in Figure 1 and list our observations as follows.

First, our proposed method achieved the best classification performance on two data sets, in terms of four evaluation metrics, followed by SGC, SCP, HOFC, and L1SVM. Specifically, our method on average improved by 2.17% and 1.71%, compared to the best comparison method SGC, respectively, on FTD and OCD, for all evaluation metrics. The possible reasons are that (i) our multi-graph fusion method takes the inter-subject variability, the heterogeneity across subjects, and the discriminative ability into account to output homogenous and discriminative representation, and (ii) our proposed method jointly selects features (i.e., brain regions) and conducts classification to avoid the influence of redundant features on high-dimensional data.

Second, L1SVM uses full FCNs to conduct classification such that outputting the worse classification performance. On the contrary, other methods use sparse FCNs. This indicates the reasonability of sparse FCNs, compared to full FCNs.

Third, the methods (e.g., HOFC, SCP and our method) design di�erent models to generate sparse FCNs, but our method achieved the best performance. This shows that our multi-graph fusion framework is feasible.

Personalized Classification

To verify the e�ectiveness of our proposed semi-supervised method, we randomly selected di�erent percentages of labeled subjects (i.e., 20%, 40%, 60%, and 80%) from the whole data set as the training set. In this case, the methods (i.e., L1SVM, HOFC, and SCP) only used labeled subjects to train the classifiers, while the methods (i.e., our method and SGC) used all subjects (i.e., labeled subjects and unlabeled subjects) to train the classifiers. We report the classification results of all methods in Figures 2 and 3.

First, our proposed method achieved the best performance, followed by SGC, HOFC, SCP and L1SVM. For example, our method on average improved by 2.31%, compared to the best comparison method SGC, in terms of accuracy, on two data sets with 80% labeled subjects for the training process.

Second, while the percentage of labeled subjects in the training set is small, all methods achieved worse performance. The main reason is that the lack of labeled subjects is di�cult to guarantee the performance of the classifiers.

Multi-graph Fusion E�ectiveness

The novelty of our method lies in the process of multi-graph fusion. In order to verify the fusion e�ect, we fed the vector representation outputted by our method to L1SVM and SGC. Note that, due to the space limitations, we only selected the best and the worst comparison methods. We reported the experimental results in Figure 4.

From Figure 4, we can see that the performance of methods (L1SVM and SGC) is better than the corresponding methods in Figure 1. This proves that sparse FCNs output by our proposed multi-graph fusion framework contains strongly discriminative ability.

Feature Selection E�ectiveness

In this section, we designed experiments to investigate the e�ectiveness of the selected features by our method. Specifically, our method selected 1270 and 898 nodes out of 4005 nodes, respectively, on FTD and OCD. We plot top selected brain regions of our method in Figure 5.

Based on the visualization of top selected brain regions, many selected regions from our method have been verified related to the neuro-diseases. Specifically, most of the nodes selected by our method occur in frontal and temporal lobes, which is consistent with the current neurobiological findings on FTD [de Haan et al., 2009]. In particular, our method finds the brain regions, such as orbital-frontal cortex, caudate, thalamus, which are included in the cortical-striato-thalamic circuits, and is considered as the theoretical neuroanatomical network of OCD [Gillan et al., 2015; Gillan et al., 2011].

4 Conclusion

In this paper, we proposed a new personalized disease diagnosis framework consisting of a multi-graph fusion method and a joint model for brain region selection and disease diagnosis. Compared with state-of-the-art methods, comprehensively experimental results on two real data sets verified the e�ectiveness of our proposed framework. In the future, we plan to conduct the brain functional connectivity analysis by considering the frequency with di�erent bands.

Acknowledgments

This work was partially supported by the Natural Science Foundation of China (Grants No: 61876046, 61836016, and 61672177); the Guangxi Collaborative Innovation Center of Multi-Source Information Integration and Intelligent Processing; the Guangxi "Bagui" Teams for Innovation and Research; the Marsden Fund of New Zealand (MAU1721); the Project of Guangxi Science and Technology (GuiKeAD17195062); and the Sichuan Science and Technology Program (No. 2019YFG0535).