Hamza Cherkaoui

Thomas Moreau

Abderrahim Halimi

Claire Leroy

Philippe Ciuciu
email: philippe.ciuciu@cea.fr

Multivariate semi-blind deconvolution of fMRI time series

Keywords: BOLD signal, HRF, sparsity, low-rank decomposition, multivariate modeling, dictionary learning, UK Biobank

Abstract

Whole brain estimation of the haemodynamic response function (HRF) in functional magnetic resonance imaging (fMRI) is critical to get insight on the global status of the neurovascular coupling of an individual in healthy or pathological condition. Most of existing approaches in the literature works on task-fMRI data and relies on the experimental paradigm as a surrogate of neural activity, hence remaining inoperative on resting-stage fMRI (rs-fMRI) data. To cope with this issue, recent works have performed either a two-step analysis to detect large neural events and then characterize the HRF shape or a joint estimation of both the neural and haemodynamic components in an univariate fashion. In this work, we express the neural activity signals as a combination of piece-wise constant temporal atoms associated with sparse spatial maps and introduce an haemodynamic parcellation of the brain featuring a temporally dilated version of a given HRF model in each parcel with unknown dilation parameters. We formulate the joint estimation of the HRF shapes and spatio-temporal neural representations as a multivariate semi-blind deconvolution problem in a paradigm-free setting and introduce constraints inspired from the dictionary learning literature to ease its identifiability. A fast alternating minimization algorithm, along with its efficient implementation, is proposed and validated on both synthetic and real rs-fMRI data at the subject level. To demonstrate its significance at the population level, we apply this new framework to the UK Biobank data set, first for the discrimination of haemodynamic territories between balanced groups (n = 24 individuals in each) patients with an history of stroke and healthy controls and second, for the analysis of normal aging on the neurovascular coupling. Overall, we statistically demonstrate that a pathology like stroke or a condition like normal brain aging induce longer haemodynamic delays in certain brain areas (e.g. Willis polygon, occipital, temporal and frontal cortices) and that this haemodynamic feature may be predictive 1. Introduction

Context

Functional magnetic resonance imaging (fMRI) non-invasively records brain activity by dynamically measuring the blood oxygenation level-dependent (BOLD) contrast. The latter reflects the local changes in the deoxyhemoglobin concentration in the brain [START_REF] Ogawa | Intrinsic signal changes accompanying sensory stimulation: functional brain mapping with magnetic resonance imaging[END_REF] and thus indirectly measures neural activity through the neurovascular coupling. This coupling is usually characterized as a linear and time-invariant system and thus summarized by its impulse response, the so-called haemodynamic response function (HRF) [START_REF] Bandettini | Processing strategies for time-course data sets in functional MRI of the human brain[END_REF][START_REF] Boynton | Linear systems anlaysis of functional magnetic resonance imaging in human V1[END_REF]. The estimation of the response is of a primary interest: a change in the haemodynamic response could be linked to the pharmacological mechanism of a drug [START_REF] Do | Hemodynamic modeling of aspirin effects on bold responses at 7t[END_REF], the effect of healthy aging [START_REF] West | BOLD hemodynamic response function changes significantly with healthy aging[END_REF] or the consequence of a neuropathological process, for example Alzheimer's disease [START_REF] Asemani | Effects of ageing and alzheimer disease on haemodynamic response function: a challenge for event-related fMRI[END_REF]. Thus, the HRF could be considered as a precious bio-marker to investigate the neurovascular function of the brain in a healthy or pathological condition. Moreover, its estimation also links the observed BOLD signal to the underlying neural activity, which can in turn be used to better understand cognitive processes in the healthy brain or to uncover functional alteration in pathological condition.

Related works

Several methods have been designed to estimate this haemodynamic response in the case of task-related fMRI (tfMRI). In this setup, the participant is engaged in an experimental paradigm (EP) during the imaging session, which alternates between rest and taskperiods [START_REF] Friston | Eventrelated fMRI: characterizing differential responses[END_REF][START_REF] Ciuciu | Unsupervised robust nonparametric estimation of the hemodynamic response function for any fMRI experiment[END_REF][START_REF] Lindquist | Validity and power in hemodynamic response modeling: a comparison study and a new approach[END_REF][START_REF] Pedregosa | Data-driven HRF estimation for encoding and decoding models[END_REF]. Commonly, supervised HRF estimation methods fit a model to explain the observed BOLD signal from the EP [START_REF] Goutte | Modeling the haemodynamic response in fMRI using smooth FIR filters[END_REF][START_REF] Ciuciu | Unsupervised robust nonparametric estimation of the hemodynamic response function for any fMRI experiment[END_REF][START_REF] Lindquist | Validity and power in hemodynamic response modeling: a comparison study and a new approach[END_REF][START_REF] Vincent | Spatially adaptive mixture modeling for analysis of fMRI time series[END_REF][START_REF] Chaari | Hemodynamic-informed parcellation of fMRI data in a joint detection estimation framework[END_REF][START_REF] Degras | A hierarchical model for simultaneous detection and estimation in multi-subject fMRI studies[END_REF][START_REF] Rosa | On the distinguishability of hrf models in fMRI[END_REF][START_REF] Pedregosa | Data-driven HRF estimation for encoding and decoding models[END_REF][START_REF] Eickenberg | Gaussian processes for hrf estimation for bold fmri[END_REF]. A limitation of these approaches is that the EP is used as a surrogate for the neural activity. Therefore they do not account for possible delays in the subject's responses compared to the task onsets, thus yielding a biased HRF estimate. Moreover, these methods cannot be used on resting-state fMRI data (rs-fMRI), where the participant is laying still in the MRI scanner and where no EP is available to serve as surrogate for neural activity. A potential solution would consist in uncovering neural activity by thresholding large values in the BOLD signal and then constructing a sparse binary sequence [START_REF] Wu | A blind deconvolution approach to recover effective connectivity brain networks from resting state fMRI data[END_REF] as an artificial EP to be subsequently used in supervised HRF estimation techniques. However, due to this two-step procedure, any error in neural activity detection tends to bias the recovery of HRF shape, especially in regard to the haemodynamic delay.

On the other hand, a long-standing literature on fMRI deconvolution methods has emerged since the late 90s to uncover the underlying activity-inducing signal at the fMRI timescale of seconds, see [START_REF] Glover | Deconvolution of impulse response in event-related bold fMRI[END_REF] for Wiener filtering and smooth estimation of activity-inducing signals and [START_REF] Gitelman | Modeling regional and psychophysiologic interactions in fMRI: the importance of hemodynamic deconvolution[END_REF] for its generalization to study psychophysiologic interactions at the neuronal level. Importantly, within a temporal fMRI deconvolution framework, the neural activity signal was estimated by involving the 2 -norm regularization (Caballero-Gaudes et al., 2011). Interestingly, most recent works have relied on sparse regularization to recover a limited number of spike-like activations directly in the time domain [START_REF] Hernandez-Garcia | Neuronal event detection in fMRI time series using iterative deconvolution techniques[END_REF]Caballero-Gaudes et al., 2011, 2012[START_REF] Watanabe | A widely applicable Bayesian information criterion[END_REF][START_REF] England | Briefing document: First incidence of stroke estimates for england 2007 to 2016[END_REF] or imposing sparsity in the activelet-domain, which is a wavelet basis that is tailored to the haemodynamic properties [START_REF] Khalidov | Activelets and sparsity: a new way to detect brain activation from fmri data[END_REF][START_REF] Khalidov | Activelets: Wavelets for sparse representation of hemodynamic responses[END_REF].

Next, a foundational work, called Total Activation [START_REF] Karahanoglu | Total activation: fMRI deconvolution through spatio-temporal regularization[END_REF], has proposed a spatio-temporal model of the underlying activity-inducing signal including both temporal and spatial sparsity-based regularization. In doing so, the recovered neural activity profiles are used to define functional networks from deconvolved BOLD signals. Such approaches have thus revealed the existence of neurally-related events in the BOLD signal during periods of rest [START_REF] Petridou | Periods of rest in fMRI contain individual spontaneous events which are related to slowly fluctuating spontaneous activity[END_REF][START_REF] Karahanoglu | Transient brain activity disentangles fmri resting-state dynamics in terms of spatially and temporally overlapping networks[END_REF][START_REF] Gonzalez-Castillo | Imaging the spontaneous flow of thought: Distinct periods of cognition contribute to dynamic functional connectivity during rest[END_REF]. However, Total Activation [START_REF] Karahanoglu | Total activation: fMRI deconvolution through spatio-temporal regularization[END_REF] suffers from a main limitations:

The HRF shape in this model was fixed and kept constant throughout the brain, making this approach primarily suitable for normal populations. Without disentangling the neural and vascular components, but taking into consideration the spatial variability in the neurovascular coupling, [START_REF] Caballero-Gaudes | Structured sparse deconvolution for paradigm free mapping of functional MRI data[END_REF] have proposed a voxelwise HRF modeling.

Alternatively some recent works have suggested either a separate or concurrent estimation of the neural activity and the HRF profile [START_REF] Wu | A blind deconvolution approach to recover effective connectivity brain networks from resting state fMRI data[END_REF]Cherkaoui et al., 2019;[START_REF] Farouj | Bold signal deconvolution under uncertain haemodynamics: A semi-blind approach[END_REF]. These approaches are often referred to as semiblind deconvolution schemes of the BOLD signal as they do not fully estimate an unknown HRF shape from scratch but instead try to uncover a transformation from a pre-existing HRF profile such as the canonical shape. Both (Cherkaoui et al., 2019) and [START_REF] Farouj | Bold signal deconvolution under uncertain haemodynamics: A semi-blind approach[END_REF] rely on the hypothesis of a block-type property for the neural activation signal as initially proposed in [START_REF] Karahanoglu | Total activation: fMRI deconvolution through spatio-temporal regularization[END_REF] but are able to infer the haemodynamic profile for instance with a limited parameterization to deal with the magnitude and delay ambiguities between the neural input the haemodynamic filter.

Goals and contributions

This paper extends preliminary contributions (Cherkaoui et al., 2019) and offers a new algorithm that aims to fit a rich multivariate decomposition of the BOLD data using a semi-blind deconvolution and low-rank sparse decomposition.

The model distinguishes two major parts in the BOLD signal: the neurovascular coupling and the neural activity signal. One of its main features is to tackle the intra-subject haemodynamic variability by introducing regionwise HRF profiles over a brain parcellation, encoding the spatial variations of the neurovascular coupling. Also, in regard to neural activity, we follow ideas developed in the convolutional dictionary learning literature to develop a low-rank approximation of a signal [START_REF] Grosse | Shift-invariant sparse coding for audio classification[END_REF][START_REF] Dupré La Tour | Multivariate convolutional sparse coding for electromagnetic brain signals[END_REF]. We thus model the neural activation signals as a combination of a limited number of piece-wise constant temporal profiles associated with spatial maps that capture the magnitude of functional networks. As such, the neural input signals are represented with their own spatial representation that departs from the neurovascular parcellation. We of course present a scalable optimization algorithm that is able to fit all parameters of interest in the whole brain in a reasonable computing time.

The rest of the paper is organized as follows. Section 2 introduces our modeling of the BOLD data and presents our semi-blind blind deconvolution algorithm.

Next, our technique is validated on simulated data and on real rs-fMRI data at the individual level in Section 3. In Section 4 we illustrate the proposed framework at the population level on the large UK Biobank database. Two main applications were targeted, namely the neurovascular discrimination of patients with stroke episodes as compared to healthy controls and the prediction of brain age. Section 5 discusses the potential impact of our method as well as its limitations. Conclusion and outlook are synthesized in Section 6.

Multivariate low-rank decomposition of the BOLD signal

In this section, we present our modeling of the BOLD signal and derive an efficient algorithm to estimate its parameters.

Notation

In what follows, y i denotes the i th entry in vector y.

Let T = T -L + 1, the convolution of two signals a ∈ R 1× T and v ∈ R 1×L is denoted by v * a ∈ R 1×T . For A ∈ R P × T , v * A ∈ R P ×T is obtained by convolving each line of A with v.
For V ∈ R P ×L , V * A ∈ R P ×T is obtained by convolving each line of A with the corresponding line of V . Let ∇ be the first-order difference operator such that et al., 1996). This model is illustrated in Fig. 1. For each voxel, the measured time series, denoted by y j ∈ R 1×T , is the convolution of a neural activation signal a j ∈ R 1× T , with a given HRF, v ∈ R 1×L and e j ∈ R 1×T refers to an additive white Gaussian noise, which leads to:

∀x ∈ R T , ∇x ∈ R T -1 with (∇x) i = x i -x i-1 , ∀i ∈ {2, . . . , T }.
y j = v * a j + e j . (1)
Although the noise that contaminates the BOLD effect is serially correlated in time [START_REF] Woolrich | Temporal autocorrelation in univariate linear modeling of fMRI data[END_REF], here for the sake of simplicity we assume that the noise corrupting the fMRI time series is white and Gaussian. As shown hereafter in Section 2.1.2, this assumption is tenable as our model will be flexible enough to capture and segregate noise-related components. However, for the sake of completeness, additional justification can be found in the supplementary material (cf.

Section A.8). In the latter, we studied on simulations the impact of a first-order autoregressive Gaussian noise on the estimation accuracy and showed that for SNR ( -8 dB) and autoregressive parameter ( 0.3) values met on real fMRI data, the reconstruction error due to the model mismatch remains quite low.

Typically, the HRF v has a restricted support in time of about 20 s. The challenge with HRF modeling is to find a fair trade-off between a flexible model that is able to capture the true haemodynamic response in each brain area and a reliable one that limits overfitting by reducing the number of degrees of freedom (diminishing the variance). Since our approach estimates the neural activity along with the HRF, reducing the number of degrees of freedom is critical to avoid the aforementioned overfitting. In this paper, we assume M different HRF with v m being the HRF corresponding to the m th region Θ m . Numerous approaches have

Figure 2: Illustration of two haemodynamic response functions (HRF) denoted h 1 , h 2 with the full-width-at-half-maximum (FWHM) and the time-to-peak (TTP) (of h 1 ) depicted. Here we illustrate these two HRFs with different dilation parameters (δ 1 , δ 2 ) such as δ 2 ≤ δ 1 , leading to TTP(h 1 ) < TTP(h 2 ) and FWHM(h 1 ) < FWHM(h 2 ).

been proposed to model this haemodynamic response v m . The Finite Impulse Response (FIR) [START_REF] Dale | Optimal experimental design for event-related fMRI[END_REF][START_REF] Glover | Deconvolution of impulse response in event-related bold fMRI[END_REF] model does not assume any particular shape for the HRF which make it very flexible but prone to data overfitting in the presence of noise. Regularization has thus been introduced to constrain the overall HRF shape in FIR models and limit their tendency to overfitting, see for instance penalization over the second-order derivative to end up with physiologically plausible smooth HRF estimates [START_REF] Ciuciu | Unsupervised robust nonparametric estimation of the hemodynamic response function for any fMRI experiment[END_REF][START_REF] Casanova | The impact of temporal regularization on estimates of the bold hemodynamic response function: a comparative analysis[END_REF]. Alternatively, the HRF has been modeled as a linear decomposition of predefined atoms such as B-splines [START_REF] Zhang | A comparative study of one-level and two-level semiparametric estimation of hemodynamic response function for fMRI data[END_REF][START_REF] Vakorin | Characterizing the functional MRI response using tikhonov regularization[END_REF], wavelets [START_REF] Khalidov | Activelets: Wavelets for sparse representation of hemodynamic responses[END_REF], a sensitivity-selected set [START_REF] Woolrich | Constrained linear basis sets for hrf modelling using variational bayes[END_REF] or more physiologically informed patterns such as the canonical HRF and its derivatives in time and with respect to the dispersion parameter [START_REF] Friston | Eventrelated fMRI: characterizing differential responses[END_REF]. All these methods intend to capture fluctuations in haemodynamic delay or shape with the minimum number of parameters. Last, to constrain even more the parameter values and reduce variance estimates, parametric models such as the inverse logit transform [START_REF] Lindquist | Validity and power in hemodynamic response modeling: a comparison study and a new approach[END_REF] have been proposed and successfully tested when benchmarking over multiple fMRI datasets [START_REF] Lindquist | Modeling the hemodynamic response function in fMRI: efficiency, bias and mis-modeling[END_REF]. In this work, we propose the time dilation HRF model (Cherkaoui et al., 2019), which captures the haemodynamic delay by dilating a reference HRF pattern:

v δ = v ref (δt)
where v ref is a reference temporal profile, here the canonical HRF. This approach is efficient while simple as it encodes delay fluctuations through a single scalar parameter δ (one degree of freedom). One limitation of this choice is that it leads to the simultaneously variation of the full width-at-half-maximum (FWHM) of the HRF and its time-to-peak (TTP) (see also Fig. 2). This coupling between TTP and FWHM partly alleviates time-shift ambiguities between the HRF and the neural activity signal during the blind deconvolution process. Although this model is quite simple, it remains attractive as it provides a low variance HRF estimate while capturing the main effects on the neurovascular coupling as described in Section 4. Our voxelwise model reads as follows:

y j = v δ * a j + e j .
(2)

Multivariate modeling

To better account for the spatial structure of the neurovascular system and the intrinsic organization of functional networks, we extend this model to the multivariate setting. Our multivariate model reads as follows:

Y = v δ * A + E , (3)
where A = ( a j ) P j=1 ∈ R P × T and E = (e j ) P j=1 ∈ R P ×T . One limitation of this straightforward approach is that it constrains the haemodynamic response to be the same across the whole brain. As the HRF shape depends on the neurovascular coupling, its features vary in space over different brain areas and between individuals [START_REF] Handwerker | Variation of bold hemodynamic responses across subjects and brain regions and their effects on statistical analyses[END_REF][START_REF] Badillo | Group-level impacts of within-and betweensubject hemodynamic variability in fMRI[END_REF]. This suggests that, for a given subject, the HRF should be modeled locally in the brain. An appropriate approach for doing so is to rely on existing brain parcellation [START_REF] Varoquaux | Learning and comparing functional connectomes across subjects[END_REF]. Ideally to accurately fit the real haemodynamic response function in a subject, we would favor a large number of regions. However, the larger this number, the smaller the number of voxels per region, which could impair the stability of HRF estimation. In practice, using a minimum of 50 voxels per region leads to stable HRF estimates. In contrast, aggregating too many voxels (e.g. > 10 3 ) results in identical HRFs over the whole brain. Hence, in all our experiments, we used parcellations where each region consisted of a few hundred voxels. We refer the reader to [START_REF] Vincent | Sensitivity analysis of parcellation in the joint detection-estimation of brain activity in fMRI[END_REF] for more information on this point. In this work, we rely on the Harvard-Oxford probabilistic brain atlas [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF]. We threshold the probabilities to obtain a fine brain parcellation that offers enough flexibility to adapt to the true haemodynamic system. In what follows, we mathematically introduce a brain parcellation with M regions with (Θ m ) M m=1 ∈ {0, 1} p 1 if the i th voxel belongs to the m th region and 0 if not. This allows us to extend Eq. (3) as follows:

Y = M m=1 Θ m v δm * A + E , (4)
The activation signals A capture for each voxel the periods of time during which any voxel is involved either in task performance or in spontaneous BOLD signal fluctuations. This model remains univariate as P independent neural activation signals ( a j ) P j=1 ∈ R P × T are estimated.

In our work, we rather introduce a low-rank constraint and learn both K temporal atoms (with K P ) and corresponding spatial maps. These maps encode various functional networks, each of them being summarized by specific neural activation profile. Mathematically, this can be modeled by replacing each vector a j in Eq. ( 4) with a linear combination of neural activation patterns Z = (z k ) K k=1 ∈ R K× T , with

z k ∈ R 1× T , modulated in space by the spatial maps U = (u k ) K k=1 ∈ R K×P , with
u k ∈ R 1×P , such that: A = K k=1 u k z k .
In other words, the spatial configuration u k encodes which voxels are linked to a given neural activation profile z k ∈ R 1× T . Note that a voxel may belong to different functional networks as overlapping can occur [START_REF] Karahanoglu | Total activation: fMRI deconvolution through spatio-temporal regularization[END_REF][START_REF] Najafi | Overlapping communities reveal rich structure in large-scale brain networks during rest and task conditions[END_REF]. Additionally, the sub-components of a functional network might also not activate simultaneously [START_REF] Allan | Functional connectivity in MRI is driven by spontaneous BOLD events[END_REF][START_REF] Betzel | Dynamic fluctuations coincide with periods of high and low modularity in resting-state functional brain networks[END_REF] as they can belong to different spatial maps u k . This is coherent with the fact that a voxel contains about one hundred thousands neurons for a typical spatial resolution (1.5 mm isotropic), and thus the underlying neural activation signals are combined with possibly different temporal fingerprints. Finally, our forward model for BOLD fMRI data is given by Eq. (5) (see also Fig. 3):

Y = M m=1 Θ m v δm * K k=1 u k z k + E .
(5)

