In standard within-subject analyses of event-related fMRI data, two steps are usually performed separately: detection of brain activity and estimation of the hemodynamic response. Because these two steps are inherently linked, we adopt the so-called region-based Joint Detection-Estimation (JDE) framework that addresses this joint issue using a multivariate inference for detection and estimation. JDE is built by making use of a regional bilinear generative model of the BOLD response and constraining the parameter estimation by physiological priors using temporal and spatial information in a Markovian model. In contrast to previous works that use Markov Chain Monte Carlo (MCMC) techniques to sample the resulting intractable posterior distribution, we recast the JDE into a missing data framework and derive a Variational Expectation-Maximization (VEM) algorithm for its inference. A variational approximation is used to approximate the Markovian model in the unsupervised spatially adaptive JDE inference, which allows automatic fine-tuning of spatial regularization parameters. It provides a new algorithm that exhibits interesting properties in terms of estimation error and computational cost compared to the previously used MCMC-based approach. Experiments on artificial and real data show that VEM-JDE is robust to model mis-specification and provides computational gain while maintaining good performance in terms of activation detection and hemodynamic shape recovery.

Functional Magnetic Resonance Imaging (fMRI) is a powerful tool to non-invasively study the relationship between a sensory or cognitive task and the ensuing evoked neural activity through the neurovascular coupling measured by the BOLD signal [1]. Since the 90’s, this neuroimaging modality has become widely used in brain mapping as well as in functional connectivity study in order to probe the specialization and integration processes in sensory, motor and cognitive brain regions [2–4]. In this work, we focus on the recovery of localization and dynamics of local evoked activity, thus on specialized cerebral processes. In this setting, the key issue is the modeling of the link between stimulation events and the induced BOLD effect throughout the brain. Physiological non-linear models [5–8] are the most specific approaches to properly describe this link but their computational cost and their identifiability issues limit their use to a restricted number of specific regions and to a few experimental conditions. In contrast, the common approach, being the focus of this paper, rather relies on linear systems which appear more robust and tractable [2, 9]. Here, the link between stimulation and BOLD effect is modelled through a convolutive system where each stimulus event induces a BOLD response, via the convolution of the binary stimulus sequence with the Hemodynamic Response Function (HRF ¹ ). There are two goals for such BOLD analysis: the detection of where cerebral activity occurs and the estimation of its dynamics through the HRF identification. Commonly, the estimation part is ignored and the HRF is fixed to a canonical shape which has been derived from human primary visual area BOLD response [10, 11]. The detection task is performed by a General Linear Model (GLM), where stimulus-induced components are assumed to be known and only their relative weighting are to be recovered in the form of effect maps [2]. However, spatial intra-subject and between-subject variability of the HRF has been highlighted [12–14], in addition to potential timing fluctuations induced by the paradigm (e.g. variations in delay [15]). To take this variability into account, more flexibility can be injected in the GLM framework by adding more regressors. In a parametric setting, this amounts to adding a function basis, such as canonical HRF derivatives, a set of gamma or logistic functions [15, 16]. In a non-parametric setting, all HRF coefficients are explicitly encoded as a Finite Impulse Response (FIR) [17]. The major drawback of these GLM extensions is the multiplicity of regressors for a given condition, so that the detection task becomes more difficult to perform and that statistical power is decreased. Moreover, the more coefficients to recover, the more ill-posed the problem becomes. The alternative approaches that aim at keeping a single regressor per condition and add also a temporal regularization constraint to fix the ill-posedness are the so-called regularized FIR methods [18–20]. Still, they do not overcome the low signal-to-noise ratio inherent to BOLD signals, and they lack robustness especially in non-activated regions. All the issues encountered in the previously mentioned approaches are linked to the sequential treatment of the detection and estimation tasks. Indeed, these two problems are strongly linked: on the one hand, a precise localization of brain activated areas strongly depends on a reliable HRF model; on the other hand, a robust estimation of the HRF is only possible in activated areas where enough relevant signal is measured [21]. This interdependence and retroactivity has motivated the idea to jointly perform these two tasks [22–24] (detection and estimation) in a Joint Detection-Estimation (JDE) framework [25] which is the basis of the approach developed in this paper. To improve the estimation robustness, a gain in HRF reproducibility is performed by spatially aggregating signals so that a constant HRF shape is locally considered across a small group of voxels, i.e. a region or a parcel. The procedure then implies a partitioning of the data into functionally homogeneous parcels, in the form of a cerebral parcellation [26].

As will be recalled in more detail in Section II, the JDE approach rests upon three main elements: i.) a non-parametric or FIR parcel-level modeling of the HRF shape; ii.) prior information about the temporal smoothness of the HRF to guarantee its physiologically plausible shape; and iii.) the modeling of spatial correlation between the response magnitudes of neighboring voxels within each parcel using condition-specific discrete hidden Markov fields. In [22, 23, 25], posterior inference is carried out in a Bayesian setting using a Markov Chain Monte Carlo (MCMC) method, which is computationally intensive and requires fine tuning of several parameters.

In this paper, we reformulate the complete JDE framework [25] as a missing data problem and propose a simplification of its estimation procedure. We resort to a variational approximation using a Variational Expectation Maximization (VEM) algorithm in order to derive estimates of the HRF and stimulus-related activity. Variational approximations have been widely and successfully employed in the context of fMRI data analysis: i.) to model auto-regressive noise in the context of a Bayesian GLM [27]; ii.) to characterize cerebral hierarchical dynamic models [28]; iii.) to model transient neuronal signals in a Bayesian dynamical system [29] or iv.) to perform inference of spatial mixture models for the segmentation of GLM effect maps [30]. As in our study, the primary goal of resorting to variational approximations is to alleviate the computational burden associated with stochastic MCMC approaches. Akin to [30], we aim at comparing the stochastic and variational-based inference schemes, but on the more complex matter of detecting activation and estimating the HRF whereas [30] treated only a detection (or segmentation) problem.

