Discovering hidden brain network responses to naturalistic stimuli via tensor component analysis of multi-subject fMRI data

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY-NC-ND 4.0 https://creativecommons.org/licenses/by-nc-nd/4.0/ Discovering hidden brain network responses to naturalistic stimuli via tensor component analysis of multi-subject fMRI data © 2022 the Authors Published version Hu, Guoqiang; Li, Huanjie; Zhao, Wei; Hao, Yuxing; Bai, Zonglei; Nickerson, Lisa D.; Cong, Fengyu Hu, G., Li, H., Zhao, W., Hao, Y., Bai, Z., Nickerson, L. D., & Cong, F. (2022). Discovering hidden brain network responses to naturalistic stimuli via tensor component analysis of multi-subject fMRI data. Neuroimage, 255, Article 119193. https://doi.org/10.1016/j.neuroimage.2022.119193 2022 NeuroImage 255 (2022) 119193 Contents lists available at ScienceDirect NeuroImage journal homepage: www.elsevier.com/locate/neuroimage Discovering hidden brain network responses to naturalistic stimuli via tensor component analysis of multi-subject fMRI data Guoqiang Hu a , Huanjie Li a , Wei Zhao a , Yuxing Hao a , Zonglei Bai b , Lisa D. Nickerson c , d , ∗ , Fengyu Cong a , e , f , g , ∗ a School of Biomedical Engineering, Faculty of Electronic Information and Electrical Engineering, Dalian University of Technology, Dalian, China b School of Electronics Engineering and Computer Science, Peking University, Beijing, China c Brain Imaging Center, Mclean Hospital, Belmont, MA, USA d Department of Psychiatry, Harvard Medical School, Boston, MA, USA e School of Artificial Intelligence, Faculty of Electronic Information and Electrical Engineering, Dalian University of Technology, Dalian, China f Key Laboratory of Integrated Circuit and Biomedical Electronic System, Liaoning Province. Dalian University of Technology, Dalian, China g Faculty of Information Technology, University of Jyvaskyla, Jyvaskyla, Finland a r t i c l e i n f o Keywords: Tensor components analysis Naturalistic stimuli fMRI Inter-subject correlation a b s t r a c t The study of brain network interactions during naturalistic stimuli facilitates a deeper understanding of human brain function. To estimate large-scale brain networks evoked with naturalistic stimuli, a tensor component analysis (TCA) based framework was used to characterize shared spatio-temporal patterns across subjects in a purely data-driven manner. In this framework, a third-order tensor is constructed from the timeseries extracted from all brain regions from a given parcellation, for all participants, with modes of the tensor corresponding to spatial distribution, time series and participants. TCA then reveals spatially and temporally shared components, i.e., evoked networks with the naturalistic stimuli, their time courses of activity and subject loadings of each component. To enhance the reproducibility of the estimation with the adaptive TCA algorithm, a novel spectral clustering method, tensor spectral clustering, was proposed and applied to evaluate the stability of the TCA algorithm. We demonstrated the effectiveness of the proposed framework via simulations and real fMRI data collected during a motor task with a traditional fMRI study design. We also applied the proposed framework to fMRI data collected during passive movie watching to illustrate how reproducible brain networks are evoked by naturalistic movie viewing. 1. Introduction There is growing interest in studying brain function in response to naturalistic stimuli, for example, viewing film clips or listening to spoken narratives or music, as naturalistic experimental paradigms may evoke human cognition and behavior that more closely resembles “real-world ”brain function ( Hasson et al., 2004 ; Huth et al., 2016 ; Meer et al., 2020 ; Nishimoto et al., 2011 ; Sonkusare et al., 2019 ; Spiers and Maguire, 2007 ). Naturalistic stimulus paradigms during functional magnetic resonance imaging (fMRI) are emerging as a powerful tool to define brain imaging-based markers of psychiatric illness ( Eickhoff et al., 2020 ), with several advantages in comparison to unconstrained resting state. Namely, studying brain network function during naturalistic stimuli may facilitate a deeper understanding of human brain function since the passive state is better constrained, and partic- ∗ Corresponding authors. E-mail addresses: [email protected] (G. Hu), [email protected] (L.D. Nickerson), [email protected] (F. Cong) . ipant motion is reduced relative to unconstrained rest, which greatly increases the quality of the fMRI data. However, new analytical strategies are needed that assess both the shared rapid temporally evolving brain responses evoked by the naturalistic stimuli in participants, as well as idiosyncratic evoked signals in individual participants ( Simony and Chang, 2020 ). Both categorical and dimensional sources of variability can contribute to inter-subject variation in responses to naturalistic stimuli fMRI. In addition, while naturalistic stimuli paradigms provide better constraint of brain activity, there are challenges with modeling the evoked responses. Namely, evoked brain activity using conventional fMRI study designs and stimuli is relatively straightforward to model, whereas naturalistic stimuli are complex and dynamic, and it is much more difficult to generate a model of evoked activity for analyses. Datadriven methods that place no assumptions on the temporal course or https://doi.org/10.1016/j.neuroimage.2022.119193 . Received 31 May 2021; Received in revised form 23 February 2022; Accepted 6 April 2022 Available online 8 April 2022. 1053-8119/© 2022 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) G. Hu, H. Li, W. Zhao et al. NeuroImage 255 (2022) 119193 spatial pattern of brain activity obviate these challenges. Inter-subject correlation (ISC, Hasson et al., 2004 ) is one such data-driven approach for characterizing the consistency of brain responses across participants viewing the same naturalistic stimuli. For dynamic complex stimuli such as movies, ISC measures share information across brains by using each individual’s measured brain activity to model another individual’s brain activity. Using this strategy, the shared brain regions that respond to the same time-locked naturalistic stimuli across subjects can be estimated, even with stimuli that reflect complex dynamic real-life contexts ( Hasson et al., 2004 ; Kauppi et al., 2014 ; Lerner et al., 2011 ; Nastase et al., 2019 ). Modifications to ISC include the temporal inter-subject functional correlation (ISFC), which considers the correlations between time courses from all possible pair-wise combinations of brain parcels across subjects ( Simony et al., 2016 ), and the spatial ISC, which is an extension of temporal ISC to multi-voxel pattern analysis ( Haxby et al., 2014 ; Norman et al., 2006 ). Inter-subject representational similarity analysis (IS-RSA) is another technique that can be used to explore between-subject variability in the relationships between brain activity and behaviors ( Finn et al., 2020 ; Kriegeskorte et al., 2008 ; Mantel, 1967 ; Meer et al., 2020 ). Van der Meer and colleagues (2020) explored the differences between movie viewing and resting state with Hidden Markov Models (HMM) and found that subject differences in brain state dynamics were linked to subjective movie ratings using IS-RSA. Through the study of between-subject variability, different schematic events ( Baldassano et al., 2018 ) and different conditions when participants recall a movie ( Chen et al., 2017 , 2016 ) can be distinguished via the corresponding brain activity. In most cases, understanding patterns of brain states that are consistent across subjects as well as patterns that reflect inter-subject variability are of interest. In this study, we leverage a popular tensor decomposition method to simultaneously estimate spatial-temporal brain activity patterns that are shared across participants and that reflect inter-subject variability. Generally, ISC-based methods are implemented with either a leaveone-out framework, in which one subject’s time course is correlated with the average of all other subjects for each region, or a pairwise framework, in which correlation analysis is performed between each possible pair of subjects ( Finn et al., 2020 ). A limitation of this computational procedure is that the resulting correlations are highly interdependent and violate the assumption of common parametric tests ( Nastase et al., 2019 ), requiring careful attention to the inference method. In order to mitigate this limitation, we use a tensor component analysis (TCA) framework that characterizes spatio-temporal patterns that are shared across subjects as well as idiosyncratic features of evoked activity to naturalistic stimuli unique to different subjects, in a purely data-driven manner that assesses all brain networks simultaneously. Tensor Component Analysis (TCA), also known as tensor Canonical Polyadic Decomposition (CPD, Kolda and Bader, 2009 ), is a fundamental model for tensor decomposition of multidimensional data with more than two dimensions. Assuming the data meet the assumption of a mixture model, in which signal sources undergo a linear mixing process, TCA more accurately estimates sources than matrix decomposition algorithms, and without any constraints ( Williams et al., 2018 ). FMRI signals are innately multidimensional and can be naturally represented in tensor form. For example, a third-order fMRI tensor is organized as space × time ×subjects (the order of modes does not impact the estimation). In TCA of the fMRI tensor, the spatial and temporal information regarding brain network activity evoked by different stimuli that is common to all participants exists in the first two dimensions. Subject loadings that capture between-subject variability exist in the third dimension. TCA has demonstrated promise in a range of neuroimaging applications. It has been shown to have superior performance in identifying hidden signal sources when compared with 2-D matrix decomposition, e.g. principle component analysis (PCA) and independent component analysis (ICA), of multidimensional data ( Williams et al., 2018 ). TCA has also been explored for magnetoencephalography (MEG) data analysis ( Zhu et al., 2020 a), and nonnegative constraint TCA applied to electroencephalography (EEG) time-frequency domain data was able to identify event-related ( Cong et al., 2015a , 2015b ; Wang et al., 2018 ) and naturalistic stimulus-evoked EEG responses ( Zhu et al., 2020 b). Mokhtari et al. (2019) investigated how different tensor organization and tensor decomposition methods applied to fMRI data impact the interpretation of dynamic functional connectivity. In our previous study ( Hu et al., 2021 ), sparse constrained nonnegative TCA was proposed to estimate frequency specific coactivation patterns. TCA has also been applied to task fMRI data to explore additional dimensions of the data other than space and time, such as run and task condition ( Andersen and Rayens, 2004 ), and has been adapted for multi-subject fMRI data analysis by placing spatial and temporal constraints to address inter-subject variability ( Beckmann and Smith, 2005 ; Helwig and Hong, 2013 ; Kuang et al., 2020 , 2015 ; Mørup et al., 2008 ; Zhou and Cichocki, 2012 ). We advance TCA for analysis of multi-subject fMRI data collected during naturalistic stimuli viewing by proposing a pipeline that does not place any constraints on the data and that enhances the reproducibility of the results. We address three key issues that have limited the use of tensor decomposition for naturalistic stimuli fMRI ( Wolf et al., 2010 ). First, the TCA algorithm may be slow or fail to converge, or have suboptimal convergence when applied to whole brain data. To mitigate this issue, the TCA is usually constrained in some way ( Beckmann and Smith, 2005 ; Zhou et al., 2014 ). Instead of constraining the TCA algorithm, we propose instead to implement a parcellation strategy to reduce the data prior to TCA to facilitate convergence. Second, we propose a novel tensor spectral clustering method to enhance the reproducibility of the estimated components. Last, model order selection is always a challenge for tensor and matrix decomposition methods ( Abou-Elseoud et al., 2010 ; Beckmann, 2012 ; Kuang et al., 2018 ). Here we propose to select the model order based on component reproducibility assessed via our novel spectral clustering method. The effectiveness of the proposed framework is first demonstrated with simulated and traditional task fMRI, in which we know the ground truth stimulation time courses (and hence have a model of the brain activity). We then apply the proposed framework to fMRI data collected during movie watching, in which there is no a priori model of brain activity, to identify spatial brain networks engaged during the task. Analysis of both task fMRI and naturalistic stimuli fMRI results show that the proposed method has several advantages compared with the widely used ISC method. The rest of the paper is structured as follows. In section 2 , Materials and Methods, we present the TCA model that was used to estimate spatio-temporal shared components, e.g., the decomposition algorithm used in this paper, the criteria that were used to evaluate the reproducibility of estimated components and the model order selection method, and the test datasets. In section 3 , Results, we show results using simulated data and motor task fMRI. After establishing the robustness of the proposed framework from simulations and conventional task fMRI, we show results from application of our analytic approach to two different naturalistic stimuli fMRI datasets: one in which participants watched a short 82 second montage of scenes, collected from 184 participants by the Human Connectome Project (HCP), and one in which participants watched a longer movie (20 minutes), collected from 17 participants by Meer et al. (2020) . In sections 4 and 5 , the Discussion and Conclusions, respectively, we discuss the benefits and pitfalls of our proposed approach and the conclusions from our work. The mathematical foundation of tensor spectral clustering used to evaluate the stability of estimated components is provided in the Appendix. 2. Materials and methods 2.1. TCA model Our proposed application of TCA of naturalistic stimuli fMRI data proceeds in three steps as shown in Fig. 1 . First, a data reduction step is introduced in which the activity of node (or region) time courses is ex2 G. Hu, H. Li, W. Zhao et al. NeuroImage 255 (2022) 119193 Fig. 1. Tensor component analysis pipeline. (A) Estimating node timeseries with dual regression from standard preprocessed dense data. (B) Stacking all subjects’ node data to construct a tensor. (C) Extracting tensor components with tensor component analysis. The estimated components of each mode  𝒔 ,  𝒄 ,  𝝈correspond to the spatial distribution, time course and subject loadings, respectively. tracted from the naturalistic stimuli fMRI data using an open-access parcellation scheme derived from independent component analysis (ICA) of resting state fMRI; next, the time courses from all participants are stacked to construct a third-order fMRI tensor; last, TCA is applied to estimate spatio-temporal patterns and their corresponding subject loadings. In the first step, rather than applying TCA to voxel-wise fMRI data, a whole brain parcellation scheme is applied to extract node time courses for the TCA. This is done for several reasons. First, BOLD signals are expected to be correlated across neighboring voxels and, with an appropriate parcellation method, node time courses will have a higher signal noise ratio (SNR) compared with the SNR of time courses from dense data ( Glasser et al., 2016 ). Parcellation schemes based on ICA of fMRI data with high model orders ( > 100 to several hundred) will be comprised of components that feature individual small brain regions, bilateral brain regions, or sparse sub-networks that may have regions overlapping other components (e.g., reflecting hubs such as posterior cingulate cortex), and can thus be considered as nodes for use in network analysis ( Smith, 2012 ). Several studies have demonstrated that brain parcellation with spatial ICA demonstrates better performance for network modeling compared with other parcellation methods (Smith et al., 2011; Arslan et al., 2018), with higher model orders providing better performance ( Pervaiz et al., 2020 ). In this study, a brain parcellation scheme derived from ICA of the Human Connectome Project (HCP, Van Essen et al., (2013) ) resting state fMRI data with model order of 300 (provided by the HCP, Smith et al., 2011) is used to extract node time courses for TCA. The time courses are extracted via the first stage of dual regression ( Nickerson et al., 2017 ) in which the 3 G. Hu, H. Li, W. Zhao et al. NeuroImage 255 (2022) 119193 full set of ICA components, 𝐒 ICA , are regressed against each participant’s 4D fMRI data (e.g., multivariate spatial regression) to extract the time courses. This is different from how conventional binary parcellation masks are used to extract average time courses from each node in that the multivariate spatial regression accurately handles any potential spatial overlap among the ICA maps due to individual brain regions participating in more than one brain network (or sub-network) represented in the ICA maps. The analysis code for the tensor component analysis framework is available at https://github.com/GHu-DUT/ TensorcomponentsanalysisfornaturalisticstimulifMRI . 2.1.1. Multilinear mixing model For naturalistic fMRI we assume, similar to the ISC model ( Finn et al., 2020 ; Nastase et al., 2019 ), that time courses 𝐂 ∈ℝ N ×J and spatial patterns 𝐒 ∈ℝ M ×J of nodes in the parcellation are stimulus-evoked responses that are consistent across subject. N is the number of timepoints, J is the number of patterns, M is the number of nodes. Thus, in our TCAbased model, for each pattern 𝑗, the loading 𝝈𝑖,𝑗 for subject 𝑖 is different from other subjects. In this case, for subject 𝑖 , the node time courses, 𝐗 𝑖 ∈ℝ M ×N , can be represented as: 𝐗 𝑖 = 𝐒 ×𝐝𝐢𝐚𝐠 (𝝈𝑖, ∶ )×𝐂 𝑇 + 𝐈 𝐝 𝑖 + 𝜺 𝑖 , (1) where 𝐈 𝐝 𝑖 are stimulus-evoked responses that are idiosyncratic to each subject and 𝜺 𝑖 corresponds to noise, which may reflect spontaneous neural activity and non-neural physiological and scanner signal sources. 𝐝𝐢𝐚𝐠 ( 𝝈𝑖, ∶ ) is a square matrix of order J with 𝝈𝑖, ∶ on the diagonal and other elements of the matrix equal to zero. 2.1.2. Tensor construction A tensor is constructed by stacking the multi-subject data, 𝐗 𝑖 , to create 𝐗 ∈ℝ M ×I ×N , with I equal to the total number of subjects. Although the spatial and temporal patterns are assumed to be shared across participants, the loadings on the components in each participant are different. Of note, the model also does not place any assumptions on the distribution of time courses or spatial distributions. 2.1.3. TCA unmixing model TCA is a basic model for tensor decomposition that is unconstrained. With different constraints, different algorithms can be derived, including non-negative canonical polyadic decomposition (NCPD) ( Zhou et al., 2014 ) and independent constrained CPD (e.g. tensor-ICA Beckmann and Smith, 2005 ). In the TCA model, a third-order tensor 𝐗 can be represented as the sum of several rank-1 tensors and a residual tensor 𝐄 ( Hitchcock, 1927 ), which is illustrated in Fig.1 (C). The mathematical formula is as follow: 𝐗 = 𝑅 ∑ 𝑟 =1  𝐬 𝑟 ◦ 𝝈𝑟 ◦ 𝐜 𝑟 + 𝐄 = 𝑅 ∑ 𝑟 =1 𝐗 𝑟 + 𝐄 , (2) where vectors  𝐬 𝑟 ∈ℝ M ×1 ,  𝝈𝑟 ∈ℝ I ×1 and  𝐜 𝑟 ∈ℝ N ×1 are 𝑟 𝑡ℎ estimated tensor component (TC) spatial distribution 𝐒 , subject loadings 𝝈and time courses 𝐂 respectively. The operator ◦represents the outer product of vectors. 𝐗 𝑟 represents the rank-1 tensor that constructed by the corresponding components of each mode. The idiosyncratic stimulus-evoked responses, spontaneous signals and noise signals are contained in the residual tensor 𝐄 and 𝑅 is the number of extracted patterns. Ideally, the number of tensor components, or model order, should be equivalent to the number of patterns, J , shared across subjects. However, in real-world applications, the number of spatio-temporal patterns is unknown. In the present study, we propose a novel method for determining model order according to stability of tensor spectral clustering. Once a solution is identified, for each tensor component, the spatial distributions are  𝐬 𝑟 , time courses are  𝐜 𝑟 . These are the same across subjects, with subject loadings,  𝝈𝑟 , reflecting between-subject variation in the strength of the patterns. 2.2. TCA estimation algorithm The alternating least-squares (ALS) algorithm ( Cichocki et al., 2015 ; Kolda and Bader, 2009 ) is used to estimate factor matrices 𝐒 , 𝐂 and 𝝈. The ALS algorithm proceeds by fixing two of the factor matrices to optimize over the third factor matrix. For example, while time courses, 𝐂 , are being estimated, the spatial patterns, 𝐒 , and subject loadings, 𝝈, are fixed. The time courses are updated with the following rule: 𝐂 ← argmi n  𝐂 1 2 |||||| 𝐗 − 𝑅 ∑ 𝑟 =1  𝐬 𝑟 ◦ 𝝈𝑟 ◦ 𝐜 𝑟 |||||| 2 𝐹 , (3) where 𝐹 represents the Frobenius norm. The updating rule is solved as a linear least-squares problem that is convex and has a closed-form solution. The other factor matrices are solved with the same updating rule and the three factor matrices are updated in an alternating fashion until the stop criteria is met, in this case when the absolute difference of the fits between two adjacent iterations is less than 1e-8, or when the number of iterations exceeds 1000. The ALS algorithm is provided open access from the tensor toolbox ( https://www.tensortoolbox.org ). 2.3. Model order selection via tensor spectral clustering Similar to ICA, model order (number of extracted components) selection is a significant methodological concern when applying these data driven algorithms for fMRI data analysis ( Abou-Elseoud et al., 2010 ; Beckmann, 2012 ; Kuang et al., 2018 ). Information-theoretic criteria (ITC) have been used in numerous signal processing applications to estimate model order, including minimum code length based minimum description length (MDL) criterion ( Rissanen, 1978 ), Akaike information criterion (AIC) ( Akaike, 1998 ), and Bayesian information criterion (BIC) ( Rissanen, 1978 ). However, studies have demonstrated that estimations based on different criteria may be different and these criteria may lose efficacy when the signal noise ratio is low ( Cong et al., 2011 ; Hu et al., 2020 ). In this study, model order is selected based on reproducibility of the decomposition results. We use the reproducibility of the estimated rank-1 tensor 𝐗 𝑟 =  𝐬 𝑟 ◦ 𝐜 𝑟 ◦ 𝝈𝑟 via a novel tensor spectral clustering approach to select the model order. In this case, the tensor decomposition is done for a range of model orders, the algorithm stability under each model order is evaluated via tensor spectral clustering, and the model order with the highest algorithm stability index is selected as the appropriate decomposition. In the stability analysis, for the given dataset, the same algorithm with the same parameters is ran 𝐾times, each with different initial conditions. For each model order, 𝑅 , 𝑅 ×𝐾components are estimated across runs for each mode (temporal, spatial, subject). The similarity matrices (or adjacency matrices calculated by correlating each pair of components for a mode) for each mode 𝐖 ( 𝐒 ) , 𝐖 ( 𝝈) , 𝐖 ( 𝐂 ) ∈ℝ 𝑅𝐾×𝑅𝐾 are then fed into tensor spectral clustering, which is a co-clustering method that enables fusing and assessing the stability information of different modes simultaneously. Details of the formulation for tensor spectral clustering are provided in the Appendix. In tensor spectral clustering, the number of clusters is set to equal number of extracted components 𝑅 . Stable components will form a tight cluster, with a stability index quantified as the average of the intra-cluster similarities. Ideally, if the estimated components for a given model order are stable, the intra-cluster similarity of the corresponding cluster is close to 1, whereas the stability index for unstable components will be close to 0. The algorithm stability is defined as the average of component stability indices. 2.4. Simulations We demonstrate the effectiveness of the proposed framework using numerical simulations performed in MATLAB. The simulated spatial maps and time courses for 29 ICA components (SimBT ( Erhardt et al., 2012 ) http://mialab.mrn.org/software ) and time courses 4 G. Hu, H. Li, W. Zhao et al. NeuroImage 255 (2022) 119193 Fig. 2. Ground truth spatio-temporal signal sources used in simulations. The left figure shows the spatial distribution of 29 ICA components. The right column shows time courses for consistent components across subjects (the first row, 𝑪 ), idiosyncratic components (the second row, 𝑰 𝒅 ) and for spontaneous components (the third row, 𝜺 ). ( http://mlsp.umbc.edu/simulated_fmri_data.html ) are shown in Fig. 2 . Four components were simulated to represent consistent spatiotemporal signals across subjects (with time courses shown in the first row of Fig. 2 ). There are also three idiosyncratic components for each subject, similar to scanner or motion artifacts present in fMRI data (second row of Fig. 2 ). These components were given a random circular shift for each subject. Time courses in the third row of Fig. 2 represent spontaneous brain activity included in the simulations. Different spontaneous components were generated independently for each subject following a Gaussian distribution. The total signal noise ratio (SNR) was fixed at 2dB. We generated a random node weight matrix 𝐒 ∈ℝ 29 ×8 and subject loading matrix 𝝈∈ℝ 10 ×8 to generate 10 participant datasets using equation (1) . All subject’s data were then stacked to construct a three-order tensor 𝐗 Simulation ∈ℝ 29 ×10 × 100 . This tensor was then decomposed with model orders ranging from 2 to 10. Under each model order, the TCA algorithm was run 20 times and the stability index was calculated with tensor spectral clustering. The correlation coefficients between the estimate components for the model order with the highest stability index and the corresponding ground truth were used as a criterion to evaluate the performance of the proposed framework. 2.5. Motor task fMRI experiment The effectiveness of the proposed framework is also demonstrated with conventional motor task fMRI in which we know the ground truth stimulation time courses that can be used to model the brain task activity. Motor task fMRI data from 100 healthy unrelated subjects (2236 years) were utilized from the WU-Minn Human Connectome Project (HCP; Van Essen et al., 2013 ). During the motor task, participants were presented with visual cues to either tap their left or right fingers, squeeze their left or right toes, or move their tongue. Blocks for each movement type lasted 12 seconds and were preceded by a 3 second cue. For each subject, a total of 284 timepoints were collected with TR = 0.72s. Additional details of the motor task for the HCP are provided in Barch et al. (2013) . 2.5.1. Preprocessing The fMRI data went through the HCP’s standard minimal preprocessing pipeline (motion correction, distortion correction, high-pass filtering (200s), and nonlinear alignment to MNI template space ( Glasser et al., 2013 )). 2.5.2. GLM In order to evaluate the performance of TCA, we used the standard GLM as the gold standard, which was applied to the same data. A conventional mixed-effects GLM analysis was conducted to derive activation estimates for each reference function, as described in Barch et al. (2013) . Five predictors were included in the Motor model —right hand, left hand, right foot, left foot, and tongue. Time courses of visual cues for each type of movement were convolved with Hemodynamic Response Function (HRF) and low-pass filtered with a highfrequency cutoff of 0.5 Hz. Spatial maps of linear contrasts (Cues, Foot, Tongue, Hand) of the parameter estimates were computed to compare each condition to baseline. Higher level GLMs, including age and sex as covariates, were used to estimate the spatial distributions of tongue, foot and hand movement and total movement evoked brain activity via one sample t-tests of the contrast of parameter estimate maps, with a correction for multiple tests to control the family-wise error (FWE). Since the evoked response to the motor mapping task was so strong ( Barch et al., 2013 ), a higher threshold (p < 0.001) was used. 2.5.3. TCA Following the recommendations of Pervais et al. (2020) , a higher model order functional parcellation based on spatial ICA components was used as the parcellation scheme for our study. Spatial ICA components provided by the HCP with a model order of 300 were used as templates with dual regression ( Nickerson et al., 2017 ) to derive the subject-specific time courses of each node, as shown in Fig. 1 (A). Postprocessing of time courses included: (1) removal of the first 10 timepoints (2) demeaning, (3) detrending linear, quadratic and cubic trends, and (4) variance normalization. The time courses for each node for each subject were then stacked in the subject dimension to create a tensor, 𝐗 Motor ∈ℝ 300 × 100 × 284 . The tensor was then decomposed with TCA, with model orders ranging from 2 to 20. The decomposition was repeated 20 times for each model order to calculate the reproducibility of the estimated components and the model order via our proposed tensor spectral clustering strategy. Task-related TCA components were identified via similarity of component time courses to the “ground truth ”time courses from the GLM analyses. 2.5.4. ISFC We compared our approach with another popular approach for naturalistic fMRI data analysis, ISFC ( Simony et al., 2016 ). ISFC was calculated with a sliding window of 20s (28 TRs). At each time point,t, the pairwise correlations between all nodes were calculated over the window interval (t, t + 20) to construct a correlation matrix for time point t. Correlation matrices were computer for each window shifted by 0.72s (1 TR) in time. The resulting sliding window correlation matrices were then analyzed with k-means clustering, with the number of clusters set 5 G. Hu, H. Li, W. Zhao et al. NeuroImage 255 (2022) 119193 from 2 to 20. For each k-means run, cluster indices of each timepoint were used to represent the occurrence of each task-related brain state (cluster). The spatial distributions of each brain state were derived by calculating the average connectivity of each node with all other nodes in the cluster. The time course of each state is calculated as the average occurrence across subjects. For ISFC, “ground truth ”time courses were assumed to be the GLM stimulation time courses convolved with the HRF and the sliding window. 2.6. Naturalistic stimuli fMRI experiment We used a naturalistic fMRI dataset that was also collected by the HCP ( Van Essen et al., 2013 ). Healthy adult participants (aged 22-35 years) underwent fMRI using a 7 T Siemens Magnetom scanner (voxel size = 1.6mm 3 , TR = 1s) during a movie-watching paradigm. The sample used here (n = 184) reflects all available data for this paradigm. Each subject underwent four runs of 15-min movie watching (MOVIE1MOVIE4). Each run comprised five video clips presented in a fixed order. FMRI data were also collected during a validation clip (the fifth video clip), consisting of a brief montage (1min 22s) of moving scenes depicting people and landscapes. The music played along with the movie scenes contained features to evoke consistent brain activity across participants ( Alluri et al., 2012 ). Five music features (Fluctuation Centroid, Fluctuation Entropy, Key Clarity, Mode and Pulse Clarity) were extracted using the MIRToolbox ( Lartillot, O., 2007 ). These music features facilitate the interpretation of the estimated tensor components. In addition, this movie included a social scene. To demonstrate relationships between inter-subject variability in TCA components with behavior and facilitate interpretation of results, the subject loadings on social scene-related TCA components were assessed for correlations with behavioral measures from the Semi-Structured Assessment for the Genetics of Alcoholism (SSAGA; BUCHOLZ et al., 1994 ; Hesselbrock et al., 1999 ) and NIH Toolbox. We selected seven measures related to social function, including Friendship, Loneliness, Perceived Hostility, Perceived Rejection, Emotional Support, Instrumental Support and Antisocial Personality Problems Raw Score. All fMRI analyses utilized the FIX-denoised data that underwent standard HCP minimal preprocessing (motion correction, distortion correction, high pass filtering, and nonlinear alignment to MNI template space ( Glasser et al., 2013 )) plus regression of 24 framewise motion estimates (six rigid-body motion parameters and their derivatives and the squares of those 12) and regression of confound components identified via ICA ( Griffanti et al., 2014 ; Salimi-Khorshidi et al., 2014 ). Details of data acquisition and basic data preprocessing can be found in previous studies ( Glasser et al., 2013 ; T. Vu et al., 2016 ; Van Essen et al., 2012 ). Both TCA and ISC were applied to same FIX-denoised dataset. TCA was implemented exactly the same way as for the analysis of the motor task fMRI data (shown in Fig. 1 ), with the resulting tensor, 𝐗 Nat uralist ic ∈ ℝ 300 × 184 ×82 , decomposed with our TCA approach. ISC was also applied to the HCP naturalistic stimuli fMRI to compare the performance of TCA and ISC. In ISC, the same parcellation strategy used for the TCA framework was also applied. In ISC framework, IS-RSA was used to evaluate the relationship between shared brain networks across subjects and subject traits. For each subject pair, brain similarity was calculated as the correlation of activity time courses. Behavioral similarity was calculated according to absolute value of the difference in behavioral data of subject pair. Representational similarity was assessed by calculating correlation between the brain and behavioral similarity matrices. Details of IS-RSA can be found in the previous study ( Finn et al., 2020 ; Kriegeskorte et al., 2008 ; Mantel, 1967 ; Meer et al., 2020 ). In order to demonstrate that our proposed framework is also suitable for naturalistic stimuli with longer durations, we used the naturalistic stimuli fMRI data provided by Meer et al. (2020) , which implemented a paradigm with 20 min movie stimuli (The Butterfly Circus). FMRI data were collected with identical movie stimuli on two separate occasions (Movie view A and Movie view B) with an interval of three months between sessions. 17 participants completed both movie viewing sessions. Two participants were excluded because of in-scanner head motion. The fMRI data was collected with TR of 2200msec and 535 volumes were acquired. Standard preprocessing was performed using fMRIPrep ( Esteban et al., 2019 ) and included slice-timing, motion correction, co-registration to the structural image, spatial normalization to MNI space, and spatial smooth with a 6 mm Gaussian kernel. ICAAROMA was subsequently performed using non-aggressive denoising. Details of the data collection, pre-processing, and brain state dynamics can be found in Meer et al. (2020) . The Butterfly Circus narrates an intense, emotionally evocative story of a man born without limbs who is encouraged by the showman of a renowned circus to overcome obstacles of self-worth and reach his own potential. Annotations for (i) the use of language, (ii) change of scenes, (iii–v) Positive/Negative Faces, and (vi–viii) Positive/Negative Scenes during the movie are also provided. Meer et al. (2020) applied a hidden Markov model (HMM) to identify brain states dynamics for 14 canonical brain networks (BN) derived with group ICA ( Shirer et al., 2012 ). In order to compare the results with the findings of Meer et al. (2020) , we used the same set of BNs for parcellating the fMRI data. The first 5 time points for each BN average temporal course were discarded. Thus, for each subject we obtained a matrix ( 530 ×14 ) for each movie viewing session. The matrices for movie viewing A and movie viewing B were then stacked and concatenated in the subject dimension to construct the tensor, 𝐗 ∈ℝ 30 × 530 ×14 , that was decomposed via our TCA approach. 3. Results 3.1. Simulations Fig. 3 shows the simulation results. Fig. 3 (A) shows the algorithm stability indices for the different model orders. We found that the algorithm stability curve reached its peak at model orders of 3 and 4. To extract more information about the simulated brain responses, results at model order 4 were further analyzed. Note that the algorithm stability index at this model order is 1, which means all components estimated at that model order are the same for all 20 runs, indicating that all the components are reproducible. The spatial distributions, time courses and subject loadings for each component are shown in Fig. 3 (B-E). Based on the correlation coefficient between estimated component time courses and spatial maps with the ground truth time courses and spatial maps, all four components are successfully estimated. In order to highlight the activated and deactivated brain regions, the spatial distribution of the estimated component is shown with a proportional threshold of 5% ( Garrison et al., 2015 ), with warm colors representing activation (relative to the within-component global average) and cool colors representing deactivation (relative to the within-component global average) within each component. 3.2. Motor task fMRI experiment For motor task fMRI data, for TCA, the algorithm stability curve reached a peak at model order of 4. Hence, the estimated TCA components with model order of 4 were selected for further analysis. For ISFC, the number of clusters is selected as 6. With this number of clusters, the estimated time courses have the greatest similarity to the corresponding ground truth time courses. The spatial distributions and time courses of the estimated components, and the results from applying the GLM and ISFC for analysis are shown in Fig. 4 , with display and thresholding at |Z| > 2.3 after transformation into Z-scores across the spatial domain, and warm colors representing activation and cool colors representing deactivation. For each component, after normalization, both stimulus timing and estimated time courses are demonstrated in the same subfigure to identify of the estimated components. In the experiment, the embedded components that correspond to the tongue, foot and cues are successfully 6 G. Hu, H. Li, W. Zhao et al. NeuroImage 255 (2022) 119193 Fig. 3. Simulation results. (A) Algorithm stability index for different model orders. The algorithm reaches the greatest stability at model order of 4, which is exactly the same as the number of consistent components. (B) - (E) estimated components and ground truths signals. 𝑟 is the correlation coefficient between the estimate component spatial map/time course and the corresponding ground truths. For spatial distributions, a proportional threshold of 5% ( Garrison et al., 2015 ) was used, with warm colors representing activation (relative to the within-component global average) and cool colors representing deactivation (relative to within-component global average) within each component. For time courses, the red represents ground truth and the black indicates estimations. estimated. However, the component corresponding to hand movements failed to be detected by TCA, which is due to poor temporal consistency across subjects (see Supplementary Fig.1). E.g., components with worse consistency across subjects will reduce algorithm stability and, in turn, the model order that may reveal such components will not be selected as the final model order. Inspection of the spatial maps and time courses in Fig. 4 shows that the estimation of spatial distributions of activation are very consistent across methods. Comparison of the estimated time courses of the visual cues from TCA and ISFC shows that the TCA time courses resemble the visual cues after being convolved with hemodynamic response function. The time courses of ISFC are similar to the stimulation time course convolved with the sliding window, which results in a distorted representation for the visual cues ( Fig. 4 (A)). Such distortions may hinder the investigation of the relationship between brain activity and movie features in naturalistic stimuli. However, ISFC successfully estimated the spatial map and time course of the hand movements whereas TCA did not at this model order (see Supplementary Fig. 2). 3.3. Naturalistic stimuli fMRI experiment For the short naturalistic stimuli fMRI from the HCP, the TCA algorithm reaches the highest stability at model order equal to 3 and all three tensor components estimated under this model order are highly reproduced (see Supplementary Fig. 3). The estimated components under this model order were selected for further analysis. The spatial distribution of the first tensor component estimated with TCA and the corresponding map estimated by ISC analyses of these data are shown in Fig. 5 (A). Bilateral occipital fusiform gyrus, lingual gyrus and superior temporal gyrus are identified as activated regions in both the TCA component and the ISC spatial map. However, bilateral postcentral gyrus, superior parietal lobule and left precentral gyrus are also identified as activated in the ISC map, but not in the TCA component. Bilateral cerebellum and lateral occipital cortex are activated in the TCA component. Fig. 5 (B) shows a bar plot of the average correlation coefficients between the subject time courses extracted from the lingual gyrus and temporal gyrus, and the temporal course of the TCA component. This shows that the time course of lingual gyrus is negatively correlated with the estimated TCA component time course, whereas the time course of the temporal gyrus is positively correlated with the estimated time course. This opposite response to the movie stimuli of the lingual gyrus and temporal gyrus is not identified by ISC. Fig. 5 (C) shows the estimated time course of the first tensor component with the onset of the annotated stimuli in the movie to illustrate the relationship between them. Inspection of the timing suggests that landscape scenes result in increases in BOLD signal for this mode, while scenes with social content are associated with decreases in BOLD signal. In an exploratory analysis, we applied the proposed framework on a social cognition task fMRI from the HCP dataset. Then we compared the spatial distribution between social task evoked brain network (Supplementary, Fig. 8) and the movie stimuli evoked component, and found that the social task network was significantly correlated (r = 0.27, P = 10 − 6 ) with the TCA component. The network estimated with ISC was also significantly correlated with the social task evoked brain network (r = 0.28, P = 10 − 6 ). Correlations between the subject loadings for the TCA component and the behavioral measures related to social function showed a significant correlation (FDR corrected, P < 0.05) between TCA component loadings and antisocial personality score ( Fig. 5 (D)). A significant correlation (r = 0.02, P = 0.01) between antisocial personality scores and shared brain activities across subjects estimated with ISC is also identified with IS-RSA but the correlation coefficient is only 0.02. The time course of the second tensor component was significantly correlated with the music feature Pulse Clarity (r = 0.39, P = 10 − 4 ), as shown in Fig. 5 (E). But the time courses estimated with ISC failed to correlated with the music feature (r = 0.19, P = 0.10). All three estimated tensor component spatial maps show effects in bilateral cerebellum, occipital fusiform gyrus, and lateral occipital cortex. The lingual gyrus is deactivated in the first TCA component, but not in the second or third components. Left frontal pole, which does not appear in the first TCA component, is deactivated in both the second and third TCA components, and the bilateral superior temporal gyrus is activated in the first and second TCA components, but not the third. Thus, TCA is able to identify overlapping brain activation modes. Fig. 6 demonstrates that our proposed framework is also suitable for fMRI data with long-duration movie stimuli. Three components were estimated with TCA. The spatial distributions show a strong correlation with the spatial distributions reported in Meer et al. (2020) using HMM (correlation indicated above each map). The first TCA component corresponds to the language state estimated with the HMM (r = 0.86, P = 10 − 5 ). The time course of the component was significantly correlated with annotation of the use of language (convolved with HRF, r = 0.12, P = 0.006), as shown in Fig. 6 (A). The second component corresponds to an interoception state (r = 0.72, P = 0.004) and the third component corresponds to task-positive network (r = 0.88, P = 10 − 5 ). 7 G. Hu, H. Li, W. Zhao et al. NeuroImage 255 (2022) 119193 Fig. 4. Comparison of TCA, GLM and ISFC. For spatial distributions, warm colors represent activation (relative to within-state global average) while cool colors represent deactivation (relative to within-state global average) within each state. In the plots of the time courses, the red lines represent estimated time courses with TCA, blue lines are the estimated time courses with ISFC, and the black dash lines represent the ground truth time courses. Note that for TCA, the ground truth time courses are the experiment stimulus time course convolved with the HRF, which is the theoretical neural response of the task. For ISFC, the ground truth time courses are the experiment stimuli time courses convolved with the sliding window (28TRs). The occupancy of these three components was shown to be higher during movie watching than during resting-state ( Meer et al., 2020 ), which demonstrates that TCA can estimate brain activity evoked with naturalistic stimuli distinct from spontaneous brain activity. We also found that subject loadings during movie viewing B are lower than for movie viewing A, so even though similar spatial and temporal patterns of activity are evoked with repeated stimuli, differences in brain activation corresponding to repeated viewing of the same stimuli can also be assessed with TCA. Last, we evaluated the relationship between subject loadings and scores from the post-movie questionnaires that were administered to participants and found that subject loadings on the first tensor component for movie viewing B were significantly correlated with boredom (P = 0.02, uncorrected), which may explain the smaller subject loadings reflecting reduced strength on this pattern during the repeat viewing. 4. Discussion In this study, an analysis framework for applying TCA to fMRI data to discover spatio-temporal components that are shared across subjects during naturalistic movie viewing was proposed. In this framework, a third-order tensor is constructed from the timeseries extracted from all brain regions from a given parcellation, for all participants. This tensor is then decomposed via TCA to identify modes corresponding to the spatial distribution (maps) and time series that are common to all participants, and subject loadings that reflect between-subject variability in the patterns of brain activation in response to naturalistic stimuli. The stability of the extracted components is evaluated with a novel clustering method, tensor spectral clustering, to guarantee the reproducibility of the results. Model order is selected based on the stability of the results across a range of model orders. Extensive testing ( Figs. 3–5 and 8