Prior information and regularizing constraints

The number of unknown parameters Z ∈ R K× T , U ∈ R K×P and δ ∈ R 1×M is lower than the number of available data Y . Indeed, the observed fMRI data has P voxels of T scans while the proposed model has to infer K temporal components of T entries, K spatial components of P voxels each, and M HRF dilation parameters.

Hence, the number of unknown parameters to be set is

M + (K × P ) + (K × T ).
Importantly, if the number of scans T is large enough and if we adequately choose M and K, then we get P × T M + (K × P ) + (K × T ).

The forward model described in Eq. ( 5) is trilinear as it depends on both unknown spatial (U ) and temporal (Z) inputs and unknown convolution filters (v δm ).

This means that any timing variation in neural activation signals Z can be symmetrically compensated by the opposite variation in the haemodynamic filter v δm or that any sign or scale variation in U can be compensated by the inverse variation in Z or v δm . To overcome these sign, scale and time-shift ambiguities, we introduce some regularizing constraints in the computation of the parameter estimates. First, to avoid any sign ambiguity in the convolution between the neural input signals and the haemodynamic filter, the HRF filter has a constant and positive maximal amplitude (see Fig. 2). Second, as there is an interplay between the spatial and temporal components in the input signals, we also impose a nonnegativity constraint over the entries of each spatial map u k , i.e. ∀j, k kj ≥ 0.

By doing this, we only allow the deactivation of a functional network to be encoded in the corresponding z k . Third, to deal with the scale ambiguity, we impose [START_REF] Karahanoglu | Total activation: fMRI deconvolution through spatio-temporal regularization[END_REF]Cherkaoui et al., 2019;[START_REF] Urunuela | Stability-based sparse paradigm free mapping algorithm for deconvolution of functional MRI data[END_REF] we will assume the neural activation signals Z to be temporally piecewise constant or equivalently consider their first-order derivative ∇z k , ∀k to be sparse. Practically speaking, we add a total variation (TV) regularization term to our model.

u k 1 = η, ∀k = 1, . . . ,

Optimization problem

In this section, we derive an algorithm for estimating the dilation parameters δ = (δ m ) M m=1 , the spatial maps U and the neural activation signals Z from the model depicted in Eq. ( 5) and the aforementioned regularizing constraints. The estimates for these parameters can be obtained by solving the following constrained minimization problem: arg min

(U ,Z,δ) 1 2 Y - M m=1 Θ m v δm * K k=1 u k z k 2 F + λ K k=1 ∇z k 1 , subject to ∀k, u k 1 = η, ∀j, u kj ≥ 0, ∀m, δ m ∈ [0.5, 2.0] . (6)
The white Gaussian noise hypothesis leads to a quadratic data fidelity term, 2019), we enforce the first-order derivative of the temporal atoms (∇z k ) K k=1 to be sparse in order to constrain each z k to a piecewise constant signal. For that purpose, we use a TV regularization term, which corresponds to the 1 norm of the gradient in time ∇z k . Interestingly, the sparsity of the firstorder derivative of the temporal components (Dz k ) k tends to impose independence between these temporal components [START_REF] Daubechies | Independent component analysis for brain fMRI does not select for independence[END_REF]. This property is preserved by the discrete integration operator. Thus, the complete low-rank neural components should tend to be de-coupled due to the independence of the temporal processes (z k ) k . Importantly, this modeling of the neural activation signals allows us to fully adapt to task-fMRI and rs-fMRI experiments and to perform paradigm-free fMRI data analyses [START_REF] Deneux | EEG-fMRI fusion of paradigm-free activity using kalman filtering[END_REF]Caballero-Gaudes et al., 2011). In the first case, as the task-related BOLD signal is classically modeled as the convolution of an input block signal, representing the experimental paradigm (the onsets of the stimulus trials) with a HRF filter, we can recover neural activation signals close to the experimental paradigm. The neural activation profiles by being inferred from the data, allows estimation of input signals corresponding either to block or event-related paradigms, the second scenario being seen as a single time point block. More interestingly, the proposed framework is even more appealing for processing rs-fMRI data and uncover spontaneous and time-varying fluctuations of brain activity as the block duration may change from one instance to the next (see also [START_REF] Petridou | Periods of rest in fMRI contain individual spontaneous events which are related to slowly fluctuating spontaneous activity[END_REF][START_REF] Karahanoglu | Transient brain activity disentangles fmri resting-state dynamics in terms of spatially and temporally overlapping networks[END_REF]).

Moving to the technical aspects for solving the constrained optimization problem (6), it is worth mentioning that it is not globally convex. However, when U and δ are fixed, problem (6) becomes convex in Z and similarly when Z and δ are fixed, it becomes convex in U . Our minimization strategy of Eq. ( 6) thus relies on a block-coordinate descent algorithm, where we alternate the minimization between the two convex problems in U and Z followed by the non-convex one involving V . Also, the non-negativity constraints are activated when solving for the spatial maps U and the boundary constraints over δ are handled in parallel for each m, i.e. each HRF pattern when solving for V . Algorithm 1 details these three main steps.

Algorithm 1: Multivariate deconvolution and low-rank decomposition of the BOLD signal. -1) andδ (i-1) : (z

Input: BOLD signal Y , 1 initialization: ∀k, z (0) k = 0 T , u (0) k = u (init) k , δ (0) = δ (init) , i = 1 ; 2 repeat 3 Estimate the temporal atoms Z (i) with fixed U (i
(i) k ) k = arg min (z k ) k 1 2 Y - M m=1 Θ m v (i-1) δm * K k=1 u (i-1) k z k 2 F +λ K k=1 ∇z k 1 .
4

Estimate the spatial maps U (i) with fixed Z (i) and δ (i-1) :

(u

(i) k ) k = arg min (u k ) k 1 2 Y - M m=1 Θ m v (i-1) δm * K k=1 u k z (i) k 2 F , subject to {∀k, u k 1 = η and ∀j, u kj ≥ 0} .
Estimate the HRFs δ (i) with fixed U (i) and Z (i) :

(v (i) m ) m = arg min (δm)m 1 2 Y - M m=1 Θ m v δm * K k=1 u (i-1) k z (i) k 2 F , subject to δ m ∈ [0.5, 2.0] . 5 until J((z (i-1) k ) k ,(u (i-1) k ) k ,(v (i-1) m )m)-J((z (i) k ) k ,(u (i) k ) k ,(v (i) m )m) J((z (i-1) k ) k ,(u (i-1) k ) k ,(v (i-1) δm )m) ≤ ;
In regards to the (z k ) K k=1 step, we performed the minimization using an adaptiverestart accelerated forward-backward algorithm [START_REF] O'donoghue | Adaptive restart for accelerated gradient schemes[END_REF].

In regards to the (u k ) K k=1 step,we noticed that it can be cast as a standard quadratic program of the form arg min vec(U ) vec(U ) Avec(U ) -2vec(U ) b with linear and non-negativity constraints on vec(U ) , with vec(•) the usual matrixto-vector transformation. Although the matrix A and constraints are sparse, the problem is high-dimensional (P K variables, K linear constraints) and may not be efficiently solvable by standard mathematical programming solvers. Consequently, we benchmarked various algorithms in the dictionary learning literature and selected the most efficient, namely the one used to update the dictionary in [START_REF] Mairal | Online dictionary learning for sparse coding[END_REF]. Last, for the minimization with respect to (δ m ) M m=1 we used the accelerated forward-backward algorithm [START_REF] Beck | A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[END_REF][START_REF] Bioucas-Dias | A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF] after checking that it leads to a correct estimation of δ. The reader can find all details of the gradient computation w.r.t. (z k ) K k=1 and (u k ) K k=1 in Cherkaoui et al. (2019), for the gradient w.r.t. δ the computation is reported in the Supplementary Material, see Section A.1. For each step, we implemented the corresponding gradient (i.e. forward move) in an efficient manner to limit multiple computations over the iterations.

Critical steps for the efficiency of this algorithm are the computation of proxi-mal operators for the non-smooth regularizers. In regards to the neural activation patterns (Z step), as the minimization is sequentially performed over the K components, we only need to compute the proximal operator of the TV norm, i.e.

g z (z k ) = λ ∇z k 1 .
This remains a challenging issue as this operator is not closed form. A seminal contribution has been done in the literature for TV minimization [START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF]. Here, we rather use the Taut-String algorithm proposed by [START_REF] Barbero | Modular proximal optimization for multidimensional totalvariation regularization[END_REF] for which we use an efficient Python implementation

available in an open source package1 . In regards to the constraints on the spatial maps (U ), we also proceed separately on the K components: the proximity oper-

ator of g u (u k ) = 1 u k 1 =η + 1 u kj ≥0
where 1 stands for the indicator function2 , is

given by:

prox g u (u k ) = [(u kj -µ) + ] 1≤j≤P (7)
where µ is defined as P j=1 max{0, u kj -µ} = η and an efficient implementation has been proposed by [START_REF] Condat | Fast projection onto the simplex and the l1 ball[END_REF]. We propose to set λ as a fraction of λ max which is the minimal value of λ for which a constant Z is solution of Eq. ( 6),we detail the computation of λ max in Supplementary Material, cf. Appendix A.2. For the rest of the paper, we will refer to λ as the fraction of λ max , such as

λ = λ f λ max , with λ f ∈ [0, 1].
Algorithm 1 converges to a local minimizer Eq. ( 6) when each main iteration does not decrease sufficiently the cost function. In practice less than 50 iterations of the main loop are needed to converge. To initialize the spatial maps (u k ) K k=1 , we apply an Independent Component Analysis (ICA), implemented in scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF], on the BOLD signals Y and retain the produced spatial maps, we initialize the (z k ) k=1...K to zero and each entry of (δ m ) m∈{1..M } to 1. As the current optimisation problem is non-convex, we investigated different initialization strategies to assess how much the recovered parameters depend on the initial values (see Supplementary Material, Section A.4).

For the sake of reproducibility of the current results and others not shown here, we offer an efficient software implementation of our fast algorithm and with numerous code optimizations. The underlying hemolearn Python package is open source3 and available to the neuroimaging community.

Model Validation

We first validate the proposed approach on numerical simulations to illustrate the gain achieved by jointly estimating the neural activity profile and the HRF shape compared to a single deconvolution scheme. Next, we will demonstrate the usefulness of the proposed framework on real rs-fMRI data at the individual level.

