Functional Brain Network Estimation with Time Series Self-scrubbing

1). Pearson’s Correlation

As pointed out previously, PC lies on the simplest extreme of the existing FBN methods. We suppose that each brain has been parcellated into N regions of interest (ROIs) based on a certain atlas, and the fMRI time series associated with the ith ROI is represented by x_i ∈ R^T, i = 1, …, N, where T is the number of volumes in each series. Then, the edge weights of the FBN based on PC can be calculated as follows:

By re-defining ${\mathbf{x}}_i\triangleq ({\mathbf{x}}_i-{\overline{\mathbf{x}}}_i)∕\sqrt{{({\mathbf{x}}_i-{\overline{\mathbf{x}}}_i)}^T({\mathbf{x}}_i-{\overline{\mathbf{x}}}_i)}$, a centralized and normalized counterpart of the original x_i, we can simplify PC as W_ij = x_i^Tx_j, which corresponds to the solution of the following optimization problem:

Here, X = [x₁, x₂, ⋯, x_N] ∈ R^T×N is the data matrix, W is the edge weight matrix, and ‖∙‖_F denotes the F-norm of a matrix. Since the BOLD signals commonly contain noises, the original PC tends to result in a FBN with dense connections. In practice, a threshold scheme is generally used to sparsify the PC-based FBN by filtering out the noisy or weak connections. For detailed discussion on the thresholding strategy, please refer to Section 3.2.1 in [36].

2). Partial Correlation

Partial correlation is used in FBN estimation for treating the confounding problem involved in the full correlation methods such as PC. A general approach to calculate partial correlation is based on the estimation of inverse covariance matrix [37]. However, this approach may be ill-posed due to the singularity of the estimated sample covariance matrix Σ = X^TX. Therefore, a regularization term R(W) is generally introduced into the FBN estimation model as follows.

Equivalently, it can be further simplified to the following matrix form:

In Eq. (3) or (4), the first term implies to invert the covariance matrix Σ, and λ is a regularized parameter to control the balance between the first (data-fitting) term and the second (regularization) term. The most popular regularizer is l₁-norm, i.e., $R(\mathbf{W})={\parallel \mathbf{W}\parallel }_1={\sum }_{i,j=1}^n\mid W_{\mathit{ij}}\mid$, which corresponds to the sparsity prior of FBN, and leads to LASSO or sparse representation (SR) model [11]. As discussed earlier, some other regularizers have been investigated in recent years to encode different priors, which is beyond the scope of this paper.

3). A Regularized FBN Estimation Framework

According to a recent study [10], a large family of FBN estimation models can be summarized by the following regularized framework:

where f(X, W) is the data-fitting term, aiming to capture some statistics of the data (e.g., the covariance or inverse covariance structure), and R(W) is the regularization term, aiming to encode the biological/physical priors of the FBN. Sometimes, specific constraints (e.g., symmetry or positive semi-definite) are included in Δ for shrinking the search space of W towards better FBNs.

Note that, such a regularized framework can not only stabilize the statistical estimation, but also, and more importantly, provide a general platform for designing the new FBN estimation method in this paper.

4). Scrubbing Metric

Despite a sophisticated preprocessing pipeline, the fMRI series may be not clean enough due to, for example, head micromovements [16]. In order to eliminate the influence of the instantaneous head motion, a scrubbing operation is typically implemented to remove some volumes according to the frame-wise displacement (FD) [18] or DVARS [38]. In particular, FD of the tth time point is defined as follows.

where x, y, z, α, β, γ are six head motion parameters, and 50 (with the unit of mm) denotes the assumed radius of the head [39]. In general, the volume would be scrubbed if FD exceeds 0.5 mm.

In contrast, the DVARS metric is defined in Eq. (7).

where Y_i,t is the value of the fMRI data at voxel i and time point t. The parameter I is the number of the voxels. The DVARS metric can be nonnalized as follows:

where ${\sigma }_i^2$ is the variance of Y_i,t and ρ_i is the correlation coefficient between Y_i,t and Y_i,t−1. Then, the “bad” volumes can be removed according to the p-value of DVARS2.

1). Motivation

Before presenting the novel FBN estimation model, we first introduce the motivation and the basic idea behind it by a toy example. As shown in Fig. 1, we have a set of rs-fMRI series from N = 5 ROIs, each of which includes T = 30 volumes. According to the previous discussion, we cannot guarantee all volumes in the fMRI series are clean. To illustrate this more clearly, some data points in Fig. 1 are labelled as “clean”, while the others are labelled as “dirty”. Equivalently, a binary indicator can also be used to denote the quality of the data, with “0” corresponding to “dirty”, and “1” to “clean”.

