Elsevier

Econometrics and Statistics

Available online 2 March 2021
Econometrics and Statistics

Spatial-Temporal Analysis of Multi-Subject Functional Magnetic Resonance Imaging Data

https://doi.org/10.1016/j.ecosta.2021.02.006Get rights and content

Abstract

Functional magnetic resonance imaging (fMRI) is one of the most popular neuroimaging technologies used in human brain studies. However, fMRI data analysis faces several challenges, including intensive computation due to the massive data size and large estimation errors due to a low signal-to-noise ratio of the data. A new statistical model and a computational algorithm are proposed to address these challenges. Specifically, a new multi-subject general linear model is built for stimulus-evoked fMRI data. The new model assumes that brain responses to stimuli at different brain regions of various subjects fall into a low-rank structure and can be represented by a few principal functions. Therefore, the new model enables combining data information across subjects and regions to evaluate subject-specific and region-specific brain activity. Two optimization functions and a new fast-to-compute algorithm are developed to analyze multi-subject stimulus-evoked fMRI data and address two research questions of a broad interest in psychology: evaluating every subject’s brain responses to different stimuli and identifying brain regions responsive to the stimuli. Both simulation and real data analysis are conducted to show that the new method can outperform existing methods by providing more efficient estimates of brain activity.

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 cm×2 cm×2 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 yvi(t) be the fMRI measurement of brain activity at voxel v of subject i at time t, t=κ,2κ,,Tκ, i=1,2,,n, and v=1,2,,V, where T is the total number of time points of fMRI time series, V 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 K different stimuli are presented in the fMRI experiment. Let ski(t) be the stimulus function that indicates the evoked times of the kth 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 W={Wk,k=1,,K}, d˜={di,i=1,,n

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 n=106 subjects’ fMRI time series at V=153 voxels, which are distributed over a 15×15×15 lattice grid. The data have K=4 stimuli. We let voxels in the center 7×7×7 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 n=106 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 PSSE(Θ) and PSSEc(Θ) produce similar HRF estimates, possibly because the differences between subject-specific coefficients P1 and P2 are small in comparison to the data error. With an additional sparsity penalty, PSSEc(Θ) 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)

  • O. Friman et al.

    Detection and detrending in fmri data analysis

    NeuroImage

    (2004)
  • K. Friston et al.

    Event-related fMRI: characterizing differential responses

    NeuroImage

    (1998)
  • K. Friston et al.

    To smooth or not to smooth?

    NeuroImage

    (2000)
  • R. Ge et al.

    Improved fastica algorithm in fmri data analysis using the sparsity property of the sources

    Journal of neuroscience methods

    (2016)
  • G. Glover

    Deconvolution of impulse response in event-related BOLD fMRI

    NeuroImage

    (1999)
  • M. Jenkinson et al.

    Improved optimization for the robust and accurate linear registration and motion correction of brain images

    Neuroimage

    (2002)
  • B. Knutson et al.

    Fmri visualization of brain activity during a monetary incentive delay task

    Neuroimage

    (2000)
  • N. Lange et al.

    Plurality and resemblance in fMRI data analysis

    NeuroImage

    (1999)
  • C. Liao et al.

    Estimating the delay of the fmri response

    NeuroImage

    (2002)
  • S. Makni et al.

    Bayesian deconvolution of fmri data using bilinear dynamical systems

    NeuroImage

    (2008)
  • B.A. Olshausen et al.

    Sparse coding of sensory inputs

    Current opinion in neurobiology

    (2004)
  • F. Pedregosa et al.

    Data-driven hrf estimation for encoding and decoding models

    NeuroImage

    (2015)
  • B.A. Vogt

    Midcingulate cortex: structure, connections, homologies, functions and diseases

    Journal of chemical neuroanatomy

    (2016)
  • K. Worsley et al.

    Analysis of fMRI time-series revisited again

    NeuroImage

    (1995)
  • K. Worsley et al.

    A general statistical analysis for fMRI data

    NeuroImage

    (2002)
  • L. Zhang et al.

    A spatio-temporal nonparametric bayesian variable selection model of fmri data for clustering correlated time courses

    NeuroImage

    (2014)
  • L. Zhang et al.

    A spatio-temporal nonparametric bayesian variable selection model of fmri data for clustering correlated time courses

    NeuroImage

    (2014)
  • T. Zhang et al.

    A semi-parametric model of the hemodynamic response for multi-subject fmri data

    NeuroImage

    (2013)
  • T. Zhang et al.

    A low-rank multivariate general linear model for multi-subject fmri data and a non-convex optimization algorithm for brain response comparison

    NeuroImage

    (2018)
  • D. Zhu et al.

    Fusing dti and fmri data: a survey of methods and applications

    NeuroImage

    (2014)
  • A. Agarwal et al.

    Noisy matrix decomposition via convex relaxation: Optimal rates in high dimensions

    The Annals of Statistics

    (2012)
  • J. Bai et al.

    Forecasting using princpal components from a large number of predcitors

    Journal of the American Statistical Association

    (2002)
  • J. Bai et al.

    Statistical analysis of factor models of high dimension

    The Annals of Statistics

    (2012)
  • L. Baldassarre et al.

    Structured sparsity models for brain decoding from fmri data

    2012 Second International Workshop on Pattern Recognition in NeuroImaging

    (2012)
  • I.N. Beloozerova et al.

    Activity of different classes of neurons of the motor cortex during locomotion

    Journal of Neuroscience

    (2003)
  • de Boor, C., 2001. A practical guide to splines (applied mathematical sciences vol....
  • G. Bush et al.

    Dorsal anterior cingulate cortex: a role in reward-based decision making

    Proceedings of the National Academy of Sciences

    (2002)
  • R.B. Buxton et al.

    A model for the coupling between cerebral blood flow and oxygen metabolism during neural stimulation

    Journal of Cerebral Blood Flow & Metabolism

    (1997)
  • R. Casanova et al.

    Evaluating the impact of spatio-temporal smoothness constraints on the BOLD hemodynamic response function estimation: an analysis based on Tikhonov regularization

    Physiological Measurement

    (2009)
  • L. Chaari et al.

    Adaptive hemodynamic-informed parcellation of fmri data in a variational joint detection estimation framework

    in 15th Proceedings MICCAI, LNCS

    (2012)
  • 1

    Equally contributing first authors.

    2

    Dr. Zhang’s research is supported by NSF-2048991.

    View full text