Spatial-Temporal Analysis of Multi-Subject Functional Magnetic Resonance Imaging Data
Introduction
Functional magnetic resonance imaging (fMRI) measures human brain activity non-invasively through blood-oxygen-level-dependent contrast (Ogawa et al., 1992) and records activity of the entire brain with a high spatial resolution. Therefore, fMRI is one of the most popular neuroimaging technologies used in human brain research. In many psychological and medical studies, researchers use different stimuli to evoke subjects’ brain activities during fMRI recordings. The purposes of this paper are to analyze these multi-subject, stimulus-evoked fMRI data, evaluate each subject’s brain responses to different stimuli, and identify voxels that have different responses to different stimuli.
The analysis of stimulus-evoked fMRI data usually uses a general linear model (GLM) (Friston, Holmes, Worsley, Poline, Frith, Frackowiak, 1995, Goutte, Nielsen, Hansen, 2000, Worsley, Friston, 1995), in which the fMRI time series follows a linear regression of the stimulus sequence through hemodynamic response functions (HRF). Under the GLM, the HRF describes each subject’s brain response to each stimulus type at each region and, therefore, is the estimation target of fMRI data analysis.
Several challenges are present when analyzing multi-subject, stimulus-evoked fMRI data. First, the data are massive in size. Typical fMRI data of one subject from one experimental session consist of 100,000 to 200,000 spatially indexed time series. Each time series contains measurements of brain activity at one brain voxel and at several hundred time points—with a time unit ranging from 1 second(s) to 2 s. A voxel is a small 3D cubic volume in the brain (2 cm2 cm2 cm for the data under study). Consequently, fMRI data analysis can be computationally intensive. Second, the HRF varies across regions, subjects, and stimulus types presented to the subject. Third, fMRI data have complex spatial and temporal properties. Fourth, each subject’s fMRI data usually have a weak signal-to-noise ratio (SNR) (Lindquist, 2008), which brings an additional difficulty to accurate HRF estimation.
Most existing approaches to estimating HRFs attend to one or part of the challenges described above. For example, to address the first two challenges, researchers commonly use voxel-wise analysis (Friston, Holmes, Worsley, Poline, Frith, Frackowiak, 1995, Goutte, Nielsen, Hansen, 2000, Worsley, Friston, 1995), that is, analyze fMRI time series at each voxel independently. Although this approach is computationally fast and accommodates variations of brain activity across voxels, the approach ignores spatial information of fMRI data and association among voxels. Despite the development of joint models for brain activities at different locations Zhang, Guindani, Versace, Vannucci, 2014, Zhang, Guindani, Versace, Engelmann, Vannucci, et al., 2016; Bowman et al. (2008); Degras and Lindquist (2014a); Makni, Ciuciu, Idier, Poline, 2005, Makni, Beckmann, Smith, Woolrich, 2008; Pedregosa et al. (2015), these approaches tend to use many parameters to accommodate the variation of brain responses across many voxels and subjects and require strong regularization or strong priors to reduce estimation errors. Moreover, existing methods are mainly focused on examining brain responses to each stimulus, which limits the ability to compare brain responses to different stimuli, as explained below.
Many psychological fMRI studies use the brain response to one stimulus as the reference and compare the brain responses to other stimuli with respect to the reference. This approach is to correct for the nuisance brain response caused by the fMRI data collection process, for example, the brain response due to anxiety and uncomfortableness experienced during fMRI scanning (Mazard, Mazoyer, Etard, Tzourio-Mazoyer, Kosslyn, Mellet, 2002, Szameitat, Shen, Sterr, 2009). In this situation, the brain response to the stimulus of interest beyond the response to the reference stimulus should represent more accurately the brain activity due to the stimulus of interest in reality. Existing methods for comparing brain responses are still through a two-step approach: estimating responses to different stimuli first and then performing tests on the response estimates through voxel-wise analysis. However, this two-step approach leads to many false discoveries (Zhang et al., 2018).
We develop a new joint model for multi-subject, stimulus-evoked fMRI data to address limitations in existing fMRI data analysis methods. With an assumption of a low-rank structure for HRFs of different voxels and subjects, the new model, called multi-subject, low-rank general linear model (MSLRGLM), uses much fewer parameters than nonparametric models to accommodate the variation of brain responses across different subjects and voxels. In addition to model flexibility, the new model also enjoys a better estimation efficiency of HRFs than existing methods.
We formulate the MSLRGLM estimation as an optimization problem. Specifically, we propose two optimization functions to address two research questions: estimating subject-specific and voxel-specific HRFs and identifying voxels with different responses to different stimuli. We use regularization penalties in the optimization functions to incorporate spatial and temporal information of fMRI data. We include a sparsity-inducing penalty in the optimization function to identify voxels that have different responses to different stimuli. We use the sparsity to produce stable and interpretable brain response estimates, similar to many existing approaches (Baldassarre, Mourao-Miranda, Pontil, 2012, Ge, Wang, Zhang, Yao, Zhang, Long, 2016), because sparse neural responses have been reported in the literature (Olshausen, Field, 1996, Olshausen, Field, 2004, Beloozerova, Sirota, Swadlow, 2003, Wixted, Squire, Jang, Papesh, Goldinger, Kuhn, Smith, Treiman, Steinmetz, 2014). To address the computational challenge in estimating HRFs of many subjects and voxels, we develop a new computational algorithm to minimize the optimization functions. We show that the proposed new model and computational algorithm provide more efficient HRF estimates with smaller variances than existing methods through both simulation and real data analysis. In summary, our new modeling and estimating approaches address the challenges in fMRI data analysis.
The rest of the article is organized as follows. Section 2 reviews existing HRF estimation methods under the GLM framework and introduces a new multi-subject GLM for fMRI data of different subjects and voxels. We present our estimation approach with new optimization functions in Section 3 and propose a new computational algorithm to minimize the functions in Section 4. The proposed HRF estimation method is compared with several existing ones through simulation studies in Section 5. Section 6 presents analysis results of multi-subject fMRI data collected in a psychological study on emotion. We evaluate each subject’s brain responses to emotional stimuli and identify voxels responsive to the stimuli using the proposed method. Section 7 concludes.
Section snippets
The General Linear Model
Let be the fMRI measurement of brain activity at voxel of subject at time and where is the total number of time points of fMRI time series, is the total number of voxels under study, and is the time interval between consecutive 3D images. For the fMRI data under analysis, equals 2 s.
Suppose different stimuli are presented in the fMRI experiment. Let be the stimulus function that indicates the evoked times of the th stimulus in
Model Estimation
We use the model (3) to estimate HRFs of different subjects, voxels, and stimulus types and identify brain voxels with different responses to different stimuli. We propose two optimization functions with different penalties to achieve these two goals.
Optimization Algorithms
A crucial challenge in analyzing fMRI data is computationally efficiently estimating massive HRFs of many subjects and voxels. For the presented real data analysis under study (Section 6), despite the low-rank structure for HRFs, the total number of model parameters needed for representing all the HRFs is more than 5.5 million. Consequently, we develop a computationally efficient algorithm to solve the above two optimization problems and estimating massive HRFs.
Let
Simulation Study
We generate fMRI data from the GLM (1) using the same experimental design as that of the real data. The simulated data contain subjects’ fMRI time series at voxels, which are distributed over a lattice grid. The data have stimuli. We let voxels in the center grid have different HRFs associated with the first two stimuli.
We simulate HRFs from a semi-parametric model developed by Zhang et al. (2013). This model characterizes subject-specific and voxel-specific brain
Real Data Analysis
The fMRI data under analysis were collected from a psychology study on the human brain’s emotion function (Coan, 2010, Coan, 2011, Coan, Schaefer, Davidson, 2006, Coan, Beckes, Allen, 2013). The fMRI experiment used threats of mild electric shocks as negative emotional stimuli to evoke subjects’ brain activities. The experiment protocol consisted of four stimuli: threat cues with 20% chances of a mild electric shock, safety cues with no electric shock, resting periods between cues, and
Discussions
The two optimization functions and produce similar HRF estimates, possibly because the differences between subject-specific coefficients and are small in comparison to the data error. With an additional sparsity penalty, enables identifying voxels with different responses to different stimuli. This voxel selection has slightly better accuracy than using a population-mean GLM for fMRI data (Zhang et al., 2018). We suppose that the larger variation of HRFs across
References (80)
Twenty years of functional mri: the science and the stories
Neuroimage
(2012)- et al.
A bayesian hierarchical framework for spatial modeling of fmri data
NeuroImage
(2008) - et al.
Cognitive and emotional influences in anterior cingulate cortex
Trends in cognitive sciences
(2000) - et al.
Fmri analysis with the general linear model: removal of latency-induced amplitude bias by incorporation of hemodynamic derivative terms
NeuroImage
(2004) - et al.
The impact of temporal regularization on estimates of the BOLD hemodynamic response function: a comparative analysis
NeuroImage
(2008) - et al.
Human ocular dominance columns as revealed by high-field functional magnetic resonance imaging
Neuron
(2001) - et al.
Childhood maternal support and neighborhood quality moderate the social regulation of neural threat responding in adulthood
International Journal of Psychophysiology
(2013) - et al.
A hierarchical model for simultaneous detection and estimation in multi-subject fmri studies
NeuroImage
(2014) - et al.
A hierarchical model for simultaneous detection and estimation in multi-subject fmri studies
NeuroImage
(2014) - et al.
Functional responses and structural connections of cortical areas for processing faces and voices in the superior temporal sulcus
Neuroimage
(2013)