(Task-based) Functional Magnetic Resonance Imaging?

• fMRI is a noninvasive neuroimaging method that measures the blood-oxygen level-dependent (BOLD) signals in the brain.
  • Large-dimensional
  • Low signal-to-noise ratios (SNRs)
  • Sophisticated spatio-temporal dependence of voxels.

Why Need Complex-valued Models

• Raw fMRI data are complex-valued (CV) after Fourier transform and inverse Fourier transform image reconstruction.
• Most fMRI studies use magnitude-only (MO) data and phase information is discarded.

Complex-valued fMRI (CV-fMRI) Data

Why Bayesian Models of CV-fMRI

• Sophisticated Bayesian spatiotemporal models have been developed for MO-fMRI (Zhang Guindani, et al., 2016).

... the most common software packages for fMRI analysis can result in false-positive rates of up to 70% --- Eklund et al. (PNAS 2016)

• Start with the most common objective: detecting brain activations in task-based fMRI.

Goal: Computationally Effcient Models for CV-fMRI

Propose

• computationally efficient models/algorithms that
  • jointly analyze real and imaginary parts of CV-fMRI data
  • better detect activation at the voxel level
  • infer brain connectivity (in preparation)

• Complex-valued Expectation-Maximization Variable Selection with Autoregressive Processes (CV-EMVS-AR)
• (Multi-subject) Bayesian Spatiotemporal Model via Kernel Convolution and Autoregressive Processes (CV-KC-AR)

Background: Rowe-Logan Constant Phase Model

From Rowe and Logan (2004), for time $t=1:T$, voxel $v=1:V$ and $p$ tasks, $$\begin{array}{cc}{y}_t^v& ={\rho }_t^v\cos \left({\varphi }^v\right)+i{\rho }_t^v\sin \left({\varphi }^v\right)+{\eta }_t^v,\\ {\rho }_t^v& ={\beta }_0^v+{\beta }_1^v{x}_{1,t}+\cdots +{\beta }_p^v{x}_{p,t}\end{array}$$

$$\left[\begin{array}{c}{y}_1^v\\ ⋮\\ y_T^v\end{array}\right]=\left[\begin{array}{cccc}1& x_{11}& \cdots & x_{p1}\\ ⋮& ⋮& \ddots & ⋮\\ 1& x_{1T}& \cdots & x_{pT}\end{array}\right]\left[\begin{array}{c}{\beta }_0^v\cos \left({\varphi }^v\right)+i{\beta }_0^v\sin \left({\varphi }^v\right)\\ ⋮\\ \underset{{\gamma }_{Re}^v}{\underset{}{{\beta }_p^v\cos \left({\varphi }^v\right)}}+i\underset{{\gamma }_{Im}^v}{\underset{}{{\beta }_p^v\sin \left({\varphi }^v\right)}}\end{array}\right]+\left[\begin{array}{c}{\eta }_1^v\\ ⋮\\ {\eta }_T^v\end{array}\right]$$

• $y^v=X{\gamma }_{Re}^v+iX{\gamma }_{Im}^v+{\eta }^v$
• $y^v=X{\gamma }^v+{\eta }^v$, ${\eta }^v\sim C{N}_T(0,{\Gamma }_v,C_v)$. ${\Gamma }_v\in C^{T\times T}$, $C_v\in C^{T\times T}$
• ${\eta }^v$ is circular normal if $C_v=0$

Background: Real-valued Representation

$y^v=X{\gamma }^v+{\eta }^v$, ${\eta }^v\sim C{N}_T(0,{\Gamma }_v,C_v)$

or equivalently, $y_r^v=X^r{\gamma }_r^v+{\eta }_r^v$, where ${\eta }_r^v\sim N_{2T}(0,{\Sigma }_v)$

Brain Activation as Variable Selection

• $y^v=X{\gamma }^v+{\eta }^v$

• ${\gamma }_j^v\ne 0$ iff voxel $v$ at task $j$ is activated (Xia, Liang, and Wang (2009); Zhang, Guindani, and Vannucci (2015)).

• Complex normal spike-and-slab prior $${\gamma }_j^v\sim (1-{\psi }_j^v)\underset{spike}{\underset{}{CN(0,{\omega }_0,{\lambda }_0)}}+{\psi }_j^v\underset{slab}{\underset{}{CN(0,{\omega }_1,{\lambda }_1)}},$$ ${\omega }_0<{\omega }_1\in C,{\lambda }_0<{\lambda }_1\in C$.

Brain Activation as Variable Selection

