A canonical polyadic tensor basis for fast Bayesian estimation of multi-subject brain activation patterns

Task-evoked functional magnetic resonance imaging studies, such as the Human Connectome Project (HCP), are a powerful tool for exploring how brain activity is influenced by cognitive tasks like memory retention, decision-making, and language processing. A fast Bayesian function-on-scalar model is proposed for estimating population-level activation maps linked to the working memory task. The model is based on the canonical polyadic (CP) tensor decomposition of coefficient maps obtained for each subject. This decomposition effectively yields a tensor basis capable of extracting both common features and subject-specific features from the coefficient maps. These subject-specific features, in turn, are modeled as a function of covariates of interest using a Bayesian model that accounts for the correlation of the CP-extracted features. The dimensionality reduction achieved with the tensor basis allows for a fast MCMC estimation of population-level activation maps. This model is applied to one hundred unrelated subjects from the HCP dataset, yielding significant insights into brain signatures associated with working memory.

1 Introduction

Task-evoked functional magnetic resonance imaging (fMRI) studies, such as the Human Connectome Project (HCP), are a powerful tool for exploring how brain activity is influenced by cognitive tasks such as memory retention, decision-making, and language processing. In the HCP, a core set of tasks was devised aiming to delineate the functionally unique nodes within the human brain to establish connections between tasks and activation signatures in critical network nodes (Barch et al., 2013). These data have motivated the work presented here.

Brain imaging data are often represented in the form of a multi-dimensional array, called a tensor. These arrays can be seen as functions of space and/or time. In multi-subject studies, it is of interest to relate these functional response data to a set of scalar predictors through a function-on-scalar model. Due to the high-dimensionality of the functional response, most models in the literature will employ a form of regularization. Regularization involves borrowing strength across observations within a function. Although regularization can be achieved through various methods such as penalized regression methods or/and sparse-inducing priors, one very important approach is by the use of basis functions. Each basis function defines a linear combination among the locations within the function, establishing a framework for borrowing strength between observations. Commonly used basis functions are splines, Fourier series, wavelets, and principal components (PCs; Morris, 2015). Some specific basis functions tailored to fMRI data are also available in Miranda and Morris (2021) and in Haxby et al. (2020). The former is based on a composite-hybrid basis obtained from a principal component analysis and the latter is based on the generalized Procrustes analysis.

Some efforts have been devoted to analyzing neuroimaging as tensor data. In one line of research, the CP decomposition is used as a descriptive tool to represent multisubject and multisession fMRI data. The resulting spatial patterns or modes of the tensor representation are, in turn, validated by visual inspection. The extracted components in each spatial, temporal, runs, among other dimensions, can be further inspected to describe variability in each of these directions (Andersen and Rayens, 2004). Although this descriptive decomposition is useful, it is often of interest to relate the tensor neuroimaging data with other variables such as age, gender, and sex, disease status, through a regression analysis framework.

Previous research that tackles tensor regression are models in which the response is a scalar and the predictor is a tensor (Zhou et al., 2013; Guhaniyogi et al., 2017; Li et al., 2018; Miranda et al., 2018). Commonly in these approaches, the goal is to understand the changes of the response outcome as the tensor variable changes. Often, a tensor decomposition is implemented as tool for dimensionality reduction of tensor covariates. For example, the authors in Miranda et al. (2018) use this approach to predict Alzheimer's Disease status as a function of magnetic resonance imaging data through a CP decomposition of the imaging covariate.

Meanwhile, tensor response regression aims to study the change of the image as the predictors such as the disease status and age vary. Research in this direction often propose a low-rank tensor decomposition of the coefficients associated with covariates (Rabusseau and Kadri, 2016; Spencer et al., 2020). In a different application context, Lock and Li (2018) propose a likelihood-based latent variable tensor response regression, where the tensor data corresponds to latent variables informed by additional covariates. The authors apply their model to a 2D facial image database with facial descriptor.

Although aiming to address the same type of problem, our proposal is methodologically distinct from these previous approaches. Here, we use the CP tensor decomposition to construct a tensor basis function that not only determines the mechanisms for borrowing information across voxels but also effectively reduces the dimensionality of the brain data, which is projected onto the basis space. The projected data are in turn modeled as functions of covariates of interest using a Bayesian MCMC algorithm. This approach is inspired by the wavelet functional mixed models proposed in Morris and Carroll (2006) and further developed in Morris et al. (2011), Zhu et al. (2011, 2018), and Zhang et al. (2016). The novelty of the methodology here is that the CP decomposition captures the 3D spatial brain intricacies across different subjects and concatenation of the spatial directions does not happen as in wavelet basis or principal components. We apply this approach to a second-level analyzes (as explained in Section 2) of fMRI data from the one hundred unrelated subjects of the HCP data. The method, however, extends to a vast range of imaging modalities such as structural MRI images, and electroencephalography data.

