fMRI hemodynamic response function estimation in autoregressive noise by avoiding the drift

Abstract

The measured functional magnetic resonance imaging (fMRI) time series is typically corrupted by instrumental drift and physiological noise due to respiration and heartbeat giving rise to temporal correlation in the signals. Most methods proposed so far for nonparametric hemodynamic response function (HRF) estimation in fMRI data do not account for these confounding effects, and thus produce biased and inefficient estimates. The aim of this paper is to address this issue by modeling the noise in fMRI time series using an autoregressive model of order p ($AR(p)$). Making use of a semiparametric model to characterize the fMRI time series and the $AR(p)$ to model the temporally correlated noise, a generalized least squares (GLS) estimator for voxelwise consistent nonparametric HRF estimation is derived. The proposed estimation method is a three-stage procedure that relies on first-order differencing to remove drift and a novel structured covariance estimator for the AR noise based on Cholesky decomposition to derive the best linear unbiased estimator (BLUE) of the HRF. We also establish the asymptotic consistency of the proposed estimator. Simulation results show that the proposed method generates more accurate HRF estimates compared to existing methods. When applied to real fMRI data, it demonstrates the effectiveness in uncovering the brain response temporal dynamics for both event-related and block-design paradigms. Our approach removes the two types of noise in fMRI data simultaneously, thus providing efficient estimation of brain hemodynamic responses, while allowing for flexible characterization of the shape and timing of the voxelwise HRF.

Introduction

Functional magnetic resonance imaging (fMRI) has been widely used in neuroscience studies to analyze cognitive functions in human brain. Besides the use in detecting activated brain areas in response to a specific stimulus or task, finding the temporal dynamics of the brain response from fMRI during activations is also an important problem to address. The blood oxygen level-dependent (BOLD) fMRI relies on the coupling between increases in neuronal activity and increases in blood flow and volume that accompany the local increase in oxygen demand to offer a proxy of the underlying local neuronal activity [1]. Identifying the temporal dynamics of brain responses during activation is among reasons why development of HRF estimation methods for the BOLD fMRI signal in a particular voxel has attracted a lot of attention. The BOLD response is usually, as a first approximation, modeled as a convolution of a stimulus function based on the experimental paradigm, with a hemodynamic response function (HRF) characterizing the temporal dynamics of the brain response. The key assumptions related to the convolution model are the stationarity and linearity of the neurovascular system for which some evidences have been provided in [2], [3], [4]. Despite the flexibilities of non-linear modeling [5], the ability to assume linearity is important to allow for the use of a linear time invariant model, which can provide more robust estimation and interpretable characterizations in noisy systems [6]. This is related to critical issues that HRF estimation is typically complicated by the confounding effect of various noise sources in the fMRI data which cause bias in the estimates, mainly from (1) low-frequency drift due to instrumental instabilities, and (2) oscillatory noise due to respiration and cardiac pulsation. In this paper, we address the important question of how to obtain the most efficient estimates of the HRF with the least amount of bias and misspecification, particularly in the presence of baseline drift and temporally-correlated physiological noise in fMRI data.

Common approach to estimating the HRF relies on parametric modeling of the HRF and the drift. This approach imposes a priori shapes on the HRF with some families of statistical distributions such as Gamma function, and on the drift with slowly-varying parametric models such as weighted discrete cosine transform (DCT) basis functions [7], [8], [9], [10], polynomials [11] or splines [12], which often assumes unrealistically the same model to fit at each voxel, under all conditions and subjects. While parametric models are very useful in providing a parsimonious description, they have drawbacks of introducing modeling biases [13] and lack of flexibility. Moreover, many fitting procedures for parametric models are computationally expensive. The alternative nonparametric approach uses a set of basis functions without the constraints of any arbitrary modeling structures, and thus allowing for a more flexible representation of a broader class of HRF and drift signals compared to using parametric models in HRF estimation [7], [8], [9], [10], [14], [15]. Allowing flexibility in both the HRF and the drift across different regions, conditions and subjects, will help reduce the bias due to model misspecification and hence more accurate HRF estimates. Moreover, this approach avoids the selection of nuisance covariates and the estimation of their coefficients, and may also capture physiologically-meaningful parameters in the HRF. We adopt a non-parametric approach proposed in our recent works [16], [17], [18], which estimates the HRF directly as unknown parameters in the linear model, by least-square fitting of signals with first-order differencing to remove the drift effect, without using any basis functions. The use of nonparametric component for the drift leads to a flexible semi-parametric estimation of the underlying HRF [19], [20].