• Complex normal spike-and-slab prior $${\gamma }_j^v\sim (1-{\psi }_j^v)\underset{spike}{\underset{}{CN(0,{\omega }_0,{\lambda }_0)}}+{\psi }_j^v\underset{slab}{\underset{}{CN(0,{\omega }_1,{\lambda }_1)}},$$ ${\omega }_0<{\omega }_1\in C,{\lambda }_0<{\lambda }_1\in C$.
• Non-active voxel: ${\psi }_j^v=0\Rightarrow {\gamma }_j^v=0$

• Active voxel: ${\psi }_j^v=1\Rightarrow {\gamma }_j^v\ne 0$

• Activation is inferred by borrowing information across voxels through a Bernoulli prior on ${\psi }_j^v$ with a common probability of activation for all voxels: ${\psi }_j^v\sim Bernoulli({\theta }_j^v={\theta }_j)$, i.e., $Pr({\psi }_j^v=1\mid {\theta }_j)={\theta }_j$.

CV-EMVS (Yu, Prado, Ombao, and Rowe, 2018)

$$\begin{array}{cc}{y}^v& =X{\gamma }^v+{\eta }^v,{\eta }^v\sim C{N}_T(0,2{\sigma }_v^{2}I,0),v=1:V,\\ {\gamma }_j^v\mid {\psi }_j^v& \stackrel{indep}{\sim }(1-{\psi }_j^v)C{N}_1(0,{\sigma }_v^2{\omega }_0,{\sigma }_v^2{\lambda }_0)+{\gamma }_j^{v}C{N}_1(0,{\sigma }_v^2{\omega }_1,{\sigma }_v^2{\lambda }_1),\\ {\omega }_0& <<{\omega }_1,{\lambda }_0<<{\lambda }_1,j=1:p,\\ {\psi }_j^v\mid {\theta }_j& \stackrel{IID}{\sim }Bernoulli\left({\theta }_j\right),{\theta }_j\stackrel{IID}{\sim }Beta(a_{\theta },b_{\theta }),{\sigma }_v^2\sim IG(a_{\sigma },b_{\sigma })\end{array}$$

• We determine if the real and imaginary parts of ${\gamma }_j^v$ are zero jointly: ${\gamma }_j^v\ne 0$ if $Pr({\psi }_j^v=1\mid {\gamma }^{\ast },{\theta }^{\ast },{\sigma }^{\ast },y)>\delta$.
  • E-step: derive

    $$Pr({\psi }_j^v=1\mid {\gamma }^{\left(l\right)},{\theta }^{\left(l\right)},{\sigma }^{\left(l\right)},y)$$

  • M-step: obtain maximum a posteriori

    $${\gamma }^{(l+1)},{\theta }^{(l+1)},{\sigma }^{(l+1)}=max\arg E_{\psi |\cdot }[\log \pi (\gamma ,\psi ,\theta ,\sigma \mid y)\mid {\gamma }^{\left(l\right)},{\theta }^{\left(l\right)},{\sigma }^{\left(l\right)},y]$$

CV-EMVS (Yu, Prado, Ombao, et al., 2018)

$$\begin{array}{cc}{y}^v& =X{\gamma }^v+{\eta }^v,{\eta }^v\sim C{N}_T(0,2{\sigma }_v^{2}I,0),v=1:V,\\ {\gamma }_j^v\mid {\psi }_j^v& \stackrel{indep}{\sim }(1-{\psi }_j^v)C{N}_1(0,{\sigma }_v^2{\omega }_0,{\sigma }_v^2{\lambda }_0)+{\gamma }_j^{v}C{N}_1(0,{\sigma }_v^2{\omega }_1,{\sigma }_v^2{\lambda }_1),\\ {\omega }_0& <<{\omega }_1,{\lambda }_0<<{\lambda }_1,j=1:p,\\ {\psi }_j^v\mid {\theta }_j& \stackrel{IID}{\sim }Bernoulli\left({\theta }_j\right),{\theta }_j\stackrel{IID}{\sim }Beta(a_{\theta },b_{\theta }),{\sigma }_v^2\sim IG(a_{\sigma },b_{\sigma })\end{array}$$

• Under circular normal prior on ${\gamma }_j^v$, i.e., ${\lambda }_0={\lambda }_1=0$,

  $$\begin{array}{cc}{\hat{\gamma }}_{Re}^v& =(X^{\prime }X+2D_v{)}^{-1}(X^{\prime }{y}^v{)}_{Re},\\ {\hat{\gamma }}_{Im}^v& =(X^{\prime }X+2D_v{)}^{-1}(X^{\prime }{y}^v{)}_{Im}.\end{array}$$

