# Mixture Components Inference for Sparse Regression: Introduction and
Application for Estimation of Neuronal Signal from fMRI BOLD

Anna Pidnebesna  Iveta Fajnerová  Jiří Horáček  and Jaroslav Hlinka This work
was supported by Czech Health Research Council project No. NV17-28427A, the
Czech Science Foundation projects No. 19-11753S, 19-07164S and 21-32608S and
by Supported by Ministry of Health Czech Republic – DRO (NUDZ, 00023752)
(Corresponding author: Jaroslav Hlinka). A. Pidnebesna and J. Hlinka are with
the Institute of Computer Science of the Czech Academy of Sciences, Pod
Vodárenskou věží 271/2, 182 07 Prague, Czech Republic (email:
<EMAIL_ADDRESS>hlinka@cs.cas.cz).A. Pidnebesna and Iveta Fajnerová and
Jiří Horáček and J. Hlinka are with the National Institute of Mental Health,
Topolová 748, 250 67 Klecany, Czech Republic<EMAIL_ADDRESS>jiri.horacek@nudz.cz).

###### Abstract

Sparse linear regression methods including the well-known LASSO and the
Dantzig selector have become ubiquitous in the engineering practice, including
in medical imaging. Among other tasks, they have been successfully applied for
the estimation of neuronal activity from functional magnetic resonance data
without prior knowledge of the stimulus or activation timing, utilizing an
approximate knowledge of the hemodynamic response to local neuronal activity.
These methods work by generating a parametric family of solutions with
different sparsity, among which an ultimate choice is made using an
information criteria. We propose a novel approach, that instead of selecting a
single option from the family of regularized solutions, utilizes the whole
family of such sparse regression solutions. Namely, their ensemble provides a
first approximation of probability of activation at each time-point, and
together with the conditional neuronal activity distributions estimated with
the theory of mixtures with varying concentrations, they serve as the inputs
to a Bayes classifier eventually deciding on the verity of activation at each
time-point. We show in extensive numerical simulations that this new method
performs favourably in comparison with standard approaches in a range of
realistic scenarios. This is mainly due to the avoidance of overfitting and
underfitting that commonly plague the solutions based on sparse regression
combined with model selection methods, including the corrected Akaike
Information Criterion. This advantage is finally documented in selected fMRI
task datasets.

fMRI, BOLD, Deconvolution, Mixtures with varying concentrations, Neuronal
signal estimation.

##  1 Introduction

The most common approach to fMRI measurement is using the Blood Oxygenation
Level Dependent (BOLD) contrast, see [1]. Soon after its discovery, it was
shown that it reliably reflects changes in local brain activity related to
spontaneous or stimulation-triggered cognitive processing [2]. However, the
BOLD signal is only indirectly related to the underlying neuronal activity. As
its name suggests, it reflects the changes in the blood oxygenation level;
these happen after the triggering neuronal activity and constitute a
considerably delayed and smoothed response that can be in the first
approximation modelled as the output of a linear filter represented by a
hemodynamic response function. Combined with a range of noise sources in fMRI
[3], statistical inference of neuronal activity from fMRI BOLD signal poses a
serious signal analysis challenge.

The common fMRI analysis approach is to impose strong prior on the timing of
the tentative neuronal activity, assuming that it amounts to a predefined
time-course reflecting the presence of particular controlled experimental
manipulation. This expected activation time course is convolved with a
template hemodynamic response and used as a predictor in a linear regression
model of the observed BOLD fMRI signal.

Yet there are still many scenarios when this approach is not suitable,
typically when the precise timing of the presumed neuronal activity is not
known as it is not perfectly locked to a predefined task paradigm, but rather
inter- and intra-individually variable, self-paced or ultimately spontaneously
occurring. The importance of having tools available to analyse such less
constrained fMRI BOLD data (compared to regular block and event-related
designs) is increasingly recognized as a valuable complement of the standard
analysis methods [4].

This approach is opposite to the common analysis – instead of convolving the
expected neuronal activity with the hemodynamic response to obtain an
approximate expected BOLD fMRI signal (and comparing it with the observed BOLD
fMRI via linear regression), the observed signal is _deconvolved_ to obtain an
estimate of the underlying neuronal signal. This ’deconvolution’ task can be
formulated as a linear regression problem, solvable by various approaches;
among them the so-called sparse regression, that attempts to preclude
overfitting of noise by inducing penalties for marking too many activations in
the data, is most common. The most used sparse regression methods are the
LASSO (least absolute shrinkage and selection operator), see [5], and the
Dantzig Selector, proposed in [6]. These methods output a parametric family of
fitted models (parameterized by a free parameter determining the strength of
the regularization and therefore the sparsity of the obtained solution). Thus
a suitable model selection procedure must be used, typically the Akaike
Information Criteria (AIC) [7] or Bayes Information Criteria (BIC) [8].

While these model selection criteria have been extensively and successfully
applied in many fields, they have their drawbacks. In particular, for small
sample size AIC tends to overfit the data, giving a too high model order,
whereas the BIC-selected model may underfit data, see [9]. Therefore the use
of the adjusted Akaike Information Criterion (AICc) [10], [11] was recommended
instead of AIC.

The use of sparse linear regression (specifically the Dantzig Selector) for
neuronal signal estimation is described in [4]. The authors compare AIC and
BIC for the neuronal signal estimation in a simulation study and show that AIC
fails to estimate the sparse stimulus signal. Therefore they suggest using the
BIC. However, the simulated stimulus was very sparse in the study, suggesting
that in the case of less sparse activations, BIC might be underfitting the
data.

We assessed the performance of fMRI deconvolution methods including the
canonical LASSO and the Dantzig Selector regularization combined with the most
common criteria (AIC, BIC and AICc), arriving at a conclusion that their
performance is far from perfect under realistic scenarios, possibly due to
nonrobustness of the model selection procedure. In this work we thus propose
an alternative approach to the estimation, based on the Bayes classifier and
theory of mixtures with varying concentrations, that provides more accurate
and robust results both in simulations and on a real data sample. The
described algorithm has been implemented into an updated version of software
BRAD [12]. The work is organized as follows. The model of BOLD signal, is
described in Section 2. The practical motivation of the study and theory
overview is given in Section 3. In Section 4 the new approach is proposed, and
methodology and simulation study are described. Finally, the real data
application is found in Section 6.

##  2 Data Model

We model the measured BOLD signal $y$ as a result of a hemodynamic response to
local brain activity $s$. The neuronal activity is represented by (a sparse)
vector of activations; we assume there is no a priori information about
stimulus timing. In line with the common fMRI BOLD modelling approaches, we
assume that the measured signal is a result of the convolution of the neuronal
activity and the hemodynamic response function (HRF), in mathematical
notation:

$y=h\star s+e,\quad y_{k}=\sum_{i=0}^{N_{h}-1}h_{i}s_{k-i}+e_{k},\quad
k\in\\{1,..,N\\},$ (1)

where $k\in\\{1,..,N\\}$ indexes the time, $\star$ denotes the convolution,
$y,s,e$ are $N\times 1$ vectors denoting the measured signal, the underlying
brain activity signal and the observational noise, respectively. The vector
$h$ denotes the HRF, where the length $N_{h}$ is smaller than the length of
the measured signal: $N_{h}<N$. We assume that $e$ is white noise, independent
and identically distributed in time. The model can be rewritten as

$y=Hs+e,$ (2)

where $H$ is the Toeplitz convolution matrix of size $N\times N$ corresponding
to the HRF $h$ (see [13]).

