LEICA: Laplacian eigenmaps for group ICA decomposition of fMRI data

Highlights

• Introduction of a non-linear model for analyzing multi-subject fMRI datasets.

• Model constructs a sparse correlation matrix, applies Laplacian eigenmaps and ICA.

• Evaluated on two fMRI datasets. Results are favorable compared to other models.

• Our components are more functionally cohesive and reproducible than other models.

• The algorithm is robust to noisy data and computationally efficient.

Abstract

Independent component analysis (ICA) is a data-driven method that has been increasingly used for analyzing functional Magnetic Resonance Imaging (fMRI) data. However, generalizing ICA to multi-subject studies is non-trivial due to the high-dimensionality of the data, the complexity of the underlying neuronal processes, the presence of various noise sources, and inter-subject variability. Current group ICA based approaches typically use several forms of the Principal Component Analysis (PCA) method to extend ICA for generating group inferences. However, linear dimensionality reduction techniques have serious limitations including the fact that the underlying BOLD signal is a complex function of several nonlinear processes. In this paper, we propose an effective non-linear ICA-based model for extracting group-level spatial maps from multi-subject fMRI datasets. We use a non-linear dimensionality reduction algorithm based on Laplacian eigenmaps to identify a manifold subspace common to the group, such that this mapping preserves the correlation among voxels’ time series as much as possible. These eigenmaps are modeled as linear mixtures of a set of group-level spatial features, which are then extracted using ICA. The resulting algorithm is called LEICA (Laplacian Eigenmaps for group ICA decomposition). We introduce a number of methods to evaluate LEICA using 100-subject resting state and 100-subject working memory task fMRI datasets from the Human Connectome Project (HCP). The test results show that the extracted spatial maps from LEICA are meaningful functional networks similar to those produced by some of the best known methods. Importantly, relative to state-of-the-art methods, our algorithm compares favorably in terms of the functional cohesiveness of the spatial maps generated, as well as in terms of the reproducibility of the results.

Introduction

Functional magnetic resonance imaging (fMRI) studies have shed light on the overall functional organization of the brain. For example, resting state studies have revealed the existence of a number of intrinsic functional networks, including the “default mode” and “salience” networks Biswal et al., 1995, Cole et al., 2014, Greicius et al., 2003. These studies have also led to the discovery of specific patterns related to brain disorders Jafri et al., 2008, Greicius et al., 2004, Garrity et al., 2007, age Damoiseaux et al. (2008), and gender Weissman-Fogel et al. (2010).

The determination of large-scale patterns of brain activity is a challenging problem due to several factors, including the high-dimensionality of the data, various noise sources, and inter-subject variability. The most widely used approach to capture group-level functional connectivity is based on Independent Component Analysis (ICA) McKeown et al. (1998), which is a data-driven approach that assumes the existence of statistically independent latent features, called sources, which can generate the data through linear mixtures. However, ICA does not naturally provide a suitable method for drawing inferences about groups of subjects Erhardt et al. (2011), which has led to the development of several multi-subject approaches that are partly based on ICA. These approaches revolve around a few strategies. One method is to concatenate the data from multiple subjects along the temporal dimension Calhoun et al. (2001) followed by the use of Principal Component Analysis (PCA) Jolliffe (2002) to reduce the dimensionality of the data, and ending with the application of ICA to extract spatial features common at the group level. The tensorial extension to ICA Beckmann and Smith (2005) uses a three-dimensional tensor to estimate shared spatial patterns and time courses between subjects, and subject-specific loadings of the components to capture the multidimensional structure of the data. In Guo and Pagnoni (2008), a general framework is proposed for using the expectation maximization algorithm to estimate an unconstrained mixing matrix for the group-level independent components. Finally, canonical ICA Varoquaux et al. (2010) uses Canonical Correlation Analysis (CCA) to find components that maximize cross-correlation before applying the ICA step.