Activation and Strength Maps: Low SNR

• The CV model outperforms other MO models in terms of activation detection and strength.
• CV: CV-EMVS
• MO: MO-EMVS
• ALA: Adaptive Lasso

Human CV-fMRI (Karaman, Bruce, and Rowe, 2015)

• Unilateral finger tapping data of dimension 96 x 96 x 510.
• KBR-CV: false positives outside the brain area.
• KBR-MO: low detecting power.
• DeTeCT-ING: Nonlinear model using the MR signal equation.

CV-EMVS with AR Noise

$$\begin{array}{ccc}{y}_t^v& =& {\gamma }_1^v+{\gamma }_2^{v}t/T+{\gamma }_3^v{x}_t+{\eta }_t^v,\\ {\eta }_t^v& =& {\phi }_v{\eta }_{t-1}^v+{\zeta }_t^v,{\zeta }_t^v\stackrel{iid}{\sim }C{N}_1(0,2{\sigma }_v^2,0),\\ {\phi }_v& \sim & Uniform(-1,1).\end{array}$$

General Bayesian spatio-temporal model

• CV-EMVS-AR much improves detection but does not explicitly model the spatial dependence across voxels.
• The execution of biological tasks involves populations of neurons spanning across several voxels, rather than a single voxel.
• Propose Bayesian variable selection tools coupled with a spatial kernel convolution (KC) structure as well as temporal autoregressive processes for CV-fMRI data. (Yu, Prado, Ombao, and Rowe, 2022)

Kernel Convolution: Definition

• Given a kernel $k(z;\varphi ),z\in S$ and a white noise process $w\left(u\right),u\in S$, the kernel convolution is $$S\left(z\right)={\int }_{S}k(z-u;\varphi )w\left(u\right)du.$$

• In practice, for sites $u_1,...,u_D$, the process is defined as $$S\left(z\right)=\sum _{d=1}^{D}k(z-u_d;\varphi )w\left(u_d\right)$$

• The spatial process $S\left(z\right)$ governs the spatial dependence of voxels in the image, and affects the probability of voxels being activated.

Kernel Convolution: Properties

• Dimension reduction: A small number $D$ of parameters $w\left(u_1\right),...,w\left(u_D\right)$ governs the entire process $S\left(z\right)$ that may have tons of measurements at location $z_1,...,z_V$. $(V>>D)$.

• $S=\left(S\right(z_1),...,S(z_V){)}^{\prime }$ is the voxel-level spatial effects

• $w=\left(w\right(u_1),...,w(u_D){)}^{\prime }$ is the $D$-dimensional latent spatial effect formed by the selected $D$ sites.

Kernel convolution (KC) vs. Gaussian process (GP)

Advantages of KC over GP

Complex-valued Bayesian Spatiotemporal Model

• Likelihood: With the indicators ${\psi }_j^v$ such that ${\gamma }_j^v\ne 0$ if ${\psi }_j^v=1$ and ${\gamma }_j^v=0$ if ${\psi }_j^v=0$, $$\begin{array}{cc}{y}^v& =X\left({\psi }^v\right){\gamma }^v\left({\psi }^v\right)+{\eta }^v\\ {\eta }^v& \sim C{N}_T(0,2{\sigma }_v^2{\Lambda }_v,0),\end{array}$$

where ${\Lambda }_v$ is the AR(1) correlation matrix.

• Empirical Bayes estimator ${\hat{\rho }}_v$ $\left({\hat{\Lambda }}_v\right)$ for AR coefficient for computation efficiency.

• Complex-valued g-prior on ${\gamma }^v$: $$\begin{array}{cc}{\gamma }^v\left({\psi }^v\right)\mid {\psi }^v,{\sigma }_v^2& \stackrel{ind}{\sim }C{N}_p\left({\hat{\gamma }}^v\right({\psi }^v),2T{\sigma }_v^2(X^{\prime }\left({\psi }^v\right){\hat{\Lambda }}_v^{-1}X\left({\psi }^v\right){)}^{-1},0)\\ {\hat{\gamma }}^v\left({\psi }^v\right)& =\left(X^{\prime }\right({\psi }^v\left){\hat{\Lambda }}_v^{-1}X\right({\psi }^v\left){)}^{-1}{X}^{\prime }\right({\psi }^v)y^v\end{array}$$

• Integrate ${\gamma }^v$ out for faster computation.

Spatial Priors

Choice of Kernels

