Abstract
Background
We present a simple, datadriven method to extract haemodynamic response functions (HRF) from functional magnetic resonance imaging (fMRI) time series, based on the Fourierwavelet regularised deconvolution (ForWaRD) technique. HRF data are required for many fMRI applications, such as defining regionspecific HRFs, effciently representing a general HRF, or comparing subjectspecific HRFs.
Results
ForWaRD is applied to fMRI time signals, after removing lowfrequency trends by a waveletbased method, and the output of ForWaRD is a time series of volumes, containing the HRF in each voxel. Compared to more complex methods, this extraction algorithm requires few assumptions (separability of signal and noise in the frequency and wavelet domains and the general linear model) and it is fast (HRF extraction from a single fMRI data set takes about the same time as spatial resampling). The extraction method is tested on simulated eventrelated activation signals, contaminated with noise from a time series of real MRI images. An application for HRF data is demonstrated in a simple eventrelated experiment: data are extracted from a region with significant effects of interest in a first time series. A continuoustime HRF is obtained by fitting a nonlinear function to the discrete HRF coeffcients, and is then used to analyse a later time series.
Conclusion
With the parameters used in this paper, the extraction method presented here is very robust to changes in signal properties. Comparison of analyses with fitted HRFs and with a canonical HRF shows that a subjectspecific, regional HRF significantly improves detection power. Sensitivity and specificity increase not only in the region from which the HRFs are extracted, but also in other regions of interest.
1 Background
In functional magnetic resonance imaging (fMRI), local brain activation temporally changes blood oxygenation, which induces a blood oxygenation level dependent (BOLD) contrast in MR images [1]. Given a model of the BOLD response to a stimulus pattern, statistics can be used to quantify the match between the predicted and measured signals in a voxel, and significant activation is assessed via hypothesis testing. Statistical parametric mapping (SPM) estimates parameters of the noise distribution in every voxel to determine a threshold for the computed statistic [2]. It uses the general linear model (GLM): the response to a stimulus pattern is modelled as the output of a linear, time invariant (LTI) system [3]. The response to one type of stimulus is modelled by convolving its temporal distribution with the response function of that type of stimulus. The total response is the sum of responses to all individual stimulus types. The stimulus pattern is known from the experimental setup. However, the haemodynamic response function (HRF), i.e., the temporal BOLD impulse response, is unknown. Essentially, it is a smooth curve that starts to rise two seconds after the stimulus, peaks after six seconds, and returns to baseline within 30 seconds.
The resolution of discrete fMRI time samples yields a coarse description of the HRF, which is a problem for designs with short and/or randomised stimuli. In slicewise MRI acquisition, there is no optimal sampling for all slices, and slight differences in temporal onset between voxels in different slices must be accurately modelled to obtain good sensitivity. This requires a continuoustime HRF model. It is sometimes preferable to use a common HRF for many experiments and many regions. In many cases though, a common HRF cannot be assumed [4,5]. Two ways to solve this problem are (i) including a HRF model with more basis functions [5], and (ii) estimating the shape of the HRF in the statistical analysis [4,6]. A drawback of the first and simple approach is that it covers only limited variation in the HRF shape, while more advanced approaches suffer from high complexity and long computation times. For previous work on the variability of the BOLD response between brain areas and even across repeated measurements of the same subject, see Neumann et al. [7].
Extracting a HRF from fMRI data requires assumptions about its shape, and is computationally expensive. A simple method described in the literature is selective averaging with a long interstimulus interval (ISI), assuming nonoverlapping responses [8,9]. Trials with overlapping responses are also averaged, ignoring the fact that the overlaps introduce errors [3,10,11]. Other studies use a functional description of the HRF, whose parameters are determined by curve fitting (examples in [1214]), one of which [12] uses frequencydomain deconvolution. Another technique based on the GLM expands the HRF into a set of basis functions [15,16].
Ciuciu et al. [4] use a Bayesian method to extract the HRF. The method can simultaneously extract multiple HRFs from multiple experiments. Woolrich et al. [6] use a fully Bayesian approach to fMRI modelling, including the HRF. Both these approaches use certain prior assumptions for the shape of the HRF: in the first case, causality, smoothness, and starting and ending at baseline level, and Gaussian temporal autocorrelations; in the second case, a set of predefined prior HRF shapes, and many priors for modelling the noise distribution and autocorrelation.
The HRF extraction method described here is datadriven instead of modeldriven. It is based on Fourierwavelet regularised deconvolution, ForWaRD for short, developed originally for denoising and deblurring of images [17]. ForWaRD combines frequencydomain deconvolution for identifying overlapping signals, frequencydomain regularisation for suppressing noise, and waveletdomain regularisation for separating signal and noise. It is related to recent waveletbased deconvolution techniques [1820], with the advantage that the roles of signal (for fMRI: sparse, high frequency) and response (for fMRI: smooth, low frequency) can be interchanged: unlike other wavelet and vaguelettebased deconvolution methods, ForWaRD does the first step entirely in the frequency domain.
The method described in this paper uses an fMRI data set and the stimulus time pattern, and produces a time series of image volumes, which contains the HRF in each voxel. Compared to the simple extraction methods above, this method has the advantage of taking overlapping responses into account. On the other hand, it is much simpler than the Bayesian ones described above, in that it does not rely on shape assumptions/priors of the HRF: the extracted time points are determined only by the fMRI signal and the stimuli: the method is datadriven. This important property means that the extracted HRF is not biased by any a priori model. The only requirements for ForWaRD are (i) LTI responses and (ii) separability of signal and noise in the frequency/wavelet domain. Studies of the linearity of the BOLD response indicate that (i) is a reasonable assumption, and given the smoothness of the HRF and the availability of smooth wavelet basis functions, (ii) is satisfied as well. Another important advantage is the fact that waveletbased signal processing methods are much better in preserving the 1/ftype temporal autocorrelations that are found in fMRI time signals than timedomain methods [21]. Models that do not preserve this structure yield biased results, leading to an increase in false positives [21,22]. Finally, our method is fast: HRF extraction from a single fMRI data set takes about the same time as spatial resampling.
To avoid possible misunderstandings we would like to stress that our method only concerns the relation between stimulus and haemodynamic response, but does not make any statement about the underlying neural response. Although it is known from the work by Logothetis et al. [23] that the BOLD response in the visual cortex of the monkey brain roughly corresponds to a convolution of the neuronal local field potential (LFP) with some neural response function (NRF), the exact form of this coupling is still unknown and remains an area of active research. In our method, we describe the haemodynamic response as a convolution of the stimulus (not the LFP) with an impulse response function, i.e., the HRF. Given the stimulus pattern our method can recover the HRF, but we do not claim to be able to extract the NRF via deconvolution.
The output time series can be used to determine subjectspecific, groupspecific, and regionspecific HRFs, e.g., by averaging time signals in the corresponding regions. We demonstrate this in an fMRI experiment using two acquisitions, where the HRF extracted from the first fMRI time series is used to predict responses in a second experiment. A functional description for the HRF is found by fitting a combination of two gamma density functions with variable parameters for height, dilation, peak location and lag. We show that contrast and localisation improve by using the extracted HRF instead of the canonical HRF from the SPM2 program (a sum of two gamma density functions (GDFs) with fixed parameters).
The remainder of this paper is organised as follows. Section 2.2 reviews the ForWaRD method for regularised deconvolution. Sect. 2.2.4 describes how ForWaRD is used for HRF extraction, and it describes how the method was tested on simulated noisy time series and demonstrates a possible application using eventrelated fMRI data. Sect. 4 contains some general conclusions.
2 Methods
This section briefly reviews the standard General Linear Model (GLM) in fMRI, and describes the ForWaRD method [17] for extracting the HRF.
2.1 The General Linear Model
Statistical parametric mapping (SPM) is a common analysis method for fMRI. It (a) computes a statistic for every voxel using the GLM, (b) chooses a threshold based on the parameters of the noise and multiple testing correction; (c) thresholds the statistic map.
The GLM describes the response in an fMRI experiment as a weighted sum of explanatory signals. An explanatory signal is the response of one stimulus type. Let the matrix Y_{[T×N] }denote the fMRI data set, whose elements y_{ij }have time index i = 1, ..., T and position index j = 1, ..., N . In the GLM,
X_{[T×M] }is the design matrix, whose columns are the explanatory signals, and are multiplied by the weights in matrix β_{[M×N]}. Matrix e_{[T×N] }is the residual signal for each y_{ij}. A leastsquares estimate b for β is given by (X^{T}X)^{1}X^{T}Y . Assuming a Gaussian temporal distribution of the residuals (this follows from the GLM if the temporal noise in Y is Gaussian), standard hypothesis testing can be used to assess the significance of the elements of b.
2.2 Deconvolution, ForWaRD, and HRF extraction
In the GLM, a single response g to a pattern f of stimuli of the same type, is a convolution of f with the impulse response h, plus an additive term representing noise and other confounding effects.
2.2.1 Deconvolution
Discrete circular convolution, denoted by '*', corresponds to multiplication in the frequency domain:
capital letters denoting Fourier transforms of the corresponding lowercase signals. In the absence of noise, and given g and f, it is possible to compute an estimate of h through deconvolution. In the frequency domain, the Fourier transform of is obtained by pointwise division:
Deconvolution of a noisy signal is an illposed problem. A unique solution may not exist, be meaningless, or at best unstable: if noise is amplified at frequencies k where F(k) is close to zero, parts of the unmodelled signal e may appear in the extracted response (see Fig. 2i). An illposed problem can be regularised by adding extra information (or constraints) to the problem, so that its solution approximates the noisefree case [24]. Examples of such constraints are: minimising the norm between the solution and the data (minimumnorm), and making the solution more stable (smooth around optimum). ForWaRD uses regularisation methods in the frequency and wavelet domains (see Fig. 1) to overcome this problem.
Figure 1. ForWaRD HRF extraction scheme. Fourier shrinkage (determined by λ) is applied to partially attenuate noise amplified during the inversion step. Subsequent wavelet shrinkage (determined by κ) effectively attenuates the residual noise.
Figure 2. Stages of ForWaRD. HRF coeffcients after: frequency domain inversion (+); frequency domain shrinkage (*); wavelet domain Wiener shrinkage (×).
2.2.2 Fourier shrinkage
Frequencydomain shrinkage attenuates the noise after the pointwise division, by multiplying each frequency coeffcient (k) by a factor λ (k). Two popular methods are Wiener shrinkage and Tikhonov shrinkage [17]:
Here is the variance of the noise e(n) in (2). The remarks at the end of this subsection explain how to estimate and determine optimal regularisation parameters α and τ. An estimate of the spectrum H(k)2 of the unknown response function needed for Wiener shrinkage is obtained by the iterative algorithm of Hillery and Chin [[25], Sect. 4]. The estimated response function is the inverse Fourier transform of (see Fig. 2ii).
2.2.3 Wavelet shrinkage
Shrinkage in the wavelet domain is done using waveletdomain Wiener shrinkage (WDWS), which reduces the noise and preserves details in the signal [17,26]. The discrete wavelet transform describes a sampled signal c^{0 }of length N as a sum of localised basis functions. A discrete wavelet transform with J levels of decomposition recursively splits the signal into an approximation part c_{J }and detail signals d_{1}, d_{2}, ..., d_{J}, which are weighted sums of shifted and dilated versions of a scaling function φ and an associated wavelet ψ, respectively. The fast wavelet transform [27], which performs downsampling at each level, is not shiftinvariant. ForWaRD uses a shiftinvariant discrete wavelet transform (SIDWT), which uses a polyphase decomposition (subsampling for all shifts) [28]. Reconstruction is possible via the inverse transform, denoted by SIIDWT.
WDWS is performed via two wavelet transforms (for details, see [[17], section IV.C]). First, a wavelet transform of is computed using (φ_{1}, ψ_{1}). A pilot estimate is obtained via scaledependent thresholding of the detail coeffcients , resulting in thresholded detail coeffcients , n =1,...,N. Another wavelet transform of with wavelet basis (φ_{2}, ψ_{2}) yields detail coeffcients . These are shrunk by waveletWiener filtering coeffcients:
Here is the noise variance at level j. Finally, the ForWaRD estimate is the SIIDWT of the shrunk coeffcients using wavelet basis (φ_{2}, ψ_{2}), see Fig. 2iii.
Remarks:
• The variance σ_{e }is estimated based on the median absolute detail (MAD) coeffcient of a onelevel wavelet transform [[29], §5.4] of the the noise e(n).
• In the frequencydomain shrinkage step, the parameters τ (Tikhonov) and α (Wiener) can be tuned to minimise the mean squared error (MSE) between the estimate and the exact response function h. The exact MSE cannot be computed, but it can be approximated very accurately, see [[17], section VII].
• A complete Matlab implementation of the original ForWaRD algorithm, on which our method is based, is available ^{1}
2.2.4 Using ForWaRD for HRF extraction
The ForWaRDbased HRF extraction scheme (see Algorithm 1) works as follows: The BOLD response to a stimulus pattern f is relative to the baseline. Lowfrequency trends in the baseline may contaminate the extracted HRF. Trends are removed beforehand by a standard waveletbased technique [30]: transform each time signal (with length N) to the wavelet domain, using a fast wavelet transform (FWT) of log_{2}(N)  3 levels; remove the detail coeffcients, and subtract the lowscale signal from the time series. The influence of the detrending on the ForWaRD estimates is minimal, because the lowpass trends are in a much lower frequency band than the HRF. Processing takes place on a voxelbyvoxel basis, so images can be partitioned to reduce the computation load. The output is a time series of volumes which, inside activated brain areas, contain the HRF. Implementation has been performed in Matlab (The Mathworks, USA).
Algorithm 1 ForWaRDbased HRF extraction scheme.
1: for all Voxels do
2: load the discrete time series g and stimulus pattern f;
3: compute and subtract the time series mean;
4: remove lowfrequency trends;
5: apply ForWaRD to the meancorrected and detrended signal g and stimulus pattern f to obtain the estimate of the HRF coeffcients.
6: end for
2.3 Tests on simulated fMRI time signals
The routine presented in Sect. 2.2.4 was tested on signals with varying properties (SNR, temporal resolution, lowfrequency trends) with different settings of the routine itself (decomposition level, wavelet filters, etc.). Figures 3a–d shows the test setup: (a) convolve a stimulus pattern with a known HRF to obtain the activation signal; (b) add noise and a lowpass trend to it to obtain the total response; (c) recover the HRF from this total response, and (d) reconstruct the activation signal by convolving the stimulus pattern with the recovered HRF. The MSE between the activation signal and the reconstructed version was used to measure the reconstruction accuracy.
Figure 3. Test setup. (a) The fMRI response as an LTI signal: HRF (i), stimulus pattern (ii), and the activation signal (iii) as (i) convolved with (ii). (b) In the GLM, confounds such as trends (ii) and noise (iii) are added to the response (i). (c) The total response: activation signal + confounds + noise. (d) The ForWaRDreconstructed HRF (i) compared to the original (ii). (e) The activation signal using the extracted (i) and original (ii) HRF. (f) The canonical HRF, HRF_{spm}, see Eq. (6).
2.3.1 Simulation of fMRI time signals
128 EPI scans of a subject at rest were acquired on a 3T Intera system (Philips Medical Systems, The Netherlands), with repetition time (TR) 3 s, image size 64 × 64 × 46 voxels of 3.5 × 3.5 × 3.5 mm^{3}. A total of 512 noise time signals were collected from a region of 8 × 8 × 8 voxels (see Fig. 4). A randomised stimulus pattern f(n), n = 1, ..., N was obtained by thresholding a vector of random values. The stimulus f(n) n = 1, ..., N was convolved with an impulse response function h(n), in this case the canonical haemodynamic response function HRF_{spm}, to describe the activation signal s(n). HRF_{spm}(t) (see Fig. 3f) is the difference of two gamma density functions (GDF),
Figure 4. Area sampled from EPI data: (a) transverse, (b) sagittal, (c) coronal.
where the GDF γ_{m, l}(t) has the form:
Four types of lowfrequency trends: flat, linear, sinusoidal and quadratic (see Fig. 5a) were added to the signal. The SNR of the time signals was set by choosing the standard deviations σ_{s }of the signal and σ_{e }of the noise, as well as the scalar m_{n }such that
Figure 5. (a) Lowfrequency trends: (i) flat, (ii) linear, (iii) sinusoidal, (iv) quadratic. (b) log_{10}(MSE) of noisy (×) and reconstructed (o) signals given the input SNR.
Trends with standard deviation σ_{t }> 0 were scaled by m_{t }so that m_{t}σ_{t }= m_{n}σ_{e}. The time signal in each voxel was the sum {EPI noise} + {activation} + {trend} (Fig. 3b). Tests included time signals with varying (a) input SNR, (b) lowfrequency trends, (c) repetition time (TR), and (d) response onset.
2.3.2 Reconstruction of activation signals
The HRF was extracted from each time signal by the ForWaRDbased routine, and the mean HRF was used to reconstruct the activation signal by convolving it with the stimulus pattern. The following settings of the ForWaRDroutine were varied: (a) type of frequency shrinkage, (b) levels of the wavelet transform, (c) waveletdomain threshold level, and (d) the wavelet basis. The default values were: SNR 0 dB, no trend, TR 2 s, onset delay 0 s, Tikhonov shrinkage, τ 0.1, decomposition level 3, θ 3, φ_{1 }Daubechies4, φ_{2 }Daubechies3 [31].
2.4 EventRelated fMRI experiments
An application of the HRF extraction method is demonstrated in an eventrelated fMRI experiment. HRF coeffcients were extracted from the fMRI data set, and HRFs were computed by fitting a model function to HRF coeffcients (both wholebrain and regionofinterest). These HRFs were then used in a subsequent GLMbased analysis.
2.4.1 FixedISI experiment
The subject in the MRI scanner had to make a fist on the appearance of a visual stimulus, and relax after 1 s. Cues were given on a white screen placed inside the scanner: a red disc was a cue to make a fist, a white disc meant that the subject had to rest. The experiment consisted of 156 scans, acquired as described in Sect. 3.1. Cues were given every 24 s (8 scans × 3 s, no jittering), starting at scan 2. Increased taskrelated activity was expected in the motor cortex, the premotor cortex, the supplementary motor area and the cerebellum.
The EPI data were denoised with a waveletbased technique [32], using SUREShrink in the wavelet domain [29]. Realignment, normalisation, and statistical analysis were done with SPM2 [2]. Slice timing correction was not applied. The design matrix X contained a set of 6 Fourier basis functions (3 sines, 3 cosines) modulated by a Hanning window, in the time interval of 8 scans after each stimulus, and a constant signal modelling the time series mean. The Fourier basis was used to model a large class of signals using only few assumptions.
The variance ratio [33] was computed in each voxel, and an Ftest was used to determine significance. False discovery rate (FDR) control [34] with q = 0.05 was used to correct for multiple hypothesis testing.
2.4.2 RandomISI experiment
A second experiment used a random stimulus sequence, created by thresholding a sequence of random values (39 stimuli, ISI mean: 6.05 scans, σ : 4.31 scans, no jittering). The number of scans was 256, scanning parameters and preprocessing were as in Sect. 3.1.
Due to overlapping responses, the windowed Fourier basis set could not be used. It was replaced by the canonical HRF from the SPM2 program and its time and dilation derivatives. An Ftest of the variance ratio with FDR q = 0.05 was used to assess significance.
2.4.3 Extracting and modelling HRFs from the time series
A stricter FDRcorrected threshold (q = 0.0001) was applied to the variance ratio maps, sampling the HRF in a smaller number of voxels.
A continuous HRF was obtained by fitting a generalised version of HRF_{spm }(6) to the coeffcients returned by ForWaRD and selective averaging. We use HRF_{gam }to denote the difference of two GDFs with 8 parameters, i.e., H(eight), D(ilation), P(eak location) and L(ag) of both GDFs:
LevenbergMarquardt's nonlinear curvefitting algorithm was used to determine these parameters and 95% prediction intervals for the fitted functions.
2.4.4 Using fitted HRFs to model responses
HRFs measured from an fMRI data set cannot be used to test for activation in that same data set: a model must be specified a priori, so that inferences are not made from models that are determined by the data itself. We tested for activation in the randomISI experiment with the HRFs fitted to the coeffcients extracted from the fixedISI data, and vice versa. The GLM was estimated using the stimulus times and the fitted HRFs, and their correlation the fMRI time signals was computed. Significance was determined via a onesample ttest, and FDR control with q = 0.05 was used to correct for multiple hypothesis testing.
3 Results
3.1 Simulated fMRI time signals
The properties of the signal and parameters of the extraction routine that noticeably influenced the original activation signal s(n) and its reconstruction r(n), are listed below.
3.1.1 Output MSE as a function of input SNR
Figure 5b shows the MSE for various input SNR values. For input SNRs up to 9 dB the MSE decreases; it increases above 9 dB.
3.1.2 Choice of shrinkage type and τ parameter
Wiener shrinkage uses an iterative algorithm [25] to estimate H^{2 }(see (4)), the number of iterations was limited to 10. Figure 6 shows the MSE for both types of frequency domain shrinkage, varying SNR and regularisation parameter τ. For low SNR and heavy regularisation (τ ≥ 1), Tikhonov regularisation outperforms Wiener shrinkage. For higher SNRs or mild regularisation, Wiener shrinkage performs better. The best value for τ depends on the shrinkage type, the SNR and the TR. For short TR, mild regularisation (τ ≤ 0.1) is preferable. A long TR requires heavy regularisation.
Figure 6. Output MSEs of Tikhonov (left) and Wiener (right) shrinkage with a TR of 0.5 s (top) and 3 s (bottom), for different input SNRs and a varying regularisation parameter τ × τ = 0.01, o τ = 0.1, □: τ = 1, *: τ = 10.
3.1.3 Different response delays
Figure 7 shows that HRFs with different onset delays can equally well be extracted with ForWaRD. The MSE hardly changes with different delays, indicating that the shape of the response is preserved. The increased MSE for negative shifts is caused by the fact that the HRF was only sampled poststimulus.
Figure 7. Output MSE for varying response onset delays, for Tikhonov (left) and Wiener (right) shrinkage, SNR = 2 dB (×), 0 dB (o), 2 dB (□), 4 dB (*), and 6 dB (+).
3.1.4 Different wavelet filters
We tested 15 different wavelet filters for (φ_{1}, ψ_{1}), as well as for (φ_{2}, ψ_{2}): Daubechies wavelets 1 ... 5 (the filter number indicates the number of vanishing moments), Daubechies' symmetric wavelets 2 ... 6 [35] (filter 1 corresponds to the Daubechies1 filter), and Coiflets 1 ... 5 [35]. Different filters did not yield large differences in performance.
3.1.5 Decomposition level and noise threshold
Figure 8 shows the MSE for different SNRs, different θ and different decomposition levels. We find that twolevel decompositions produce the smallest errors for the lower SNRs, and threelevel decompositions perform best for the higher. Fourlevel and fivelevel decompositions yield higher errors. A higher θ often yields a lower MSE.
Figure 8. Output MSE for different SNRs, with various levels of decomposition and threshold levels. Left: θ = 2, right: θ = 3, with a twolevel transform (×), threelevel(o), fourlevel (□), fivelevel(*).
3.1.6 Different lowfrequency trends
Tests with lowfrequency trends (see Fig. 5a) show that the type of frequency domain shrinkage changes the result. (see Fig. 9). The MSE is higher with Tikhonov than with Wiener shrinkage, especially for lower SNRs. Trends are not removed perfectly with either shrinkage type, but the extra information about f and the power spectral density of h in Wiener shrinkage make it less sensitive to trend residuals.
Figure 9. Output MSE for different SNRs and lowpass trends in the data. Left: Tikhonov shrinkage, right: Wiener shrinkage. No trend (×), a linear trend (o), a sinusoidal trend (□), or a quadratic trend (*).
3.1.7 Summary
Even though the default parameters (see Sect. 2.3.2) make the algorithm robust for a wide range of signals (different SNR, sampling frequency, etc.), the method is also robust changing its parameters, such as the wavelet filters and decomposition level. The MSE of the reconstructed signal was lower than the input MSE in most of the tested situations.
The robustness to onset delay makes ForWaRD an attractive alternative to other delay correction methods such as including a temporal derivative of the standard HRF in the model [16]: these only correct for small synchronisation errors, and are not usable when the HRF is not well modelled by the standard function.
3.2 Eventrelated fMRI: fixedISI
Activation maps are shown in Fig. 10a as 'glass brain' maximum intensity projections (MIP) in two orthogonal directions. Low statistic values are shown in grey, high values in black. The highest value is indicated with a '<' sign. Activation was found in the expected areas, predominantly in the motor cortex. A wholevolume HRF was extracted from the poststimulus volumes using selective averaging [9], by taking the mean response of each volume, weighted by the map of significant statistic values. A regionspecific HRF in the 7 × 7 × 7voxel neighbourhood of a selected voxel (see the '<' in Fig. 10a) was computed by using only the time signals from that region. Figure 11a–b shows the extracted HRFs.
Figure 10. Maps of the variance ratio using the Ftest with FDR q = 0.05. The highest value is indicated with a '<' sign. (a) FixedISI experiment, range 3.86–20.63 and (b) randomISI experiment, range 5.43–41.63 (b). The indicated areas are: l.motor cortex (1), supplementary motor area (2), premotor area (3), r.cerebellum (4).
Figure 11. HRFs extracted from the fixedISI data by selective averaging (top row) and ForWaRD (bottom row). Left: wholevolume, right: regionspecific. The extracted coeffcient are the × at each TR. Dotted lines: fits of HRF_{gam }to the coeffcients. Error bars show the 95% confidence intervals for the fitted function.
The ForWaRD algorithm used 128 scans of the experiment, starting with scan 2 (first stimulus). Wholevolume HRF and regionspecific HRF time points were computed using the poststimulus time series (see Fig. 11c–d). The HRFs extracted by ForWaRD are similar to the results from selective averaging, except that the baseline of the ForWaRDextracted HRF decreases. This is because the HRF does not return to baseline within the sampled interval (24 s), so in the GLM the response decreases at the end of every stimulus. In some cases (see Fig. 11c–d) this was fixed manually by setting the baseline to the last HRF coeffcient.
3.3 Eventrelated fMRI: randomISI
The contrast of the randomISI experiment was generally weaker than in the fixedISI experiment, but localisation was better, see Fig. 10b.
A poststimulus time series was made with ForWaRD, using all 256 scans of the experiment. Selective averaging could not be used because of overlapping responses. With a random ISI, a much longer poststimulus interval can be sampled [11]. The poststimulus volumes produced by the extraction routine were used to create a wholevolume HRF and a regionspecific HRF in the same way as in the fixedISI case.
Figure 12 shows the extracted HRF coeffcients and the fitted HRFs. A comparison between Fig. 11c–d and Figure 12 shows that the ForWaRDextracted HRF coeffcients from the randomISI design agree better with the model function than those from the fixedISI design, especially in the regionspecific case. A possible explanation is that the fixedISI stimulus signal does not contain enough frequency information to do the Fourier inversion, resulting in a lowerquality estimate. In contrast, selective averaging does not work with overlapping responses and random ISIs.
Figure 12. HRFs extracted from the randomISI experiment by ForWaRD: wholevolume (a) and regionspecific (b). ×: extracted HRF coeffcients. Dashed lines: function HRF_{gam }fitted to ×. 95% prediction intervals for the fitted functions are shown as error bars.
3.4 Activation detected with fitted HRFs
Figure 13 shows the activation detected in the fixedISI dataset with the HRFs extracted by ForWaRD from the randomISI dataset. There is good correspondence between the maps of the wholevolume (a) and regionspecific (b) HRFs, respectively, and the detected activations match the expected pattern (see Fig. 10). Figure 14 shows the activation detected in the randomISI dataset with the HRFs from the fixedISI dataset, using both selective averaging (ab) and ForWaRD (cd), which are in very good agreement. Where selective averaging is possible, ForWaRD yields results very similar to those produced by selective averaging. The statistical maps for ForWaRD between the fixedISI and randomISI time series are in good correspondence.
Figure 13. Activation in the fixedISI data set, using the Ftest with FDR q = 0.05. The highest value is indicated with a '<' sign. HRFs extracted from the randomISI data set by ForWaRD and modelled with HRF_{gam}. (a) wholevolume HRF, range 3.08–10.62. (b) regionspecific HRF, range 2.99–12.72.
Figure 14. Activation maps of the randomISI data set, using the Ftest with FDR q = 0.05. The highest value is indicated with a '<' sign. HRFs modelled by HRF_{gam}, coeffcients extracted from the fixedISI data set by selective averaging with wholevolume HRF (a, range 3.10–10.11) and regionspecific HRF (b, range 3.10–10.21), and ForWaRD with wholevolume HRF (c, range 3.11–10.15) and regionspecific HRF (d, range 3.13–10.10).
An analysis was also performed performed with the general HRF_{spm}. The activation maps in Fig. 13 are in very good correspondence with Fig. 15a, indicating that both HRF_{spm }and HRF_{gam }model the data well. The values for the variance ratio in Fig. 16 for the fixedISI experiment show that all models explain the data well, and HRF_{gam }with the regional HRF yields the best fit to the data.
Figure 15. Activation maps of the fixedISI data set (a, range 3.10–10.80), and the randomISI data set (b, range 3.32–8.61), using HRF_{spm}. The highest value is indicated with a '<' sign.
Figure 16. Maximum variance ratio values found in the tests.
Figures 14 and 15b show that HRF_{spm }yields a poor fit to the data; this is also shown in Fig. 16. The higher variance ratios for HRF_{gam }indicate that a larger portion of the measured variance is explained by the model, and that the residual of the GLM (1) is small. The fitted regionspecific HRFs generally perform better than wholevolume HRFs, and the maps of detected activation indicate that the fitted HRFs do not only detect activation in the region from which they were extracted, but that they are general enough to also detect activation in other areas, corresponding to the predicted activated regions (see the blue ellipses in Fig. 10).
4 Conclusion
We have developed a technique to extract HRF coeffcients from fMRI time series based on the ForWaRD deconvolution technique. Frequencydomain deconvolution allows extraction of the HRF even when the responses to subsequent stimuli overlap, and the sensitivity to noise of frequencydomain deconvolution is compensated by Wiener or Tikhonov shrinkage in the frequency domain, followed by waveletdomain Wiener shrinkage. Before applying ForWaRD, lowfrequency trends are removed from the time signal with a standard waveletbased method. Tests of the extraction routine using simulated activation, several types of trend and noise from a real fMRI time series, demonstrate its robustness. Results show that the method is robust to trends in the data, and the performance does not differ much between the noise levels we tested. The output of our algorithm is a poststimulus time series, representing the HRF coeffcients in every voxel. At present, the extraction method is capable of recovering one HRF from one time series. The algorithm used in this paper is entirely datadriven: it does not use a priori models of the data. Other advantages of this algorithm are its simplicity, i.e., the algorithm works independently of other pre and postprocessing steps; and its speed and low computational complexity.
The default settings of the method used in the tests as well as in the fMRI experiments (see Sect. 2.3.2) lead to a good and robust performance of the extraction algorithm.
We have demonstrated the use of ForWaRDextracted HRFs in combination with continuous HRF functions to predict eventrelated fMRI responses. Given the output of the extraction routine, continuous functions (in this case gamma densities) are fitted to the average HRF coeffcients in a region of a (statistical or anatomical) map.
Results from the experiments with extracted and fitted HRFs indicate that subjectspecific and regionspecific HRFs lead to stronger contrasts and better localisation than a standard HRF. The results with the randomISI data suggest that it is possible to use ForWaRD in combination with relatively complex prior experiments to extract HRFs, and that it is beneficial to use these HRFs instead of standard functions for detection and modelling of subsequently acquired data. The advantage of using a single, tailored HRF over a model that spans several basis functions, is that the statistical analysis is more specific and more powerful, resulting in a stronger contrast and better localisation.
A possible extension of the method is the extraction of multiple HRFs from one or multiple experiments. Deconvolution in the frequency domain of multiple waveforms has already been done in ERP research [36]. The ForWaRDbased method may be extended in a similar way.
Competing interests
The author(s) declare that they have no competing interests.
Authors' contributions
AMW wrote the software, did the analyses and wrote most of the paper. JBTMR designed critical parts of the method and coauthored the paper. HH led the image acquisition.
Note
Figure 16 Maximum variance ratio values found in the tests.
Acknowledgements
This research is part of the project "Wavelets and their applications", funded by the Dutch National Science Foundation (NWO), project no. 613.006.570.
References

