Identifying Multisubject Cortical Activation in Functional MRI: A Frequency Domain Approach

: Functional magnetic resonance imaging (fMRI) has, since its description ﬁfteen years ago, become the most common in-vivo neuroimaging technique. FMRI allows the identiﬁcation of brain areas which are related to speciﬁc tasks, by statistical analysis of the BOLD (blood oxigenation level dependent) signal. Classically, the observed BOLD signal is compared to an expected haemodynamic response function (HRF) using a general linear model (GLM). However, the results of GLM rely on the HRF speciﬁcation, which is usually determined in an ad hoc fashion. For periodic experimental designs, we propose a multisubject frequency domain brain mapping, which requires only the stimulation frequency, and consequently avoids subjective choices of HRF. We present some computational simulations, which demonstrate a good performance of the proposed approach in short length time series. In addition, an application to real fMRI datasets is also presented.


Introduction
Since the description of the blood oxygenation level dependent signal (BOLD) by Ogawa et al. (1990), the number of studies based on functional magnetic resonance imaging (fMRI) has increased very rapidly.Tolias et al., 2005;Logothetis et al., 2002;Logothetis and Pfeuffer, 2004) showed that BOLD signal can be considered as an indirect measure of neuronal activity, reflecting local blood flow properties.In an fMRI session, many images are acquired at short time intervals, allowing the temporal monitoring of the relationships between the BOLD signal and stimulus presentation.Hence, BOLD signal reflects the activation level at a specific brain region (voxel).An important point is that the concept of stimuli is something related to conditions or active tasks performed by an individual.Thus, hands movements, speaking, calculating, are also considered as stimuli.
Structural magnetic resonance images are commonly acquired for medical diagnostics following lesions, tumors or strokes.In contrast, functional magnetic resonance imaging (fMRI) focuses on the observation of neuronal activity in the whole brain at short time intervals to establish the relationship of BOLD to a specific stimulus.The main advantages of fMRI analysis compared to EEG (eletroencefalography) and PET (positron emission tomography) are its noninvasive properties and also the high spatial resolution.The images in a fMRI session are acquired as multiple slices (Figure 1A) providing a 3-dimensional visualization of the whole brain (volume).In terms of data structure, each slice is a matrix composed by X × Y voxels, and each voxel represents a small brain area.Therefore, considering all the slices, the whole brain is represented by a 3-dimensional matrix of voxels.In a fMRI scanning session, several volumes are acquired through time (Figure 1B), resulting in T observations {y t , t = 1, .., T } for each voxel (a volume is a set of voxels).In conclusion, an fMRI dataset consists in many time series of BOLD, each one related to a small brain area.In a fMRI experiment, individuals lie inside a magnetic resonance scanner and are exposed to a sequence of stimuli.The evaluation of the temporal change in BOLD can then be used to infer neural dynamics, i.e., response to an experimental task.A general linear model (GLM, Graybill (1976)) is the most often used to identify brain regions activated by a particular task.SPM (statistical parametric mapping) is the most commonly used technique to describe the analysis results.In GLM analysis, a linear regression is performed considering the observed BOLD signal as dependent variable and the expected haemodynamic response function (HRF) as the regressor (Turner et al. (1998), Buchel et al. (1996), Friston et al. (1998,b)).The HRF is the expected variation or response of BOLD signal after a stimulus presentation (examples of possible choices for HRF are plotted in the left side of Figure 2).The model is given by where Y t is the value of BOLD signal observed at time t, X t is the expected predicted response assuming stimulation (HRF) and t is a random error.As the dependent variable is obtained directly from the data, the first issue in this analysis is the choice of the HRF regressor, which is subjective.Assuming that the sequence and times of stimulus are known, there are many suggestions for an appropriate HRF (Buchel et al. (1998), Friston et al. (2000)) based on haemodynamic delay, e.g., Gamma functions, Volterra kernels, and non-parametric smoothing.
Finally, the decision if a certain voxel was or not activated by the stimulus, is reached by evaluating the statistical significance of activation obtained using a Wald statistics on HRF coefficient (β).The null hypothesis of no activation is given by H 0 : β = 0.In other words, the GLM analysis provides a decision whether an observed BOLD signal is similar to an expected variation under stimulation (HRF).
Here, we propose a multisubject activation mapping in frequency domain for cases of periodic stimulation.In this approach, HRF specification is not necessary, avoiding subjective or incorrect choices.This paper is structured as follow: the proposal is presented in section 2. Section 3 contains simulation results related to statistical power and small sample approximations.Finally, an application of the proposed approach to a real fMRI dataset involving a motor task is presented in section 4.

