Suppose x_i is the multivariate random variable (random vector) associated with the i^th ROI. Then, the observed rs-fMRI signals¹ x_i = (x_1i, x_2i, ⋯, x_ni)^T, i = 1, 2, ⋯, p can be considered as a sampling of the multivariate random variable (or population) x_i, where p is the number of ROIs and n is the number of time points. Since our goal is to estimate the edge weights of the BFN, the simplest and empirically effective way is to calculate the sample PC coefficient w^ij(i,j=1,2,⋯,p) of pair-wise ROIs, as follows.

where x¯i is the mean vector corresponding to x_i. Under Gaussian assumption, Eq. 1 gives an asymptotically unbiased estimation for the population PC coefficient. Without loss of generality, we redefine xi≜(xi-x¯i)/(xi-x¯i)T(xi-x¯i). Then, the estimator of population PC coefficient can be simplified as follows:

where X = [x₁, x₂, ⋯, x_p] is the rs-fMRI data matrix whose columns are the rs-fMRI time series associated with different ROIs. Therein, W^ is the generalized estimator in matrix form.

Sparse representation is one of the commonly used methods for calculating partial correlation among ROIs. A regularization term encoding sparsity prior is introduced into the BFN

estimation model. Specifically, the mathematical model of SR is given by:

where W is the edge weight matrix of BFN.

Similar to PC, we can rewrite SR in matrix form:

where the constraint w_ii = 0 plays a role in removing x_i from X to avoid the trivial solution.

Despite its popularity and effectiveness, the traditional PC can only construct the low-order BFN. That is, the connection between two ROIs is determined by the correlation of the corresponding rs-fMRI time series. However, in practice, a connection can be described in both low-order and high-order views. For example, we can directly define a connection between ROI i and j if there is a relationship. Besides, if ROIs i and j are connected to the same brain region, we can infer with a great possibility that there is a connection between ROI i and j. The former corresponds to the correlation in the traditional low-order view, while the latter can be considered as CC in a high-order perspective. This results can be achieved by a two-step procedure. First, the low-order BFN is estimated via PC. According to the formula in Eqs 1 or 2, the adjacency matrix W^=(w^ij)p=p of the PC-based BFN can be calculated as follows.

where w^i and w^j are the i^th and j^th columns of W^, respectively. For simplicity, in Eq. 5, w^i and w^j has been centralized and normalized as the case in Eq. 2. As a result, the CC-based high-order BFN is defined as follows,

As we know, PC is a measure of the linear correlation between pair-wise rs-fMRI time series associated with the ROIs. In other words, a time series x_i = (x_1i, x_2i, ⋯, x_ni)^T can be linearly represented by other time series x_j = (x_1j, x_2j, ⋯, x_nj)^T as

where a_ij is the representative coefficient, and ε_i = (ε_1i, ε_2i, ⋯, ε_ni)^T is the random error vector. That is, for each variable, we have:

Generally, we assume that the random variable ε_ki follows a normal distribution with mathematical expectation 0, i.e., ε_ki∼𝒩(0, σ²). Therefore, given x_kj and a_ij are constants, x_ki follows the normal distribution x_ki∼𝒩(a_ijx_kj, σ²). Then, The following formula can be obtained by the maximum likelihood estimation of x_ki (see Appendix A for details).

Note that Eq. 9 can be equivalently written as the following least-squares problem:

The optimal solution to Problem (10) is given by a^ij=(xjTxj)-1xjTxi=xiTxj=w^ij, considering that all of time series x_i, i = 1, 2, ⋯, p have been normalized by xi||xi||. This means that the solution of Eqs 9, 10 is the same as PC shown in Eq. 2. Therefore, in what follows, we only use w_ij instead of a_ij for the consistency of mathematical notations.

To provide a more flexible framework for BFN estimation, we further generalize PC in Bayesian view by introducing a prior distribution on w_ij. Although various distributions can be used as the prior, here we first consider the standard normal distribution, i.e., w_ij∼𝒩(0,1), since it provides a basis for understanding more complex cases. However, in practice, the entries w_ij in W may not be apriori independent of each other, but exist a relationship. Due to w_ij∼𝒩(0,1), the edge weight w_ij of BFN has the following prior probabilistic density:

Next, a maximal posterior estimation of w_ij (see Appendix B for details) can be obtained as follows,

The above problem is equivalent to the regularized least-squares problem:

where λ = σ² in the case of standard normal distribution. In practice, λ≜σ2/σ02 is a hyper-parameter that controls the balance between the two terms in Eq. 13, where σ02 corresponds to the variance of normal distribution of the edge weight w_ij. Setting the gradient of the objective function to zero, we obtain the optimal solution of Eq. 13 as follows:

We find that it is a shrinkage of the original estimation of PC which helps remove the weak connections in BFN.

Note that Eq. 13 only considers finding one edge weight at a time. Without loss of generality, with the assumption that the variables w_ij in W are independent, we can rewrite Eq. 13 in the following matrix form:

where tr(WWT)=∑i,jwij2 is the trace operator of WW^T. As a result, we achieve the estimation of BFN in a batching way as follows,

This formula is essentially a generalization of Eq. 2 with a shrinkage factor (1+λ)⁻¹. When λ+0, Eq. 16 reduces to the traditional PC.

In Section “Pearson’s Correlation-Based Brain Functional Network Learning Framework in Bayesian View,” we reformulate PC and then generalize it in Bayesian view by introducing a standard normal prior w_ij∼𝒩(0,1) for each pair of ROIs (i, j), i, j = 1, 2, ⋯, p. However, in practice, the entries w_ij in W may not be apriori independent of each other, but exist a relationship. Even so, Section “Pearson’s Correlation-Based Brain Functional Network Learning Framework in Bayesian View” provides a flexible probabilistic framework to develop new brain network estimation methods. Inspired by this point, we lay down theoretical support for CC from the Bayesian perspective. More importantly, instead of assuming that w_ij in W are independent, we propose a Bayesian high-order model (BHM) for BFN estimation by introducing the prior of matrix-variate normal distribution to the low-order BFN W. BHM learns the high-order relationship from the data automatically, rather than manually define as the case in CC. Interestingly, such a scheme can simultaneously learn low-order and high-order BFNs by considering the spatial structure of network connections. To distinguish the low-order and high-order correlations, an illustration is shown in Figure 1. The connections among w_ij can be considered as a high-order correlation h_ij,kl, while w_ij denotes the traditional low-order linear correlation between ROIs.

To evaluate the effectiveness of the proposed BHM, we conduct experiments on Autism Brain Imaging Data Exchange (ABIDE) database. The objective is to identify subjects with ASD from typical controls (TCs). Considering the heterogeneity of multi-site data, we only use data from the NYU site in our study. The dataset includes 184 subjects (79 ASD patients and 105 TCs). The detailed scan procedures and protocols are described on the ABIDE website.² The demographic information of all participants is summarized and displayed in Table 1.

All rs-fMRI images were acquired using a standard echo-planar imaging sequence on a clinical routine 3T Siemens Allegra scanner. During the 6-min rs-fMRI scanning procedure, most subjects were required to relax with their eyes focusing on a white fixation cross in the middle of the black background screen projected on a screen. A few participants close their eyes. The functional scan parameters include the flip angle = 90°, 33 slices, TR/TE = 2000/15 ms with 180 volumes, FOV = 240 mm and voxel size = 3 ×3 ×4 mm³. The rs-fMRI data were preprocessed by DPARSF³ software. Specifically, to avoid the interference of early signal instability, the first 5 rs-fMRI volumes of each subject were discarded. The remaining volumes were calibrated as follows: (1) Slice timing correction and head motion correction; (2) Regression of nuisance signals (ventricle, white matter) and head-motion with Friston 24-parameter model (Friston et al., 1996); (3) Normalization and register to MNI space with resolution of 3 3 3 mm³; (4) Segmentation using DATTEL; (5) Spatial smoothing by a kernel of 6 mm. After that, since our focus is functional connectivity, the rs-fMRI time series signals were partitioned into 160 ROIs, according to the functional atlas Dosenbach 160 (Dosenbach et al., 2010). Finally, the mean time series of the ROI were put into a data matrix X ∈ R^175=160, which will be used for the subsequent BFN estimation.

