FAST Adaptive Smoothing and Thresholding for Improved Activation Detection in Low-Signal fMRI

TL;DR

This paper introduces FAST, an automated adaptive smoothing and thresholding method that enhances activation detection in low-signal fMRI data, improving reliability in identifying brain regions during cognitive tasks.

Contribution

The paper presents a novel FAST algorithm combining smoothing and extreme value theory for better thresholding in low-signal fMRI analysis.

Findings

Effective in low-signal scenarios

Performs well in identifying sensory-specific brain regions

Shows promising results across diverse experiments

Abstract

Functional Magnetic Resonance Imaging is a noninvasive tool for studying cerebral function. Many factors challenge activation detection, especially in low-signal scenarios that arise in the performance of high-level cognitive tasks. We provide a fully automated fast adaptive smoothing and thresholding (FAST) algorithm that uses smoothing and extreme value theory on correlated statistical parametric maps for thresholding. Performance on experiments spanning a range of low-signal settings is very encouraging. The methodology also performs well in a study to identify the cerebral regions that perceive only-auditory-reliable or only-visual-reliable speech stimuli.

Tables

Table 1. TABLE I: ϕ bold-italic-ϕ \boldsymbol{\phi} s for the AR( p 𝑝 p ) scenarios used in our simulations.

$p$	Decreasing	Increasing
2	(0.6, 0.3)	(0.3, 0.6)
3	(0.4, 0.3, 0.2)	(0.2, 0.3, 0.4)
4	(0.3, 0.25, 0.20, 0.15)	(0.15, 0.20, 0.25, 0.3)
5	(0.3, 0.25, 0.20, 0.10, 0.05)	(0.05, 0.10, 0.20, 0.25, 0.30)
$p$	Decreasing-Increasing	Increasing-Decreasing
2	(0.6,0.3)	(0.3,0.6)
3	(0.4,0.1,0.4)	(0.1,0.7,0.1)
4	(0.4,0.05,0.05,0.4)	(0.05,0.4,0.4,0.05)
5	(0.25,0.15,0.1,0.15,0.25)	(0.1,0.15,0.4,0.15,0.1)

Table 2. TABLE S1: Number of SPMs (out of 1000) simulated under null conditions where we found even one activated pixel using (a) ALL-FAST, (b) AM-FAST and (c) AR-FAST using different values of α 𝛼 \alpha and for different values of ρ 𝜌 \rho (the assumed radially symmetric first order autocorrelations between neigboring voxels in each simulated SPM).

						$\rho$
		0	0.01	0.025	0.05	0.075	0.1	0.25	0.5	0.75	0.99
$\alpha$	0.001	0	0	0	0	0	0	0	0	0	0
	0.01	0	0	0	0	0	0	0	0	0	0
	0.025	0	1	0	0	0	0	0	0	0	0
	0.05	0	2	2	1	0	0	0	0	1	0
	0.075	1	4	4	1	0	0	0	0	2	2
	0.1	3	5	6	1	0	0	0	0	2	3

Table 3. (a)

						$\rho$
		0	0.01	0.025	0.05	0.075	0.1	0.25	0.5	0.75	0.99
$\alpha$	0.001	0	0	0	0	0	0	0	0	0	0
	0.01	0	0	0	0	0	0	0	0	0	0
	0.025	0	1	0	0	0	0	0	0	0	0
	0.05	0	2	2	1	0	0	0	0	1	0
	0.075	1	4	4	1	0	0	0	0	2	2
	0.1	3	5	6	1	0	0	0	0	2	3

Table 4. (b)

						$\rho$
		0	0.01	0.025	0.05	0.075	0.1	0.25	0.5	0.75	0.99
$\alpha$	0.001	0	0	0	0	0	0	0	0	0	0
	0.01	0	0	0	0	0	0	0	0	0	0
	0.025	0	0	0	0	0	0	0	0	0	0
	0.05	0	1	0	0	0	0	0	0	0	0
	0.075	0	1	0	0	0	0	0	0	0	0
	0.1	0	1	0	0	0	0	0	0	0	0

Table 5. (c)

						$\rho$
		0	0.01	0.025	0.05	0.075	0.1	0.25	0.5	0.75	0.99
$\alpha$	0.001	0	0	0	0	0	0	0	0	0	7
	0.01	0	0	0	0	0	0	0	0	0	14
	0.05	0	0	0	0	0	0	0	0	0	15
	0.05	0	0	0	0	0	0	0	0	0	20
	0.075	0	0	0	0	0	0	0	0	0	22
	0.1	0	0	0	0	0	0	0	0	0	22

Equations

$${\mathbf{Y}}_i=\mathbf{X}{\mathbf{\beta }}_i+{\mathbf{ϵ}}_i,$$

$$\begin{array}{rl}{\Phi }_{(n)}(x/\varrho )& \le \mathbb{P}[{\mathbf{R}}^{1/2}\mathbf{Z}\le x{\mathbf{R}}^{1/2}1/\varrho ]\\ & =\mathbb{P}[\mathbf{X}\le x1],\text{mbox}where\mathbf{X}\sim N_n(0,\mathbf{R})\\ & =\mathbb{P}[X_{(n)}\le x]=F_{(n)}(x).\end{array}$$

$$\underset{n\to \infty }{\lim }[F_n(a_{n}x+b_n)]\ge \exp \{-\exp (-x)\},$$

$${\varphi }_{\eta }^{\bullet }(y;\varrho )=\frac{\varphi (\frac{y}{\varrho })}{\varrho \Phi (\frac{\eta }{\varrho })}I(y<\eta );{\Phi }_{\eta }^{\bullet }(y;\varrho )=\frac{\Phi [\frac{\min (y,\eta )}{\varrho }]}{\Phi (\frac{\eta }{\varrho })},$$

$$\underset{n\to \infty }{\lim }[{\Phi }_{\eta ,(n)}^{\bullet }(a_n^{\bullet }x+b_n^{\bullet };\varrho )]=\exp \{-(-x{)}^{-\nu }\}I(x\le 0).$$

$$\underset{y\to \eta }{\lim }\frac{(\eta -y){\varphi }^{\bullet }(y;\varrho )}{1-{\Phi }^{\bullet }(y;\varrho )}=\nu$$

$$\begin{array}{rl}\underset{y\to \eta }{\lim }\frac{(\eta -y){\varphi }_{\eta }^{\bullet }(y;\varrho )}{1-{\Phi }_{\eta }^{\bullet }(y;\varrho )}& =\underset{y\to \eta }{\lim }\frac{(\eta -y)\frac{d}{dx}{\varphi }_{\eta }^{\bullet }(y;\varrho )-{\varphi }_{\eta }^{\bullet }(y;\varrho )}{-{\varphi }_{\eta }^{\bullet }(y;\varrho )}\\ & =\underset{y\to \eta }{\lim }\frac{(\eta -y){\varphi }^{\prime }(\frac{y}{\varrho })/\varrho -\varphi (\frac{y}{\varrho })/\varrho }{-\varphi (\frac{y}{\varrho })/\varrho }=1.\end{array}$$