In particular, we will highlight the impact of the hyper-parameter selection on the decomposition, describe the interpretation of its component and its stability. The usefulness of our method in the context of large cohorts will be investigated in Section 4. both activated regions (in bottom panels) while standard deviation across activated voxels is encoded by transparency around these mean curves. To make a comparison with respect to the state of the art, we deconvolved the observed signal using the method proposed in [START_REF] Wu | A blind deconvolution approach to recover effective connectivity brain networks from resting state fMRI data[END_REF]. This approach first detects the main neurally-related events by thresholding the observed BOLD signal. A binary sparse sequence associated with these spikes is thus constructed to be used as a proxy of a neural input signal. Then, for each voxel, the HRF profile is fitted to the data using this sequence of neural events. As a result, the BOLD signal can then be deconvolved voxelwise using the estimated HRF.

Numerical results

In a first step, we only estimated the pair (Z, U ) from the synthetic fMRI time series Y and kept the HRF profile v constant. The results are reported in Fig. 4(a).

The HRF shape used in this deconvolution process is shown in green in Fig. 4(a)

and actually differs from the true shape used for simulating the data. Because of this discrepancy in terms of haemodynamic delay and peak magnitude, the neural activation signals are not properly recovered (orange traces in Fig. 4(a)).

The magnitude of the estimates Z is much larger than the true one. This is partly due to compensate for the smaller magnitude of the HRF (green trace in large. However, we noticed that in both spatial maps, the non-negative magnitudes U are very well estimated. This is a direct consequence of using non-overlapping activating regions for the two neural traces. Additionally, we report, using a dashed green line, the case where the HRF profile is set to the true shape. In this case, we observe that the neural temporal components are accurately recovered.

In a second step, on the same data set Y we jointly estimated (Z, U ) and v using our full semi-blind deconvolution scheme. We kept the same initialization for the HRF shape for the sake of consistency. The results are reported in Fig. 4(b).

The HRF estimate v is shown in orange (dashed line) and actually matches the true curve. Consequently, the neural activation signals Z are properly estimated both in time and in magnitude and the corresponding RMSEs are one order of magnitude smaller than those reported in the previous simulated results (see Fig. 4(a)). This second synthetic setting did not impact the spatial maps, which are still well estimated.

Finally, for comparison purposes, we report in Fig. 4(c), the estimates computed by [START_REF] Wu | A blind deconvolution approach to recover effective connectivity brain networks from resting state fMRI data[END_REF]. We noticed an average error of 3s on the haemodynamic delay estimate. Moreover, for a selected voxel localized in the first activated region, we display the estimated neural events as well as the deconvolved BOLD signal and compare it to the true first temporal atom, its first-order derivative and the observed noisy BOLD signal. It is worth noticing that [START_REF] Wu | A blind deconvolution approach to recover effective connectivity brain networks from resting state fMRI data[END_REF] provided a poor alignment of the estimated neural spikes with the true ones as compared to our approach.

These results on synthetic data confirmed the good expected behavior of the proposed method. From a computational viewpoint, the estimation with constant HRF ran in 0.5s while the full estimation took 1s approximately on a machine with 15 GB of RAM and an Intel processor i7-7600U (2 physical cores, 2.80 GHz).

3.2. Single-subject analysis on rs-fMRI data

Dataset and parameter setting

To illustrate the proposed semi-blind deconvolution algorithm, we analyzed a single subject extracted from the UKBB resting-state fMRI data set. More investigation on a larger cohort of this data set will be presented in Section 4. The rs-fMRI data was 6 min10s long with TR = 0.735s. The first ten seconds were discarded (dummy scans) so that we end up with T = 490 scans (6min). The data was collected on a 3T Skyra Siemens MAGNETOM MR system at an isotropic resolution of 2.4×2.4×2.4 mm 3 using the multi-band GRE sequence (mb = 8). [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF].

In this subsection, we manually set the temporal regularization parameter to

λ f = 0.8 (remember λ f ∈ [0, 1]
). This setting achieves a bias-variance trade-off between two extreme situations, namely data overfitting (λ f = 0) on one hand and entirely sparse neural activation signals (λ f = 1 as λ = λ max ) on the other hand. The question of the unsupervised tuning of λ f is critical. It could be driven either from a statistical viewpoint (e.g. using the maximum likelihood criterion) that characterizes how likely the measured time series may be observed or using an external task and its corresponding metric such as classification performance (e.g. accuracy in prediction). Because the former approach does not admit a closed form solution, we explored in the Supplementary Material (see Section A) the impact of changing the temporal regularization (i.e. amount of sparsity in the activation neural signals) on the spatial map of haemodynamic delays. The spatial consistency we reported across regularization levels gave us confidence on the haemodynamic parameter estimates to be further used in subsequent classification tasks (see Section 4). In this setting, we implemented a cross-validation step with a leave-one-out loop to tune λ max in an unsupervised way.

In regard to the number of spatio-temporal atoms K, we set it using the explained variance (or R 2 -score averaged over the time series) as target metric in a preliminary study. For this set of parameters, the model estimation took around 1 minute on a machine with 15 GB of RAM and an Intel processor i7-7600U (2 physical cores, 2.80 GHz).

Results

Model selection. The first question we addressed on real rs-fMRI data was to optimally set the number of spatio-temporal atoms K and to find the best compromise between model complexity and model accuracy. For this purpose, we looked at two complementary criteria. The first one is standard and corresponds to the R 2 score that quantifies the variance explained by model ( 5) over the total sum of squares whereas the second one is given by the determinant of the correlation matrix between the neural activation signals.

The R2-score is defined as follows: R 2 = 1 -SSres SStot where SS tot quantifies the variance of the data Y and SS res the variance of the residuals after fitting model ( 5) by minimizing the cost function described in Eq. ( 6). The R 2 -score may vary from -∞ in pathological cases to 1 for a perfectly matching model. A good model is normally associated with R 2 > 0 and means that the L 2 norm of the residual is lower than the variance of the data. We therefore ran multiple model fitting for K in a range of {2, 3, 4, . . . , 10, 15, 20, . . . , 50}. The results are shown in Fig. 5(a)

and illcrustrate that the model accuracy first increases as a function of K up to reaching a plateau around R 2 0.55 for K = 20. So adding more spatio-temporal components no longer improves its ability to capture variability in the data while it becomes more complex.

The second information measure we used to help us select K was based on the determinant of the correlation matrix

Σ K = (E[(z k -m k )(z -m ) T /σ 2 k σ 2 ]) k,
between the temporal atoms (z k ) K k=1 . The quantities σ 2 k and σ 2 define the variance of the neural activation signals z k and z . As Σ K is semi-positive definite with entries between 0 and 1, its eigenvalues are positive or null and so its determinant varies between 0 and 1: det Σ K = 1 when matrix Σ K defines a basis, which means that all atoms are orthogonal and decorrelated like in a PCA decomposition. In contrast, det Σ K = 0 when matrix Σ K is not of full rank so at least one atom could be obtained as a linear combination from the others. Therefore, as before we ran multiple model fitting for K in a range of {2, 3, 4, . . . , 10, 15, 20, . . . , 50}

and we plotted in Fig. 5(b) the evolution of the determinant of Σ K as a function of K. The results show us that beyond K ≥ 20, we get a correlation matrix with det Σ K ≤ 10 -10 which tends to zero. According to this criterion, from the we should therefore not exceed 20 temporal atoms. Thus, from the R2-score criterion and this det Σ K criterion, in our following experiments we will keep K = 20.

Analysis of spatial decomposition. Fig. 6 shows the spatial maps of this spatio-temporal decomposition for this individual and Tab. 1 summarizes the list of functional networks retrieved in this setting. It is worth mentioning that the lates with components 10 and 20 pretty strongly6 between 0.46 and 0.54. Overall, this analysis shows that the proposed approach does not separate RSN in single components. However, it still achieves a meaningful decomposition. Although we do not report here the spatial decomposition for K ∈ {5, 8, 10, 15, 30, 40, 50}, the latter can be easily obtained using the HemoLearn toolbox7 . Interestingly, in such decomposition some RSNs are either not recovered or mixed together for small K ≤ 15 whereas for large values of K ∈ {25, 30, . . . , 50} each RSN is split in multiple components. Increasing the number of components extends the range of variation of the correlation coefficients in both positive and negative senses. This confirms that our approach does not have statistical independence guaranties like PCA and ICA do when decomposing the neural activity. In this subsection, we manually set the temporal regularization parameter to λ f = 0.8 as in the prior rs-fMRI analysis and the number of spatio-temporal atoms to K = 30 to offer a richer spatio-temporal decomposition with possibly more focal components, compared to resting-state analysis. We compared our approach to [START_REF] Pedregosa | Data-driven HRF estimation for encoding and decoding models[END_REF] in which a voxelwise HRF estimate is yielded using a rank-1 GLM model. This permits a fair comparison of both haemodynamic delays and spatial activation maps corresponding to a given condition from the EP.

Results

For illustrative purposes, we focused on the left-hand motor task. Fig. 9(a)

represents the 100 most activated voxels estimated by [START_REF] Pedregosa | Data-driven HRF estimation for encoding and decoding models[END_REF] in response to this condition. These voxels were obtained by thresholding the regression coefficients of the model displayed in 2015)'s method due to the recovered activity in the visual cortex vs longer TTP for our method that picks up the activation in the primary right motor cortex).

haemodynamic estimation stability over time

The shape of the HRF is controlled by the neurovascular coupling including both neural and non neural factors such as glial cell activity, cerebral energy metabolism, and the cerebral vasculature. Abnormalities in the local vascular system or cell communication due to pathological state or changes in cerebral blood flow upon psychoactive drugs could influence this haemodynamic response. As we expect the HRF estimate to be stable if none of those events took place, we propose to study the within-subject stability of HRF estimates over time, namely between consecutive time periods. For doing so, we compare the within-subject variability of the HRF whole brain dilation parameter vector δ to the inter-subject variability of the same quantity. We thus introduce two reference 2 distances, namely the within-subject distance WS(δ s 1 ,

δ s 2 ) = δ s 1 -δ s 2 2
2 where (δ s 1 , δ s 2 ) correspond to the vectors of spatially aggregated HRF dilation parameters that were estimated over two periods of time T 1 and T 2 in the same individual s. Similarly, for any pair of subjects (s 1 , s 2 ) and a given period T , we measure the between-subject distance between (δ s 1 T , δ s 2 T ) vectors as follows: BS(δ

s 1 T , δ s 2 T ) = δ s 1 T -δ s 2 T 2 2 .
The goal is then to compare the within-and between-subject distances across individuals and show that the intra-subject variability is significantly lower than the inter-subject one over a sufficiently large population.

Data set and numerical analysis