With the preprocessed rs-fMRI data, we estimate the low- and high-order BFNs using the proposed method, i.e., Bayesian low-order Network W and Bayesian high-order Network Ω, respectively. For comparison, we also choose PC, SR, and traditional CC as baseline methods to construct BFNs.

Once the BFNs are constructed, the next step is feature selection and classification. In our study, we directly use edge weights of the estimated BFN as features for ASD identification. Despite its simplicity (without complex feature design), such a scheme easily causes the curse of dimensionality due to limited sample size. As described previously, the number of ROIs is 160 and thus the estimated feature edges are 160=(160−1)/2=12720, which is far greater than the sample size (i.e., the number of subjects 184). To alleviate the problem of small sample size, we adopt a two-sample t-test with an empirically fixed p values to select features before ASD classification. In our experiments, we evaluate five candidate parametric values of p, that is [0.001,0.005,0.01,0.05,0.1]. The specific parameter analysis results are given in Section “Sensitivity to Network Modeling Parameters.”

To perform the following classification task, we use a linear support vector machine (SVM) (Chang and Lin, 2011) with default C = 1 as the classifier. To evaluate the model, we adopt leave-one-out cross-validation (LOOCV) in our experiments due to the limited data samples. Specifically, a LOOCV works in each run and only one sample is used to test while the rest are used to train a classifier. The final performance is obtained by the averaged results of all the runs. Note that the model parameters are involved in certain methods, including SR and BHM. Therefore, we additionally adopt an inner LOOCV procedure on the training data to obtain the optimal parametric value. Specifically, for SR, the regularization parameter λ is set to [2⁻², 2⁻¹, 2⁰, 2¹, 2²]. For the proposed BHM, the regularization parameter λ is set to [0.0001, 0.001, 0.01, 0.1, 1]. To be consistent with the number of parameters in other methods, the coefficient δ of the perturbation involved in Eq. 20 is set to 0.1 empirically.

To evaluate the classification results of different methods, we use accuracy (ACC), sensitivity (SEN), specificity (SPE) as performance metrics. The definition of these quantities are reported in Table 2. Note that, in this work, we treat ASD patients as the positive class while the NCs as the negative class.

In Table 3, we report the ASD classification results of five methods. It can be observed that the Bayesian low-order network (BHM-W) and Bayesian high-order network (BHM-Ω) constructed by the proposed BHM perform better than the BFNs constructed by the traditional PC and CC, respectively. Moreover, BHM-Ω achieves the best performance. Besides, the high-order BFNs (traditional CC and BHM-Ω) are associated with better recognition performance when they are compared with the baseline methods PC, SR, and BHM-W. This means that the high-order network structure can provide more helpful information for BFN analysis to some extent. Furthermore, for two corresponding low-order methods, the performance of the traditional PC and BHM-W are approximately similar whereas BHM-W has slightly better accuracy than the traditional PC. This may benefit from the guidance information provided by the Bayesian high-order network Ω in the optimization process.

To evaluate the BFNs estimated by different methods, we randomly select a subject and visualize the BFNs constructed by PC, SR, CC, BHM, as shown in Figure 2. Specifically, the different colors of Figure 2 indicate different weights of the edge weights matrix (i.e., the BFN), ranging from −1 to 1. It is observed that: (1) Compared with SR, the BFNs estimated by the correlation-based methods (i.e., PC, CC, BHM) are denser since the sparsity prior is introduced into SR. (2) There are fewer areas of cold colors in the PC network heatmap, implying that the edges with the negative weights are less. (3) Compared with PC, CC’s network heatmap shows a sharper distinction between the areas of warm and cold colors, indicating that the positive edge weights of CC’s BFN are larger and the negative edge weights are smaller. (4) The BHM-W estimated from the Bayesian perspective has a greater distinction between positive and negative edge weights than PC. (5) The BHM-Ω as a Bayesian version of the traditional CC tends to produce a greater distinction between positive and negative edge weights than CC. Combined with the fact that the high classification accuracy of the BHM method, we can infer that the negative edge weights of BFN also have important information for classification. (6) The BFNs based on the correlation methods show a degree of consistency., as shown in the black box in Figure 2. Similar structures appear in the four BFNs estimated by PC, CC, BHM, which can provide certain support for the reliability of the three correlation methods. In addition, we cluster the brain regions using spectral clustering (Ng et al., 2001) and visualize the clustered adjacency matrices of these 5 methods in Figure 3. It is observed that the BFN estimated by BHM-Ω shows a more significant modular structure.

