An Anisotropic 4D Filtering Approach to Recover Brain Activation From Paradigm-Free Functional MRI Data

1. Introduction

Functional MRI (fMRI) is a non-invasive imaging technique that indirectly probes brain function by providing a measure of the metabolic activity consequent to an increased neural activation. Two experimental setups are commonly used to acquire fMRI data, task-fMRI, and resting-state fMRI (rs-fMRI). In the former, the subject is asked to follow an experimental paradigm, whereas in the latter the subject is asked to rest in the scanner and not do or think of anything in particular. Standard approaches for the analyses of task-fMRI data are based on the well-known general linear model (GLM) adapted by Friston and colleagues in 1998 in the context of fMRI data analyses (Friston et al., 1998). This approach requires prior knowledge of the task parameters and timing of events, as well as assumptions about the neural and hemodynamic responses. Therefore, the GLM can be used only for task-fMRI experiments, where the expected stimulus response is given by experimental paradigm. In contrast to task-fMRI, where the focus is on the response to a specific stimulus, rs-fMRI provides insight on brain function in the absence of stimuli. It also allows us to map the brain activity of patients whose condition does not allow them to perform tasks or follow an experimental paradigm. Furthermore, there are also brain activations that cannot be modeled and expected, such as seizures in epileptic patients (Karahanoğlu et al., 2013). For these data, which represent unpredictable brain activity, the GLM approach is not suitable (Gusnard and Raichle, 2001).

Data-driven methods have been proposed to analyze images obtained in resting-state, when no information about the occurrence of the activation is available. They include blind source separation approaches such as the independent component analysis (ICA) (McKeown et al., 1998; Beckmann and Smith, 2004; Calhoun and Adali, 2006), the principal component analysis (PCA) (Andersen et al., 1999; Baumgartner et al., 2000), the temporal clustering analysis (Liu et al., 2000; Morgan et al., 2008), and clustering methods (Cordes et al., 2002; Salvador et al., 2005; Golland et al., 2008; Lee et al., 2012). These methods are of interest if the aim is to group voxels showing the same spatial or temporal features but they cannot be used if the goal is identifying activations at the voxel level. Indeed, they do not consider including any hemodynamic effect. They are also limited by the necessity of choosing a priori the number of components or clusters and by their interpretation (Gaudes et al., 2011).

To overcome these limitations, deconvolution approaches have been developed to address the problem of studying and uncovering brain activations hidden within fMRI time series at the voxel level. fMRI deconvolution was introduced by Glover in Glover (1999), who investigated the performance of Wiener deconvolution for deblurring the fMRI response and reduce image distortions. This approach resulted in smooth recovered activation (Karahanoğlu et al., 2013) and required an independent measurement of the noise spectral density (Gitelman et al., 2003). Gitelman et al. (2003) developed an approach based on linear deconvolution and modeled the interplay between areas as effects of psycho-physiological interactions. In addition, dynamical filter methods, such as Kalman and Bayesian filtering, and local linearization filters have been developed and applied to fMRI (Riera et al., 2004; Friston et al., 2008; Havlicek et al., 2011). However, because these approaches are based on non-linear models in continuous time, they are limited by the high computational cost and convenient only for the analysis of localized regions of interests (ROIs) (Karahanoğlu et al., 2013). Other approaches make spatial and/or temporal assumptions on the underlying signals, thus adding priors in the optimization problems. In particular, sparse regularization on the recovered activation maps was exploited to force to zero the weights of regressors, which did not contribute to the activation (Flandin and Penny, 2007; Smith and Fahrmeir, 2007; Harrison et al., 2008). L1-norm regularization approaches have also been developed to exploit sparse temporal features of the hidden neural activation. This was done by means of the majorization–minimization of a cost function to find an optimal solution to the inverse problem (Hernandez-Garcia and Ulfarsson, 2011). Caballero Gaudes et al. developed a ridge-regression regularization by minimizing both the variance of the residuals and the power of the resulting estimate of the input signal representing the brain activity (Gaudes et al., 2011) and a sparse regression (Caballero Gaudes et al., 2013) by assuming short neuronal activations. Temporal regularized optimization problems based on wavelets were also explored (Khalidov et al., 2011). These methods exploit the temporal features of the hemodynamic response function (HRF) (Karahanoğlu et al., 2013) attempting to deconvolve it from the fMRI data, rather than using any information on the timing of the neural activations which translate on the blood oxygen level dependent (BOLD) response. Recently, by supposing the brain activates in constant blocks, Karahanoğlu et al. (2013), later revisited by Farouj et al. (2017), developed a deconvolution approach which involves both spatial and temporal regularization called total activation (TA). These approaches split the optimization problem into two decoupled spatial and temporal regularizations that increases the number of parameters to be set and requires the solver to alternate between the constraints. The temporal regularization of the latter was then improved by Costantini et al. (2018) by proposing a joint approach using the least angle regression (LARS) algorithm and the L-curve, which overcame the limitation of having to choose a priori the regularization parameter, and reduced the computation time. Recently, deconvolution algorithms have also incorporated a multivariate formulation to perform spatiotemporal deconvolution in 2D (Uruńuela et al., 2020, 2021). The first allows the estimation of blocks of activations coupled with a subsampling approach based on stability selection to avoid the choice of the regularization parameter. The latter enables a method that finds global fluctuations due to motion artifacts or physiological signals as well as neuronal activity. Also, Bolton et al. (2019) proposed a spatiotemporal deconvolution approach that includes structurally informed regularization.

Our approach is based on the idea that brain activations occur in spatially and temporally coherent patterns with sharp boundaries. These patterns are projected to fMRI data via the HRF and corrupted by physiological and motion artifacts. To recover the underlying brain activations, we propose a novel method based on partial differential equations (PDEs), named anisotropic 4D filtering fMRI (A4D-fMRI). The A4D-fMRI applies an anisotropic diffusion process whose diffusivity is steered by derivatives of the evolving image to smooth the fMRI image and simultaneously enhance important features such as spatial edges and temporal functional activations. Regularization methods have been enriched by the use of non-linear PDEs in several contexts for the last 30 years. First applied to physics and fluid mechanics, it has been shown that non-linear PDEs allow smoothing the data while preserving large global features, such as discontinuities of the signal (Tschumperle and Deriche, 2005), which can be found, for example, in image contours and corners (Tschumperlé and Deriche, 2007). This approach is based on the isotropic diffusion equation, i.e., heat flow, and has subsequently been extended to other theoretical contributions. Among them there are the anisotropic smoothing (Weickert, 1998; Sapiro, 2006) and the PDEs-based gradient descent used to solve energy functionals minimizations (Rudin et al., 1992; Chambolle and Lions, 1997; Charbonnier et al., 1997; Kimmel et al., 2000; Aubert and Kornprobst, 2006). The pioneering work that employed anisotropic diffusion PDEs for the restoration of noisy and blurred digital data was proposed by Perona and Malik (1990) overcoming the limitations associated with linear filtering approaches (Tschumperle and Deriche, 2002). To date, PDEs-based regularization algorithm has been applied to 2-D scalar images (Perona and Malik, 1990; Nielsen et al., 1997; Weickert, 1998; Aubert and Kornprobst, 2006) and vector-valued images (Tschumperle and Deriche, 2002).

The A4D-fMRI proposed in this work has been conceived for the geometrical regularization of 4D fMRI images (3D space × 1D time) based on PDEs. The A4D-fMRI acts concurrently in space and time, thus overcoming the limitation of previous deconvolution approaches that consider the two problems of spatial and temporal regularization as decoupled processes. In our method, the regularization flow is performed according to the time and to the local geometry of the image to evaluate the presence of an edge and its local strength.

The rest of the article is organized as follows. We first introduce the PDEs regularization theoretical framework, illustrate the mathematical problem, and how we solved it. Next, we validate the A4D-fMRI on phantom data and on task-fMRI data from 51 subjects using the experimental paradigm as ground truth and positively compare its performance to the TA approach. We also applied the A4D-fMRI approach to rs-fMRI data. Finally, the advantage of our approach of not requiring any knowledge of the experimental paradigm is shown by recovering activity that agrees with the GLM, which uses and need this a priori knowledge.