$$\begin{array}{rl}{\Phi }_{\frac{\eta }{\varrho },(n)}^{\bullet }(x/\varrho ;1)=\mathbb{P}[\mathbf{U}\le x1/\varrho ]& \le \mathbb{P}[{\mathbf{R}}^{1/2}\mathbf{Y}\le x{\mathbf{R}}^{1/2}1/\varrho ]\\ & =\mathbb{P}[\mathbf{X}\le x1]\equiv F_{\eta }^{\bullet }(x),\end{array}$$

$$\underset{n\to \infty }{\lim }[F_{\eta }^{\bullet }(a_n^{\bullet }x+b_n^{\bullet })]\ge \exp \{-(-x{)}^{-\nu }\}I(x\le 0).$$

$$\ell (\mathbf{h}\mid \sigma ,{\mathbf{\Gamma }}_{(-\mathbf{h})})=constant-\frac{1}{2}\log |{\sigma }^2{\mathbf{S}}_{\mathbf{h}}|-\frac{1}{2{\sigma }^2}{\mathbf{\Gamma }}_{(-\mathbf{h})}^{\prime }{\mathbf{\Gamma }}_{(-\mathbf{h})}.$$

$${\eta }_k=\{\begin{array}{ll}{a}_{n_0}{\iota }_{\alpha }^{\mathcal{G}}+b_{n_0}& \text{mbox}fork=1\\ a_{n_{k-1}}^{\bullet }{\iota }_{\alpha }^{\mathcal{W}}+b_{n_{k-1}}^{\bullet }& \text{mbox}otherwise.\end{array}$$

$${\tilde{s}}_{(k,w)}=\frac{{\left(n\right)}^{\frac{1}{2}}{\left[{\sum }_{|e_{i(k),j}|<1}{\Gamma }_{i(k)}^2{\left(1-e_{i(k)}^2\right)}^4\right]}^{\frac{1}{2}}}{|{\sum }_{|e_{i(k)}<1|}\left(1-e_{i(k)}^2\right)\left(1-5e_{i(k)}^2\right)|}$$

$${\mathbf{\Lambda }}_{\mathbf{\beta }}^{-1}={\beta }_1{\mathbf{J}}_{n_1}\otimes {\mathbf{I}}_{n_1}\otimes {\mathbf{I}}_{n_1}+{\beta }_2{\mathbf{I}}_{n_1}\otimes {\mathbf{J}}_{n_2}\otimes {\mathbf{I}}_{n_2}+{\beta }_3{\mathbf{I}}_{n_1}\otimes {\mathbf{I}}_{n_2}\otimes {\mathbf{J}}_{n_3},$$

$$\left(\begin{array}[]{rrrrrr}1&-1&0&0&0&\cdots\\ -1&2&-1&0&0&\cdots\\ 0&-1&2&-1&0&\cdots\\ &&&\ddots&&\\ \cdots&0&0&-1&2&-1\\ \cdots&0&0&0&-1&1\\ \end{array}\right)\\$$

$$\hat{\mathbf{\beta }}=\underset{\mathbf{\beta }}{\mathrm{argmin}}\left[\log |{\mathbf{\Lambda }}_{\mathbf{\beta }}+\mathbf{R}|+{\mathbf{\Gamma }}^{\prime }({\mathbf{\Lambda }}_{\mathbf{\beta }}+\mathbf{R}{)}^{-1}\mathbf{\Gamma }\right].$$

$${\sigma }^2+{\left\{{\beta }_1\left[2-2\cos \left(\frac{(i-1)\pi }{n_1}\right)\right]+{\beta }_2\left[2-2\cos \left(\frac{(j-1)\pi }{n_2}\right)\right]+{\beta }_2\left[2-2\cos \left(\frac{(k-1)\pi }{n_3}\right)\right]\right\}}^{-1}$$

$$\begin{array}{rl}\hat{\mathbf{\beta }}=\underset{\sigma ,\mathbf{\beta }}{\mathrm{argmax}}\sum _{i,j,k}& [\log \left({\sigma }^2+{\left\{{\beta }_1\left[2-2\cos \left(\frac{(i-1)\pi }{n_1}\right)\right]+{\beta }_2\left[2-2\cos \left(\frac{(j-1)\pi }{n_2}\right)\right]+{\beta }_2\left[2-2\cos \left(\frac{(k-1)\pi }{n_3}\right)\right]\right\}}^{-1}\right)\\ & +\frac{w_{ijk}^2}{{\sigma }^2+{\left\{{\beta }_1\left[2-2\cos \left(\frac{(i-1)\pi }{n_1}\right)\right]+{\beta }_2\left[2-2\cos \left(\frac{(j-1)\pi }{n_2}\right)\right]+{\beta }_2\left[2-2\cos \left(\frac{(k-1)\pi }{n_3}\right)\right]\right\}}^{-1}}].\end{array}$$

Full text

I Introduction

Functional Magnetic Resonance Imaging (fMRI) [1, 2, 3, 4, 5, 6, 7, 8, 9] studies the spatial characteristics and extent of brain function while at rest or, more commonly, while performing tasks or responding to external stimuli. The latter scenario is the setting for this paper. Here, the imaging modality acquires voxel-wise Blood-Oxygen-Level-Dependent (BOLD) measurements [10, 11] at rest and during stimulation or performance of a task. After pre-processing, a general linear or other statistical model [4, 12] is fit to the time course sequence against the expected BOLD response [13, 14, 15]. Statistical Parametric Mapping (SPM) [16] techniques provide voxel-wise test statistics summarizing the association between the time series response at each voxel and the expected BOLD response [3]. The map of test statistics is then thresholded to identify significantly activated voxels [17, 18, 19]. The analysis of fMRI datasets is challenged [20, 21, 22, 23] by factors such as scanner, inter- and intra-subject variability, voluntary/involuntary or stimulus-correlated motion and also the several-seconds delay in the BOLD response as the neural stimulus passes through the hemodynamic filter [24, 23, 25]. Pre-processing  [22, 26] mitigates some of these effects, but additional challenges are presented by the fact that an fMRI study is expected to have no more than 1-3% activated voxels [27, 8]. Also, many activation studies involving high-level cognitive processes have low contrast-to-noise ratios (CNR), throwing up additional challenges as illustrated next.

I-A Activation Detection during Noisy Audiovisual Speech

The most important form of human communication is speech [28, 29, 30], which the brain is adept at understanding even in noisy surroundings. This ability may be due [31] to the brain’s capacity for multisensory integration of independently-acquired visual and auditory input information which reduces noise and allows for more accurate perception [32, 33]. Recently, [31] studied the role of the superior temporal sulcus (STS) in perceiving noisy speech, through fMRI and behavioral experiments, and established increased connectivity between the STS and the auditory or the visual cortex depending on whichever modality was more reliable, that is, less noisy.

[31] provided results on regions of interest (ROIs) drawn on the STS and the auditory and visual cortices. However, the full benefit of fMRI can be realized only if we move beyond assessing cerebral function at the ROI level to understanding it at the voxel level. Reliable voxel-wise activation detection in individual subjects may increase the adoption of fMRI in a clinical setting. All these are potential scenarios with low CNRs where accurate activation detection methods are needed.

I-B Background and Current Practices

Many thresholding methods [34, 35, 36, 37, 38] in fMRI address multiple testing issues in determining significance of test statistics but ignore spatial resolution. Acquired images are instead often spatially smoothed prior to analysis, but such non-adaptive smoothing reduces both the adaptive spatial resolution and the number of available independent tests for activation detection [39]. There are also iterative adaptive smoothing and segmentation methods such as propagation-separation (PS) [39] and adaptive-segmentation (AS) [40] that essentially segment the SPM into activated and inactivated voxels. PS approximately yields a random $t$-field and uses Random Field Theory for segmentation while AS uses multi-scale testing. [40] argued for AS because of its more general development and fewer model assumptions. AS also requires no heuristic corrections for spatial correlation, provides decisions at prescribed significance levels and showed [40] better performance over PS in an auditory experiment. However, AS requires pre-specified bandwidth sequences and ignores correlation within the SPM. So Section of this paper develops theory and methodology for fully automated Fast Adaptive Smoothing and Thresholding (FAST) algorithms that account for correlation and obviate the need for setting all but one parameters. Performance evaluations on real datasets and large-scale simulation experiments are in Section III. Section IV revisits the dataset of Section I-A, while Section V provides discussion. A supplement with sections, figures and tables referenced using the prefix “S” is available.

II Theory and Methods

II-A Preliminaries

Let $\boldsymbol{Y}_{i}$ be the time series vector of the observed BOLD response at the $i$th voxel obtained after preprocessing for registration and other corrections. It is common to relate $\boldsymbol{Y}_{i}$ to the expected BOLD response via the general linear model

[TABLE]

where $\boldsymbol{\epsilon}_{i}$ is a $p$th-order auto-regressive (AR) Gaussian error vector with AR coefficients $\boldsymbol{\phi}_{i}\!\!=\!\!(\phi_{i1},\phi_{i2},\ldots,\phi_{ip})$ and marginal variance $\sigma_{i}^{2}$. Without loss of generality (w.l.o.g.), assume that the design matrix $\boldsymbol{X}$ has the intercept in the first column, the expected BOLD response for the $k$ stimulus levels in the next $k$ columns, and polynomial terms for the drift parameter in the remaining $m$ columns. Therefore, $\boldsymbol{\beta}$ is a coefficient vector of length $d\!=\!k\!+\!m\!+\!1$. We assume that the image volume has $n$ voxels, so $i\!=\!1,2,\ldots,n$ in (1). The parameters $(\hat{\boldsymbol{\beta}}_{i},\hat{\sigma}_{i}^{2},\hat{\boldsymbol{\phi}}_{i})$s are usually estimated via generalized least squares or restricted maximum likelihood. A typical analysis approach then applies (voxel-wise) hypothesis tests with the null hypothesis specifying no activation owing to the stimulus or task. SPMs of the form $\boldsymbol{\Gamma}=\{c^{\prime}\hat{\boldsymbol{\beta}}_{i}\}_{i\in V}$ with appropriate contrasts $c^{\prime}{\boldsymbol{\beta}}_{i}$ are then formulated at each voxel.

Many researchers use models that assume independent or AR(1) errors, while others pre-whiten the time series before fitting (1) under independence. Misspecified models can yield less precise SPMs [41, 42, 43, 44] so here we assume AR($p$) errors, with $p$ assessed by the Bayes Information Criterion (BIC) [45, 46] that trades a fitted model’s complexity against its fidelity to the data. Tests on the SPM $\boldsymbol{\Gamma}$ identify voxels that are activated with the application of the stimulus. Our objective is to develop an approach that adaptively and automatically smooths and thresholds the SPM while incorporating spatial correlation and the fact that the sequential thresholding results in SPMs from truncated distributions. Before detailing our methods, we provide some theoretical development.

II-B Theoretical Development

We assume $t$-distributed SPMs with degrees of freedom large enough for them to be approximately standard normally distributed under the hypothesis of no activation. The SPM has a homogeneous correlation structure, a reasonable assumption with our use of radially symmetric smoothing kernels. We have

Theorem 1**.**

Let ${\boldsymbol{X}}\sim\mathcal{N}_{n}(\boldsymbol{0},{\boldsymbol{R}})$ where $\boldsymbol{X}=(X_{1},\ldots,X_{n})^{\prime}$ and ${\boldsymbol{R}}$ is a circulant correlation matrix with only nonnegative elements such that $\boldsymbol{R}^{1/2}$ also has no negative entries. Writing $\boldsymbol{1}=(1,1,\ldots,1)^{\prime}$, we let $\varrho$ be the sum of the elements in any row of ${\boldsymbol{R}}^{\frac{1}{2}}$. Further, let $X_{(n)}$ be the maximum value of ${\boldsymbol{X}}$, that is, $X_{(n)}=\max{\{X_{1},X_{2},\ldots,X_{n}\}}\equiv\max\boldsymbol{X}$. Then the cumulative distribution function (CDF) $F_{(n)}(x)$ of $X_{(n)}$ is given by $F_{(n)}(x)=P(X_{(n)}\leq x)\geq[\Phi(x/\varrho)]^{n}$, where $\Phi(\cdot)$ is the CDF of the standard normal random variable. The equality holds when $\boldsymbol{R}^{-1/2}$ also has no negative entries.

Proof.

Let $\boldsymbol{Z}\!\sim\!\mathcal{N}_{n}(\boldsymbol{0},\boldsymbol{I}_{n})$ and $Z_{(n)}=\max{\{Z_{1},Z_{2},\ldots,Z_{n}\}}$ have CDF $\Phi_{(n)}(z)\equiv[\Phi(z)]^{n}$. Then $\Phi_{(n)}(x/\varrho)=\mathbb{P}[Z_{(n)}\leq x/\varrho]=\mathbb{P}[\boldsymbol{Z}\leq x\boldsymbol{1}/\varrho]$ so that

[TABLE]

Now $\boldsymbol{R}^{-1/2}$ is also circulant and $\boldsymbol{R}^{1/2}\boldsymbol{R}^{-1/2}\boldsymbol{1}=\tilde{\varrho}\boldsymbol{R}^{1/2}\boldsymbol{1}=\tilde{\varrho}\varrho\boldsymbol{1}$ where $\tilde{\varrho}$ is the sum of the elements of any row of $\boldsymbol{R}^{-1/2}$ and $\tilde{\varrho}=1/\varrho$. If $\boldsymbol{R}^{-1/2}$ has no negative elements, $F_{(n)}(x)=\mathbb{P}[\boldsymbol{X}\leq x\boldsymbol{1}]\leq\mathbb{P}[{\boldsymbol{R}}^{-1/2}\boldsymbol{X}\leq x{\boldsymbol{R}}^{-1/2}\boldsymbol{1}]=\mathbb{P}[\boldsymbol{Z}\leq x\tilde{\varrho}\boldsymbol{1}]=\Phi_{(n)}(x/\varrho)$ and then equality holds in (LABEL:th1). ∎

Corollary 2**.**

For $\boldsymbol{X}$ and $X_{(n)}$ as in Theorem 1, the limiting distribution of $X_{(n)}$ is bounded below by one that lies in the domain of attraction of the Gumbel distribution, and satisfies

[TABLE]

where $a_{n}=\varrho/[n\phi(b_{n}/\varrho)]$ and $b_{n}=\varrho\Phi^{-1}(1-1/n)$, with $\phi(\cdot)$ the standard normal probability density function (PDF).

Proof.

Each element of $\boldsymbol{W}\sim\mathcal{N}_{n}(\boldsymbol{0},\varrho^{2}\boldsymbol{I}_{n})$ has CDF $\Phi(\frac{w}{\varrho})$ and PDF $\phi(\frac{w}{\varrho})/\varrho$. Then $W_{(n)}\!\doteq\!\max\boldsymbol{W}$ has CDF $G_{(n)}(\cdot)$ that satisfies $\lim_{n\rightarrow\infty}G_{(n)}(a_{n}x+b_{n})\!=\!\exp\{-\exp(-x)\}$ with $a_{n}\!=\!\varrho/[n\phi(b_{n}/\varrho)]$ and $b_{n}\!=\!\varrho\Phi^{-1}(1-1/n)$ [47]. The result follows from Theorem 1. ∎

This paper uses $\boldsymbol{R}$ for which $\boldsymbol{R}^{1/2}$ can be shown, using Theorem 2 of [48], to have no negative entries. Then Corollary 2 provides a conservative bound for the quantiles of the limiting distribution of $\max\boldsymbol{X}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{R})$ with the conservatism determined by the negative elements of $\boldsymbol{R}^{-1/2}$.

