Modeling the Hemodynamic Response Function Using EEG-fMRI Data During Eyes-Open Resting-State Conditions and Motor Task Execution

Abstract

Being able to accurately quantify the hemodynamic response function (HRF) that links the blood oxygen level dependent functional magnetic resonance imaging (BOLD-fMRI) signal to the underlying neural activity is important both for elucidating neurovascular coupling mechanisms and improving the accuracy of fMRI-based functional connectivity analyses. In particular, HRF estimation using BOLD-fMRI is challenging particularly in the case of resting-state data, due to the absence of information about the underlying neuronal dynamics. To this end, using simultaneously recorded electroencephalography (EEG) and fMRI data is a promising approach, as EEG provides a more direct measure of neural activations. In the present work, we employ simultaneous EEG-fMRI to investigate the regional characteristics of the HRF using measurements acquired during resting conditions. We propose a novel methodological approach based on combining distributed EEG source space reconstruction, which improves the spatial resolution of HRF estimation and using block-structured linear and nonlinear models, which enables us to simultaneously obtain HRF estimates and the contribution of different EEG frequency bands. Our results suggest that the dynamics of the resting-state BOLD signal can be sufficiently described using linear models and that the contribution of each band is region specific. Specifically, it was found that sensory-motor cortices exhibit positive HRF shapes, whereas the lateral occipital cortex and areas in the parietal cortex, such as the inferior and superior parietal lobule exhibit negative HRF shapes. To validate the proposed method, we repeated the analysis using simultaneous EEG-fMRI measurements acquired during execution of a unimanual hand-grip task. Our results reveal significant associations between BOLD signal variations and electrophysiological power fluctuations in the ipsilateral primary motor cortex, particularly for the EEG beta band, in agreement with previous studies in the literature.

Appendices

Appendix

Identification of Hammerstein and Wiener-Hammerstein Models

The Hammerstein model (Fig. 1a; Eq. (1)) can be estimated efficiently from the input–output data using orthonormal basis functions for the representation of the LTI block (Gómez and Baeyens 2004), which is given by

$${\text{h}}\left( {\text{m}} \right) = \sum\limits_{{{\text{j}} = 0}}^{{{\text{L}} - 1}} {{\text{b}}_{{\text{j}}} } {\text{B}}_{{\text{j}}} \left( {\text{m}} \right)$$

(6)

where $\left\{{\text{B}}_{\text{j}}\left(\text{n}\right);\text{j}=0,\dots ,\text{ L}-1;\text{n}=0,\dots ,\text{M}\right\}$ is a set of $\text{L}$ orthonormal basis functions, and ${\text{b}}_{\text{j}}$ is the unknown expansion coefficient of the $\text{j}$-th order basis function. The use of orthonormal bases reduces the number of required free parameters in the model and allows parameter estimation using least-squares regression. This leads to increased estimation accuracy in the presence of noise even from short experimental data-records. Combining Eqs. (1, 2) and (6), the input–output relationship can be written as