We selected 100 healthy subjects from the Human Connectome Project (HCP) data set [START_REF] Van Essen | The wu-minn Human Connectome Project: An overview[END_REF] at random. We used this data set because of the availability of a 12-min long rs-fMRI run with a short time of repetition (T R=0.753s), see [START_REF] Glasser | The minimal preprocessing pipelines for the human connectome project[END_REF] for a full description of the acquisition parameters and the pre-statistics processing steps. In this rs-fMRI run for each individual, we extracted two segments of 4 minutes each, denoted as T 1 and T 2 hereafter, the first and last parts of the recording. We then applied the proposed multivariate spatio-temporal decomposition to each segment using K = 8 spatio-temporal atoms (z k , u k ) K k=1 and a brain atlas Θ = (Θ m ) M m=1 [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF] composed of M = 96 regions of interest (ROIs). This haemodynamic brain parcellation thus yields 96 HRF dilation parameters δ = (δ m ) M m=1 for each individual. In practice, in the definition of WS(•, •) and BS(•, •), the true vectors δ s i T j (i = 1, . . . , 100, j = 1, 2) have been replaced by their estimates δ s i T j computed by solving Eq. ( 6) for the two 4-min rs fMRI data sets (T 1 and T 2 ). To make sure that our conclusions hold for a large scale of temporal regularization parameters, we spanned the range λ f ∈ [0, 1] and repeated the same procedure over 6 discrete values of λ f within this interval.

Results

In Fig. 10, the box plot in blue shows the within-subject distance WS(

δ s i T 1 , δ s i T 2 )
between the two 4-min rs-fMRI segments for all individuals and across 6 values of λ f covering the whole interval [0, 1]. The orange and green box plots in Fig. 10 depict the between-subject distances computed over the first and second segments respectively, namely BS(

δ s i T 1 , δ s j T 1 ) and BS( δ s i T 2 , δ s j
T 2 ) with i = j. We observed that the within-subject (i.e. inter-segment) variability is systematically lower than the between-subject variability and that all metrics remain stable across regularization levels. To go further, we performed a statistical analysis (paired t-test) by comparing the mean of the WS and BS distributions and we obtained significant pvalues (p < 10 -8 ) showing that the within-subject haemodynamic variability is significantly lower than the between-subject fluctuations. In contrast, the statistical inter-individual comparison between the two segments is not significant (p 10 -2 ).

These results are valid for all tested regularization levels indicating a minor impact of the regularization parameter onto the haemodynamic parameter estimate.

In sum, this analysis demonstrates that the whole brain characterization of the vascular system remains stable in a given individual between two periods shortly spaced in time, compared to the same analysis between individuals and so that the haemodynamic response discriminates each subject from the others.

Clinical validation at the population level

In the previous section, the numerical experiments were devoted to demonstrate the meaningfulness and reliability of the proposed multivariate spatio-temporal within-subject decomposition of fMRI data, especially in resting-state experiments.

In this section, our main objective is to showcase the application of this approach to clinical diagnosis. For this purpose, we leverage the functional features (haemodynamic delays, neural activation signals, etc.) output by our approach to first characterize patients with history of stroke compared to healthy controls and then to discriminate middle-age vs elderly subjects. In both analyses, we again used the 6-min long rs-fMRI data from the UK Biobank database.

Characterization of patients with an history of stroke

Stroke is a medical condition in which the blood supply to is interrupted or reduced in a brain area, resulting in ischemic brain tissue and neuronal damage.

This pathology is considered as a major health issue nowadays [START_REF] England | Briefing document: First incidence of stroke estimates for england 2007 to 2016[END_REF].

In this field, multiple studies [START_REF] Min | Development of an algorithm for stroke prediction: A national health insurance database study in korea[END_REF] have proposed approaches to better estimate the stroke risk for patients. However, in these attempts a major issue is the precise estimation of the brain damage that occurs in the neurovascular system during and after a stroke episode. To that purpose, we tested our approach to characterize the effect of stroke on the haemodynamic response in the brain.

We considered 24 patients of both genders and various ages who suffered from a stroke in the past from the UK Biobank database. For comparison purposes, we selected 24 healthy controls matched in age and gender from the same database.

We applied the same decomposition (K = 20, M = 96, same λ f ) to each patient and healthy control. all dilation parameter values by their within-subject average, namely δHC and δSP respectively. We first observed that the dilation parameters were larger in average in the healthy condition compared to stroke ( δHC > δSP ). This corresponds to shorter and more homogeneous TTPs in the brain in the healthy condition. The shortest TTP found in the healthy control was actually located in the primary visual cortex (axial slice, z=-2, left hemisphere), a result consistent with the literature on fastest haemodynamic responses often detected in visual areas [START_REF] Handwerker | Variation of bold hemodynamic responses across subjects and brain regions and their effects on statistical analyses[END_REF][START_REF] Badillo | Group-level impacts of within-and betweensubject hemodynamic variability in fMRI[END_REF]. In contrast, Fig. 11(b) illustrates that the haemodynamic dilation parameters δ SP are smaller (so the TTPs longer) in the stroke patient Also, we found less variability in the healthy subject since the difference between the maximum and minimum TTPs were smaller (∆ HC TTP = 1.25 s) compared to the stroke patient (∆ SP TTP = 2.25 s). Importantly, Fig. 11(a) illustrates the relative symmetry of haemodynamic territories that exists in normal subjects between both hemispheres [START_REF] Raemaekers | Knowing left from right: asymmetric functional connectivity during resting state[END_REF].

On the contrary, Fig. 11(b) reveals a wider asymmetry between the two hemispheres in the stroke patient. Interestingly, in this patient we noticed the presence of larger TTPs in the middle left precentral gyrus and left motor cortex (resp.

z=44 and z=60), namely the brain regions supposedly impacted by the stroke episode.

To go one step further, we quantified the spatial asymmetry of the haemodynamic structure within each individual. To this end, we computed the interhemispheric haemodynamic 2 distance (IHD) between the HRF dilation parameter vectors estimated over the left and right hemispheres in laterally matched brain regions, respectively denoted δ L and δ R . This normalized intra-subject distance is defined as follows:

IHD(δ s R , δ s L ) = δ s L -δ s R 2 δ s L+R 2
, ∀s = 1, . . . , 24.

A zero-valued distance thus reflects a perfect symmetry of the estimated haemodynamic responses. In contrast, we expect to uncover asymmetry between haemodynamic territories respectively localized in the ischemic and normal hemispheres.

By pulling down the values of IHD across all individuals within each group (HC vs SP), we estimated the IHD distributions for the two populations of interest, as shown in Fig. 11(c). In the latter graph, we illustrate how different the two cohorts are in terms of neurovascular asymmetry. The group of 24 stroke patients exhibit larger haemodynamic differences between the ischemic and normal hemispheres.

We statistically assessed such difference between the two distributions using a twosample Kolmogorov-Smirnov test and found a significant p-value (p = 9.1 10 -5 ).

This quantification thus confirmed preliminary visual assessment. We report this p-value with a temporal regularization defined such as λ f = 0.001. However, we obtain similar p-value results when using the 5 others levels of temporal regularization (λ f ∈ [0.001, 0.9]).

In summary, this analysis has shown that the proposed framework is instrumental in discriminating healthy subjects from stroke patients, both at the individual and group-levels, using haemodynamic features and an neurovascular asymmetry index, which allowed us to localize pathological haemodynamic delays.

Middle-age vs elderly subjects classification

