Fast Bayesian estimation of brain activation with cortical surface fMRI data using EM

Consider cortical surface fMRI (cs-fMRI) data from a scan, a time series of length $T$ represented as $\mathbf{y}_{v}\in\mathbb{R}^{T}$ for locations $v=1,\ldots,N$ for $N$ locations in order to model data from a single subject. These data are gathered while a subject completes $K$ different tasks at specific time intervals. The expected BOLD response to task $k$, $\mathbf{x}_{k}\in\mathbb{R}^{T}$, is generated by convolving the binary stimulus timing function (on/off) with a canonical haemodynamic response function (HRF) representing the time-lagged BOLD response that follows neuronal firing. After preprocessing the data to reduce noise, eliminating the baseline signal, and inducing residual independence as outlined in A, the general linear model can be written in the form

where $\beta_{v,k}$ is the activation effect in terms of percent signal change at location $v$ for task $k$. Note that the task-based regressors $x_{k,v}$ vary across locations $v$. This is because our preprocessing includes location-specific (or spatially variable) prewhitening to account for differences in residual autocorrelation across the cortex (Parlak et al., 2022). It is typical to fit this model in a massive univariate manner which we refer to as the classical GLM, which is able to be fitted quickly, even to high-resolution data. However, this advantage comes at the cost of ignoring the spatial dependence in brain activation, which reduces model power significantly when compared to the Bayesian GLM, as shown in Spencer et al. (2022).

In order to facilitate spatial Bayesian modeling, we rewrite the model in Equation (1) to represent data across all locations and times at once in Equation (2). Here, the preprocessed response $\mathbf{y}=(\mathbf{y}_{1}^{\prime},\ldots,\mathbf{y}_{N}^{\prime})^{\prime}\in\mathbb{R}^{TN}$ at times $t=1,\ldots,T$ and locations $v=1,\ldots,N$ is explained by linear effects $\boldsymbol{\beta}_{k}\in\mathbb{R}^{N}$ for each of $K$ tasks ($\mathbf{X}_{k}=\text{block-diagonal}(\mathbf{x}_{1,k},\ldots,\mathbf{x}_{N,k})\in\mathbb{R}^{TN\times N}$) and an error term $\mathbf{e}\in\mathbb{R}^{TN}$.

This linear model facilitates incorporating spatial dependence, in contrast to the $N$ separate linear models constituting the classical GLM, which does not include spatial dependence at the model-fitting stage. We fit this model across an entire cortical hemisphere, and the two cortical surface hemispheres are analyzed separately for computational reasons, which does not affect the results because they are represented by separate surface meshes that are physically distinct.

The SBSB GLM incorporates spatial dependence in the activation amplitudes $\boldsymbol{\beta}_{k}$ using a special class of Gaussian Markov random field (GMRF) processes known stochastic partial differential equation (SPDE) priors (Mejia et al., 2020; Spencer et al., 2022; Lindgren et al., 2011; Bolin and Lindgren, 2013). The specific construction of the prior is

where $\boldsymbol{\Psi}_{k}\in\mathbb{R}^{N\times n}$ maps the original $N$ data locations to a triangular mesh consisting of $n$ vertices. $\mathbf{C}$ is a fixed diagonal matrix describing the relative mesh vertex precisions. $\mathbf{G}$ is a fixed sparse matrix describing the neighborhood structure of the mesh vertices, in which entries corresponding to neighboring locations take nonzero values. The mesh is usually created to maximize the minimum interior angle across the triangles, and extra vertices may be added along the boundaries to avoid adverse boundary effects. An example of this mesh structure can be seen in Figure 7. The default setting for the INLA implementation used in Mejia et al. (2020) sets the priors for $\tau_{k}$ and $\kappa_{k}$ to be log-Normal with mean 2 and variance 2. Under the SPDE representation in $\mathbb{R}^{d}$, the variance of the random field can be represented as $\phi_{k}=\Gamma(\nu)\left(\Gamma(\alpha)(4\pi)^{d/2}\kappa_{k}^{2\nu}\tau_{k}^{2}\right)^{-1}$ where $\nu=\alpha-d/2$. Spectral theory shows that an integer must be picked for $\alpha$ to obtain a Markov field. The form of the prior in Equation (3) assumes $\alpha=2$, which makes $\nu=1$ for a surface, which has dimension $d=2$. Therefore, the variance of the random field simplifies to

