Spatiotemporal wavelet analysis for functional MRI

Characterizing the spatiotemporal behavior of the BOLD signal in functional Magnetic Resonance Imaging (fMRI) is a central issue in understanding brain function. While the nature of functional activation clusters is fundamentally heterogeneous, many current analysis approaches use spatially invariant models that can degrade anatomic boundaries and distort the underlying spatiotemporal signal. Furthermore, few analysis approaches use true spatiotemporal continuity in their statistical formulations. To address these issues, we present a novel spatiotemporal wavelet procedure that uses a stimulus-convolved hemodynamic signal plus correlated noise model. The wavelet fits, computed by spatially constrained maximum-likelihood estimation, provide efficient multiscale representations of heterogeneous brain structures and give well-identified, parsimonious spatial activation estimates that are modulated by the temporal fMRI dynamics. In a study of both simulated data and actual fMRI memory task experiments, our new method gave lower mean-squared error and seemed to result in more localized fMRI activation maps compared to models using standard wavelet or smoothing techniques. Our spatiotemporal wavelet framework suggests a useful tool for the analysis of fMRI studies.     

Bayesian spatiotemporal model of fMRI data

This research describes a new Bayesian spatiotemporal model to analyse block-design BOLD fMRI studies. In the temporal dimension, we parameterise the hemodynamic response function's (HRF) shape with a potential increase of signal and a subsequent exponential decay. In the spatial dimension, we use Gaussian Markov random fields (GMRF) priors on activation characteristics parameters (location and magnitude) that embody our prior knowledge that evoked responses are spatially contiguous and locally homogeneous. The result is a spatiotemporal model with a small number of parameters, all of them interpretable. Simulations from the model are performed in order to ascertain the performance of the sampling scheme and the ability of the posterior to estimate model parameters, as well as to check the model sensitivity to signal to noise ratio. Results are shown on synthetic data and on real data from a block-design fMRI experiment.     

Spatiotemporal nonlinearity in resting-state fMRI of the human brain

In this work, the spatiotemporal nonlinearity in resting-state fMRI data of the human brain was detected by nonlinear dynamics methods. Nine human subjects during resting state were imaged using single-shot gradient echo planar imaging on a 1.5T scanner. Eigenvalue spectra for the covariance matrix, correlation dimensions and Spatiotemporal Lyapunov Exponents were calculated to detect the spatiotemporal nonlinearity in resting-state fMRI data. By simulating, adjusting, and comparing the eigenvalue spectra of pure correlated noise with the corresponding real fMRI data, the intrinsic dimensionality was estimated. The intrinsic dimensionality was used to extract the first few principal components from the real fMRI data using Principal Component Analysis, which will preserve the correct phase dynamics, while reducing both computational load and noise level of the data. Then the phase-space was reconstructed using the time-delay embedding method for their principal components and the correlation dimension was estimated by the Grassberger-Procaccia algorithm of multiple variable series. The Spatiotemporal Lyapunov Exponents were calculated by using the method based on coupled map lattices. Through nonlinearity testing, there are significant differences of correlation dimensions and Spatiotemporal Lyapunov Exponents between fMRI data and their surrogate data. The fractal dimension and the positive Spatiotemporal Lyapunov Exponents characterize the spatiotemporal nonlinear dynamics property of resting-state fMRI data. Therefore, the results suggest that fluctuations presented in resting state may be an inherent model of basal neural activation of human brain, cannot be fully attributed to noise.     

Multivariate analysis of neuronal interactions in the generalized partial least squares framework: simulations and empirical studies