##  3 Motivation/Theoretical background

###  3.1 Sparse Regression Methods

We use the two most prominent sparse regression methods - the LASSO (least
absolute shrinkage and selection operator) and the Dantzig Selector as the
benchmarks. LASSO estimator was proposed in [5] as a regression analysis
method that performs variable selection. The aim (in notation congruent with
the deconvolution introduction) is to estimate a sparse vector of linear
regression coefficients $s$ in linear regression:

$y=Hs+e,$

where $y$ is the dependent variable, $H$ is a matrix of regressors, and $e$ is
an additive i.i.d. Gaussian noise. Then the LASSO problem is formulated as

$\hat{s}=\arg\min_{s}\Big{[}\lambda\|s\|_{1}+\|y-Hs\|^{2}_{2}\Big{]},$ (3)

where $\lambda\geq 0$ is a tuning parameter and $\|.\|_{p}$ is a norm in
$L_{p}$-space, i.e. $\|x\|_{p}=\big{(}\sum_{i=1}^{T}|x_{i}|^{p}\big{)}^{1/p}$.
This is an instance of convex optimization and quadratic programming.

The Dantzig Selector was proposed in [6] as:

$\hat{s}=\arg\min_{s}\|s\|_{1},\quad\|H^{T}(y-Hs)\|_{\infty}<\lambda,$ (4)

where $\lambda\geq 0$ is a tuning parameter, $\|.\|_{\infty}$ is a norm in
$L_{\infty}$-space: $\|x\|_{\infty}=\sup_{i}|x_{i}|$. Relation of LASSO and
Dantzig Selector estimators was discussed e.g. in [14], [15], [16].

###  3.2 Selection Criteria

In the model selection context, we assume that for given data, we obtain a set
of candidate models. The selection criteria assess the relative quality of the
candidate models by finding a compromise between accuracy of the data fit and
model parsimony. The model with the lowest value of the selection criteria
among the tested models is considered the best.

The Akaike Information Criteria (AIC) is defined as

$AIC=-2\ln L+2k,$

where $L$ is the maximised value of the likelihood function of the model, and
$k$ is the number of estimated model parameters.

The Bayes Infromation Criteria (BIC) is defined as

$BIC=-2\ln L+k\ln N,$

note that it also takes into account the sample size $N$.

The adjusted Akaike Information Criterion (AICc) is recommended instead of AIC
(see [11], [9]) when $k$ is large relative to sample size $N$. It is defined
as

$AIC_{c}=AIC+\frac{2k(k+1)}{N-k-1}.$

For LASSO and Dantzig Selector, AIC and BIC give

$\lambda^{AIC}=\arg\min_{\lambda}\Big{[}N\ln(\lVert
y-H\hat{s}(\lambda)\lVert_{2}^{2})+2k(\lambda)\Big{]},$ (5)
$\lambda^{BIC}=\arg\min_{\lambda}\Big{[}N\ln(\lVert
y-H\hat{s}(\lambda)\lVert_{2}^{2})+k(\lambda)\ln N\Big{]},$ (6)

where $\hat{s}(\lambda)$ is the solution of (3) and (4), that corresponds to
$\lambda$, $y$ is the measured signal, $H$ is the matrix of regressors,
$k(\lambda)$ is the number of non-zero parameters in $\hat{s}(\lambda)$ (see
[17]). Comparison of using the AIC and BIC to Dantzig Selector in
deconvolution tasks in neuroimaging was discussed in [4].

##  4 MCI

###  4.1 Basic idea: Bayes classifier

The idea of our proposed approach is to utilize the information available in
the observed data and the derived signals (namely the ordinary least square
estimate of the brain activity and the family of regularized estimates) within
a Bayesian inference tool that separates the true activations from false
positives. The starting point is a basic estimate of the brain activity
$\xi=(\xi_{1},\dots,\xi_{N})$, $\xi_{j}\in\mathbb{R}$ (think of e.g. the OLS
solution $\hat{s}_{OLS}$). We assume that its non-zero elements correspond to
apparent activations due to one of two cases: neuronal response (true
activation) or noise (false activation), and apply a naive Bayes classifier to
separate the two classes.

Let us discuss the Bayes classifier in more detail. Let
$g:\mathbb{R}\rightarrow\\{1,\dots,M\\}$ be the function (classifier), which
returns the number of mixture component (class) $g(\xi)\in\\{1,\dots,M\\}$ for
each possible value of the observed characteristics $\xi\in\mathbb{R}$. In our
case, there are two classes, i.e. $M=2,\quad g(\xi_{j})\in\\{1,2\\},$ where
$g(\xi_{j})=1$ corresponds to the conclusion ’$j$-th element falls to a true
activation component’, and $g(\xi_{j})=2$ corresponds to the conclusion
’$j$-th element falls to a false detection component’. The Bayes classifier is
of the form

$g^{B}(\xi_{j})=\arg\max_{m=1,\dots,M}\,p_{j}^{m}f^{m}(\xi_{j}),$ (7)

where $j\in\\{1,\dots,N\\}$ indexes the classified elements of the vector,
$\xi_{j}\in\mathbb{R}$ are the basic (OLS) estimate values, $p_{j}^{m}$ are
the probabilities for $j$-th activation to fall to the $m$-th class, and
$f^{m}$ denotes the probability densities of the characteristic $\xi$,
provided it falls within the class $m$. Thus, the classifier $g^{B}(.)$ for
each $j$-th activation compares the posterior probabilities of the element
belonging to a neuronal activation or noise. For details and properties of
such a classifier see [18] and [19].

However, true values of the probabilities $p^{m}$ and the densities $f^{m}(x)$
are usually unknown a priori. Instead, we may need to obtain some estimates
$\hat{p}^{m}$ and $\hat{f}^{m}(x)$ from data. We suggest two approaches to
obtain such estimates. First, by using the Gaussian mixture model (GMM),
described in [20]. It is based on the assumption of the normal distribution of
all mixture components. This approach is well-known, and its main advantage is
a developed theory and existence of software tools for the automatic
estimation of all the parameters of the model. For our purposes, it is
important to estimate both of the needed elements, probability vectors
$\hat{p}^{m}$ and component densities $\hat{f}^{m}(x)$. We test this approach
in Section 5.

The second approach is based on utilizing a model with varying concentrations
(MVC) for estimating the component densities, and a heuristic procedure based
on LASSO/Dantzig Selector algorithms for approximating the probabilities.
Using this concept enables us, not only to disregard the assumption about the
normality of the component distributions, but it also allows for utilizing
element-specific prior probabilities – that are in the present scenario
informed by the regularized regressions solutions. The detailed description of
this concept is provided below and its block-scheme is shown in the left panel
of Fig. 1.

Figure 1: Left panel: Block-scheme of the compared neural signal estimation
approaches. The BOLD signal is first used to estimate the OLS signal estimate
$s_{OLS}$ and the family of sparse estimates $s_{DS}(\lambda)$ parametrized by
the regularization parameter $\lambda$. Then, either [left branch] the GMM can
be applied to the OLS estimate $s_{OLS}$, to provide an estimate for $p^{m}$,
$f^{m}(x)$ and ultimately a signal estimate, or [right branch] one of the
AIC/AICc/BIC selection criteria can be applied to the family of sparse
estimates $s_{DS}(\lambda)$ to select a particular sparse estimate, or finally
[middle branch], MCI uses $s_{DS}(\lambda)$ to estimate ${p_{j}}^{m}$, this is
used together with the OLS estimate $s_{OLS}$ to estimate $f^{m}(x)$, and
these are then subjected to a Bayesian Classifier to define the final neural
signal estimate. Right panel: scheme of simulation-based comparisons. From the
simulated neuronal activity $s$ the BOLD signal $y$ is generated by
convolution and adding noise. Then, different sparse methods are applied:
LASSO/DS with AIC/AICc/BIC selection criteria, or OLS estimate followed by
removing noisy activations using GMM or MCI. These three types of estimates
are then compared against the ground truth.