The product $\kappa_{k}^{2}\tau_{k}^{2}$ in the variance remains the same if $\kappa_{k}^{2^{\prime}}=c\kappa_{k}^{2}$ and $\tau_{k}^{2^{\prime}}=\tau_{k}^{2}/c$ for any constant $c$, and early experiments showed that this can cause the EM algorithm to produce values for $\kappa_{k}^{2}$ and $\tau_{k}^{2}$ in which one approaches infinity while the other approaches zero. Therefore, we rewrite the precision in the SPDE prior in terms of the variance of the random field as

In this form, $\phi_{k}$ controls the variance and $\kappa_{k}^{2}$ controls the spatial correlation of the random field.

To avoid blurring across structural boundaries, cortical surface data are already represented on a triangular mesh.

For mathematical convenience, the model in Equation (2) can be rewritten using the notation in Equation (3) to include all $k$ task covariates in matrix form:

where $\mathbf{X}\in\mathbb{R}^{TN\times NK}$ is a column-bound matrix of the covariates, $\boldsymbol{\Psi}\in\mathbb{R}^{NK\times nK}$ is a block diagonal matrix projecting the covariate values onto the triangular mesh, and $\mathbf{w}\in\mathbb{R}^{nK\times 1}$ is the vector of coefficients corresponding to the covariates on the mesh.

In order to provide a computational alternative to INLA for the SBSB GLM, here we derive an EM algorithm (Dempster et al., 1977; Gelman et al., 2013) to find the mode of the posterior distribution of $\mathbf{w}$, and by extension $\boldsymbol{\beta}$, along with maximum likelihood parameter estimates for $\boldsymbol{\theta}=(\kappa_{1}^{2},\ldots,\kappa_{K}^{2},\phi_{i},\ldots,\phi_{K},\sigma^{2})^{\prime}$. This algorithm could be easily extended to incorporate priors on these hyperparameters to account for uncertainty in their values and avoid underestimation of the posterior variance of the latent fields. In addition, we obtain an estimate of the posterior precision of $\boldsymbol{\beta}_{k}$, which can be used in conjunction with the excursions method developed by Bolin and Lindgren (2015, 2017, 2018) to determine areas of activation in the latent fields.

Mejia et al. (2020) develops a joint multi-subject model that combines results across multiple single-subject models by taking a weighted average of their hyperparameters ($\boldsymbol{\theta}_{G}$), which results in the multi-subject posterior distribution for the SPDE hyperparameters. Our own multi-subject model differs because each single-subject result returns a point estimate of the maximum a posteriori values for $\boldsymbol{\theta}_{m}$ for each subject $m$, rather than a distribution. For our multi-subject model, we assume that the values found in the single-subject analysis are noisy estimates of the fixed true values of $\boldsymbol{\theta}_{G}$. Therefore, in order to find the maximum a posteriori values for $\boldsymbol{\theta}_{G}$ across all $M$ subjects, the multi-subject estimate of the SPDE hyperparameters is found as

where $\hat{\boldsymbol{\theta}}_{m}$ is a vector of the maximum a posteriori values found for each subject $m$ and $\lambda_{m}$ is a weight assigned to each subject’s observation, which can be assigned based on the amount of data in each subject’s scan and their respective noise levels. We then adjust the posterior distributions of $\mathbf{w}_{m}$ found for each of the $M$ subjects to reflect these updated hyperparameter values, calculating $\mathbf{Q}_{G}$ as in Equation (10) using $\hat{\boldsymbol{\theta}}_{G}$:

Next, $H$ samples are taken from these updated posterior distributions for the subjects, and a contrast vector $\mathbf{c}\in\mathbb{R}^{KM\times 1}$ is specified in order to find the multi-subject posterior distribution for a linear combination of tasks and subjects. For example, if an equally-weighted average of the task activation amplitudes across all subjects for the first task is desired, then the contrast vector would take the form