Identification of spatiotemporal interactions within/between neuron populations is critical for detection and characterization of large-scale neuronal interactions underlying perception, cognition, and behavior. Univariate analysis has been employed successfully in many neuroimaging studies. However, univariate analysis does not explicitly test for interactions between distributed areas of activity and is not sensitive to distributed responses across the brain. Multivariate analysis can explicitly test for multiple statistical models, including the designed paradigm, and allows for spatial and temporal model detection. Here, we investigate multivariate analysis approaches that take into consideration the 4D (time and space) covariance structure of the data. Principal component analysis (PCA) and independent component analysis (ICA) are two popular multivariate approaches with distinct mathematical constraints. Common difficulties in using these two different decompositions include the following: classification of the revealed components (task-related signal versus noise), overall signal-to-noise sensitivity, and the relatively low computational efficiency (multivariate analysis requires the entire raw data set and more time for model identification analysis). Using both Monte Carlo simulations and empirical data, we derived and tested the generalized partial least squares (gPLS) framework, which can incorporate both PCA and ICA decompositions with computational efficiency. The gPLS method explicitly incorporates the experimental design to simplify the identification of characteristic spatiotemporal patterns. We performed parametric modeling studies of a blocked-design experiment under various conditions, including background noise distribution, sampling rate, and hemodynamic response delay. We used a randomized grouping approach to manipulate the degrees of freedom of PCA and ICA in gPLS to characterize both paradigm coherent and transient brain responses. Simulation data suggest that in the gPLS framework, PCA mostly outperforms ICA as measured by the receiver operating curves (ROCs) in SNR from 0.01 to 100, the hemodynamic response delays from 0 to 3 TR in fMRI, background noise models of Guassian, sub-Gaussian, and super-Gaussian distributions and the number of observations from 5, 10, to 20 in each block of a six-block experiment. Further, due to selective averaging, the gPLS method performs robustly in low signal-to-noise ratio (<1) experiments. We also tested PCA and ICA using PLS in a simulated event-related fMRI data to show their similar detection. Finally, we tested our gPLS approach on empirical fMRI motor data. Using the randomized grouping method, we are able to identify both transient responses and consistent paradigm/model coherent components in the 10-epoch block design motor fMRI experiment. Overall, studies of synthetic and empirical data suggest that PLS analysis, using PCA decomposition, provides a stable and powerful tool for exploration of fMRI/behavior data.     

Spectral spatiotemporal imaging of cortical oscillations and interactions in the human brain

This paper presents a computationally efficient source estimation algorithm that localizes cortical oscillations and their phase relationships. The proposed method employs wavelet-transformed magnetoencephalography (MEG) data and uses anatomical MRI to constrain the current locations to the cortical mantle. In addition, the locations of the sources can be further confined with the help of functional MRI (fMRI) data. As a result, we obtain spatiotemporal maps of spectral power and phase relationships. As an example, we show how the phase locking value (PLV), that is, the trial-by-trial phase relationship between the stimulus and response, can be imaged on the cortex. We apply the method to spontaneous, evoked, and driven cortical oscillations measured with MEG. We test the method of combining MEG, structural MRI, and fMRI using simulated cortical oscillations along Heschl's gyrus (HG). We also analyze sustained auditory gamma-band neuromagnetic fields from MEG and fMRI measurements. Our results show that combining the MEG recording with fMRI improves source localization for the non-noise-normalized wavelet power. In contrast, noise-normalized spectral power or PLV localization may not benefit from the fMRI constraint. We show that if the thresholds are not properly chosen, noise-normalized spectral power or PLV estimates may contain false (phantom) sources, independent of the inclusion of the fMRI prior information. The proposed algorithm can be used for evoked MEG/EEG and block-designed or event-related fMRI paradigms, or for spontaneous MEG data sets. Spectral spatiotemporal imaging of cortical oscillations and interactions in the human brain can provide further understanding of large-scale neural activity and communication between different brain regions.     

Hierarchical Bayesian estimation for MEG inverse problem

Source current estimation from MEG measurement is an ill-posed problem that requires prior assumptions about brain activity and an efficient estimation algorithm. In this article, we propose a new hierarchical Bayesian method introducing a hierarchical prior that can effectively incorporate both structural and functional MRI data. In our method, the variance of the source current at each source location is considered an unknown parameter and estimated from the observed MEG data and prior information by using the Variational Bayesian method. The fMRI information can be imposed as prior information on the variance distribution rather than the variance itself so that it gives a soft constraint on the variance. A spatial smoothness constraint, that the neural activity within a few millimeter radius tends to be similar due to the neural connections, can also be implemented as a hierarchical prior. The proposed method provides a unified theory to deal with the following three situations: (1) MEG with no other data, (2) MEG with structural MRI data on cortical surfaces, and (3) MEG with both structural MRI and fMRI data. We investigated the performance of our method and conventional linear inverse methods under these three conditions. Simulation results indicate that our method has better accuracy and spatial resolution than the conventional linear inverse methods under all three conditions. It is also shown that accuracy of our method improves as MRI and fMRI information becomes available. Simulation results demonstrate that our method appropriately resolves the inverse problem even if fMRI data convey inaccurate information, while the Wiener filter method is seriously deteriorated by inaccurate fMRI information.     