2. Theory

The fMRI BOLD signal is a combination of brain activation, the HRF, physiological artifacts, and noise. More specifically, fMRI data are modeled as a brain activation that is convolved with the HRF and corrupted by noise. So, starting from the corrupted fMRI images, the purpose of this work is to remove from this signal the noise and the hemodynamic effect and to keep only the underlying activation signal that is at its origin. To achieve this goal, in this paper we propose to solve this regularization problem using diffusion theory. Inspired by the physics of fluids, many authors assimilated the process of image regularization with the diffusion of chemical concentrations and proposed to apply the following diffusion PDE process (Weickert, 1998, 1999; Tschumperle and Deriche, 2002, 2005; Tschumperlé and Deriche, 2007),

$\begin{array}{ll}\frac{\partial \text{I}}{\partial \text{t}}=\text{div}(\text{D}\nabla \text{I})& (1)\end{array}$

where I is the input image, ∇ is the gradient operator, t is the time, div(·) is the divergence operator, and

$\begin{array}{ll}\text{D}={\lambda }_1\text{u}{\text{u}}^T+{\lambda }_2\text{v}{\text{v}}^{\text{T}}& (2)\end{array}$

is the diffusion tensor of the image I, also called structure tensor (Förstner, 1986; Förstner and Gülch, 1987; Tschumperle and Deriche, 2005). The diffusion tensor D is a symmetric and positive definite matrix, and has λ₁, λ₂ as positive eigenvalues and u, v as corresponding orthogonal eigenvectors that drive the regularization process; the amount of diffusion in the directions u and v will be weighted by λ₁ and λ₂, respectively. PDEs smooth the image at each step with a notion of scale-space (Perona and Malik, 1990; Nielsen et al., 1997; Lindeberg, 2013). The scale-space technique involves generating coarser resolution images by convolving the original image with a Gaussian kernel. Therefore, starting from an original image, according to one parameter, i.e., the size of the smoothing kernel, a set of smoothed images are produced. At each iteration, the image is smoothed and fine-scale properties, such as noise in our case, are gradually suppressed. To apply the diffusion process to fMRI data, Equation (1) will need to be reformulated to handle 4D images. Let us define a scalar-valued image as a function I:Ω⊂ℝ⁴, where Ω is the domain of the 4D (3D space × 1D time) image and let us assume Neumann boundary conditions on δΩ, specifying the values in which the derivative of the solution is applied within the boundary of the domain. Let us now define a structure tensor D as a 4 × 4 symmetric and positive-definite matrix. By definition, D has four positive eigenvalues (λ₁≥λ₂≥λ₃≥λ₄≥0) and their associated four orthogonal eigenvectors (θ₁, θ₂, θ₃ and θ₄) explain the distribution and orientation of the gradient ∇I = (I_x, I_y, I_z, I_t) of the image I in a given neighborhood. A structure tensor can distinguish between anisotropic and isotropic diffusion. If λ₁>>λ₂, λ₃ and λ₄, the structure tensor has a principal orientation (in this case θ₁) and the diffusion is anisotropic. It can be represented with an ellipsoid oriented along θ₁. On the other hand, if λ₁≈λ₂≈λ₃≈λ₄, the structure tensor is not oriented in a main direction and θ₁, θ₂, θ₃, and θ₄ are eigenvectors of D with equal weight. In this situation, the diffusion is isotropic and the structure tensor can be represented with a sphere.

Inspired by the physical process of diffusion, we link the diffusion to fMRI image regularization and we propose the A4D-fMRI for the enhancement of coherent structures found in fMRI data. The A4D-fMRI recovers brain activations and smooths small variations while preserving large variations via a regularization that is applied on the 4D fMRI image, acting simultaneously in the 3D space and the 1D time dimensions. Using this approach corresponds to performing minimization of image variations as well as a blind image deconvolution. To do this, we propose a regularization process based on a gradient descent computed with PDEs, such that

$\begin{array}{ll}\frac{\partial \text{I}}{\partial t}=(1-\alpha )\frac{{\text{H}}^T({\text{I}}_0-\text{HI})}{\Vert {\text{I}}_0{\Vert }_2}+\alpha \frac{\text{div}(\tilde{D}\nabla \text{I})}{\Vert \text{div}({\tilde{D}}_0\nabla {\text{I}}_0){\Vert }_2}& (3)\end{array}$

where the term on the left, ∂I/∂t, is the regularization flow, the first term on the right is the data fitting term, also called fidelity term, that prevents the solution from straying far from the input data, and the second term on the right minimizes image variations. The parameter α∈[0, 1] is the user-defined regularization parameter that balances the fidelity and the regularization terms. Starting from the initial image I₀, the restored image I is regularized as t increases reducing noise and extracting coherent space and time variations. At the same time, no new structures are introduced in the image (Aubert and Kornprobst, 2006).

Going more into the details of the fidelity term

$\begin{array}{ll}F(\text{I})=\frac{{\text{H}}^T({\text{I}}_0-\text{HI})}{\Vert {\text{I}}_0{\Vert }_2}& (4)\end{array}$

where I₀ and I are the original and the regularized image, respectively, ||I₀||₂ in the denominator is the normalization factor, H is the HRF (Khalidov et al., 2011) operator, and H^T is its transpose. The multiplication of H^T with (I₀−HI) corresponds to a correlation and can be implemented via convolution with the time-reversed HRF. Note that this product is computed only along the time dimension. The HRF considered in this work is the linearized time-HRF operator proposed by Friston et al. (2000) and Khalidov et al. (2011).

The regularization term is defined as

$\begin{array}{ll}R(\text{I})=\frac{\text{div}(\tilde{D}\nabla \text{I})}{\Vert \text{div}({\tilde{D}}_0\nabla {\text{I}}_0){\Vert }_2}& (5)\end{array}$

where $\Vert \text{div}({\tilde{D}}_0\nabla {\text{I}}_0){\Vert }_2$ is the normalization term and $\tilde{D}$ is the regularization tensor, distinct from the structure tensor D. In order to elucidate the regularization term, let us start by the definition of the operator

$\begin{array}{ll}\text{D}=\frac{\nabla I\nabla I^T}{\Vert \nabla I{\Vert }^2}*G& (6)\end{array}$

that is the 4D structure tensor of I smoothed by the Gaussian kernel G with standard deviation σ_G via the convolution operator *. The matrix D being the diffusion tensor of the image I, its eigendecomposition gives a set of eigenvalues and eigenvectors such that, if the gradient in one direction is large, the eigenvalue associated with that direction is large, whereas the eigenvalues associated with the other three directions are relatively small. Since we are processing fMRI images with the aim of saving activations and contours that occur concomitant to a large gradient in a certain direction, we aim at reversing the diffusion process, therefore at reversing the effect of D into $\tilde{\text{D}}$ to enhance and simultaneously simplify coherent structures of the fMRI image. Here, we propose to compute the operator $\tilde{D}$ by modifying the eigenvalues of the operator D in Equation (6). Specifically, we defined the directions of the image variations by an eigendecomposition of D such that

$\begin{array}{ll}\text{D}=\text{Q}\Lambda {\text{Q}}^T& (7)\end{array}$

where Q contains the orthogonal eigenvectors (θ₁, θ₂, θ₃, θ₄) of D and Λ contains their associated eigenvalues (λ₁≥λ₂≥λ₃≥λ₄). We then recomputed the matrix

$\begin{array}{ll}\tilde{D}=\text{Q}\tilde{\Lambda }{\text{Q}}^T& (8)\end{array}$

where $\tilde{\Lambda }$ is a diagonal matrix with entries ${\tilde{\lambda }}_1$, ${\tilde{\lambda }}_2$, ${\tilde{\lambda }}_3$, and ${\tilde{\lambda }}_4$ such that for each voxel the highest eigenvalue is given by

$\begin{array}{ll}{\tilde{\lambda }}_1=\text{exp}(-\frac{{\lambda }_1^2}{\text{max}{(\Lambda )}^2}\frac{1}{2{\sigma }_D^2})& (9)\end{array}$