Compared to the JDE MCMC implementation, the proposed approach does not require priors on the model parameters for inference to be carried out. However, such priors may be injected in the adopted model for more robustness and to make the proposed approach fully auto-calibrated. Experiments on artificial and real data demonstrate the good performance of our VEM algorithm. Compared to the MCMC implementation, VEM is more computationally efficient, robust to mis-specification of the parameters, to deviations from the model, and adaptable to various experimental conditions. This increases considerably the potential impact of the JDE framework and makes its application to fMRI studies in cognitive and clinical neuroscience easier and more valuable. This new framework has also the advantage of providing straightforward criteria for model selection.

The rest of this paper is organized as follows. In Section II, we introduce the hierarchical Bayesian model for the JDE framework in the within-subject fMRI context. In Section III, the VEM algorithm based on variational approximations for inference is described. Evaluation on real and artificial fMRI datasets are reported in Section IV and the performance comparison between the MCMC and VEM implementations is carried out in Section V. Finally, Section VI discusses the pros and cons of the proposed approach and areas for further research.

Matrices and vectors are denoted with bold upper and lower case letters (e.g. P and μ). A vector is by convention a column vector. The transpose is denoted by ^t. Unless stated otherwise, subscripts j, m, i and n are respectively indexes over voxels, stimulus types, mixture components and time points. The Gaussian distribution with mean μ and covariance matrix Σ is denoted by (μ, Σ). The main acronyms used in the paper are defined in Table I. Table II gathers definitions of the main variables and parameters. A graphical representation of the model is given in Fig. 1.

A. The parcel-based model

We first recast the parcel-based JDE model proposed in [23, 25] in a missing data framework. Let us assume that the brain is decomposed in = ( ) _{γ =1: ϒ} parcels, each of them containing J_γ voxels and having homogeneous hemodynamic properties. The fMRI time series y _j is measured in voxel j ∈ at times (t_n ) _{n =1… N} , where t_n = nTR, N being the number of scans and TR, the time of repetition. The number of different stimulus types or experimental conditions is M. For a given parcel containing a group of connected voxels, a unique BOLD signal model is used in order to link the observed data Y = {y _j ∈ ℝ ^N , j ∈ } to the unknown HRF h _γ ∈ ℝ ^{D +1} specific to , and also to the unknown response amplitudes A = {a ^m , m = 1 … M} with being the magnitude at voxel j for condition m. More specifically, the observation model at each voxel j ∈ is expressed as follows [23]:

(1)

where S _j h _γ is the summation of the stimulus-induced components of the BOLD signal. The binary matrix is of size N × (D + 1) and provides information on the stimulus occurrences for the m-th experimental condition, Δt < TR being the sampling period of the unknown HRF h _γ = {h_{d Δ t} , d = 0 … D} in . This hemodynamic response is a consequence of the neuronal excitation which is commonly assumed to occur following stimulation. The scalars ’s are weights that model the response magnitude evoked by the stimuli, whose occurrences are informed by the X _m matrices (m = 0 … M). They model the transition between stimuli and the vascular response informed by the filter h _γ . It follows that the ’s are generally referred to as Neural Response Levels (NRL). The rest of the signal is made of matrix P, which corresponds to physiological artifacts accounted for via a low frequency orthonormal function basis of size N × O. With each voxel j is associated a vector of low frequency drifts ℓ _j ∈ ℝ ^O which has to be estimated. Within parcel , these vectors may be grouped into the same matrix L = {ℓ _j , j ∈ }. As regards observation noise, the b _j ’s are assumed to be independent with at voxel j (see Section II-B.1 for more details). The set of all unknown precision matrices (inverse of the covariance matrices) is denoted by Γ = {Γ _j , j ∈ }. The forward BOLD model expressed in Eq. (1) relies on the classical assumption of a linear and time invariant system which is adopted in the GLM framework [2]. Indeed, it can easily be recast in the same formulation where the response magnitudes ’s and drift coefficient ℓ _j ’s are equivalent to the effects associated with stimulus-induced and low frequency basis regressors, respectively. However, the JDE forward model generalizes the GLM model since the hemodynamics filter is unknown.

Finally, detection is handled through the introduction of activation class assignments Q = {q ^m , m = 1 … M} where and represents the activation class at voxel j for condition m. The NRL coefficients will therefore be expressed conditionally to these hidden variables. In other words, the NRL coefficients will depend on the activation status of the voxel j, which itself depends on the activation status of neighbouring voxels thanks to a Markov model used as a spatial prior on Q (cf Section II-B.2.c). Without loss of generality, we consider here two activation classes akin to [25] (activated and non-activated voxels). An additional deactivation class may be considered depending on the experiment as proposed within the JDE context in [31]. In the following developments, all provided formulas are general enough to cover this case.

B. A hierarchical Bayesian Model

In a Bayesian framework, we first need to define the likelihood and prior distributions for the model variables (Y, A, h _γ , Q) and parameters (Θ). Using the hierarchical structure between Y, A, h _γ , Q and Θ, the complete model is given by the joint distribution of both the observed and unobserved (or missing) data: p(Y, A, h _γ , Q; Θ) = p(Y | A, h _γ ; Θ) p(h _γ ; Θ)p(A | Q; Θ) p(Q; Θ). To fully define the hierarchical model, we now specify each term in turn.