The thresholding steps yield truncated (and correlated) random variables for potential thresholding in subsequent steps. We account for this added complication by deriving the limiting distribution of the maximum of a correlated sample from a right-truncated normal distribution.

Suppose that $Y_{1},Y_{2},\ldots Y_{n}$ are independent identically distributed (IID) random variables from $\mathcal{N}(0,\varrho^{2})$ but truncated at $\eta$, then each $Y_{i}$ has PDF $\phi_{\eta}^{\bullet}(y;\varrho)$ and CDF $\Phi_{\eta}^{\bullet}(y;\varrho)$, where

[TABLE]

where $I(\cdot)$ is the indicator function. Then $Y_{(n)}=\max{\{Y_{1},Y_{2},\ldots,Y_{n}\}}$ has CDF $\Phi_{\eta,(n)}^{\bullet}(y;\varrho)\!=\!\left[{\Phi(\frac{\min(y,\eta)}{\varrho})}/{\Phi(\frac{\eta}{\varrho})}\right]^{n}$ with limiting distribution as follows.

Theorem 3**.**

Let $Y_{1},Y_{2},\ldots,Y_{n}$ be a sample from (4). Then the limiting distribution of $Y_{(n)}$ belongs to the domain of attraction of the reverse Weibull distribution and satisfies