Time series analysis
A time series is a set of observations y t , each one recorded at a specified time t.The obvious correlation introduced by the sampling of adjacent points in time can often restrict the applicability of many conventional statistical methods, which rely on the assumption that these adjacent observations are independent and identically distributed.Time series analysis is the systematic approach by which one goes about answering the mathematical and statistical questions posed by these temporal correlations.In general, we observe one or more realizations {y t , t = 1, • • • , T } of a stochastic process {Y t }, describe its properties and make inferences.In time series analysis, there are two common approaches: time domain and frequency domain analysis.An interesting focus of several studies has been the class of second order stationary processes.The main properties of these processes are: a) the time-invariant unconditional expectation and b) variance and covariance between observations in different time points depending only on the delay between them (Brillinger, 1980).
In the time domain, the correlation structure of a second order stationary process is described by the autocovariance function (acf) γ y (k) defined as where µ is the unconditional expectation of the process.This function measures the linear dependence between two points separated by lag k in a series.
In the frequency domain, the correlation structure is represented by the spectral density, where λ is measured in cycles per unit of time.The spectral density describes the properties of the process in terms of periodic components at different frequencies in any given realization.In addition, the spectral density may also be interpreted as a variance decomposition in orthogonal components (sines and cosines).

Spectral analysis of periodic stimulus designs
In the analysis of fMRI voxel time series, considering for periodic designs, the stimulation frequency is defined as the fundamental frequency of activation (λ a ).Our aim is the identification of voxels which have BOLD signal oscillating at the stimulus frequency.As in the case of GLM, the analysis is performed in order to detect responses to the stimulus.We consider the model where y t represents the BOLD signal, K is the number of components, R j is the amplitude, λ j is the frequencies of interest and t is a white noise series.The variance of a voxel time series, {y t , t = 1, • • • , T }, attributable to an oscillation of frequency λ j is obtained through the spectral density (1) with λ = λ j , which can be estimated by the periodogram defined as where d y (λ j ) is the discrete Fourier transform of y 1 , . . ., y T , at the Fourier frequencies and defined as (2.4) In practice, the discrete Fourier transform is obtained using the Fast Fourier Transform (fft) algorithm, which is computationally efficient.Most of the frequencies will contain information solely about the correlation structure of the underlying stochastic process at the chosen voxel.The value of the periodogram at each frequency represents the amount of time series variance related to this frequency (power).Figure 2 shows some haemoresponse functions and their respective periodograms.The amount of variance explained by a frequency is defined as the power (or energy) in this frequency.Note that the power (Eq.4) at the fundamental frequency is so large, that visually, it makes the power in other frequencies seems to be zero.
The large value of the periodogram I y (λ j ) at the fundamental frequency is indicative of response to the stimulus.Hence, for a periodic design, we are only interested in the spectral density at a fourier frequency λ a (the fundamental frequency of activation).
The asymptotic sampling properties of the periodogram are well-known.For stationary series, (I) I y (λ j ) and I y (λ k ) are asymptotically independent, for all j = k; independently, where χ 2 2 denotes a chi-squared random variable with 2 degrees of freedom (Shumway and Stoffer (2000), Brillinger (1980Brillinger ( , 1981))), and D means convergence in distribution.