###  4.2 MVC-based density estimates

In the classical finite mixture model, the mixing probabilities and thus
probability densities of all features $\xi_{j}$ are supposed to be the same
for all the observed objects $j$ [20]:

$f_{j}(\xi_{j})=f(\xi_{j})=\sum_{m=1}^{M}p^{m}f^{m}(\xi_{j}).$ (8)

As a generalization of this approach, the model with varying concentrations
(MVC) was proposed [21]. According to this model, mixing probabilities
$p^{m}_{j},\,m\in\\{1,\dots,M\\},\,j\in\\{1,\dots,N\\}$, vary for different
objects $j$, although still satisfying the condition: $\forall j:0\leq
p^{m}_{j}\leq 1$, $\sum_{m=1}^{M}p_{j}^{m}=1.$ Then the distribution for the
$j$-th object can be modelled as

$f_{j}(\xi_{j})=\sum_{m=1}^{M}p^{m}_{j}f^{m}(\xi_{j}).$ (9)

In MVC, the mixing probabilities $p^{m}_{j}$ are supposed to be known, whereas
the densities of components $f^{m}$ are unknown and estimated from the data
using the kernel estimate

$\hat{f}^{m}(x)=\frac{1}{b_{m}N}\sum_{j=1}^{N}a^{m}_{j}K\Big{(}\frac{x-\xi_{j}}{b_{m}}\Big{)},$
(10)

where $K(.)$ is a kernel function and $b$ is the bandwidth and $a^{m}_{j}$ are
the weights, aimed to distinguish the $m$-th component and suppress the
influence of other components on it (see [21]). The parameters are obtained
using estimates developed in [22] and [23]. Classification of components of
mixture are discussed in [24] and in [25]. In particular, let us show the
formulas for weights’ computation. We denote the inner product of two vectors
$a=(a_{1},...,a_{N})$, and $b=(b_{1},...,b_{N})$ by $\langle
a,b\rangle=\frac{1}{N}\sum_{i=1}^{N}a_{i}b_{i}$ and define the matrix

$\Gamma_{N}=(\langle p^{l},p^{m}\rangle_{N})_{l,m=1}^{M}.$

Note that this matrix is assumed to be invertible, with $\gamma_{km}$ being
the $km$-minor of $\Gamma_{N}$. Then the weights are given by:

$a^{k}_{j}=\frac{1}{\det\Gamma_{N}}\sum_{m=1}^{M}(-1)^{m+k}\gamma_{km}p^{m}_{j}.$

For the case of two-components, the weights $a^{m}_{j}$ are:

$a^{1}_{j}=\frac{(1-\langle p^{1},1\rangle)p^{1}_{j}+(\langle
p^{1},p^{1}\rangle-\langle p^{1},1\rangle)}{\langle
p^{1},p^{1}\rangle-(\langle p^{1},1\rangle)^{2}},$ $a^{2}_{j}=\frac{\langle
p^{1},p^{1}\rangle-\langle p^{1},1\rangle p^{1}_{j}}{\langle
p^{1},p^{1}\rangle-(\langle p^{1},1\rangle)^{2}},\quad j=1,...,N.$

The kernel width $b_{m}$ is estimated using the Silverman’s rule-of-thumb
(see, for example, [26]):

$b_{m}\approx 1.06\hat{\sigma}_{m}N^{-1/5},$

where $\hat{\sigma}_{m}$ and $\hat{\mu}_{m}$ are the estimated standard
deviation and mean of the $m$-th component given by
$\hat{\sigma}_{m}^{2}=1/(N-1)\sum_{j=1}^{N}a_{j}^{m}(\xi_{j}-\mu_{m})^{2}$,
and $\hat{\mu}_{m}=1/N\sum_{j=1}^{N}a_{j}^{m}\xi_{j}$.

###  4.3 Probability estimates

As described above, for the classification of signal values $s_{j}$, we need
to set the probability $p^{m}_{j}$ that it belongs to the true or false
activation class. In this section, we propose a procedure based on the
numerical solutions of the optimisation problems with a regularisation
parameter, in particular the homotopy algorithms for LASSO (3) and Dantzig
Selector (4) estimators.

Let us denote $s(\lambda)$ the solution of (3) (or (4)) corresponding to the
value $\lambda$. The support of $s(\lambda)$ and its sign vector are

$\mathcal{B(\lambda)}=\text{supp}\,s(\lambda)=\\{j|1\leq j\leq
N,\,s(\lambda)_{j}\neq 0\\},$
$sgn(\lambda)=\\{sgn(s(\lambda)_{j}),\,j=1,\dots,N\\}.$

For both mentioned estimators, the particular solution $s(\lambda)$, its
support $\mathcal{B}(\lambda)$ and sign vector $sgn(\lambda)$ are specified by
the parameter value $\lambda$. For the LASSO problem, it was shown in [27],
that for a given response vector $y$, there is a finite sequence of
$\lambda$’s, called transition points (see [17]),

$\lambda_{0}>\lambda_{1}>\dots>\lambda_{K}=0,\quad\text{such that:}$

  * •

for all $\lambda>\lambda_{0}$, solution of the optimization problem (3)
corresponds to the trivial solution, i.e. its support is the empty set
$\mathcal{B}(\lambda)=\emptyset$;

  * •

the value $\lambda_{K}=0$ corresponds to the ordinary least squares solution;

  * •

for $\lambda$ values from the interval $I_{m}:=(\lambda_{m+1},\lambda_{m})$,
the solutions of (3) are constant with respect to $\lambda$ in terms of their
supports and sign vectors;

$\mathcal{B}_{m}:=\mathcal{B}(\lambda),\quad
sgn_{m}:=sgn(\lambda)\quad\forall\lambda\in I_{m}.$

In other words, there is a sequence $\mathcal{B}_{m}$, $m=0,\dots,K$ of
subsets of indices $j,j=1,\dots,N$ that correspond to nonzero coordinates of
$s(\lambda)$ for $\forall\lambda\in I_{m}$. Thus, for each $j=1,\dots,N$ there
is a set of intervals $I_{m}$, where the $j$ lies in the corresponding
solution support $\mathcal{B}_{m}$:

$\mathcal{I}_{j}:=\cup_{\\{m:j\in\mathcal{B}_{m}\\}}I_{m}=\cup_{\\{m:j\in\mathcal{B}_{m}\\}}(\lambda_{m+1},\lambda_{m}).$

Let us denote the union of all intervals $I_{m}$ as
$\mathcal{I}:=\cup_{m=0}^{K}I_{m}=(0,\lambda_{0})$. We propose to approximate
the probability $\omega_{j}$ that $j$-th activation comes from the true
neuronal signal as

$\hat{p}^{1}_{j}=\frac{|\mathcal{I}_{j}|}{|\mathcal{I}|}=\frac{\sum_{m=0}^{K-1}(\lambda_{m}-\lambda_{m+1})\mathbb{I}\\{j\in\mathcal{B}_{m}\\}}{\lambda_{0}}.$
(11)