and the other eigenvalues are ${\tilde{\lambda }}_2={\tilde{\lambda }}_3={\tilde{\lambda }}_4=1$. The rationale is that, if λ₁ is large, the current voxel may be located on a edge or activation and the diffusion tensor $\tilde{D}$ is steered to be anisotropic, by setting ${\tilde{\lambda }}_1\ll {\tilde{\lambda }}_2,{\tilde{\lambda }}_3,{\tilde{\lambda }}_4$. Since we aim at performing a smoothing only along the other three directions to smooth preferably along the coherence directions, the three eigenvalues λ₂≈λ₃≈λ₄ are set to 1. On the other hand, if λ₁ is small, the diffusion will be isotropic in the four directions because ${\tilde{\lambda }}_1\approx 1$ and ${\tilde{\lambda }}_2={\tilde{\lambda }}_3={\tilde{\lambda }}_4=1$. Using the function in Equation (9) corresponds to reassigning to each voxel different eigenvalues constituting the matrix $\tilde{\Lambda }$, before recomputing the operator $\tilde{D}$ as in Equation (8). In fact, if λ₁/max(Λ) is large, the highest eigenvalue $\widetilde{{\lambda }_1}$ of the considered voxel will tend to zero. This steers the geometrical regularization to be anisotropic, because the smoothing will apply equally in the remaining three directions but it will be negligible in the normal to the detected contour. Otherwise, if λ₁/max(Λ) is small, the greatest eigenvalue $\widetilde{{\lambda }_1}$ will tend to 1 and this leads to an isotropic regularization almost in all the four directions (x, y, z, t). In both cases, $\widetilde{{\lambda }_2}$, $\widetilde{{\lambda }_3}$, $\widetilde{{\lambda }_4}$ are set to 1. This procedure is applied to each voxel of the entire 4D image such that at each iteration the image I computed in Equation (3) is removed from the image at the previous iteration.

In practice, the gradient and divergence operator are implemented using finite differences. The gradient was implement using a left to right scheme and the divergence using right to left, leading to a discrete Laplacian if the diffusion tensor was omitted. To maintain stability across iterations used to numerically solve the PDE, the step size must be small, but reducing it increases computation times. We empirically set it to 0.1, which provided stability and acceptable computation times.

In this way, supposing the brain activates in constant blocks, we regularized the image together in space and time. We were able to keep large image variations occurring during brain activations or spatial edges, and to gradually remove small variations, corresponding to noise, while conserving and enhancing coherent structures of the fMRI image. The theoretical framework of the A4D-fMRI explained above has been implemented in a Python package that can be downloaded at the following link: gitlab.inria.fr/cobcom/a4dfmri. The validation illustrated in the following sections is performed using the A4D-fMRI package.¹

3. Methods

3.1. Simulation of fMRI Data

To reproduce the acquired fMRI signals, we started by simulating the fMRI activation as a piece-wise constant function as proposed by Farouj et al. (2017). For each voxel v, we modeled the activity-inducing signal as a boxcar function

$\begin{array}{l}u(t)=H(t-a)-H(t-b)\end{array}$

where H(t) is the Heaviside step function. We added noise to u(t) representing the random intrinsic electrical fluctuations within neuronal networks that are not associated with encoding a response to internal or external stimuli. To do this, we corrupted the activity-inducing signal u(t) with an additive random Gaussian noise with zero mean and standard deviation σ_m that we called “model noise” ϵ_m. The noisy activity-inducing signal is

$\begin{array}{l}{u}_n(v,t)=u(t)+{ϵ}_m.\end{array}$

We modeled the activity-related signal x(t), consequent to the neural activation as the convolution of u_n(t) with the HRF, modeled as the linear time-invariant system h(t) (Khalidov et al., 2011):

$\begin{array}{l}x(v,t)=u_n(v,t)*h(t).\end{array}$

Real time series acquired using the fMRI technique is corrupted by different kinds of noise and artifacts given by mechanisms that do not reflect any neurophysiological function, such as heart rate, respiratory fluctuations, motion artifacts, thermal noise, and scanner drifts (Lund et al., 2006). For this reason, we added noise to x(t) thus obtaining the acquired fMRI signals

$\begin{array}{l}y(v,t)=x(v,t)+{ϵ}_a=u_n(v,t)*h(t)+{ϵ}_a\end{array}$

where ϵ_a is the additive random Gaussian noise with zero mean and standard deviation σ_a. A scheme representing the model of the phantom fMRI data is shown in Figure 1.

FIGURE 1

Figure 1. fMRI model. From left to right, u(t) is the activity inducing signal that represents the neural activation as a piece-wise constant signal. To this signal, a model noise ϵ_m with a Gaussian distribution was added, thus leading to the signal u_n(t) that was then convolved with the hemodynamic response function operator ${H}$, thus obtaining x(t), also called activity-related signal. Adding the noise ϵ_a to the signal x(t) gives the simulated acquired fMRI data denoted as y(t).

3.2. Validation of Synthetic Data

To test and validate the A4D-fMRI, we scaled a 3D activation map computed with the FMRIB Software Library (FSL²) Physics-Oriented Simulated Scanner for Understanding MRI (POSSUM³) in the range [0, 3], with a 2-mm isotropic resolution (Figure 2A). We multiplied it by a piece-wise constant signal u(t) of 100 s, with one onset of 40 s, from 20 to 60 s (Figure 2B). Starting from u(t), we simulated the acquired fMRI time-courses y(t) as explained in the previous section. We tested the A4D-fMRI on several simulated images obtained by adding different amount of noise for each experiment. We regularized the whole image using the A4D-fMRI as shown in Section 2, and we recovered the voxel-wise activity-inducing signals û(t). To evaluate the results, we compared the simulated activity-inducing signal u(t) with the recovered û(t). To do this, for each voxel's time course belonging to the GM, we first computed the mean square errors (MSEs) between û(t) and u(t). Second, we computed the roots of the mean (RMSE) and the standard deviation (STD) of this array set. To be able to compare these values to the TA results, which restrict computations to the Gray Matter (GM) voxels, we only considered values obtained among the voxels belonging to the GM mask. Similarly, the Pearson correlation (r) was computed for each GM voxel between the simulated activity-inducing signal u(t) and the recovered û(t). The mean and the STD was than computed on this set of r values. We compared our results with those obtained using the TA approach (Farouj et al., 2017), implemented in the TA toolbox.⁴ Note that because of the L1-norm regularization used by TA which tends to underestimate signal amplitudes, comparing the results based only on amplitudes is not fair (see Supplementary Materials for additional discussion). For this reason also, we compared the performance of both algorithms based on the correlation between the simulated and the recovered activation. Pearson's correlation coefficient is invariant to the scale of the input signals and therefore ignores the signal shrinkage associated with the L1-norm and quantifies only the shape of reconstruction. We did not explicitly set any regularization parameters for the TA approach and instead relied on the automatic selection heuristics provided by the implementation of the TA toolbox⁴ (Karahanoğlu et al., 2013; Farouj et al., 2017). The results presented were obtained by providing only our data with its specific repetition time (TR). The temporal regularization parameter is initialized in the TA toolbox based on a pre-estimated noise level value of the data fit, derived from the median absolute deviation of fine-scale wavelet coefficients, then updated at each iteration. The spatial temporal regularization parameter is already preset in the TA toolbox for the gray matter constrained total variation algorithm and therefore follows the recommendation of the authors. To investigate the performance of our algorithm to event-related designs, we also simulated brain activations as spikes and applied our approach on the noisy synthetic fMRI images.

FIGURE 2

Figure 2. Spatial map (A) and time series (B) of the activation considered as ground truth for functional MRI simulated data. The time course of the activation u(t) was simulated with a repetition time of 1 s.

3.3. Validation of Real Experimental Data

3.3.1. Task-fMRI Data

In addition to validating the A4D-fMRI on simulated data, we evaluated its performance on real data. Although our approach was conceived to be applied to rs-fMRI data, or in any situation where no experimental paradigm is given, we test it here on task-fMRI data. This testing strategy provides a ground truth, i.e., the task timing, which can be used to assess the performance of the algorithm.

