Total activation: fMRI deconvolution through spatio-temporal regularization

Abstract

Confirmatory approaches to fMRI data analysis look for evidence for the presence of pre-defined regressors modeling contributions to the voxel time series, including the BOLD response following neuronal activation. As more complicated questions arise about brain function, such as spontaneous and resting-state activity, new methodologies are required. We propose total activation (TA) as a novel fMRI data analysis method to explore the underlying activity-inducing signal of the BOLD signal without any timing information that is based on sparse spatio-temporal priors and characterization of the hemodynamic system. Within a variational framework, we formulate a convex cost function—including spatial and temporal regularization terms—that is solved by fast iterative shrinkage algorithms. The temporal regularization expresses that the activity-inducing signal is block-type without restrictions on the timing nor duration. The spatial regularization favors coherent activation patterns in anatomically-defined brain regions.

TA is evaluated using a software phantom and an event-related fMRI experiment with prolonged resting state periods disturbed by visual stimuli. The results illustrate that both block-type and spike-type activities can be recovered successfully without prior knowledge of the experimental paradigm. Further processing using hierarchical clustering shows that the activity-inducing signals revealed by TA contain information about meaningful task-related and resting-state networks, demonstrating good abilities for the study of non-stationary dynamics of brain activity.

Introduction

Conventional analysis of functional magnetic resonance imaging (fMRI) data heavily relies on approaches based on general linear models (GLMs) where prior knowledge about the experimental paradigm; i.e., onsets and durations of stimuli, is used to construct temporal regressors, which are then fitted to the time course of every voxel. The subsequent statistical hypothesis testing, for a given contrast of fitted weights, is a mass univariate approach leading to an “activation map” that highlights brain regions for which sufficient evidence is present to be related to the experimental paradigm. The relationship between the measured blood oxygen level dependent (BOLD) and the experimental paradigm can be modeled as a linear shift-invariant system with the hemodynamic response function (HRF) as an impulse response. In standard GLM approached (Friston et al., 1998), the HRF is predefined using gamma functions; more flexible techniques estimate HRF components in a subject- or time-dependent way in order to deal with inter- and intra-subject variability (Aguirre et al., 1998). More notably, parcel-based HRF estimation methods through joint detection estimation (JDE) framework are studied based on Bayesian approaches (Chaari et al., 2013, Makni et al., 2008, Vincent et al., 2010). Recently, an adaptive parcel identification driven by the hemodynamics is proposed using JDE (Chaari et al., 2012, Thirion et al., 2006). These HRF identification methods are mainly combined with GLM analysis to explore the parcel/subject/group/task specific hemodynamic models.

Not all brain activity can be modeled beforehand using stimulus functions; e.g., interictal epileptic discharges occur spontaneously. In addition, mere resting-state, which was neglected before as background noise, is known (Raichle, 2006) to produce characteristic patterns of brain activity referred to as resting-state networks. Such unpredictable activity cannot be inferred from traditional GLM analysis approaches (Gusnard and Raichle, 2001). Therefore, there is an increasing need for methodologies that enable the exploration of hemodynamic brain activity without predefined responses (Cole et al., 2010). Data-driven methods have been proposed for that purpose such as fuzzy clustering (Baumgartner et al., 2000), temporal clustering analysis (TCA) (Liu et al., 2000, Morgan et al., 2008), seed correlation analysis (Biswal et al., 1995), or subspace decomposition methods such as independent component analysis (ICA) (Beckmann and Smith, 2004, Calhoun and Adali, 2006), canonical correlation analysis (CCA) (Afshin-Pour et al., 2012) and agnostic canonical variates analysis (agnostic-CVA) (Evans et al., 2010). ICA is probably the most commonly used data-driven method. It provides a bilinear decomposition of the data into components that consist out of a spatial map with an associated time course. Its application to fMRI typically relies on (a surrogate for) spatial statistical independence between the components. However, such criterion is not directed specifically to identify “activation like” components since no knowledge is taken into account about the hemodynamics or about the type of activity-driven signal (e.g., spikes versus sustained activity).