[TABLE]

for some $\nu>0$. Here $a_{n}^{\bullet}=\eta-{\Phi_{\eta}^{\bullet}}^{-1}(1-1/n;\varrho)$ and $b_{n}^{\bullet}=\eta$.

Proof.

Note that $\eta=\sup\{x\mid\Phi_{\eta}^{\bullet}(x)<1\}$. From Theorem 10.5.2 in [49], for $Y_{(n)}$ to be in the domain of attraction of the reverse Weibull, it sufficient to show that

[TABLE]

for some $\nu>0$ [50]. In our case, the limit holds because $\eta<\infty$. Then upon using L’Hôpital’s rule, we have

[TABLE]

Thus the right-truncated normal distribution satisfies the reverse Weibull condition and converges to the reverse Weibull distribution with $\nu\!\equiv\!1$ in (5). The constants in the theorem are as per extreme value theory [49, 47]. ∎

Theorem 4**.**

Let $\boldsymbol{X}$ be a random vector from the $\mathcal{N}_{n}(\boldsymbol{0},\boldsymbol{R})$ density but that is right-truncated in each coordinate at $\eta$, with ${\boldsymbol{R}}$ and $\varrho$ as in Theorem  1. Then the CDF $F^{\bullet}_{\eta,(n)}(x)$ of $X_{(n)}$ is $F^{\bullet}_{\eta}(x)=\mathbb{P}[X_{(n)}\leq x]\geq\Phi_{\eta,(n)}^{\bullet}(x;\varrho)$.

Proof.

For IID random variables $Y_{1},Y_{2},\ldots,Y_{n}$ from (4), $Y_{(n)}$ has CDF $\Phi^{\bullet}_{\eta,(n)}(x;\varrho)=\mathbb{P}[Y_{(n)}\leq x]=\mathbb{P}[U_{(n)}\leq\frac{x}{\varrho}]=\Phi^{\bullet}_{\frac{\eta}{\varrho},(n)}(\frac{x}{\varrho};1)$, where $U_{(n)}=\max\boldsymbol{U}$ with $\boldsymbol{U}$ a vector of $n$ IID random variables from $\boldsymbol{\Phi}^{\bullet}_{\frac{\eta}{\varrho}}$. Then

[TABLE]

proving the statement of the theorem. ∎

Corollary 5**.**

Let $\boldsymbol{X}$ and $X_{(n)}$ be as in Theorem 4. Then the limiting distribution of $X_{(n)}$ belongs to the domain of attraction of the reverse Weibull distribution, and satisfies:

[TABLE]

where $a_{n}^{\bullet}=\eta-{\Phi_{\eta}^{\bullet}}^{-1}(1-1/n;\varrho)$ and $b_{n}^{\bullet}=\eta$.

Proof.

The result is immediate from Theorems 3 and 4. ∎

II-C Fast Adaptive Smoothing and Thresholding

