Salima Makni

Jérôme Idier

Thomas Vincent

Bertrand Thirion

Ghislaine Dehaene-Lambertz

Philippe Ciuciu
email: philippe.ciuciu@cea.fr

A fully Bayesian approach to the parcel-based detection-estimation of brain activity in fMRI

Keywords: Bayesian modelling, fMRI, gamma-Gaussian mixture model, detection-estimation, Markov Chain Monte Carlo methods, Bayes factor, model comparison

. Then, in each parcel, prior information is embedded to constrain this estimation. Detection is achieved by modelling activating, deactivating and non-activating voxels through mixture models within each parcel. From the posterior distribution, we infer upon the model parameters using Markov Chain Monte Carlo (MCMC) techniques. Bayesian model comparison allows us to emphasize on artificial datasets first that inhomogeneous gamma-Gaussian mixture models outperform Gaussian mixtures in terms of sensitivity/specificity trade-off and second that it is worthwhile modelling serial correlation through an AR(1) noise process at low signal-to-noise (SNR) ratio. Our approach is then validated on an fMRI experiment that studies habituation to auditory sentence repetition. This phenomenon is clearly recovered as well as the hierarchical temporal

Introduction

Since the first report of the BOLD effect in human [START_REF] Ogawa | Brain magnetic resonance imaging with contrast dependent on blood oxygenation[END_REF], functional magnetic resonance imaging (fMRI) has represented a powerful tool to non-invasively study the relation between cognitive stimulus and the hemodynamic (BOLD) response. Within-subject analysis in fMRI is usually addressed using a hypothesis-driven approach that actually postulates a model for the HRF and enable voxelwise inference in the General Linear Model (GLM) framework. In this formulation, the modelling of the BOLD response i.e., the definition of the design matrix is crucial. In its simplest form, this matrix relies on a spatially invariant temporal model of the BOLD signal across the brain meaning that the expected response to each stimulus is modelled by a single regressor. Assuming the neurovascular system as linear and timeinvariant (LTI), this regressor is built as the convolution of a sparse spike train representing the stimulation signal and the canonical HRF, i.e., a composition of two gamma functions which reflects the BOLD signal best in the visual and motor cortices [START_REF] Glover | Deconvolution of impulse response in event-related BOLD fMRI[END_REF].

Intra-individual differences in the characteristics of the HRF have been exhi--bited between cortical areas in [START_REF] Aguirre | The variability of human BOLD hemodynamic responses[END_REF][START_REF] Miezin | Characterizing the hemodynamic response: effects of presentation rate, sampling procedure, and the possibility of ordering brain activity based on relative timing[END_REF][START_REF] Neumann | Bayesian second-level analysis of functional magnetic resonance images[END_REF][START_REF] Handwerker | Variation of BOLD hemodynamic responses across subjects and brain regions and their effects on statistical analyses[END_REF]. Although smaller than inter-individual fluctuations, this regional variability is large enough to be regarded with care. To account for these spatial fluctuations at the voxel level, one usually resorts to hemodynamic function basis. For instance, the canonical HRF can be supplemented with its first and second derivatives to model differences in time [START_REF] Friston | Imaging neuroscience: Principles or maps?[END_REF][START_REF] Henson | Detecting latency differences in event-related BOLD responses: application to words versus nonwords and initial versus repeated face presentations[END_REF]. To make the basis spatially adaptive, Woolrich et al. (2004a) have proposed a half-cosine parameterisation in combination to the selection of the best basis set. Although powerful and elegant, the price to be paid for such a flexible modelling lies in a loss of sensitivity of detection: the larger the number of regressors in the basis, the smaller the number of effective degrees of freedom in any subsequent statistical test. Crucially, in a GLM involving several regressors per condition, the Student-t statistic can no longer be used to infer on differences between experimental conditions. Rather, an unsigned Fisher statistic has to be computed, making direct interpretation of activation maps more difficult. Indeed, the null hypothesis is actually rejected whenever any of the contrast components deviates from zero and not specifically when the difference of the response magnitudes is far from zero.

In this paper, to facilitate cognitive interpretations, we argue in favour of a spatially adaptive GLM in which a local estimation of the HRF is performed. This allows us to factorise the expected BOLD response with a single regressor attached to each experimental condition and to enforce direct statistical comparisons based on response magnitudes. However, to conduct the analysis in an efficient and reliable manner, local estimation is performed at the scale of several voxels.

As mentioned earlier, the localisation of brain activation strongly depends on the modelling of the brain response and thus of its estimation. Of course, the converse also holds: HRF estimation is only relevant in voxels that elicit signal fluctuations correlated with the paradigm. Hence, detection and estimation are intrinsically linked to each other. The key point is therefore to tackle the two problems in a common setting, i.e., to set up a formulation in which detection and estimation enter naturally and simultaneously. This setting cannot be the classical hypothesis testing framework. Indeed, the sequential procedure which first consists in estimating the HRF on a given dataset and then building a specific GLM upon this estimate for detecting activations in the same dataset, entails statistical problems in terms of sensitivity and specificity: the control of the false positive rate actually becomes hazardous due to the use of an erroneous number of degrees of freedom. We rather propose a Bayesian approach that provides an appropriate framework to address both detection and estimation issues in the same formalism.

The literature on Bayesian fMRI methods offers several approaches to adequately choose priors for detection. As introduced in [START_REF] Everitt | Mixture model mapping of brain activation in functional magnetic resonance images[END_REF][START_REF] Vaever Hartvig | Spatial mixture modeling of fMRI data[END_REF]Penny and Friston, 2003) and further developed in [START_REF] Smith | Assessing brain activity through spatial Bayesian variable selection[END_REF][START_REF] Woolrich | Mixture models with adaptive spatial regularization for segmentation with an application to fMRI data[END_REF][START_REF] Ou | From spatial regularization to anatomical priors in fMRI analysis[END_REF][START_REF] Woolrich | Variational Bayes inference of spatial mixture models for segmentation[END_REF][START_REF] Flandin | Bayesian fMRI data analysis with sparse spatial basis function priors[END_REF], prior mixture models define an appropriate way to perform the classification or the segmentation of statistical parametric maps into activating, non-activating or deactivating brain regions. The pioneering contributions related to mixture modelling in fMRI [START_REF] Everitt | Mixture model mapping of brain activation in functional magnetic resonance images[END_REF][START_REF] Vaever Hartvig | Spatial mixture modeling of fMRI data[END_REF] have proposed a voxel by voxel classification to decide whether the estimated effect is analogous to signal or noise in each voxel. Yet, the use of mixture modelling in a joint detection-estimation problem introduces specific concerns in comparison to the usual "hypothesis testing framework". Indeed, our data are not the voxelwise z-statistics but rather the raw fMRI time courses, which are required for the estimation step.

As regards HRF estimation, various priors may be thought of depending on the underlying HRF model. Basically, three classes of models coexist. Parametric models appeared first in the literature [START_REF] Friston | Statistical parametric mapping[END_REF][START_REF] Lange | Empirical and substantive models, the Bayesian paradigm, and meta-analysis in functional brain imaging[END_REF][START_REF] Cohen | Parametric analysis of MRI data using linear systems methods[END_REF][START_REF] Rajapakse | Modeling hemodynamic response for analysis of functional MRI time-series[END_REF][START_REF] Kruggel | Modeling the hemodynamic response in single-trial functional MRI experiments[END_REF]. In this setting, the estimation problem consists in minimising some criterion with respect to (w.r.t.) some parameters of a precise function (e.g., Gaussian, gamma ,...). However, parametric models tend to introduce some bias in the HRF estimate, since it is unlikely that they capture the true shape variations of the brain dynamics. Moreover, the objective function to be minimized is often non-convex making the parameter estimates unreliable and dependent of the initialisation. Hence, more flexible semi-parametric approaches have been proposed later to capture these variations [START_REF] Genovese | A Bayesian time-course model for functional magnetic resonance imaging data (with discussion)[END_REF][START_REF] Gössl | Bayesian spatio-temporal modeling of the hemodynamic response function in BOLD fMRI[END_REF]Woolrich et al., 2004a). In a semi-parametric framework, the HRF time course is decomposed into different periods (initial dip, attack, rise, decay, fall, ...), each of them being described by specific parameters. In the same time, non-parametric approaches or Finite Impulse Response (FIR) models have emerged in the fMRI literature as a powerful tool to infer on the HRF shape [START_REF] Nielsen | Comparison of two convolution models for fMRI time series[END_REF][START_REF] Goutte | Modeling the haemodynamic response in fMRI using smooth filters[END_REF][START_REF] Marrelec | Robust Bayesian estimation of the hemodynamic response function in event-related BOLD MRI using basic physiological information[END_REF][START_REF] Ciuciu | Unsupervised robust non-parametric estimation of the hemodynamic response function for any fMRI experiment[END_REF][START_REF] Marrelec | Estimation of the hemodynamic response function in event-related functional MRI: Bayesian networks as a framework for efficient Bayesian modeling and inference[END_REF]. Most of these works take place in the Bayesian setting and constrain the HRF to be temporally smooth, which warrants a stable estimation in case of ill-posed identification.

Whatever the model in use, most methods are massively univariate and therefore neglect the spatial structure of the BOLD signal. Early investigations have shown that estimating the HRF using regularised FIR models over a functionally homogeneous region of interest provides more reliable results [START_REF] Gössl | Bayesian spatio-temporal modeling of the hemodynamic response function in BOLD fMRI[END_REF][START_REF] Ciuciu | Outlier detection for robust region-based estimation of the hemodynamic response function in event-related fMRI[END_REF]). In the following, a region-based non-parametric model of the HRF is therefore adopted. Then, the critical issue arising consists in exhibiting a functionally homogeneous clustering of the fMRI datasets over the whole brain. To that end, the grey matter's mask is segregated into a few hundreds of connected Regions of Interest (ROIs), called parcels. The parcels are derived using the parcellation procedure proposed by [START_REF] Thirion | Dealing with the shortcomings of spatial normalization: Multi-subject parcellation of fMRI datasets[END_REF], according to a compound criterion balancing spatial and functional homogeneity. The second step of our analysis solves for the detection-estimation problem over each parcel.

The rest of this paper is organised as follows. Section 2 details how anatomical information is handled and how parcels are built up. Then, the forward parcelbased model of the BOLD signal is derived and priors over the unknown parameters are specified. In Section 3, we explain the key steps of our inferential procedure based on MCMC methods, posterior mean (PM) HRF estimation and marginal Maximum A Posteriori (MAP) classification for detection. On artificial datasets, Section 4 reports the performance of our approach in terms of sensitivity-specificity trade-off depending on the mixture prior and the noise modelling. In Section 5, our joint detection-estimation approach is tested on real fMRI data acquired during an habituation study to auditory sentence repetition. On the same datasets, we also performed a classical GLM analysis employing the widely used Statistical parametric mapping (SPM) software1 . The two approaches are then compared and the main differences are exhibited. The pros and cons of the proposed method are discussed in Section 6 and some future extensions are envisaged.

Methodology

Definition of functionnally homogeneous brain regions

Anatomical representation

The segmentation of the grey-white matter interface is performed on an anatomical T 1 -weighted MRI image using the BrainVisa software2 [START_REF] Mangin | From 3D magnetic resonance images to structural representations of the cortex topography using topology preserving deformations[END_REF]. It provides us with the anatomical representation of the cortex. To accommodate the coarser spatial resolution of fMRI data (typically, 3.5mm along each direction), the grey matter mask M a is dilated using a sphere as structural element, with a radius equal to the resolution of functional images.

Concurrently, a functional mask M f was computed from the motion-corrected3 BOLD EPI volumes. Also, an average EPI volume was created. Then, we carried out a histogram analysis of this volume: a Gaussian density N (µ, σ 2 ) was fitted on the main mode of the EPI signal of interest. A threshold defined as µ -3σ was used to obtain the functional mask. Finally, the mask of interest where activation most likely occurs was built as M s = M a ∩ M f .

Parcellation of the grey-matter

The volume in mask M s was then divided in K functionally homogeneous parcels or ROIs using the parcellation technique proposed in [START_REF] Thirion | Dealing with the shortcomings of spatial normalization: Multi-subject parcellation of fMRI datasets[END_REF]. The goal of this procedure is to segregate the brain into connected and functionally homogeneous components. For doing so, the parcellation algorithm relies on the minimisation of a compound criterion reflecting both the spatial and functional structures and hence the topology of the dataset. The spatial similarity measure favours the closeness in the Talairach coordinates system. The functional part of this criterion is computed on parameters that characterise the functional properties of the voxels. These parameters can be chosen either as the fMRI time series themselves or as the β-parameters estimated during a first-level SPM analysis. The latter choice is nothing but a projection onto a subspace of reduced dimension, i.e., the feature space. Typically, the feature space is defined from a F-contrast in a SPM analysis.

The number of parcels K needs to be set by hand. The larger the number of parcels, the higher the degree of within-parcel homogeneity. Of course, there exists a trade-off between the within-parcel homogeneity and the signal-tonoise ratio (SNR). If the number of voxels is too small in a given parcel, the HRF estimation may become inaccurate, specifically in regions where no voxel elicits a specific response to any experimental condition. To objectively choose an adequate number of parcels, [START_REF] Thyreau | Anatomofunctional description of the brain: a probabilistic approach[END_REF] have used the Bayesian information criterion (BIC) and cross validation techniques on an fMRI study of ten subjects. They have shown converging evidence for K ≈ 500 for a whole brain analysis and recommend K = 200 as a fair setting for a restricted analysis to the grey matter's mask leading to typical parcel sizes around a few hundreds voxels.