1. Likelihood The definition of the likelihood depends on the noise model assumptions. In [23, 32], an autoregressive (AR) noise model has been adopted to account for serial correlations in fMRI time series. It has also been shown in [23] that a spatially-varying first-order AR noise model helps control the false positive rate. In the same context, we will assume such a noise model with where Λ _j is a tridiagonal symmetric matrix which depends on the AR(1) parameter ρ_j [23]: (Λ _j )_1,1 = (Λ _j ) _N,N = 1, for n = 2 : N − 1 and (Λ _j ) _{n +1 ,n} = (Λ _j ) _{n,n +1} = −ρ_j for n = 1 : N − 1. These parameters are assumed voxel-specific due to their tissue-dependence [33, 34]. Denoting and ȳ _j = y _j − Pℓ _j − S _j h _γ , the likelihood can be factorized over voxels as follows:

(2)

2. Model priors

a. Hemodynamic response function Akin to [23, 25], we introduce constraints in the HRF prior that favor smooth variations in h _γ by controlling its second order derivative: h _γ ~ (0, v _hR) with where D ₂ is the second-order finite difference matrix and v _h is a parameter to be estimated. Moreover, boundary constraints have also been fixed on h _γ as in [23, 25] so that h ₀ = h_{D Δ t} = 0. The prior assumption expressed on the HRF amounts to a smooth FIR model introduced in [18] and is flexible enough to recover any HRF shape.

b. Neural response levels Akin to [23, 25], the NRLs are assumed to be statistically independent across conditions: where θ_a = {θ _m , m = 1 … M } and θ _m gathers the parameters for the m-th condition. A mixture model is then adopted by using the assignment variables to segregate non-activated voxels from activated ones ( ). For the m-th condition, and conditionally to the assignment variables q ^m , the NRLs are assumed to be independent: . If then . It is worth noting that the Gaussianity assumed for a given NRL is similar to the assumption of Gaussian effects in the classical GLM context [10]. The Gaussian parameters θ _m = {μ_im , v_im , i = 0, 1} are unknown. For the sake of conciseness, we rewrite θ_a = (μ, v) where μ ={μ _m , m = 1 … M} with μ _m = {μ _{0 m} , μ _{1 m} } and v = {v _m , m = 1 … M} with v _m = {v _{0 m} , v _{1 m} }. More specifically, for non-activated voxels we set for all m, μ _{0 m} =0.

c. Activation classes As in [25], we assume prior independence between the M experimental conditions regarding the activation class assignments. It follows that where we assume in addition that p(q ^m ; β_m ) is a Markov random field prior, namely a Potts model. Such prior modeling assumption is consistent with the physiological properties of the fMRI signal where the activity is known to be correlated in space [33, 35]. Here, the prior Potts model with interaction parameter β_m [25] is expressed as:

(3)

and where Z(β_m ) is the normalizing constant and for all (a, b) ∈ ℝ², I(a = b) = 1 if a = b and 0 otherwise. The notation j ~ k means that the summation is over all neighboring voxels. Moreover, the neighboring system may cover a 3D scheme through the brain volume. The unknown parameters are denoted by β = {β_m , m = 1 … M}. In what follows, we will consider a 6-connexity 3D neighboring system.

For the complete model, the whole set of parameters is denoted by Θ = {Γ, L, θ_a, v _h, β} and belong to a set Θ.

We propose to use an Expectation-Maximization (EM) framework to deal with the missing data namely, A ∈ , h _γ ∈ , Q ∈ .

A. Variational Expectation-Maximization principle

Let be the set of all probability distributions on × × . EM can be viewed [36] as an alternating maximization procedure of a function on , for all q ∈ ,

(4)

where E _q [.] denotes the expectation with respect to q and (q) = −E _q [log q(A, h _γ , Q)] is the entropy of q. This function is called the free energy functional. It can be equivalently expressed in terms of the log-likelihood as (q, Θ) = log p(Y |Θ)–KL(q ||p(A, h _γ , Q | Y, Θ)) where KL(q || p(A, h _γ , Q | Y, Θ)) is the Kullback-Leibler (KL) divergence between q and p(A, h _γ , Q | Y, Θ):

(5)

Hence maximizing the free energy with respect to q amounts to minimizing the Kullback-Leibler divergence between q and the posterior distribution of interest p(A, h _γ , Q | Y, Θ). Since the KL divergence is always non-negative, and because the KL divergence of the posterior distribution to itself is zero, it follows easily that the maximum free energy over all q ∈ is the log-likelihood. The link to the EM algorithm follows straightforwardly. At iteration (r), denoting the current parameter values by Θ ^{( r −1)}, the alternating procedure proceeds as follows:

(6)

(7)

However, the optimization step in Eq. (6) leads to , which is intractable for our model. Hence, we resort to a variational EM (VEM) variant in which the intractable posterior is approximated by constraining the space of possible q distributions in order to make the maximization procedure tractable. In that case, the free energy optimal value reached is only a lower bound on the log-likelihood. The most common variational approximation consists of optimizing over the distributions in that factorize as a product of three pdfs on , and respectively.