$$\begin{array}{*{20}l} {{\text{y}}\left( {\text{n}} \right) = \sum\limits_{{{\text{m}} = 0}}^{{\text{M}}} {\left( {\sum _{{{\text{j}} = 0}}^{{{\text{L}} - 1}} {\text{b}}_{{\text{j}}} {\text{B}}_{{\text{j}}} \left( {\text{m}} \right)} \right)} \left( {\sum _{{{\text{i}} = 1}}^{{\text{Q}}} \sum _{{{\text{p}} = 1}}^{{\text{P}}} {\text{a}}_{{{\text{i}},{\text{p}}}} {\text{g}}_{{\text{i}}}^{{\left( {\text{p}} \right)}} \left( {{\text{u}}_{{\text{i}}} \left( {\text{n}} \right)} \right)} \right) + \varepsilon \left( {\text{n}} \right)} \hfill \\ \quad\quad\; { = \sum\limits_{{{\text{j}} = 0}}^{{{\text{L}} - 1}} {\sum\limits_{{{\text{i}} = 1}}^{{\text{Q}}} {\sum\limits_{{{\text{p}} = 1}}^{{\text{P}}} {{\text{a}}_{{{\text{i}},{\text{p}}}} } } } {\text{b}}_{{\text{j}}} \cdot \sum _{{{\text{m}} = 0}}^{{\text{M}}} {\text{B}}_{{\text{j}}} \left( {\text{m}} \right){\text{g}}_{{\text{i}}} \left( {{\mathbf{u}}\left( {{\text{n}} - {\text{m}}} \right)} \right) + \upvarepsilon \left( {\text{n}} \right)} \hfill \\\quad\quad\; { = \sum\limits_{{{\text{j}} = 0}}^{{{\text{L}} - 1}} {\sum\limits_{{{\text{i}} = 1}}^{{\text{Q}}} {\sum\limits_{{{\text{p}} = 1}}^{{\text{P}}} {{\text{a}}_{{{\text{i}},{\text{p}}}} } } } {\text{b}}_{{\text{j}}} \cdot {\text{v}}_{{{\text{i}},{\text{j}}}}^{{\left( {\text{p}} \right)}} \left( {\text{n}} \right) + \upvarepsilon \left( {\text{n}} \right)} \hfill \\\quad\quad\; \end{array}$$

(7)

where ${\text{v}}_{\text{i},\text{j}}^{\left(\text{p}\right)}\left(\text{n}\right)$ denotes the convolution between the $\text{p}$-th polynomial power of the $\text{i}$-th input with the $\text{j}$-th basis function. Equation (7) can be re-expressed as a linear regression problem

$$\mathbf{y}=\mathbf{V}\mathbf{c}+{\varvec{\upvarepsilon}}$$

(8)

where $\mathbf{V}=\left[{\text{v}}_{\text{1,0}}^{\left(1\right)},\dots , {\text{v}}_{\text{Q},0}^{\left(1\right)}, \dots ,{\text{v}}_{1,\text{L}-1}^{\left(1\right)},\dots ,{\text{v}}_{\text{Q},\text{L}-1}^{\left(1\right)} ,\dots \dots , {\text{v}}_{\text{1,0}}^{\left(\text{P}\right)},\dots , {\text{v}}_{\text{Q},0}^{\left(\text{P}\right)}, \dots ,{\text{v}}_{1,\text{L}-1}^{\left(\text{P}\right)},\dots ,{\text{v}}_{\text{Q},\text{L}-1}^{\left(\text{P}\right)}\right]$, and ${\mathbf{c}} = \left[ {{\text{a}}_{{1,1}} {\text{b}}_{0} , \ldots ,{\text{a}}_{{1,1}} {\text{b}}_{{{\text{L}} - 1}} , \ldots ,{\text{a}}_{{{\text{Q}},1}} {\text{b}}_{0} , \ldots ,{\text{a}}_{{{\text{Q}},1}} {\text{b}}_{{{\text{L}} - 1}} , \ldots \ldots ,{\text{a}}_{{1,{\text{P}}}} {\text{b}}_{0} , \ldots ,{\text{a}}_{{1,{\text{P}}}} {\text{b}}_{{{\text{L}} - 1}} , \ldots ,{\text{a}}_{{{\text{Q,P}}}} {\text{b}}_{0} , \ldots ,{\text{a}}_{{{\text{Q,P}}}} {\text{b}}_{{{\text{L}} - 1}} } \right]^{{\text{T}}}$ is a vector of the unknown model parameters.

Power fluctuations within distinct EEG frequency bands are highly correlated as previously reported in the literature (de Munck et al. 2009). As a result, the columns in $\mathbf{V}$ are strongly collinear, which makes estimation of $\mathbf{c}$ using ordinary least-squares numerically unstable due to ill-conditioning of the Gram matrix $\left[{\mathbf{V}}^{\mathbf{T}}\mathbf{V}\right]$. Therefore, to obtain a numerically more stable estimate of the unknown parameter vector $\mathbf{c}$, we employed partial least-squares regression (PLSR) (Rospiral and Kramer 2006).