Ogawa S, Lee TM, Kay AR, Tank DW: Brain Magnetic Resonance Imaging with Contrast Dependent on Blood Oxygenation.
Proceedings of the National Academy of Sciences 1990, 87:98689872. Publisher Full Text

Friston KJ, Holmes AP, Worsley KJ, Poline JP, Frith CD, Frackowiak RSJ: Statistical Parametric Maps in Functional Imaging: A General Linear Approach. [http://www.fil.ion.ucl.ac.uk/spm] webcite
Human Brain Mapping 1995, 2:189210. Publisher Full Text

Boynton GM, Engel SA, Glover GH, Heeger DJ: Linear Systems Analysis of Functional Magnetic Resonance Imaging in Human V1.
The Journal of Neuroscience 1996, 16(13):42074221. PubMed Abstract  Publisher Full Text

Ciuciu P, Poline JB, Marrelec G, Idier J, Pallier C, Benali H: Unsupervised robust nonparametric estimation of the hemodynamic response function for any fMRI experiment.
IEEE Transactions on Medical Imaging 2003, 22(10):12241234. PubMed Abstract  Publisher Full Text

Handwerker DA, Ollinger JM, D'Esposito M: Variation of BOLD hemodynamic responses across subjects and brain regions and their effects on statistical analyses.
NeuroImage 2004, 21(4):16391651. PubMed Abstract  Publisher Full Text

Woolrich MW, Jenkinson M, Brady JM, Smith SM: Fully Bayesian spatiotemporal modeling of FMRI data.
IEEE Transactions on Medical Imaging 2004, 23(2):213231. PubMed Abstract  Publisher Full Text

Neumann J, Lohmann G, Zysset S, von Cramon DY: Withinsubject variability of BOLD response dynamics.
Neuroimage 2003, 19(3):784796. PubMed Abstract  Publisher Full Text

Bandettini PA, Cox RW: EventRelated fMRI Contrast When Using Constant Interstimulus Interval: Theory and Experiment.
Magnetic Resonance in Medicine 2000, 43:540548. PubMed Abstract  Publisher Full Text

Buckner RL, Bandettini PA, Craven KMO, Savoy RL, Petersen SE, Raichle ME, Rosen BR: Detection of cortical activation during averaged single trials of a cognitive task using functional magnetic resonance imaging.
Proceedings of the National Academy of Sciences 1996, 93:1487814883. Publisher Full Text

Zarahn E, Aguirre G, D'Esposito M: A TrialBased Experimental Design for fMRI.
NeuroImage 1997, 6:122138. PubMed Abstract  Publisher Full Text

Burock MA, Buckner RL, Woldorff MG, Rosen BR, Dale AM: Randomized EventRelated Experimental Designs Allow for Extremely Rapid Presentation Rates Using Functional MRI.
NeuroReport 1998, 9:37353739. PubMed Abstract  Publisher Full Text

Glover GH: Deconvolution of Impulse Response in EventRelated BOLD fMRI.
NeuroImage 1999, 9:416429. PubMed Abstract  Publisher Full Text

Miezin FM, Macotta L, Ollinger JM, Petersen SE, Buckner RL: Characterizing the Hemodynamic Response: Effects of Presentation Rate, Sampling Procedure and the Possibility of Ordering Brain Activity Based On relative timing.
NeuroImage 2000, 11:735759. PubMed Abstract  Publisher Full Text

Hinrichs H, Scholz M, Tempelmann C, Woldorff MG, Dale AM, Heinze HJ: Deconvolution of EventRelated fMRI Responses in FastRate Experimental Designs: Tracking Amplitude Variations.
Journal of Cognitive Neuroscience 2000, 12(6(suppl 2)):7689. PubMed Abstract  Publisher Full Text

Josephs O, Turner R, Friston KJ: EventRelated fMRI.
Human Brain Mapping 1997, 5:243248. Publisher Full Text

Friston KJ, Fletcher P, Josephs O, Holmes A, Rugg MD, Turner R: Characterizing Differential Responses.
NeuroImage 1998, 7:3040. PubMed Abstract  Publisher Full Text

Neelamani R, Choi H, Baraniuk RG: ForWaRD: FourierWavelet Regularized Deconvolution for IllConditioned Systems.
IEEE Transactions on Signal Processing 2004, 52(2):418433. Publisher Full Text

Johnstone I, Kerkyacharian G, Picard D, Raimondo M: Wavelet deconvolution in a periodic setting.
Journal of the Royal Statistical Society 2004, 66(3):547573. Publisher Full Text

Kalifa J, Mallat S, Rouge B: Deconvolution by thresholding in mirror wavelet bases.
IEEE Transactions on Image Processing 2003, 12(4):446457. PubMed Abstract  Publisher Full Text

SanchezAvila C: Wavelet domain signal deconvolution with singularitypreserving regularization.

Bullmore ET, Long C, Suckling J, Fadili J, Calvert G, Zelaya F, Carpenter TA, Brammer M: Colored Noise and Computational Inference in Neurophysiological (fMRI) Time Series Analysis: Resampling Methods in Time and Wavelet Domains.
Human Brain Mapping 2001, 12:6178. PubMed Abstract  Publisher Full Text

Zarahn E, Aguirre GK, D'Esposito M: Empirical Analyses of BOLD fMRI Statistics: I. Spatially Unsmoothed Data Collected under NullHypothesis Conditions.
NeuroImage 1997, 5:179197. PubMed Abstract  Publisher Full Text

Logothetis NK, Pauls J, Augath M, Trinath T, Oeltermann A: Neurophysiological Investigation of the Basis of the fMRI Signal.
Nature 2001, 412:150157. PubMed Abstract  Publisher Full Text

Kilmer ME, O'Leary DP: Choosing Regularization Parameters in Iterative Methods for IllPosed Problems.
SIAM Journal on Matrix Analysis and Applications 2000, 22(4):12041221. Publisher Full Text

Hillery AD, Chin RT: Iterative Wiener Filters for image restoration.
IEEE Transactions on Signal Processing 1991, 39:18921899. Publisher Full Text

Ghael SP, Sayeed AM, Baraniuk RG: Improved Wavelet Denoising via Empirical Wiener Filtering.
Proc SPIE: Wavelet Applications in Signal and Image Processing 1997, 3169:389399.

Mallat SG: A Theory for Multiresolution Signal Decomposition: The Wavelet Representation.
IEEE Transactions on Pattern Analysis and Machine Intelligence 1989, 11(7):674693. Publisher Full Text

Holschneider M, KronlandMartinet R, Morlet J, Tchamitchian P: A realtime algorithm for signal analysis with the help of the wavelet transform.
In Wavelets, TimeFrequency Methods and Phase Space Edited by Combes JM, Grossmann A, Tchamitchian P. 1989, 286297.

Donoho DL, Johnstone IM: Adapting to Unknown Smoothness by Wavelet Shrinkage.
Journal of the American Statistical Association 1995, 90:12001224. Publisher Full Text

Meyer FG: Waveletbased estimation of a semiparametric generalized linear model of fMRI timeseries.
IEEE Transactions on Medical Imaging 2003, 22(3):315322. PubMed Abstract  Publisher Full Text

Daubechies I: Orthonormal Bases of Compactly Supported Wavelets.
Communications on Pure and Applied Mathematics 1988, 41:909996. Publisher Full Text

Wink AM, Roerdink JBTM: Denoising functional MR images: a comparison of wavelet denoising and Gaussian smoothing.
IEEE Transactions on Medical Imaging 2004, 23(3):374387. PubMed Abstract  Publisher Full Text

Josephs O, Henson RNA: Eventrelated Functional Magnetic Resonance Imaging: Modelling, Inference and Optimization.
Philosophical Transactions of the Royal Society 1999, 354:12151228. Publisher Full Text

Genovese CR, Lazar NA, Nichols TE: Thresholding of Statistical Maps in Functional Neuroimaging Using the False Discovery Rate. [http://www.sph.umich.edu/~nichols/FDR] webcite
NeuroImage 2002, 15:772786. PubMed Abstract  Publisher Full Text

Daubechies I: Orthonormal Bases of Compactly Supported Wavelets: II. Variations on a Theme.
SIAM Journal on Mathematical Analysis 1993, 24(2):499519. Publisher Full Text

Hansen JC: Separation of Overlapping Waveforms Having Known Temporal Distributions.
Journal of Neuroscience Methods 1983, 9:127139. PubMed Abstract  Publisher Full Text
Prepublication history
The prepublication history for this paper can be accessed here: