A total of 143 participants from the ADNI dataset (Wang et al., 2019b) were used in this work, including 95 MCIs and 48 NCs. Each participant was scanned at one or more visits in this study, and the interval between two visits is at least 6 months, resulting in a total of 299 resting-state fMRI scans. These 299 scans include 154 NC cases and 145 MCI cases. The scanning parameters were listed as follows. The slice thickness is 3.31 mm, TE (echo time) is 30 ms, TR (repetition time) is 2.2 − 3.1 s, the in-plane image resolution is 2.29 − 3.31 mm, and the scanning time for each subject is 7 min (resulting in 140 volumes). Table 1 shows the demographic information of these 299 scans.

To reduce the influence of nuisance signals, the preprocessing pipeline in FSL FEAT is used in this work. Specifically, for each subject, we discard the first 3 volumes in the fMRI time course for signal stabilization. And then the remaining 137 volumes were processed via the standard pipeline. The main steps in preprocessing include slice timing correction, head motion estimation, bandpass filtering, and regression of nuisance covariates (i.e., white matter, cerebrospinal fluid, and motion parameters) by mean regression. Note that subjects were removed if their head motion > 2.0 mm of maximal translation or 2.0^o of maximal rotation. After that, we performed skull stripping based on T₁-weighted MRI and aligned the skull-stripped fMRIs onto the Montreal Neurological Institute space. Then, we used a Gaussian kernel with full-width-at-half-maximum of 6 mm to spatially smooth the volumes. The subjects with more than 2.5 min of frame-wise displacement (FD > 0.5) were excluded from further analysis. Finally, every brain was divided into 116 ROIs based on the Automated Anatomical Labeling (AAL) template (Tzourio-Mazoyer et al., 2002), and the mean time series (with band-pass filtered 0.015 − 0.15 Hz) of each ROIs were extracted as the input data for FBN estimation.

Due to the crucial role in exploring the neurodegenerative diseases, many FBN estimation methods have been proposed in the past decades. In this section, we briefly review several PC-based methods that are closely related to our study.

As pointed out earlier, PC is the simplest and most popular method for FBN estimation (Smith et al., 2013). The edge weight of PC-based method for FBN estimation is defined as follows:

Where xi∈Rm is the blood oxygen level-dependent (BOLD) signals associated with the i^th ROI, and xi¯∈Rm is a vector containing the mean of the elements in x_i.

Without loss of generality, we redefine xi=(xi−x¯i)/(xi−x¯i)T(xi−x¯i). Then, Equation (1) can be simplified as wij=xiTxj, which is exactly the solution of the following optimization problem:

or its matrix form shown as follows:

Where W is the adjacency matrix to be estimate by PC, X=[x1,x2,⋯,xn]∈Rm×n is the data matrix containing the fMRI time courses, and n is the number of ROIs. We will note shortly that such an optimization view of PC can help us to develop new and more flexible FBN estimation methods.

Despite its popularity, PC aims to measure the full correlation between signals of all ROIs in the brain, thus generating dense networks where all nodes/ROIs are densely connected. This is not consistent with the empirical finding that sparsity has been proven to be the most popular property of FBN (Sporns, 2016). Therefore, in practice, a threshold is generally used to sparsify the originally estimated FBN by removing the edges with weak weights. An alternative to sparsify FBN is the l₁-regularized PC (PC_Sparsity) (Li et al., 2017) whose model is given as follows:

where ||W||₁ is the l₁-norm of a matrix for encoding the sparsity prior of FBN, and λ is a regularized parameter for controlling the sparsity of W.

Besides sparsity, FBNs usually have more richer structures (Sporns, 2016) such as modularity, and some of these structures may guide us to estimate more reasonable FBNs. Recently, Zhou et al. (2018b) developed an M-FBN method to estimate FBN by further introducing a trace norm regularizer into Equation (4), resulting in the following model:

Where ||W||_* is the trace norm of the matrix W, and λ₁ and λ₂ are regularized parameters for controlling the balance between the three terms in the objective function. With the combination of l₁ and trace norms, Equation (5) has been verified to be able to discover modular structures of the estimated FBN. However, since this model is an approximation of the sparse and low-rank matrix, the network estimated by M-FBN could not include a notable modular structure.