The study was conducted on the motor task-fMRI data from 51 subjects from the Human Connectome Project (HCP) database (Van Essen et al., 2013). The data underwent a minimal pre-processing pipeline (Glasser et al., 2013), which includes correction of gradient-non linearity-induced distortions, registration of each image frame to the signal-band reference image to achieve motion correction, phase-encoding distortion correction, EPI image distortion correction, registration of the fMRI volumes to the structural data, coregistration of the fMRI data to the Montreal Neurological Institute (MNI) space, masking and fMRI image intensity normalization to the 4D whole global mean of 10³. As additional pre-processing steps, each voxels' time course was detrended to remove linear trends and normalized to 0 mean and unit standard deviation. The motor task is initiated by a visual cue followed by the movement of the left and right foot, the left and right hand, and the tongue. The tasks starting points were considered equal for each subject and inter-subjects differences of the order of milliseconds were neglected.

After applying the A4D-fMRI on the entire brain images of each subject, we recovered the reconstructed activity-inducing signals û(t) without prior knowledge on the onset/offset times and location of the evoked stimuli. The regularization parameter α was set experimentally to 0.9997, σ_G was set to 1, σ_D was set to 0.2, and we computed up to 40 iterations. The value of α is driven by two factors: the divergence of the expected solution and the amount of noise in the data. One is very small converging toward zero and a second one stabilizing at the variance of the noise in the data. This stopping criteria was chosen because it was the minimal number of iterations required for the Pearson correlation to stabilize for all tasks fMRI datasets.

To highlight the ability of the A4D-fMRI to recover brain activations without knowledge of the experimental paradigm, we qualitatively compared brain regions recovered using the A4D-fMRI to those recovered using the GLM as implemented in the FSL library. The GLM model requires as input the exact occurrence of the tasks that the subject is asked to perform. We run the A4D-fMRI on the whole brain thus blindly recovering brain activations, without prior knowledge on the intervals of the evoked stimuli. To estimate the results obtained using the A4D-fMRI, we computed the voxel-wise correlation maps, by estimating the Pearson correlation coefficient (r) between the recovered activations and the five tasks. The tasks were simulated as a piece-wise constant signal with unit amplitude when the subject is performing the task and zeros elsewhere. As for the GLM, we included the five tasks in a design matrix and we estimated the regressors' weights with FSL. Results showing differences and similarities of both approaches were qualitatively assessed.

Subsequently, we quantitatively compared the results obtained using the A4D-fMRI with the ones given by the TA, on the sample data composed by 51 subjects. We first defined four ROIs located in brain regions that are involved in the five considered tasks. The ROI related to the tongue was bilateral, whereas for the hands and the feet we defined separate ROIs for the left and the right side of the brain. To do this, we started by defining the ROIs from the work proposed by Roux et al. (2018), who mapped the somatosensory homunculus MNI coordinates using the electrostimulation. For each coordinate center, we built a spherical 3-mm-radius ROI and we grouped the multiple ROIs related to each task into a unique ROI. Coordinates' centers are shown in Table 1.

TABLE 1

Table 1. Montreal Neurological Institute (MNI) coordinates centers of the brain areas found to respond to the somatosensory stimulation.

After defining the ROIs, similarly to the comparison between the A4D-fMRI and the GLM, we computed the whole-brain voxel-wise correlation maps between the time course related to each task and the recovered activity-inducing signals û(t) obtained with the A4D-fMRI and the TA. For each subject, we first computed the average of the Pearson correlation coefficients (r) inside the GM-masked ROIs and then calculated the mean and the standard deviation of these averaged correlation values across the 51 subjects belonging to the sample data.

Furthermore, to show that the A4D-fMRI is able to differentiate between a region that is active and one that is not, the time courses û(t) of one representative subject (100307) were averaged in two ROIs of 6 × 6 × 6mm³: one that is expected to be active during the task, and one located in a brain area which is not involved in the task. We selected the task related to the tongue, and we chose one ROI centered in the Brodmann Area 4p (rBA4p; MNI coordinates: 62, −14, 30) that is activated during a tongue motor task, and another centered in the primary auditory cortex (TE1.2; MNI coordinates: 56, 4, 10; Kiviniemi et al., 2003), that is not involved in the tongue movement. The Pearson correlation (r) was computed between the tongue activation and the recovered û(t) for each voxel, and then averaged among the voxels belonging to the two GM-masked ROIs. Again the tongue activation was simulated as a piece-wise constant signal with unit amplitude during the task and zeros elsewhere. We compared results obtained using the A4D-fMRI with the ones obtained using the TA toolbox.

3.3.2. Resting-State fMRI Data

Finally, we applied the A4D-fMRI on the rs-fMRI image of one subject (100307) from the HCP database. The data were acquired with a SIEMENS MAGNETOM Connectome Syngo MR D11 using a gradient-echo EPI sequence (TR = 720 ms; TE = 33.1 ms; flip angle = 52°; FOV = 208 × 180 mm; slice thickness 2.0 mm; number of slices = 72; 2.0 mm isotropic voxels; multiband factor = 8). The subject was asked to lay in the scanner without thinking about anything in particular. The number of acquired frames was 1,200 and the duration of the acquisition was 14:33 min. In the case of rs-fMRI data, the task paradigm is unavailable since the subject does not perform any task in the scanner. The data underwent the same minimal preprocessing of the task fMRI data as proposed in the HCP pipeline (Glasser et al., 2013). In addition, the time series were detrended to remove linear drifts and normalized to zero mean and unit standard deviation. We applied the A4D-fMRI algorithm on the entire rs-fMRI sample and we observed the dynamics of the recovered activation maps across time.

4. Results

4.1. Validation on Synthetic Data

Figure 3 shows examples of regularized spatial maps (Figure 3A) and time series (Figure 3B) using the A4D-fMRI (û_A4D−fMRI) and the TA (û_TA). Both approaches do not require any prior knowledge of the paradigm timing. The regularized spatial maps in Figure 3A represented in the axial plane show how the regularized fMRI image recovered with the A4D-fMRI (û_A4D−fMRI) is closer to the ground truth (u) in terms of signal amplitude with respect to those obtained using the TA û_TA. This is verified for different peak-SNRs (pSNRs), i.e., 6.54, 5.99, and 3.93 dB. In Figure 3B, we show examples of regularized time courses and that we recover an amplitude closer to the ground truth when compared to the method implemented in the TA tool. We also show smoother recovered signals when compared to the TA. Additional tests and results are shown in Supplementary Figure 1.

FIGURE 3

Figure 3. (A) From left to right, spatial maps of the simulated fMRI image y, ground truth activation u, recovered activation using the total activation (TA) approach (û_TA) and our approach (û_A4D−fMRI). Each row corresponds to a different peak-SNR (pSNR): 6.54, 5.99, and 3.93 dB from the top to the bottom. (B) Reconstructed time series û(t) obtained with our approach [û_A4D−fMRI(t), light-blue] and the TA approach [û_TA(t), magenta] superimposed on the activation [u(t), black] and functional MRI (fMRI) signal [y(t), orange].

Figure 4 shows that the roots of MSEs ± STDs computed between the simulated activation u(t) and the recovered one û(t) change for different pSNRs. We show lower errors with lower standard deviations than the ones obtained using TA.

FIGURE 4

Figure 4. On the left, the graph shows, for different pSNRs, the mean Pearson correlation coefficients (μ_r) and related standard deviation (σ_r) computed between u(t) and û(t) and averaged among the voxels belonging to the GM. On the right, The graph shows, for different pSNRs, the roots of the mean squared errors (MSEs) and standard deviations (STDs) between u(t) and û(t) averaged among the voxels belonging to the GM.

Figure 4 also shows that the activation recovered with the A4D-fMRI is more correlated with the ground truth (r≈1) for different pSNRs. Although the results obtained with TA are more sensitive to noise and show better performances for less noisy data: the mean correlation increases according to the pSNR while the standard deviation decreases. Results related to event-related designs, where activations were simulated as spikes, are illustrated in Figure 5. The results show that we are able to recognize a spike activation and remove the noise (see Figure 5 between 30 and 50 time points). However, the obtained activity-inducing signals are smoothed and do not match the exact lengths of the simulated spikes. We refer the reader to Supplementary Figures 2, 3 to see the results for different spike trains and model and additive amount of noise.