As stated earlier, some methods including SR and BHM involve optional model parameters. Different parameter values may have a significant impact on the results. Therefore, we calculate the accuracy of different methods under different parameter values, as shown in Figure 4. It is worth noting that, traditional PC and CC models do not involve optional parameters. However, for comparison, we fix their values with the final classification accuracy in Figure 4 for visualization. We can observe that BHM-Ω and BHM-W are quite sensitive to the parameters. When λ in BHM is set to a large real number, the accuracy of BHM-Ω decreases, which may be because the large value of λ, the algorithm has difficulties to converge. Besides, the accuracy of BHM-W increases as λ increases. For this, we empirically tested a larger lambda range [2⁻⁵, 2⁻⁴, ⋯, 2⁰, ⋯, 2⁴, 2⁵] and find that as λ continues to increase, the accuracy decreases, which is consistent with the performance of BHM-Ω. Moreover, SR is not sensitive to different parameter values, but its accuracy performs average in general.

Considering that different p-values significantly influence the results, we show the classification accuracies of 5 methods under different p-values in Figure 5. Note that all 5 methods are sensitive to different p-values. We selected the optimal parameter value for feature selection, so that different methods can get the best classification performance.

In this work, for the ASD classification task, we use the edge weights of the estimated BFN as features. With the empirically optimal parameter, we construct the BFNs using the proposed BHM, then apply a two-sample t-test to rearrange the features according to the p-values. Particularly, we choose the BHM-Ω since it outperforms the BFNs estimated by the other methods. As a result, we obtain the discriminative edge connections with a threshold value p < 0.001 as shown in Figure 6. Here, the thickness of each arc represents the discriminative power that is inversely proportional to the corresponding p-value. The colors of each arc are assigned randomly for better visualization.

In Figure 6, the top discriminative features and the corresponding brain regions, that may contribute to ASD identification include occipital lobe, post-cingulate, dorsal frontal cortex, inferior parietal lobule, precuneus, anterior prefrontal cortex, lateral cerebellum, temporal lobe, fusiform gyrus, mid insula, etc. in order of discriminant ability. The findings are consistent with previous studies (Nickl-Jockschat et al., 2012; Hashem et al., 2020; Lau et al., 2020). We visualize the ROIs using the Brainmesh of Ch2 with Cerebellum in Figure 7, where the size of node spheres depends on the original value in the node file provided by the Dosenbach 160 template.

As described before, we first give an equivalent probability explanation for PC by introducing a normal distribution for the rs-fMRI signal values. Then we reformulate PC with Bayesian rule, thus getting two perspectives of PC. This provides a platform for generalizing PC to CC by assuming that the edge weight matrix W follows the MVN distribution prior. As a result, we derive a probabilistic explanation of CC and develop a high-order BFN estimation framework that allows the introduction of different priors (or regularizers).

Besides the introduced normal distribution prior on W for BFN estimation, we can also introduce other priors on W. For example, considering Laplacian distribution prior for w_ij, e.g., P(wij)=12βe-|wij|β where β is a scale parameter. In this way, the regularized least square problem is

We can find that the Laplacian distribution generates sparse BFN due to the regularizer w_ij. Besides, we can get the optimal solution by the soft thresholding, as follows:

Although different prior distributions can be tried to introduce the proposed probabilistic framework, we do not formulate their models in detail since this paper focuses on the formulation of CC.