Parcel-based modelling of the BOLD signal

Vectors and matrices are displayed in lower and upper cases, respectively, both in bold font (e.g., y and P ). Unless stated otherwise, subscripts i, j, m and n are respectively indexes over mixture components, voxels, stimulus types and time points. We refer the reader to appendix A for the definitions of the non-standard probability density functions (pdf). Also, the pdf families are denoted using calligraphic letters (e.g., N and G for the Gaussian and gamma densities).

The regional forward model of the BOLD signal introduced in [START_REF] Makni | Joint detectionestimation of brain activity in functional MRI: a multichannel deconvolution solution[END_REF] is used to account for voxel-dependent and task-related fluctuations of the magnitudes of the BOLD response. Hereafter, the latter magnitudes are called Neural Response Levels (NRLs). In short, this time-invariant model characterises each and every parcel by a single HRF shape and a NRL for each voxel and stimulus type. As shown in Fig. 1, this means that although the HRF shape is assumed constant within a parcel, the magnitude of the activation can vary in space and across experimental conditions. Let P = (V j ) j=1:J be the current parcel and V j a voxel in P. Then, the generative BOLD model reads:

y j = M m=1 a m j X m h + P j + b j , (1)
where

• y j = (y j,t n ) n=1:N denotes the BOLD fMRI time course measured in voxel V j at times (t n ) n=1:N (N is the number of scans) with t n = nTR and TR is the time of repetition, • X m = (x m t n -d∆t ) n=1:N,d=0:D is a N × (D + 1) binary matrix corresponding to the arrival times for the mth condition. ∆t is the sampling period of the HRF, usually smaller than TR. The onsets of the stimuli are put on the ∆t-sampled grid by moving them to the nearest time points on this grid. Note that X m can be adapted to paradigms having trial varying stimulus magnitudes or durations.

• Vector h = (h d∆t ) d=0:D represents the unknown HRF shape in parcel P (D+ 1 is the number of HRF coefficients). It actually seems reasonable to assume a single HRF shape in homogeneous parcels. • a m j stands for the NRL in voxel V j for condition m (M is the number of experimental conditions in the paradigm). Hence, the activation profile associated to the mth stimulus type in voxel V j is computed as the product h × a m j .

• P = [p 1 , . . . , p Q ] is the low frequency orthogonal matrix of size N × Q. It
consists of an orthonormal basis of functions p q = (p q (t n )) n=1:N . To each voxel is attached an unknown weighting vector j that has to be regressed in order to estimate the trend in V j . We denote l = ( j ) j=1:J the set of low frequency drifts involved in P. • b j ∈ R N is the noise vector in voxel V j . In [START_REF] Woolrich | Temporal autocorrelation in univariate linear modelling of fMRI data[END_REF][START_REF] Worsley | A general statistical analysis for fMRI data[END_REF] an autoregressive (AR) noise model has been introduced to account for the serial correlation of the fMRI time series in the detection analysis. Importantly, when this temporal correlation is correctly estimated, the number of false positives decreases, yielding more conservative detection results. Similarly, in a joint detection-estimation framework, [START_REF] Makni | Joint detectionestimation of brain activity in fMRI using an autoregressive noise model[END_REF] have shown that the introduction of a spatially-varying first-order AR noise in model (1) provides a lower false positive rate. Hence, b j is defined by b j,tn = ρ j b j,t n-1 + ε j,tn , ∀j, t, with ε j ∼ N (0 N , σ 2 ε j I N ), where 0 N is a null vector of length N , and I N is the identity matrix of size N .

Although the noise structure is correlated in space (Woolrich et al., 2004b) and non-stationary across tissues [START_REF] Worsley | A general statistical analysis for fMRI data[END_REF], we do not essentially account for this correlation for two reasons. First, it is likely that a large part of the noise may be due to misspecification of the HRF. Second, we actually assume that the spatial correlation of the signal of interest is more important. Hence, the fMRI time series y = (y j ) j=1:J are supposed to be statistically conditionally independent and identically distributed (iid). The Fig. 1. Summary of the proposed regional BOLD model. The size of each parcel P varies typically between by a few tens and a few hundreds: 80 J 350. The number M of experimental conditions involved in the model usually varies from 1 to 5. In our example, M = 2, the NRLs (a 1 j , a 2 j ) corresponding to the first and the second conditions are surrounded by circles and squares, respectively. Note that our model accounts for asynchronous paradigms in which the onsets do not necessarily match acquisition time points. As illustrated, the NRLs take different values from one voxel to another. The HRF h can be sampled at a period of 1s and estimated on a range of 20 to 25s (e.g., D = 25). Most often, the LFD coefficients j are estimated on a few components (Q = 4). likelihood then factors over voxels:

p(y | h, a, l, θ 0 ) = J j=1 p(y j | h, a j , j , θ 0,j ) ∝ J j=1 |Λ j | 1/2 σ -N ε j exp - J j=1 ỹt j Λ j ỹj 2σ 2 ε j (2)
where θ 0,j = (ρ j , σ 2 ε j ), θ 0 = (θ 0,j ) j=1:J and ỹj = y jm a m j X m h -P j . Note that σ -2 ε j Λ j defines the inverse of the autocorrelation matrix of b j . According to Kay (1988, Chap VI, p. 177)

, Λ j is tridiagonal, with (Λ j ) 1,1 = (Λ j ) N,N = 1, (Λ j ) n,n = 1 + ρ 2 j and (Λ j ) n+1,n = (Λ j ) n,n+1 = -ρ j ∀ n = 2 : N -1. Its determinant is given by |Λ j | = 1-ρ 2 j .
In what follows, we do not approximate Eq. ( 2) by dropping the term |Λ j | 1/2 , as done in previous works [START_REF] Roberts | Variational Bayes for generalized autoregressive models[END_REF]Penny et al., 2003;Woolrich et al., 2004b). Indeed, when the AR parameter ρ significantly departs from zero (e.g., ρ ≥ 0.4), this approximation is biased and potentially far from the exact likelihood.

On the sole basis of the likelihood function (2), it seems impractical to identify the pair (h, a). Indeed, Maximum likelihood (ML) estimation of (h, a) is a bilinear inverse problem since (1) is linear w.r.t. h when a is fixed and viceversa. Therefore, the ML solution (h * , a * ) is not unique. For instance, every couple (h * /s, a * ×s) defines another pair of solutions in the ML sense whatever the scale parameter s > 0. To get rid of identifiability problems and reach a more reliable estimation, in the Bayesian formalism we introduce suitable prior distributions attached to the unknown quantities (h, a).

Priors

The Hemodynamic response function

Akin to [START_REF] Buxton | A model for the coupling between cerebral blood flow and oxygen metabolism during neural stimulation[END_REF][START_REF] Goutte | Modeling the haemodynamic response in fMRI using smooth filters[END_REF][START_REF] Marrelec | Robust Bayesian estimation of the hemodynamic response function in event-related BOLD MRI using basic physiological information[END_REF], the HRF is characterised as follows: (i) its variations are smooth; (ii) it is causal and returns to a baseline after a given time interval T (h t = 0, ∀ t < 0 and t > T ); T is fixed by the user according to the experimental paradigm (usually 25 seconds).

Condition (i) may be fulfilled using an approximation of the second-order derivative ∂ 2 h 2 :

(∂ 2 h) dτ ≈ (h (d+1)τ -2h dτ + h (d-1)τ )/τ 2 , ∀ d = 1 : D -1.
In matrix form, we get ∂ 2 h = D 2 h. Condition (ii) is ensured with a HRF h whose magnitude vanishes at first and last time points (h 0 = h D = 0). Hence, D 2 is the truncated second-order finite difference matrix of size (D -1) × (D -1) and

∂ 2 h 2 = h t R -1 h with R = (D t 2 D 2
) -1 a symmetric positive definite matrix. The prior on h thus reads h ∼ N (0, σ 2 h R). To overcome the scale ambiguity problem mentioned earlier, we constrain the HRF to be of unitary norm ( h = 1). Alternative constraints such as setting the value of the peak could be considered.

The "neural" response levels

Mixture models are often used as a second stage to segment the SPMs (i.e., the statistical maps) resulting from a first-level temporal analysis of fMRI time series [START_REF] Vaever Hartvig | Spatial mixture modeling of fMRI data[END_REF][START_REF] Everitt | Mixture model mapping of brain activation in functional magnetic resonance images[END_REF][START_REF] Woolrich | Mixture models with adaptive spatial regularization for segmentation with an application to fMRI data[END_REF]. This means that the data to be classified correspond to some normalised effect c t β/ std (c t β), where vector c defines a contrast of interest (typically a comparison between two experimental conditions) and β is the vector of parameter estimates after fitting a GLM against the fMRI data.

In the present paper, as well as in [START_REF] Makni | Joint detectionestimation of brain activity in functional MRI: a multichannel deconvolution solution[END_REF], prior mixture models are used in a different way, closer to that proposed by [START_REF] Svensen | Probabilistic modeling of single-trial fMRI data[END_REF]. In the same spirit, a mixture model is introduced on the NRLs for every experimental condition m and not specifically on the linear combination c t β.

Table 1

Model definition and notations. Here, AR(1) stands for a first-order autoregressive noise in time whose parameters vary in space. In this respect, it is referenced as a spatially-varying AR(1) noise. The second noise model under study is a spatiallyvarying white noise. The columns describe the different NRLs priors: GMM stands for a two-class Gaussian mixture model (a zero-mean Gaussian density (G) for nonactivating voxels and a Gaussian density (G) for activating voxels). GaGMM stands for a two-class gamma-Gaussian mixture model (a centred Gaussian density (G) for non-activating voxels and a gamma density (Ga) for activating voxels). GaGGaMM stands for a three-class mixture model composed of a zero-mean Gaussian density and two gamma densities (a gamma density (Ga) for activating voxels and a flipped gamma density (Ga) for deactivating voxels).

GaGMM GMM GaGGaMM AR(1) noise M 1 M 3 white noise M 2 M 4 M 5
In (Makni et al., 2006a), it was stressed that a two-class Gaussian mixture model (GMM) may be inadequate for segregating noise from true activations.

In particular, it can be shown that this kind of independent mixture may degenerate in the sense that the two probability density functions (pdf) overlap almost entirely if there are not enough activating voxels in the current parcel (see (Makni et al., 2005, §VII.)). For this reason, we have rather adopted an inhomogeneous prior mixture model. Among several possibilities (Gaussianlognormal MM, Gaussian-truncated Gaussian MM, ...), a gamma-Gaussian mixture model (GaGMM) has been retained for technical reasons that will become clearer in what follows. The non-activating voxels are still modelled using a zero-mean Gaussian pdf while a gamma distribution is used to enforce positivity of activating voxels. Akin to [START_REF] Vaever Hartvig | Spatial mixture modeling of fMRI data[END_REF][START_REF] Woolrich | Mixture models with adaptive spatial regularization for segmentation with an application to fMRI data[END_REF], a three-class mixture prior model is actually considered to account for deactivating voxels. Since we assume that deactivation corresponds to a negative BOLD response, we use a flipped gamma density defined on the left part of the real line leading to define the GaGGaMM extension (see Table 1).

In our model, different stimulus types are supposed to induce statistically independent hemodynamic magnitudes or NRLs i.e., p(a

| θ a ) = m p(a m | θ m )
with a = (a m ) m=1:M , a m = (a m j ) j=1:J and θ a = {θ 1 , . . . , θ m }. Vector θ m denotes the set of unknown hyper-parameters related to the mth stimulus type. Because our mixture model is voxelwise, the prior pdf factors over voxels:

p(a m | θ m ) = j p(a m j | θ m ).
Importantly, the hyper-parameters are kept constant for all voxels in a given parcel because of the within-parcel homogeneity. These parameters may actually vary from one parcel to another. Let q m j be the allocation variable (the label) that indicates whether voxel V j is activating (q m j = 1), deactivating (q m j = -1) or non-activating (q m j = 0) in condition m.

The marginal density p(a m j | θ m ) thus reads:

p(a m j | θ m ) = 1 i=-1 Pr(q m j = i | λ m )f (a m j | q m j = i, θ m ) = 1 i=-1 λ i,m f i (a m j | θ m ), (3)
with

λ m = (λ -1,m , λ 0,m , λ 1,m ) and f 0 (a m j | θ m ) = N (0, v 0,m ), f ±1 (a m j | θ m ) = G(α ±1,m , β ±1,m
). The λ i,m parameters define the prior probabilities of the threeclass mixture on the NRLs ( i λ i,m = 1). For instance, λ 1,m gives us the prior probability of being activated in response to condition m. Since the mixture is independent in space, we have λ 1,m = Pr(q m j = 1 | λ m ), ∀j. Note that q m j | λ m follows a multinomial distribution over the 3-dimensional probability simplex, i.e., q m j ∼ MN 3 (1; λ m ) (see Appendix A). Hence, seven hyper-parameters are necessary to describe the prior mixture for each experimental condition m:

θ m = {λ ±1,m , α ±1,m , β ±1,m , v 0,m }.
Compared to [START_REF] Woolrich | Mixture models with adaptive spatial regularization for segmentation with an application to fMRI data[END_REF], we set the mean of the non-activating class to zero (µ 0,m = 0, ∀ m), while we do not need to place restrictions on the mode of the activation and deactivation gamma classes.

The nuisance variables

We assume that l is a random process independent of h such that p(l; σ 2 ) = j p( j ; σ 2 ) and j ∼ N (0, σ 2 I Q ).

The hyper-parameters

All the hyper-parameters are concatenated into the overall parameter vector Θ = {θ 0 , σ 2 h , σ 2 , θ a }. Without informative prior knowledge, the following priors are retained for (σ 2 h , σ 2 , θ 0 ):

p(σ 2 h , σ 2 ) = (σ h σ ) -1 , p(θ 0 ) = J j=1 p(ρ j , σ 2 ε j ) = J j=1 σ -1 ε j 1 (-1,1) (ρ j ), (4)
to ensure stability of the AR(1) noise process [START_REF] Kay | Modern Spectral Estimation[END_REF].

Mixture parameters. As regards variances v 0,m , an improper Jeffreys' prior p(v 0,m ) = v -1/2 0,m is considered because we do expect non-activating voxels in a given parcel. Hence, class 0 should never be empty a priori. However, to avoid emptiness and subsequent degeneracy problems making the sampling of the posterior distribution of v 0,m unfeasible, a conjugate prior could also be chosen, that is, an inverse gamma density IG(v 0,m , a v 0 , b v 0 ), where (b v 0 , c v 0 ) are fixed values chosen in an appropriate way to make the prior flat enough.

The non-negativity of parameters α i,m is guaranteed through the use of an exponential density E(α i,m ; s i ) ≡ G(α i,m ; 1, s i ) as prior distribution (see Appendix A). For parameters β i,m we resort to the conjugate prior, given by a gamma density G(β i,m ; b i , c i ) for i = ±1.

Mixture probabilities. As regards mixture probabilities λ m ∈ [0, 1] 3 , a Dirichlet prior distribution is used as it is conjugate to the multinomial distribution used for labels, i.e., MN 3 (q m j | λ m ). More exactly, a symmetric Dirichlet

density D 3 (λ m | δ) is selected with δ = δ1 3 and δ > 0 (see Appendix A).
The full prior density p(θ m ) thus reads:

p(θ m ) = v -1/2 0,m Γ(3δ) 3Γ(δ) i=±1 λ δ-1 i,m s i c b i i Γ(b i ) β b i -1 i,m exp(-s i α i,m -c i β i,m ). (5)
Values of (a ±1 , b ±1 , c ±1 , s ±1 , δ) are fixed empirically but do not really influence the results in most cases4 . These parameters make the sampling steps of (α ±1,m , β ±1,m ) always possible even when one of the two classes ±1 is empty, because the hyper-prior densities have been chosen proper.

The full posterior distribution

Combining data-driven information in each parcel with prior knowledge using Bayes' rule, we get the full posterior distribution, which is the keystone both for localising activations and deactivations as well as for estimating the corresponding parcel-based HRF:

p(h, a, l, Θ | y) ∝ p(y | h, a, l, θ 0 ) p(a | θ a ) p(h | σ 2 h ) p(l | σ 2 ) p(Θ) ∝ σ -D h σ -JQ J j=1 (1 -ρ 2 j ) 1/2 σ N +1 ε j 1 (-1,1) (ρ j ) × • • • × exp - h t R -1 h 2σ 2 h - J j=1 1 2σ 2 ε j ỹt j Λ j ỹj + 1 2σ 2 j 2 × • • • × M m=1 p(θ m ) J j=1 p(a m j | θ m ) (6)
where p(a m j | θ m ) and p(θ m ) are defined by ( 3) and ( 5), respectively.

Note that the parcel-based HRF h can be identified if at least one voxel elicits activation in response to one or several experimental conditions involved in model( 1). In addition, other identifiability problems may occur on hyperparameters such as the mean and variance parameters. It is necessary that at least two voxels belong to each class in order to properly estimate the variances attached to the mixture components. In practice, there is no numerical problem because of the choice of proper priors for the hyper-parameters; see Subsection 3.1.2 for practical details.

Inference scheme

Our objective is to obtain an estimate of the joint posterior distribution of all unknown parameters, given the observed data. Exact and analytical approaches are not feasible with non-Gaussian models such as ( 6). Several competing inferential schemes are possible. For instance, approximations to the full posterior distribution can be derived in the Variational Bayes (VB) framework or using Taylor series expansion. In our context, given the bilinear structure of the generative BOLD model ( 6), the VB formulation would be feasible only at the expense of separability assumptions between a and h in the approximation of the posterior distribution. Further work is required to decide whether or not this hypothesis is tenable. Instead, we resort to a more computationally demanding but exact approach to simulate realisations of the full posterior distribution.

Gibbs sampling algorithm

To draw realisations of the full posterior distribution, a Gibbs sampler is implemented. This consists in building a Markov chain, whose stationary distribution is the joint posterior pdf (6), by sequentially generating random samples from the full conditional pdfs of all the unknown parameters and hyper-parameters ; see [START_REF] Liu | Monte Carlo strategies in scientific computing[END_REF][START_REF] Robert | The Bayesian Choice[END_REF] for a general introduction to MCMC.

As shown in appendix B, direct sampling according to the full conditional distributions is only feasible for the HRF h, the labels q, the nuisance variables l, the noise variances σ ε , the mixture probabilities λ m , and part of the hyper-parameters (scales σ h and σ , class 0 variances v 0 and shape parameters β i for i = ±1). In contrast, direct simulation is not tractable for the other parameters, i.e., the NRLs a corresponding to classes ±1, the AR parameters ρ and the scale parameters α i of the gamma densities for i = ±1. Therefore, single-component Metropolis-Hastings jumps [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF] are specifically designed. More precisely, separate jumps are proposed for each of the parameters in turn. To this end, suitable instrumental distributions regarding the parameters of interest are designed (see Appendix B for details).

Initialisation

Parameters are uniformly initialised. This means that we set up all voxels with the same noise statistical parameters (θ 0,j = θ 0 , ∀j) and that we use the same starting values of mixture hyper-parameters (θ m = θ , ∀m). In the first parcel, the HRF is initialised to the canonical shape [START_REF] Glover | Deconvolution of impulse response in event-related BOLD fMRI[END_REF]. In the next ones, the HRF is set up using the mean of the estimates computed over the already processed neighbouring parcels. We resort to the same strategies for the labels and the corresponding NRLs when the parcel sizes match approximately. We have checked that this strategy provides shorter burn-in periods5 and thus reduces the computation load.

identifiability issues

To cope with these identifiability problems, we have carried out the following three steps procedure over the first iterations of our MCMC algorithm:

• initialise each parcel-specific HRF with a fixed shape in order to obtain a first estimate of labels q; • check that the class of activating voxels is effectively not empty for at least one experimental condition:

--If ∃ m ∈ N * M = {1, .
. . , M } such that in every parcel P ∃ j ∈ {V 1 , . . . , V J } ∈ P | q m j = 1 then release the HRF constraint to estimate the complete model i.e., including the HRF shape; --otherwise, discard the current parcel: the HRF estimate is not reliable in P. Since the corresponding NRLs are close to zero in that case there is no evoked activation due to the experimental paradigm.

Convergence diagnosis

We use a burn-in period of 500 iterations, followed by 1000 subsequent jumps and compute PM and MAP estimates every two jumps. Observations of the chain with different initial conditions confirmed that a burn-in of 500 jumps was sufficient. In addition, convergence has been checked by monitoring online the behaviour of the estimated values of some scalar parameters (e.g., noise variances, AR parameters, ...) from one iteration to another. These observations confirmed also that 1000 iterations were sufficient.

Computational load and parallel implementation

Our current implementation (PyHRF package) is in Python, while the most intensive computations (e.g., computation of the inverse covariance matrix of h) have been coded in C-language and interfaced with the Gnu Scientific Library (GSL)6 . This allows us to take advantage of a parallel computing system available through the Seppo library (Simple Embarassingly Parallel Python7 ) and the Pyro (Python Remote Object) server. Using such a system, all the parameter estimates are obtained in about 2 mn for a parcel of mean size (250 voxels) for two experimental conditions (M = 2). Since about 200 parcels are necessary to cover the grey matter's mask, a complete withinsubject analysis takes about 2 hours when running four processes on a dual core bi-processors Pentium IV (2.7 GHz). PyHRF will be available in the next release of BrainVisa8 in March 2008.

Derivation of parcel-based summaries

After convergence of the MCMC algorithm in each parcel P, the samples of the quantities of interest are averaged over iterations to compute approximations of marginal posterior expectations:

x PM P = L 1 k=L 0 x (k) /L, L = L 1 -L 0 + 1, ∀ x ∈ {h, a, l, Θ} ⊂ P, ( 7
)
where L 0 stands for the length of the burn-in period and L the effective number of iterations. For classification purpose, we proceed in two steps:

1. Compute the PM estimates (p m j ) i of Pr(q m j = i | y j ) for i = -1, 0 using the following expression:

(p m j ) i = L 1 k=L 0 I (q m j ) (k) = i /L, ( 8
)
where I stands for the identity function. Then, deduce (p m j ) 1 from the constraint of unitary probability mass:

(p m j ) 1 = 1 -(p m j ) -1 -(p m j ) 0 .
2. Sort the probabilities (p m j ) i and select the MAP estimate:

( q m j ) MAP = arg max i Pr(q m j = i | y j ) ≈ arg max i (p m j ) i . (9)
whatever the number of components in the mixture. The MAP estimator is easily obtained in the two-class mixture case: V j is non-activating (( q m j ) MAP = 0) for the mth condition if (p m j ) 1 < 0.5.

In combination with these PM estimates, one can attach uncertainty measures to the NRLs. More precisely, the error bars are derived as follows:

e m j = 1 L L 2 k=L 1 (σ m i,j ) (k) with i = ( q m j ) MAP . ( 10
)
Interestingly, σ m 0,j is directly given by v m 0,j since the full conditional posterior distribution of the zero class is Gaussian, i.e., N (µ m 0,j , v m 0,j ). In contrast, the standard deviations (SD) σ m ±1,j require further computation since these full conditional densities are gamma-Gaussian (see Section A.5). As derived in Eq. (A.13), the variance of a gamma-Gaussian density admits a closed form expression, which gives σ m ±1,j after taking the square root. These SD estimates are then plugged into (10) to get corresponding error bars e m j .

The stochastic algorithm is summarised in Table 2.

Statistical comparisons for cognitive interpretation

Akin to the contrast definition in any GLM-based approach, statistical comparison between our task-related NRL estimates can be addressed in the proposed formalism. One might be interested in assessing unsigned or signed differences like using Fisher or Student-t tests, respectively in the classical hypothesis testing framework.

Let m and m be the indexes of the conditions we plan to contrast across the brain. This contrast can be assessed by measuring how close the voxelwise marginal distributions (p m j , p m j ) of the NRLs (a m j , a m j ) are in every voxel V j . Since these densities write as posterior mixtures, say p m j = i π i f m i,j , we start with identifying the MAP estimates ( q m j , q m j ) and then we compare the full conditional posterior densities (f m q m j ,j , f m q m j ,j ) instead of computing a distance between p m j and p m j . Hence, three different (respectively, six) situations may arise depending on the mixture prior in use (two or three-class mixture, respectively). The different cases correspond to all possible combinations of the pair ( q m j , q m j ):

Table 2 Gibbs sampling algorithm in a given parcel P. The parameters of the sampled distributions are derived in appendix B.

• Setting up: choose h 0 , a 0 , l 0 , θ 0 , θ 0 a .

• Iteration k: draw samples h k , a k , λ k , ( 2 ) k , θ k a from the conditional posterior pdfs:

--HRF:

h k ∼ N (µ h , Σ h ), --HRF variance: (σ 2 h ) k ∼ IG D/2, h t R -1 h/2 , --NRLs: for every condition m and every voxel j, --(u m j ) k ∼ U[0, 1] ; if (u m j ) k λ m -1,j , then (q m j ) k = -1 else if (u m j ) k λ m -1,j + λ m 0,j then (q m j ) k = 0, otherwise (q m j ) k = 1. --(a m j ) k | (q m j ) k = 0 ∼ N ((µ m 0,j ), v m 0,j ). (a m j ) k | (q m j ) k = ±1 ∼ GN (a m j | α ±1,m , µ m ±1,j , v m ±1,j ). --drift coefficients: ∀ j, ( j ) k ∼ N (µ j , Σ j ) --Noise variances: ∀ j, (σ 2 j ) k ∼ IG((N + 1)/2, ỹj 2 Λ j /2).
--AR parameters:

∀ j, (ρ j ) k ∼ 1 -ρ 2 j exp - A j 2σ 2 ε j ρ j - B j A j 2 1 (-1,1) (ρ j ).
--Mixture parameters: for every condition m, --Weighting probabilities λ m :

(λ m ) k ∼ D(δ ), with δ i = δ + Card C i,m = j ∈ 1 : J | κ m j = i =J i,m , ∀i = -1 : 1.
--Variance of NRLs for non-activating voxels:

(v 0,m ) k ∼ IG(η k 0,m , ν k 0,m ). --Shape parameters: (α ±1,m ) k ∼ exp(J i,m τ i,m α i,m )/Γ(α i,m ) J i,m I R + (α i,m ). --Scale parameters: (β ±i,m ) k ∼ G(J i,m α i,m + b i + 1, j∈C i,m a m j + c i ).
• Iterate until convergence is achieved. PMEs of {h, a, l, θ a } are computed using (7).

• Classification is performed according to the MAP criterion using ( 8)-(9).

a. if q m j = q m j = -1, voxel V j generates deactivations for both conditions. Comparing the NRLs (a m j , a m j ) is achievable by measuring how close (f m -1,j , f m -1,j ) are. This comparison therefore answers the question of deciding whether or not the deactivation is stronger for one condition w.r.t. the other (signed comparison) or if there is any difference between the two conditions (unsigned comparison). b. if q m j = q m j = 0, voxel V j is non-activating for both conditions. Comparing the NRLs (a m j , a m j ) amounts to computing a criterion between (f m 0,j , f m 0,j ). The interesting comparison consists in deciding whether or not there is some difference in the non-activating profile. c if q m j = q m j = 1, both conditions elicit activations in V j . By measuring how close (f m 1,j , f m 1,j ) are, we hope to know if activation occurring for condition m or m is stronger or if there is any difference irrespective of its sign. d. if q m j = -1 and q m j = 0, V j is deactivating in response to the mth stimulus type but is non-activating in response to the m th condition. To quantify this decision, one can measure a signed or unsigned criterion between (f m -1,j , f m 0,j ). By symmetry this case is equivalent to q m j = 0 and q m j = 1. e. if q m j = -1 and q m j = 1 or vice-versa, V j is activating in condition m and deactivating in condition m. To quantify this decision, one can measure a signed or unsigned distance between (f m -1,j , f m 1,j ). f. if q m j = 0 and q m j = 1 or vice-versa, V j is activating in condition m and non-activating in condition m. To quantify this decision, one can measure a signed or unsigned criterion between (f m 0,j , f m 1,j ).

Due to the use of mixture models, these comparisons can allow us to assess the null hypothesis (H 0 : a m j = a m j ) or the alternative one (e.g., H 1 : a m j = a m j or H 1 : a m j < a m j ) depending on the computed criterion. The question is now to define what kind of signed or unsigned criteria we can implement to quantitatively discriminate the two underlying distributions (f m i,j , f m i ,j ).

Unsigned task comparison

Unsigned comparison between f m i,j and f m i ,j can be computed using the Kullback-Leibler (KL) divergence i.e.,

D(f m i,j f m i ,j ) = R f m i,j (a) log f m i,j (a) f m i ,j (a)
da.

In the present case, its exact computation is only feasible when the two distributions are Gaussian i.e., when i = i = 0 (case b); see (A.3) in appendix A for details. Otherwise, an approximation of D(• •) has to be derived. For doing so, we proceed as follows. In cases (a, c, e), the sampling step of the NRLs (a m j , a m j ) relies on two Metropolis jumps, one for each NRL. The corresponding instrumental laws are truncated normal distributions (see (A.4) in appendix B). Therefore, we approximate f m i,j and f m i ,j by these positive Gaussian distributions which mean and variance parameters are given in (A.5)-(A.6). We end up by applying the KL divergence formula (Eq. (A.3)) to these truncated Gaussian approximations. In cases (d, f ), we proceed similarly for the single activating or deactivating component.

Signed task comparison

To go one step further and recover a sign information regarding the difference d m-m j = a m j -a m j , we need to estimate its posterior probability distribution f m-m j from a histogram H B j (.) with B time bins (β b ) b=1:B constructed over the last 500 iterations (i.e., the generated values (d m-m j

) (k) = (a m j ) (k) -(a m j ) (k)
) in any voxel of the mask M f . The posterior cumulative distribution function (cdf) F (•) can then be easily estimated from H B j (.). Contrast-based posterior probability maps (PPMs) are thus given by looking at differences d m-m j above a given threshold α:

P (d m-m j > α) = 1 -F (d m-m j α) = 1 - α -∞ f m-m j (t) dt (11) ≈ 1 - d n=1 H B j β n + β n+1 2 ∆β d < α d + 1. ( 12
)
where ∆β = β n+1 -β n . Setting α = 0, we actually find the voxels where a m j > a m j . Finally, we can threshold P (d m-m j > α) at level η to retain the voxels which make the comparison significant at this level (e.g., η = 0.95). Formally, the thresholded PPMs are given by P (d m-m j > α) > η. Note that this only provides uncorrected PPMs for multiple comparisons. The control of the familywise error is an open issue in the Bayesian formalism and is beyond the scope of this paper.

Results on synthetic data

Goal of the study

A comparison between two different prior mixture models has been done in (Makni et al., 2006a). In short, it has been shown that the gamma-Gaussian mixture model (GaGMM ) introduced on the NRLs is more efficient than a two-class Gaussian mixture model (GMM ) in terms of specificity: it provides a better control of the false positive rate. Similar conclusions have been drawn in [START_REF] Makni | Joint detectionestimation of brain activity in fMRI using an autoregressive noise model[END_REF] when considering an AR(1) noise model instead of a white Gaussian one in combination with a GMM prior. As the two changes induce higher computation time, it is worth assessing which modelling effort is preferable i.e. leads to the more significant improvement: the introduction of an inhomogeneous prior mixture or the consideration of serial correlation. For doing so, the models described in Table 1 are tested on the same artificial fMRI dataset.

Artificial fMRI dataset

These data were obtained by first generating two sets of trials, each of them corresponding to a specific stimulus (M = 2). These binary time series were then multiplied by a stimulus-dependent scale factor. Here, the functionally homogeneous region P consisted of J = 60 voxels. The number of activating voxels J 1,m was varied with the stimulus type m according to (J 1,1 , J 1,2 ) = (22, 30). Positive NRLs corresponding to activating voxels were simulated according to gamma pdfs: activating voxels :

a 1 j ∼ G(α 1 = 3, β 1 = 1), a 2 j ∼ G(α 2 = 10, β 2 = 2), non-activating voxels : a 1,2 j ∼ N (0, v 0,m = 0.1).
Remark that the chosen gamma parameter values yields a lower SNR for condition 1 ((µ 1 , v 1 ) = (3, 3) vs. (µ 2 , v 2 ) = (5, 2.5)). For all voxels, the binary stimulus sequence was convolved with the canonical HRF h c , whose exact shape appears in Fig. 2(a) in -line. An AR(1) noise b j was then added to the stimulus-induced signal m a m j X m h in every voxel V j . All AR parameters were set to the same value: (ρ j ) j=1:J =0.4, which is compatible with the serial correlation observed on actual fMRI time series. Also, a low SNR (SNR = 0.3) was considered in our simulations, in conformity with the real situation. Space-varying low-frequency drifts P j (generated from a cosine transform basis with coefficients j drawn from a normal distribution) were also added to the fMRI time courses according to (1).

General comments

As shown in Figs. Since the artificial fMRI time courses were synthetised using a GaGMM prior and some correlated noise, it is not surprising that the estimation performed under model M 1 provides the most accurate NRL estimates. Let us remark that the NRL estimates have a small but not negligible amount of bias, which is due to the bias/variance trade-off arising in the Bayesian approach in the non-asymptotic case. Nonetheless, we have checked that the bias tends to zero when the SNR increases. ), the noise model has only little influence on the recovered shape, as already advocated in [START_REF] Marrelec | Robust Bayesian estimation of the hemodynamic response function in event-related BOLD MRI using basic physiological information[END_REF]. As regards AR parameters, the estimated first order coefficients (ρ i ) i are close to the true values in every voxel for both models M 1 -M 3 (results not shown).

Influence of the noise model

We also assessed the sensitivity and the specificity of the four models. Figs. [2-5](d)-(e) show the posterior mean estimates (p m j ) 1 of deciding that voxel V j lies in class 1, i.e., is activating for models M 1 -M 4 and conditions 1 and 2, respectively. These results confirm our expectations: the modelling of the temporal correlation significantly improves both the sensitivity and the specificity. A higher/lower value of (p m j ) 1 is obtained with M 1 -M 3 when V j is truly activating/non-activating (compare Figs. . This means that models M 1 -M 3 provide lower false positive (FP) and false negative (FN) rates than models M 2 -M 4 , respectively. This effect is stronger in condition 1. This is in agreement with the idea that the precision of the noise model plays a more important role at a lower SNR.

Influence of the mixture prior

Not surprisingly, the estimated NRLs are recovered more accurately using the true prior mixture (M 1 -M 2 ): compare Figs. c) for a white noise model. This effect is much more important at low SNR, i.e., for condition 1. However, we have checked that when the true NRLs of the activating voxels follow a Gaussian distribution, the estimated shape and scale parameters of the gamma density in the GaGMM mixture provide close estimates of the mean and variance parameters of the uncentered Gaussian distribution (results not shown).

We are now interested in assessing the differences between M 2 and M 3 . The purpose of such a comparison is to decide whether or not a good mixture type provides more accurate and sensitive results than a precise noise modelling. Contrasting Figs. [3][4](b) allows us to note that M 2 outperforms M 3 in terms of accuracy of estimation for the first experimental condition. The NRLs attached to the non-activating voxels are over-estimated, leading to a much larger bias. In case of high SNR arising for the second condition, the comparison of Figs. [3-4](b) is less clear. The small NRLs are still over-estimated but the large ones are better estimated using M 2 in some cases (e.g., voxels 27, 58, 60). In terms of detection, Fig. 3(d) shows that a single false negative (voxel 3) is retrieved by model M 2 for condition 1, while five FNs are found by model M 2 , as shown in Fig. 4(d) (voxels 3, 6, 19, 26, 32). Hence, model M 2 achieves better results in terms of sensitivity and specificity. There- fore, we conclude that introducing an inhomogeneous prior mixture is more powerful than modelling the serial correlation as regards both estimation and detection.

Receiver-operator-characteristic (ROC) curves have been also computed to quantitatively evaluate the differences between models M 1 -M 4 . Fig. 6 illustrates and confirms that model M 1 provides the most sensitive detection when specificity is fixed and a better specificity at a given sensitivity. These ROC curves also validate that model M 4 is the less sensitive and the less specific out of the four models. Finally, model M 2 outperforms M 3 and provides better results in terms of sensitivity and specificity, irrespective of the stimulus type. pact in detection at low SNR since the distance between continuous and dotted lines is larger in Fig. 6(a) than in Fig. 6(b), except at very low specificity (0.1). This holds whatever the mixture type.

Deactivation modelling

Our purpose was to compare an inhomogeneous two-class mixture model with its three class extension. In the latter case, a third class is used to account for putative deactivation phenomenon arising for instance during sustained bursts of interictal epileptiform activity [START_REF] Bagshaw | Analysis of the EEG-fMRI response to prolonged bursts of interictal epileptiform activity[END_REF][START_REF] Bénar | EEG-fMRI of epileptic spikes: Concordance with EEG source localization and intracranial EEG[END_REF]. Suitable artificial fMRI datasets were simulated accordingly. We considered a ROI consisting of J = 60 voxels. Let J -1,m , J 0,m J 1,m be respectively the number of deactivating, non-activating and activating voxels in response to condition m. We set (J 1,1 , J -1,1 ) = (28, 19) and (J 1,2 , J -1,2 ) = (32, 12), so that (J 0,1 , J 0,2 ) = (13, 16). We simulated the NRLs as follows:

activating voxels : a 1 j ∼ G(3, 1), a 2 j ∼ G(5, 2) non-activating voxels : a 1,2 j ∼ N (0, 0.1) deactivating voxels : -a 1 j ∼ G(5, 4), -a 2 j ∼ G(5, 4)

The same procedure as before (see §4.1.1) was applied to simulate artificial fMRI time series. The only difference concerns the noise type, which is white, Gaussian and homogeneous in space to save computation time (∀j, σ 2 j = 0.3). Hence, model M 2 and M 5 (see Table 1) were tested and compared in terms of estimation, detection performance and evidence.

The HRF estimates corresponding to models M 2 -M 5 are shown in Figs. [7][8](a), respectively. These estimated time courses appear very close to the true HRF shape. Figs. [7-8](b)-(c) show the NRL estimates related to conditions 1 and 2, computed for model M 2 and M 5 , respectively. First, we observe that M 2 provides under-estimated NRLs for activating voxels but over-estimated parameters for deactivating ones, irrespective of the stimulus type. The estimated error bars also appear significantly larger when deriving from M 2 . In contrast, model M 5 provides more reliable NRL estimates with smaller error bars, as illustrated in Fig. 8-(b)(c). Also, the mean square error is decreased for the NRLs corresponding to deactivating and non-activating voxels. Fig. 7(d)-(e) demonstrates that model M 2 reports a few FN voxels (see upward arrows). All these voxels have small NRL coefficients, inducing their assignment to class 0. More importantly, we observe that the truly non-activating and deactivating voxels are mixed in class 0, irrespective of the condition. Fig. 8(d)-(e) reports the posterior mean estimates (p m j ) i (see ( 8)), which are then combined to get the final classification according to the MAP criterion ( q m j ) MAP (see ( 9)). As indicated on these graphs, model M 5 produces an accurate classification. Fig. 8(d)-(e) respectively show the presence of three FN voxels for condition 1 and only two FNs for condition 2. These classification errors could be explained by the low values taken by the true NRL coefficients in these voxels, making likely the assignment to class 0.