Estimation of the intrinsic dimensionality of fMRI data

A new method based on an autoregressive noise model of order 1 is introduced to the problem of detecting the number of components in fMRI data. Unlike current information-theoretic criteria like AIC, MDL, and related PPCA, which do not incorporate autocorrelations in the noise, the new method leads to more consistent estimates of the model order, as illustrated in simulated and real fMRI resting-state data.     

Estimating the global order of the fMRI noise model

One of the major issues in GLM-based fMRI analysis techniques is the presence of temporal autocorrelations in the residual signal after regression. A possible correction method is that of prewhitening, which fits an autoregressive (or other) model to the residual and uses the expected temporal autocorrelations of the model to transform the data and design matrix such that the residual becomes white noise. In this article, a method is introduced to estimate the global autoregressive model order of a data set, based on the residuals after regression. The proposed global standardized partial autocorrelation (SPAC) method tests whether the spatial profile of partial autocorrelations at a certain lag is random, and uses random field theory to account for the spatial correlations typical for fMRI data. It is tested both on synthetic and fMRI data, and is compared to two traditional techniques for model order estimation.     

Robust Bayesian General Linear Models

We describe a Bayesian learning algorithm for Robust General Linear Models (RGLMs). The noise is modeled as a Mixture of Gaussians rather than the usual single Gaussian. This allows different data points to be associated with different noise levels and effectively provides a robust estimation of regression coefficients. A variational inference framework is used to prevent overfitting and provides a model order selection criterion for noise model order. This allows the RGLM to default to the usual GLM when robustness is not required. The method is compared to other robust regression methods and applied to synthetic data and fMRI.     

Improved decoding of neural activity from fMRI signals using non-separable spatiotemporal deconvolutions

The goal of most functional Magnetic Resonance Imaging (fMRI) analyses is to investigate neural activity. Many fMRI analysis methods assume that the temporal dynamics of the hemodynamic response function (HRF) to neural activation is separable from its spatial dynamics. Although there is empirical evidence that the HRF is more complex than suggested by space-time separable canonical HRF models, it is difficult to assess how much information about neural activity is lost when assuming space-time separability. In this study we directly test whether spatiotemporal variability in the HRF that is not captured by separable models contains information about neural signals. We predict intracranially measured neural activity from simultaneously recorded fMRI data using separable and non-separable spatiotemporal deconvolutions of voxel time series around the recording electrode. Our results show that abandoning the spatiotemporal separability assumption consistently improves the decoding accuracy of neural signals from fMRI data. We compare our findings with results from optical imaging and fMRI studies and discuss potential implications for classical fMRI analyses without invasive electrophysiological recordings.     

Spatiotemporal clustering of fMRI time series in the spectral domain

We propose a new method for the analysis of functional magnetic resonance images (fMRI). The decision that a voxel v0 is activated is based not solely on the value of the fMRI signal at v0, but rather on the comparison of all time series s(v)(t) in a small neighborhood Nv0 around v0. Our approach explicitly takes into account the intrinsic spatiotemporal correlations that exist in the data. We focus on experimental designs with periodic stimuli, and therefore we can capture most of the features of the BOLD signal with a low dimensional subspace in the frequency domain. The presence of activated time series can be detected by partitioning the time series in this low dimensional space. Experiments with simulated data, and experimental fMRI data, demonstrate that our approach can outperform standard methods of analysis, such as the t-test.     

fMRI Signal Restoration Using a Spatio-Temporal Markov Random Field Preserving Transitions