We propose our FAST algorithm that adaptively and, in sequence, smooths and identifies activated regions by thresholding. We estimate the amount of smoothing robustly or from the correlation structure that we assume is well-approximated by an ellipsoidally-symmetric 3D Gaussian kernel oriented along the three axes and with parameters $\boldsymbol{h}=(h_{1},h_{2},h_{3})$. That is, under the null hypothesis (of no activation anywhere), we assume that the SPM $\boldsymbol{\Gamma}\sim\mathcal{N}_{n}(\boldsymbol{0},\sigma^{2}\boldsymbol{R})$ where $\boldsymbol{R}\!=\!\boldsymbol{S}_{\boldsymbol{h}}$, with $\boldsymbol{S}_{\boldsymbol{h}}$ a circulant smoothing matrix [51]. Let $\boldsymbol{\Gamma}_{(-{\boldsymbol{h}})}\sim\boldsymbol{S}_{\boldsymbol{h}}^{-\frac{1}{2}}\boldsymbol{\Gamma}$. We estimate $\boldsymbol{h}$ and $\sigma$ by maximizing the loglikelihood function

[TABLE]

Both $\boldsymbol{\Gamma}_{(-\boldsymbol{h})}$ and $|\boldsymbol{S}_{\boldsymbol{h}}|$ are speedily computed using Fast Fourier Transforms (FFTs). Starting with the SPM $\boldsymbol{\Gamma}$, obtained as discussed in Section II-A, we propose the algorithm:

Initial Setup. At this point, assume that $\zeta_{i}\equiv 0$ $\forall$ $i$, where $\zeta_{i}$ is the activation status of the $i$th voxel. That is, assume that all voxels are inactive. Set $\zeta_{i}^{(0)}\equiv\zeta_{i}$. Also denote $\boldsymbol{\Gamma}^{(0)}=\boldsymbol{\Gamma}$, and $n_{0}=n$, where $n_{k}$ denotes the number of voxels for which $\zeta_{i}^{(k)}=0$. 2. 2.

Iterative Steps, For $k\!=\!1,2,\ldots,$ iterate as follows:

1. (a)

Smoothing. Smooth $\boldsymbol{\Gamma}^{(k-1)}$ in one of three ways:

1. i.

Adaptive Likelihood maximization (ALL-FAST, pronounced ôl-fast): Maximize (7) given $\boldsymbol{\Gamma}^{(k-1)}$ to obtain ${\boldsymbol{h}}^{(k)}$. Smooth $\boldsymbol{\Gamma}^{(k-1)}$ with $\boldsymbol{S}_{{\boldsymbol{h}}^{(k)}}$ to get $\boldsymbol{\Gamma}^{(k)}$. 2. ii.

Adaptive Model-based smoothing (AM-FAST, pronounced ăm-fast): Use a Markov Random Field (MRF) prior model with parameters estimated using empirical Bayes methods as described in Section S1 to get $\boldsymbol{\Gamma}^{(k)}$ from $\boldsymbol{\Gamma}^{(k-1)}$. 3. iii.

Adaptive Robust smoothing (AR-FAST, ahr-fast): Use [52] to smooth $\boldsymbol{\Gamma}^{(k-1)}$ and get $\boldsymbol{\Gamma}^{(k)}$. 2. (b)

Standardization. Maximize (7) given $\boldsymbol{\Gamma}^{(k)}$ to obtain $\sigma_{(k)}$, ${\boldsymbol{h}}^{(k)}$ and $\varrho_{k}={\boldsymbol{S}}_{{\boldsymbol{h}}^{(k)}}^{\frac{1}{2}}{\boldsymbol{1}}$. Standardize $\boldsymbol{\Gamma}^{(k)}$ by scaling with a robust version of $\sigma_{(k)}$, say $\tilde{\sigma}_{(k)}$. 3. (c)

Adaptive Thresholding. This consists of two steps:

1. i.

For $k=1$, use Corollary 2 to obtain $(a_{n},b_{n})$, otherwise (i.e. for $k>1$) use Corollary 5 to get $(a^{\bullet}_{n_{k-1}},b_{n_{k-1}}^{\bullet})$. In both cases, use $\boldsymbol{R}=\boldsymbol{S}_{{\boldsymbol{h}}^{(k)}}$. 2. ii.

From the Gumbel (for $k=1$) or reverse Weibull distributions (for $k>1$), get

[TABLE]

where $\iota_{\alpha}^{\mathcal{G}}$ and $\iota_{\alpha}^{\mathcal{W}}$ are the upper-tail $\alpha$-values for the Gumbel and the reverse Weibull (with $\nu=1$) distributions, respectively. 3. iii.

Set $\zeta^{(k)}_{i}=1$ if $\zeta^{(k-1)}_{i}=0$ and if the $i$th coordinate of $\boldsymbol{\Gamma}^{(k)}$ exceeds $\eta_{k}$. Let $n_{k}=\sum_{i=1}^{n}\zeta_{i}^{(k)}$. 3. 3.

Termination.

Declare no activation and terminate if $\boldsymbol{\zeta}^{(1)}\equiv 0$. Otherwise, let $J(\boldsymbol{\zeta}^{(k)},\boldsymbol{\zeta}^{(k-1)})$ be the Jaccard Index [53, 54] of the activation maps in the $k$th and $(k-1)$th iterations. If $J(\boldsymbol{\zeta}^{(k)},\boldsymbol{\zeta}^{(k-1)})\geq J(\boldsymbol{\zeta}^{(k+1)},\boldsymbol{\zeta}^{(k)})$, the algorithm terminates – the final activation map is $\boldsymbol{\zeta}^{(k)}$.

Comments

A few comments are in order:

Correlation structure

A circulant correlation structure allows for spatial context in the association between the values of the voxel-wise test statistics, while having the added benefit of speedy computations via the use of FFTs.

Robust Estimation of $\sigma_{(k)}$

The estimate $\hat{\sigma}_{(k)}$ from (7) assumes no activation in $\boldsymbol{\Gamma}_{(k)}$. Ignoring the activation can inflate the estimate, so we obtain $\tilde{\sigma}_{(k)}=\hat{\sigma}_{(k)}\tilde{s}_{(k;w)}/s_{(k)}$ where $s_{(k)}$ is the estimated SD of $\boldsymbol{\Gamma}_{(k)}$, and $\tilde{s}_{(k;w)}$ is its biweight-estimated SD [55], both assuming a zero mean. Specifically, if $\boldsymbol{\Gamma}_{(k)}$ has $i$th component $\Gamma_{i(k)}$, we calculate

[TABLE]

where $e_{i(k),j}={\Gamma_{i(k)}}/(w\tilde{s}_{(k)})$, $\tilde{s}_{(k)}$ the median absolute deviation of $\boldsymbol{\Gamma}_{(k)}$ from 0, and $w=\operatorname*{argmin}_{w\in(0,6)}\tilde{s}_{(k,w)}$.

Comparison with AS