A contrast matrix is then created using the Kronecker product such that $\mathbf{C}=\mathbf{I}_{N}\otimes\mathbf{c}^{\prime}\in\mathbb{R}^{N\times NKM}$ and the multi-subject posterior draws for the average activation amplitude for task 1 is found to be $\boldsymbol{\beta}^{*}=\mathbf{C}\boldsymbol{\beta}^{(1:H)}$, where $\boldsymbol{\beta}^{(1:H)}=(\boldsymbol{\beta}_{1}^{*}\cdots\boldsymbol{\beta}_{H}^{*})\in\mathbb{R}^{NKM\times H}$, and $\boldsymbol{\beta}_{h}^{*}\in\mathbb{R}^{NKM\times 1}=(\boldsymbol{\beta}_{1,h}^{*\prime}\cdots\boldsymbol{\beta}_{M,h}^{*\prime})^{\prime}$ is a column-bound matrix constructed where $\boldsymbol{\beta}_{m,h}^{*^{\prime}}\in\mathbb{R}^{NK\times 1}=\boldsymbol{\Psi}\mathbf{w}_{m,h}^{*}$ and $\mathbf{w}_{m,h}^{*}\in\mathbb{R}^{nK\times 1}$ is the vector of the $h$th draw from the posterior distribution of $\mathbf{w}_{m}^{*}$. This principled combination of single-subject results provides the distribution of a multi-subject average activation amplitude effect, which can be used with the excursions method to determine population-level areas of activation, as discussed in the next section.

In addition to assessing an estimate of an effect size for brain activation, it is also desirable to determine which regions of the brain are significantly activated above a given threshold. The Bayesian GLM proposed by Mejia et al. (2020) and validated by Spencer et al. (2022) takes advantage of the posterior distributions for the activation effects at all data locations on the cortical surface or subcortical volume through the use of the excursions set method proposed by Bolin and Lindgren (2015). The excursions method examines the joint distribution of the activation $\boldsymbol{\beta}_{k}=(\beta_{1,k},\ldots,\beta_{V,k})^{\prime}$ for task $k$ and location $v=1,\ldots,V$. Specifically, the excursions method finds the largest set $\mathcal{V}=(a_{1},a_{2},\ldots,a_{p})^{\prime}$ of data locations such that $P(\beta_{a_{1},k}>\gamma,\beta_{a_{2},k}>\gamma,\ldots,\beta_{a_{p},k}>\gamma)\geq 1-\alpha$ for a given threshold $\gamma$ and credible level $1-\alpha$. This is done using the posterior distributions, either at the single- or multi-subject, for the activation amplitudes $\boldsymbol{\beta}_{k}$. Implementation of the excursions method has been done in the excursions package in R (Bolin and Lindgren, 2015, 2017, 2018), which is used by the BayesfMRI package.

We begin with simulations for a single subject, simulating data under nine different conditions. These conditions are the combination of resampling resolution at $n=\{2,000;5,000;10,000\}$ cortical surface vertices per hemisphere, and setting the number of tasks to be $K=\{2,5,8\}$. Under these conditions, the length of the scan is set to $T=300$, and the maximum value for the spatial coefficient is set to 2. Under these conditions, the average value for the nonzero true coefficients is approximately 0.5 due to spatial smoothing in the data generation. The error term was set to have a variance of 1, and the error was assumed to be temporally independent.

Each of the nine conditions was used to generate 10 datasets, which were all fit to the classical GLM, and the SBSB GLM using the INLA and EM implementations. The INLA implementation allows for the use of the PARDISO linear algebra library (Alappat et al., 2020; Bollhöfer et al., 2020, 2019), which allows for multiple processor threads to be used in expensive matrix solve computations. In these trials, the INLA implementation is allowed to use 6 threads as an advantageous setting for analysis. We compare the results in terms of the amount of time to perform the analysis after preprocessing (Figure 1(a)) and in terms of their root mean squared error (RMSE) (Figure 1(b)). These results show that the EM implementation is faster in terms of computation time than the INLA implementation when the number of tasks ($K$) is relatively small. The speeds are similar when $K=5$, but the INLA implementation is currently faster than the EM implementation due to faster solve times stemming from its use of the PARDISO library. Implementing the EM algorithm to take advantage of a solver that allows the use of multiple cores is likely to further increase its speed, potentially outstripping the INLA implementation in terms of computation time across all settings. However, no parallelization libraries are currently available across UNIX and Windows systems that are also compatible with the RcppEigen library. Using the RMSE measure, the EM implementation’s accuracy is similar to that of the INLA implementation, though it is important to note that the INLA implementation failed to converge in the analysis of 13 of the 90 datasets, and thus did not produce estimates in these cases. Twelve of these datasets had the condition that $K=2$, and of those, 7 had $n=10000$. This suggests that the EM algorithm is more stable than the INLA implementation when evaluating experiments with few tasks as resolution increases. A visualization of an estimated and true coefficient surface for the sampling condition when $n=5000$ and $K=5$ can be seen in Figure 2(a).