Several previous Studies (Grone et al., 1990; Zhou et al., 2018b) proposed to jointly minimize the l₀-“norm” and the rank of the edge weighting matrix W, since a sparse (via minimizing l₀-“norm”) and low-rank (via minimizing the rank) matrix tends to result in modular structures. However, these two regularizers are both non-convex with respect to W, making the optimization problem intractable. In practice, they are generally relaxed to l₁-norm ||W||₁and trace norm ||W||_* respectively. Such relaxation only achieves an approximation of a sparse low-rank matrix, and thus fails to guarantee the modular structure of FBN.

Motivated by the theorem (Heider, 1946) that the multiplicity k of the eigenvalues 0 of a Laplacian matrix is equal to the number of connected components in a graph, we propose a new method to estimate FBNs with an accurate number of modularity. Denote n (n = 116 in this work) as the number of ROIs and k (k = 8 in this work) as the number of modules. Specifically, based on the l₁-regularized PC method defined in Equation (4), we propose to constrain the rank of L_W as n−k to generate a FBN with an accurate number of modules. Mathematically, the proposed AM-PC model is given as follows:

Due to the nonnegative constraint in Equation (6), the l₁-norm can be replaced by the sum of the elements in the matrix. In fact, the nonnegative assumption for edge weights is supported by the structural equilibrium theory (Fan, 1949; Cartwright and Harary, 1956), and can simplify the subsequent analysis for FBNs. With Equation (6), we can explicitly construct an FBN with a total of k modules for each subject.

Given the data matrix X and a matrix A = X^TX, we can rewritten Equation (6) as follows:

In general, it is not straightforward to solve (Equation 7) since rank(L_W) = n−k is a strict constraint. In what follows, we derive an efficient optimization algorithm to solve this challenging problem.

Denote σ_i(L_W) as the i-th smallest eigenvalue of L_W. Accordingly, Equation (7) can be equivalently converted to the following problem:

Where L_W is positive semi-definite to guarantee σ_i(L_W)≥0, and a large α enables Σi=1kσi(LW)=0 to meet the constraint of rank(L_W) = n−k.

Furthermore, our proposed optimization method is based on the Ky Fan's theorem (Grant and Boyd, 2014) as follows:

By combining Equations (8), (9), we have the problem:

Which can be efficiently solved by the following alternating optimization algorithm.

Step 1: When W is fixed, Equation (10) becomes

Whose the optimal solution is formed by the k eigenvectors corresponding to the k smallest eigenvalues of L_W.

Step 2: When F is fixed, Equation (10) becomes

which is equivalent to the following problem:

For simplicity, we denote ||fi-fj||22 as f_ij. Then, by expanding and combining the like terms, Equation (13) can be rewritten as follows:

Which is a quadratic programming problem and can be easily solved by, for example, the CVX toolbox (Grant and Boyd, 2014). We summarize the optimization algorithm to solve Equation (10) in Algorithm 1.

Based on the learned edge weighting matrix W, we can construct a specific FBN with k modules for each subject. Given FBNs for all subjects, the subsequent task is to classify MCI and NC based on the estimated FBNs. Since the given FBNs are symmetry matrices, the upper and lower triangle matrix features are the same, and for ease of calculation, we take the upper triangle features including (116 × (116 − 1))÷2 = 6, 670 elements. Although the edge weight matrix contains the full information of the network, it typically causes the curse of dimensionality, since the number of feature dimension (Kuo and Sloan, 2005), i.e. 6670, is far greater than the number of subjects.

To address this issue, we propose to employ the t-test algorithm to select the most informative features by a fixed p values (p = 0.05 in this work). As pointed out in Wee et al. (2014), the classifier design has a big influence on the ultimate accuracy. For this reason, we employ a linear support vector (SVM) with the default parameter (i.e., C = 1) as classifier (Chang and Lin, 2011) in this work, considering that it is simple and widely used in neuroimaging-based brain disorder classification.

In the experiments, we use a 5-fold subject-level cross-validation (CV) strategy (Li et al., 2017) to evaluate different methods, to ensure that fMRI scans of the same subject will not appear in both training and testing sets. To reduce the biased introduced by random partition cross validation, the 5-fold CV process was repeated 100 times for all methods. Besides, for the fair comparison, all the competing methods use the same linear SVM (with C = 1) as the classifier.