PLSR is performed in three phases. In phase 1, the algorithm finds projections of $\mathbf{V}$ and $\mathbf{y}$ to a new co-ordinate system such that the covariance of these projections is maximized. This is achieved using a linear decomposition of both $\mathbf{V}$ and $\mathbf{y}$ into a set of orthonormal latent variables (scores) and loadings given by

$$\mathbf{V}=\mathbf{T}{\mathbf{P}}^{\mathbf{T}}+{\mathbf{e}}_{1}$$

(9)

$$\mathbf{y}=\mathbf{U}{\mathbf{W}}^{\mathbf{T}}+{\mathbf{e}}_{2}$$

(10)

where $\mathbf{T}$ and $\mathbf{U}$ are matrices of latent variables associated with $\mathbf{V}$ and $\mathbf{y}$, respectively. $\mathbf{P}$ and $\mathbf{W}$ are the corresponding loadings for each latent variable matrix, and ${\mathbf{e}}_{1,2}$ are error terms. The decomposition of $\mathbf{V}$ and $\mathbf{y}$ is performed such that the covariance between $\mathbf{T}$ and $\mathbf{U}$ is maximized. In phase 2, the algorithm performs ordinary least-squares regression analysis between the latent variables $\mathbf{T}$ and system output $\mathbf{y}$

$$\mathbf{y}=\mathbf{T}{\varvec{\uptheta}}+\mathbf{e}$$

(11)

$$\therefore {\widehat{{\varvec{\uptheta}}}}_{\mathbf{L}\mathbf{S}\mathbf{E}}={\left[{\mathbf{T}}^{\mathbf{T}}\mathbf{T}\right]}^{-1}{\mathbf{T}}^{\mathbf{T}}\mathbf{y}$$

(12)

where ${\varvec{\uptheta}}={\left[{\uptheta }_{1},\dots , {\uptheta }_{\text{Q}\times \text{L}\times \text{P}}\right]}^{\text{T}}$ is a vector of the regression coefficients. Note that in this case the Gram matrix $\left[{\mathbf{T}}^{\mathbf{T}}\mathbf{T}\right]$ is well-conditioned since the columns in $\mathbf{T}$ are orthonormal. In phase 3, the estimated ${\widehat{{\varvec{\uptheta}}}}_{\mathbf{L}\mathbf{S}\mathbf{E}}$ coefficients are projected back to the original parameter space yielding unbiased estimates of the original model parameters ${\widehat{\mathbf{c}}}_{\mathbf{P}\mathbf{L}\mathbf{S}}$ (die Jong 1993).

To uniquely identify the unknown parameters of the Hammerstein model described by (1, 2) and (6, 7), the bilinear parameter vector $\mathbf{c}$ needs to be dissociated into its constituent ${\text{a}}_{\text{i},\text{p}}$ and ${\text{b}}_{\text{j}}$ parameters. The parameter vector $\mathbf{c}$ can be reshaped into a block-column matrix ${\mathbf{c}}_{\mathbf{a}\mathbf{b}}$, such that