This approach corresponds to the idea, that more frequent appearance of non-
zero activation estimate in the solution sets means higher probability for
this moment to correspond to a true neuronal activation. Although the
probabilistic interpretation of this ’relative occurence in the solution set’
is clearly heuristic, it is reasonably motivated as the LASSO operator is
assumed to detect true variables among the noise.

Numerically, the set of transition points and corresponding sets of solutions
for LASSO can be obtained by the homotopy algorithm LARS, presented in [27].
For Dantzig Selector, there is an available homotopy algorithm Primal-Dual
Pursuit, which was presented in [16]. Note that the Primal-Dual Pursuit
produces a finite sequence of so-called critical values of the parameters,
that includes the points of changes in the support or sign vectors of the
solutions; without loss of generality, these critical values can be considered
as the transition points.

Thus, the probability $p_{j}^{1}$ for $j$-th activation to come from true
activation, can be approximated based on the numerical solution of the
optimisation problem with regularisation parameters. Then, the probability
$p_{j}^{2}=1-p_{j}^{1}$, and the densities $f^{1}(.),f^{2}(.)$ can be
estimated. Finally, this information is enough to build the classifier and
separate the components of vector $\xi$ on true and noisy activations.

##  5 Simulation study

To assess the newly proposed methods against the standard methods we first use
numerical simulations. We generate a signal $y=s\star h+e$ using HRF $h$,
while parametrically varying the the properties of the signal $s$ and noise
$e$. Note that the Balloon-Windkessel continuous model can be used for
generation of even more realistic BOLD signals; we compared it with our data
model from Section 2, and found out that the correlation between convolutional
BOLD and the signal generated using the Balloon-Windkessel model is around
0.99.

After the BOLD signal is generated, LASSO and Dantzig Selector with
AIC/AICc/BIC selection criteria are used to estimate the input signal $s$.
These estimates are compared with two approaches based on Bayes Classifier:
GMM (components estimated by (8)) and the proposed method Mixture Components
Inference (MCI, components estimated by (9)). To assess the quality of the
obtained estimates, several measures are used. General scheme is shown in the
right panel of Fig. 1.

The HRF kernel is obtained using the function spm_hrf(.) in the SPM toolbox
(The Wellcome Dept. of Cognitive Neurology, University College London),
setting scan repeat time (RT) to 2.5s. To define the level of added noise $e$,
we use the Signal-to-Noise Ratio (SNR), $SNR=\frac{\sigma_{s\star
h}}{\sigma_{e}}$, where $\sigma_{s\star h}$ and $\sigma_{e}$ are standard
deviation of BOLD signal and noise, respectively. Note that there is a high
variability of presented SNR definitions, moreover the values of observed SNR
for real data in the literature [28] likely depend on many aspects of data
acquisition, preprocessing and analysis techniques. For our simulations we
used values of SNR in the range $[0.8,3.5]$. This rangewas chosen based on SNR
estimated for the datasets presented in the current paper as well as in [12].
The input signal $s$ was simulated as a vector of the length $N$, that
consists of $K$ non-zero activations (value 1) and $N-K$ zero values.
Positions of non-zero elements are chosen uniformly in the range
$[1,\dots,N]$. Fig. 2 shows an example of a simulated signal and the obtained
estimate. The value of the selection criteria depends on the length of the
signal $N$ and number of non-zero elements $K$. Thus, we explore the estimates
for several different noise levels and multiple numbers of true non-zero
coefficients. To assess the model quality we use the Jaccard index,
sensitivity, and specificity.

Figure 2: Example of a simulated and estimated signal. The simulated neuronal
signal is shown in red, the simulated BOLD signal is shown in blue, the
estimated neuronal signal is shown with thick black lines.

###  5.1 Numerical results

In this section, we present the numerical results of the simulation study.
First, our simulations show that DS/LASSO with AIC/BIC tends to mark almost
all times as activated, in other words, extremely over-fits the data. This
observation is illustrated by the figures in the Supplementary material,
Section D. Thus, the rest of this section presents only the results for DS
AICc, MCI DS, and Bayes classifier based on GMM components (BC GMM).

Fig. 3 shows the averaged results for 100 simulations for varying noise levels
and activation density. Clearly, BC GMM works well only for very sparse
signals, and moreover tends to over-estimate the signal for lower SNR.

Figure 3: Comparison of the performance of the deconvolution methods using
simulated data in terms of the number of detected events, Jaccard index (JI),
sensitivity (SNS), and specificity (SPC). Time series length is $N=300$, the
number of non-zero peaks is shown on the $x$-axis, the noise level on the
$y$-axis. Estimation methods are shown in different rows: DS AICc, MCI DS and
BC GMM.

The other two methods show generally better results. DS AICc tends to slightly
over-estimate the number of detected events for higher SNR and less sparse
signals, whereas it severely underestimates for low SNR. Otherwise, the method
performs relatively well and gives stable results for the presented
parameters. The Jaccard index is typically about 0.6, which is much better
than the results of BIC or uncorrected AIC.

The MCI DS shows generally the best results as illustrated by the almost
uniformly dominant Jaccard index, with balanced performance in terms of
sensitivity and specificity. Indeed it does fail only in the case of extremely
sparse signal (simulated by the case of a single activation from 300 points) –
this happens due to inability to estimate the distribution of true activations
using only one point.

##  6 Real data examples

For validation of the proposed deconvolution approach on real data, we use two
different datasets with distinct experimental designs and data acquisition
parameters. In particular, we use an internal dataset of subjects performing a
memory task and a publicly available dataset originating from the Human
Connectome Project involving motor task fMRI data.

###  6.1 Working memory task

The first used dataset has been collected in the National Institute of Mental
Health (Czech Republic). Only the measured fMRI BOLD signals are used for the
deconvolution procedure, and information about subjects’ responses are used
for quality assessment of MCI in comparison with LASSO/DS.

####  6.1.1 Working memory task: Subjects

The real data example includes the data of 56 subjects (24 males and 32
females, age = $33.66\pm 10.97$) that performed fMRI recording with n-back
task as a part of a complex study with repeated psychiatric, neuroanatomical
and neuropsychological evaluation. All subjects have been recruited for the
study based on the following exclusion criteria: history of mental disorder,
history of neurological disorder, presence of any artificial objects that
would interfere with the magnetic resonance imaging. The study was conducted
in accordance with the Declaration of Helsinki and approved by the local
Ethics Committee of the National Institute of Mental Health (Czech republic).
All participants were informed about the purpose of the study, the
experimental procedures, as well as the fact that they could withdraw from the
study at any time, and provided written informed consent prior to their
participation. The data set of 89 MRI sessions in total was analyzed,
including data from repeated sessions (6 months apart) in 25 subjects (2
visits for 17 subjects and 3 visits for 8 subjects).

####  6.1.2 Working memory task: Experimental design

The n-back task is a continuous performance task commonly used for assessment
of working memory performance.The subject is presented with a sequence of
stimuli, and required to indicate when the current stimulus matches the one
from n steps earlier in the sequence. The individual button press timings are
used as a proxy for true activations of the motor cortex for individual
analysis. As all subjects had the same experimental design, we also carry out
a group level analysis, setting as group ground-truth as button presses
present in at least 80% of subjects (71 of 89) and averaging the measured BOLD
signal in order to decrease noise. On average the subjects responded by button
press 44.63 times (SD = 6.18, min = 27, max = 58). Most of the responses were
correct, on average 5 of the responses were incorrect.

####  6.1.3 Working memory task: fMRI acquisition and preprocessing