Previous attempts to use variational inference [37, 38] and in particular in fMRI [27, 30] have been successful, with this type of approximations usually validated by assessing its fidelity to its MCMC counterpart. In Section IV, we will also provide such a comparison. The actual consequences of the factorization may vary with the models under study. Some couples of latent variables may capture more dependencies that would then need to be kept whereas others may induce only weak local correlation at the expense of a long-range correlation which to first order can be ignored (see [39] for more details on the consequences of the factorization for particular models). The way it may affect inference is that often variational approximations are shown to lead to underestimated variances and consequently to confidence intervals that are too narrow. Note that [38] suggested that nonparametric bootstrap intervals whenever possible may alleviate this issue. Also, when the concern is the computation of maximum a posteriori estimates, the required ingredients for designing an accurate variational approximation lie in the shape of the optimized free energy. All that is needed for our inference to work well is that the optimized free energy have a similar shape (mode and curvature) to the target log-likelihood whenever the likelihood is relatively large. As a matter of fact, there are cases, such as mixtures of distributions from the exponential family, where the variational estimator is asymptotically consistent. The experiments in [38] even report very accurate confidence intervals. Unfortunately no general theoretical results exist that would include our case to guarantee the accuracy of estimates based on the variational approximation. The other few cases for which the variational approach has shown good theoretical properties are to the best of our knowledge simpler than our current setting.

The fact that the HRF h _γ can be equivalently considered as a missing variable or a random parameter induces some similarity between our Variational EM variant and the Variational Bayesian EM algorithm presented in [37]. Our framework varies slightly from the case of conjugate exponential models described in [37] and more importantly, our presentation offers the possibility to deal with extra parameters Θ for which prior information may not be available. This is done in a maximum likelihood manner and avoids the use of non-informative priors that could be problematic (difficulties with non informative priors are listed in [40, pp. 64–65]). As a consequence, the variational Bayesian M-step of [37] is transferred into our E-step while our M-step has no equivalent in the formulation of [37].

B. Variational Joint Detection-Estimation

We propose here to use a EM variant in which the intractable E-step is instead solved over , a restricted class of probability distributions chosen as the set of distributions that factorize as p̃ _{A,H γ Q} = p̃ _A p̃ _{H γ} p̃_Q where p̃_A , p̃ _{H γ} and p̃_Q are probability distributions on , and , respectively. It follows then that our E-step becomes an approximate E-step, which can be further decomposed into three stages that consist of updating the three pdfs, p̃ _{H γ} , p̃_A and p̃_Q in turn using three equivalent expressions of when p factorizes as in . At iteration (r) with current estimates denoted by and Θ ^{( r −1)}, the updating rules become:

In other words, the factorization is used to maximize the free energy by alternately maximizing it with respect to p _{H γ} , p_A and p_Q while keeping the other distributions fixed. The steps above can then be equivalently written in terms of minimizations of some Kullback-Leibler divergences. The properties of the latter lead to the following solutions (see Appendix A for details):

(8)

(9)

(10)

The corresponding M-step is (since Θ and are independent, see Eq. (4)):

(11)

These steps lead to explicit calculations for and the parameter set . Although the approximation by a factorized distribution may have seemed initially drastic, the equations it leads to are coupled. More specifically, the approximation consists of replacing stochastic dependencies between latent variables by deterministic dependencies between moments of these variables (see Appendices A to C as mentioned below).

• VE-H step: From Eq. (8) standard algebra enables to derive that is a Gaussian distribution whose parameters are detailed in Appendix B. The expressions for and are similar to those derived in the MCMC case [23, Eq. (B.1)] with expressions involving the ’s replaced by their expectations with respect to (wrt) .
• VE-A step : Using Eq. (9), standard algebra rules allow to identify the Gaussian distribution of which writes as with . More detail about the update of is given in Appendix C. The relationship with the MCMC update of A is not straightforward. In [23, 25], the ’s are sampled independently and conditionally on the ’s. This is not the case in the VEM framework but some similarity appears if we set the probabilities either to 0 or 1 and consider only the diagonal part of .
• VE-Q step: Using the expressions of p(A|Q) and p(Q) in Section II, Eq. (10) yields which is intractable due to the Markov random field prior. To overcome this difficulty, a number of approximation techniques are available. To decrease the computational complexity of our VEM algorithm and to avoid introducing additional variables as done in [30], we use a mean-field like algorithm which consists of fixing the neighbours to their mean value. Following [41], can be approximated by a factorized density such that if ,

  (12)

  where q̃ ^m is a particular configuration of q ^m updated at each iteration according to a specific scheme, ~ j denotes neighbouring voxels to j on the brain volume and . Hereabove, and denote the m and (m, m′) entries of the mean vector ( ) and covariance matrix ( ), respectively. The Gaussian distribution with mean μ_im and variance v_im is denoted by ( . ; μ_im, v_im ), while . More details are given in Appendix D.
• M step: For this maximization step, we can first rewrite Eq. (11) as

  (13)

  The M-step can therefore be decoupled into four sub-steps involving separately (L, Γ), (μ, v), v _h and β. Some of these sub-steps admit closed-form expressions, while some other require resorting to iterative or alternate optimization. For more details about the related calculations, the interested reader can refer to Appendix E.

This section aims at validating the proposed variational approach. Synthetic and real contexts are considered respectively in sub-sections IV-A and IV-B. To corroborate the effectiveness of the proposed method, comparisons with its MCMC counterpart, as implemented in [25], will also be conducted throughout the present section. The two approaches have been tuned at best so as to make them as close as possible. This reduces essentially to moderate the effect of the MCMC priors. The hyper-priors have been parameterized by a set of hyper-parameters that have a limited impact on the priors themselves. For instance, for the mean and variance parameters involved in the mixture model (μ, v), conjugate Gaussian (0, w_μ ) and inverse-gamma (a_v , b_v ) hyper-prior distributions have been considered whose parameters have been tuned by hand so as to make them flat while proper (eg, w_μ < +∞). For doing so, the relationship between the hyper-parameters (eg, a_v , b_v ) and the statistical moments of the distribution have been carefully studied to guarantee a large variance or a large entropy in the latter. Moreover, given the above mentioned flatness of the hyper-priors, the MCMC approach is fairly robust to hyper-parameter setting.