In fMRI studies, Gaussian filtering is usually applied to improve the detection of activated areas. Such lowpass filtering enhances the signal to noise ratio. However, undesirable secondary effects are a bias on the signal shape and a blurring in the spatial domain. Neighboring activated areas may be merged and the high resolution of the fMRI data compromised. In the temporal domain, activation and deactivation slopes are also blurred. We propose an alternative to Gaussian filtering by restoring the signal using a spatiotemporal Markov Random Field which preserves the shape of the transitions. We define some interaction between neighboring voxels which allows us to reduce the noise while preserving the signal characteristics. An energy function is defined as the sum of the interaction potentials and is minimized using a simulated annealing algorithm. The shape of the hemodynamic response is preserved leading to a better characterization of its properties. We demonstrate the use of this approach by applying it to simulated data and to data obtained from a typical fMRI study.     

Bayesian fMRI time series analysis with spatial priors

We describe a Bayesian estimation and inference procedure for fMRI time series based on the use of General Linear Models (GLMs). Importantly, we use a spatial prior on regression coefficients which embodies our prior knowledge that evoked responses are spatially contiguous and locally homogeneous. Further, using a computationally efficient Variational Bayes framework, we are able to let the data determine the optimal amount of smoothing. We assume an arbitrary order Auto-Regressive (AR) model for the errors. Our model generalizes earlier work on voxel-wise estimation of GLM-AR models and inference in GLMs using Posterior Probability Maps (PPMs). Results are shown on simulated data and on data from an event-related fMRI experiment.     

Determining significant connectivity by 4D spatiotemporal wavelet packet resampling of functional neuroimaging data

An active area of neuroimaging research involves examining functional relationships between spatially remote brain regions. When determining whether two brain regions exhibit significant correlation due to true functional connectivity, one must account for the background spatial correlation inherent in neuroimaging data. We define background correlation as spatiotemporal correlation in the data caused by factors other than neurophysiologically based functional associations such as scanner induced correlations and image preprocessing. We develop a 4D spatiotemporal wavelet packet resampling method which generates surrogate data that preserves only the average background spatial correlation within an axial slice, across axial slices, and through each voxel time series, while excluding the specific correlations due to true functional relationships. We also extend an amplitude adjustment algorithm which adjusts our surrogate data to closely match the amplitude distribution of the original data. Our method improves upon existing wavelet-based methods and extends them to 4D. We apply our resampling technique to determine significant functional connectivity from resting state and motor task fMRI datasets.     

Wavelets and functional magnetic resonance imaging of the human brain

The discrete wavelet transform (DWT) is widely used for multiresolution analysis and decorrelation or "whitening" of nonstationary time series and spatial processes. Wavelets are naturally appropriate for analysis of biological data, such as functional magnetic resonance images of the human brain, which often demonstrate scale invariant or fractal properties. We provide a brief formal introduction to key properties of the DWT and review the growing literature on its application to fMRI. We focus on three applications in particular: (i) wavelet coefficient resampling or "wavestrapping" of 1-D time series, 2- to 3-D spatial maps and 4-D spatiotemporal processes; (ii) wavelet-based estimators for signal and noise parameters of time series regression models assuming the errors are fractional Gaussian noise (fGn); and (iii) wavelet shrinkage in frequentist and Bayesian frameworks to support multiresolution hypothesis testing on spatially extended statistic maps. We conclude that the wavelet domain is a rich source of new concepts and techniques to enhance the power of statistical analysis of human fMRI data.     

Multivariate dynamical systems models for estimating causal interactions in fMRI