FIGURE 5

Figure 5. Reconstructed time series û(t) obtained with the anisotropic 4D-fMRI (A4D-fMRI) [û_A4DfMRI(t), light-blue] superimposed on the spike activation [u(t), black] and functional MRI (fMRI) signal [y(t), orange]. Peak-SNR = 13.23 dB.

4.2. Validation on Real Data

4.2.1. Task-fMRI Data

Figure 6 shows the structure tensors $\tilde{D}$ (Equation 8) in a coronal brain slice for a representative subject of the HCP database. Note that the image refers to the spatial maps, and the temporal dimension is not represented. The presence of ellipsoids oriented in different directions rather than spheres shows the anisotropic nature of the regularization.

FIGURE 6

Figure 6. Representation of the structure tensors computed with the anisotropic 4D-fMRI (A4D-fMRI). The images show for a 3D spatial map, superimposed to the standard MNI template, how the structure tensors look like ellipsoids or spheres meaning that the geometric regularization is applied anisotropically. Images were made using MRview (Tournier et al., 2012).

As for the real data analyses, and specifically the comparison between the A4D-fMRI and the GLM, we show that the correlation maps related to each task computed with the A4D-fMRI were well overlapped to the values of the regressors coefficients obtained using the GLM as shown for one illustrative subject (100307) in Figure 7. The GLM shows results that follow the GM, while the activations found with the A4D-fMRI, which again were performed across the whole brain, and not masked with the GM mask, cover also voxels across the white matter. Interestingly, the found activations overlap the areas found to be active in the motor homunculus brain (Penfield and Boldrey, 1937).

FIGURE 7

Figure 7. Qualitative comparison between the general linear model (GLM) and the anisotropic 4D-fMRI (A4D-fMRI). On the left column, in a blue-lightblue color-map, superimposed to the standard MNI template, the β-regressors map obtained using the GLM implemented in FMRIB Software Library (FSL). On the right column, in a red-yellow color map, the whole-brain voxel-wise correlation maps obtained using the A4D-fMRI superimposed to the standard MNI brain. The Pearson correlation (r) was computed voxel-wise across the whole brain, between the reconstructed activity inducing signals û(t) and the five motor tasks simulated as piece-wise constant signals with ones in the time points where the subject is executing the task and zeros elsewhere. The values r of the correlations are indicated by the color bars. Each row corresponds to a specific motor task, from the top to the bottom: the tongue, the right and left hand, and the right and left foot. A, anterior; P, posterior; S, superior; I, inferior; R, right; L, left. The unthresholded images are available in Supplementary Materials.

Quantitative comparison between the activity-inducing signal recovered using the A4D-fMRI and the TA is shown in Figure 8. Results show that the mean Pearson correlation values estimated for each ROI across the data sample increase while increasing the number of iterations, until it converges after 25 iterations for the hands, 5 iterations for the feet, and about 35 iterations for the tongue. Moreover, starting from the first iteration, we show statistically significant higher correlation values compared to the ones obtained using the TA. In particular for the comparison between the A4D-fMRI and the TA, Figure 9 shows the reconstructed signals û(t) (Figure 9A) and the correlations values (Figure 9B) for a single subject (100307). We show a clear difference between the correlation values estimated in the area involved in the task and the one that is not involved. In fact, we show higher correlation between the tongue activation and the recovered activation û(t) in the ROI rBA4p, which was expected to be involved in the motor task, while a low correlation is shown in the ROI rTE12 that instead is not involved. The TA was not able to clearly distinguish between an active and an inactive region since it showed low correlation values for both ROIs. Figure 9 clearly show another advantage of using the A4D-fMRI over TA in the fact that the recovered amplitude is higher for A4D-fMRI than TA. If a single threshold is applied on the output of the A4D-fMRI to define either the region is active or not, we are clearly able to make this distinction using the A4D-fMRI. This is not the case if a threshold is applied on the TA output, because for a low threshold both regions would be considered active, whereas for a high threshold both would result inactive. We compared the performance of the implementations of TA and A4D-fMRI by running both algorithms on the same hardware. We used the fMRI motor task data analysis with a size of 109 × 91 × 109 × 284 where 284 is the number of volumes and 109 × 91 × 109 is the MNI space dimension. TA, as implemented in the TA toolbox, and constrained in the GM voxels, ran in for approximately 9–9.5 h while our algorithm processed the whole dataset in 7 min per iteration, with a total computation time of 4.5–5 h for 40 iterations.

FIGURE 8

Figure 8. Barplots of the mean (μ_r) ± standard deviations (σ_r) of the Pearson correlation coefficients (r) computed on the sample data (51 subjects) in five regions of interests (ROIs) related to the tasks of the left and right hand, the tongue, and the left and right foot. For each task, the bars in light-blue represent the results using the anisotropic 4D filtering fMRI (A4D-fMRI) for an increasing number of iterations (from 1 to 40). The bars in magenta represent the results using the total activation (TA) toolbox. (lHAND, left hand; rHAND, right hand; lFOOT, left foot; rFOOT, right foot).

FIGURE 9

Figure 9. (A) Reconstructed signals û(t) obtained with the anisotropic 4D filtering fMRI (A4D-fMRI) (light-blue) and the total activation (TA) tool (magenta) superimposed on the real acquired functional MRI (fMRI) signals (orange) and the simulated tongue activation (black). The plot on the left is related to the regions of interest (ROI) located on the Brodmann Area 4p (rBA4p), the plot on the right is associated to the ROI positioned on the primary auditory cortex (rTE1.2). All the signals were averaged across the voxels belonging to the GM-masked ROIs. The gray areas represent the occurrence and the duration of the tongue movements. (B) Mean Pearson correlation coefficients (μ) and their associated standard deviations (σ) computed between the tongue activation and the recovered signals û(t) averaged across the voxels belonging to the GM-masked ROIs (rBA4p on the left, rTE1.2 on the right).

4.2.2. Resting-State fMRI Data

As for the resting-state data analysis, in Figure 10 we show few example spatial maps extracted from the 1,200 analyzed time points. The figure shows a dynamic between the regions composing the default mode network, i.e., the posterior cingulate cortex, the medial prefrontal cortex, the lateral parietal lobules, and the temporal cortex (Buckner et al., 2008). As shown in the figure, these areas are not all active simultaneously, but they assemble and disassemble over time in different combinations.

FIGURE 10

Figure 10. Spatial maps obtained with the anisotropic 4D filtering fMRI (A4D-fMRI) spanned across 150 repetition times (TRs) (from 200 to 350 TRs). The maps show the dynamics between the areas composing the default mode network. These regions are not all active at the same time but dynamically assemble and disassemble over time. The maps show how the lateral parietal lobules are not always active simultaneously as well as the medial prefrontal cortex that can activate selectively with the posterior cingulate cortex and the right temporal lobe.

5. Discussion

In this paper, we have described an innovative method to analyze fMRI images and recover the location and the occurrence in time of the functional neural activations. The proposed approach achieves this blindly, without the necessity of a priori knowledge of timing, duration, and location of the underlying activations. The approach we proposed, namely the A4D-fMRI, geometrically regularizes the fMRI image such that it saves and highlights large image variations as they are present at the occurrence of a brain activation or in the presence of a spatial edge with respect to small image variations that instead are removed to reduce noise. To do this, we used the PDEs in a iterative algorithm and exploited the 4-D image structure tensor that defines the directions of the gradient in the neighborhood of a voxel and directs toward an anisotropic or isotropic regularization. This gradient contains all the four principal directions of the fMRI image that is composed of a 3D spatial image repeatedly acquired in time, which suggests that the whole 4D fMRI image was smoothed contemporaneously in space and time at once.