A. Artificial fMRI datasets

In this section, experiments have been conducted on data simulated according to the observation model in Eq. (1) where P has been defined from a cosine transform basis as in [23]. The simulation process is illustrated in Fig. 2.

Different studies have then been conducted in order to validate the detection-estimation performance and robustness. For each of these studies, some simulation parameters have been changed such as the noise or the paradigm properties. Changing these parameters aims at providing for each simulation context a realistic BOLD signal while exploring various situations in terms of Signal-to-Noise Ratio (SNR).

1. Detection-Estimation performance The first artificial data analysis was conducted on data simulated with a Gaussian white noise (I _N is the N-dimensional identity matrix). Two experimental conditions have been considered (M = 2) while ensuring stimulus-varying Contrast-to-Noise Ratios (CNR) ² achieved by setting μ ₁₁ = 2.8, μ ₁₂ = 1.8, μ ₀₁ = μ ₀₂ = 0, and v ₁₁ = v ₁₂ = v ₀₁ = v ₀₂ = 0.5, so that a higher CNR is simulated for the first experimental condition (CNR _{m =1} = 15.68) compared to the second one (CNR _{m =2} = 6.48). For each of these conditions, the initial artificial paradigm comprised 30 stimulus events. The simulation process finally yielded time-series of 268 time-points. Condition-specific activated and non-activated voxels (Q values) were defined on a 20×20 2D slice as shown in Fig. 2[left]. No parcellation was performed. Simulated NRLs and BOLD signal were thus assumed to belong to a single parcel of size 20 × 20.

The Posterior Probability Maps (PPM) obtained using MCMC and VEM are shown in Fig. 3[middle] and Fig. 3[right]. PPMs here correspond to the activation class assignment probabilities . These figures clearly show the gain in robustness provided by the variational approximation. This gain consists of lower miss-classification error (a lower false positive rate) illustrated by higher PPM values, especially for the experimental condition with the lowest CNR (m = 2).

For a quantitative evaluation, the Receiver Operating Characteristic (ROC) curves corresponding to the estimated PPMs using both algorithms were computed. As shown in Fig. 4, they confirm that both algorithms perform well at high CNR (m = 1) and that the VEM scheme outperforms the MCMC implementation for the second experimental condition (m = 2).

Fig. 5 shows the NRL estimates obtained by the two methods. Although some differences are exhibited on the PPMs, both algorithms report similar qualitative results wrt the NRLs. However, the difference between NRL estimates (VEM-MCMC) in Fig. 5[right] points out that regions corresponding to activated areas for the two conditions present positive intensity values, which means that VEM helps retrieving higher NRL values for activated area compared to MCMC.

Quantitatively speaking, the gain in robustness is confirmed by reporting the Sum of Squared Error (SSE) ³ values on NRL estimates which are slightly lower using VEM compared to MCMC for the first experimental condition (m = 1: SSE_MCMC = 0.012 vs. SSE_VEM = 0.010), as well as for the second experimental condition (m = 2: SSE_MCMC = 0.010 vs. SSE_VEM = 0.009). These error values indicate that, even though the MCMC algorithm gives the most precise PPMs for the high CNR condition (Fig. 4, m = 1), the VEM approach is more robust than its MCMC alternative in terms of estimated NRLs. These values also indicate slightly lower SSE for the second experimental conditions (m = 2) compared to the first one (m = 1) with higher CNR. This difference is explained by the presence of larger non-activated areas for m = 2 where low NRL values are simulated, and for which the SSE is very low.

In addition, a test for equality of means has been conducted to test whether the quadratic errors means over voxels obtained with VEM and MCMC were significantly close. Very low p-values of 0.0377 and 0.0015 were obtained respectively for condition m = 1 and m = 2, which means that the null hypothesis (the two means are equal) is rejected for the usual 5% threshold. In other words, although the obtained error difference is small, it is significant from a statistical viewpoint. Interestingly, As m = 2 corresponds to a lower CNR, this highlights the significance of the gain in robustness we got with the VEM version in a degraded CNR context.

As regards HRF estimates, Fig. 6 shows both retrieved shapes using MCMC and VEM. Compared to the ground truth (solid line), the two approaches yield very similar results and preserve the most important features of the original HRF like the peak value (PV), time-to-peak (TTP) and time-to-undershoot (TTU).

2. Estimation robustness Since estimation errors may be caused by different perturbation sources, in the following, various Monte Carlo analyses have been conducted by varying one-at-a time several simulation parameters, namely the stimulus density, the noise parameters and the amount of spatial regularization. The differences in the obtained estimation errors have all been tested for statistical significance using tests for equality of means for the errors over 100 runs. For all reported comparisons, the obtained p-values were very low (p < 0.001) indicating that the performance differences, even when small, were statistically significant.

