Eneko's Digital Garden

A Low Rank and Sparse Paradigm Free Mapping Algorithm for Deconvolution of fMRI Data

Motivation

The performance of existing deconvolution approaches can be hampered considerably in presence of global, widespread signal changes due to head jerks, hardware artifacts, or prominent non-neuronal physiological events (e.g. deep breaths) (Power et al. (2017)).

These global events are difficult to compensate during data preprocessing (Power et al. (2018)) and can be misinterpreted as neuronally related since their temporal signature can closely resemble the hemodynamic response function (HRF) assumed by the deconvolution model to describe neurovascular coupling.

We propose a new approach (LR+MV-SPFM) for the spatio-temporal deconvolution of fMRI data that can simultaneously estimate global signal fluctuations and neuronal-related activity based on the low-rank plus sparse matrix decomposition method (Otazo et al. (2015)).

Signal model

The whole-brain fMRI data $\mathbf{Y} \in \R^{N \times V}$ , where $N$ is the number of volumes and $V$ is the number of voxels, can be decomposed into three terms: the neuronal-related component ( $\mathbf{HS}$ ), low-rank components ( $\mathbf{L}$ ) and additional Gaussian noise ( $\mathbf{N}$ ).

$\mathbf{Y} = \mathbf{HS} + \mathbf{L} + \mathbf{N}$

$\mathbf{HS}$ is the convolution of voxel-specific activity-inducing signal ( $\mathbf{S}$ ) with the Toeplitz matrix $\mathbf{H} \in \R^{N \times N}$ with shifted HRFs in its columns (Caballero Gaudes et al. (2013)).

global fluctuations can be captured as the sum of $P$ spatially-widespread (i.e., global) low-rank components $\mathbf{L} = \sum_{p=1}^P v_p a_p^T$ , where $v_p \in \R^{N}$ and $a_p \in \R^V$ denote their corresponding spatial and temporal signatures.

Algorithm

In order to simultaneously estimate the activity-inducing signal $\mathbf{S}$ and low-rank components $\mathbf{L}$ , the following multivariate regularized least-squares problem must be solved:

$\hat{\mathbf{L}}, \hat{\mathbf{S}} = \argmin_{\mathbf{L, S}} | \mathbf{Y} - \mathbf{HS} \mathbf{L} |_F^2 + \lambda_L |\mathbf{L}|_* + (1-\rho) | \mathbf{SD}_S |_{2,1} + \rho |\mathbf{SD}_S |_1$

$| \cdot |_F^2$ denotes the Frobenius norm.

The $\ell_{2,1} + \ell_1$ norm promotes temporal sparsity and spatial structure to the estimates of the activity-inducing signal.

$\rho = 0.8$ (empirically set) balances the trade-off between temporal sparsity and spatial structure.

$\mathbf{D}_S = diag(\lambda_{S_1}, ..., \lambda_{S_V})$ is a diagonal matrix with voxel-specific regularization parameters to promote the sparsity of $\mathbf{S}$ .

For each voxel, we set $\lambda_S$ to the median absolute deviation estimate of the noise standard deviation from the fine-scale wavelet coefficients of the voxel time series (Daubechies, order 3).

The nuclear norm $| \cdot |_*$ promotes the detection of low-rank, global components.

We solve the regularization problem by means of monotone FISTA with variable acceleration (Zibetti et al. (2018)).

Results

Simulated data

We simulated 1000 voxels including two groups of 50 voxels with a known BOLD signal, whereas the remaining voxels did not contain any BOLD signal. For each voxel, we added noise of different sources (motion-related, thermal and physiological noise) (Caballero Gaudes et al. (2013)), as well as two global low-rank components with a random voxelwise amplitude simulating widespread signal changes due to two deep breaths (Power et al. (2018)) and motion-related spikes, respectively.

Regardless of the simulated SNR and the BOLD/total voxels ratios, the ROC curves demonstrate the proposed LR+MV-SPFM algorithm is more specific and provides higher sensitivity in comparison with the original SPFM method, except with the highest BOLD/total number of voxels ratio and highest SNR where it seems that the multivariate nature of the model prevents the algorithm from fitting accurately at the voxel level.

The estimation of the low-rank components improves when the BOLD-only/total voxels ratio decreases for $\rho = 0.8$ .

Experimental data

For this subject and session, the proposed approach estimated $P = 3$ low-rank components (LRC). Among participants, the number of low-rank components ranged between 1 and 5.

The Euclidean norm of the motion displacements (E-norm), DVARS, and global signal (GS) time series.

Grayplots of the preprocessed data (RAW), estimated low-rank components and estimated neuronal-related components in grey matter (GM) and white matter (WM) voxels.

LRC1 captures signal fluctuations related to head movements and susceptibility artefacts during the "moving the tongue" condition. LRC2 has a time series that closely follows the global signal and its spatial map actually delineates major arteries and draining veins that strongly contribute to the global signal. LRC3 is clearly related to global physiological fluctuations.

The estimated activity-inducing signal correctly matches the timings of the tongue condition.

LR+MV-SPFM maps reveal clusters of activity in similar regions to those inferred with general linear model analyses ( $p < 0.001$ ). LR+MV-SPFM maps still depict the tongue areas of the motor cortex bilaterally despite the timing of the first low-rank component followed the tongue condition.

LRC1 illustrates motor areas that show a pattern similar to the activity maps in the tongue condition, while inferior regions of the brain show patterns that match with LRC3. LRC2 spatial map delineates major arteries and draining veins that strongly contribute to the global signal. LRC3 map shows inferior parts of the brain that are highly susceptible to physiological fluctuations.

ROC curves of the five motor task conditions show that MV-SPFM and LR+MV-SPFM provide higher sensitivity at the cost of reduced specificity.The higher complexity involved in estimating the low-rank components does not reduce the accuracy in deconvolving the neuronal-related component of the signal.

Take-home messages

We present a novel formulation for the deconvolution of BOLD fMRI data that captures global fluctuations due to motion artifacts or physiological signals.

Capturing global fluctuations does not reduce accuracy of the estimated neuronal-related activity.

Future work:

Once demonstrated that

Capturing global fluctuations does not reduce accuracy of the estimated neuronal-related activity.

, we will evaluate the algorithm on resting-state data., where the low-rank components should only describe motion artifacts and physiological signals.

We will explore robust approaches for selecting $\lambda$ .

The selection of the regularization parameter $\lambda$ is critical for the precise estimation of the neuronal-related and low-rank components.

We will work on extending the proposed formulation for multi-echo fMRI acquisitions.

We will work on describing the neuronal-related signal in terms of its innovation signals (the temporal derivative of the activity-inducing signal), which is sparser and better fits the sparse nature of the regularized least-squares problem proposed*.

Referenced in

my contributions

E. Uruñuela, S. Moia, C. Caballero-Gaudes. A Low Rank and Sparse Paradigm Free Mapping Algorithm for Deconvolution of fMRI Data. IEEE 18th International Symposium on Biomedical Imaging (ISBI 2021).

A Low Rank and Sparse Paradigm Free Mapping Algorithm for Deconvolution of fMRI Data