Our FAST algorithms are similar to AS [40] in that they also smooth and threshold iteratively. But there are a few fundamental differences. The AS approach has a set user-specified sequence of bandwidths that smooths $\boldsymbol{\Gamma}^{(k)}$ at each step. In contrast, ALL-, AM- and AR-FAST use likelihood, empirical Bayes and robust methods to optimally determine $\boldsymbol{h}$ at each step. AS also thresholds but uses a general Fréchet extreme value distribution that ignores spatial context and the correlated truncated nature of the random variables that arise from the smoothing and thresholding at each iteration. Our development represents the procedure more accurately because we account for both the correlation structure (with the initial cut-off decided as per the Gumbel distribution) and the truncation (with subsequent cut-offs determined by the reverse Weibull distribution). Our more general $\boldsymbol{h}$ allows for different amounts of smoothing in each axis. Finally, our method is entirely data-driven, with termination declared only if there is no initial activation or when $\mathcal{J}$ between subsequent activation maps decreases.

Two-sided alternatives

Our development here builds from one-sided tests where large values are the extreme values of the SPM. For two-sided alternatives, we use the algorithm individually on the SPM and its negative, but replacing $\alpha$ in (8) with $\alpha/2$. This provides two (disjoint) activation maps, the union of which is the two-sided activation map.

III Performance Evaluations

We studied performance of FAST relative to the most popular and relevant alternatives. Our evaluations were on real and simulated datasets and compared FAST with cluster thresholding (CT) applied with $\alpha=0.001$ per [38], a second-order neighborhood and number of voxels in cluster determined by  [56]’s 3dClustSim function, threshold-free cluster enhancement (TFCE) [37], permutation-based testing (PBT) [57], AS and PS (applied as AWS or adaptive-weighted smoothing [58]). We used $\alpha=0.01$ and $\alpha=0.05$ in FAST to obtain insight into the role of $\alpha$. We used R packages RFASTfMRI for FAST, fMRI for AS and AWS, AnalyzeFMRI for CT and permuco for TFCE and PBT.

III-A Finger-Tapping Experiments

Our first set of evaluations used the 12 replicated SPMs [24, 25] from the right-hand (RH) and left-hand (LH) finger tapping study of a RH-dominant male. For each method, Figure 1 displays

the summarized Jaccard similarity ($\ddot{\omega}$) between the activation maps from the 12 replications. For the RH experiments, AR-FAST showed the highest reliability of detected activation with $\alpha=0.05$. AR-FAST at $\alpha=0.01$ and TFCE were marginally behind and AS and AWS also doing reasonably. TFCE was a bit better than AR-FAST for the LH experiments. The generally low $\ddot{\omega}$ for all methods points to potential issues in data quality and processing [54].

III-B Experiments on Simulated Phantom Data

Our next set of examples evaluated performance on phantom data simulated using (1) under different conditions.

III-B1 Motif and Stripes

We first study performance using the simulation setup of [40]. We thank K. Tabelow for readily sharing code that created the motif and three striped ($16\times 16$, $32\times 32$, $64\times 64$) phantoms of Figures 2LABEL:sub@motif-LABEL:sub@str64. The phantoms have 14, 50, 28 and 47% truly activated voxels, or more than the 1-3% expected in typical fMRI experiments. We used $\beta$s as per [40] and CNRs of between 0.75 to 2.68 for the motif and 1 to 2 for the stripes. These examples are of two-sided alternatives. All simulations had AR(1) errors with $\rho=0.3$. For AS and AWS, we adopted the maximum bandwidth sequence values ($h_{k}^{*}=3.06$, 1, 2 and 3) in [40] for the four respective phantoms as the best-case specific choices. Figure 2LABEL:sub@JI-Tabelow summarizes performance. There is no clear overall winner but AM-FAST, PBT and TFCE find no activation at all (Figure S1). ALL-FAST and AR-FAST, in that order, perform creditably in some situations but not in others.

III-B2 Large-scale study with modified Hoffman phantom

The phantoms in [40], with uniform underlying structure ($\beta_{0}$) and no drift, but varying CNR, are not particularly representative of cerebral activation and do not provide much insight into performance of different activation methods. So we performed a large simulation study using a more realistic phantom and experimental setup that matches (1). We used a modified version of the digitized $128\!\times\!128$ 2D Hoffman phantom [59] of Positron Emission Tomography, with 3465 in-brain pixels, representing two types (say, A and B) of anatomic structures – the latter has 138 deemed truly activated pixels in two distinct regions (Figure 3).

The $i$th pixel in the phantom had values $\boldsymbol{\beta}_{i}\!=\!(\beta_{i0},\beta_{i1},\beta_{i2})$ in (1) as per its location (see Figure 3).

As in (1), the design matrix $\boldsymbol{X}$ had the intercept in the first column. The second column had the hemodynamic response function (HRF) [7] convolved with the input stimulus time series that alternated as 16 on-off blocks of 6 time-points each. The third column of $\boldsymbol{X}$ represented linear drift and was set to $t$ ($t\!=\!1,2,\ldots,96$). As per (1), AR($p$) Gaussian errors were simulated for different $p$ and at each pixel. Specifically, for each $p$, we considered AR coefficients for a range of $\boldsymbol{\phi}\!\equiv\!\boldsymbol{\phi}_{i}$s with coefficients $(\phi_{1},\phi_{2},\ldots,\phi_{p})$ that were, with lag, (a) all equal, (b) decreasing, (c) increasing, (d) first decreasing, then increasing and (e) first increasing and then decreasing. We restricted $\sum^{p}_{j=1}\phi_{j}=0.9$ to ensure stationary solutions. So $\phi_{1}\equiv 0.9$ for all AR(1) cases. For $p\!=\!2,3,4$, $\phi_{i}\equiv 0.9/p$ for the equal AR coefficients scenario and as per Table I for the other cases.

Finally, $\sigma_{0}=1200,800,600$ to correspond to very low to moderate CNR = 0.50, 0.75, and 1.0. By design, our SNRs were 10 times our CNRs. Time series images were simulated using (1) and the setup of Fig. 3 and Table I.

For each pixel, we fit (1) with $\hat{p}$ chosen from {0,1,2,3,4,5} using BIC. SPMs were generated as per Section II-A. Figure 4 provides sample SPMs and activation maps with the three top-performing methods: AR-FAST, AM-FAST and AS, for the three CNR settings with AR(2) errors and coefficients decreasing with order.

(See Figures S2 for activation maps using all methods and Figures S2–S5 for cases with other $p$ and/or decreasing order.) All methods do poorly for CNR=0.5 in this example, but AS correctly finds a few activated voxels. AM- and AR-FAST perform very well for CNR=0.75 and 1.0. Other methods – in particular ALL-FAST, CT, TFCE and PBT – barely find activation on this example.