a. Varying the stimulus density In this experiment, simulations have been conducted by varying the stimulus density from 5 to 30 stimuli in the artificial paradigm, which leads to decreasing Inter Stimuli Intervals (ISI) (from 47 s to 9 s, respectively). At each stimulus density, 100 realizations of the same artificial dataset have been generated so as to evaluate the estimation bias and variance of NRLs. Here, the stimuli are interleaved between the two conditions so that the above mentioned ISIs correspond to the time interval between two events irrespective of the condition they belong to. A second order autoregressive noise (AR(2)) has also been used for the simulation providing a more realistic BOLD signal [23]. The rest of the simulation process is specified as before. In order to quantitatively evaluate the robustness of the proposed VEM approach to varying input SNR ⁴ , results (assuming white noise in the model used for estimation for both algorithms) are compared while varying the stimulation rate during the BOLD signal acquisition. Fig. 7 illustrates the error evolution related to the NRL estimates for both experimental conditions wrt the ISI (or equivalently the stimulus density) in the experimental paradigm. Estimation error is illustrated in terms of Mean Squared Error (MSE) ⁵ , which splits into the sum of the variance (Fig. 7[top]) and squared bias (Fig. 7[bottom]). This figure shows that at low SNR (or high ISI), VEM is more robust in terms of estimation variance to model mis-specification irrespective of the experimental condition. At high SNR or low ISI, the two methods perform similarly and remain quite robust. As regards estimation bias, Fig. 7[bottom] shows less monotonous curves for m = 2, which may be linked to the lower CNR of the second experimental condition compared to the first one. However, the two methods still perform well since squared bias values are very low, meaning that the two estimators are not highly biased. As reported in Section IV-A.1, error values on NRL estimates remain comparable for both experimental conditions and all ISI values, although PPM results present some imprecisions for the low CNR condition (m = 2).

As regards hemodynamic properties, Fig. 8[left] depicts errors on HRF estimates inferred by VEM and MCMC in terms of variance (Fig. 8[left-top]) and squared bias (Fig. 8[left-bottom]) wrt the ISI (or equivalently the stimulus density). The VEM approach outperforms the MCMC scheme over the whole range of ISI values, but the bias and variance remain very low for both methods. When evaluating the estimations of the key HRF features (PV, TTP and TTU), it turns out that the TTP and TTU estimates remain the same irrespective of the inference algorithm, which corroborates the robustness of the developed approach (results not shown). As regards PV estimates, Fig. 8[right] shows the error values wrt the ISIs. The VEM algorithm outperforms MCMC in terms of squared bias over all ISI values. However, the performance in terms of estimation variance remains similar with very low variance values over the explored ISIs range. For more complete comparisons, similar experiments have been conducted while changing the ground truth HRF properties (PV, TTP, TTU), and similar results have been obtained. In contrast to NRL estimations, variance and squared bias are here comparable. We can even note higher squared bias for PV estimates compared to the variance.

b. Varying the noise parameters In this experiment, several simulations have been conducted using an AR(2) noise with varying variance and correlation parameters in order to illustrate the robustness of the proposed VEM approach to noise parameter fluctuation. For each simulation, 100 realizations were generated so as to estimate the estimation bias and variance.

Performance of the two methods are evaluated in terms of MSE (splits into the sum of the variance and squared bias). Fig. 9 illustrates in a semi-logarithmic scale for NRLs the variance and squared bias of the two estimators plotted against the input SNR when varying the noise variance. This figure clearly shows that the bias introduced by both estimators is very low compared to the variance. Moreover, the bias introduced by VEM is very low compared to MCMC. As regards the variance, our results illustrate that VEM also slightly outperforms MCMC.

Fig. 10 depicts the variance and squared bias plotted in semi-logarithmic scale against input SNR when varying the noise autocorrelation. Overall, the same conclusions as for Fig. 9 hold. Moreover, as already observed in [42] at a fixed input SNR value, the impact of a large autocorrelation is stronger than that of a large noise variance irrespective of the inference scheme. Comparing Figs 9 and 10, this property is mainly visible at low input SNR (as usually observed on real BOLD signals).

Although a slight advantage is observed for the VEM approach in terms of estimation error and for both experimental conditions, the two methods perform generally well with a relatively low error level. The slightly better performance of VEM compared to MCMC is likely due to a better fit under model mis-specification.

c. Varying the spatial regularization parameter This section is dedicated to studying the robustness of the spatial regularization parameter estimation. For doing so, the synthetic activation maps of Fig. 2 are replaced with maps obtained as simulations of a 2-class Potts model with interaction parameter β that varies from 0.5 to 1.4. When positive, this parameter (β) favors spatial regularity across adjacent voxels, and hence smoother activation maps. Fig. 11 shows the estimated mean value β̂ and standard deviations for β = {β_m , m = 1 … M} over 100 simulations using both algorithms and for the two experimental conditions. Three main regions can be distinguished for both experimental conditions. The first one corresponds to β ≤ 0.8, which approximatively matches the phase transition critical value for the 2-class Potts model. For this region, Fig. 11 shows that the VEM estimate (green curve) appears to be closer to the Ground truth (black line) than the MCMC one (blue curve). Also, the proposed VEM approach yields more accurate estimation, especially for the first experimental condition (m = 1) having relatively high CNR. The second region corresponds to β ∈ (0.8, 1.1), where MCMC inference becomes more robust than the VEM one. The third region is identified by β ≥ 1.1, where both methods give less robust estimation than for the first two regions. Based on these regions, we conclude that the variational approximation (mean-field) improves the estimation performance up to a given critical value. It turns out that such an approximation is more valid for low β values, which usually correspond to observed β values on real fMRI data.

When comparing estimates for the two conditions, the curves in Fig. 11 show that both methods generally estimate more precise β values for the first experimental condition (m = 1) having higher input CNR. For both cases, and across the three regions identified hereabove, the error bars show that the VEM approach generally gives less scattered estimates (lower standard deviations) than the MCMC one, which confirms the gain in robustness induced by the variational approximation.