Another key ingredient of fMRI time series models besides the drift is specification of the additive noise. The convenient assumption of temporally independent noise is inappropriate, as the residuals of fMRI signals were shown to be non-white, but colored noise exhibiting temporal correlation arising from the aliased physiological artifacts [21], [22], [23]. Ignoring this autocorrelation in fMRI data renders the ordinary least squares (LS) estimate of HRF parameters asymptotically inefficient, and worse introduces negative bias in the estimated parameter standard errors, resulting in invalid statistical inference. To overcome this problem, most fMRI analyses adopted the pre-whitening strategy which first estimates the autocorrelations by assuming parametric models for the colored noise, e.g. autoregressive process of order one, $AR(1)$ [24], [25], [26], [27] or higher orders, $AR(p)$ [11], [21], [28], [29], and then use the estimated parameters to ‘pre-whiten’ or decorrelate the errors to produce an efficient HRF estimate. However, these studies still relied on the parametric estimation of HRF and drift with the drawbacks discussed above.

In this paper, we propose a unified nonparametric framework for consistent estimation of hemodynamic response by simultaneously dealing with the temporally correlated noise and avoiding drift present in fMRI data. More precisely, we develop a novel nonparametric estimator of HRF by first-order differencing of the fMRI signals to remove the drift, and then derive a generalized LS (GLS) estimator to account for the temporal autocorrelations in the noise. The proposed GLS estimator is an extension of the LS fitting of differenced fMRI signals for the white noise case considered in our previous approach [16], [17], [18]. It is obtained by incorporating the covariance structure of the colored noise to further improve the HRF estimation. The GLS estimator is a best linear unbiased estimator (BLUE) (with minimum variance) given a known noise covariance. However, it is typically unknown and current methods are unsatisfactory in providing an accurate estimate to reduce the bias. In order to construct the BLUE of HRF, we further derive a consistent estimator for the noise covariance based on drift-free residuals and a parameterization of the covariance structure by Cholesky decomposition assuming an $AR(p)$ noise model. The entries in the Cholesky decomposition of the covariance matrix can be interpreted as AR parameters and innovation variances. Therefore, the proposed estimation procedure consists of three steps: [Step 1.] The drift-corrected, unbiased residuals are generated from the LS fitting of the non-parametric fMRI model to the differenced signals. [Step 2.] The covariance estimator is then constructed by substituting in the Cholesky factors with the AR parameters estimated based on the differenced residuals in Step 1. [Step 3.] A novel GLS estimator of HRF with improved asymptotic efficiency (lower bias and variance) can be computed from the necessary statistics obtained in the first two steps. To our knowledge, there were limited studies on the asymptotic properties of HRF estimation. We provide the asymptotic consistency of the proposed estimator which is important for theoretical inferences.

Our proposed GLS estimator which integrates the estimated noise covariance directly in the LS fitting, is more general than using it separately in an initial stage for pre-whitening. A similar GLS method was proposed for non-parametric HRF estimation in [30] which, however, has not taken into account the drift effect and only used an $AR(1)$ model for noise. In contrast, our approach utilizes the drift-corrected signals and a consistent noise covariance estimate based on the Cholesky factorization of $AR(p)$ noise structure, thus producing more efficient HRF estimates. It is also computationally more efficient than the voxel-wise expectation maximization (EM) estimation in [31], [27], since the estimates of the large covariance matrix and its inverse are simple constructs of the AR Cholesky factors. Besides, the derivation of our GLS estimator is not straightforward. It relies on the covariance of a differenced $AR(p)$ noise process equivalent to an AR-moving average ($ARMA$) process, instead of an $AR(p)$ process in the standard GLS.

The rest of the paper is organized as follows. The proposed voxel-wise semiparametric fMRI time series model is presented in Section 2. The three-step procedure to obtain the novel GLS HRF estimator is described in Section 3. The performance evaluation of the proposed HRF estimator via simulations and real fMRI data for both block-design and event-related paradigms are given in Sections 4 and 5 respectively. The conclusion is given in Section 6.