Neighborhood structure-guided brain functional networks estimation for mild cognitive impairment identification

Pearson’s correlation

As mentioned earlier, PC is a relatively simple, efficient and popular method for establishing BFNs (Smith et al., 2013). We assume that each brain, due to some atlas, is segmented into n ROIs. And the edge weight matrix W = (w_ij) ∈ R^n×n of BFN is given as follows: (1)${\displaystyle w_{ij}=\frac{{\left(x_i-{\overline{x}}_i\right)}^T\left(x_j-{\overline{x}}_j\right)}{\sqrt{{\left(x_i-{\overline{x}}_i\right)}^T\left(x_i-{\overline{x}}_i\right)}\sqrt{{\left(x_j-{\overline{x}}_j\right)}^T\left(x_j-{\overline{x}}_j\right)}}}$ where w_ij represents the estimated functional connectivity between the ith and jth ROIs, x_i ∈ R^m×1 $\left(i=1,\dots ,n\right)$ represents the observed time series (i.e., BOLD signal) in relation to the ith ROIs, m represents the number of the time volumes in every series. All entries in ${\overline{x}}_i\in R^{m\times 1}$ are the average values of elements in x_i, $x_i-{\overline{x}}_i$ represents the centralization of signal x_i. Without loss of generality, we redefine $x_i=\left(x_i-{\overline{x}}_i\right)/\sqrt{{\left(x_i-{\overline{x}}_i\right)}^T\left(x_i-{\overline{x}}_i\right)}$. Consequently, Eq. (1) can be simplified as $w_{ij}=x_i^T{x}_j$, or an equivalent matrix form: (2)${\displaystyle W=X^{T}X}$ where X = [x₁, x₂, …, x_n] ∈ R^m×n represents the rs-fMRI data matrix, W = (w_ij) ∈ R^n×n represents the edge weight matrix of the estimated BFN.

Practically, Eq. (2) can serve as a solution to the following optimization problem (Li et al., 2017): (3)${\displaystyle {min}_W{∥W-X^{T}X∥}_F^2}$ where ||⋅||_F stands for the Frobenius norm of a matrix.

Sparse Pearson’s correlation

As we all know, because noises are usually included in BOLD signals, the original PC are more prone to bringing about a BFN with intensive links (Fornito, Zalesky & Bullmore, 2016). To alleviate this problem, a threshold structure is usually applied to sparse the PC BFNs in practical applications. However, the selection of thresholds is usually empirical and lacks elegant mathematical representations. Therefore, Li et al. introduced the l₁-norm as regularization term (Li et al., 2017), resulting in the PC_Sparsity model as below: (4)${\displaystyle {min}_W{∥W-X^{T}X∥}_F^2+\lambda \left|\right|W|{|}_1}$ where λ represents a regularized parameter to controll the sparsity of W. Apparently, the PC_Sparsity decreases to the original PC as λ = 0. Here, we select the l₁-norm to sparse the network since it is simple and popular.

Sparse representation

Unlike PC-based methodologies that measure the total correlation between ROIs, through degenerating the hybrid influence of other ROIs, partial correlations may cause more correlated BFN estimation, which depicts a relatively reliable interaction between brain regions (Peng et al., 2009). A universal method for calculating partial correlation is founded on the inverse covariance matrix estimation. Nevertheless, because of the singularity of the covariance matrix, the estimation of partial correlation may be maladaptive. For acquiring a steady estimate of BFN, an l₁-norm regularization term is usually presented in the partial correlation model. For instance, partial correlation estimation based on SR is a classic statistical method. The mathematical model of SR is represented as below:

${\displaystyle {min}_{w_{ij}}\sum _{i=1}^n\left(\parallel x_i-\sum _{j\ne i}{w}_{ij}{x}_j{\parallel }^2+\lambda \sum _{j\ne i}\left|w_{ij}\right|\right)}$ (5)${\displaystyle s.t.w_{ii}=0,\forall i=1,\dots ,n}$

Then, Eq. (5) can be further simplified to:

${\displaystyle \underset{W}{min}{\parallel X-XW\parallel }_F^2+\lambda {\parallel W\parallel }_1}$ (6)${\displaystyle s.t.w_{ii}=0,\forall i=1,\dots ,n}$ where ||⋅||_F and ||⋅||₁ represent the F-norm and l₁- norm of a matrix. ${∥X-XW∥}_F^2$ represents a adapting term, ${∥W∥}_1$ represents a l₁-regularized term used to get sparse solutions of W, λ represents the regularized parameter to control the equilibrium of two terms in Eq. (6). The constraint w_ii = 0 is applied for removing x_i from X to stay away from the problem of trivial solution.

Motivation

As described above, the classification performance of the brain disease recognition task critically depends on the quality of the constructed BFNs. Therefore, many BFNs estimation schemes have been proposed over the last few decades, and the three methods (i.e., PC, PC_Sparsity and SR) introduced in the previous section are the most common statistical methods. These three methods, whether using a hard-threshold to remove weak connection edges or introducing l₁-norm as a regularization term to generate an elegant mathematical representation to build the network, can ultimately obtain sparse BFNs, which is consistent with the scant prior information on brain network (Sporns, 2011). Although the above methods have been widely used in neurological diseases, the topological structure of brain networks is ignored in estimating BFNs. Besides, recent studies have indicated that due to the principle of simplicity in wiring costs, ROI closer in space has higher connectivity prospects than ROI farther in space (Fornito, Zalesky & Bullmore, 2016). Due to all the viewpoints above, we present a new BFN estimation method, which uses spatial distance information to reconstruct the PC model. The details of our method will be introduced in the following section.

Model and algorithm

In order to obtain sparse BFNs and incorporate the topology structure information, we propose a neighborhood structure-guided BFN construction scheme in this article. The whole framework is divided into the following steps.

Step1: Calculate the Euclidean distance between ROIs. There are many methods to measure distance, including Euclidean distance, Cosine similarity, Manhattan distance, and so on. Cosine similarity is more suitable for measuring the directionality and correlation between vectors, without considering their absolute value size. The Manhattan distance, on the other hand, is not as direct when considering the distances of various dimensions as the Euclidean distance, in some cases, it may not be possible to capture the feature relationships of the data. So, we choose the most suitable and commonly used Euclidean distance as the measurement strategy of spatial distance information between ROIs. The mathematical calculation formula for Euclidean distance is as follows:

${\displaystyle d_{ij}=\sqrt{{\left({\mathrm{x}}_{1,\mathrm{i}}-{\mathrm{x}}_{1,\mathrm{j}}\right)}^2+{\left(x_{2,i}-x_{2,j}\right)}^2+\cdots {\left(x_{m,i}-x_{m,j}\right)}^2}}$ (7)${\displaystyle =\sqrt{\sum _{l=1}^m{\left(x_{l,i}-x_{l,j}\right)}^2}}$ where x_i = (x_1,i, x_2,i, …, x_m,i)^T represents the mean time series extracted from the ith ROIs, d_ij represents the distance between the ith and jth ROIs. The distance matrix D = (d_ij) ∈ R^n×n is obtained by calculating the distance between all ROIs using Eq. (7).

Step2: Construct the adjacency matrix with topological structure. The distance between each ROIs and other ROIs is sorted from smallest to most significant according to the distance matrix D. Then, the KNN algorithm is used to find the k ROIs closest to the current ROIs (that is, the top k ROIs with the smallest distance) as neighbors and connect them with edges to generate an adjacency graph. The details are shown in Fig. 1. This graph can be represented by the following adjacency matrix S = (s_ij) ∈ R^n×n: (8)${\displaystyle s_{ij}=\left\{\begin{array}{cc}1,\hfill & e_{ij}\text{exists}\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}\right.}$ where e_ij is the edge between the ith and jth ROIs. Eq. (8) indicates that if an edge is connected between the ith and jth ROIs, the edge weight s_ij is regarded as 1. Otherwise, the edge weight s_ij is 0. In this way, the adjacency matrix S that only distinguishes between connected and unconnected networks is obtained. This is in line with the previous introduction that ROIs with close spatial distance are more inclined to be connected by edges. The selection of k value in KNN algorithm will be discussed in ‘Discussions’.

Step3: Estimating BFNs based on adjacency matrix S. According to matrix S, we propose a new sparse PC model for estimating BFNs. Specifically, the calculation formula of this model is as follows: (9)${\displaystyle w_{ij}=\left\{\begin{array}{cc}ll\frac{{\left(x_i-{\overline{x}}_i\right)}^T\left(x_j-{\overline{x}}_j\right)}{\sqrt{{\left(x_i-{\overline{x}}_i\right)}^T\left(x_i-{\overline{x}}_i\right)}\sqrt{{\left(x_j-{\overline{x}}_j\right)}^T\left(x_j-{\overline{x}}_j\right)}},\hfill & \text{if}{s}_{ij}=1\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}\right.}$ where w_ij represents the edge weight between the ith and jth ROIs. According to Eq. (9), the correlation coefficients between all ROIs with edges connected can be calculated as the edge weights w_ij. As a result, the sparse BFN containing topology can be finally generated. For facilitating the subsequent description, the BFN constructed by the KNN is abbreviated as K-BFN.

Data acquisition and preprocessing

For evaluating the performance of the presented approach, we conduct MCI recognition task in accordance with the benchmark database, i.e., Alzheimer’s Disease Neuroimaging Initiative (ADNI) dataset in the experiment. Concretely, 137 subjects (68 MCIs, 69 NCs) were participated as in the researches (Zhou et al., 2018; Van der Haar et al., 2023; Sun et al., 2021). All these data are publicly available in the ADNI database (http://adni.loni.usc.edu/). All subjects were scanned through a 3.0T Philips scanner, with the scanning time of 7 min, equivalent to 140 volumes. Table 1 shows the demographic information of these subjects.

For decreasing the impact of interference signals, the preprocessing pipeline in SPM8 toolbox (http://www.fil.ion.ucl.ac.uk/spm/) was applied to enhance the quality of fMRI data before BFN estimation. Especially, the first 3 volumes of one subject in the fMRI time course were removed to keep the signal stability. Next, the volumes leave behind were handled based on a popularly-used framework. The major stages were (1) modifying slice timing and head movement (Friston et al., 1996); (2) fitting the ventricle and white matter signals with Friston 24-parameters of head movement; (3) enrolling the modified rs-fMRI graphics to Montreal Neurological Institute (MNI) standard space (Brett et al., 2001); (4) performing spatial smoothing on the signal with a maximum half width of four mm and using a bandpass filter (0.015–0.15 Hz) on it. In the end, each brain was divided into 116 ROIs due to the Automated Anatomical Labeling (AAL) template (Tzourio-Mazoyer et al., 2002). The average time series of each ROI for each subject were extracted and placed in the data matrix X ∈ R^137×116 which is just used to estimate BFNs.

Experimental setting

Feature selection and classification

After gaining the BFNs, the follow-up mission is identifying MCI and NC due to the estimated BFNs. Now, the issue becomes the decision of which features and classifiers are going to be selected for MCI recognition. Given that feature selection and classifier design have significant effects on the identification performance, it’s ambiguous to tell whether the BFN estimation approaches or the follow-up steps in the follow-up pipeline is beneficial for the outcome. Hence, we merely take a relatively simple and efficient feature selection way (t-test with four accepted p-values) and a very popular used support vector machine (SVM) classifier (Chang & Lin, 2011), as our primary focus is BFN estimation. Figure 2 shows the basic procedure to classify MCI and NC using our designed pipeline. In the experiment, weights between ROIs are regarded as features for identification. As the AAL atlas contains 116 ROIs, there will be 116 × (116 − 1)/2 = 6670 connections/attributes for each BFN, which is far greater than the sample size (the number of subjects applied in the experiment), resulting in the collapse of dimensionality. For overcoming the issue, we apply two-sample t-test for feature selection in advance of MCI identifications.

Further, to make full use of the limited samples, i.e., the estimated BFNs of distinct subjects, during the training phase, we use 5-fold cross validation 100 times to investigate the effectiveness of the relevant approaches (Qiu et al., 2020; Diogo, Ferreira & Prata, 2022). Specifically, the 137 subjects of ADNI are first divided into five folds, each containing nearly the equivalent amount of connections the samples. Then, four folds are applied for training, and the other one is for testing. Note that the same linear SVM with default parameter C = 1 is applied in every method used in this work for fair. As shown in Fig. 2, an inner 5-fold CV is operated through grid-search to choose the optimal parameters. The accuracy is computed by averaging the outcomes of all subjects.

Network visualization

Using the preprocessed rs-fMRI data, we get the modified BFNs with the presented K-PC methods. Moreover, we construct BFNs for comparison according to the baselines including PC, PC_Sparsity and SR. Figure 3 shows the adjacency matrices of BFNs constructed by four various approaches for MCI recognition (the K = 40). From Fig. 3, we can observe that BFNs estimated based on PC_Sparsity, SR and K-PC methods are very clear compared with PC. The adjacency matrices of PC based BFNs are too divergent, and compared to them, the adjacency matrices of the other three methods are more clustered in a certain type of linear regression. And the BFN is not sparse, as pairwise complete correlation is applied for modelling the network adjacency matrix. Consequently, it results in new incorrect connections, affecting the final classification accuracy. Conversely, the networks on the basis of PC_Sparsity and SR are sparse because of the introduction of the l₁-norm regularizer. Besides, it is worth noting that the BFN estimated through K-PC exhibits a more pronounced module topological architecture than other methods. To evaluate the modularity of BFNs, the signed modularity maximization method is applied for calculating the modularity scores of BFNs constructed through various approaches (Gómez, Jensen & Arenas, 2009; Rubinov & Sporns, 2011). Moreover, for keeping off the randomness of the results, 28 subjects are selected equidistantly. As Fig. 4 shown, the BFN estimated by K-PC is significantly more modular than the other three baseline methods, indicating that K-PC tends to achieve cleaner networks. Additionally, modularization of BFN estimated by PC_Sparsity, SR and K-PC acquire higher scores than that by PC in most cases, suggesting that moderate sparsity promotes the modularity architecture of the network.

Classification performance

An array of quantitative mensuration is conducted to test the performance of different approaches, like classification accuracy (ACC), sensitivity (SEN), specificity (SPE), positive predictive value (PPV) and negative predictive value (NPV). The definitions of each method are shown in Table 2, where TP, TN, FP, and FN represent true positive, true negative, false positive and false negative.

Table 3 reports the MCI recognition outcomes of different approaches under five evaluation metrics. As shown in Table 3, our presented K-PC exceeds the baselines according to the average ACC score under the circumstance that p-values are 0.01, 0.05 and 0.001 in accordance with 5-fold CV. Note that our K-PC gets the ACC of 81.75%, SEN of 83.09%, SPE of 80.43%, PPV of 80.70% and NPV of 82.87% when p-value is 0.05, which are significantly better than the baselines. Regarding the situation where the p-value is 0.005, our K-PC, i.e., even if the ACC, SEN and NPV performances are not satisfactory, also gets prominently better results than the other compared approaches on SPE and PPV. The outcome of the experiments evidence that the presented K-PC framework can enhance the baselines, i.e., at least the discrimination of the estimated BFNs.

For verifying the validity of the presented approach, we compared it with the most commonly used time series methods in deep learning. We chose classic BiGRU, BiLSTM, and Transformer for comparison. From Table 4, it can be seen that among these methods, Transformer performs the best in terms of five performance indicators. The time series methods such as Transformer typically directly model raw data without an accurate feature selection process. They rely on the learning ability of the model to learn the most relevant features from the raw data. When selecting the appropriate p-value, the K-PC method outperforms the traditional time series prediction methods. The common algorithms for deep learning exhibit unstable stability, with significant performance differences between each fold, indicating poor generalization ability and excellent performance in capturing complex relationships in MCI data. Among the five evaluation metrics (i.e., ACC, SEN, SPE, PPV and NPV), the degree of variation in the K-PC method was 1 / 7 to 1 / 5 that of other traditional time-series approaches. In the field of MCI identification outcome, we can consider the K-PC method superior to deep learning traditional temporal methods.

Parameter analysis

It is well known that feature selection is a widely used dimensionality reduction technique, which can effectively reduce redundant features in data. Therefore, it is also generally applied in the research of neurological diseases as per fMRI data. Now, there are many proposed feature selection approaches. In this article, we only focus on two-sample t-test, since it is relatively simple and efficient in practice. Specifically, we conducted comparative analysis based on the t-test with four various p-values, i.e., 0.01, 0.05, 0.001, 0.005, k = 40, in the experiment, and then show the classification results in Fig. 5. Compared with the outcomes, we can see the following facts: when using t-test under the p = 0.05, our method can obtain the best classification accuracy, of 86.13%. However, we note that compared with other methods, the SR-based method can achieve a better final accuracy of 93.43% with p-values of 0.01, and even the classification performance is much higher than our K-PC. The phenomenon suggests that the selection of the P value is crucial to the classification outcomes, so in the following discussion, we choose p = 0.05 as the parameter of t-test for feature selection. In conclusion, appropriate p-value selection has a positive impact on achieving better identification outcomes.

Sensitivity to network modelling parameters

Virtually, in addition to the p-values applied in the feature selection stage, the network modeling parameters usually have an essential impact on the structure of BFN and the subsequent classification accuracy. For studying the sensitivity of distinct approaches to relevant parameters, we calculate the evaluation metric with various parameters (we set p = 0.05 in this part for better performance according to the result from 6.1).

As the experimental outcomes shown in Fig. 6, the variation in classification accuracy corresponds to that in parametric values, suggesting that the majority of the BFN estimation approaches are easily affected by the free parameters. In terms of PC, it can be found that moderate sparsity inclines to lead to higher classification accuracies. Interestingly, the dense BFNs without thresholding operation (parameter value = 0) can achieve high accuracy as well, which may be due to the use of a t-test to remove weak discriminative connections before classification. Compared with PC, PC_Sparsity can obtain relatively high recognition results of brain diseases in the parameter interval [2⁻⁵, 2⁻¹], indicating that the sparsified BFNs can effectively remove the redundant connections in PC, thus improving the recognition accuracy of patients with brain diseases. However, when the parameter is 2⁰, the accuracy of PC_Sparsity is reduced to 50.36%, and the accuracy of PC_Sparsity does not change when the parameter is more than 2⁰, which proves that being too sparse can actually lead to a decrease in classification accuracy. For SR, l₁-regularizer in the BFN estimation model is applied to remove the weak connections and control the estimated BFN density. From Fig. 6C, it can be found that the best value of λ is 2³ for SR, and the corresponding accuracy is 85.43%. In particular, according to Fig. 6D, we can find that although the classification performance of K-PC fluctuates with the number of neighbors k, it is least affected by the parameters compared with the three baseline methods. More importantly, K-PC obtains no less than 78% results in the entire parameter range [110, 100, …, 20, 10], which indicates that its recognition performance is more advantageous relative to other methods.

Top discriminative features

Beyond the evaluation metrics, it is more enlightening to find which features/connections in BFN will result in better recognition performance. For BFNs based on K-PC, we study some favorable clues to recognize subjects with MCI from NCs. Particularly, we choose 80 connections whose p-values are lower than 0.005. As can be seen from Fig. 7, the gauge of the arc is in reverse proportion to the p-value, suggesting the discriminative ability of the corresponding edge. In terms of MCI identification, we sort the brain regions according to the chosen discriminative features and identify the first three, hippocampus, left cerebellum 6, and amygdala. They are generally considered to have biological significance for MCI identification, according to the researches (Greicius, 2008; Albert et al., 2011) for the last few years, which further validates the effectiveness of our K-PC method in constructing reliable functional brain networks for MCI identification.

However, this method also has its limitations because the p-value indicates the degree of influence of this feature. If the p-value is smaller, the significance of the feature’s influence is higher. Therefore, features with a p-value less than 0.005 are selected. However, the smaller the threshold value, the fewer features will be selected, which may impose certain restrictions on forming a more accurate, clearer, and continuous boundary highlighting MCI features.

From the perspective of algorithm optimization, deep neural network (DNN) offers notable advantages over KNN in terms of feature learning, handling high-dimensional data, generalization ability, and computational efficiency. However, in the context of this study involving small-scale datasets and interpretability requirements, KNN remains a viable option. Nevertheless, further optimization is necessary to strike a balance in defining feature boundaries.