Note here that estimated β values in the experiment of Section IV-A.1 lie in the first region for the first experimental condition ( and ). For the second condition (m = 2), and because input SNR is relatively low, no clear conclusion can be made since MCMC and VEM give relatively different values ( and ) and no ground truth is available since the activation maps have been drawn by hand and not simulated according to the Markov model.

In this section, the computational performance of the two approaches are compared on both artificial and real fMRI datasets. Both algorithms were implemented in Python and fully optimized by resorting to the efficient array operations of the Numpy library 7 as well as C-extensions for the computationally intensive parts (e.g. NRL sampling in MCMC or the VE-Q step for VEM). Moreover, our implementation handled distributed computing resources as the JDE analysis consists of parcel-wise independent processings which can thus be performed in parallel. This code is available in the PyHRF package ⁸ .

For both the VEM and MCMC algorithms, the same stopping criterion was used. This criterion consists of simultaneously evaluating the online relative variation of each estimate. In other words, for instance for the estimated ĥ _γ , one has to check whether . By evaluating a similar criterion c_A for the NRLs estimates, the algorithm is finally stopped once c_H ≤ 10⁻⁵ and c_A ≤ 10⁻⁵. For the MCMC algorithm, this criterion is only computed after the burn-in period, when the samples are assumed to be drawn from the target distribution. The burn-in period has been fixed manually based on several a posteriori controls of simulated chains relative to different runs (here 1000 iterations). More sophisticated convergence monitoring techniques [46] should be used to stop the MCMC algorithm, but we chose the same criterion as for the VEM to carry out a more direct comparison.

Considering the artificial dataset presented in Section IV-A.1, Fig. 14 illustrates the evolution of c_H and c_A with respect to the computational time for both algorithms. Only about 18 seconds are enough to reach convergence for the VEM algorithm, while the MCMC alternative needs about 1 minute to converge on the same Intel Core 4 - 3.20 GHz - 4 Gb RAM architecture. The horizontal line in the blue curve relative to the MCMC algorithm corresponds to the burn-in period (1000 iterations). It can also be observed that the gain in computational efficiency the VEM scheme brings is independent of the unknown variables (e.g. NRLs or HRF) we look at.

To illustrate the impact of the problem dimensions on the computational cost of both methods, Fig. 15 shows the evolution of the computational time of one iteration when varying the number of voxels (left), the number of experimental conditions (center) and the number of scans (right). The three curves show that the computational time increases almost linearly (see the blue and red curves) for both algorithms, but with different slopes. Blue curves (VEM) have steeper slopes than red ones (MCMC) in the three plots showing that the computational time of one iteration increases faster with VEM than with MCMC wrt the problem dimensions.

As regards computational performance on the real fMRI data set presented in Section IV-B and comprising 600 parcels, the VEM also appeared faster as it took 1.5 hours to perform a whole brain analysis whereas the MCMC version took 12 hours. These analysis timings were obtained by a serial processing of all parcels for both approaches. When resorting to the distributed implementation, the analysis durations boiled down to 7 mins for VEM and 20 mins for MCMC (on a 128-core cluster). To go further, we illustrate the computational time difference (Δt = t _MCMC −t _VEM) between both algorithms in terms of parcel size which ranged from 50 to 580 voxels. As VEM vs. MCMC efficiency appears to be influenced by the level of activity within the parcel, we resorted to the same criterion as in Section IV-B to distinguish non-activated from activated parcels and tag the analysis durations accordingly in Fig. 16.

Fig. 16(a) clearly shows that the differential timing Δt between the two algorithms is higher for non-activated parcels (blue dots) and increases with the parcel size, which confirms the utility of the proposed VEM approach especially in low CNR/SNR circumstances. To further investigate the gain in terms of computational time induced by VEM, Fig. 16(b) illustrates the gain factor (G_f = t _MCMC /t _VEM) for activated and non-activated parcels. This figure shows that the VEM algorithm always outperforms the MCMC alternative since G_f ≥ 1 in all parcels (see horizontal line in Fig. 16[b]). Moreover, the gain factor G_f is clearly higher for non-activated parcels for which the input SNRs and CNRs are relatively low, and we found G_f ∈ [2.7, 80].

In this paper, we have proposed a new intra-subject method for parcel-based joint detection-estimation of brain activity from fMRI time series. The proposed method relies on a Variational EM algorithm as an alternative solution to intensive stochastic sampling used in previous work [23, 25]. Compared to the latter formulation, the proposed VEM approach does not require priors on the model parameters for inference to be carried out. However, to achieve gain in robustness and make the proposed approach completely auto-calibrated, the adopted model can be extended by injecting additional priors on some of its parameters as detailed in Appendix E for β and v _h estimation.

Illustrations on simulated and real datasets have been deeply conducted in order to assess the robustness of the proposed method compared to its MCMC counterpart in different experimental contexts. Simulations have shown that the proposed VEM algorithm retrieved more accurate activation labels and NRLs especially at low input CNR, while yielding similar performance for HRF estimation. VEM superiority in terms of NRLs and labels estimation has been confirmed by mean equality statistical tests over estimation errors obtained through 100 Monte Carlo runs. Conducted tests showed very low p-values, which means that differences in terms of obtained error means were statistically significant. Simulations have also shown that our approach was more robust to any decrease of stimulus density (or equivalently to any increase of ISI value). Similar conclusions have been drawn wrt noise level and autocorrelation structure. In addition, our VEM approach provided more robust estimation of the spatial regularization parameter and more compact activation maps that are likely to better account for functional homogeneity. These good properties of the VEM approach are obtained faster than using the MCMC implementation. Simulations have also been conducted to study the computational time variation wrt the problem dimensions, which may significantly vary from one experimental context to another.

