Identification and comparison of stochastic metabolic/hemodynamic models (sMHM) for the generation of the BOLD signal

Abstract

This paper extends a previously formulated deterministic metabolic/hemodynamic model for the generation of blood oxygenated level dependent (BOLD) responses to include both physiological and observation stochastic components (sMHM). This adds a degree of flexibility when fitting the model to actual data by accounting for un-modelled activity. We then show how the innovation method can be used to estimate unobserved metabolic/hemodynamic as well as vascular variables of the sMHM, from simulated and actual BOLD data. The proposed estimation method allowed for doing model comparison by calculating the model’s AIC and BIC. This methodology was then used to select between different neurovascular coupling assumptions underlying sMHM. The proposed framework was first validated on simulations and then applied to BOLD data from a motor task experiment. The models under comparison in the analysis of the actual data considered that vascular response was coupled to: (I) inhibition, (II) excitation, (III) both excitation and inhibition. Data was best described by model II, although model III was also supported.

Appendices

Appendix A: Expression for the drift

The drift term of system (8) is given by:

$${\mathbf{a}} = {\left( {\begin{array}{*{20}c} {{y_{2} {\left( t \right)}}} \\ {{\theta _{1} {\left( {u_{e} {\left( t \right)} - 1} \right)} - 2\theta _{2} y_{2} {\left( t \right)} - \theta ^{2}_{2} {\left( {y_{1} {\left( t \right)} - 1} \right)}}} \\ {{y_{4} {\left( t \right)}}} \\ {{\theta _{3} {\left( {u_{i} {\left( t \right)} - 1} \right)} - 2\theta _{4} y_{4} {\left( t \right)} - \theta ^{2}_{4} {\left( {y_{3} {\left( t \right)} - 1} \right)}}} \\ {{y_{6} {\left( t \right)}}} \\ {{\theta _{8} {\left( {u_{e} {\left( t \right)} - 1} \right)} + \theta _{{26}} {\left( {u_{i} {\left( t \right)} - 1} \right)} - 2\theta _{9} y_{6} {\left( t \right)} - \theta ^{2}_{9} {\left( {y_{5} {\left( t \right)} - 1} \right)}}} \\ {{\frac{1}{{\theta _{{10}} {\left( {\theta _{7} + 1} \right)}}}{\left( {y_{3} {\left( t \right)} + \theta _{7} \frac{{\lambda - \frac{1}{{1 - e^{{ - \theta _{5} {\left( {y_{1} {\left( t \right)} - \theta _{6} } \right)}}} }}}}{{\lambda - \frac{1}{{1 - e^{{ - \theta _{5} {\left( {1 - \theta _{6} } \right)}}} }}}}y_{1} {\left( t \right)}} \right)} - \frac{{{\left( {\theta _{{10}} {\left( {y_{8} {\left( t \right)}} \right)}^{{\frac{1}{{\theta _{{12}} }}}} + \theta _{{11}} y_{5} {\left( t \right)}} \right)}y_{7} {\left( t \right)}}}{{\theta _{{10}} {\left( {\theta _{{10}} + \theta _{{11}} } \right)}y_{8} {\left( t \right)}}}}} \\ {{\frac{1}{{\theta _{{10}} + \theta _{{11}} }}{\left( {y_{5} {\left( t \right)} - {\left( {y_{8} {\left( t \right)}} \right)}^{{\frac{1}{{\theta _{{12}} }}}} } \right)}}} \\ \end{array} } \right)}$$

(19)

where $\lambda = 2 - \frac{1}{{1 - e^{\theta _5 \theta _6 } }}$

Appendix B: Local linearization filter algorithm

Let a state-space model be defined by the continuous state equation:

$$d{\mathbf{y}}\left( t \right) = {\mathbf{a}}\left( {t,{\mathbf{y}}\left( t \right),{\mathbf{\theta }}} \right)dt + {\mathbf{b}}\left( {\mathbf{\theta }} \right)d{\mathbf{w}}\left( t \right){\text{, for }}t \geqslant t_0 $$

(20)

and the discrete observation equation:

$$z_{t_k } = {\mathbf{c}}\left( {\mathbf{\theta }} \right)^T {\mathbf{y}}\left( {t_k } \right) + {\mathbf{e}}_{t_k } ,\quad {\text{for }}k = 0,1,....,N{\text{.}}$$

(21)

Assuming that θ is given and starting with $\widehat{\mathbf{y}}_{{{t_0 } \mathord{\left/ {\vphantom {{t_0 } {t_0 }}} \right. \kern-\nulldelimiterspace} {t_0 }}} = {\mathbf{y}}_{{{t_0 } \mathord{\left/ {\vphantom {{t_0 } {t_0 }}} \right. \kern-\nulldelimiterspace} {t_0 }}}$, $\widehat{{\mathbf{P}}}_{{{t_{0} } \mathord{\left/ {\vphantom {{t_{0} } {t_{0} }}} \right. \kern-\nulldelimiterspace} {t_{0} }}} = {\mathbf{P}}_{{{t_{0} } \mathord{\left/ {\vphantom {{t_{0} } {t_{0} }}} \right. \kern-\nulldelimiterspace} {t_{0} }}}$, the LL filter for the model (20)–(21) is defined by the recursive computation of the following estimates (Jiménez and Ozaki 2003; Riera et al. 2007):

For all k = 0,...,M − 1.