To more fully understand performance, we replicated our experiment 25 times for each simulation setting. Figure S6 shows performance in estimating $p$, with over-estimation and mild under-estimation for large and small values of the true AR order. The pattern broadly holds at all CNRs and $\phi$s. We now discuss performance of activation detection methods on SPMs obtained upon fitting AR($\hat{p}$). Figure 5 displays overall performance, in terms of $\omega$, of all methods in the decreasing-$\phi$ cases. Performances with other types of $\phi$s are similar (Figure S7). Figures S8–S12 display the number of activated voxels and the false and true positive rates (FPRs and TPRs). FAST methods are among the top performers at all CNRs: AM- and AR-FAST do very well for experiments with orders other than AR(1). (All methods do poorly at CNR=0.50 in the AR(1) case: we surmise that is because of highly-correlated noise obtained with $\phi=0.9$.) AS and AWS are the next best performers but CT, TFCE and PBT perform very poorly with TPRs of below 25%. FAST methods have very high TPRs but FPRs of up to 0.2, 0.3 and 0.5% for ALL-, AM- and AR-FAST for $\alpha=0.05$, low CNR and $p$. Overall, the slightly higher FPRs of the FAST methods are overwhelmed by their vastly higher TPRs, leading to their having the highest $\mathcal{J}$s.

Our threshold $\alpha$ has a role, with smaller values performing better at higher CNRs and conversely. We suggest $\alpha\approx 0.05$ for low-CNR tasks and $\alpha\approx 0.01$ for high-CNR tasks. We suggest determining low or high CNR scenarios accordingly as whether the upper percentile of the estimated voxel-wise CNRs is less than the standard normal upper percentile (2.33) or not: the upper percentile of the estimated CNRs is chosen to include an activated voxel (if such exists) in the CNR determination. In our studies, ALL-FAST required more Step 2 iterations, but was computationally faster than AM-FAST AR-FAST and had lower TPR, FPR and $\omega$, especially at low CNRs. Regardless, our methods were the fastest among all methods.

III-C Performance in Null Activation Scenarios

III-C1 Resting-State Dataset

A reviewer’s suggestion led us to apply our FAST algorithms on SPMs obtained upon fitting (1) to a resting state dataset [60, 61], with no activation identified even at $\alpha=0.05$. This zero FPR (when CNR=0) as opposed to the small FPR in low-CNR experiments may be due to Step 3 of our algorithm correctly allowing more Step 2 iterations to attenuate stray high-valued SPM voxels – termination is earlier in the low CNR cases given the spatially located weaker-signal peaks. For higher CNRs, Step 3 again adaptively admits more smoothing iterations that dampen stray high values in the SPM without substantially degrading the true high-signal peaks.

III-C2 Null-simulated SPMs

Another reviewer was concerned about multiple significance. Our use of $\alpha$ is to set a threshold and not as a significance level: still, experiments on simulated null SPMs (Section S2) indicate no practical concerns.

IV Activation During Perception of Noisy Speech

The dataset, provided as data6 in the AFNI tutorial [56], is originally from an fMRI study [31] where

a subject heard and saw a female volunteer speak words, separately, in two different formats. The audio-reliable setting had the subject clearly hear the spoken word but see a degraded image of the speaker while the visual-reliable case had the subject clearly see the speaker vocalize the word but the audio was of reduced quality. There were three experimental runs, each consisting of a randomized design of 10 blocks, equally divided into blocks of audio-reliable and visual-reliable stimuli. $\mbox{T}_{2}^{*}$-weighted images with volumes of $80\times 80\times 33$ (with voxels of dimension $2.75\times 2.75\times 3.0\ mm^{3}$) from echo-planar sequences (TR=2s) were obtained over $152$ time-points. Our interest was in determining activation corresponding to the audio ($H_{0}:\beta_{a}=0$) and visual ($H_{0}:\beta_{v}=0$) tasks. At each voxel, we fitted AR models for $p=0,1,2,3,4,5$ and chose $p$ with the highest BIC. Figure 6 uses AFNI and Surface Mapping (SUMA) to display activated regions obtained using AR-FAST on the SPM: see Figure S13 for maps drawn from ALL-FAST, AS, AWS and CT. We used $\alpha=0.01$ because of the high (greater than 4) upper percentile of the voxel-wise estimated CNRs. Most of the activation occurs in Brodmann areas 18 and 19 (BA18 and BA19) which comprise the occipital cortex and the extrastriate (or peristriate) cortex. In humans with normal sight, this area is for visual association where feature-extraction, shape recognition, attentional and multimodal integrating functions occur. We also see increased activation in the STS, which recent studies [62] have related to distinguishing voices from environmental sounds, stories versus nonsensical speech, moving faces versus moving objects, biological motion and so on. ALL-FAST performs similarly as AR-FAST, while the other methods also identify the same regions but they identify a lot more activated voxels, some of which appear to be false positives. Although a detailed analysis of the results of this study is beyond the purview of this paper, we note that AR-FAST finds interpretable results even when applied to a single subject high-level cognition experiment.

V Discussion

We propose a new fully automated fast adaptive smoothing and thresholding algorithm suite called FAST with the ability to detect activation in low-signal settings. Three variants – ALL-FAST, AM-FAST and AR-FAST – are proposed with AR-FAST generally recommended because of its consistent good performance across a range of low-CNR experiments and real datasets. AM-FAST’s performance, while good, is more variable, while ALL-FAST appears to undersmooth but performs better in two-sided activation detection scenarios. Our methodology realistically accounts for both spatial correlation structure and is also developed under more accurate extreme value theory. Our algorithm suite is implemented in a R package RFASTfMRI available at https://github.com/ialmodovar/RFASTfMRI and is fully automated with one threshold choice for which we provide easily-implemented guidance. This contrasts with AS and AWS that require setting maximum smoothing bandwidths related to the expected diameter of activated regions [40] – a determination that may require considerable dexterity and is ambivalent when different-sized activation regions are expected.

A reviewer has pointed to the joint detection-estimation literature [63, 64] where estimation of the HRF and activation detection occur jointly. The FAST, AS and AWS algorithms can be placed in a related framework, with the distinction that the estimation step is of a more spatially consistent (smoothed) SPM. We also agree with another reviewer on other ways of ensuring spatial contiguity such as through Markov Random Field priors [65] and on the need to incorporate approaches also allowing for nonhomogeneous smoothing. Our algorithms converge by construction and are guaranteed to terminate. They are also seen to have good overall performance but, as observed by a reviewer, establishing the optimality properties and conditions and assumptions governing such properties may provide more solid theoretical grounding for FAST and improve its understanding and widen its applicability. Developing FAST for more sophisticated time series and spatial models, including in the context of complex-valued fMRI [66] as well as increased use of diagnostics in understanding activation and cognition are other important research areas and directions that would benefit from further attention.

Acknowledgment

The authors sincerely thank four anonymous reviewers and an Associate Editor whose helpful and insightful comments on an earlier version of this article greatly improved its content.

Supplementary Materials

S1 Model-based smoothing