Other approaches have been proposed to analyze fMRI data. Among them are the (i) the GLM that fits a linear model to the fMRI time series, but it assumes prior knowledge of the tasks (Friston et al., 1998); (ii) deconvolution methods, which are used to uncover brain activations from the BOLD response without prior information on the underlying activity (Gaudes et al., 2011; Caballero Gaudes et al., 2013; Karahanoğlu et al., 2013; Farouj et al., 2017). The deconvolution approach in Karahanoğlu et al. (2013) splits the problem into a spatial and temporal regularization problems, meaning that the user has to specify two regularization parameters and the two weights used to have a solution that is given by a weighted sum of the two separate regularization processes. In contrast, the A4D-fMRI overcomes several limitations found in the previous literature. When comparing the regions recovered using the GLM and A4D-fMRI, we noted overall very good agreement between the two methods. Moreover, the activations found using our approach are bigger than those found using the GLM. Both approaches found active voxels outside the GM. The A4D-fMRI was run on the entire brain volume and did not rely on a GM mask. To reduce the size of the recovered action, a GM mask can be applied to the recovered brain activation. However, in this paper, whose aim was presenting the novel algorithm, we omitted this step to show the potentials of our approach, capable of running on a whole 91 × 109 × 91 × 1,200 sized image taking only few minutes per iteration. Future developments could employ the application of a GM-mask for a functional activation analysis, or not if the interest is to investigate white matter functional activation. It should again be emphasized that, while the GLM requires knowledge of the experimental paradigm, the A4D-fMRI does not. These results highlight that the A4D-fMRI can be used to recover brain activity in the absence of an experimental paradigm, such as rs-fMRI. When comparing correlation maps obtained for the A4D-fMRI and the TA, correlation values obtained with the A4D-fMRI were significantly higher than those obtained with the TA suggesting an improved recovery of brain activity.

Compared to previous approaches, A4D-fMRI considers both the spatial and temporal nature of the fMRI signals simultaneously. This not only leads to an elegant formulation, but also allows the algorithm to leverage the notion smoothness across dimensions. For example, a brief activation of a volume of the brain or a long activation of a small area would both be detected as they correspond to a 3D subset of the 4D data. The algorithm also allows efficient implementation, which leads to reduced computation times when compared to existing strategies. Finally, while we only considered the notion of connectivity the natural representation of fMRI leading to a 4D algorithm, long-range connectivity could potentially be added by increasing the number of dimensions. This could, for example, lead to non-local filtering tied to white matter connectivity obtained from diffusion MRI and embed the notion of brain network naturally.

The A4D-fMRI can be used for different purposes, for example to recover brain activations in a task experimental paradigm as well as in a rs-fMRI study, where the subject is asked not to perform any task while lying in the MRI scanner. Hence, the A4D-fMRI could be useful to analyze the functional brain activity for those subject affected by neurological diseases that make them unable to perform a task or to analyze unexpected brain activities, for example in the case of epilepsy, thus improving the recovery of brain dynamics for future clinical applications. The A4D-fMRI could help in the recovery of time series and spatial maps that could be post-processed afterwards to perform statistical analyses. Another application of A4D-fMRI is the detection of innovation-driven co-activation patterns (iCAPs) (Karahanoğlu and Van De Ville, 2015). Indeed, iCAPs identifies spatially overlapping activation maps that reveal transients in spontaneous neural temporal activity on rs-fMRI data. It has been demonstrated that decomposing rs-fMRI data using iCAPs reveals the spatiotemporal dynamics below the so-called resting-state networks. Within this context, A4D-fMRI could be used to recover the innovation signal, which is derived from the recovered brain activations. A second application of interest is structure–function mapping (Deslauriers-Gauthier et al., 2020b) where A4D-fMRI would serve as a preprocessing step. This could potentially enhance the signal to noise ratio of the functional connectivity matrices and improve the training of models used to predict function from structure. It has been shown that including fMRI information from the A4D-fMRI approach to find active regions to be exploited as priors on cortical regions allows to select plausible structural connections to yield a tractable optimization problem to infer white matter information flow (Deslauriers-Gauthier et al., 2020a). Results, and specifically the ones shown in Figure 10, show the great potentials of the A4D-fMRI to analyze the dynamics of rs-fMRI data (Preti et al., 2017) on a larger sample in healthy control as well as in a patient group. The A4D-fMRI could also be exploited in a framework that employs the dynamic causal modeling (Friston et al., 2003). It must be clarified that the A4D-fMRI does not assume any hypothesis on the interactions between brain regions. Once the activations are inferred using the A4D-fMRI, the dynamic causal modeling could be exploited on the recovered signals to reveal possible causality between brain regions' activity. Hence, the information flow can be inferred between cortical regions known to be active using the A4D-fMRI that allows to blindly identify active regions, without requiring a manual selection of the regions of interest. A drawback of A4D-fMRI is its tendency to smooth the boundary of activity inducing signal as seen in Figure 3, particularly when compared to l1-norm approaches such as TA. It is caused by the local smoothing of the 4D tensor embodied in the parameter σ_G. However, this is a necessary step of the algorithm as it improves the local estimation of the tensor and ensures it is full rank. A second source of smoothing is caused by our choice to set three eigenvalues of the modified structure tensor to 1. Filtering can therefore only be prevented in one direction, which is ideal for many signal boundaries. However, this introduces smoothing in certain cases where ideally fewer eigenvalues should be set to 1 (one such example is the spike train presented in Section 3). We are currently investigating strategies that would allow anisotropic filtering in an adapted number of directions and reduce over-smoothing.

Future works could include into the solution of the inverse problem the information given by the diffusion MRI data that would provide us with a more complex neighborhood defined by white matter connectivity. In this way, the neighborhood would no longer be given only by the surrounding voxels, but also by the voxels which are anatomically segregated and therefore functionally connected to achieve the same function. In addition, the proposed A4D-fMRI approach could be exploited to investigate possible fMRI activations in the white matter, which is an emerging debated topic in the neuroimaging field (Huang et al., 2018).

6. Conclusion

In this article, we proposed and validated a new method to blindly regularize the fMRI images and recover the brain activity from fMRI signals without prior knowledge. Our findings show that the A4D-fMRI enabled us to solve an important problem, coupling the spatial and the temporal dimension and to recover brain activations overlapping the ones obtained with the GLM. Our results also show higher correlations of the recovered time courses with the ground truth compared to the TA. This opens a new channel for the analyses of rs-fMRI data and the recovery of paradigm-free neural activity to be used for investigations in future clinical applications.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found at: https://db.humanconnectome.org/app/template/Login.vm;jsessionid=55C7D675B7D4CA21400384E01478FD3A.

Author Contributions

IC, RD, and SD-G: study conception and design, data collection, analysis and interpretation of results, and manuscript preparation. All authors contributed to the article and approved the submitted version.

Funding

This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (ERC Advanced Grant agreement No 694665: CoBCoM—Computational Brain Connectivity Mapping).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Data were provided [in part] by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research and by the McDonnell Center for Systems Neuroscience at Washington University.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnimg.2022.815423/full#supplementary-material

Footnotes

1. ^gitlab.inria.fr/cobcom/a4dfmri

2. ^https://fsl.fmrib.ox.ac.uk/fsl/fslwiki

3. ^https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/POSSUM

4. ^https://miplab.epfl.ch/index.php/software/total-activation

References

Andersen, A. H., Gash, D. M., Avison, M. J. (1999). Principal component analysis of the dynamic response measured by fMRI: a generalized linear systems framework. Magn. Reson. Imaging 17, 795–815. doi: 10.1016/S0730-725X(99)00028-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Aubert, G., Kornprobst, P. (2006). Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations, Vol. 147. New York, NY: Springer Science and Business Media. doi: 10.1007/978-0-387-44588-5

CrossRef Full Text | Google Scholar