$${\mathbf{c}}_{\mathbf{a}\mathbf{b}}=\left[\begin{array}{cc}{\text{a}}_{\text{1,1}}{\text{b}}_{0}& {\text{a}}_{\text{2,1}}{\text{b}}_{0}\\ \vdots & \vdots \\ \begin{array}{c}\vdots \\ {\text{a}}_{\text{1,1}}{\text{b}}_{\text{L}-1}\end{array}& \begin{array}{c}\vdots \\ {\text{a}}_{\text{2,1}}{\text{b}}_{\text{L}-1}\end{array}\end{array} \begin{array}{cc}\cdots & {\text{a}}_{\text{Q},1}{\text{b}}_{0}\\ \ddots & \vdots \\ \begin{array}{c}\ddots \\ \cdots \end{array}& \begin{array}{c}\vdots \\ {\text{a}}_{\text{Q},1}{\text{b}}_{\text{L}-1}\end{array}\end{array} \begin{array}{cc}\cdots \cdots & {\text{a}}_{1,\text{P}}{\text{b}}_{0}\\ \vdots & \vdots \\ \begin{array}{c}\vdots \\ \cdots \cdots \end{array}& \begin{array}{c}\vdots \\ {\text{a}}_{1,\text{P}}{\text{b}}_{\text{L}-1}\end{array}\end{array} \begin{array}{cc}\cdots & {\text{a}}_{\text{Q},\text{P}}{\text{b}}_{0}\\ \ddots & \vdots \\ \begin{array}{c}\ddots \\ \cdots \end{array}& \begin{array}{c}\vdots \\ {\text{a}}_{\text{Q},\text{P}}{\text{b}}_{\text{L}-1}\end{array}\end{array}\right]=\mathbf{b}{\mathbf{a}}^{\text{T}}$$

(13)

where $\mathbf{a}={\left[{\text{a}}_{\text{1,1}},\dots ,{\text{a}}_{\text{Q},1},{\text{a}}_{\text{1,2}},\dots \dots , {\text{a}}_{\text{Q},\text{P}-1}, {\text{a}}_{1,\text{P}},\dots ,{\text{a}}_{\text{Q},\text{P}}\right]}^{\text{T}}$ and $\mathbf{b}={\left[{\text{b}}_{0},\dots ,{\text{b}}_{\text{L}-1}\right]}^{\text{T}}$. An optimal, in the least-squares sense, estimate of the model parameters ${\widehat{\mathbf{a}}}_{\mathbf{L}\mathbf{S}\mathbf{E}}$ and ${\widehat{\mathbf{b}}}_{\mathbf{L}\mathbf{S}\mathbf{E}}$ can be obtained solving the following constrained minimization problem

$$\begin{array}{c}\left({\widehat{\mathbf{a}}}_{\mathbf{LSE}},{\widehat{\mathbf{b}}}_{\mathbf{LSE}}\right)=\underset{{\mathbf{a}}^{\mathbf{^{\prime}}},{\mathbf{b}}^{\mathbf{^{\prime}}}}{\text{argmin}}\left\{{\Vert {\widehat{\mathbf{C}}}_{\mathbf{PLS}}-{\mathbf{b}}{\mathbf{a}}^{\text{T}}\Vert }_{2}^{2}\right\}\\ {\text{s.t.}} {\Vert {\mathbf{a}}\Vert }_{2}=1, {\mathop{\arg\max}\limits_{\text{m}}}\left\{\left|{\text{h}}\left({\text{m}}\right)\right|\right\}>0\end{array}$$

(14)

where $\text{h}\left(\text{m}\right)$ is given by (6). Note that as a result of normalizing the polynomial coefficients $\mathbf{a}$, the estimate ${\widehat{\mathbf{a}}}_{\mathbf{L}\mathbf{S}\mathbf{E}}$ reflects a relative rather than absolute contribution of individual EEG bands to the BOLD signal variance. A solution to (A9) is provided by the singular value decomposition (SVD) of matrix ${\mathbf{c}}_{\mathbf{a}\mathbf{b}}$ (Gómez and Baeyens 2004). Specifically,

$$\begin{array}{c}{\widehat{\mathbf{a}}}_{\mathbf{L}\mathbf{S}\mathbf{E}}=\\ {\widehat{\mathbf{b}}}_{\mathbf{L}\mathbf{S}\mathbf{E}}=\end{array}\begin{array}{c}{\mathbf{U}}_{1}\\ {{\mathbf{V}}_{1}\bullet{\varvec{\Sigma}}}_{1}\end{array}$$

(15)