In this section, we describe the model-based smoothing used in the AM-FAST algorithm. The starting point of our methodology is a SPM $\boldsymbol{\Gamma}$ and our smoothed estimator is given by $\hat{\boldsymbol{\Gamma}}=\boldsymbol{H}_{\boldsymbol{\beta}}\boldsymbol{\Gamma}$ where $\boldsymbol{H}_{\boldsymbol{\beta}}$ is a smoothing matrix with parameter $\boldsymbol{\beta}$. We write $\boldsymbol{H}_{\boldsymbol{\beta}}$ = $\boldsymbol{\Lambda}_{\boldsymbol{\beta}}(\boldsymbol{\Lambda}_{\boldsymbol{\beta}}+\boldsymbol{R})^{-1}$ where

[TABLE]

with $\otimes$ denoting the Kronecker product between matrices, $\boldsymbol{I}_{q}$ is a $q\times q$ identity matrix and $\boldsymbol{J}_{q}$ is a $q\times q$ tridiagonal matrix of the form

[TABLE]

with $i$th eigenvalue given by $j_{i}=2-2\cos((i-1)\pi/q)$ for $i=1,2,\ldots,q$. In the above representation, $\boldsymbol{\Lambda}_{\boldsymbol{\beta}}^{-1}$ is split into three components to represent the three axes of the volume image. Indeed $\boldsymbol{\Lambda}_{\boldsymbol{\beta}}$ represents the inverse of the dispersion matrix of a 3-D first-order Gaussian Markov Random Field with a neighborhood structure for each interior voxel given by the six nearest neighbors (one nearest in each of the two directions in each plane). Also $\boldsymbol{\beta}=(\beta_{1},\beta_{2},\beta_{3})$ with $\beta_{i}$ measuring the strength of the interaction between the neighbors in the $i$th plane. For the edge voxels, the neighborhood structure is similar and given by the nearest neighboring voxels in each of the two directions in each plane, provided they exist in the imaging grid.

Our smoothed SPM is a penalized maximum likelihood estimator and is the same as the maximum a posterioiri (MAP) or posterior mean estimator of $\boldsymbol{\mu}$ when we have $\boldsymbol{\Gamma}\mid\boldsymbol{\mu}\sim N(\boldsymbol{\mu},\boldsymbol{R})$ and $\boldsymbol{\mu}\sim N(\boldsymbol{0},\boldsymbol{\Lambda}_{\boldsymbol{\beta}})$. So, the regularization parameter $\boldsymbol{\beta}$ may be chosen using an empirical Bayes estimate. This estimate is given by

[TABLE]

We consider $\boldsymbol{R}=\sigma^{2}\boldsymbol{I}_{n_{1},n_{2},n_{2}}\equiv\sigma^{2}\boldsymbol{I}_{n_{1}}\otimes\boldsymbol{I}_{n_{2}}\otimes\boldsymbol{I}_{n_{3}}$ in our application. Then the eigenvalues of $\boldsymbol{\Lambda}_{\boldsymbol{\beta}}+\sigma^{2}\boldsymbol{I}$ are all real and given by

[TABLE]

Let $\boldsymbol{C}_{q}^{\prime}$ be the $q\times q$ Type-2 1D Discrete Cosine Transform (DCT) matrices and $\boldsymbol{C}_{q}$ be the inverse DCT. Then $\boldsymbol{C}^{\prime}=\boldsymbol{C}^{\prime}_{n_{1}}\!\otimes\boldsymbol{C}^{\prime}_{n_{2}}\!\otimes\!\boldsymbol{C}^{\prime}_{n_{3}}$ is the Type-2 DCT in the 3D volume of $n_{1}\times n_{2}\times n_{3}$ voxels and $\boldsymbol{C}=\boldsymbol{C}_{n_{1}}\!\otimes\boldsymbol{C}_{n_{2}}\!\otimes\!\boldsymbol{C}_{n_{3}}$ is the inverse 3D DCT. Then $\boldsymbol{W}=\boldsymbol{C}^{\prime}\boldsymbol{\Gamma}$ is the 3D DCT of $\boldsymbol{\Gamma}$ with entries $w_{ijk}$. Then (E11) reduces to

[TABLE]

Further, the smoothed estimator $\boldsymbol{\Gamma}$ can be obtained by the inverse DCT of $(\boldsymbol{L}-\sigma^{2}\boldsymbol{I}_{n})\boldsymbol{L}^{-1}\boldsymbol{W}$, where $\boldsymbol{L}$ is the diagonal matrix with the diagonal element corresponding to the $(i,j,k)$th voxel given by (E12).

S2 Performance of FAST Algorithms on Simulated Null-Activated SPMs

The second reviewer expressed continued concerns regarding multiple significance. Our use of $\alpha$ is only to set a threshold: its connections with testing in deciding on this threshold notwithstanding, our algorithm is not built on, and does not use, the actual nominal significance level. Further, performance of FAST on the resting state dataset of Section III-C1 did not raise concerns of identifying false positives. Nevertheless, we also performed a comprehensive and thorough 2D simulation study to evaluate the issue of false positives. We generated SPMs under the null hypothesis of no activation, and used ALL-, AM- and AR-FAST for different $\alpha$ and evaluated the incidence of even one pixel in an image (falsely) identified as activated. Specifically, we generated SPMs $\boldsymbol{\Gamma}\sim N(\boldsymbol{0},\boldsymbol{R})$ where $\boldsymbol{R}$ is a 2D circulant matrix, with first order autocorrelation structure indexed in each dimension by $\rho$. (For simplicity, we used a radial correlation structure in the simulation, but not in the estimation and activation detection.) We used $\rho\in\{0,0.01,0.025,0.05,0.075,0.1,0.25,0.5,0.75,0.99\}$ to provide a comprehensive evaluation even though in our experience most SPMs from real datasets have estimated $\rho$s substantially lower than 0.1. Our SPMs were of dimensions $128\times 128$. For each $\rho$, we simulated 1000 SPMs to account for simulation variability in drawing our conclusions and then applied both the FAST algorithms for a range of $\alpha\in\{0.001,0.01,0.025,0.05,0.075,0.1\}$.

Table S1 represents the number of cases (out of 1000) where ALL-FAST, AM-FAST or AR-FAST for different $\rho$ and $\alpha$ found (falsely) even one activated voxel. For $\rho\leq 0.75$, AR-FAST had no false positives. while for AM- and AR-FAST, there were (out of several thousands) only a handful of cases where even one pixel was falsely identified as activated and, in the case of AM-FAST, these all happened at values of $\alpha>0.05$. (Indeed, even with ALL-FAST, for $\alpha\leq 0.075$, these cases with even one false positive identification mostly had 1-2 pixels in that image falsely identified as activated.) Thus, as with the case of resting state data of Section III-C1, and for all pratical fMRI settings, FAST performs very well in scenarios of no true activation.