• Gaussian kernel $k(z_v-u_d;\varphi )=\exp (-\frac{\parallel z_v-u_d{\parallel }^2}{2\varphi })$
• Bezier kernel $$k(s_v-u_d;\nu ,\varphi )=\{\begin{array}{cc}{(1-\frac{\parallel z_v-u_d{\parallel }^2}{{\varphi }^2})}^{\nu },& \parallel z_v-u_d\parallel <\varphi \\ 0,& \text{otherwise}\end{array},$$ where $\nu$ is the smooth parameter and $\varphi$ is the range parameter.

• The Bezier kernel has a compact support.
• To capture neighboring spatial dependence
  • this avoids unrealistically relating any two voxels together
  • lets the model learn the neighboring structure of a given voxel by learning $\varphi$

Bezier Kernel: $\nu =2$; $\varphi =2,4,6$

• Here shows you what the Bezier kernel looks like with different value of parameters.

Bezier Kernel: $\nu =0.5,2,5$; $\varphi =3$

Simulated data with AR coefficient 0.5

• Estimate $Pr({\psi }_j^v=1|y)$ by computing the number of 1s of ${\psi }_j^v$ in the posterior sample divided by the number of MCMC iterations.

Simulated data with AR coefficient 0.5

• Estimate $Pr({\psi }_j^v=1|y)$ by computing the number of 1s of ${\psi }_j^v$ in the posterior sample divided by the number of MCMC iterations.

Simulated data with AR coefficient 0

KC is Less Sensitive to Dimension Reduction

Computing Time

• Average computing time over 10 times of running the MCMC algorithms.
• The time unit is seconds per 1000 MCMC iterations

Model	$16$	$25$	$100$
CV-KC	0.51	0.59	1.48
CV-GP	0.30	0.41	3.36

• The example shows the KC model can use just about 15% to 20% of the computing time of the GP model per 1000 MCMC iterations to reach the similar detecting activation performance

CV vs. MO

• The CV model is better than the MO model especially when signals are noisy.
• The MO models improve more when an explicit spatial prior is included.

MO-GP Has a Region Probability Inflation

If a spatial region contains true activated voxels, other nonactivated voxels in the region will also have relatively high probability of activation $⟹$ increases false positives.

Spatial Models Encourage Activation in Clusters

Zoom In

Multi-resolution

• $D$ coarse-resolution sites (white dots) and $H$ fine-resolution sites (green dots)

$$S^v=\sum _{d=1}^D{k}_c(z_v-s_d;{\varphi }_c)w^d+\sum _{h=1}^H{k}_f(z_v-s_d;{\varphi }_f)b^h$$

Take-home Message

• Complex-valued modeling improves activation.

• CV-EMVS and CV-KC-AR are computationally efficient.

• Spatial models encourage activation in clusters.

• The MO models need a sophisticated spatial structure to reach good performance as CV models.

• The KC models are flexible and outperform the baseline GP models.

  • Any valid kernel can be used.
  • Avoid dealing with problems of shape of regions.
  • How a voxel is influenced by its neighboring voxels is model-based.
  • Less affected by dimension reduction, leading to fast computation.

References

[1] M. Karaman, I. P. Bruce, and D. B. Rowe. "Incorporating relaxivities to more accurately reconstruct MR images". In: Magnetic Resonance Imaging 33 (2015), pp. 374-384.

[2] D. B. Rowe and B. R. Logan. "A complex way to compute fMRI activation". In: NeuroImage 23 (2004), pp. 1078-1092.

[3] J. Xia, F. Liang, and Y. Wang. "FMRI analysis through Bayesian variable selection with a spatial prior". In: Proceedings of the 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro. Boston, MA, USA, 2009, pp. 714-717.

[4] C. Yu, R. Prado, H. Ombao, et al. "A Bayesian variable selection approach yields improved detection of brain activation from complex-valued fMRI". In: Journal of American Statistical Association: Applicaiton and Case Studies 113 (2018), pp. 1395-1410.

[5] C. Yu, R. Prado, H. Ombao, et al. "Bayesian spatiotemporal modeling on complex-valued fMRI signals via kernel convolutions". In: Biometrics (2022), pp. 1-13.

[6] L. Zhang, M. Guindani, and M. Vannucci. "Bayesian models for fMRI data analysis". In: WIREs Computational Statistics 7 (2015), pp. 21-41.

[7] L. Zhang, M. Guindani, F. Versace, et al. "A spatio-temporal non-parametric Bayesian model of multi-subject fMRI data". In: Annals of Applied Statistics (2016).