The rest of the paper is organized as follows. In Section 2 we introduce the CP decomposition and estimation method, including a brief review of tensor concepts and properties (Section 2.1) and a review of the method used to obtain subject-level brain maps for fMRI task data (Section 2.2, and the design the CP tensor basis functions with the Bayesian modeling approach (Section 2.3). In Section 3 we apply the proposed approach to the Working Memory data from HCP to find population-level maps of a particular Working Memory contrast. In Section 4 we provide some concluding remarks and directions for future work.

2 Methods

The proposed approach involves a two-stage procedure. In the first stage, the 4D data for each subject are modeled through the Bayesian composite-hybrid basis model proposed in Miranda and Morris (2021) and described in Section 2.2. In the second stage, posterior means of contrast brain maps are modeled by the CP tensor basis model. This modeling strategy is flexible, as subject-level maps can be estimated by any known single-subject estimation technique and combined into coefficients maps for any contrasts of interest. Furthermore, these maps are assembled into a 4D tensor (3D volume and subjects) to be associated with population-level covariates, e.g., age and sex.

2.1 Background

We review a few important concepts and operations necessary for the comprehension of this manuscript. For an exhaustive review on tensor decomposition and applications (see Kolda and Bader, 2009).

2.1.1 Working with tensors

A tensor is a multidimensional array with order given by the number of its dimensions, also known as modes. For example, a matrix is a tensor of order 2. A tensor $\mathcal{X}\in {ℝ}^{I_1\times I_2\dots \times I_K}$ of order K is said to be a rank-one tensor if it can be written as an outer product of K vectors, i.e.,

$X=a^{(1)}\circ a^{(2)}\circ \dots \circ a^{(K)},$

where ° is the outer product of vectors. This is equivalent to writing each element of $x_{i_1{i}_2\dots i_K}=a_{i_1}^{(1)}{a}_{i_2}^{(2)}\dots a_{i_K}^{(K)}$, for 1 ≤ i_k ≤ i_K.

Matricization or unfolding of a tensor is the process of reordering the tensor into a matrix. We will denote the mode-n matricization of the tensor $\mathcal{X}$ above as X_(n) and arrange the mode-n dimension to be the columns of X_(n). Please see Kolda and Bader (2009) for how the observations are arranged for a general tensor order. As an example, consider the matricization of a tensor $\mathcal{X}\in {ℝ}^{2\times 3\times 4}$. The mode-3 matricization X₍₃₎ is a matrix of size 6 × 4 where, for a fixed column k in the matrix, we have elements x_11k, x_21k, x_12k, x_22k, x_13k, x_23k in the kth row.

Tensor decompositions are higher-order generalizations of the singular value decomposition and principal component analysis. Although there are various tensor decomposition methods, we focus here on the Canonical Polyadic (CP). The CP decomposition factorizes a tensor as a sum of rank-one tensors, as illustrated in Figure 1 for a tensor of order four.

Figure 1

Figure 1. CP decomposition of a 4D array into the sum of rank-one tensors. An element of the tensor on the left can be written as the sum of the products of the elements in vectors a,b,c, and d, e.g., $x_{j_1{j}_2{j}_3{j}_4}=\sum _{r=1}^R{\lambda }_r{a}_{r{j}_1}{b}_{r{j}_2}{c}_{r{j}_3}{d}_{r{j}_4}$.

Mathematically, we write the CP decomposition above as

$X\approx 〚\lambda ;A,B,C,D〛\equiv {\displaystyle \sum _{r=1}^R{\lambda }_r}{a}_r\circ b_r\circ c_r\circ d_r$

The factorization means that every element of the four-dimensional array $\mathcal{X}\in {ℝ}^{J_1\times J_2\times J_3\times J_4}$ can be written as $x_{j_1{j}_2{j}_3{j}_4}=\sum _{r=1}^R{\lambda }_r{a}_{j_{1}r}{b}_{j_{2}r}{c}_{j_{3}r}{d}_{j_{4}r}$, where the vectors $a_r\in {ℝ}^{J_1}$, $b_r\in {ℝ}^{J_2}$, $c_r\in {ℝ}^{J_3}$, and $d_r\in {ℝ}^{J_4}$ form the columns of the factor matrices A,B,C, and D.

The factor matrices refer to the combination of the vectors from the rank-one components, i.e., A = [a₁a₂…a_R] and likewise for B, C, and D. We assume that the columns of A, B, C, and D are normalized to length one, with the weights absorbed into the vector λ∈ℝ^R. The approximation refers to the estimation of the factor matrices that minimize the Frobenius norm of the difference between $\mathcal{X}$ and ⟦λ;A,B,C,D⟧ (Carroll and Chang, 1970).

2.2 Subject-level modeling

For each voxel v = 1, …, N_v in a volumetric image, the task fMRI response function is a time series of T time points based on the Blood Oxygen Level Dependent (BOLD) contrast (Ogawa et al., 1990). These measurements are an indirect and non-invasive measure of brain activity and can be modeled as a function of the external stimuli:

$\begin{array}{ll}{y}_t(v)={\displaystyle \sum _{j=1}^p}{b}_j(v)x_j(t)+e_t(v),& (1)\end{array}$

where x_j(t) is the covariate for stimulus j, representing the indicator function of stimulus j at time t convolved with the hemodynamic response of a neural event at time t (Lindquist, 2008), and e_t(v) is a measurement error. The coefficients b_j(v) characterize the relationship of brain activity during stimulus j. The primary goal of activation studies is to identify which voxels are differentially activated by specified stimuli.

There are many ways to estimate the coefficient maps b_j(v), including the popular general linear model available through the software SPM (Friston et al., 2006). The model in Equation 1 can be written in matrix form as:

$\begin{array}{l}Y=XB+E\end{array}$

To estimate the brain maps $B\in {ℝ}^{p\times N_v}$, we will use the Bayesian approach proposed in Miranda and Morris (2021), which consists of fitting a basis-space model:

$\begin{array}{l}{Y}^{*}=X^{*}{B}^{*}+E^{*}\end{array}$

This can be written in vectorized form as:

$\begin{array}{l}\text{vec}(Y^{*})=(I_S\otimes X^{*})(I_S\otimes B^{*})+\text{vec}(E^{*}),\end{array}$

with vec(𝕐^*) = (Ψ^TΦ^T ⊗ W)vec(𝕐) and vec(𝔼^*) = (Ψ^TΦ^T⊗W)vec(𝔼). In summary, the authors' tensor basis approach consists of choosing the set of bases {Ψ,Φ,W} to account for global and local spatial correlations and long-memory temporal correlations. A Bayesian MCMC is proposed and samples from the coefficient maps are obtained in the basis space, and the inverse transform is applied to obtain the MCMC samples in the voxel space. For more details, refer to Miranda and Morris (2021).

2.3 Multi-subject modeling

Let $\mathcal{Y}$ be the 4D tensor representing each coefficient map or contrast map for all N subjects. Then $\mathcal{Y}\in {ℝ}^{p_1\times p_2\times p_3\times N}$. We can model $\mathcal{Y}$ as a function of p subject-level covariates Z∈ℝ^p as follows:

$\begin{array}{ll}Y=Z\gamma +E,& (2)\end{array}$

where Y is a N × N_v matrix representing the mode four unfolding of the tensor $\mathcal{Y}$ with N_v = p₁ × p₂ × p₃ being the number of voxels; Z is a N×p design matrix of p subject-level covariates; γ is a p × N_v matrix of coefficients; and E is a N×N_v matrix of random errors.

Properly modeling the data as in Equation 2 is challenging due to its very large number of voxels, typically around a million, and the spatial dependency inherited from the brain volume data. Most Bayesian models rely on either downsampling the volume data to reduce the number of voxels or rely on simplified assumptions on the prior of γ and/or E. The solution proposed here involves projecting the data into a basis space through a tensor basis obtained from the CP decomposition of the tensor $\mathcal{Y}$. This decomposition is low rank and effectively extracts brain volume features that are common to all subjects. By extracting the common features, the remaining information is subject-specific and can be then related to the covariates of interest.

The CP decomposition of $\mathcal{Y}$ is given by:

$\begin{array}{ll}Y\approx {\displaystyle \sum _{r=1}^R}{\lambda }_r{a}_r^{(1)}◦a_r^{(2)}◦a_r^{(3)}◦g_r\equiv ⟦\lambda ;A^{(1)},A^{(2)},A^{(3)},G⟧,& (3)\end{array}$

where the vectors $a_r^{(1)}\in {ℝ}^{p_1}$, $a_r^{(2)}\in {ℝ}^{p_2}$, $a_r^{(3)}\in {ℝ}^{p_3}$, and $g_r\in {ℝ}^N$. The matrices A⁽¹⁾, A⁽²⁾, A⁽³⁾, and G are called factor matrices and their columns are associated with these vectors. From Kolda and Bader (2009), we can write Y from model (Equation 2) as:

$\begin{array}{l}Y=G\Lambda {(A^{(3)}\odot A^{(2)}\odot A^{(1)})}^T+ϵ,\end{array}$

where G, A⁽¹⁾, A⁽²⁾, and A⁽³⁾ are estimated to minimize the Frobenius norm of ϵ (Kolda and Bader, 2009). Let $L=\Lambda {(A^{(3)}\odot A^{(2)}\odot A^{(1)})}^T\in {ℝ}^{N_v\times R}$, then P =L^†, the Moore-Penrose pseudoinverse of L, provides a projection matrix formed by spatial basis that can be used to effectively reduce data from N_v to R spatial locations.

Based on the projection matrix, we multiply both sides of Equation 2 by P to obtain the basis space model:

$\begin{array}{ll}G=Z{\gamma }^{*}+E^{*},& (4)\end{array}$

where γ^* =γP is the p × R matrix of coefficients in the basis space, and E^* =EP is a N × R matrix of random errors.

2.3.1 Basis-space Bayesian modeling

The columns of G are correlated due to the nature of the CP decomposition, and model assumptions should reflect that. For each row of E^*, we assume $E_i^{*}~N_R(0,{\Sigma }_{ϵ})$ and priors:

$\begin{array}{l}p({\gamma }^{*},{\Sigma }_{ϵ})=p({\Sigma }_{ϵ})p({\gamma }^{*}|{\Sigma }_{ϵ}),\end{array}$

with ${\gamma }^{*}~\mathrm{\text{MN}}(g_0,L_0^{-1},{\Sigma }_{ϵ})$ and ${\Sigma }_{ϵ}~{\mathrm{\text{W}}}^{-1}(V_0,{\nu }_0)$, where MN indicates the matrix normal distribution and W⁻¹ indicates the Inverse Wishart distribution. The joint posterior distribution of ${\gamma }^{*},{\Sigma }_{ϵ}|G,Z$ factorizes as a product of an Inverse Wishart and a Matrix Normal distributions as follows:

$\begin{array}{c}{\Sigma }_{ϵ}|G,Z~{\text{W}}^{-1}(V_n,{\nu }_n),\\ {\gamma }^{*}|G,Z,{\Sigma }_{ϵ}~\text{MN}(g_n,L_n^{-1},{\Sigma }_{ϵ}),\end{array}$

with:

$\begin{array}{c}{V}_n=V_0+{(Y-Z{g}_n)}^T(Y-Z{g}_n)+{(g_n-g_0)}^T{L}_0(g_n-g_0),\\ {\nu }_n={\nu }_0+N,\\ g_n={(Z^{T}Z+L_0)}^{-1}(Z^{T}G+L_0{g}_0),\\ L_n=Z^{T}Z+L_0.\end{array}$

Sampling from these distributions is straightforward and done through Gibbs sampling.

2.3.2 Rank selection

To choose the appropriate rank for the CP decomposition, we propose a 10-fold cross-validation procedure with the following steps for each of the training/test set combination:

1. Obtain the factor matrices of the CP decomposition for the training data,

2. Proceed to estimate γ^* from model Equation 4 following the MCMC algorithm detailed in Section 2.3.1,

3. Using $\widehat{{\gamma }^{*}}$ from the previous step and Z_Test, estimate $G_{\mathrm{\text{Test}}}=Z_{\mathrm{\text{Test}}}\widehat{{\gamma }^{*}}$,

4. Obtain ${\hat{\mathcal{Y}}}_{\mathrm{\text{Test}}}=⟦\lambda ;A^{(1)},A^{(2)},A^{(3)},G_{\mathrm{\text{Test}}}⟧$ as in Equation 3,

5. Choose the rank that minimizes the Frobenius norm of the average of the ${\mathcal{Y}}_{\mathrm{\text{Test}}}-{\hat{\mathcal{Y}}}_{\mathrm{\text{Test}}}$, obtained across the training/test sets.

2.3.3 Inference on the brain space

For every MCMC sample of γ^*, we obtain γ =γ^*L in the original voxel space. To control for the experimentwise error rate in the voxel space, joint credible bands as in Ruppert et al. (2003) can be calculated. For each voxel v, a joint 100(1−α)% credible interval for any group-level contrast C is given by

$\begin{array}{l}{I}_{\alpha }(v)=Ĉ(v)\pm q_{(1-\alpha )}\left[\widehat{\text{std}}(Ĉ(v))\right],\end{array}$

where the quantile q_(1−α) is obtained from z⁽¹⁾, …, z^(M) with $z^{(m)}=\underset{v\in \mathcal{V}}{max}\left|\frac{C^{(m)}(v)-Ĉ(v)}{\widehat{\mathrm{\text{std}}}(Ĉ(v))}\right|$, for m = 1, …, M and M is the total number of MCMC samples after burn-in and thinning (Miranda and Morris, 2021).

These joint credible bands account for multiple testing across the voxel space and can be inverted to flag activation signatures based on voxels for which the 100(1−α)% joint credible bands exclude zero (Meyer et al., 2015; Miranda and Morris, 2021). Denote by P_SimBaS the minimum α at which each interval excludes zero (P_SimBaS = min{α:0∉I_α(v)}). P_SimBaS can be directly computed by

$\begin{array}{ll}{P}_{SimBaS}(v)=\frac{1}{M}{\displaystyle \sum _{m=1}^M}1\left\{\left|\frac{Ĉ(v)}{\widehat{\text{std}}(Ĉ(v))}\right|\le z^{(m)}\right\}.& (5)\end{array}$

These are probabilities of the simultaneous band scores. For an specific value of α, we can identify v for which P_SimBaS(v) < α as significant. This is equivalent to checking if the joint credible intervals, I_α(v), cover zero at a specific α-level (Meyer et al., 2015).

3 Application to the HCP data

The proposed model is applied to the 100 unrelated subjects data from the Working Memory task of the Human Connectome Project. The volumes considered were collected from the right-left phase. The experiment consists of a total of eight blocks, half corresponding to the 2-back task and the other half the 0-back task. In each block, the participants were shown a series of 10 images. Images corresponding to 4 stimuli (tools, places, body parts, and faces) were embedded within the memory task. Each image was shown for 2.5 s, followed by a 15 s inter-block interval. fMRI volumes were obtained every 720 ms (TR). Each volume consisted of images of size 91 × 109 × 91 for a total of 405 time frames.

3.1 Subject-level analysis

We follow the analysis described in Section 2.2. To obtain the local spatial basis, we partitioned the brain into ROIs using a digitalized version of the original Talairach structural labeling that was registered into the MNI152 space. The atlas can be obtained from the FSL atlas library with information found at: https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Atlases (Jenkinson et al., 2012). We considered ROIs that had at least 125 voxels and obtained a total of 298 regions. Within each region, the local basis Φ was formed by the eigenvectors of the SVD decomposition for the matrix data in each ROI as described in Miranda and Morris (2021). Next, the global spatial basis Ψ was obtained by first projecting the brain data using the local basis and then calculating the SVD decomposition of the projected data, exactly as in Miranda and Morris (2021). The temporal basis W was chosen to be the Haar wavelet. The estimation was performed at the cluster Cedar managed by the Digital Research Alliance of Canada (former Compute Canada). The jobs were submitted in parallel for the 100 subjects and results took an hour to finish.

For each subject, after estimation, we constructed a contrast map of the 2 vs. 0 back working memory task. The maps were obtained from the posterior means of the coefficients, and assembled into a 4D array $\mathcal{Y}\in {ℝ}^{91\times 109\times 91\times 405}$.

3.2 Multi-subject analysis

We take the tensor $\mathcal{Y}$ containing the posterior means of the contrast 2 vs. 0 back for the maps obtained in the single-level analysis. Next, we follow the method proposed in Section 2.3 to estimate the coefficients maps γ for the mean of 2 vs. 0 back and the covariate Sex. We are interested in finding common activation maps for all subjects and aim to investigate if there are differences in working memory task based on sex. Next, we obtain the CP decomposition of the tensor $\mathcal{Y}$, taking R = 30 selected based on the cross-validation procedure described in Section 2.3.2. Figure 2 shows the Frobenius norm as a function of the rank, for rank values R = 10, 15, 20, 25, …, 80. Following the decomposition, we aim to identify significant locations in the brain for the Working Memory task 2 vs. 0 back, we proceed as follows.

(i) Obtain MCMC samples for γ^*, as detailed in Section 2.3.1.

(ii) Use the inverse transform of P, $L=\Lambda {(A^{(1)}\odot A^{(2)}\odot A^{(3)})}^T\in {ℝ}^{N_v\times R}$, to obtain MCMC samples for γ =γ^*L in the voxel space.

(iii) For each subject-level covariate, use the MCMC samples of γ to obtain P_SimBas as described in Equation 5 to control for the experimentwise error rate on multiple testing.

(iv) For P_SimBas < α = 0.01 flag the voxel as different than zero, indicating activation.

(v) Obtain clusters of significant voxels that are at least 125 voxels large.

Figure 2

Figure 2. Cross-validation to choose the rank R. We want R to minimize the Frobenius norm.

The computational time to run this analysis is 1.62 minutes on a MacBook Pro 2.8 GHz Quad-Core Intel Core i7 with 16 GB 2,133 MHz LPDDR3 memory. For the cross-validation procedure, it took 38 minutes for a grid of 15 rank values. For faster computational time on rank selection, one could run a 5-fold cross validation instead and/or a less fine grid of ranks.

We consider two maps of coefficients, the first is the intercept, which is the average activation for the contrast 2 vs. 0 back of the Working Memory task. The second is the map associated with the variable sex. We did not find any significant differences in the working memory contrast between females and males. Figure 2 shows the results for the intercept revealing 4,445 active voxels. The larger cluster in red shows 1,268 voxels encompassing 24 distinct brain regions as subdivided by the Tailarach Atlas, (labels can be found at: http://www.talairach.org/about.html#Labels). Within this cluster, larger active regions are in the right parietal lobe, including inferior, superior, precuneous parietal regions. The cluster in green is formed by 1182 voxels, with active regions on the left parietal lobe, including the inferior, superior, precuneous, and cuneous regions. The parietal lobe has been previously identified as playing a key role in working memory (Berryhill and Olson, 2008; Koenigs et al., 2009; Lee et al., 2022).

The next cluster in blue shows the right frontal lobe, including the superior and middle frontal gyrus. The cluster in yellow shows 345 voxels in the right frontal lobe, including the superior and middle frontal gyrus. The cluster in light orange shows 276 voxels on the left frontal lobe in the middle frontal gyrus. These regions have also been identified as playing a critical role in working memory (Owen et al., 1998). The cluster in burgundy indicates 313 voxels on the right cerebellum at the posterior lobe region. A smaller cluster of 139 voxels was also identified on the posterior lobe at the left cerebellum. The authors in Deverett et al. (2019) found that the cerebellum can influence the accurate maintenance of working memory.

Next, we compare the results obtained from the Bayesian CP tensor basis approach with the results obtained from the Statistical Parametric Map (SPM) software (Friston et al., 2006). We ran both the first and the second-level analysis in SPM with parameter settings to reflect no spatial normalization and an AR(1) noise component. The results for the Average of 2 vs. 0 back contrast are displayed in Figure 3. The two largest clusters were sizes 64 and 58 voxels, respectively. Overlaid onto the template is a mask with all significant voxels based on a voxel-level threshold that controls for the family-wise error rate of 0.05. Inspecting Figure 4 reveals that this SPM voxelwise model only capture small parts of the parietal and frontal lobe and has significant less power than the CP tensor basis approach.

Figure 3

Figure 3. Working memory contrast of 2 vs. 0 back.

Figure 4

Figure 4. Working memory contrast of 2 vs. 0 back obtained from SPM.

4 Final remarks

In this paper, we propose a fast Bayesian model that takes tensor responses as a function of population-level covariates. The model is based on using the CP tensor decomposition approach to construct a basis function with the goal of modeling the spatial structure of the tensor data. This approach allows us to borrow information across voxels, increasing the power for the statistical tests performed. Because the estimation of the modeling components is performed in the basis space, we have a very fast Bayesian procedure. The Bayesian MCMC estimation method for the second-level analysis, including cross-validation for rank selection, estimation, and inference, takes 40 min on a regular laptop. This time can be much faster if using a cluster with more powerful computing capabilities.

We apply this approach to a second-level analysis of fMRI data from one hundred unrelated subjects of the HCP dataset. The results obtained are consistent with findings in the literature showing the role of the parietal lobe, the middle frontal gyrus, and the posterior lobe of the cerebellum in working memory.

One limitation of the proposed approach is that CP decomposition applies only to volumetric data. For more modern imaging files, such as the CIFTI files available in the Human Connectome Project, the CP decomposition would need to aggregate spatial information as one large dimension. In this case, the benefits of the CP decomposition are unclear as the CP decomposition might be similar to other basis approaches, including principal components and wavelets. The authors plan to conduct an investigation on that front.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found at: https://db.humanconnectome.org/.

Ethics statement

The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required from the participants or the participants' legal guardians/next of kin in accordance with the national legislation and institutional requirements.

Author contributions

MM: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing - original draft, Writing - review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was funded by Natural Sciences and Engineering Research Council of Canada, grant number: RGPIN-2020-06941.

Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Andersen, A. H., and Rayens, W. S. (2004). Structure-seeking multilinear methods for the analysis of fMRI data. NeuroImage 22, 728–739. doi: 10.1016/j.neuroimage.2004.02.026

PubMed Abstract | Crossref Full Text | Google Scholar

Barch, D. M., Burgess, G. C., Harms, M. P., Petersen, S. E., Schlaggar, B. L., Corbetta, M., et al. (2013). Function in the human connectome: task-fMRI and individual differences in behavior. NeuroImage 80, 169–189. doi: 10.1016/j.neuroimage.2013.05.033

PubMed Abstract | Crossref Full Text | Google Scholar

Berryhill, M. E., and Olson, I. R. (2008). The right parietal lobe is critical for visual working memory. Neuropsychologia 46, 1767–1774. doi: 10.1016/j.neuropsychologia.2008.01.009

PubMed Abstract | Crossref Full Text | Google Scholar

Carroll, J. D., and Chang, J.-O. (1970). Analysis of individual differences in multidimensional scaling via an n-way generalization of “eckart-young” decomposition. Psychometrika 35, 283–319.

Google Scholar

Deverett, B., Kislin, M., Tank, D. W., and Wang, S. S. (2019). Cerebellar disruption impairs working memory during evidence accumulation. Nat. Commun. 10:3128. doi: 10.1038/s41467-019-11050-x

PubMed Abstract | Crossref Full Text | Google Scholar

Friston, K. J., Penny, W. D., Ashburner, J., Kiebel, S. J., and Nichols, T. E. (2006). Statistical Parametric Mapping: The Analysis of Functional Brain Images. London: Academic Press.

Google Scholar

Guhaniyogi, R., Qamar, S., and Dunson, D. B. (2017). Bayesian tensor regression. J. Machine Learn. Res. 18, 1–31. Available at: https://jmlr.org/papers/v18/16-362.html

Google Scholar

Haxby, J. V., Guntupalli, J. S., Nastase, S. A., and Feilong, M. (2020). Hyperalignment: modeling shared information encoded in idiosyncratic cortical topographies. eLife 9:e56601. doi: 10.7554/eLife.56601

PubMed Abstract | Crossref Full Text | Google Scholar

Jenkinson, M., Beckmann, C. F., Behrens, T. E., Woolrich, M. W., and Smith, S. M. (2012). Fsl. NeuroImage 62, 782–790. doi: 10.1016/j.neuroimage.2011.09.015

PubMed Abstract | Crossref Full Text | Google Scholar

Koenigs, M., Barbey, A. K., Postle, B. R., and Grafman, J. (2009). Superior parietal cortex is critical for the manipulation of information in working memory. J. Neurosci. 29, 14980–14986. doi: 10.1523/JNEUROSCI.3706-09.2009

PubMed Abstract | Crossref Full Text | Google Scholar

Kolda, T. G., and Bader, B. W. (2009). Tensor decompositions and applications. SIAM Rev. 51, 455–500. doi: 10.1137/07070111X

Crossref Full Text | Google Scholar

Lee, B., Kim, J. S., and Chung, C. K. (2022). Parietal and medial temporal lobe interactions in working memory goal-directed behavior. Cortex 150, 126–136. doi: 10.1016/j.cortex.2022.02.005

PubMed Abstract | Crossref Full Text | Google Scholar

Li, X., Xu, D., Zhou, H., and Li, L. (2018). Tucker tensor regression and neuroimaging analysis. Stat. Biosci. 10, 520–545. doi: 10.1007/s12561-018-9215-6

PubMed Abstract | Crossref Full Text | Google Scholar

Lindquist, M. A. (2008). The statistical analysis of fMRI data. Stat. Sci. 23, 439–464. doi: 10.1214/09-STS282

Crossref Full Text | Google Scholar

Lock, E. F., and Li, G. (2018). Supervised multiway factorization. Electr. J. Stat. 12, 1150–1180. doi: 10.1214/18-EJS1421

PubMed Abstract | Crossref Full Text | Google Scholar

Meyer, M. J., Coull, B. A., Versace, F., Cinciripini, P., and Morris, J. S. (2015). Bayesian function-on-function regression for multilevel functional data. Biometrics 71, 563–574. doi: 10.1111/biom.12299

PubMed Abstract | Crossref Full Text | Google Scholar

Miranda, M. F., and Morris, J. S. (2021). Novel bayesian method for simultaneous detection of activation signatures and background connectivity for task fMRI data. arXiv [preprint] arXiv:2109.00160. doi: 10.48550/arXiv.2109.00160

Crossref Full Text | Google Scholar

Miranda, M. F., Zhu, H., and Ibrahim, J. G. (2018). TPRM: tensor partition regression models with applications in imaging biomarker detection. Ann. Appl. Stat. 12, 1422–1450. doi: 10.1214/17-AOAS1116

PubMed Abstract | Crossref Full Text | Google Scholar

Morris, J. S. (2015). Functional regression. Ann. Rev. Stat. Its Appl. 2, 321–359. doi: 10.1146/annurev-statistics-010814-020413

Crossref Full Text | Google Scholar

Morris, J. S., Baladandayuthapani, V., Herrick, R. C., Sanna, P., and Gutstein, H. (2011). Automated analysis of quantitative image data using isomorphic functional mixed models, with application to proteomics data. Ann. Appl. Stat. 5, 605–1125. doi: 10.1214/10-aoas407

PubMed Abstract | Crossref Full Text | Google Scholar

Morris, J. S., and Carroll, R. J. (2006). Wavelet-based functional mixed models. J. Royal Stat. Soc. Ser. B 68, 179–199. doi: 10.1111/j.1467-9868.2006.00539.x

PubMed Abstract | Crossref Full Text | Google Scholar

Ogawa, S., Lee, T. M., Kay, A. R., and Tank, D. W. (1990). Brain magnetic resonance imaging with contrast dependent on blood oxygenation. Proce. Natl. Acad. Sci. U. S. A. 87, 9868–9872.

Google Scholar

Owen, A. M., Stern, C. E., Look, R. B., Tracey, I., Rosen, B. R., and Petrides, M. (1998). Functional organization of spatial and nonspatial working memory processing within the human lateral frontal cortex. Proc. Natl. Acad. Sci. U. S. A. 95, 7721–7726.

Google Scholar

Rabusseau, G., and Kadri, H. (2016). “Low-rank regression with tensor responses,” in Advances in Neural Information Processing Systems, Vol. 29, eds. D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett (Curran Associates, Inc.), 1867–1875.

Google Scholar

Ruppert, D., Wand, M. P., and Carroll, R. J. (2003). “Semiparametric Regression,” in Cambridge Series in Statistical and Probabilistic Mathematics (Cambridge: Cambridge University Press).

Google Scholar

Spencer, D., Guhaniyogi, R., and Prado, R. (2020). Joint bayesian estimation of voxel activation and inter-regional connectivity in fMRI experiments. Psychometrika 85, 845–869. doi: 10.1007/s11336-020-09727-0

PubMed Abstract | Crossref Full Text | Google Scholar

Zhang, L., Baladandayuthapani, V., Zhu, H., Baggerly, K. A., Majewski, T., Czerniak, B. A., et al. (2016). Functional car models for large spatially correlated functional datasets. J. Am. Stat. Assoc. 111, 772–786. doi: 10.1080/01621459.2015.1042581

PubMed Abstract | Crossref Full Text | Google Scholar

Zhou, H., Li, L., and Zhu, H. (2013). Tensor regression with applications in neuroimaging data analysis. J. Am. Stat. Assoc. 108, 540–552. doi: 10.1080/01621459.2013.776499

PubMed Abstract | Crossref Full Text | Google Scholar

Zhu, H., Brown, P. J., and Morris, J. S. (2011). Robust, adaptive functional regression in functional mixed model framework. J. Am. Stat. Assoc. 106, 1167–1179. doi: 10.1198/jasa.2011.tm10370

PubMed Abstract | Crossref Full Text | Google Scholar

Zhu, H., Versace, F., Cinciripini, P. M., Rausch, P., and Morris, J. S. (2018). Robust and Gaussian spatial functional regression models for analysis of event-related potentials. NeuroImage 181, 501–512. doi: 10.1016/j.neuroimage.2018.07.006

PubMed Abstract | Crossref Full Text | Google Scholar