Finally, note that modelling the third class induces a higher computation time. In our simulations, inferring the parameters of models M 2 and M 5 takes about 6 and 11 minutes, respectively. If the ROI is large or if the experimental paradigm involves numerous conditions, it seems reasonable to start with a careful analysis of the paradigm to anticipate potential deactivations before inferring upon parameters of M 5 instead of M 2 .

Bayesian model comparison

More formally, from a statistical point of view we compare models M 1 -M 5 by computing sample-based approximations to the model evidence p(y | M m ). That allows us to derive Bayes factors BF mn as ratios of model evidence (see Appendix C for computational details). Bayes factor provides us with good statistical summary for model comparison. As reported in Table 3, there is a strong evidence in favour of Model M 1 . More interestingly, our conclusion drawn from the parameter estimates are also confirmed when comparing M 2 with M 3 using Bayes factor (line 2, Table 3). This also holds for the comparison between the two-class and the three-class mixtures, M 2 and M 5 respectively (line 5, Table 3).

Table 3 Values of the integrated log-likelihood log p(y | M m ) computed from a stabilized version of the harmonic mean identity [START_REF] Raftery | Estimating the integrated likelihood via posterior simulation using the harmonic mean identity[END_REF] 5 Results on real fMRI data

fMRI experiment

MRI settings

The experiment was performed on a 3-T whole-body system (Bruker, Germany) equipped with a quadrature birdcage radio frequency (RF) coil and a head-gradient coil insert designed for echo planar imaging (EPI). Functional images were obtained with a T2*-weighted GE-EPI sequence with an acquisition matrix at the 64 × 64 in-plane spatial resolution and 32 slices. A highresolution (1 × 1 × 1.2 mm 3 ) anatomical image was also acquired for each subject using a 3-D gradient-echo inversion-recovery sequence.

Experimental paradigm and contrast selection

The reader may refer to [START_REF] Dehaene-Lambertz | Functional segregation of cortical language areas by sentence repetition[END_REF] for details about this fMRI experiment. In short, the motivation of this study was to measure the reduction in the neural activity subserving a cognitive representation when this representation is accessed twice (the so-called "repetition suppression" effect), resulting in a detectable adaptation of the measurable signal in fMRI [START_REF] Grill-Spector | fmri-adaptation: a tool for studying the functional properties of human cortical neurons[END_REF][START_REF] Naccache | Unconscious semantic priming extends to novel unseen stimuli[END_REF]. The experiment consisted of a single session of N = 216 scans lasting TR = 2.4 seconds each. Sixty sentences presented in a slow event-related design (SOA=14.4s) were recorded. Each sentence (S 1 ) could be repeated two (S 2 ), three (S 3 ) or four (S 4 ) times in a row. The main goal of our subsequent analysis was twofold. First, our primary interest was to exhibit regions which activation to a given sentence either decrease with repetition or keep a constant magnitude across the repetitions. Second, we were interested in inferring the hierarchical temporal organisation from the parcel-based HRF estimates along the superior temporal sulcus (STS).

Since the most significant habituation effect occurs between the first and second sentence repetitions, we modelled the four conditions S 1 -S 4 but we only studied the contrast S 1 -S 2 .

Pre-processings

As explained in Subsection 2.1.1, the grey matter's mask was first computed (see Fig. 9(a)) and then dilated using a 4mm-radius sphere to account for the width of the cortical ribbon. We checked that for nine out of ten subjects the raw fMRI data were motionfree approximately. All T 1 -weighted MRI images were normalised onto the MNI template and functional images were transformed accordingly. fMRI volumes were also spatially smoothed using a Gaussian kernel with F W HM = 6 mm along each direction. A first level analysis was performed for each subject using SPM2. The GLM modelled the four presentations of a given sentence with two regressors (the canonical HRF and its time derivative). Then, the parcellation was computed from the parameter estimates of this analysis. We chose a relevant

F-contrast c = [1, 0, -1, 0, • • • ; 0, 1, 0, -1, • • • , .
..] to study the habituation effect between the first and second presentations of a given sentence (S 1 -S 2 ). Fig. 9(c) depicts an axial view of this parcellation for the same slice (z = -4 mm).

Our approach strongly relying on a functional homogeneity assumption, we started by comparing the results using increasing parcel numbers (from K = 100 to K = 500 parcels in M s = M a ∪ M f ). We checked that K = 200 is large enough to guarantee a higher and sufficient degree of homogeneity. Here, the smallest and the largest parcels contained 44 and 190 voxels, respectively. Within each parcel, the degree of functional homogeneity was measured by computing a correlation matrix over the parameter estimates of the GLM.

Note that this could also be done over the fMRI signals attached to each parcel.

Results

Our method was tested on the nine datasets. Here, we only report results for Subject 1. Although the habituation effect and brain dynamics (i.e., the HRF shape) are subject to inter-individual variability in terms of spatial localisation and activation delay, the conclusions drawn for Subject 1 remain quite valid for the others.

In what follows, the proposed joint detection-estimation algorithm was applied to each parcel of Subject 1's brain. Activating voxels appear in red colour.

Probing for putative deactivation

This first analysis was devoted to looking at putative deactivations, that is the presence of negative NRLs. We actually performed tests on all parcels to assess differences between the GaGMM prior and its 3-class extension. For illustrative purpose, Fig. 11 depicts the results of such a comparison on two parcels P 1 and P 2 , composed of 129 and 135 voxels, respectively. Interestingly, the vast majority of voxels in P 1 elicit a coherent activation in response to the first presentation of a sentence (S 1 ), while in P 2 , most voxels are nonactivating. As shown in Fig. 11(a)-(b), the same S 1 -based classification map is obtained in P 1 irrespective of the mixture model. The same conclusion holds with respect to S 2 (results not shown). In P 2 , Fig. 11(c)-(d) illustrates that a few voxels move from the non-activating state to the deactivating one. However, the corresponding NRL estimates are of small magnitudes indicating that this new classification may arise by chance. Bayesian model comparison statistically confirms our result since numerical evaluation of Bayes factors gives us log BF 52 = -1.2 for P 1 and log BF 52 = -.5 for P 2 . These results show that there is less evidence in this dataset for supporting model M 5 . Therefore, the introduction of the third class in the mixture is not necessary to analyse these data, particularly in the brain regions involved in the treatment of phonological stimuli (language comprehension). In the rest of the paper, we restrict ourselves to the GaGMM prior.

Unsigned comparison between conditions

The result of our KL-distance map is shown in Fig. 12(a) in one slice of the brain (z = -4 mm). In Fig. 12(b)-(c), the KL-distance map is then compared to the standard F-test map by extracting the most significant, say MS, voxels from both volumes, according to the corresponding criteria (the largest KL divergence and the highest F value). This thresholding procedure shows common features in the activation patterns but also discrepancies in the temporal lobes that elicit different responses. On this slice, the KL-based criterion provides less activations. This may be due either to our approximation of the KL divergence or to the criterion itself.

Signed comparison between conditions

The mapping of the habituation phenomenon calls for signed comparisons since we are looking for voxels where a significant decrease of NRLs between S 1 and S 2 can be observed. Fig. 13 shows such a comparison over the whole brain between the PPM derived using the proposed methodology and the thresholded T-map obtained using SPM2. The activated regions (in blue) for the contrast S 1 > S 2 elicit therefore a higher response when a sentence is presented only once.

Statistical differences appear between the proposed PPM and the corresponding SPM. On the one hand, in the majority of slices we observe more activations on the PPM shown in Fig. 13(a). To a certain extent, these differences can be explained by the shape variations of the HRF estimate and its deviation from the canonical shape prescribed in SPM. Moreover, the correction for multiple comparisons used for the t-map may dramatically reduce the number of activating voxels. In contrast, no correction has been applied over the PPM.

On the other hand, in the temporal lobes (top row, right slices of Fig. 13(b)), sensitivity of detection seems better on the Student-t map: the activated clusters appear larger. Note also the presence of isolated activating voxels in Fig. 13(a). It is likely that this reflects the presence of false positives. To circumvent this issue and decrease the FP rate in the PPM, a spatial correlation model between neighbouring voxels can be introduced; see (Vincent et al., 2007b,a) and the discussion.

In Fig. 14, we represent the estimated HRFs in parcels corresponding to areas where we notice sensitivity differences. It clearly appears that in some parcels our HRF estimates exhibit unexpected timing properties. For instance, in central regions (x = 0, y < 0), we obtained initial dips that were difficult to predict in advance. This requires further analysis (see Discussion). In regions where the HRF estimate is very close to the canonical shape, the PPM and SPM provide similar activation patterns. Finally, some regions were also detected as activating (S 1 > S 2 ) by both methods while there is no evidence in the literature to suppose a priori that they elicit responses to auditory stimuli (see for instance Fig. 14(b) along the interhemispheric axis).

Habituation and temporal organisation

We focus on superior temporal regions ranging from the primary auditory cortex (Heschl's gyrus) to associative areas (middle and posterior STS). Although it has been shown in [START_REF] Dehaene-Lambertz | Functional segregation of cortical language areas by sentence repetition[END_REF] that repetitions affect both amplitude and delay of responses, we only model habituation effect on the NRLs by considering the different sentence presentations as different conditions. This procedure is not optimal but remains quite simple.

The first interesting region is Heschl's gyrus located in the primary auditory cortex around voxel with coordinates (-48, -12, 0) mm in the standard Talairach space. This area shows the same response magnitude each time a sentence was presented. Fig. 15(a)-(b) shows the NRL estimates in parcel P He , which is circled in black for S 1 and S 2 , respectively. In every voxel of P He , these magnitude parameters are very close to each other making the KL-distance between the marginal posterior distributions of S 1 and S 2 close to zero (see Fig. 15(c)). Hence, the measured difference between S 1 and S 2 is not statisti- cally significant. Summary statistics computed over the NRL estimates in P He are reported in Table 4 (left col.) and confirm our first analysis quantitatively.

The same study was done in P Br containing Broca's area and centered around voxel (-40, 24, 0) mm in the Talairach space. Fig. 15(d)-(e) clearly indicates a strong decrease in the NRLs between S 1 and S 2 . The higher value of the KL divergence reported in Fig. 15(f) confirms the presence of a significant habituation effect in P Br . Our quantitative analysis (see Table 4, right col.) also shows the same trend as outlined by the strong discrepancies between the mean and standard deviations of S 1 and S 2 .

In [START_REF] Dehaene-Lambertz | Functional segregation of cortical language areas by sentence repetition[END_REF], a sine-wave GLM was designed and fitted to study the speed of habituation. It allowed one to exhibit a temporal organisation of the temporal lobe with fastest responses located in Heschl's gyrus and slowest ones in temporal poles. Here, given the proposed methodology, the temporal organisation is studied more directly by measuring and sorting the timing properties of the parcels P He and P Br .

-47 In these parcels, we also investigated the temporal aspects of the STS organisation by measuring different features (time-to-peak T peak , time-toundershoot T undershoot ) on the estimated HRFs (see Fig. 15(g)-(h)). We computed these quantities over each parcel (P He and P Br and their four closest neighbours) separately before averaging them to get a mean estimate. We found that the responses in P He (T peak = 6 s. and T undershoot = 13.3 s.) occur and return to the baseline earlier than the responses in P Br (T peak = 6.4 s. and T undershoot = 13.6 s.). After averaging these quantities over the five parcels, we obtained congruent results (T peak = 6.1 s. and T undershoot = 13.2 s. around P He vs. T peak = 6.5 s. and T undershoot = 13.5 s. around P Br ) meaning that the region embedding Hesch's gyrus elicits brain activations faster than the region including Broca's area. This confirms more directly what has been already derived in [START_REF] Dehaene-Lambertz | Functional segregation of cortical language areas by sentence repetition[END_REF] although no statistical test is provided to assess the significance of this result. The next question concerns of course the putative reasons of these earlier responses in Hesch's gyrus. There is a large evidence in these datasets for supporting faster neurodynamics as the main origin of these results instead of faster hemodynamics. This can be checked for instance using complementary analysis like a phase analysis conducted in [START_REF] Dehaene-Lambertz | Functional segregation of cortical language areas by sentence repetition[END_REF].

38 0 4 (a) (b) (c) -47 38 0 4 (d) (e) (f) % ∆ BOLD signal % ∆ BOLD signal (g) (h)

Discussion

In this paper, we have proposed an original method to perform a parcel-based joint detection-estimation of brain activity from fMRI data. It has been shown on simulated datasets that a gamma-Gaussian mixture as prior pdf on the NRLs outperforms a Gaussian mixture in terms of sensitivity/specificity tradeoff. It has also been reported that the noise model has an influence over this compromise, particularly at lower SNR: a first-order AR model provides lower false positive and negative rates in comparison with a white noise model.

Our method extends previous works [START_REF] Makni | Joint detectionestimation of brain activity in functional MRI: a multichannel deconvolution solution[END_REF](Makni et al., , 2006b,a) ,a) to deal with anatomically informed whole brain analysis. As already done in [START_REF] Smith | Assessing brain activity through spatial Bayesian variable selection[END_REF][START_REF] Nieto-Castanon | Region of interest based analysis of functional imaging data[END_REF][START_REF] Flandin | Improved detection sensitivity of functional MRI data using a brain parcellation technique[END_REF], analysis was constrained to the mask of the grey matter obtained from a segmentation of the T1-weighted MRI. Our approach also relies on functional homogeneity as-sumptions at a regional scale that can be assessed either from the fMRI time series themselves or from the GLM parameter estimates. To meet these conditions, we resort to an automatic parcellation technique developed in [START_REF] Thirion | Dealing with the shortcomings of spatial normalization: Multi-subject parcellation of fMRI datasets[END_REF] but alternative clustering strategies may be thought of. Our approach therefore depends on this prior decomposition making the global within-subject analysis a two-steps procedure. The quality of the parcels will have an impact on the model fitting and a slight modification of the parcellation may generate different results especially in case of identifiability problems. By varying the number of parcels, we have checked that our results remain quite stable for different parcellations. The solution to this problem actually lies in the coupling of the parcellation procedure with our detectionestimation approach. This a very appealing direction of research but remains beyond the scope of the present work. At the expense of an increased computational complexity, the two steps could be merged in a combined approach through a hierarchical Bayesian model: one might be interested in improving the parcellation from the results of the detection-estimation stage using an iterative strategy: neighbouring parcels would be grouped if their underlying hemodynamics share similar features. The algorithm should take place in the context of reversible jumps MCMC to properly handle fusion/segragation moves between parcels [START_REF] Green | Reversible jump MCMC computation and Bayesian model determination[END_REF]; [START_REF] Richardson | On Bayesian analysis of mixtures with an unknown number of components (with discussion)[END_REF]. Besides, the parcellation identification issue could also be attacked using triplet Markov fields [START_REF] Benboudjema | Unsupervised statistical segmentation of nonstationary images using triplet Markov fields[END_REF], which seem suitable for modelling nonstationarities in image segmentation.