fMRI deconvolution methods have been proposed to uncover the underlying activity-inducing signal at the fMRI timescale of seconds. Initially, Glover (1999) introduced Wiener deconvolution filtering that is optimal for Gaussian sources and thus results in very smooth activity-inducing signals. This work was generalized by Gitelman et al. (2003) to study the psychophysiologic interactions at the neuronal level. Following recent advances in convex optimization theory, these methods can be made more sophisticated by adding sparse priors on the underlying signal and solve maximum a posteriori estimation rather than naive maximum likelihood estimation. They have been exploited as an extension to standard GLM analysis by defining spatial priors on spatial activation maps (Flandin and Penny, 2007, Harrison et al., 2008, Smith and Fahrmeir, 2007, Vincent et al., 2010). Within a temporal fMRI deconvolution framework, while the linear system assumption on the hemodynamic model is retained, regularization terms use ℓ₁-norm to favor sparse solutions in time; i.e., a limited number of spike-like activations. For example, Hernandez-Garcia and Ulfarsson (2011) use the majorization–minimization scheme, Gaudes et al., 2011, Gaudes et al., 2013 rely on ridge-regression and sparsity-promoting estimators, and (Khalidov et al., 2007, Khalidov et al., 2011) impose sparsity in the activelet-domain, a wavelet basis that is tailored to the hemodynamic properties. These methods do not require any knowledge on the timing and exploit temporal properties of the HRF. Promising results have been demonstrated for local event detection (Gaudes et al., 2011), especially in epilepsy (Lopes et al., 2012), but also for resting-state analysis (Petridou et al., 2012). Recently, non-linear models have also attracted a lot of attention for blind deconvolution to explore the network dynamics. These methods solve non-linear state-space model in continuous time, which bring along high resolution solutions, and access both hidden states and (non)dynamic parameters related to the neuronal activity via Bayesian filtering, Cubature Kalman filtering and Local Linearization filters (Friston et al., 2008, Friston et al., 2010, Havlicek et al., 2011, Riera et al., 2004). These methods, however, have high computational cost compared with linear models and are mostly applied for uncovering the hemodynamics of a priori regions of interest.

Here we propose a novel deconvolution method for fMRI data analysis, for which we coin the term “total activation” (TA). TA includes some unique features to go beyond the limitations mentioned above:

• Express temporal properties of the activity-inducing signal. The deconvolution identifies the “innovation” signal (which is spike-type) as the sparse driver of the BOLD signal. However, the activity-inducing signal can be more flexible such as block-type signals (Karahanoglu et al., 2011).

• Structured sparsity for combined temporal and spatial regularization. Spatial regularization is incorporated using mixed-norms based on anatomical priors of brain regions (Baritaux et al., 2011); i.e., time courses of voxels in the same brain regions are favored to be coherent.

• Take advantage of efficient optimization schemes. We employ the efficient generalized forward–backward splitting algorithm (Raguet et al., 2012), which is a fast iterative shrinkage algorithm that alternates between temporal and spatial domain solutions until convergence to the final estimate of the underlying activity-inducing signal.

The paper is organized as follows. We first introduce the TA theoretical framework. Next, the feasibility of TA is demonstrated on both synthetic and experimental data. The simulation study allows validating the performance for block-type activity-inducing signals with different block durations. The experimental study is based on fMRI data acquired in three healthy subjects while being at rest but with several (unexpected) visual stimuli. The TA deconvolved signals indicate strong and short periods of activity in the primary visual regions that match with the experimental timing. The dynamic “activation maps” show coordinated activation (and “de-activation”) in large-scale networks, mostly with much larger average block lengths. We also show that resting-state networks can be retrieved using hierarchical clustering of the region-averaged activity-inducing signals.