In the previous part, we performed group-level statistical analysis in the classical way. In this part, we intend to assess the prediction power of the proposed framework in order to classify middle-age vs elderly subjects using standard machine learning tools [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. The reason for choosing this classification task between middle-age and elderly individuals lies first in the fact that multiple studies have pointed out the modification of the haemodynamic system with healthy aging [START_REF] Ances | Effects of aging on cerebral blood flow, oxygen metabolism, and blood oxygenation level dependent responses to visual information stimulation[END_REF][START_REF] Li | Agingassociated changes in cerebral vasculature and blood flow as determined by quantitative optical coherence tomography angiography[END_REF][START_REF] West | BOLD hemodynamic response function changes significantly with healthy aging[END_REF] and second in a regain of interest in the literature for brain age analysis using multiple neuroimaging techniques [START_REF] Engemann | Combining magnetoencephalography with magnetic resonance imaging enhances learning of surrogate-biomarkers[END_REF]. We thus intend to assess whether our approach is able to capture the effect of aging and if so, whether the haemodynamic features reflect more brain aging compared to the neural activity signals.

For that purpose, we still relied on the UK Biobank database as in the previous experiments as the short TR (TR=0.735 s) in the rs-fMRI data set (6 min long) provides a suitable setting to investigate the HRF evolution with aging. Here, we selected 459 healthy subjects of both genders and divided them in two balanced groups: the middle-age (MA: 40-44 yo) and elderly (E: 64-70 yo) groups. We applied the decomposition (6) to each subject using 5 levels of temporal regularization (λ f ∈ [0.001, 0.9]), K = 20 temporal components and we used the same brain parcellation of M = 96 ROIs as before to segregate the HRFs in space.

First, we analyzed the haemodynamic differences between the two populations by computing a two-sample t-test on the distributions of dilation parameters. We used the temporal regularization parameter λ f = 0.675 which is the one selected through cross validation for our classification model in the subsequent paragraph.

The results were first quantified with t-scores to compare the dilation parameters in each region. The results are presented in Fig. 12(a) and illustrate large differences between the two populations. More specifically, higher dilation parameters or shorter TTP were retrieved in middle-age subjects as the t-scores were mostly positive for the comparison δMA m > δE m . This is notably visible in the Willis Polygon, temporal cortices, angular gyri, the medial prefrontal cortices and the superior frontal cortices. To assess the statistical significance, we also computed the log-transformed p-values, i.e. -log 10 p val (shown in Fig. 12(b)), after correcting for multiple comparisons using the Bonferroni correction across the M = 96 ROIs.

We noticed first that a large majority of significant brain regions appear bilaterally indicating larger haemodynamic dilation parameters or shorter TTPs in younger individuals. Second, the negative t-values reported in the cerebellum are not statistically significant after correcting for multiple comparisons.

Then, in an attempt to be exhaustive we constructed three different logistic regression (LR) models based either on the individual (i) neural activity signals,

(ii) HRF shapes and (iii) haemodynamic dilation parameters. We did not consider the spatial activity maps as input features in these models as they do not permit to perform dimension reduction. We trained these LR models using the Scikit-Learn software [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] to predict the age label for each subject (1 for the elderly people, 0 for middle-aged people). A 2 -norm regularization was used in the estimation of the LR model parameters with an hyper-parameter β > 0. We grid-searched the temporal regularization parameter λ f and the classifier hyperparameter β (see Fig. E in Supplementary Materials for the stability of the setting).

We chose the accuracy as the classification metric and cross-validated the score to provide an estimation of the generalization error with a 10 times repeated stratified 4-fold split.

Fig. 13(a) shows that the haemodynamic properties have an improved prediction power to discriminate the age compared to the neural activation signals (i.e.

temporal components). The latter actually reaches an average accuracy score of 0.557, whereas the mean accuracy associated with the HRF shape and haemodynamic dilation parameter estimates respectively goes up to 0.741 and 0.743.

Also, the distribution of accuracy scores across trials is more concentrated for the HRF dilation parameters compared to the whole HRF shape. This is likely due to the dimension reduction operated to extract this parameter which fluctuates less than the complete profile of the haemodynamic response. This analysis thus demonstrates that our decomposition is able to capture the brain age based on neurovascular information. However, due to the large between-subject variability even within each class of age the neural activation signals do not define a good feature for the brain age prediction. Complementary to that, Fig. 13 Overall, this experiment has permitted to demonstrate that haemodynamic features are a good biomarker of the normal aging, as already reported in the literature [START_REF] Grady | Understanding variability in the bold signal and why it matters for aging[END_REF][START_REF] West | BOLD hemodynamic response function changes significantly with healthy aging[END_REF]. Moreover, it highlighted that the inter-hemispheric asymmetry in neurovascular coupling brings key information to discriminate middle-age from elderly people.

Discussion

Separating vascular and neuronal components on fMRI BOLD signals. Both neural and non neural factors such as glial cell activity, cerebral energy metabolism and the cerebral vasculature contribute independently and synergistically to the fMRI BOLD signal. A mis-estimation of individual and regional HRFs may lead to an interpretation that haemodynamic changes as neural variations could have considerable implications for the interpretability and reliability of findings in fMRI studies. Previous literature supports the notion that HRF variability corrupts fMRI data deconvolution [START_REF] Gitelman | Modeling regional and psychophysiologic interactions in fMRI: the importance of hemodynamic deconvolution[END_REF][START_REF] Roebroeck | The identification of interacting networks in the brain using fMRI: model selection, causality and deconvolution[END_REF][START_REF] Seth | Granger causality analysis of fMRI BOLD signals is invariant to hemodynamic convolution but not downsampling[END_REF]. Deconvolution-based approaches such as the one presented here allows discrimination of cerebrovascular components from neural activaty and minimizes the confound of HRF variability in the exploration of brain physiology, functional connectivity and cognitive processes.

In this work, we developed a new algorithm that proposed the joint estimation of the HRF and neural activation signal as a semi-blind deconvolution multivariate problem in a paradigm-free setting. Synthetic and real resting-state fMRI data allowed us to demonstrate that this approach is able to faithfully capture the individual's haemodynamic response function and intrinsic functional networks with low intra-subject variability and relative minimal impact of hyper-parameters on the reliability of HRF estimation. Aims of this paper were also to verify that these methodological developments have practical impacts as related to detection and classification. Further demonstrations using well defined research protocols will be required to refine the tool for use in clinical applications.

Related works. Disentangling the neurovascular coupling from the neural activity is well documented in the literature [START_REF] Goutte | Modeling the haemodynamic response in fMRI using smooth FIR filters[END_REF][START_REF] Ciuciu | Unsupervised robust nonparametric estimation of the hemodynamic response function for any fMRI experiment[END_REF][START_REF] Handwerker | Variation of bold hemodynamic responses across subjects and brain regions and their effects on statistical analyses[END_REF][START_REF] Lindquist | Validity and power in hemodynamic response modeling: a comparison study and a new approach[END_REF][START_REF] Vincent | Spatially adaptive mixture modeling for analysis of fMRI time series[END_REF][START_REF] Chaari | Hemodynamic-informed parcellation of fMRI data in a joint detection estimation framework[END_REF][START_REF] Pedregosa | Data-driven HRF estimation for encoding and decoding models[END_REF], however many features present in our model are pivotal to making a significant step forward. First, unlike the previous works our approach is paradigm-free and can be applied to both task-related and resting-state fMRI data. As such, the proposed methodology is not restricted to a few brain regions (e.g. visual, auditory and sensorimotor cortices) that are typically involved in a specific activation. Second, in contrast to [START_REF] Wu | A blind deconvolution approach to recover effective connectivity brain networks from resting state fMRI data[END_REF][START_REF] Farouj | Bold signal deconvolution under uncertain haemodynamics: A semi-blind approach[END_REF]Cherkaoui et al., 2019) we propose a multivariate approach to reliably summarize the neurovascular coupling and the neural activity over the whole brain. Third, to the best of knowledge, the proposed semi-blind deconvolution approach is the first to perform such an extensive experimental validation including i) simulation-based numerical experiments for model validation, ii) performance assessment on taskrelated and resting-state fMRI data at the individual level, including comparisons to [START_REF] Wu | A blind deconvolution approach to recover effective connectivity brain networks from resting state fMRI data[END_REF][START_REF] Pedregosa | Data-driven HRF estimation for encoding and decoding models[END_REF], and iii) two meaningful cohort-level evaluations in different clinical settings. It is worth mentioning that we used a common brain parcellation -the Harvard-Oxford probabilistic atlas [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF] -across all individuals, to ensure group-level analysis and facilitate between-group comparisons. Hence, our main findings on the asymmetries in haemodynamic features between stroke patients and healthy controls on one hand, and the differences related to normal aging on the other hand, are dependent on this atlas and could slightly differ with another parcellation. This question is left for future research.

Interest in identifying the vascular component for clinical applications. Non-invasive imaging modalities are of undeniable relevance for the diagnostic and prognostic work-up in patients and this present study constitutes a proof of concept in terms of the interest and feasibility of the proposed approach. This type of analysis was made possible due to the fast convergence of our algorithm and the numerous optimizations on the implementation side. By making the HemoLearn Python package open source, we offer a unique tool to the neuroimaging community that will permit complementary investigation in other clinical contexts.

The current challenge for diagnostic imaging methods is to find metrics that capture relevant information including biomarkers and the present work on haemodynamic blind deconvolution might help uncover these measurable indicators. In that context, in agreement with various imaging studies on aging [START_REF] Bangen | Differential age effects on cerebral blood flow and BOLD response to encoding: associations with cognition and stroke risk[END_REF][START_REF] Shafto | The cambridge centre for ageing and neuroscience (Cam-CAN) study protocol: a cross-sectional, lifespan, multidisciplinary examination of healthy cognitive ageing[END_REF], our algorithm, tested on a large cohort of UK BioBank rs-fMRI data sets (459 subjects) has proven its sensitivity to classify middle-age vs elderly individuals with respect to the estimated neurovascular coupling. Vascular aging, as characterized by progressive deterioration in the cellular structure of the blood vessel wall [START_REF] Handwerker | Variation of bold hemodynamic responses across subjects and brain regions and their effects on statistical analyses[END_REF][START_REF] Havlicek | A dynamical model of the laminar BOLD response[END_REF], undoubtedly impacts both resting-state cerebral blood flow and the ability to modulate it during neuronal activity. Our study confirmed a regional variability in the HRF features with a significantly slower neurovascular coupling in elderly people. Importantly, we demonstrated that the dilation parameter δ, tightly coupled to the haemodynamic delay, achieves better accuracy scores (0.74) compared to the neural activation signals (0.56, just above chance level), suggesting a causal role of cerebral micro-vascular dysregulation due to cognitive aging.

Based on the proposed haemodynamic asymmetry index, our analysis tested in patients with a history of stroke has also proven its utility to individually detect very slow haemodynamic delays in a restricted brain territory probably related to local ischemic tissue consecutive to stroke. This finding is perfectly consistent with the literature [START_REF] Altamura | The longitudinal changes of BOLD response and cerebral hemodynamics from acute to subacute stroke. A fMRI and TCD study[END_REF] showing that the delay in peak latency that arises as patients advance from the acute to the subacute stroke phase is related to the deterioration of cerebral haemodynamics. Diffusion-weighted MRI and perfusion imaging (e.g. ASL) remain the reference imaging modalities as a noninvasive method to perform the diagnosis in the acute episode of stroke.

However, in the post-acute period, rs-fMRI acquisition which does not require the patient's engagement would be extremely valuable to measure the re-establishment of haemodynamic response function in stroke patients as a prediction of subsequent recovery of cerebral function.