A strong feature of our approach is the possibility to derive parcel-based HRF time courses throughout the brain. It allows us to assess the spatial variability of the HRF shape and to check that this shape greatly fluctuates across parcels.

Since the parcellation procedure is derived at the group level, one is able to compare subject specific HRFs in a given parcel. Doing so, we have noticed that the between-subject variability in the HRF shape seems to be larger than the within-subject spatial variability, as already suggested in [START_REF] Handwerker | Variation of BOLD hemodynamic responses across subjects and brain regions and their effects on statistical analyses[END_REF].

Our results also suggest that the modelling of spatial HRF fluctuations is important to segregate brain regions involved in the experimental paradigm (Heschl's gyrus, Broca's area, ...). The adaptation effect was particularly evident between the first and second sentences. The pattern of adaptation was different across regions with a set of regions demonstrating the same response each time a sentence was presented ( i.e., Heschl's gyrus) and regions showing a more or less strong decrease between the first and second presentation (eg, Broca's area, superior STS).

Nonetheless, there exist fMRI experiments for which the proposed approach may fail because of the inhomogeneity of the HRF shape both in space and across conditions at a regional scale. This may occur when neurodynamic and hemodynamic fluctuations intermix. Indeed, recent studies of the fine structure of the fusiform face area (FFA) [START_REF] Grill-Spector | High-resolution imaging reveals highly selective nonface clusters in the fusiform face area[END_REF] have shown that the FFA is actually highly heterogeneous. It appears that the FFA is composed in reality of several small scale subregions that respond strongly, not only to faces, but also to cars and sculptures. The subregions discovered in [START_REF] Grill-Spector | High-resolution imaging reveals highly selective nonface clusters in the fusiform face area[END_REF] were associated with very distinct HRFs. For such studies, the best way to capture these HRF fluctuations is to perform a voxel-based non-parametric analysis on the basis of which statistical comparisons can be done across conditions, as demonstrated in [START_REF] Ciuciu | Unsupervised robust non-parametric estimation of the hemodynamic response function for any fMRI experiment[END_REF][START_REF] Marrelec | Estimation of the hemodynamic response function in event-related functional MRI: Bayesian networks as a framework for efficient Bayesian modeling and inference[END_REF].

We have also distinguished some differences between our PPMs and the SPMs derived from a classical GLM-based analysis. Our results confirm the interest of a simultaneous procedure for detecting and estimating brain activity. The proposed procedure actually improves the sensitivity of detection in some regions where the temporal characteristics (time-to-peak, time-to-undershoot, ...) of our HRF estimate deviate from those of the canonical shape. Unfortunately, this effect is not systematic on the datasets since in other regions a loss of sensitivity was observed. The reasons underlying this unexpected decrease have to be identified. To elucidate this issue, future work will be devoted to the comparison of a degraded version of our joint detection-estimation procedure with the actual one. The degraded version corresponds to an inference scheme where the HRF is maintained fixed and is not sampled at every iteration of the Gibbs sampler. We will quantify statistically this sensitivity difference on simulated data where we know exactly the "ground" truth. For doing so, we could generate fMRI data over the whole brain using the fMRI simulator developed by [START_REF] Drobnjak | Development of a functional magnetic resonance imaging simulator for modeling realistic rigid-body motion artifacts[END_REF].

In the present paper, we get rid of the scale ambiguity problem due to the bilinearity of (1) w.r.t. the pair (h, a) by imposing a unitary norm constraint over h. This could have a dramatic impact on the convergence of our sampling scheme and then on the recovered HRF shape by distorting the target distribution of the MCMC scheme. Alternative strategies may first consist in cancelling this normalisation step. However, the resulting sampling scheme is too slow to converge in a reasonable amount of time [START_REF] Veit | Rééchantillonnage de l'échelle dans les algorithmes MCMC pour les problèmes inverses bilinéaires[END_REF]). An efficient alternative has been proposed in [START_REF] Veit | Rééchantillonnage de l'échelle dans les algorithmes MCMC pour les problèmes inverses bilinéaires[END_REF]. It has been applied to the joint detection-estimation of brain activity in [START_REF] Ciuciu | Application du rééchantillonnage stochastique de l'échelle en détection-estimation de l'activité cérébrale par IRMf[END_REF]. It consists in adding to the MCMC procedure a sampling step of a positive scalar parameter s coding for the HRF scale. It can be shown that its sampling is fast, follows a generalized inverse Gaussian distribution in case of Gaussian mixtures and guarantees the theoretical convergence of the generated Markov chain to the posterior distribution. Deriving the target distribution of this scale parameter for inhomogeneous mixtures is beyond the scope of this paper.

To conclude about the real impact of our within-subject analysis, inference should take place at the group level. In other words, we should compare the results of two random effect analyses (RFX) based on the same group statistics (e.g., mean effect) and statistical test (Student t-test). The first RFX analysis would correspond to the gold standard, in which the input data are given by the normalised effects of a standard individual SPM analysis. The second analysis would take the results of our algorithm for each subject as inputs. Of course, this is only feasible in case of multi-subjects parcellation, what is currently obtained using the procedure described in [START_REF] Thirion | Dealing with the shortcomings of spatial normalization: Multi-subject parcellation of fMRI datasets[END_REF][START_REF] Thirion | Structural analysis of fMRI data revisited: Improving the sensitivity and reliability of fMRI group studies[END_REF].

From a methodological point of view, we have shown that our joint detectionestimation technique is able to identify deactivations in the brain. This is owing to the introduction of a third class in the prior mixture model associated to the NRLs. Nonetheless, we did not exhibit real deactivations on the analysed datasets. In the future, we should therefore validate the 3-class extension on specific datasets. A good candidate could be a dataset acquired during an event-related auditory paradigm in which silence events are presented randomly to compare activations to a baseline derived from such events. As already shown in [START_REF] Ciuciu | Unsupervised robust non-parametric estimation of the hemodynamic response function for any fMRI experiment[END_REF], silence events may generate deactivations in the temporal lobe if they are presented to the subject when the gradients of the scanner are switched off. This will be the subject of further work.

Smoothing the data spatially provides a reliable manner for recovering clusters of activation instead of isolated spots, at the expense of a loss of resolution.

To avoid this preprocessing, the proposed method could be extended by introducing spatial correlation in the prior model. This could be done either on the NRLs (a) or on the underlying states (labels q). We argue in favour of the second solution for simplicity reasons. As already derived for Gaussian mixtures in (Vincent et al., 2007b,a), it is quite simple to sample from an Ising (2-class model) or Potts (3-class model) Markov random field (MRF) that enforce neighbouring voxels to be classified in the same state (e.g., activating). This approach actually seems more reasonable in terms of computational load than considering edge-preserving MRF based on non-quadratic potentials [START_REF] Green | Bayesian reconstructions from emission tomography data using a modified EM algorithm[END_REF][START_REF] Geman | Statistical methods for tomographic image reconstruction[END_REF]. Also, for computational reasons this extension has been developed in a supervised framework meaning that the hyper-parameter encoding spatial regularity of the hidden MRF is set by hand. Future work will be focused on an spatially adaptive extension in which this parameter is estimated as well, as already done in [START_REF] Woolrich | Mixture models with adaptive spatial regularization for segmentation with an application to fMRI data[END_REF][START_REF] Woolrich | Variational Bayes inference of spatial mixture models for segmentation[END_REF].

Another extension that could be introduced at little expense concerns the analysis on the cortical surface, as proposed in [START_REF] Andrade | Detection of fMRI activation using cortical surface mapping[END_REF]. This will probably improve the sensitivity of detection by constraining the analysis to the cortical surface. Such study needs first a segmentation of the anatomical MRI, then an extraction of the grey-white matter interface (e.g., as a mesh), and finally requires an interpolation of the fMRI signal on the nodes of the mesh (see for instance [START_REF] Grova | Anatomically informed interpolation of fMRI data on the cortical surface[END_REF] for a suitable approach).

Finally, the model presented here assumes that the NRLs are constant in time.

Hence, to account for putative habituation effects, it requires to model repetitions of the same stimulus as different experimental conditions, what may be not optimal in terms of sensitivity of detection. To account for trial-varying NRLs due to adaptation or learning effects arising either as a direct consequence of the paradigm or as a alteration of subject's arousal, the proposed model can be generalised in a way that makes the number of unknown parameters not too large. As proposed in [START_REF] Ciuciu | Modeling non-linear and nonstationary effects of the BOLD response using mixture models in fMRI[END_REF], habituation can be modelled at the voxel level by a pair of parameters: the NRL to the first trial of the stimulus and a mean habituation speed across the consecutive trials that follows a hyperbolic parametric model depending on the inter-stimulus intervals.

Hopefully, all these additional points will induce improvements in the detectionestimation results and will help to a better comprehension of brain functions.

A.1 Multivariate Normal density

The multivariate Normal density for d-dimensional variable x with mean vector µ and covariance matrix Σ is given by

N (x | µ, Σ) = |2πΣ| -1 2 exp - 1 2 (x -µ) T Σ -1 (x -µ) . (A.1)
The general formula of Kullback-Leibler (KL) divergence between a test density, say q(x) and a reference density p(x) is

KL(q||p) = q(x) log q(x) p(x) dx. (A.2)
For multivariate normal densities q(x) = N (x | µ q , Σ q ) and p(x) = N (x | µ p , Σ p ), Eq (A.2) takes the following form:

KL(q||p) = 1 2 log |Σ p | |Σ q | + 1 2 tr (Σ -1 p Σ q ) + 1 2 (µ q -µ p ) t Σ -1 p (µ q -µ p ) - d 2 .
For univariate normal densities q(x) = N (x | µ q , σ q ) and p(x) = N (x | µ p , σ 2 p ), the KL (A.2) distance becomes:

KL(q||p) = = σ 2 p -σ 2 q 2 + (µ p -µ q ) 2 σ 2 p + σ 2 q 4σ 2 p σ 2 p . (A.3)

A.2 Positive normal density

The truncated normal distribution for a scalar variable x ∈ R + with parameters (m, v) expresses as follows:

     N + (x | m, v) = C -1 exp -(x -m) 2 /2v I R + (x), C = √ πv 1 + erf m/ √ 2v / √ 2, (A.4)
where erf is the error function (Abramowitz and Stegun, 1970, p. 297): erf(z) = 2/ √ π z 0 e -t 2 dt. Parameter m defines the mode of the density if m > 0. Note that the knowledge of C is not required for simulating a realisa-tion of the density. Indeed, its mean µ and variance σ 2 are given by

µ = m + 2v π exp(-m 2 /2v) 1 + erf(m/ √ 2v) , (A.5) σ 2 = v + m 2 4 -µ - m 2 2 . (A.6)
Hence, the standard inversion technique of the cumulative distribution function [START_REF] Devroye | Non-Uniform Random Variate Generation[END_REF][START_REF] Gelfand | Bayesian analysis of constrained parameter and truncated problems using Gibbs sampling[END_REF]) may be used. First, it consists in simualting an uniform variate u ∼ U ([0, 1]) and then in calculating:

x = m + √ 2v -1 erf u + erf m/ √ 2v (u -1) ,
The erf function is approximated numerically in practice. In cases where the approximation error becomes important (i.e., when |m| is too large), this simulation method is inefficient. Instead, we use efficient alternatives which are based on accept-reject algorithms [START_REF] Robert | Simulation of truncated normal variables[END_REF], the most powerful relies on multiple instrumental distributions (Mazet et al., 2005) 9 .

