This paper proposes a simple yet effective approach for detecting activated voxels in fMRI data by exploiting the inherent sparsity property of the BOLD signal in temporal and spatial domains. In the time domain, the approach combines the General Linear Model (GLM) with a Least Absolute Deviation (LAD) based regression method regularized by the pseudonorm
The medical imaging modality known as functional Magnetic Resonance Imaging (fMRI) using blood oxygen leveldependent (BOLD) contrast is a noninvasive technique widely accepted as a standard tool for localizing brain activity [
Since its development in the early 1990s, a variety of univariate and multivariate methods for analyzing fMRI data have been developed [
In addition to these drawbacks, there are also a variety of research studies that have pointed out the existence of a common underlying principle involved in sensory information processing referred to as “sparse coding” [
In the last decade, there has been a growing interest in the use of sparse signal representation techniques for analyzing fMRI data based on the assumption that the components of each voxel’s fMRI signal are sparse and the neural integration of those components is linear. These studies include, among others, sparse representation of nonparametric components of partially linear models [
In the context of sparse regression over the coefficients of the GLM, the solution of the regularized inverse problem induced by the GLM is obtained by using
Motivated by recent developments in sparse signal representation and the biological findings of ‘sparse coding’ in the brain, in this paper, we propose a simple yet effective approach based on the sparsity of underlying BOLD signal in fMRI data that exploits both temporal and spatial sparse properties of the fMRI images. This paper has two aims, the first to explore the ability of a new robust regression method named
Our rationale in spatial domain is that the set of activated voxels in response to a given task is spatially sparse in the sense that only a reduced number of voxels that comprise the brain volume are activated. In fact, the behavior exhibited by different sets of spatial parameters associated with different stimuli of the
The assumption of a Laplacian distribution model for explaining the set of spatial parameters has been proposed [
A very preliminary version of this manuscript was presented at the 33rd Annual International Conference of the IEEE Engineering in Medicine and Biology Society which focused mainly on exploiting the sparsity of the BOLD signal in the time domain whereas in the spatial domain the activation maps were generated by selecting the 300 most significant spatial parameters [
This paper is organized as follows. In Section
The GLM for the observed response variable
In this work, the solution to the inverse problem (
An important component in the GLM is the design matrix
Given the linear model,
The
Once the regression parameters are determined following the approach described above, each voxel is represented by a
As was mentioned above, the reasons that drive us to use a Laplacian model as a sparsityinducing model are based on the following considerations.
Assuming that the probability distribution of the entries of
Given the critical lack of ground truth, a major difficulty with fMRI research concerns the validation of any analysis of experimental
In all experiments the dictionary is designed according to the structure
The synthetic data were generated by blending a simulated activation in nonactivated time series that were extracted from the PBAIC2007 fMRI database. Thus, knowledge of the ground truth is assured while the noise is representative for real data. According to this, the synthetic time series
(a) Synthetic boxcar signal
Likewise, in order to reduce the number of voxels to be analyzed, instead of trying to model the cerebral cortex a synthetic volume was designed. This volume (depicted in Figure
A synthetic volume consists of 4 slices. White regions denote the set of voxels that has been activated by the stimulus
The synthetic fMRI database was generated from 32000 nonactivated time series in the following steps. First, the 32000 nonactivated time series were divided into 20 sets of 1600 time series each. Then, each set of 1600 time series was inserted into the synthetic volume. Finally, for each one of the 108 selected voxels (white regions) the synthetic activation was added to their BOLD signals, while remaining voxels retain their original real BOLD signals. Using 4 different levels of activation strength
Signaltonoise ratio of the synthetic dataset.
Group  1  2  3  4 

SNR  0.1419  0.2838  0.4256  0.5675 
SNR (dB) 




Regarding the design of the dictionary, in this experiment
The PBAIC2007 fMRI dataset was collected with a Siemens 3T Allegra scanner using an EpiBOLD sequence, with imaging parameters TR and TE being set to 1.75 s and 25 ms, respectively. Subjects were engaged in a Virtual Reality task during which they had to perform a number of tasks in a hypothetical neighborhood. The tasks included, among others, the acquisition of pictures of neighbors with piercing, the gathering of specific objects (e.g., fruits, weapons) and the avoidance of a growling dog (for a more detailed description of the experiment, see [
Given the acceptance of the SPM software in practice [
In this experiment the dictionary is designed following the framework used by SPM [
In the second experiment, predictors that model confounds are incorporated to the dictionary. Specifically, the 13 DCT basis functions defined by (
In order to identify the voxels activated by the stimulus
Given the knowledge of the ground truth, the proposed method is evaluated by comparing the set of voxels detected as activated against the set of true activation. From this information, it is determined the number of true detections, false alarms, and missed activated voxels. Figure
Average of activated voxels, true detections, false alarms, and missed activations for synthetic data base.
Group  Activated voxels  True detections  False alarms  Missed activated voxels 


72.40  67.45  4.95  40.55 

102.80  102.25  0.55  5.75 

106.25  106.25  0  1.75 

107.20  107.20  0  0.8 
(a) Activation maps (left to right: slice 1 to slice 4; top to bottom: group 1 to group 4) for stimulus
To illustrate the performance of the proposed method in detecting activation we choose the data of subject 14, run 1 and two representative tasks: Instructions and Faces. The selection of these particular tasks is based on the a priori knowledge of the anatomical and functional localization of auditory and visual cortex. With regard to the parameters
Figures
Activation maps (left to right, top to bottom: slices 11 to 22) for the Instructions task obtained with (a) proposed method with reduced dictionary (14 atoms) and (b) SPM software with FWE rate.
Activation maps (left to right, top to bottom: slices 9 to 19) for the Faces task obtained with (a) proposed method with reduced dictionary (14 atoms) and (b) SPM software with FWE rate.
Under these experimental conditions, for the Instructions task, the proposed method activates 676 voxels contained in 26 slices (2 to 26 and 28) and SPM activates 708 voxels contained in 22 slices (2 to 23). Although the activated slices are not exactly the same, the number of matching activated slices is high, to be more precise 22 (slices 2 to 23) for a matching percentage of 84.61% at slices level. At voxels level, however, the percentage of matching is 57.99%. It can be seen in Figure
Figures
Results of experiments 1 and 2 indicating number of activated voxels (av), number of activated slices (as), and the matching percentages at voxel level (mpvl) and at slice level (mpsl) with respect to SPM software.
Stimulus  SPM  Reduced dictionary  Extended dictionary  

av  as  av  as  mpvl  mpsl  av  as  mpvl  mpsl  
Instructions  708  22  676  26  57.99%  84.62%  717  22  63.18%  77.27% 
Faces  91  11  95  17  54.96%  64.71%  120  17  47.50%  64.71% 
Activation maps (left to right, top to bottom: slices 11 to 22) for the Instructions task generated by (a) the proposed method with extended dictionary (26 atoms), and (b) SPM software with FWE rate.
Activation maps (left to right, top to bottom: slices 9 to 19) for the Faces task generated by (a) the proposed method with extended dictionary and (b) SPM software with FWE rate.
Activation maps for the Instructions task generated with the proposed method by using the reduced dictionary (a) and extended dictionary (b). The activated regions contained in the ellipses show the effect of “filling” observed when the parameter estimation is performed using the extended dictionary.
Activation maps for the Faces task generated with the proposed method by using the reduced dictionary (a) and extended dictionary (b). The activated regions contained in the ellipses show the effect of “filling” observed when the parameter estimation is performed by using the extended dictionary.
In this paper, a new method for detecting activation in fMRI data that exploits the sparseness of the BOLD signal in both time and spatial domains is presented. By adopting a sparse approach in both domains two objectives are achieved. First, in time domain, the
As mentioned in Section
Data analysis through histograms strongly suggests as a plausible model the Laplacian distribution mainly due to the concentration of mass centered around zero and the heaviness of the tails of the empirical distribution. Nevertheless, two probability distribution models are fitted to data: specifically, the Gaussian model which has been the standard model assumed for spatial analysis of fMRI data and the Laplacian model. For both models, maximum likelihood (ML) estimators are used to compute the respective distribution’s parameters (i.e., the location and dispersion parameters) using the spatial parameters
Statistical behavior of spatial parameters set
QQplots for Gaussian and estimated distributions of stimuli: (a) Dog, (b) Faces, (c) Instructions, and (d) WeaponsTools.
Since graphical tools lack rigor, a quantitative measure of the distance between the estimated distribution and the proposed models based on the AIC is used. The AIC is defined by
Sign of the difference
Previous analysis shows that the tails of the estimated distribution are generally heavier than those of a Gaussian distribution, but not as heavy as those of a Laplacian distribution. This finding suggests the feasibility of adopting a more general statistical model: the Generalized Gaussian (GG) distribution, which contains the Gaussian and Laplacian distributions as special cases. The GG density has the following form [
In case of synthetic database, the shape parameter estimated by the function
Results of analyzing in vivo data with the normalp package indicating the shape parameter
Stimuli 




Arousal  1  0.00979738  0.00979711 
Dog  1  0.1483275  0.14082884 
Faces  1  1.24216047  1.24212600 
FruitsVegetables  1  0.1328233  0.13281962 
Hits  1  0.39250583  0.39249493 
Instructions  1.3637  0.35744665  0.41064940 
InteriorExterior  1  0.13646107  0.13645728 
SearchFruit  1  0.04498026  0.04497901 
SearchPeople  1  0.01944126  0.01944072 
SearchWeapons  1  0.03387324  0.03387230 
Valence  1  0.04011155  0.04011044 
Velocity  1.7865  0.1221627  0.18050036 
WeaponsTools  1  0.40075143  0.40074030 
A partial version of this manuscript was presented as a plenary at the 33rd Annual International Conference of the IEEE Engineering in Medicine and Biology Society in Boston, in 2011 [
The authors declare no conflicts of interest.
The authors would like to acknowledge ECOSNORDFONACIT under Grant PI20100000299, the Council of Scientific, Humanistic, Technological and Artistic Development of the Universidad de los Andes (CDCHTAULA) under Project I13361202B, and the Universidad Nacional Experimental del Táchira (UNET) for their financial support.