Assessment of cerebrovascular function through approaches that employ haemodynamic deconvolution --as presented in this paper --bring new opportunities for the exploration of brain plasticity and pathogenesis. For example, modeling the haemodynamic response function improves sensitivity of fMRI data to delineate epileptogenic area [START_REF] Storti | Modelling hemodynamic response function in epilepsy[END_REF]. Introducing some flexibility in the HRF shape is of undeniable benefit to consider 'silent' vascular risk factors in the pathogenesis or exacerbation of neurological diseases (D 'Esposito et al., 2003;[START_REF] Ogaki | Vascular abnormalities and the role of vascular endothelial growth factor in the epileptic brain[END_REF] or in pharmacological fMRI studies (phMRI) where drugs, used as stimuli, do not always elicit predictable changes in neurovascular coupling (Cherkaoui et al., 2021).

Limitations and extensions. Some limitations of our tool do exist. First, there are free parameters in the proposed modeling (K, λ f , M ) that need to be set in an appropriate manner. We explored two model selection criteria for setting K, namely the R 2 score and the determinant of the correlation matrix between the neural activation signals. Based on these metrics, we found a fair compromise between accuracy and model complexity for K = 20. We thus constantly used this value hereafter in the individual decomposition. Of course, other model selection approaches might be envisaged to optimize K and λ f using for instance a (widely)

Bayesian information criterion [START_REF] Neath | The Bayesian information criterion: background, derivation, and applications[END_REF][START_REF] Watanabe | A widely applicable Bayesian information criterion[END_REF], or the log-likelihood in the standard classical framework. The selected model would thus be the one associated with the lowest BIC value or largest log-likelihood. More recently, the concept of bi-level optimization [START_REF] Bennett | Model selection via bilevel optimization[END_REF] erator [START_REF] Cherkaoui | Learning to solve tv regularised problems with unrolled algorithms[END_REF] for instance could also be directly plugged into the current algorithm. Fourth, as in standard multivariate data-driven methods, the inter-subject comparison of spatial maps is currently difficult in the proposed formulation. In the same spirit as group-ICA [START_REF] Calhoun | A review of group ICA for fMRI data and ICA for joint inference of imaging, genetic, and erp data[END_REF], canonical ICA [START_REF] Varoquaux | A group model for stable multi-subject ICA on fMRI datasets[END_REF] or multi-subject dictionary learning [START_REF] Varoquaux | Multi-subject dictionary learning to segment an atlas of brain spontaneous activity[END_REF], the current within-subject decomposition could be extended to the group-level to become more stable. One possibility would be to impose the same spatial maps across all individuals like in [START_REF] Calhoun | A review of group ICA for fMRI data and ICA for joint inference of imaging, genetic, and erp data[END_REF] while another more flexible approach would permit spatial variations around a group-level spatial template [START_REF] Varoquaux | Multi-subject dictionary learning to segment an atlas of brain spontaneous activity[END_REF]. In this context, the neural activation signals could remain subject-specific with large fluctuations both in timings and magnitudes.

This kind of extension will be investigated in the near future.

Fifth, we experimentally observed both on numerical simulations and real fMRI data (ADHD cohort [START_REF] Milham | The adhd-200 consortium: a model to advance the translational potential of neuroimaging in clinical neuroscience[END_REF]) that a TR larger than 1 s may be detrimental to a precise estimation of the haemodynamic dilation parameter. For that reason, all analyses were performed on fMRI acquisitions with short TR. This type of data is usually collected using simultaneous multi-slice imaging [START_REF] Feinberg | Ultra-fast MRI of the human brain with simultaneous multi-slice imaging[END_REF][START_REF] Hesamoddin | Advantages of short repetition time resting-state functional MRI enabled by simultaneous multi-slice imaging[END_REF] to keep this parameter below 1s.

Sixth, because the proposed HRF model relies solely on a time dilation parameter, its magnitude is fixed and the fluctuations of the BOLD signal across the brain are thus captured through the neural activity atoms (z k ) K k=1 on one hand and the spatial maps (u k ) K k=1 ) on the other hand. However, the norm of the spatial maps being constrained, the real BOLD signal amplitude is captured by the neural activation signals. A recent work [START_REF] Tsvetanov | The effects of age on resting-state bold signal variability is explained by cardiovascular and neurovascular factors[END_REF] has shown that the resting-state fluctuation amplitude is crucial to predict brain age in healthy subjects. One possible enhancement of the current model would be to add a magnitude parameter to each HRF. In that case, we should fix the scale ambiguity issue by setting the amplitude of the temporal atoms (z k ) K k=1 .

This modification would significantly increase the computational complexity due to the calculation of the proximal operator associated with the new regulariza-

tion term g z ((z k ) k ) = λ K k=1 ∇z k 1 + I z k ∞=α .
Last, thus far we have used the canonical HRF as the reference shape in v ref . This setting could be easily updated to perform investigations in specific populations (e.g. newborns) where the true haemodynamic response function is known to deviate from the canonical shape [START_REF] Arichi | Development of BOLD signal hemodynamic responses in the human brain[END_REF].

Conclusion

In this paper, we have presented a semi-blind deconvolution approach to jointly estimate the haemodynamic response function and the neural activity signals across the whole brain. As the proposed methodology is paradigm-free, it enables the analysis of resting-state fMRI data in an semi-supervised manner as the regularization parameters (K, λ f ) may be tuned using a trade-off between model accuracy and complexity. Beyond the model validation on synthetic and real fMRI data, we have demonstrated the interest of the proposed approach in two applications in neuroscience. Both aimed at characterizing cerebral haemodynamic delays in specific populations, namely stroke patients and elderly people by contrasting them with healthy and younger controls. Most importantly, we proposed an haemodynamic asymmetry index to lateralize the stroke episode while confirming the presence of a prolonged haemodynamic delay in these patients. We also demonstrated that haemodynamic properties are predictable of brain age. Finally, this new framework opens the door to new research avenues for functional connectivity analysis based on the neural input signals instead of the BOLD signal themselves. In contrast to existing techniques [START_REF] Wu | A blind deconvolution approach to recover effective connectivity brain networks from resting state fMRI data[END_REF], our approach would be less biased by a constant haemodynamic response shape across the whole brain. We depict the axial slices of the corresponding activation maps for each approach ((c) refers to our competitor, (d) to our method). We observe that [START_REF] Pedregosa | Data-driven HRF estimation for encoding and decoding models[END_REF] mainly identify the activation in the cerebellum whereas our approach localizes the activation in the right motor cortex. Bottom (e): We display the neural activation time courses estimated by our approach in orange for the selected voxels. The standard deviation across voxels is encoded by transparency around the mean curve. Similarly, we report the observed BOLD signal in blue and display the temporal profile of the motor task in light gray. Bottom (f ): We display the HRF estimates for each method. The green HRF corresponds to [START_REF] Pedregosa | Data-driven HRF estimation for encoding and decoding models[END_REF] with the between-voxel standard deviation encoded by transparency around the mean curve. Similarly, we depict our HRF estimate in orange. We notice larger haemodynamic delays estimated by our method. ) with i = j. These distributions are assessed for 6 levels of temporal regularization (λ f ∈ {10 -6 , 10 -3 , 10 -2 , 10 -1 , 5 10 -1 , 1 -10 -6 }) and remain stable. Statistical analysis (paired t-test) was conducted to assess the significance of the difference between the mean of the within-and between-subject 2 2 norm distributions. Significant differences are marked with a * . k=1 , (ii) HRF shape estimates ( v δm ) 96 m=1 or (iii) the haemodynamic dilation parameters ( δ) 96 m=1 . From top to bottom, the distribution of the classification scores is shown from for the predictive features (i)-(iii), respectively. The best accuracy scores (average 0.74) are reached using the haemodynamic parameters and the smallest variability in the prediction using specifically the dilation parameter estimates. (b): Learning curve of accuracy scores as a function of the number of individuals (middle-age vs elderly subjects) used for the training stage both for the haemodynamic dilation parameter (blue curve) and the temporal components (orange curve). As a plateau is reached for 459 people, we presented the corresponding performances in panel (a).

Figure 1 :

 1 Figure 1: Illustration of the voxelwise BOLD signal modeling: a is a time series encoding the neural activation signal, v being the haemodynamic response function (HRF) and e the additive white Gaussian noise. The measured fMRI signal is denoted y and obtained via a * v +e.