A.3 Gamma density

The gamma density for variable x ∈ R + with shape parameter α > 0 and scale parameter β > 0, is defined by

G(x | α, β) = β α Γ(α) x α-1 exp (-βx) I R + (x) (A.7)
where Γ(x) is the gamma function defined as

Γ(x) = (x -1)Γ(x -1) = ∞ 0 t x-1 e -t dt. (A.8)
The mean and variance are respectively given by E [x] = α/β and var [x] = α/β 2 , while the mode of the distribution is given by M 0 = (α-1)/β. Particular cases of the gamma distribution are the Erlang distribution G(x | α, 1), the

exponential distribution E(x | β) = G(x | 1, β) and the chi-squared distribution G(x | ν/2, 1/2) denoted by χ 2 ν .
The Erlang distribution may be used in pratice as a preliminary step for simulating a gamma variate x ∼ G(α, β) since we get a right sample if x = u/β and u ∼ G(α, 1).

Note also that the inverse gamma distribution, denoted by IG(•α, β) throughout the paper is the distribution of x -1 when x ∼ G (•α, β).

where K α is the normalising constant:

K α = +∞ 0 x α-1 exp -(x -µ) 2 /2v d x = exp -µ 2 /4v v α 2 Γ(α) D α (-µ/ √ v). (A.11)
The last equation follows from (Gradshteyn and Ryzhik, 1994, p. 337, Eq 3.462-1) that relies upon the gamma function (see (A.8)) and the Parabolic Cylinder functions D ν (Gradshteyn and Ryzhik, 1994, p. 885, Eq 7.711-1 & p. 1065)10 . Importantly, the first two centered moments of a gamma-Gaussian random variable x ∼ GN (x | α, µ, v) can be computed analytically:

E x | α, µ, v = 1 K α +∞ 0 x α exp -(x -µ) 2 /2v d x = K α+1 K α (A.12) var x | α, µ, v = E x 2 | α, µ, v -E x | α, µ, v 2 = K α+2 K α - K 2 α+1 K 2 α . (A.13)
Simulating gamma-Gaussian random variables x ∼ GN (• | α, µ, v) is not straightforward and thus more effortful as compared to sampling from standard laws. Following [START_REF] Moussaoui | Separation of non-negative mixture of non-negative sources using a Bayesian approach and MCMC sampling[END_REF], to solve for this problem, we resort to a Metropolis-Hastings algorithm which needs the specification of an instrumental distribution q [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF][START_REF] Robert | The Bayesian Choice[END_REF]. To avoid high rejection rate, this instrumental pdf has to be chosen to fit the target distribution f = GN (• | α, µ, v) at best. In this regard, expression (A.10) is useful to characterisef in terms of mode, mean or variance from which the instrumental distribution q may be adjusted. Calculating the first-order derivative of (A.10) w.r.t. x and equating to zero, the mode of f is obtained as the solution of the following second order equation:

x 2 -µx -v (α -1) = 0, subject to x 0.

Let us denote ∆ = µ 2 + 4v (α -1).The mode of f , which is non-negative by definition, expresses as follows:

ν =      0 if ∆ < 0 max (µ + √ ∆)/2, 0 otherwise. (A.14)
Hence, the instrumental density q is taken as a truncated normal distribution N + (x | ν, v), which is easier to sample from as detailed in Appendix A.2.

A.6 Multinomial density

The density of a multinomial discrete distribution for variable x = {x 1 , . . . , x N } with parameters π = {π 1 , . . . , π N } is defined by

MN M (x | π) = M ! N i=1 x i ! N i=1 π x i i I i x i =M
where x i 0, π i > 0 and N i=1 π i = 1.

A.7 Dirichlet density

Let δ = {δ 1 , . . . , δ N } be some positive parameters. The probabilistic density function of the N -state Dirichlet distribution for variable π = {π 1 , . . . , π N } satisfying π i 0 with parameters δ, is defined by

D N (π | δ) = Γ( N i=1 δ i ) N i=1 Γ(δ i ) N i=1 π δ i -1 i I i π i =1 ,
where Γ(x) is, as before, the gamma function (A.8). Parameters δ i are prior observation counts for events governed by π i . The Dirichlet distribution is the conjugate prior of the parameters of a multinomial distribution. One special case is the symmetric Dirichlet distribution where δ i = δ 0 ∀i. In this case, the density becomes

D N (π | δ 0 ) = Γ(N δ 0 ) Γ(δ 0 ) N N i=1 π δ 0 -1 i .
The real vector (X 1 /S, . . . , X N /S) follows a Dirichlet distribution denoted as

D N (• | δ) if X i ∼ G(δ i , β)
are independent, and S = i X i . This holds true for any β, so in practice we choose β = 1. This result is very useful in practice for simulating realisations of a Dirichlet process.

B Computational details for the MCMC procedure

In this section, we derive the full conditional distributions of the quantities (h, (a, q), and Θ) to be sampled. When the sampling procedure w.r.t. a given parameter cannot be implemented as a Gibbs sampling step, we provide the reader with the derivations of some relevant instrumental distribution needed in the corresponding Metropolis Hastings move.

B.1 The HRF h and its scale σ

2 h Let us denote S j = m a m j X m . h is N (µ h , Σ h )-distributed with: Σ -1 h = σ -2 h R -1 + J j=1 S t j Λ j S j and µ h = Σ h J j=1 σ -2 ε j S t j Λ j y j -P j . (B.1) Variance σ 2 h is simulated according to p(σ 2 h | h) = IG D/2, h t R -1 h/2 .

B.2 The nuisance variables and their scale σ 2

Vectors j being independent (j = 1 : J), they can be sampled in parallel according to N (µ j , Σ j ) where

Σ -1 j = σ -2 I Q + σ -2 ε j P t Λ j P and µ j = σ -2 ε j Σ j P t Λ j (y j -S j h). (B.2)
Variance σ 2 is simulated according to IG((QJ + 1)/2, j || j || 2 /2).

B.3

The voxelwise mixtures (q, a)

Although we do not introduce any spatial correlation between the NRLs, the latter are sampled one-at-a time since the distribution p(a m j | rest)11 depends on a \m j = [a 1 j , , . . . , a m-1 j , a m+1 j , . . . , a M j ] due to the linearity of model ( 1) with respect to a j . The sampling of the NRLs is therefore implemented through two nested loops, the inner corresponding to the stimulus types (e.g., index m) and the outer to voxels (e.g., index j). Since the full conditional posterior p(a m j | rest) is a mixture, its sampling can be achieved in two steps. The first one consists in drawing a realisation of class q m j while the second one proceeds conditionally on class q m j . To carry out the first step, we need to identify the posterior mixture in voxel j and for condition m: p(a m j | rest) ∝ exp -

1 2σ 2 ε j ||e j,m -a m j g m || 2 Λ j 1 i=-1 λ i,m f i (a m j | θ i,m ) (B.3)
where g m = X m h, and e j,m = y j -P jn =m a n j g n .

After straightforward calculations including the normalisation of (B.3), we get the following expression: The mixing probabilities are given by

λ m i,j = 1 + i =i λm i ,j / λm i,j
-1 , ∀ i = -1 : 1 λm 0,j = λ 0,m v m 0,j /v 0,m 1/2 exp (µ m 0,j ) 2 /2v m 0,j ,

λm i,j = λ i,m β α i,m i,m
Γ(α i,m ) K m i,j exp (µ m i,j ) 2 /2v m i,j , for i ∈ {-1, 1} , and depend on the normalising constants K m i,j , for i = ±1, of the corresponding gamma-Gaussian densities; see (A.11) for its closed-form expression.

The parameters of the conditional posterior densities f i (a m j | θ m i,j ) are given by:

     v m 1,j = v m -1,j = σ 2 ε j (g t m Λ j g m ) -1 , v m 0,j = v -1 0,m + (v m 1,j ) -1 -1 µ m i,j = v m i,j σ -2 ε j g t m Λ j e m,j -i β i,m , ∀ i = -1 : 1. (B.8)
Sampling the class q m j first amounts to generating u m j ∼ U([0, 1]) and then to applying the following rules:

q m j              -1 if u m j λ m -1,j , 0 if λ m -1,j < u m j λ m -1,j + λ m 0,j , 1 otherwise.
Once q m j is correctly set, it remains to sample from the conditional distribution f i (• | θ m i,j ) as suggested by (B.4). If q m j = 0, this operation is straightforward because f 0 (• | µ m 0,j , v m 0,j ) is Gaussian (cf. (B.5)). However, if q m j = ±1, this operation is computationally more expensive since f ±1 (• | θ m ±1,j ) are gamma-Gaussian; see Appendix A.5 for details. Strictly speaking, the sampling of f -1 (• | θ m 1,j ) consists first in simulating a realisation of (-a m j ) using a wellsuited positive normal density and then negating that realisation.

Interestingly, when α i,m = 1, which corresponds to taking an exponential prior for the NRL distribution, the conditional posterior density f i (• | θ m i,j ) is exactly a truncated normal distribution with parameters equal to those of p i .

In that case, the Metropolis-Hastings is not necessary since all proposals are accepted (the acceptation rate equals to 1). The sampling of the truncated normal density can be achieved efficiently as detailed in Appendix A.

B.4 Mixture probabilities

Since we have Pr(q m j = i | λ m ) = λ i,m , for i = -1 : 1 and because the prior is conjugate, i.e., a symmetric Dirichlet distribution D(λ m | δ) with δ > 0, the full conditional posterior distribution is also Dirichlet and reads: p(λ m | q m j = i, δ) ∝ Pr(q m j = i | λ m ) p(λ m | δ) ∼ D(δ ), with δ i = δ + 1, and δ l = δ, ∀ l = i.

The spatial correlation being not modelled in (6), we may write Pr(q m = κ | λ m ) = j Pr(q m j = κ j | λ m ), with κ j = -1 : 1 and the joint posterior distribution p(λ m | κ, δ) is given by

p(λ m | κ, δ) ∝ j Pr(q m j = κ j | λ m )p(λ m | δ) ∼ D(δ ) (B.9) with δ i = δ + Card C i,m = j ∈ 1 : J | κ m j = i =J i,m
.

B.5 Mixture hyper-parameters

Variance v 0,m is very easy to sample because p(v 0,m | z m ) = IG (J 0,m -1)/2, ν 0,m /2 , where ν 0,m = j∈C 0,m (a m j ) 2 .

For the two other classes, we proceed to the sampling of the scale and shape gamma distribution parameters α ±1,m and β ±1,m , respectively. Following [START_REF] Moussaoui | Separation of non-negative mixture of non-negative sources using a Bayesian approach and MCMC sampling[END_REF], we use a Metropolis-Hastings step for α ±1,m parameters with a gamma instrumental density defined below. Simulating parameters β ±1,m is easier since they follow a gamma distribution. (a m j -s i )/J i,m .

This posterior distribution does not belong to a known family, so its simulation requires a MH jump. Akin to [START_REF] Moussaoui | Separation of non-negative mixture of non-negative sources using a Bayesian approach and MCMC sampling[END_REF], to obtain a good instrumental law q(α i,m ), we propose to approximate function g(α i,m ) using a gamma density G(t i,m , u i,m ). More precisely, parameters (u i,m , t i,m ) of this density are determined in order for its mode and inflexion points match those of function g(α i,m ). After some simple manipulations, we obtain:

t i,m = 1 + α 2 mode /(α mode -α infl ) 2 , u i,m = α mode /(α mode -α infl ) 2 , (B.11)
where α mode and α infl are the mode and the superior inflexion point (α infl > α mode ) of g(α i,m ). Calculating the first and second derivatives of g(α i,m ) yields these two non-linear equations that implicitly define α mode and α infl :

ψ(α mode ) =τ i,m and ψ (1) (α infl ) = (ψ(α infl ) -τ i,m ) 2 , (B.12)

where ψ is the digamma function defined by ψ(x) = d dx log Γ(x) and ψ (1) is its first derivative (trigamma function). Details about these functions are provided in (Abramowitz and Stegun, 1970, p. 253). The resolution of the two equations (B.12) is done using a root finding numerical method (cf. (Press et al., 1992, Ch. 9)). Finally, the posterior density (B.10) is simulated using a Metropolis-Hastings algorithm with a instrumental density q(α i,m ) chosen as a gamma distribution G(t i,m , u i,m ) whose parameters are given by t i,m = J i,m (t i,m -1) + 1, u i,m = J i,m u i,m .

(B.13)

The sampling of shape parameters β i,m is done according to the full conditional posterior distribution