Since the parameters may significantly affect the structure of the constructed networks and the classification results, we select optimal parameters through grid search via inner 5-fold CV based on only training data. Specifically, we uniformly utilize 11 candidate values [2⁻⁵, 2⁻⁴, ⋯ , 2⁵] for the regularization parameters (λ, λ₁ and λ₂) in the four competing methods (i.e., PC_Sparsity, PC_Sparsity+, M-PC, and M-PC+). The proposed AM-PC has three parameters, i.e., k, λ and α. The module number k is empirically set as 8 based on the prior knowledge (Wong, 2015). The optimal value of α can be determined by a heuristic approach. That is, we first initialize α with a small value. Then, in each iteration, we compute the number of zero eigenvalues in L_W, if it is larger than k, then divide α by 2; if it is smaller than k, then multiply α by 2; otherwise stop the iteration. Thus, in AM-PC, we only need to tue the parameter λ. In our experiment, the elements in the network tend to zero when λ is equal to the number in [2⁰, 2¹, ⋯ , 2⁵], and the network loses its discriminative capability. Therefore, the optimal value of λ in our AM-PC is selected from [2⁻¹¹, 2⁻¹⁰, ⋯ , 2⁻¹] via inner 5-fold CV.

In this section, we perform MCI identification (i.e., MCI vs. NC classification) based on the FBNs estimated by seven different methods (including PC, PC+, PC_Sparsity, PC_Sparsity+, M-PC, M-PC+, AM-PC). For seven methods, we employ the same t-test algorithm for feature selection and the linear SVM for classification to ensure the fair comparison. Four metrics are used to evaluate the classification performance, including accuracy (ACC), sensitivity (SEN), specificity (SPE) and AUC. Denote TP, TN, FP and FN as true positive, true negative, false positive and false negative, respectively. These four metrics are defied as follows: ACC=TP+TNTP+FP+TN+FN, SEN=TPTP+FN, SPE=TNTN+FP, and AUC is the area under the ROC (receiver operating characteristic) curve.

In Table 2, we report the classification results of MCI vs. NC classification achieved by seven different ways. The term marked by “*” denotes that the result of proposed AM-PC is significantly better than that of all six competing methods (with p < 0.05). It can be seen from this table that our proposed AM-PC method consistently outperforms the six competing methods in terms of four evaluation metrics. And the AM-PC is significantly better (with p < 0.05) than the six methods in terms of ACC, SPE and AUC values. These results imply that the modularity prior introduce by our method plays an important role in modeling the rich topological structure of functional brain networks, thus helping boost the identification performance of subjects with MCI.

We visually compare the FBNs estimated by our AM-PC and six competing methods. In Figures 3A–G, we take a randomly selected subject from ADNI as an example to visualize the FBNs estimated by these seven methods. Here, the corresponding parametric values are λ = 2⁻² for PC_Sparsity and PC_Sparsity+, λ1=23, λ2=22 for M-PC and M-PC+, and λ = 2⁻¹ for our AM-PC.

It can be seen from Figures 3A–G that the FBNs estimated by PC_Sparsity, PC_Sparsity+, M-PC, M-PC+ and AM-PC are sparse. This is due to the introduction of the l₁-norm regularizer in these three methods. Besides, the FBN estimated by our proposed AM-PC method shows clearer modular structure than others. The underlying reason could be that AM-PC apply a low-rank constraint to L_W, thus yielding more clear modules in the estimated functional brain network.

To show the modular structure of the network constructed by AM-PC more clearly, we use scatter plots to illustrate the networks estimated by M-PC and our AM-PC in Figures 3H, I, respectively. From Figures 3H, I, we can observe that our AM-PC can generate an accurate number (i.e., k = 8) modules in the estimated FBN, compared with the M-PC method that aims to generate FBNs with approximate modules. These results further validate the effectiveness of the proposed method in generating FBNs with clear modular structures.