Two methods involving temporal concatenation of the data include GIFT Calhoun et al. (2001) (http://icatb.sourceforge.net/) and MELODIC Beckmann and Smith (2004) (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/MELODIC/). GIFT approximates the PCA of the group temporally concatenated data by first reducing each subject's dataset to dimension m, and then running the PCA on the concatenation of the principal components across all subjects to further reduce temporal dimensionality. MELODIC, however, uses an incremental approach called MELODIC's Incremental Group PCA (MIGP) Smith et al. (2014) to implement the approximation. Although MELODIC and GIFT have been broadly used for group level fMRI analysis, they are not without limitations as will be explained next.

These ICA-based strategies involve the use of PCA to reduce the dimensionality of the data. PCA is basically a linear projection of the original data onto a low-dimensional subspace in such a way as to retain most of the variance of the data. However, PCA does not distinguish between variance due to noise versus variance due to inherent underlying signal variations. Moreover, there is no principled justification for using the variance as the basis for data reduction based solely on a linear transformation. The BOLD signal is a complex function of neural activity, oxygen metabolism, cerebral blood volume, cerebral blood flow, and other physiological parameters Wang et al. (2003). The dynamics underlying neural activity and hemodynamic physiology involve multiple nonlinearities Birn et al., 2001, Miller et al., 2001, Xie et al., 2008, a fact that has motivated the exploration of nonlinear approaches for fMRI analysis Thirion and Faugeras, 2004, Mannfolk et al., 2010. To the best of our knowledge, nonlinear models have not been proposed to extract common spatial patterns for multi-subject fMRI analysis.

Several nonlinear dimensionality reduction approaches have been successfully used to address a wide range of applications Van Der Maaten et al. (2009). Kernel PCA Schölkopf et al. (1998) maps data to a feature space using a nonlinear transformation (kernel) and performs linear PCA in the feature space. Isomap Tenenbaum et al. (2000) is based on computing a low dimensional representation that tries to preserve the pairwise geodesic distances, which are approximated by the pairwise shortest paths between corresponding points. Maximum variance unfolding (MVU) Weinberger and Saul (2006) computes a low-dimensional manifold such that local distances and angles are preserved. Locally linear embedding (LLE) Roweis and Saul (2000) is based on the geometric intuition of local linearity, which assumes that each point and its neighbors lie on an approximately linear patch of a low-dimensional manifold. The mapping is computed by trying to preserve linearity locally, such that each point has the same neighborhood structure as in the high-dimensional space. A technique of a different flavor is the t-distributed stochastic neighbor embedding (t-SNE) Maaten and Hinton (2008), which captures similarities between pairs of high-dimensional points through a probability distribution, followed by determining a mapping to a lower-dimensional space that achieves a similar distribution. Laplacian eigenmaps Belkin and Niyogi (2003) project data into a low dimensional Euclidean space such that local proximity relations are preserved as much as possible, mapping close points in the original space to close points in the low-dimensional space. Of these techniques, Laplacian eigenmaps appear to be the most suitable for fMRI data since its overall objective function can be expressed as minimizing

$$\sum _{ij}{\Vert {\mathit{y}}_{\mathit{i}}-{\mathit{y}}_{\mathit{j}}\Vert }^2{W}_{ij}$$

where ${\mathit{y}}_{\mathit{i}},{\mathit{y}}_{\mathit{j}}$ are the low-dimensional vectors corresponding to points $i,j$ in the original space, and $W_{ij}$ is a proximity measure of points i and j in the original space. If we choose $W_{ij}$ to be the zero-thresholded correlation between voxels i and j in the original space, this nonlinear projection maps voxels with strong correlations to nearby points in the low-dimensional Euclidean space. Correlations between voxels in fMRI data reveal functional relationships and hence contain the information that can be used to construct functional networks.

We note that Laplacian eigenmaps have been successfully used in analyzing fMRI data Marquand et al., 2017, Haak et al., 2016. Our approach differs from earlier studies in that they are ROI-based methods focusing on specific networks using prior brain-related knowledge, while our framework makes use of ICA to extract intrinsic networks for general fMRI datasets.

In this paper, we propose a group-level model that uses Laplacian eigenmaps as the main data reduction step, which preserves the correlation information in the original data as best as possible. The nonlinear map is robust relative to noise in the data and to inter-subject variability. After the Laplacian eigenmaps transformation, ICA is applied to the reduced data to estimate spatial features, followed by thresholding to extract group-level independent spatial maps. We call our algorithm LEICA (Laplacian Eigenmaps for group ICA decomposition).

The rest of the paper is organized as follows. In section 2, we present the details of our algorithm and inference method, and describe several tests used to validate our method and compare the results with state-of-the-art methods. In section 3, we show the corresponding test results on both resting state and working memory task fMRI datasets for 100 subjects from the Human Connectome Project. Finally, in section 4, we discuss our results.