As a general comment, this performance of VEM compared to MCMC may be counterintuitive. The sampled chain in MCMC is supposed to converge to the true target distribution after the burn-in period. However, since non-informative priors are used in the MCMC model such as the uniform prior for Potts parameters β, imprecisions in the target distribution may occur. In this case, VEM may outperform MCMC as observed in our simulations. Also, in some of the experiments, the model assumes a 2-class Potts model for the activation classes while we used more realistic synthetic images instead. The images we used (e.g. the house shape) are more regular and realistic than would be a typical realization of a Potts model. Then, our results show lower MSE (variance and squared bias) for VEM when the noise model is mis-specified and we suspect that VEM is less sensitive to model mis-specification. In a mis-specification context, the variational approach may be favored by the fact that the factorization assumption acts as an extra regularizing term that smoothes out the solutions in a more appropriate manner. There exists a number of other studies in which the variational approach is compared to its MCMC counterpart and provides surprisingly accurate results [30, 38, 47, 48]. It is true that results showing the superior performance of VEM over MCMC are more seldom. For instance, the results reported in [38] show that their variational approach is highly accurate in approximating the posterior distribution. These authors show like us, smaller MSE for the variational approach vs MCMC on simulations and point out a MCMC sensitivity to initialization conditions. On a computational efficiency point of view, most works on variational methods lead to EM-like algorithms in which one iteration consists of updating sufficient statistics (e.g. means and variances) that characterize the distribution approximating the target posterior. In our setting the approximating distribution is made of Gaussian parts fully specified by their mean and variance. In contrast, sampling-based methods like MCMC are not focused on sufficient statistics computation, but rather simulate realizations from the full posterior. Approximating a limited number of moments is less complex than approximating a full distribution, which may also be an ingredient that explains improved VEM efficiency and performance.

Regarding real data experiments, VEM and MCMC showed similar results with a higher specificity for the former. Compared to the classical GLM approach, the JDE methodology yields similar results in the visual areas where the canonical HRF is well recovered, whereas in the areas involved in the CS contrast, the estimation of more adapted HRFs that strongly differ from the canonical version enables higher sensitivity in the activation maps. These results further emphasize the interest of using VEM for saving a large amount of computational time. From a practical viewpoint, another advantage of the proposed algorithm lies in its simplicity to track convergence even if the latter is local: the VEM algorithm only requires a simple stopping criterion to achieve a local minimizer in contrast to the MCMC implementation [25]. It is also more flexible to account for more complex situations such as those involving higher AR noise order, habituation modeling [49] or considering three instead of two activation classes with an additional deactivation class.

To confirm the impact of the proposed inference, comparisons between the MCMC and VEM approaches should also take place at the group level. The most straightforward way would be to compare the results of random effect analyses (RFX) based on group-level Student t-test on averaged effects, the latter being computed either by a standard individual SPM analysis or by the two VEM and MCMC JDE approaches. In this direction, a preliminary study has been performed in [14] where group results based on JDE MCMC intra-subject analyses provide higher sensitivity than results based on GLM based intra-subject analyses. Ideally, the JDE framework could be extended to perform group-level analysis and yield group-level NRL maps as well as group-level HRFs. To break down the complexity, this extension could operate parcel-wisely by grouping subject-dependent data into group-level functionally homogeneous parcels. This procedure would result in a hierarchical mixed effect model and encode mean group-level values of NRLs and HRFs so that subject-specific NRL and HRF quantities would be modeled as fluctuations around theses means.

Such group-level validations would also shed the light on the impact of the used variational approximation in VEM. In fact, no preliminary spatial smoothing is used in the JDE approach in contrast to standard fMRI analyses where this smoothing helps retrieving clearer activation clusters. In this context, the used mean field approximation especially in the VE-Q step should help getting less noisy activation clusters compared to the MCMC approach. Eventually, akin to [23], the model used in our approach accounts for functional homogeneity at the parcel scale. These parcels are assumed to be an input of the proposed JDE procedure and can be a priori provided independently by any parcellation technique [26, 50]. In the present work, parcels have been extracted from functional features we obtained via a classical GLM processing assuming a canonical HRF for the entire brain. This assumption does not bias our HRF local model estimation since a large number of parcels is considered (600 parcels) with an average parcel size of 250 voxels.

However, on real dataset, results may depend on the reliability of the used parcellation technique. A sensitivity analysis has been performed in [51] on real data and for MCMC JDE version, that assesses the reliability of this parcellation against a computationally heavier approach which tends, by randomly sampling the seed positions of the parcels, to identify the parcellation that retrieved the most significant activation maps. Still, it would be of interest to investigate the effect of the parcellation choice in the VEM context, and more generally to extend the present framework to incorporate an automatic online parcellation strategy to better fit the fMRI data while accounting for the HRF variability across regions, subjects, populations and experimental contexts. The current variational framework has the advantage to be easily augmented with parcel identification at the subject-level as an additional layer in the hierarchical model. Automatically identifying parcels raises a model selection problem in the sense of getting sparse parcellation (reduced number of parcels) which guarantees spatial variability in hemodynamic territories while enabling the reproducibility of parcel identification across fMRI datasets. More generally, a model selection approach can be easily carried out within the VEM implementation as variational approximations of standard information criteria based on penalised log-evidence can be efficiently used [52].