2. 1 .

 1 Linear and time-invariant modeling 2.1.1. Univariate modeling A common model for the multivariate (P voxels, T scans) BOLD data Y ∈ R P ×T with Y = (y j ) P j=1 is the linear and time-invariant model (LTI) (Boynton

Figure 3 :

 3 Figure 3: Illustration of the low-rank multivariate BOLD signal model (the colors are here for illustrative purposes). Y stands for the observed BOLD data, v m the m th HRF, Θ m represents the predefined brain parcellation, U = (u k ) K k=1 defines the K spatial components and Z = (z k ) K k=1 the corresponding "neural" activation signals. Last, E represents the additive white Gaussian noise.

  to measure how well we reconstruct the observed fMRI signals. Moreover, in alignment with with Caballero-Gaudes et al. (2012); Karahanoglu et al. (2013); Cherkaoui et al. (

  We generated two temporal Dirac signals of length T = 500 with a fixed sparsity level. The first generated Dirac signal is composed of randomly drawn signed spikes, with location chosen uniformly in time and intensity drawn from a Gaussian distribution N (0, 1) to simulate spontaneous (i.e. resting-state) neural activity fluctuations. To produce the second temporal atom, we generated 4 spikes (2 with positive amplitudes, 2 with negative amplitudes) to simulate a taskrelated response. To produce the corresponding block signals Z -shown in blue in Fig. 4[right panels] -we integrated over time these signals and convolved them with a predefined HRF v δ to yield two corresponding pure BOLD time series. The chosen HRF has length L = 25 and is shown in blue in Fig. 4[left panel]. For the sake of simplicity, we considered a single HRF profile (M = 1) in this synthetic setting, so the haemodynamic properties were supposed constant in space. We then assigned these BOLD signals to spatial locations. Hence, we defined K = 2 corresponding 2D maps U = (u k ) 2 k=1 (20 × 20, i.e. P = 200). Each spatial map has a single activating squared region consisting of 4 × 4 = 16 pixels. Each active pixel has a randomly drawn non-negative magnitude, the other ones being set to zero. Then, we normalized each map by its 1 -norm. Finally we added white Gaussian random noise to produce observed, i.e. noisy BOLD signals Y of length T = T +L-1 = 524 scans (T R = 1s) with a signal-to-noise-ratio (SNR) of -10 dB. The mean synthetic BOLD signals are reported in black traces in Fig. 4(a)-(b) for

Fig. 4

 4 Fig. 4(a)) used for deconvolving the BOLD signals. Consequently, the residual mean square errors (RMSEs) computed on the neural activation signals are pretty

  4 Standard preprocessing steps were applied: motion correction usingMCFLIRT Jenkinson et al. (2002); grand-mean intensity normalisation of the entire 4D data set by a single multiplicative factor; high-pass temporal filtering (Gaussian-weighted least-squares straight line fitting); EPI unwarping; gradient distortion correction unwarping. Structured artefacts are removed by ICA processing, see the documentation 5 for a full description. Finally, we down-sampled the voxels to reduce the spatial dimension of our decomposition to 61, 605 voxels to match the resolution of 2 mm isotropic of the Harvard-Oxford probabilistic atlas

  Fig.8depicts the correlation matrix between the corresponding neural activation signals. It is then insightful to notice that the correlation between the multiple components in a given network are quite strong. For the visual network we observed a correlation coefficient varying between 0.31 and 0.48, the largest value being reached for areas located both in the extrastriate cortex. The same conclusion holds in the motor network with a correlation level varying between 0.1 and 0.36. In regard to the DMN, component 11 plays the role of a hub as it corre-

3. 3 .

 3 Single-subject analysis on task-fMRI data3.3.1. Data set and parameter settingTo further illustrate the proposed semi-blind deconvolution algorithm, we analyzed a single subject extracted from the HCP[START_REF] Van Essen | The wu-minn Human Connectome Project: An overview[END_REF] data set and focused on the motor task in the activation fMRI experiment. The fMRI run was 3min34s long with TR = 0.753s. The first ten seconds were discarded (dummy scans) so that we ended up with T = 284 scans collected by interleaved simultaneous multislice echo-planar images with a multi-band factor of 8 and a spatial resolution of 2 × 2 × 2mm 3 . The fMRI images were already preprocessed using a classical pipeline including realignment, coregistration, spatial normalization and smoothing (5mm isotropic) 8 . Additionally, we spatially down-sampled the fMRI images to reduce the spatial dimension of our decomposition to 61, 605 voxels.The EP was divided in two sets of motor tasks, with 15s fixation blocks at the beginning, in the middle of the acquisition and at the end of the recording. Each task set was composed of 5 blocks of 12 s each, preceded by a 3s cue indicating the task to be performed by the participant. The latter corresponded to moving the tongue, tapping the left or right finger or contracting the left or right toes. In what follows, we only consider one participant even though our results are reproducible across individuals.

  Fig. 9(e). In contrast, the neural activity recovered by the proposed multivariate

  Fig. 11(a) andFig. 11(b) show respectively the corresponding normalized maps of haemodynamic dilation parameters (δ m ) 96 m=1 in a healthy control and stroke patient, respectively. The normalization has been done by dividing

  (b) illustrates the progression of the mean accuracy score with the number of individuals involved in the LR model and clearly depicts that a plateau is reached around 459 subjects (the total size of the sampled cohort) both for the haemodynamic dilation parameter. Also, one can see the rapid progression of the mean accuracy with the number of individuals for the LR model based on haemodynamic properties compared to the one constructed from the neural activation signals.

  has emerged to set hyper-parameters. In this case, an upper-level cost function (e.g. a supervised training score on the features of the decomposition) has to be minimized with respect to the unknown hyper-parameters while staying intrinsically connected to the lower-level problem, namely the multivariate decomposition. Because of the extra-computation cost required by these approaches, such aspects are beyond the scope of this paper. Second, the proposed regional analysis is conditioned by the parcellation atlas (and the value of M ). It would be interesting to deepen this research by testing the reproducibility of the tool with some atlas variations and the creation of an atlas using subject-specific assessment of the cerebral vasculature.Third, to recover more structured spatial maps, an advanced regularization model based on TV-elastic net[START_REF] De Pierrefeu | Structured sparse principal components analysis with the TV-elastic net penalty[END_REF] or structured sparsity[START_REF] Jenatton | Multiscale mining of fMRI data with hierarchical structured sparsity[END_REF][START_REF] Baldassarre | Structured sparsity models for brain decoding from fMRI data[END_REF] could be used in space while keeping the same algorithmic structure. Recent progress in solving the TV proximity op-

Figure 4 :

 4 Figure 4: Comparison to (Wu et al., 2013) on synthetic data: (a): Firstly, we consider the case when the HRF is set to a different shape than the one used for simulating the data. (b): Secondly, we learn the HRF along with the spatio-temporal components (a/b)-top: The yellow-purple maps define the spatial activation patterns. (a/b)-bottom-left: The blue curve represents the true HRF, the green trace corresponds to the initialization and the dashed orange shape corresponds to the HRF used to deconvolve the observed BOLD signal. (a/b)bottom-right: The observed BOLD signal in the associated activation region is depicted in black, whereas the true and recovered temporal atoms are plotted in blue and orange, respectively. (a)-bottom-right: We add, in dashed green, the estimate of the temporal atoms if we choose the true HRF as the initialization. The standard deviation across voxels is encoded by transparency around mean curves. (c): We report the result of the blind-deconvolution of (Wu et al., 2013) for comparison. (c)-bottom-left:We display in green the histogram of the haemodynamic delays estimated with this concurrent approach for the activated voxels and reports the associated average delay using a green vertical line. The true haemodynamic delay in plotted in blue while the one estimated by our approach is shown in orange. (c)-bottom-right: For the voxel localized in the yellow-purple map (bottom -far right) we display the true neural activation signal in blue and the one estimated by our approach in orange along with the deconvolved BOLD signal obtained with(Wu et al., 2013) (bottom). On top of that panel, in order to quantify how accurate in time the localization of neural events is, we plot the first-order derivative of the true neural activation signal and our estimate with the same color code along with the detected events from[START_REF] Wu | A blind deconvolution approach to recover effective connectivity brain networks from resting state fMRI data[END_REF]'s method in green.

Figure 5 :Figure 6 :

 56 Figure 5: Model selection: compromise between model complexity and accuracy. (a): Evolution of the R 2 score as a function of the number of spatio-temporal atoms K in model (5) ranging from 2 to 50. (b): Evolution of the determinant of the correlation matrix Σ K between neural activation signals as a function of K ranging in the same interval as mentioned earlier.

Figure 7 :

 7 Figure 7: Single subject results from rs-fMRI semi-blind deconvolution analysis.. Top Row (a): Neural activation signal z 1 (left) and corresponding spatial map (axial view) u 1 (right), mostly involving activated voxels in the visual cortex. Second row (b): Neural activation signal z 11 (left) and corresponding spatial map (axial view) u 11 (right), mostly involving activated voxels in the default mode network (DMN). Third row, left (c): Fastest haemodynamic region. Fastest HRF estimate v δ f (left) localized in the middle temporal gyrus as shown on the parcel mask Θ f (right). Third row, right (d): Slowest haemodynamic region. Slowest HRF estimate v δs (left) localized in the frontal orbital cortex as shown on the parcel mask Θ s (right). (Bottom row (e): Voxelwise time courses. Estimate of the neural activation signal (in blue), superimposed on the denoised BOLD signal (in orange) computed as the convolution with the local HRF estimate. The observed BOLD time course in shown in black.

Figure 8 :

 8 Figure 8: Correlation matrix between neural activation signals for K = 20. Triangular inferior view of the semi-definite positive matrix Σ K for K = 20. All entries vary between -1 and +1 as they reflect correlation coefficients.

Figure 9 :

 9 Figure 9: Comparison to (Pedregosa et al., 2015) on task-fMRI data: We consider the experimental condition of the left-hand motor action. Top (a)-(b): We localize the 100 most activated voxels w.r.t. the condition. (a) corresponds to the activated voxels estimated from the thresholded regression coefficients as done in[START_REF] Pedregosa | Data-driven HRF estimation for encoding and decoding models[END_REF] and (b) corresponds to the activated voxels in our approach obtained by thresholding the spatial map and keeping the largest values. Center (c)-(d): We depict the axial slices of the corresponding activation maps for each approach ((c) refers to our competitor, (d) to our method). We observe that[START_REF] Pedregosa | Data-driven HRF estimation for encoding and decoding models[END_REF] mainly identify the activation in the cerebellum whereas our approach localizes the activation in the right motor cortex. Bottom (e): We display the neural activation time courses estimated by our approach in orange for the selected voxels. The standard deviation across voxels is encoded by transparency around the mean curve. Similarly, we report the observed BOLD signal in blue and display the temporal profile of the motor task in light gray. Bottom (f ): We display the HRF estimates for each method. The green HRF corresponds to[START_REF] Pedregosa | Data-driven HRF estimation for encoding and decoding models[END_REF] with the between-voxel standard deviation encoded by transparency around the mean curve. Similarly, we depict our HRF estimate in orange. We notice larger haemodynamic delays estimated by our method.

Figure 10 :

 10 Figure 10: Within-vs between-subject analysis of the haemodynamic variability. The box plots show respectively in blue, orange and green the distribution of WS( δ si T1 , δ si T2 ) for all subjects (i = 1, . . . , 100), BS( δ si T1 , δ sj T1 ) and BS( δ si T2 , δ sj T2) with i = j. These distributions are assessed for 6 levels of temporal regularization (λ f ∈ {10 -6 , 10 -3 , 10 -2 , 10 -1 , 5 10 -1 , 1 -10 -6 }) and remain stable. Statistical analysis (paired t-test) was conducted to assess the significance of the difference between the mean of the within-and between-subject 2 2 norm distributions. Significant differences are marked with a * .

Figure 11 :

 11 Figure 11: Haemodynamic discrimination between stroke patients (SP) and healthy controls (HC). Top (a)-Middle (b): Normalized haemodynamic dilation parameter maps in a healthy control (top row) and stroke patient (bottom row), respectively. The maps have been respectively normalized by the within-subject mean value ( δs = 1 M

Figure 13 :

 13 Figure 13: Accuracy score for classifying middle-age vs. elderly subjects (459 individuals sampled from the UK Biobank database). (a): The prediction was performed by pulling individual features either based on (i) the estimated neural activation signals ( z k ) 20k=1 , (ii) HRF shape estimates ( v δm ) 96 m=1 or (iii) the haemodynamic dilation parameters ( δ) 96 m=1 . From top to bottom, the distribution of the classification scores is shown from for the predictive features (i)-(iii), respectively. The best accuracy scores (average 0.74) are reached using the haemodynamic parameters and the smallest variability in the prediction using specifically the dilation parameter estimates. (b): Learning curve of accuracy scores as a function of the number of individuals (middle-age vs elderly subjects) used for the training stage both for the haemodynamic dilation parameter (blue curve) and the temporal components (orange curve). As a plateau is reached for 459 people, we presented the corresponding performances in panel (a).

Table 1 :

 1 Taxonomy of functional networks involved in the spatio-temporal decomposition (5) with K = 20, whose spatial maps are shown in Fig.6.

	Networks	#Comp. id k
	Default Mode	6, 11, 18, 20
	Attention	5
	Central executive	10, 12, 13
	Salience	19
	Sensori-motor	3, 7, 9, 15, 17
	Auditory-language	2, 8, 14
	Visual	1, 4, 16
	sensory (visual, auditory and motor) and resting-state networks (RSN) are quite
	well retrieved by multiple components, localized respectively in the occipital (com-
	ponents 1, 4, 16) temporal (components 2, 8 and 14), fronto-parietal (components
	5, 10, 11, 13, 17 and 20), frontal/fronto-polar (components 6, 12, 18 and 19) and

https://pypi.org/project/prox tv/

This function is zero-valued inside the constraint set and equals infinity elsewhere.

code available at https://github.com/hemolearn/hemolearn.

Acquisition details can be found at https://www.fmrib.ox.ac.uk/ukbiobank/protocol/.

Preprocessing details can be found at https://biobank.ctsu.ox.ac.uk/crystal/ crystal/docs/brain_mri.pdf

no statistical test performed at the individual level

https://github.com/hemolearn/hemolearn

Acquisition details can be found at https://humanconnectome.org/storage/app/media/ documentation/s1200/HCP_S1200_Release_Appendix_VI.pdf

Acknowledgment

This work was supported by a CEA PhD scholarship, the UK Royal Academy of Engineering under the RF/201718/17128 grant and the SRPe PECRE 1718/15 Award. We would like to thank our colleagues Bertrand Thirion, Kamalaker Reddy Dadi and Thomas Bazeille from Inria for fruitful discussions that helped us investigate the proposed data analyses. Finally we are grateful to Cécilia Garrec (CEA/DRF-Dir) for her comments and suggestions that have significantly improved the manuscript.