Brain images were obtained using Siemens Prisma 3T MR machine. Functional
T2*-weighted images with BOLD contrast were acquired with voxel size 3x3x3 mm,
slice dimensions 64x64 voxels, 32 axial slices, repetition time 2000 ms, echo
time 30 ms, flip angle $70^{\circ}$. Then, the standard fMRI preprocessing was
done using the SPM8 toolbox with spatial realignment, slice-time correction,
normalization to the standard anatomical space, and spatial smoothing. The
primary motor cortex time series was extracted by independent component
analysis using the GIFT toolbox.

####  6.1.4 Working memory task: Results

By comparing the measured signal with the theoretical one, we found out that
the data is in line with the HRF delayed by 1 TR (shifted HRF), rather than
the default HRF, see Fig. 4. We further compared our methods MCI DS and MCI
LASSO with the standard approaches DS and LASSO with the selection criteria
AIC, AICc, BIC.

Figure 4: Selection of appropriate hemodynamic response function (HRF) using a
forward model on real data from the working memory task. Left: the default and
shifted HRF. Comparison of the averaged measured BOLD signal with the
convolution of the true activations vector and the original HRF (middle) and
shifted HRF (right).

Numerical results for the combined (group) signal are shown in Table 1, and
the corresponding receiver operating characteristic (ROC) curve is shown in
Fig. 5. For DS BIC, the estimate is almost a trivial vector. It contains only
true activations, however only a fraction of them. Conversely, the DS AIC-
estimate includes more true activations (high sensitivity), but also contains
many false detections. AICc suffers from mediocre sensitivity. MCI DS then
shows the highest values of all three quality characteristics. LASSO shows a
similar picture.

Fig. 6 shows the quality of estimates of each subject’s time series. The black
and green points that mark DS AIC and DS BIC methods form similarly placed
clouds and demonstrate that almost all estimates contain too many positive
detections. The blue crosses mark the DS AICc and show a lot of estimates with
too few true positive detections. The rest of the blue crosses forms a cloud
that has an intersection with all other methods. Red color is used for MCI DS
estimates, and this cloud lies away from the DS AIC and BIC and demonstrates
the most appropriate estimates in terms of true and false positive rates. See
table 2 for summary results.

Figure 5:  ROC curve of the Dantzig Selector and LASSO solutions for the
group-level real data signal from the working memory task. The blue line of
points corresponds to the set of solutions obtained by LASSO/DS as function of
its one free parameter. The other colored points specify the particular
solution selected by the selection procedures (in case of AIC, AICc, BIC, DS),
or the newly proposed method (MCI). Table 1: Summary group-level performance:
working memory task.

| JI | SNS | SPC
---|---|---|---
DS BIC | 0.07 | 0.07 | 1.00
DS AIC | 0.35 | 0.87 | 0.41
DS AICc | 0.31 | 0.46 | 0.81
MCI DS | 0.51 | 0.61 | 0.92
LASSO BIC | 0.00 | 0.00 | 1.00
LASSO AIC | 0.42 | 0.94 | 0.45
LASSO AICc | 0.11 | 0.12 | 0.96
MCI LASSO | 0.49 | 0.78 | 0.73

JI = Jaccard index, SNS = sensitivity, SPC = specificity. The novel MCI method
results in bold.

Table 2: Summary subject-level performance: working memory task.

| JI | SNS | SPC
---|---|---|---
DS BIC | 0.14(0.12) | 0.37(0.38) | 0.72(0.38)
DS AIC | 0.30(0.04) | 0.86(0.12) | 0.30(0.11)
DS AICc | 0.14(0.11) | 0.26(0.25) | 0.86(0.24)
MCI DS | 0.33(0.05) | 0.56(0.12) | 0.72(0.10)
LASSO BIC | 0.12(0.16) | 0.32(0.42) | 0.76(0.32)
LASSO AIC | 0.30(0.05) | 0.87(0.07) | 0.28(0.11)
LASSO AICc | 0.12(0.15) | 0.19(0.30) | 0.91(0.16)
MCI LASSO | 0.34(0.08) | 0.70(0.10) | 0.58(0.11)

JI = Jaccard index, SNS = sensitivity, SPC = specificity. Mean (std) across
subjects is shown for each method. The novel MCI method results in bold.

Figure 6: Comparison of deconvolution methods performance on individual-
subject level - results from the working memory task. Left: results for the
Dantzig Selector, right: results for LASSO.

###  6.2 Motor task

As a second dataset for our analysis, we use the task-related data from the
Human Connectome Project. We choose this dataset as it is publicly available,
and allows testing the role of specific artifacts due to the acquired
physiological data.

####  6.2.1 Motor task: Subjects

We analyze preprocessed task-related fMRI data of 62 subjects from the Human
Connectome Project [29], https://humanconnectome.org/study/hcp-young-adult.
One of the goals of this analysis is to study the influence of the
physiological processes (breathing and heartbeat) on the deconvolution
quality. Thus, from the 100 downloaded subjects (Motor task), we choose only
those for whom corresponding physio-data were available. The analysis is
performed on the averaged signals from primary motor cortex (Broadman area 4)
for both hemispheres (the LR-modality). The original sample rate of the data
is $0.72s$, however, we subsample the data to $TR=2.16s$ for comparability
with our previous dataset.

####  6.2.2 Motor task: Experimental design

For our analysis, we choose the motor task. Volunteers were asked to tap their
left or right fingers, or squeeze their left or right toes, or move their
tongue to map motor areas. Each block included 10 movements (12 seconds). In
total there were 13 blocks, with 4 of hand movements (2 right and 2 left), 4
of foot movements (2 right and 2 left), and 2 of tongue movements. In
addition, there were 3 15-second fixation blocks. Detailed description of the
experimental paradigm can be found in [30].

####  6.2.3 Motor task: Physiological data

It is widely known that physiological fluctuations resulting from the
heartbeat and respiration are a substantial source of noise in fMRI, as it
perturbs BOLD by fluctuations of non-neuronal origin ([31], [32], [33], [34]).
This causes challenges in the analysis, in particular, false-positive results
may emerge if non-neuronal physiological fluctuations correlate with neuronal
activation.

A commonly used physiological noise correction techniques for fMRI data
include model-based RETROspective Image CORrection (RETROICOR), using
respiratory volume per time and heart rate variability responses (RVT/HRV)
[35]. In our analysis, we use the mentioned techniques implemented in the
PhysIO toolbox, a part of the TNU ”TAPAS” software suite [36]. PhysIO is used
for obtaining the matrix of physiological regressors, which is later used for
cleaning the BOLD signal from breath/heartbeat influence for every subject.

####  6.2.4 Motor task: Analysis and Results

For every subject, we extract two time-series corresponding to the primary
motor cortex (left and right hemisphere) and a matrix of regressors that
consists of physiological regressors obtained by PhysIO, and the motion
regressors as provided for HCP data. The mentioned BOLD signal, obtained after
standard HCP preprocessing, we refer to as the _original_ signal; also, we run
the linear regression with described physiological regressors, to get the
_cleaned_ BOLD signal. In this section we present only the results obtained by
Dantzig Selector on the right hemisphere time series.

Figure 7: Deconvolution results for the motor task data (right hemisphere)
from the HCP dataset, focusing on the influence of physiological artifacts on
BOLD. Top panel: results of deconvolution of original (left) and cleaned
(right) signal. Four different methods of deconvolution are marked in
different colors. Middle and bottom panels: Jaccard Index for original and
cleaned signal, DS MCI, and DS AIC respectively.