However, the indicator cannot be observed directly in practice. Traditionally, a data-scrubbing operation (e.g., based on FD or DVARS) is used to identify the quality of the data points (or, equivalently, determine the value of the indicator) by detecting head motion [16]. Despite its potential effectiveness, such a scheme (1) cannot remove all the “dirty” data points caused by different factors [23]; (2) tends to scrub too many time points due to the lack of a control mechanism [19]; and (3) is independent on the specific task, and thus will not necessarily benefit the ensuing FBN estimation.

2). Model

To address the above problems, in this paper, we propose a task-dependent data scrubbing method for FBN estimation. In particular, we consider the binary indicator as a latent variable v_t for denoting the quality of the tth time point (v_t = 0 for “dirty”, and v_t = 1 for “clean”), and thus the regularized FBN estimation framework in Eq. (5) can be extended to the following form with an latent indicator variable.

where X^(t) is tth row of the data matrix X. Note that, when the indicator variable v_t = 0, the tth time point in fMRI series will be removed, meaning that it has no contribution to FBN estimation; when v_t = 1 for all t = 1, 2, ⋯, T, Eq. (9) will reduce to the original FBN estimation framework given in Eq. (5).

In practice, however, Eq. (9) will always return a trivial solution that v_t = 0 for all t = 1, 2, ⋯, T. Therefore, we introduce a negative l₁-norm regularization term $-\gamma {\sum }_{t=1}^T\left|v_t\right|$ into Eq. (9) for avoiding an exceedingly sparse solution. Additionally, we relax⁴ the binary indicator variable v_t ∈ {0,1} to a real range [0,1], for simplifying the solving of the optimization problem. Them we get the following model:

where λ and γ are two regularized parameters for controlling the balance of three terms in the objective function. Especially for γ, it plays a role to determine the number of removed time points from the whole series. When the value of γ approaches 0, most of the time points will be removed. On the other hand, when γ has a large value, all the time points in the series will be kept for FBN estimation. In other words, the proposed model can scrub the data adaptively in the process of FBN estimation, by controlling the hyper-parameter γ and learning the indicator variable v_t from the data. Therefore, we name our proposed scheme as FBN estimation with self-scrubbing (SS for short).

In principle, any data-fitting/regularization term can all be adopted to realize the SS model, but, in this paper, we select the partial correlation strategy shown in Eq. (3) or (4), since it overcomes the confounding effect and is empirically verified to more effective than the full correlation [9]. Consequently, the partial correlation based self-scrubbing FBN estimation model is given as follows:

where $x_i^{(t)}$ is the tth time point of the series associated with the ith node (e.g. ROI). For simplicity, we rewrite Eq. (11) into the following matrix form:

where V = diag(v₁, v₂, ⋯, v_T) ∈ R^T×T is a diagonal matrix containing the indicator variables on its principal diagonal. The constraint, W_ii = 0, is employed only to avoid the trivial solution that leads to W being an identity matrix. For R(W), we can, in principle, use any off-the-shelf regularizers such as l₁-norm [11], l_2,1-norm [13], trace norm and their combination [10]. However, this problem goes beyond the main focus in this paper. Therefore, we only attempt l₁-norm (sparsity prior) due to its simplicity and effectiveness [11] and get the specific FBN estimation model (named SR+SS) as follows.

3). Optimization Algorithm

Considering that there are two variables V and W involved in Eq. (12), in this paper, we employ alternative convex search (ACS) [40] method to solve them alternately, by two following steps. Note that, we first initialize indicator V as an identity matrix, meaning that all time points are retained in first iteration.

Step 1: With a fixed V, Eq. (13) reduces to a traditional FBN estimation problem that can be solved by many convex optimization methods. Here, we use the proximal method [41] due to its efficiency and simplicity. Two main steps are involved in the proximal method, including gradient descent and proximal operation. First, for the data-fitting term $f(\mathbf{X},\mathbf{W})={\parallel \mathbf{VX}-\mathbf{VXW}\parallel }_F^2$, whose gradient w.r.t W is ∇_Wf(X, W) = 2X^TV^TVXW − X^TV^TVX, we have the following update formula according to the gradient descent criterion:

where α_k denotes the step size of the gradient descent⁵ [42]. Then, the proximal operator is imposed on the current W. For the sparsity regularizer λ‖W‖₁, the proximal operator is defined as follows.

where sgn(W_ij) and abs(W_ij) return the sign and absolute value of W_ij respectively.

Step 2: With a fixed W, we then update V. Now, Eq. (13) reduces to the following optimization problem,

or

Eq. (17) can be further simplified to the following problem.

where $c_t={\sum }_{i=1}^N{(x_i^{(t)}-{\sum }_{j\ne i}^n{W}_{\mathit{ij}}{x}_j^{(t)})}^2={\parallel {\mathbf{X}}^{(t)}-{\mathbf{X}}^{(t)}\mathbf{W}\parallel }^2$ (i.e., loss value of the data-fitting term) is a constant. Note that Eq. (18) is a linear programming problem, and thus we can easily get its optimal solution as follows.

This means that 1) if ‖X^(t) − X^(t)W‖² > γ, the tth time point is more likely to be removed by labelling it with 0; on the contrary, 2) if ‖X^(t) − X^(t)W‖² < γ, the tth time point will be kept by labelling it with 1. Such a formula for updating v_t coincides well with the intuition that the “dirty” time points (labelled as 0) can cause a poor fitting with high bias/residual (i.e., ‖X^(t) − X^(t)W‖² > γ ). In fact, in the experimental section, we will further illustrate this problem based on a toy dataset, and empirically verify that the proposed method can automatically detect “dirty” time points. Finally, we summarize the algorithm for solving Eq. (13) in ALGORITHM I.

1). FBN estimation

After obtaining the preprocessed fMRI data, we estimate FBNs based on three different methods, PC, SR, and the proposed SR+SS. In our experiments, the average convergence times of SR and SR+SS are 0.0861s and 0.0916s per network, respectively, based on Win 10 OS, core i7 and MATLAB 2016a. The convergence time is calculated by tic and toc tool in MATLAB. This means that the proposed method does not spend much more running time than the baseline method.

In Fig. 3, we visualize the adjacency matrices⁷ of the FBN estimated by these methods on ADNI dataset. For SR, we simply set the regularized parameter λ = 1 and for SR+SS, λ = 1 and λ = 0.5. As shown in Fig. 3, the FBN estimated by PC has a highly different topology from that of the partial correlation-based methods (i.e., SR and SR+SS), since they use different data fitting term. In contrast, SR and SR+SS lead to a similar FBN structure by using the same kind of data fitting term. The differences of the FBNs estimated by SR and SR+SS methods lie mainly in several specific brain regions.

For investigating the effect of SS on the estimated FBN, in Fig. 4, we map the most significantly changed connections between SR and SR+SS onto the International Consortium for Brain Mapping (ICBM) 152 surface template. For a better visualization, we do not provide the cerebellum part. It can be observed, in Fig. 4, that the changed connections mainly concentrate on the hippocampus, parahippocampus, cuneus and temporal lobe regions. In Part 4), we will give a further discussion on the self-scrubbing scheme.

2). Experimental Settings

Once we obtain the FBNs of all subjects, the subsequent task is to identify the subjects with MCI from NCs based on the estimated FBNs. Then, the problem turns to determine which features and classifiers should be used for identification. Considering the big influence of different steps in the classification pipeline on the final accuracy, it is difficult to conclude whether the FBN estimation methods or the ensuing feature selection and classification methods contribute to the ultimate result. Therefore, we only adopt the simplest feature selection method (t-test with p<0.01) and the most popular SVM [44] classifier (linear SVM with default parameter C=1) in our experiment.

Note that, for comparing the proposed method with the traditional FD-based scrubbing method, we conducted experiments on both scrubbed and non-scrubbed data. The scrubbing criterion is based on the commonly used scheme with FD > 0.5mm.

Due to limited samples, we test the involved methods using the leave one out cross validation (LOOCV), in which only one subject is kept out for testing while the remaining are used for training the models. This is repeated until each fold has been used once for testing. For obtaining optimal values of the regularized parameters, an inner LOOCV is further conducted on the training data by a grid-search. For the parameter λ, the candidate value is ranged in [2⁻⁵, 2⁻⁴, ⋯, 2⁰, ⋯, 2⁴, 2⁵]; for the parameter γ , the candidate value is ranged in [0.1,0.2,⋯,0.9,1]; for the threshold in PC, we use 20 sparsity levels from [5%, 10%, ⋯, 95%, 100%], where, for example, 90% means that 10% of the weak edges are filtered out from the FBN.