Analysis of dynamical interactions between distributed brain areas is of fundamental importance for understanding cognitive information processing. However, estimating dynamic causal interactions between brain regions using functional magnetic resonance imaging (fMRI) poses several unique challenges. For one, fMRI measures Blood Oxygenation Level Dependent (BOLD) signals, rather than the underlying latent neuronal activity. Second, regional variations in the hemodynamic response function (HRF) can significantly influence estimation of causal interactions between them. Third, causal interactions between brain regions can change with experimental context over time. To overcome these problems, we developed a novel state-space Multivariate Dynamical Systems (MDS) model to estimate intrinsic and experimentally-induced modulatory causal interactions between multiple brain regions. A probabilistic graphical framework is then used to estimate the parameters of MDS as applied to fMRI data. We show that MDS accurately takes into account regional variations in the HRF and estimates dynamic causal interactions at the level of latent signals. We develop and compare two estimation procedures using maximum likelihood estimation (MLE) and variational Bayesian (VB) approaches for inferring model parameters. Using extensive computer simulations, we demonstrate that, compared to Granger causal analysis (GCA), MDS exhibits superior performance for a wide range of signal to noise ratios (SNRs), sample length and network size. Our simulations also suggest that GCA fails to uncover causal interactions when there is a conflict between the direction of intrinsic and modulatory influences. Furthermore, we show that MDS estimation using VB methods is more robust and performs significantly better at low SNRs and shorter time series than MDS with MLE. Our study suggests that VB estimation of MDS provides a robust method for estimating and interpreting causal network interactions in fMRI data.     

Spatiotemporal Brain Imaging of Visual-Evoked Activity Using Interleaved EEG and fMRI Recordings

Combined analysis of electroencephalography (EEG) and functional magnetic resonance imaging (fMRI) has the potential to provide higher spatiotemporal resolution than either method alone. In some situations, in which the activity of interest cannot be reliably reproduced (e.g., epilepsy, learning, sleep states), accurate combined analysis requires simultaneous acquisition of EEG and fMRI. Simultaneous measurements ensure that the EEG and fMRI recordings reflect the exact same brain activity state. We took advantage of the spatial filtering properties of the bipolar montage to allow recording of very short (125–250 ms) visual-evoked potentials (VEPs) during fMRI. These EEG and fMRI measurements are of sufficient quality to allow source localization of the cortical generators. In addition, our source localization approach provides a combined EEG/fMRI analysis that does not require any manual selection of fMRI activations or placement of source dipoles. The source of the VEP was found to be located in the occipital cortex. Separate analysis of EEG and fMRI data demonstrated good spatial overlap of the observed activated sites. As expected, the combined EEG/fMRI analysis provided better spatiotemporal resolution than either approach alone. The resulting spatiotemporal movie allows for the millisecond-to-millisecond display of changes in cortical activity caused by visual stimulation. These data reveal two peaks in activity corresponding to the N75 and the P100 components. This type of simultaneous acquisition and analysis allows for the accurate characterization of the location and timing of neurophysiological activity in the human brain.     

A state-space model of the hemodynamic approach: nonlinear filtering of BOLD signals

In this paper, a new procedure is presented which allows the estimation of the states and parameters of the hemodynamic approach from blood oxygenation level dependent (BOLD) responses. The proposed method constitutes an alternative to the recently proposed Friston [Neuroimage 16 (2002) 513] method and has some advantages over it. The procedure is based on recent groundbreaking time series analysis techniques that have been, in this case, adopted to characterize hemodynamic responses in functional magnetic resonance imaging (fMRI). This work represents a fundamental improvement over existing approaches to system identification using nonlinear hemodynamic models and is important for three reasons. First, our model includes physiological noise. Previous models have been based upon ordinary differential equations that only allow for noise or error to enter at the level of observation. Secondly, by using the innovation method and the local linearization filter, not only the parameters, but also the underlying states of the system generating responses can be estimated. These states can include things like a flow-inducing signal triggered by neuronal activation, de-oxyhemoglobine, cerebral blood flow and volume. Finally, radial basis functions have been introduced as a parametric model to represent arbitrary temporal input sequences in the hemodynamic approach, which could be essential to understanding those brain areas indirectly related to the stimulus. Hence, thirdly, by inferring about the radial basis parameters, we are able to perform a blind deconvolution, which permits both the reconstruction of the dynamics of the most likely hemodynamic states and also, to implicitly reconstruct the underlying synaptic dynamics, induced experimentally, which caused these states variations. From this study, we conclude that in spite of the utility of the standard discrete convolution approach used in statistical parametric maps (SPM), nonlinear BOLD phenomena and unspecific input temporal sequences must be included in the fMRI analysis.