As a complex network, FBN include the hierarchical structure (Meunier et al., 2009, 2010), and the number of modules increase with the deeper of the FBN level. Thus, it is a practical problem to select the optimal values of the parameter k in the proposed objective function in Equation (6). Previous studies works (Bertolero et al., 2015; Geerligs et al., 2015; Gallen et al., 2016; Murakami et al., 2018) have shown that the number of modules in the FBN is usually less than 10. To study the influence of the number of modules on the classification performance, we report the four performance metrics (i.e., ACC, SEN, SPE and AUC) achieved by the proposed AM-PC method using different values of k in Figure 4. From Figure 4, we can see that the AM-PC method can achieve overall stable results when 6 ≤ k ≤ 9, and the best AUC value is achieved with k = 8.

We further show the most discriminative functional connections in FBNs estimated by the proposed AM-PC method. We empirically set k = 8 and use the t-test (with p = 10⁻⁶) to select the top 12 most informative connections. The identified discriminative connections are shown in Figure 5, where the discriminating power of a connection between two ROIs is represented by the thickness of an arc. As can be seen from Figure 5, the discriminative brain regions corresponding to these selected connections include several important ROIs, such as right hippocampus, right amygdala, and middle temporal gyrus. Especially, a clear discriminative functional connectivity exists between right temporal pole sup and right amygdala, and such connectivity plays an important role in cognition and emotion (Menon, 2018). The discriminative ROIs identified by our method also include right hippocampus, a brain region that is primarily associated with memory (Disouky et al., 2022). This finding is consistent with previous research (Albert et al., 2011; Kesler, 2014; Zhu et al., 2015), which further validates the effectiveness of our AM-PC method in constructing reliable functional brain networks for MCI identification.

In Figure 6, we visually show the modular structure of the average functional brain network (i.e., the average FBN of all training subjects), including 8 modules (with α = 2⁻³). The gray dotted lines represent the 12 most discriminative connections identified by our AM-PC method (see Figure 6). In the following Figure 7, we map Figure 6 on the brain template through BrainNet Viewer (Xia et al., 2013), where each node is a brain region, and each edge is the link between the brain region.

From Figures 6, 7, we have the following observations. First, we can clearly see that these eight modules are sparsely connected with each other, while nodes/ROIs within each module are densely connected, which caters to the results in previous papers (Meunier et al., 2010; Bertolero et al., 2018). Second, the ROIs contained in the Module 1 and Module 3 are associated with cognitive functions of the brain, involving the middle temporal gyrus (ROI IDs: 85 and 86), hippocampus (IDs: 37 and 38), parahippocampus (IDs: 39 and 40), precuneus (IDs: 67 and 68), amygdalae (IDs: 41 and 42), supramarginal gyrus (IDs: 63 and 64), inferior parietal lobules (IDs: 61 and 62), superior-medial frontal gyrus (IDs: 23 and 24), and medial orbitofrontal gyrus (IDs: 25 and 26). These regions are believed to be biologically associated with MCI, as reported in previous studies (Yetkin et al., 2006; He et al., 2007; Fair et al., 2008). Besides, we found that most of the 12 most discriminative connections are distributed between Module 1, Module 2 and Module 3. These results may imply that the interruption of the connections between the three modules could be used as potential biomarkers for MCI detection.

We also compare our AM-PC method with several state-of-the-art methods on FBN-based MCI identification with rs-fMRI data from ADNI. In Table 3, we briefly summarize the classification results of several previous studies as well as our method, where the top 2 best results are shown in bold. As can be seen from Table 3, our AM-PC method can achieve the overall comparable results compared with four SOTA methods. Even though the ACC and AUC results reported in Yang et al. (2019) are better than ours, their evaluation was based on a relatively smaller dataset compared with this work.

The current work has some limitations. First, in this work, we obtained the modularity of FBNs based on the PC-based networks. In fact, our method can also be combined with other methods (such as sparse representation) to estimate FBNs with modularity (Qiao et al., 2018), which will be our future work. Second, we focus on the within-module connection in this work by introducing an accurate number of modules in the estimated FBNs, without emphasizing between-module connections. In the future, we plan to incorporate both within-module and between-module connections in to the proposed framework for FBN-based brain disease analysis. Third, the small sample size could be an important limitation to the generalizability and replicability of this study. To alleviate this problem, we will utilize transfer learning (Pan and Yang, 2010; Valverde et al., 2021) or meta-learning (Finn et al., 2017; Hospedales et al., 2021) strategy to model the modularity structure of FBN.