We use blocks of movements as known stimuli. Note that we do not distinguish
activations caused by hand and foot or tongue. We deconvolve both types of
signal, original and cleaned, for both hemispheres with all discussed methods;
subject-level results can be found in Fig. 7. Similarly to the working memory
task, DS BIC and AICc show robust underestimate (basically, both of these
criteria choose almost trivial solution), AIC tends either to choose a trivial
solution or overestimate it by choosing almost all the peaks as activation.
MCI DS solutions, marked by red stars, perform relatively reasonable in terms
of both, TPR and FPR.

As another check of how much the cleaning of the signal from physiological and
motion artifacts influences deconvolution results, we compare distribution of
JI before and after cleaning. Paired $t$-test is used for comparison, showing
the significant difference only for DS AIC right hemisphere ($p<0.001$);
similar results are obtained for the LASSO approach instead of DS.

Additionally to the subject-level analysis, we do also group analysis. On a
group level, we deconvolve averaged original signal and averaged cleaned
signal. While the standard criteria of AIC, BIC and AICc fail to choose the
appropriate solution, the MCI approach shows a relatively good Jaccard Index
($0.49$ for original signal and $0.48$ for cleaned data). The method clearly
does not work perfectly, but performs much better than standard ones. Our
results, obtained on the subject- and group-level, suggest that for reasonably
good quality data with appropriate preprocessing the results are not
dominantly affected by false-positives activations due to physiological
artifacts.

Finally, we compare JI for MCI DS with SNR for original and cleaned signals.
We see a highly significant correlation of the Jaccard Index with SNR for
both, original and cleaned versions of BOLD ($r=0.41$, $p=0.0014$ and
$r=0.53$, $p=0.0025$ respectively). This is an additional confirmation of
earlier stated conclusion: deconvolution works better on good signals. In
general, results obtained on the HCP dataset are in line with the previous
ones from the working memory dataset, despite experimental design does not
correspond to exactly sparse neuronal activations, but rather continuous
blocks of activity. The proposed MCI approach shows more relevant
deconvolution results, whereas the typical behavior of classical selection
criteria AIC/BIC/AICc is either to overestimate or underestimate the signal
for Dantzig Selector/LASSO.

##  7 Discussion and conclusions

We have introduced a novel method for recovering a sparse true signal from
smoothed and noisy observations, exemplified by estimation of neuronal
activity from BOLD signal. The main idea of the approach is to overcome the
challenges of regularization combined with standard model selection criteria,
which is achieved by aggregating the information contained in the solutions,
instead of selecting a single one favoured by the criterion. Practically this
entails applying additional operations on the result of standard
regularization-based deconvolution methods, such as LASSO and Dantzig
Selector, and thus builds upon previous achievements in the field.

The advantages of this method over the application of the Bayes information
criterion and Akaike’s information criterion (and its improved variant) were
demonstrated both in numerical simulations as well as in experimental data,
where it provided more accurate estimates of timing of motor responses
occurring during simple motor or more complex cognitive tasks. This provides a
proof of principle that the method would improve the detection of timing of
neuronal activation in other experimental situations, where the brain activity
timing is not directly controlled by the experimental procedure.

Notably, the general Mixture Components Inference method was exemplified here
by application to a particular neuroimaging task – the estimation of the
underlying neuronal activity from fMRI BOLD signal. The BOLD signal character
poses several main challenges to such deconvolution task: its low sampling
rate compared to the scale of neuronal activity, its indirect and postponed
dependence on neuronal signal (hemodynamic response), and complex nature of
sources driving it. Low sampling rate is in principle a major limitation
compared to for instance electrophysiological measurements. However due to the
sluggishness of the hemodynamic response, the BOLD signal anyway corresponds
to a smoothed version of the otherwise fast-scale neuronal activations. BOLD
signal thus corresponds to relatively long-scale weighted temporal averages
(compared to the neuronal firing timescales), justifying thus using relatively
low temporal sampling in both the measurements and models. It is customary to
impose for convenience the same temporal discretization for the estimated
neural activity as for the BOLD measurements, and multiple true activations
within a single TR would thus be, even in the best scenario of perfect
estimate, be interpreted as a single joint activation in the corresponding
time point. Of course, in principle the model could be extended to include
finer temporal sampling of the estimated neural activity, however, given the
smoothness of the signal and unknown precise hemodynamic response, it might
remain challenging to provide reliable estimates with a sub-TR temporal
precision.

Future work may involve combining the presented approach with some other
directions recently taken up in the field, namely spatial regularization,
change sparsity, HRF estimation and further processing. In more detail, the
possible extensions and integrations include: the joint estimation of the HRF
and neuronal activity, such as in the fused LASSO and Tikhonov regularization
of the HRF [37], spatial structure modelling including regularization as in
the work of [38], which was further extended in [39] to remove the need for a
priori brain parcelation, or alternative penalty schemes in regularization
such as combining total variation and negativity penalties as in [40]. Of
interest could be the extension to application of such deconvolution using
more complex data models, such as the bilinear model including neuronal
activity autocorrelation, stimulation and modulation terms [41, 42],
potentially in the context of network analysis.

As the outlined approach relies on being combined with an initial method
providing a parameterized family of solutions (in this paper, we used the
LASSO and Danzig Selector to provide examples of such families), we leave open
the question of how it relates and could synergize with other families of
approaches such as using the cubature Kalman filtering [43], activelets
(wavelets designed using prior knowledge on the hemodynamic response function
form [44]), spatiotemporal HRF deconvolution [45], multivariate semi-blind
deconvolution [46], improved sparse paradigm free mapping with no selection of
the regularization parameter [47] or nonparametric hemodynamic deconvolution
of fMRI using homomorphic filtering [48]. We can generally speculate that the
key idea of combining model outputs using e.g. Bayesian approaches instead of
crisp model selection within a constrained parametric family may be applicable
to other schemes beyond classical regularization schemes such as LASSO, as
long as these involve selection among parametric model families.

The obtained information on local neural activations can be of interest on its
own such as in the application to study cortical response to unpredictable
mental events as outlined in [49], or possibly used for obtaining meta-
analytic interpretation of the brain activation patterns [50]. Importantly, it
has been proposed that functional connectivity in MRI is driven by spontaneous
BOLD events detectable by deconvolution methods [51]; and the use of
deconvolution as an intermediate step to estimate effective connectivity has
been advocated [52]. Similarly, the dynamics of deconvolution-based brain
activation patterns can be further explicitly modelled to fit sparse coupled
hidden Markov models [53] characterizing the large-scale brain dynamics.

Moreover, deconvolution techniques are used in a range of specific
applications in neuroimaging. In [54], authors are using the deconvolution for
ROI identification. Deconvolution for multi-echo functional MRI is discussed
in [55], and the authors claim that the technique is very promising in
understanding the dynamic nature of brain activity in resting-state and
paradigms with unknown timing of the BOLD events. It is also used for
computation of activity maps for getting the dynamical functional connectivity
in spontaneous flow of thought/resting-state studies [56]. Authors of [57]
mention that even though they use thresholding of the BOLD signal, there is a
place for the application of more sophisticated methods as deconvolution.
Including deconvolution is also mentioned as possible future work for
filtering of the physiological noise [58]. Improving the signal-to-noise ratio
of the data with the deconvolution is pointed out as an option in [59].
Deconvolution has also been used for obtaining the effective connectivity
matrix in an autism study [60].