Testing for a response to the stimulus
The periodogram provides a baseline to test the null hypothesis of no power peak in a frequency of interest.We define the ratio statistic at the fundamental frequency of activation (frequency of stimulus), λ a , for each voxel time series as as to obtain a test statistic for significant activation.Large value of W a indicates a large effect at the fundamental frequency, and consequently, that the BOLD signal of this brain area is correlated with the stimulus presentation.It is important to highlight that, theoretically, the W a statistic is adequate only in cases where the energy of the signal is concentrated solely in one frequency.The power of the test is reduced if there are two or more peaks.However, in some cases, the power at the frequency of interest is so large compared to the other peaks, that an approximation to the asymptotic distribution of W a statistic holds.From (I) and (II), under the hypothesis of no activation, the asymptotic distribution of W a is given by 2W when T → ∞ (Shumway and Stoffer (2000), Brillinger (1980Brillinger ( ), 1981))).Analogously, for a multisubject analysis considering N individuals, we have (2.8) Hence, we reject the hypothesis of no activation for large energy in the fundamental frequency of activation.Nevertheless, some fMRI time series are autocorrelated.In these cases, pre-whitening filters which preserve the time series periodicity should be applied.Periodicity detection is the main concern in the analysis of fMRI cyclical experiments, hence trends or other components are not important and may be discarded.

Simulations
Following the derivation of the asymptotic distribution of W * a statistic, the validity of approximation in small samples need to be evaluated.Further, the effects of the number of subjects, energy in the frequency of interest or time series length are also unknown.Hence, in this section, we present the results of Monte  Carlo simulations related to these questions.All simulations were performed using the R Statistical Software (www.r-project.org).
Firstly, consider the case of a white noise time series.We simulated 10000 Gaussian white noise processes of length 100 for 6 subjects in order to empirically estimate the W * a probability density function under the null hypothesis (frequency λ a = 0.05).The histogram, theoretical and estimated (Kernel) densities are presented in Figure 3.Note that the estimated and theoretical densities are similar, indicating a good asymptotic approximation even for short length time series.
Focusing on the power evaluation, consider the following model: where T is the time series length, n = 1, .., N is an index representing each subject, λ a is the fundamental frequency of activation, R a is a coefficient representing the energy in frequency λ a and t is a Gaussian white noise.Ten thousand simulations were performed for each evaluation and the size of the tests is α = 0.05.The effect in the power of the test by increasing the energy in the fundamental frequency of activation R a for N = 6, T = 100 and λ a = 0.05 is shown in Figure 4A.The effects of sample length (T ) and number of subjects (N ) are presented in Figures 4B and 4C, respectively.The simulations provide evidence of a satisfactory performance of the proposed statistical test.Figure 3 shows that W * a statistic asymptotic distribution is a reasonable approximation even for short length time series.This result is useful, as that in many cases of periodic designs, the fMRI time series have short length in order to avoid habituation effects.Further, Figure 4 shows a reasonable increase in the test power, as the energy at the frequency of interest (R a ), the time series length (T ) and the number of subjects (N ) increase, respectively.