Under the simulation condition with resolution $n=5000$ and number of tasks $K=5$, simulated data were generated for two runs in a single session to allow for sharing of the hyperparameters $\boldsymbol{\theta}=(\kappa_{1},\ldots,\kappa_{5},\phi_{1},\ldots,\phi_{5},\sigma^{2})^{\prime}$ across the runs. The INLA and EM implementations of the Bayesian GLM were then used to analyze the first run by itself and the two runs combined. The effect estimates and activations illustrating the difference in the 1- and 2-run analyses can be seen in Figure 2. As there is a small run effect simulated in the generation of the data, the true values of the activation amplitude are not identical between the two runs. Therefore, the true coefficient nonzero region in the single run is not centered exactly in the same location as the true nonzero coefficient region across both runs. In terms of both the coefficient estimation and activation, the EM implementation performs very similarly to the INLA implementation. The posterior standard deviations for the activation amplitudes for a single-run analysis are shown for the EM and INLA implementations in Figure 8, which shows that the posterior standard deviations are slightly lower in the INLA implementation due to the prior imposed on the range and variance hyperparameters.

In order to compare the performance of the multi-subject analysis between the EM and INLA implementations of the SBSB GLM and the classical GLM, we generated a simulated population of 100 subjects with true amplitude maps with location translations from a population amplitude map by an average of 5mm. Each subject’s data were simulated at a resolution of about 5,000 vertices per hemisphere, with a TR of 1 second for a total length of 300 seconds. Two tasks were included in the simulation, each with a maximum true signal-to-noise ratio of 2 in the simulation before the BOLD response and design matrices were scaled. Each subject’s data was analyzed separately following the single subject analysis for the EM and INLA implementations of the SBSB GLM and the classical GLM as outlined in Section 2.1. Next, 100 different subsets of size 10 were drawn from the population, and the multi-subject analysis was performed using the single-subject analysis results from the subjects within the subset. For each multi-subject analysis, the square root of the mean squared error was found with respect to the true simulated population amplitudes. The Dice coefficient (Dice, 1945) is a measure of set overlap between sets $\mathcal{A}$ and $\mathcal{B}$ found as

where $|\mathcal{X}|$ is the number of elements in set $\mathcal{X}$ and $\mathcal{X}\cap\mathcal{Y}$ is the intersection of sets $\mathcal{X}$ and $\mathcal{Y}$. The Dice coefficient can take values between 0 and 1, and higher values indicate more overlap between sets. We calculate the Dice coefficient between the set of locations with true positive activation amplitudes in the population and the set of locations determined to have amplitudes significantly greater than zero ($\gamma=0\%$) in each subset for each modeling method, as outlined in Section 2.4.

Figure 3(a) shows the true value for the multi-subject activation amplitude for a single task, as well as the locations where the true amplitudes are greater than $\gamma=\{0\%,0.5\%,1\%\}$ signal change in the analysis of a single 10-subject subset. The EM and INLA implementations of the SBSB GLM produce similar amplitude estimates and areas of activation, showing attenuation between the true values and the estimates. In contrast, the point estimates from the classical GLM show many values that are appreciably distant from zero in the true zero-valued region, while only finding one location that is significantly different from zero within the true nonzero-valued region. None of the modeling methods found any locations in which the amplitude was found to be significantly greater than 0.5% due to the attenuation expected in the Bayesian model and the low power of the classical GLM. For the Bayesian GLM implementations, the areas of activation are found using the excursions method, and are found as areas deemed jointly greater than threshold $\gamma$ with probability 0.01. For the classical GLM, the areas of activation are found by hypothesis test after fixing the false discovery rate at 0.01.

Figure 3(b) shows the calculated values for the square root of the mean squared error (RMSE) and the Dice coefficient for each of 100 subsets of ten subjects each in order to compare the accuracy between the different methods. The true amplitude values for the whole population of 100 simulated subjects was used to determine the values of the measures of accuracy, rather than the true values for each of the subsets, as population-level inference is the goal of the multi-subject analysis. The EM implementation of the SBSB GLM exhibits generally lower values for the RMSE than the INLA implementation and the classical GLM, while exhibiting only slightly lower values for the Dice coefficient. The EM estimates exhibit slightly less attenuation than in the INLA implementation (shown in Figure 9), which is expected because no prior distributions are assumed on the hyperparameters in the EM implementation. The classical GLM performs poorly in terms of the Dice coefficient compared to both implementations of the SBSB GLM. These analyses suggest that the multi-subject analysis using the EM implementation of the SBSB GLM is likely to deliver powerful results in the analysis of real cs-fMRI data.