Baumgartner, R., Ryner, L., Richter, W., Summers, R., Jarmasz, M., Somorjai, R. (2000). Comparison of two exploratory data analysis methods for fMRI: fuzzy clustering vs. principal component analysis. Magn. Reson. Imaging 18, 89–94. doi: 10.1016/S0730-725X(99)00102-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Beckmann, C. F., Smith, S. M. (2004). Probabilistic independent component analysis for functional magnetic resonance imaging. IEEE Trans. Med. Imaging 23, 137–152. doi: 10.1109/TMI.2003.822821

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolton, T. A., Farouj, Y., Inan, M., Van De Ville, D. (2019). “Structurally informed deconvolution of functional magnetic resonance imaging data,” in 2019 IEEE 16th International Symposium on Biomedical Imaging (ISBI 2019) (Venice: IEEE), 1545–1549. doi: 10.1109/ISBI.2019.8759218

CrossRef Full Text | Google Scholar

Buckner, R., Andrews-Hanna, J., Schacter, D. (2008). The brain's default network-anatomy, function, and relevance to disease. Cogn. Neurosci. 2008, 1–38. doi: 10.1196/annals.1440.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Caballero Gaudes, C., Petridou, N., Francis, S. T., Dryden, I. L., Gowland, P. A. (2013). Paradigm free mapping with sparse regression automatically detects single-trial functional magnetic resonance imaging blood oxygenation level dependent responses. Hum. Brain Mapp. 34, 501–518. doi: 10.1002/hbm.21452

PubMed Abstract | CrossRef Full Text | Google Scholar

Calhoun, V. D., Adali, T. (2006). Unmixing fMRI with independent component analysis. IEEE Eng. Med. Biol. Mag. 25, 79–90. doi: 10.1109/MEMB.2006.1607672

PubMed Abstract | CrossRef Full Text | Google Scholar

Chambolle, A., Lions, P.-L. (1997). Image recovery via total variation minimization and related problems. Numer. Math. 76, 167–188. doi: 10.1007/s002110050258

CrossRef Full Text | Google Scholar

Charbonnier, P., Blanc-Féraud, L., Aubert, G., Barlaud, M. (1997). Deterministic edge-preserving regularization in computed imaging. IEEE Trans. Image Process. 6, 298–311. doi: 10.1109/83.551699

PubMed Abstract | CrossRef Full Text | Google Scholar

Cordes, D., Haughton, V., Carew, J. D., Arfanakis, K., Maravilla, K. (2002). Hierarchical clustering to measure connectivity in fMRI resting-state data. Magn. Reson. Imaging 20, 305–317. doi: 10.1016/S0730-725X(02)00503-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Costantini, I., Filipiak, P., Maksymenko, K., Deslauriers-Gauthier, S., Deriche, R. (2018). “fMRI deconvolution via temporal regularization using a lassomodel and the LARS algorithm,” in 40th International Engineering in Medicine and Biology Conference (Honolulu, HI).

Google Scholar

Deslauriers-Gauthier, S., Costantini, I., Deriche, R. (2020a). Non-invasive inference of information flow using diffusion MRI, functional MRI, and MEG. J. Neural Eng. 17:045003. doi: 10.1088/1741-2552/ab95ec

PubMed Abstract | CrossRef Full Text | Google Scholar

Deslauriers-Gauthier, S., Zucchelli, M., Frigo, M., Deriche, R. (2020b). A unified framework for multimodal structure-function mapping based on eigenmodes. Med. Image Anal. 66:101799. doi: 10.1016/j.media.2020.101799

PubMed Abstract | CrossRef Full Text | Google Scholar

Farouj, Y., Karahanoǧlu, F. I., Van De Ville, D. (2017). “Regularized spatiotemporal deconvolution of fMRI data using gray-matter constrained total variation,” in 2017 IEEE 14th International Symposium on Biomedical Imaging (Melbourne, VIC: IEEE), 472–475. doi: 10.1109/ISBI.2017.7950563

CrossRef Full Text | Google Scholar

Flandin, G., Penny, W. D. (2007). Bayesian fMRI data analysis with sparse spatial basis function priors. NeuroImage 34, 1108–1125. doi: 10.1016/j.neuroimage.2006.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Förstner, W.. (1986). “A feature based correspondence algorithm for image matching,” in ISPRS ComIII (Rovaniemi), 150–166.

Google Scholar

Förstner, W., Gülch, E. (1987). “A fast operator for detection and precise location of distinct points, corners and centres of circular features,” in Proceedings of ISPRS Intercommission Conference on Fast Processing of Photogrammetric Data (Interlaken), 281–305.

Google Scholar

Friston, K. J., Fletcher, P., Josephs, O., Holmes, A., Rugg, M., Turner, R. (1998). Event-related fMRI: characterizing differential responses. Neuroimage 7, 30–40. doi: 10.1006/nimg.1997.0306

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K. J., Harrison, L., Penny, W. (2003). Dynamic causal modelling. Neuroimage 19, 1273–1302. doi: 10.1016/S1053-8119(03)00202-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K. J., Mechelli, A., Turner, R., Price, C. J. (2000). Nonlinear responses in fMRI: the balloon model, volterra kernels, and other hemodynamics. NeuroImage 12, 466–477. doi: 10.1006/nimg.2000.0630

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K. J., Trujillo-Barreto, N., Daunizeau, J. (2008). DEM: a variational treatment of dynamic systems. Neuroimage 41, 849–885. doi: 10.1016/j.neuroimage.2008.02.054

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaudes, C. C., Petridou, N., Dryden, I. L., Bai, L., Francis, S. T., Gowland, P. A. (2011). Detection and characterization of single-trial fMRI bold responses: paradigm free mapping. Hum. Brain Mapp. 32, 1400–1418. doi: 10.1002/hbm.21116

PubMed Abstract | CrossRef Full Text | Google Scholar

Gitelman, D. R., Penny, W. D., Ashburner, J., Friston, K. J. (2003). Modeling regional and psychophysiologic interactions in fMRI: the importance of hemodynamic deconvolution. Neuroimage 19, 200–207. doi: 10.1016/S1053-8119(03)00058-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Glasser, M. F., Sotiropoulos, S. N., Wilson, J. A., Coalson, T. S., Fischl, B., Andersson, J. L., et al. (2013). The minimal preprocessing pipelines for the human connectome project. Neuroimage 80, 105–124. doi: 10.1016/j.neuroimage.2013.04.127

PubMed Abstract | CrossRef Full Text | Google Scholar

Glover, G. H.. (1999). Deconvolution of impulse response in event-related bold fMRI1. Neuroimage 9, 416–429. doi: 10.1006/nimg.1998.0419

PubMed Abstract | CrossRef Full Text | Google Scholar

Golland, Y., Golland, P., Bentin, S., Malach, R. (2008). Data-driven clustering reveals a fundamental subdivision of the human cortex into two global systems. Neuropsychologia 46, 540–553. doi: 10.1016/j.neuropsychologia.2007.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Gusnard, D. A., Raichle, M. E. (2001). Searching for a baseline: functional imaging and the resting human brain. Nat. Rev. Neurosci. 2:685. doi: 10.1038/35094500

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrison, L. M., Penny, W., Flandin, G., Ruff, C. C., Weiskopf, N., Friston, K. J. (2008). Graph-partitioned spatial priors for functional magnetic resonance images. NeuroImage 43, 694–707. doi: 10.1016/j.neuroimage.2008.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Havlicek, M., Friston, K. J., Jan, J., Brazdil, M., Calhoun, V. D. (2011). Dynamic modeling of neuronal responses in fMRI using cubature Kalman filtering. Neuroimage 56, 2109–2128. doi: 10.1016/j.neuroimage.2011.03.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Hernandez-Garcia, L., Ulfarsson, M. O. (2011). Neuronal event detection in fMRI time series using iterative deconvolution techniques. Magn. Reson. Imaging 29, 353–364. doi: 10.1016/j.mri.2010.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, Y., Bailey, S. K., Wang, P., Cutting, L. E., Gore, J. C., Ding, Z. (2018). Voxel-wise detection of functional networks in white matter. NeuroImage 183, 544–552. doi: 10.1016/j.neuroimage.2018.08.049

PubMed Abstract | CrossRef Full Text | Google Scholar

Karahanoǧlu, F. I., Caballero-Gaudes, C., Lazeyras, F., Van De Ville, D. (2013). Total activation: fMRI deconvolution through spatio-temporal regularization. Neuroimage 73, 121–134. doi: 10.1016/j.neuroimage.2013.01.067