Last but not least, the development of deconvolution approaches is relevant
not only for processing of most common fMRI BOLD signal but also in the
context of functional near-infrared spectroscopy [61, 62] and quantitative
susceptibility mapping [63]. In general, this field demonstrates the fruitful
interaction between advanced engineering attempts to tackle experimental and
technical challenges specific to each measurement modality or even particular
type of experiment and the development of novel algorithmic and mathematical
approaches, possibly applicable across disciplines.

##  8 Acknowledgment

We thank David Tomeček for sample fMRI data and Jessica Barilone for helping
with language and style editing.

## References

  * [1] S. Ogawa, T.-M. Lee, A. R. Kay, and D. W. Tank, “Brain magnetic-resonance-imaging with contrast dependent on blood oxygenation,” _PNAS_ , vol. 87, no. 24, pp. 9868–9872, 1990.
  * [2] S. Ogawa _et al._ , “Intrinsic signal changes accompanying sensory stimulation - functional brain mapping with magnetic-resonance-imaging,” _PNAS_ , vol. 89, no. 13, pp. 5951–5955, 1992.
  * [3] M. Bianciardi _et al._ , “Sources of functional magnetic resonance imaging signal fluctuations in the human brain at rest: a 7 T study,” _Magn Reson Imaging_ , vol. 27, no. 8, pp. 1019–1029, 2009.
  * [4] C. C. Gaudes, N. Petridou, S. T. Francis, I. L. Dryden, and P. A. Gowland, “Paradigm free mapping with sparse regression automatically detects single-trial functional magnetic resonance imaging blood oxygenation level dependent responses,” _Hum Brain Mapp_ , vol. 34, no. 3, pp. 501–518, 2013\.
  * [5] R. Tibshirani, “Regression shrinkage and selection via the LASSO,” _J R Stat Soc Series B Stat Methodol_ , vol. 58, no. 1, pp. 267–288, 1996.
  * [6] E. Candes and T. Tao, “The Dantzig selector: statistical estimation when p is much larger than n,” _Ann Stat_ , vol. 35, no. 6, pp. 2313–2351, 2007.
  * [7] H. Akaike, “A New Look at the Statistical Model Identification,” _IEEE Trans Automat Contr_ , vol. AC-19, no. 6, pp. 716–723, 1974.
  * [8] G. Schwarz, “Estimating the Dimension of a Model,” _Ann Stat_ , vol. 6, no. 2, pp. 461–464, 1978.
  * [9] K. P. Burnham and D. R. Anderson, “Multimodel inference: understanding AIC and BIC in Model Selection,” _Sociol Methods Res_ , vol. 33, pp. 261–304, 2004.
  * [10] N. Sugiura, “Further analysis of the data by Akaike’ s information criterion and the finite corrections,” _Commun Stat Theory Methods_ , vol. 7, no. 1, pp. 13–26, 1978.
  * [11] C. M. Hurvich and C.-L. Tsai, “Regression and time series model selection in small samples,” _Biometrika_ , vol. 76, no. 2, pp. 297–307, 1989.
  * [12] A. Pidnebesna, D. Tomeček, and J. Hlinka, “BRAD: Software for BRain Activity Detection from hemodynamic response,” _Comput Methods Programs Biomed_ , vol. 156, pp. 113–119, 2018.
  * [13] R. M. Gray, “Toeplitz and circulant matrices: a review,” _Found Trends Somm Inf Theory_ , vol. 2, no. 3, pp. 155–239, 2006.
  * [14] G. M. James, P. Radchenko, and J. Lv, “DASSO: connections between the Dantzig selector and LASSO.” _J R Stat Soc Series B Stat Methodol_ , vol. 71, no. 1, pp. 127–142, 2009.
  * [15] P. J. Bickel, Y. Ritov, and A. B. Tsybakov, “Simultaneous analysis of Lasso and Dantzig selector,” _Ann Stat_ , vol. 37, no. 4, pp. 1705–1732, 2009\.
  * [16] M. S. Asif and J. Romberg, “Dantzig selector homotopy with dynamic measurements,” in _Comput Imaging VII_ , ser. Proceedings of SPIE, vol. 7246, 2009.
  * [17] H. Zou, T. Hastie, and R. Tibshirani, “On the ”degrees of freedom” of the LASSO,” _Ann Stat_ , vol. 35, no. 5, pp. 2173–2192, 2007.
  * [18] I. Rish, “An empirical study of the naive Bayes classifier,” IBM Research Division Thomas J. Watson Research Center, Tech. Rep., 2001.
  * [19] D. J. Hand and K. Yu, “Idiot’s Bayes - Not So Stupid After All?” _Int Stat Rev_ , vol. 69, no. 3, pp. 385–398, 2001.
  * [20] G. McLachlan and D. Peel, _Finite mixture models_ , 2000.
  * [21] R. Maiboroda and O. Sugakova, “Statistics of mixtures with varying concentrations with application to DNA microarray data analysis,” _J Nonpar Stat_ , vol. 24, no. 1, pp. 201–215, 2012.
  * [22] R. Maiboroda, _Otsiniuvannia ta klasyfikatsiia za sposterezhenniamy iz sumishi. Monohrafiia._ , R. Maiboroda and O. Sugakova, Eds. Publishing and printing center ”Kyjiv University”, 2008\.
  * [23] A. Doronin and R. Maiboroda, “Testing hypotheses on moments by observations from a mixture with varying concentrations,” _Mod Stoch Theor Applic_ , vol. 1, no. 2, pp. 195–209, 2015.
  * [24] O. Sugakova, “Classification of components of a mixture,” _Theory Probab Math.Stat_ , vol. 72, pp. 157–166, 2006.
  * [25] F. Autin and C. Pouet, “Adaptive test on components of densities mixture,” _Math Meth Stat_ , vol. 21, no. 2, pp. 93–108, 2012.
  * [26] W. K. Hardle, M. Muller, S. Sperlich, and A. Werwatz, _Nonpar Semipar Models_. Springer, 2006.
  * [27] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, “Least Angle Regression,” _Ann Stat_ , vol. 32, no. 2, pp. 407–451, 2004.
  * [28] M. Welvaert and Y. Rosseel, “On the Definition of Signal-To-Noise Ratio and Contrast-To-Noise Ratio for fMRI Data,” _PLOS ONE_ , vol. 8, no. 11, pp. 1–10, 2013.
  * [29] D. C. Van Essen, S. M. Smith, D. M. Barch, T. E. J. Behrens, E. Yacoub, and K. Ugurbil, “The WU-Minn Human Connectome Project: An overview,” _NeuroImage_ , vol. 80, pp. 62–79, 2013.
  * [30] D. M. Barch _et al._ , “Function in the human connectome: task-fMRI and individual differences in behavior,” _NeuroImage_ , vol. 80, pp. 169–189, 2013.
  * [31] G. Krüger and G. H. Glover, “Physiological noise in oxygenation-sensitive magnetic resonance imaging,” _Magn Reson Med_ , vol. 46, no. 4, pp. 631–637, 2001.
  * [32] C. Hutton _et al._ , “The impact of physiological noise correction on fMRI at 7T,” _NeuroImage_ , vol. 57, no. 1, pp. 101–112, 2011.
  * [33] C. Triantafyllou _et al._ , “Comparison of physiological noise at 1.5 T, 3 T and 7 T and optimization of fMRI acquisition parameters,” _NeuroImage_ , vol. 26, no. 1, pp. 243–250, 2005.
  * [34] T. T. Liu, “Noise contributions to the fMRI signal: An overview,” _NeuroImage_ , vol. 143, pp. 141–151, 2016.
  * [35] G. H. Glover, T.-Q. Li, and D. Ress, “Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR,” _Magn Reson Med_ , vol. 44, no. 1, pp. 162–167, 2000.
  * [36] L. Kasper _et al._ , “The PhysIO Toolbox for Modeling Physiological Noise in fMRI Data,” _J Neurosci Methods_ , vol. 276, pp. 56–72, 2017.
  * [37] P. Aggarwal, A. Gupta, and A. Garg, “Joint estimation of activity signal and HRF in fMRI using fused LASSO,” in _2015 IEEE GlobalSIP_ , 2015, pp. 829–833.
  * [38] F. I. Karahanoglu, C. Caballero-Gaudes, F. Lazeyras, and D. Van De Ville, “Total activation: fMRI deconvolution through spatio-temporal regularization,” _NeuroImage_ , vol. 73, pp. 121–134, 2013.
  * [39] Y. Farouj, F. I. Karahanoglu, and D. Van De Ville, “Regularized spatiotemporal deconvolution of fMRI data using gray-matter constrained total variation,” in _2017 IEEE Isbi 2017_ , 2017, pp. 472–475.
  * [40] L. Hernandez-Garcia and M. O. Ulfarsson, “Neuronal event detection in fMRI time series using iterative deconvolution techniques,” _Magn Reson Imaging_ , vol. 29, no. 3, pp. 353–364, 2011.
  * [41] W. Penny, Z. Ghahramani, and K. Friston, “Bilinear dynamical systems,” _Philos Trans R Soc Lond B Biol Sci_ , vol. 360, pp. 983–993, 2005.
  * [42] S. Makni, C. F. Beckmann, S. Smith, and M. Woolrich, “Bayesian deconvolution fMRI data using bilinear dynamical systems,” _NeuroImage_ , vol. 42, no. 4, pp. 1381–1396, 2008.
  * [43] M. Havlicek, J. Jan, M. Brazdil, and V. D. Calhoun, “Dynamic Granger causality based on Kalman filter for evaluation of functional network connectivity in fMRI data,” _NeuroImage_ , vol. 53, no. 1, pp. 65–77, 2010\.
  * [44] I. Khalidov, J. Fadili, F. Lazeyras, D. V. D. Ville, and M. Unser, “Activelets: Wavelets for sparse representation of hemodynamic responses,” _Signal Process_ , vol. 91, no. 12, pp. 2810–2821, 2011.
  * [45] K. M. Aquino, P. A. Robinson, M. M. Schira, and M. Breakspear, “Deconvolution of neural dynamics from fMRI data using a spatiotemporal hemodynamic response function,” _NeuroImage_ , vol. 94, pp. 203–215, 2014.
  * [46] H. Cherkaoui, T. Moreau, A. Halimi, C. Leroy, and P. Ciuciu, “Multivariate semi-blind deconvolution of fMRI time series,” _NeuroImage_ , vol. 241, no. November 2020, p. 118418, 2021.
  * [47] E. Urunuela, S. Moia, and C. Caballero-Gaudes, “A low rank and sparse paradigm free mapping algorithm for deconvolution of FMRI data,” _Proc Intern Symp Biom Imaging_ , vol. 2021-April, pp. 1726–1729, 2021.
  * [48] K. R. Sreenivasan, M. Havlicek, and G. Deshpande, “Nonparametric hemodynamic deconvolution of fMRI using homomorphic filtering,” _IEEE Trans Med Imaging_ , vol. 34, no. 5, pp. 1155–1163, 2015.
  * [49] C. C. Gaudes, N. Petridou, I. L. Dryden, L. Bai, S. T. Francis, and P. A. Gowland, “Detection and characterization of single-trial fMRI BOLD responses: paradigm free mapping,” _Hum Brain Mapp_ , vol. 32, no. 9, pp. 1400–1418, 2011.
  * [50] F. M. Tan _et al._ , “Decoding fMRI events in sensorimotor motor network using sparse paradigm free mapping and activation likelihood estimates,” _Hum Brain Mapp_ , vol. 38, no. 11, pp. 5778–5794, 2017.
  * [51] T. W. Allan _et al._ , “Functional Connectivity in MRI Is Driven by Spontaneous BOLD Events,” _PLOS ONE_ , vol. 10, no. 4, 2015.
  * [52] G.-R. Wu, W. Liao, S. Stramaglia, J.-R. Ding, H. Chen, and D. Marinazzo, “A blind deconvolution approach to recover effective connectivity brain networks from resting state fMRI data,” _Med Image An_ , vol. 17, no. 3, pp. 365–374, 2013.
  * [53] T. A. W. Bolton, A. Tarun, V. Sterpenich, S. Schwartz, and D. Van De Ville, “Interaction Between Large-Scale Functional Brain Networks are Captured by Sparse Coupled HMMs,” _IEEE Trans Med Imaging_ , vol. 37, no. 1, pp. 230–240, 2018.
  * [54] A. Veloz _et al._ , “Ant Colony Clustering for ROI Identification in Functional Magnetic Resonance Imaging,” _Comput Intell Neurosci_ , vol. 2019, 2019.
  * [55] C. Caballero-Gaudes, S. Moia, P. Panwar, P. A. Bandettini, and J. Gonzalez-Castillo, “A deconvolution algorithm for multi-echo functional MRI: Multi-echo Sparse Paradigm Free Mapping,” _NeuroImage_ , vol. 202, p. 116081, 2019.
  * [56] J. Gonzalez-Castillo, C. Caballero-Gaudes, N. Topolski, D. A. Handwerker, F. Pereira, and P. A. Bandettini, “Imaging the spontaneous flow of thought: Distinct periods of cognition contribute to dynamic functional connectivity during rest,” _NeuroImage_ , vol. 202, p. 116129, 2019.
  * [57] G. Deco, E. Tagliazucchi, H. Laufs, A. Sanjuán, and M. L. Kringelbach, “Novel Intrinsic Ignition Method Measuring Local-Global Integration Characterizes Wakefulness and Deep Sleep,” _eNeuro_ , vol. 4, no. 5, 2017\.
  * [58] Särkkä S. _et al._ , “Dynamic retrospective filtering of physiological noise in BOLD fMRI: DRIFTER.” _NeuroImage_ , vol. 60, no. 2, pp. 1517–1527, 2012.
  * [59] F. De Vico Fallani _et al._ , “Graph analysis of functional brain networks: practical issues in translational neuroscience.” _Philosophical transactions of the Royal Society of London. Series B, Biological sciences_ , vol. 369, no. 1653, p. 20130521, 2014.
  * [60] G. Deshpande, L. Libero, K. Sreenivasan, H. Deshpande, and R. Kana, “Identification of neural connectivity signatures of autism using machine learning,” _Front Hum Neurosci_ , vol. 7, p. 670, 2013.
  * [61] H. Santosa, F. Fishburn, X. Zhai, and T. J. Huppert, “Investigation of the sensitivity-specificity of canonical- and deconvolution-based linear models in evoked functional near-infrared spectroscopy,” _Neurophotonics_ , vol. 6, no. 2, pp. 1–10, 2019.
  * [62] A. Seghouane and D. Ferrari, “Robust Hemodynamic Response Function Estimation From fNIRS Signals,” _IEEE Trans Signal Process_ , vol. 67, no. 7, pp. 1838–1848, 2019.
  * [63] M. Costagli _et al._ , “Quantitative Susceptibility Mapping of Brain Function During Auditory Stimulation,” _IEEE Trans Radiat Plasma Med Sci_ , vol. 3, no. 4, pp. 516–522, 2019.