• 1. Prediction: For j = 0,…,M − 1,(M ≥ 1):

  $$s_j = t_k + j\frac{{t_{k + 1} - t_k }}{M}$$

  (22)
  $$\widehat{\mathbf{y}}_{s_{j + 1} /t_k } = \widehat{\mathbf{y}}_{s_j /t_k } + \int\limits_0^{s_{j + 1} - s_j } {e^{{\mathbf{D}}_{s_j } \left( {s_{j + 1} - s_j - s} \right)} } \left( {{\mathbf{a}}\left( {s_j ,\widehat{\mathbf{y}}_{s_j /t_k } ,{\mathbf{\theta }}} \right) + {\mathbf{d}}_{s_j } s} \right)ds$$

  (23)
  $$\widehat{\mathbf{P}}_{s_{j + 1} /t_k } = e^{{\mathbf{D}}_{s_j } \left( {s_{j + 1} - s_j - s} \right)} \widehat{\mathbf{P}}_{s_j /t_k } e^{{\mathbf{D}}_{_{s_j } }^T \left( {s_{j + 1} - s_j - s} \right)} + \int\limits_0^{s_{j + 1} - s_j } {e^{{\mathbf{D}}_{s_j } s} } {\mathbf{bb}}^{\text{T}} e^{{\mathbf{D}}_{{\text{s}}_{\text{j}} }^{\text{T}} s} ds$$

  (24)

• 2. Innovation:

  $$\widehat\vartheta _{t_{k + 1} } = z_{t_{k + 1} } - {\mathbf{c}}^T \widehat{\mathbf{y}}_{{{t_{k + 1} } \mathord{\left/{\vphantom {{t_{k + 1} } {t_k }}} \right.\kern-\nulldelimiterspace} {t_k }}} $$

  (25)
  $$\widehat{\text{ $ \Sigma $ }}_{t_{k + 1} /t_k }^\nu = {\mathbf{c}}^T \widehat{\mathbf{P}}_{{{t_{k + 1} } \mathord{\left/{\vphantom {{t_{k + 1} } {t_k }}} \right.\kern-\nulldelimiterspace} {t_k }}} {\mathbf{c}} + \theta _8 $$

  (26)

• 3. Filter:

  $$\widehat{\mathbf{y}}_{{{t_{k + 1} } \mathord{\left/{\vphantom {{t_{k + 1} } {t_{k + 1} }}} \right.\kern-\nulldelimiterspace} {t_{k + 1} }}} = \widehat{\mathbf{y}}_{{{t_{k + 1} } \mathord{\left/{\vphantom {{t_{k + 1} } {t_k }}} \right.\kern-\nulldelimiterspace} {t_k }}} + {\mathbf{k}}_{t_{k + 1} } \widehat\vartheta _{t_{k + 1} } $$

  (27)
  $$\widehat{\mathbf{P}}_{{{t_{k + 1} } \mathord{\left/{\vphantom {{t_{k + 1} } {t_{k + 1} }}} \right.\kern-\nulldelimiterspace} {t_{k + 1} }}} = \widehat{\mathbf{P}}_{{{t_{k + 1} } \mathord{\left/{\vphantom {{t_{k + 1} } {t_k }}} \right.\kern-\nulldelimiterspace} {t_k }}} - {\mathbf{k}}_{t_{k + 1} } {\mathbf{c}}^T \widehat{\mathbf{P}}_{{{t_{k + 1} } \mathord{\left/{\vphantom {{t_{k + 1} } {t_k }}} \right.\kern-\nulldelimiterspace} {t_k }}} $$

  (28)

  where ${\mathbf{k}}_{t_{k + 1} } = {\mathbf{\hat P}}_{t_{k + 1} /t_k } {\mathbf{c}}^{\text{T}} \left( {\hat \Sigma _{t_{k + 1} }^\vartheta } \right)^{ - 1}$ is the filter gain. In the above algorithm M = 1 denotes the number of points between a pair of consecutive observations t _k ,t _{k  + 1} used to compute the predictions $\widehat{\mathbf{y}}_{{{t_{k + 1} } \mathord{\left/ {\vphantom {{t_{k + 1} } {t_k }}} \right. \kern-\nulldelimiterspace} {t_k }}}$ and $\widehat{\mathbf{P}}_{{{t_{k + 1} } \mathord{\left/ {\vphantom {{t_{k + 1} } {t_k }}} \right. \kern-\nulldelimiterspace} {t_k }}}$, ${\mathbf{D}}_{s_j } = \frac{\partial }{{\partial {\mathbf{y}}}}{\mathbf{a}}\left( {s_j ,{\mathbf{y}}_{s_j /t_k } ;{\mathbf{\theta }}} \right)$ and ${\mathbf{d}}_{s_j } = \frac{\partial }{{\partial s}}{\mathbf{a}}\left( {s_j ,{\mathbf{y}}_{s_j /t_k } ,{\mathbf{\theta }}} \right)$. Explicit expressions in terms of exponential matrices for integrals (23) and (24) are given in Jiménez and Ozaki, (2003). Thus, the numerical computation of the conditional moments ${\mathbf{\hat y}}_{t_{k + 1} /t_k }$ and ${\mathbf{\hat P}}_{{{t_{k + 1} } \mathord{\left/{\vphantom {{t_{k + 1} } {t_k }}} \right. \kern-\nulldelimiterspace} {t_k }}}$ is reduced to use a convenient algorithm to compute matrix exponentials. This was achieved in the present paper by using Krylov subspace methods (Hochbruck and Lubich 1997).