Application to fMRI Data
In this section, we illustrate the usefulness of the proposed approach in the analysis of a fMRI motor dataset with a block design.An experiment with block design consists on segmented presentation of stimulus.In other words, each type of stimulus (e.g., sounds, fingertap, lights, resting, etc) is presented during long time intervals.In a session (run), the subjects are exposed to these stimuli, which may occur alternately or periodically.
Seven right handed healthy volunteers (4 male and 3 female, 36 to 75 yearsold) participated in this study.The datasets were acquired at the HC Radiology Institute -University of Sao Paulo (Brazil).The subjects performed a simple motor task: finger tapping movements of right hand.A volume contains observations of BOLD signal of each voxel (regions) in the whole brain.Each volume was composed of 15 slices acquired in a 1.5Tesla Signa LX (GE, Milwaukee, USA) magnetic resonance scanner.For each subject, one hundred volumes were acquired in different time points, resulting in BOLD signals of length 100 per voxel (T = 100).The sampling of BOLD signal occurs at fixed intervals called TR (time of repetition), which is in general specified in miliseconds (ms).Thus, for each region in the brain the BOLD signal of seven individuals (N = 7) was monitored.
The frequency domain multisubject periodicity test was then applied to these time series, resulting in statistical measures of association between stimulation and activity of each brain region.The visualization of the regions with significant periodicity (red spots) is overlayed on structural brain images.
Each block design session consisted of 5 cycles with two conditions each (30 seconds of active finger movements and 30 seconds of rest) in response to a visual cue.The BOLD sampling interval is TR=3000ms.Notice that the sequence of stimulations is periodic (rest followed by motor task), with fundamental frequency of activation, λ a , corresponding to 0.05TR/s.Here, the main interest is the identification of brain voxels which show a significant periodicity in the stimulation frequency, indicating that these regions are activated by the stimulus.
The images were preprocessed considering motion realignment, slice time correction (SPM2, http://www.fil.ion.ucl.ac.uk/spm/ ) and spatial smoothing.This first stage is necessary to remove scanning artifacts.Individuals have different brain sizes and shapes, thus, warping the volumes to a common brain template (stereotatic space of Talairach and Tornoux (1988) is necessary in multisubject analysis.
The BOLD signal was detrended (polinomial of order 2) and also pre-whitened considering an AR(1) model in order to remove autocorrelation from the data (Marchini and Smith (2003), Woolrich et al. (2001)).This filter does not change the periodicity of the signal.Finally, brain activation maps were obtained using the proposed multisubject periodicity test.The red spots show the brain areas activated during the fingertap task (Figures 5 and 6).According to activation maps (p-value< 10 −4 , less than 1 expected false activated voxel per slice), we found more activated clusters in left primary motor cortex (LM1), left primary sensitive area (LS1) and supplementary motor area (SMA).Small clusters were found in right primary motor cortex (RM1).The left primary motor area is classicaly involved in right hand movements and SMA is often related to movements programming.The activation in LS1 is expected as fingertap tasks involve also finger touching.The interhemispheric communication is thought to be responsible for a reduction of the oxygenation, and consequently BOLD effect in this area, probably in order to give priority to the motor learning in the RM1 (Li et al. (1996), Tinazzi and Zanette (1998), Allison et al. (2000)).This pattern of activation has also been found in classical fMRI motor studies with finger-tapping, which is a simple task but with motor and sensory components (Bandettini et al. (1992), Kwong et al. (1992), Allison et al. (2000)).Figure 7 presents the average time series of all subjects corresponding to an active voxel in primary motor area.
In conclusion, we found that all activated regions detected using the proposed approach are in total agreement with the literature related to this topic.This provides preliminary evidence for the reliability of the method.

Conclusion
The number of neuroscientific studies based on fMRI has been increasing rapidly.However, the quality of results relies on the choice of haemodynamic response function (HRF).In this paper, we propose a frequency domain multisubject approach, which is based solely in stimulus periodicities, avoiding subjective specification of HRF.The power and usefulness of the new approach in small samples are evaluated by simulations.We also present an application to fMRI dataset involving a motor experiment.The results provide evidence of the reliability of the proposed approach in normal subjects.Future works involve the application of this approach to detect differences of activation between healthy and unhealthy groups.

Figure 1 :
Figure 1: A: A brain volume of multiple slices acquired at one time point.B: A fMRI time series is a set of volumes (BOLD in each voxel) observed at different time points.

Figure 2 :
Figure 2: Some haemodynamic response functions and their respective periodograms.

Figure 3
Figure 3: W * a estimated probability density function (full line).The theoretical function is represented by the dashed line.

Figure 4 :
Figure 4: a) Test power as a function of energy in the fundamental frequency of activation(R a ).b) Test power as a function of time series length(T ).c) Test power as a function of the number of subjects(N )

Figure 5 :
Figure 5: Block design data analysis: the maps show brain slices and red spots indicating the areas (voxels) related to the motor stimulus (radiological notation).

Figure 6 :Figure 7 :
Figure 6: Block design data analysis: the figure shows a 3D-view of activated voxels (red spots) in different angles.