p(β i,m | rest) ∝ β (J i,m α i,m +b i ) i,m exp -β i,m ( j∈C i,m a m j + c i ) ∼ G(J i,m α i,m + b i + 1,

B.7 AR parameters

For each voxel V j , we have: p j (ρ j ) = p(ρ j | rest) ∝ 1 -ρ 2 j exp -a j (ρ j -m j ) 2 1 (-1,1) (ρ j ), (B.15) where a j = A j /2σ 2 ε j and m j = B j /A j , with A j = N -1 n=2 ỹ2 j,n and B j =

N -1 n=1 ỹj,n ỹj,n+1 .

The density p j is log-concave, unfortunately it does not seem to belong to a referenced family of pdf, from which an efficient sampling technique would be available. Here, we propose to resort to a Metropolis-Hastings independence algorithm that uses a beta pdf g j ∼ Be(ζ j , κ j ) defined over (-1, 1)12 as the instrumental distribution:

g j (ρ) ∝ (1 + ρ) ζ j -1 (1 -ρ) κ j -1 , ∀ |ρ| < 1. (B.16)
The parameters ζ j and κ j have to be tuned in an appropriate way, so that g j approximates p j as closely as possible. Here, ζ j and κ j are chosen in such a way that log g j and log p j have the same curvature around a common maximizer over (-1, 1). Let us first remark that the maximizer r j of log p j is uniquely defined by 2a j (r j -m j )(1 -r 2 j ) + r j = 0, |r j | < 1. Moreover, r j takes an explicit expression, as the root of a polynomial of degree three. Then, (ζ j , κ j ) can be found by solving

    
(log g j ) (r j ) = 0 (log g j ) (r j ) = (log p j ) (r j ) which is a linear system. After some straightforward calculations, the solution can be expressed as follows:

ζ j = a j (1 + r j ) 2 (1 + m j -2r j ) + 3/2 κ j = a j (1 -r j ) 2 (1 -m j + 2r j ) + 3/2 It can be practically checked that g j (ρ) and p j (ρ) take very similar values for all ρ ∈ (-1, 1). Therefore, the proposal ρ j (sampled from g j ) has a high acceptation probability min 1, p j (ρ j )g j (ρ j )/p j (ρ j )g j (ρ j ) . In practice, the worse acceptation ratio that we observed was about 0.92.

  [2-5](a), all HRF estimates obtained using the four different models match the canonical time course h c pretty well. Figs. [2-5](b) show the corresponding NRL estimates that we obtained from models M 1 -M 4 , respectively in response to condition 1 while Figs. [2-5](c) summarize the same results for condition 2.

Figs

  Figs. [2-3](b)-(c) illustrate the impact of the noise model: a more accurate estimation of the NRLs, with smaller error bars and lower mean square error, is observed in Fig. 2(b)-(c) compared to Fig. 3(b)-(c), that is for model M 1 compared to model M 2 . This is a direct consequence of accounting for serial correlation in M 1 . The same conclusion holds when looking at Figs. [4-5](b)-(c), so irrespective of the prior mixture type. As regards the HRF estimate (compare Figs. [2-3](a)), the noise model has only little influence on the recovered shape, as already advocated in[START_REF] Marrelec | Robust Bayesian estimation of the hemodynamic response function in event-related BOLD MRI using basic physiological information[END_REF]. As regards AR parameters, the estimated first order coefficients (ρ i ) i are close to the true values in every voxel for both models M 1 -M 3 (results not shown).

  [2-3](d)-(e) for GaGMM priors or Figs.[ 4-5](d)-(e) for GMM priors)

  [2-4](b)-(c) one to another for an AR(1) noise model and observe the difference in Figs. [4-5](b)-(

Fig. 2 .

 2 Fig. 2. Estimation results on the simulated data using model M 1 . (a) HRF results: Symbols and • represent the true h c and its corresponding HRF estimate, respectively. (b)-(c): NRL estimates for conditions 1 and 2, respectively. True values appear on the x-axis and estimated values on the y-axis. The error bars follow Eq. (10). (d)-(e): PM estimates of activation probabilities pm j (• symbols) for the conditions 1 and 2, respectively. Symbols * depict the true class attached to each voxel.

Fig. 3 .

 3 Fig. 3. Simulation results using model M 2 . The same legend as in Fig. 2 holds. Only one FN voxel is present, indicated with an upward arrow.

Fig. 4 .

 4 Fig. 4. Simulation results using model M 3 . The same legend as in Fig. 2 holds. FN voxels are indicated by upward arrows.

Fig. 5 .

 5 Fig. 5. Simulation results using model M 4 . The same legend as in Fig. 2 holds. FP and FN voxels are indicated by downward and upward arrows, respectively.

Fig. 6 .

 6 Fig. 6. (a)-(b): ROC curves associated to the four different models for condition 1 (a) and condition 2 (b), respectively. Continuous line, interrupted line with •, continuous line with • and interrupted line with * represent the ROC curves for models M 1 -M 4 , respectively.

Fig. 7 .

 7 Fig. 7. Simulation results using model M 2 . FN voxels are indicated by upward arrows. The truly deactivating voxels that have been mixed in class 0 are surrounded by rectangle.

Fig. 8 .

 8 Fig. 8. Simulation results using model M 5 . The same legend as in Fig.2 holds. FP and FN voxels are surrounded by rectangles (3 FNs in (d) and 2 FPs in (e)). •, • and * symbols represent (p m j ) 1 , (p m j ) -1 and (p m j ) 0 , respectively.

  for models M m with m ∈ N * 5 . Model comparison based on the computation of Bayes factors log BF mn = log p(y | M m ) -log p(y | M n ) for every pair (m, n). NR stands for Not Relevant. Model m log p(y | M m ) Fig. # log BF mn , n = 1

Fig. 9 .

 9 Fig. 9. (a): Slice of Subject 1's anatomical mask (z = -4 mm). (b): its dilated version M a to match the functional resolution. (c): corresponding parcellation in the same slice. Each colour codes for a different parcel.

  Fig 10(a)-(b) shows the maps of the NRL estimates corresponding to conditions S 1 and S 2 , respectively, in a given slice of the brain. In the same slice, Fig 10(c)-(d) shows the activation probability maps attached to S 1 and S 2 (see (9) for details) our algorithm provides.

Fig

  Fig. 10. (a)-(b) NRL estimates in one slice of subject 1's brain (at z = -4 mm) in response to S 1 (a) and S 2 (b). Values correspond to the NRL coefficients only for voxels belonging to M s . Otherwise they are equal to 0. (c)-(d): detection results in the same slice for S 1 and S 2 , respectively. Voxels colour-coded in red are detected as activating. In black are the parcel borders that are superimposed to the different map results to show the parcellisation influence on the estimation of such parameters.

Fig. 12 .

 12 Fig. 11. Comparison of S 1 -based classification maps between the GaGMM prior (left column, (a)-(c)) and its 3-class extension (right column, (b)-(d)) in two different parcels, P 1 (top row, at y = -16 mm) and P 2 (bottom row, at y = -52 mm). Square spots in white, orange and brown match with activating, non-activating and deactivating voxels, respectively.

Fig. 13 .

 13 Fig. 13. From left to right and from top to bottom, brain slices organised along increasing axial axis (the bottom of the brain appear first in the superior left corner). Study of the contrast S 1 > S 2 . (a): thresholded PPM (see (11)) at η = 0.999; (b): thresholded SPM t-map at T=3.09 (corrected for multiple comparisons) obtained using SPM2.

Fig. 14 .

 14 Fig. 14. (a): Single-slice comparison (z = 0 mm, slice 7 in Fig.13) for the contrast S 1 > S 2 . left panels: thresholded PPM at η = 0.999; right panels: thresholded SPM t-map at T=3.09 (corrected for multiple comparisons) obtained from the SPM analysis. (b): same comparison in slice z = 24 mm (slice 13 in Fig.13). HRF estimates (in blue) in these regions compared to the canonical hemodynamic response function used in SPM (in red).

Fig. 15 .

 15 Fig. 15. (a)-(b) & (d)-(e): NRL estimates in the sagittal slices located at x = -48 mm and x = -40 mm, respectively. (a)-(b) provide the magnitudes in response to S 1 while (d)-(e) give us the NRLs in response to S 2 . (c)-(f ): KL-distance maps between the corresponding NRLs. Images are superimposed to functional data. The parcels containing the Heschl's gyrus (top row) and Broca's area (bottom row) are surrounded in black and referenced as P He and P Br in the text. (g)-(h): HRF estimates in P He and P Br , respectively.

  j ) = -GN (-a m j | α -1,m , -µ m -1,j , v m -1,j ) (B.7)

  Sampling the noise variances σ 2 ε can be performed in parallel. Drawing a noise variance is straightforward because p(σ 2ε j | rest) = IG((N + 1)/2, ỹj2Λ j /2).

Table 4

 4 Summary statistics in parcels containing Heschl's gyrus and Broca's area.

	Statistics	P He , Heschl's gyrus	P Br , Broca's area
		S 1 (m = 1) S 2 (m = 2) S 1 (m = 1) S 2 (m = 2)
	max j ( a m j ) PM	39.45	31.29	46.57	22.33
	min j ( a m j ) PM	0.35	-0.88	-4.21	-0.32
	Mean ( a m j ) PM	7.96	3.76	19.81	9.91
	Median ( a m j ) PM	5.09	1.44	19.95	10.49
	Std ( a m j ) PM	7.62	6.72	8.95	4.92

  The posterior density of each hyper-parameter α i,m takes the formp(α i,m | rest) ∝ i,m ) = exp(τ i,m α i,m )/Γ(α i,m ), τ i,m = ln β i,m +

	j∈C i,m	β i,m α i,m Γ(α i,m )	(a m j ) α i,m -1 p(α i,m | s j∈C i,m

i ) ∝ g(α i,m ) J i,m I R + (α i,m ), (B.10)

where g(α

http://www.fil.ion.ucl.ac.uk/spm/

http://brainvisa.info

We applied the SPM2 motion-correction algorithm.

except potentially when the corresponding class is empty: J i,m = 0 for i = ±1.

The burn-in period is the starting part of the Markov chain built by any MCMC algorithm which is used to ensure that the subsequent samples follow the equilibrium target distribution, i.e., the posterior law.

http://www.gnu.org/software/gsl

see http://www.its.caltech.edu/ ∼ astraw/seppo.html

http://brainvisa.info.

This code available on-line at http://www.iris.cran.uhp-nancy.fr/francais/si/ Personnes/Perso Mazet/rpnorm-fr.htm.

See also http://mahieddine.ichir.free.fr for implementation of Parabolic Cylinder functions.

rest stands for the "remaining variables".

If x ∈ (0, 1) and x ∼ Be(ζ, κ) then ρ = 2x -1 is said Be(ζ, κ)-distributed over (-1, 1).

Acknowledgements

The authors are grateful to Guillaume Flandin, Tim Behrens and Mark Woolrich for fruitful discussion and helpful comments. This work has been partly supported by the ANR french agency under grant number ANR-05-MMSA-0014-03 (OPTIMED project).

Appendix A Densities

We give the definitions of the densities used throughout this paper. We also provide numerical recipes for efficient simulations according to complex distributions when necessary.

A.4 Beta density

The Beta density for variable x ∈ [0, 1] with shape parameter α > 0 and scale parameter β > 0, is defined by

where B (α, β) is the Beta function:

The mean and variance are respectively given by

.

The mode of the Beta distribution evolves w.r.t. the domains of (α, β):

Importantly, for simulating Beta-distributed random variate x ∼ Be(• | α, β), we proceeds as follows. First, generate two exponential variables

x as e 1 /(e 1 + e 2 ).

A.5 gamma-Gaussian density

The non-standard gamma-Gaussian density for variable x ∈ R + with shape α > 0, mean µ > 0 and variance parameters v > 0, is defined by

C Bayesian model comparison

We have introduced Bayesian model comparison through the computation of Bayes factors:

These quantities are thus computed as the ratios of integrated likelihoods or model evidence of the different models. Of course, they are approximated from the MCMC outputs using the metholodogy proposed in [START_REF] Kass | Bayes factors[END_REF][START_REF] Chib | Marginal likelihood from the Gibbs output[END_REF][START_REF] Chib | Marginal likelihood from the Metropolis-Hastings output[END_REF] and further developed in [START_REF] Raftery | Estimating the integrated likelihood via posterior simulation using the harmonic mean identity[END_REF]. The latter relies on the harmonic mean identity:

This suggests that the model evidence can be approximated by the harmonic mean of the likelihoods p(y | θ (l) , M * ) based on L draws θ (1) , . . . , θ (L) from the posterior distribution p(θ | y, M * ):

with L = L 1 -L 0 + 1. These sample might come out of a standard MCMC implementation. Although p(y | M * ) is consistent as the sample size L increases, its precision is not guaranteed: it may have an infinite variance. Therefore, we have implemented a stabilized version of p(y | M * ) which is presented in detail in [START_REF] Raftery | Estimating the integrated likelihood via posterior simulation using the harmonic mean identity[END_REF]. In short, it consists in replacing p(y | θ (l) , M * ) by p(y | f (θ (l) ), M * ) in (C.1) such that f is a measurable function of θ and a dimension reduction transformation. Since the voxels are assumed independent in space, one may proceed separately for each voxel:p(y | f (θ (l) ), M * ) = j p(y j | f (θ (l) j ), M * ). More precisely, this means that f can be derived by integrating out analytically the NRLs a j and the noise variance σ 2 ε j for each voxel V j . This seems sufficient to ensure that var (p(y | f (θ), M * )) -1 | y < ∞. Hence, in practice we consider the following estimator:

Doing so, we have computed the log-evidence log p(y | M m ) of models M 1 -M 5 once they have been fitted against the first set of artificial data. Then, the logarithms of Bayes factors have been derived: log BF mn = log p(y | M m ) -log p(y | M n ), ∀(m, n) ∈ N * 4 .

The same procedure has been applied to models M 2 -M 5 with artificial data eliciting deactivations.