In addition, we evaluate the classification performance of different methods by a set of quantitative measures, including accuracy, sensitivity and specificity, which are defined as follows:

where TP, TN, FP and FN denote the number of True Positive, True Negative, False Positive and False Negative, respectively.

3). Results

The classification results on ADNI dataset (with and without traditional scrubbing⁸) are reported in TABLE I. The corresponding classification results on NITRC dataset are given in TABLE II. Based on the results, we can observe that the proposed SR+SS method achieves the best discriminability on both scrubbed and non-scrubbed datasets. The proposed methods significantly out perform the baseline under the 95% confidence interval based on the DeLong’s non-parametric statistical test [45]. The p-value of the DeLong test is summarized in TABLE III. The receiver operating characteristic (ROC) is given in Fig. 5. It can be clearly found that the proposed method performs better than the traditional scrubbing method and can further improve the performance on the scrubbed data.

It is worth noting that compared to PC and SR, the proposed SS scheme has an extra hyper-parameter γ for controlling the number of the removed volumes. In Fig. 6, we count the removed volumes across subjects for different values of the hyper-parameter and find no significant differences between two groups (i.e., MCI and NC). This means that the self-scrubbing operation itself does not introduce bias into MCI and NC data. In other words, the improvement of the discriminability may mainly benefit from the proposed FBN estimation method.

4). Further Discussion on the Self-Scrubbing Scheme

In this section, we designed two experiments for investigating the scrubbed volumes by our proposed method. First, we conduct an experiment to investigate the removed volumes by the traditional scrubbing (DVARS) and our proposed self-scrubbing scheme. Based on the result shown in Fig. 7, we note that they can remove some common volumes with higher possibility than the random case. For example, if 30 volumes were removed from the time series with totally 137 volumes, nearly 13 volumes are common, on the average, by DVARS and the proposed method. However, the number of common volumes is only 5 for the random case. This result illustrates that the scrubbed volumes by the proposed criterion contain a part of stereotyped structured noisy points such as head motion.

Second, in order to further explore the removed volumes by our scheme, we estimate FBNs by the time series with and without SS⁹, and then determine the brain regions (ROIs) with connections influenced significantly by SS. For each brain region, we add up their corresponding connection weights, and the sum for the ith brain region is calculated as follows:

where N is the number of ROIs, R is the number of subjects, W_ij is the connection weight value between ROIs i and j in the estimated FBN and k is the k-th subject. For each brain region, we have a series with R elements.

Then, we investigate the statistical differences of brain regions with and without SS by pair-wise t-test. The result is mapped onto the ICBM 152 template as shown in Fig. 8, where the color bar represents the p-value of t-test. In particular, under the ICBM 152 template and AAL atlas, we simply set the value of corresponding ROI as the p-value (p<0.05). In Fig.8, we note that the differences of the estimated FBNs mainly concentrate on hippocampus, temporal lobe, cuneus and frontal regions. These brain regions have been shown to be closely related to memory, emotion and auditory sense [46–49]. Therefore, we assume that those removed points may be related to some abnormal resting-state processes such as mind-wandering (memory), worrying (emotion), or distraction due to scanner sounds [23]. In other words, the time points that do not fit the FBN estimation may be related to some “disturbances” from the scanning conditions, memory, emotion or mind wandering. That is, compared with the traditional motion-based scrubbing methods, the proposed method can further scrub some “abnormal resting-states” points.

5). Discriminative Connection

For FBNs based on SR+SS, we investigate the most discriminative connections for identifying subjects with MCI from NCs. Specifically, we select 64 connections whose p-value are less than 0.01. As visualized in Fig. 9, the thickness of the arc is inversely proportional to the p-value, indicating the discriminative power of the corresponding edge. It is worth pointing out that the significantly changed connections or regions, as shown in Fig. 4 and 8, are contained in these discriminative connections. Therefore, it can support the classification accuracy improvement of the proposed method to some extent. Additionally, we note that several regions, including superior-medial frontal gyrus, medial orbitofrontal gyrus, hippocampus, para-hippocampus and precuneus, are selected in our proposed method. These regions are generally involved in the default mode network [50], and believed to be biologically associated with MCI identification, according to previous studies [51,52].