PubMed Abstract | CrossRef Full Text | Google Scholar

Karahanoǧlu, F. I., Van De Ville, D. (2015). Transient brain activity disentangles fMRI resting-state dynamics in terms of spatially and temporally overlapping networks. Nat. Commun. 6:7751. doi: 10.1038/ncomms8751

PubMed Abstract | CrossRef Full Text | Google Scholar

Khalidov, I., Fadili, J., Lazeyras, F., Van De Ville, D., Unser, M. (2011). Activelets: wavelets for sparse representation of hemodynamic responses. Sign. Process. 91, 2810–2821. doi: 10.1016/j.sigpro.2011.03.008

CrossRef Full Text | Google Scholar

Kimmel, R., Malladi, R., Sochen, N. (2000). Images as embedded maps and minimal surfaces: movies, color, texture, and volumetric medical images. Int. J. Comput. Vis. 39, 111–129. doi: 10.1023/A:1008171026419

CrossRef Full Text | Google Scholar

Kiviniemi, V., Kantola, J.-H., Jauhiainen, J., Hyvärinen, A., Tervonen, O. (2003). Independent component analysis of nondeterministic fMRI signal sources. Neuroimage 19, 253–260. doi: 10.1016/S1053-8119(03)00097-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, M. H., Hacker, C. D., Snyder, A. Z., Corbetta, M., Zhang, D., Leuthardt, E. C., et al. (2012). Clustering of resting state networks. PLoS ONE 7:e40370. doi: 10.1371/journal.pone.0040370

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindeberg, T.. (2013). Scale-Space Theory in Computer Vision, Vol. 256. Boston, MA: Springer Science and Business Media.

Google Scholar

Liu, Y., Gao, J.-H., Liu, H.-L., Fox, P. T. (2000). The temporal response of the brain after eating revealed by functional MRI. Nature 405:1058. doi: 10.1038/35016590

PubMed Abstract | CrossRef Full Text | Google Scholar

Lund, T. E., Madsen, K. H., Sidaros, K., Luo, W.-L., Nichols, T. E. (2006). Non-white noise in fMRI: does modelling have an impact? Neuroimage 29, 54–66. doi: 10.1016/j.neuroimage.2005.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

McKeown, M. J., Makeig, S., Brown, G. G., Jung, T.-P., Kindermann, S. S., Bell, A. J., et al. (1998). Analysis of fMRI data by blind separation into independent spatial components. Hum. Brain Mapp. 6, 160–188. doi: 10.1002/(SICI)1097-0193(1998)6:3<160::AID-HBM5>3.0.CO;2-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Morgan, V. L., Li, Y., Abou-Khalil, B., Gore, J. C. (2008). Development of 2DTCA for the detection of irregular, transient bold activity. Hum. Brain Mapp. 29, 57–69. doi: 10.1002/hbm.20362

PubMed Abstract | CrossRef Full Text | Google Scholar

Nielsen, M., Florack, L., Deriche, R. (1997). Regularization, scale-space, and edge detection filters. J. Math. Imaging Vis. 7, 291–307. doi: 10.1023/A:1008282127190

CrossRef Full Text | Google Scholar

Penfield, W., Boldrey, E. (1937). Somatic motor and sensory representation in the cerebral cortex of man as studied by electrical stimulation. Brain 60, 389–443. doi: 10.1093/brain/60.4.389

CrossRef Full Text | Google Scholar

Perona, P., Malik, J. (1990). Scale-space and edge detection using anisotropic diffusion. IEEE Trans. Pattern Anal. Mach. Intell. 12, 629–639. doi: 10.1109/34.56205

CrossRef Full Text | Google Scholar

Preti, M. G., Bolton, T. A., Van De Ville, D. (2017). The dynamic functional connectome: state-of-the-art and perspectives. Neuroimage 160, 41–54. doi: 10.1016/j.neuroimage.2016.12.061

PubMed Abstract | CrossRef Full Text | Google Scholar

Riera, J. J., Watanabe, J., Kazuki, I., Naoki, M., Aubert, E., Ozaki, T., et al. (2004). A state-space model of the hemodynamic approach: nonlinear filtering of bold signals. NeuroImage 21, 547–567. doi: 10.1016/j.neuroimage.2003.09.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Roux, F.-E., Djidjeli, I., Durand, J.-B. (2018). Functional arctschumperle2007anisotropichitecture of the somatosensory homunculus detected by electrostimulation. J. Physiol. 596, 941–956. doi: 10.1113/JP275243

PubMed Abstract | CrossRef Full Text | Google Scholar

Rudin, L. I., Osher, S., Fatemi, E. (1992). Nonlinear total variation based noise removal algorithms. Phys. D 60, 259–268. doi: 10.1016/0167-2789(92)90242-F

CrossRef Full Text | Google Scholar

Salvador, R., Suckling, J., Coleman, M. R., Pickard, J. D., Menon, D., Bullmore, E. (2005). Neurophysiological architecture of functional magnetic resonance images of human brain. Cereb. Cortex 15, 1332–1342. doi: 10.1093/cercor/bhi016

PubMed Abstract | CrossRef Full Text | Google Scholar

Sapiro, G.. (2006). Geometric Partial Differential Equations and Image Analysis. Cambridge: Cambridge University Press.

Google Scholar

Smith, M., Fahrmeir, L. (2007). Spatial Bayesian variable selection with application to functional magnetic resonance imaging. J. Am. Stat. Assoc. 102, 417–431. doi: 10.1198/016214506000001031

PubMed Abstract | CrossRef Full Text | Google Scholar

Tournier, J.-D., Calamante, F., Connelly, A. (2012). Mrtrix: diffusion tractography in crossing fiber regions. Int. J. Imaging Syst. Technol. 22, 53–66. doi: 10.1002/ima.22005

PubMed Abstract | CrossRef Full Text | Google Scholar

Tschumperle, D., Deriche, R. (2002). Diffusion PDES on vector-valued images. IEEE Sign. Process. Mag. 19, 16–25. doi: 10.1109/MSP.2002.1028349

PubMed Abstract | CrossRef Full Text | Google Scholar

Tschumperle, D., Deriche, R. (2005). Vector-valued image regularization with PDES: a common framework for different applications. IEEE Trans. Pattern Anal. Mach. Intell. 27, 506–517. doi: 10.1109/TPAMI.2005.87

PubMed Abstract | CrossRef Full Text | Google Scholar

Tschumperlé, D., Deriche, R. (2007). Anisotropic diffusion partial differential equations for multichannel image regularization: framework and applications. Adv. Imaging Electr. Phys. 145, 149–209. doi: 10.1016/S1076-5670(06)45004-7

CrossRef Full Text | Google Scholar

Uruńnuela, E., Jones, S., Crawford, A., Shin, W., Oh, S., Lowe, M., et al. (2020). “Stability-based sparse paradigm free mapping algorithm for deconvolution of functional MRI data,” in 2020 42nd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (Montreal, QC: IEEE), 1092–1095. doi: 10.1109/EMBC44109.2020.9176137

PubMed Abstract | CrossRef Full Text | Google Scholar

Uruńnuela, E., Moia, S., Caballero-Gaudes, C. (2021). “A low rank and sparse paradigm free mapping algorithm for deconvolution of fMRI data,” in 2021 IEEE 18th International Symposium on Biomedical Imaging (Nice: IEEE), 1726–1729. doi: 10.1109/ISBI48211.2021.9433821

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Essen, D. C., Smith, S. M., Barch, D. M., Behrens, T. E., Yacoub, E., Ugurbil, K., et al. (2013). The Wu-minn human connectome project: an overview. Neuroimage 80, 62–79. doi: 10.1016/j.neuroimage.2013.05.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Weickert, J.. (1998). Anisotropic Diffusion in Image Processing, Vol. 1. Stuttgart: Teubner Stuttgart.

Google Scholar

Weickert, J.. (1999). Coherence-enhancing diffusion of colour images. Image Vis. Comput. 17, 201–212. doi: 10.1016/S0262-8856(98)00102-4

CrossRef Full Text | Google Scholar