where ${\mathbf{U}}_{1}$ is the first left singular vector, ${\mathbf{V}}_{1}$ the first right singular vector, and ${{\varvec{\Sigma}}}_{1}\in {\mathbb{R}}$ the first singular value of the SVD of ${\widehat{\mathbf{c}}}_{\text{ab}}$.

Hammerstein-Wiener Model Identification

The MISO Hammerstein-Wiener (HW) described by Eq. (3) can be re-expressed in a compact matrix form as

$$\mathbf{y}=\mathbf{F}\mathbf{z}+{\varvec{\upvarepsilon}}$$

(16)

where $\mathbf{F}$ denotes a matrix the columns of which are polynomial powers of ${\text{y}}_{\text{H}}$, and $\mathbf{z}$ is a vector of the unknown polynomial coefficients, which can be estimated using ordinary least-squares

$${\widehat{\mathbf{z}}}_{\text{LSE}}={\left[{\mathbf{F}}^{\mathbf{T}}\mathbf{F}\right]}^{-1}{\mathbf{F}}^{\mathbf{T}}\mathbf{y}.$$

(17)

Orthonormal Basis Functions

There are several sets of orthonormal basis functions that can be used for modeling the impulse response function of the LTI block in the Hammerstein and Wiener-Hammerstein model configuration (Heuberger et al. 2005). The selection of the appropriate basis set depends on the dynamic behavior of the system to be modelled. One basis set that has been extensively used in the literature for modeling of physiological systems is the Laguerre basis. Laguerre basis functions exhibit exponentially decaying structure and constitute an orthonormal set in $\left[0,\infty \right)$, which makes them suitable for modeling causal systems with finite memory (Marmarelis 1993).

In this work we employ a smoother variant of the Laguerre basis functions that is based on the spherical Laguerre basis functions (Leistedt and McEwen 2012), which allow obtaining robust HRF estimates in single voxels even during resting conditions where the signal-to-noise ratio (SNR) is particularly low. The $\text{j}$-th spherical Laguerre basis function ${\text{b}}_{\text{j}}\left(\text{n}\right); \text{j}=0,\dots ,\text{L}-1;\text{ n}=1,\dots ,\text{M}$ is given by

$${\text{b}}_{\text{j}}\left(\text{n}\right)=\sqrt{\frac{\text{j}!}{\left(\text{j}+2\right)!}}\frac{{\text{e}}^{\frac{\text{n}}{2\alpha}}}{\sqrt{\alpha^{3}}}{\cdot \text{K}}_{\text{j}}\left(\frac{\text{n}}\alpha\right)$$

(18)

where $\alpha\in {\mathbb{R}}_{+}$ is a parameter that determines the rate of exponential asymptotic decline of ${\text{b}}_{\text{j}}\left(\text{n}\right)$, and ${\text{K}}_{\text{j}}\left(\text{n}\right)$ is the $\text{j}$-th generalized Laguerre polynomial of order two, defined as

$${\text{K}}_{\text{j}}\left(\text{n}\right)=\sum_{\text{r}=0}^{\text{j}}\left(\begin{array}{c}\text{j}+2\\ \text{j}-\text{r}\end{array}\right)\frac{{\left(-\text{n}\right)}^{\text{r}}}{\text{r}!}.$$

(19)

The $\text{j}$-th spherical Laguerre basis function ${\text{b}}_{\text{j}}\left(\text{n}\right)$ were convolved with a Gaussian kernel G(μ,τ), with μ > 0 and τ = 1. The μ > 0 parameter controls how late ${\text{b}}_{\text{j}}\left(\text{n}\right)$ will start to fluctuate. The range of this parameter was set to 2.5 < μ < 5 corresponding to HRF estimates with time-to-peak ranging between 3 and 10 s (Fig. S12). The orthonormal properties of the derived basis set was ensured by entering the result of this convolution into an orthogonalization process based on the Gram-Schmidt orthogonalization algorithm.