Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3696586

Formats

Article sections

Authors

Related links

Med Image Comput Comput Assist Interv. Author manuscript; available in PMC 2013 June 30.

Published in final edited form as:

Med Image Comput Comput Assist Interv. 2009; 12(0 1): 1009–1017.

PMCID: PMC3696586

NIHMSID: NIHMS136681

See other articles in PMC that cite the published article.

The standard general linear model (GLM) for rapid event-related fMRI design protocols typically ignores reduction in hemodynamic responses in successive stimuli in a train due to incomplete recovery from the preceding stimuli. To capture this adaptation effect, we incorporate a region-specific adaptation model into GLM. The model quantifies the rate of adaptation across brain regions, which is of interest in neuroscience. Empirical evaluation of the proposed model demonstrates its potential to improve detection sensitivity. In the fMRI experiments using visual and auditory stimuli, we observed that the adaptation effect is significantly stronger in the visual area than in the auditory area, suggesting that we must account for this effect to avoid bias in fMRI detection.

Rapid event-related (ER) functional magnetic resonance imaging (fMRI) is one of the most popular imaging methods in cognitive neuroscience. In the rapid ER fMRI studies, individual stimuli are presented every few seconds or faster. Although less efficient for localizing activation, rapid ER fMRI has several advantages over the traditional block design, including the ability to randomize trial types and to sort data based on behavioral responses.

The standard analysis for rapid ER fMRI models activation as a linear system [2, 5, 9]; the hemodynamic response to multiple input stimuli is assumed to be a superposition of the responses to individual stimuli. This approach estimates the impulse response function, also known as the hemodynamic response function (HRF), of this linear system via de-convolution, and compares the estimates to the null hypothesis or to estimates from other experimental conditions.

fMRI signals commonly do not comply with the linear assumption. Independent studies have demonstrated a substantial adaptation effect in the hemodynamic response [1, 11, 12, 16, 17], i.e., if two stimuli are presented within the adaptation period, the amplitude of the response to the second stimulus is reduced. Furthermore, the adaptation effect strengthens as the inter-stimulus interval (ISI) decreases. Several studies demonstrated that when a pair of visual stimuli is presented less than 1 sec apart, the response amplitude to the second stimulus is approximately 55% of that to the first stimulus, with recovery to 90% at a 6 sec ISI [12, 16, 17]. This evidence suggests that the adaptation effect must be modeled in the analysis, especially when stimuli are presented frequently.

The adaptation effect is expected to vary spatially due to differences in neural and hemodynamic properties of functional areas in the brain [1, 13, 16]. While physiological mechanisms for adaptation are not clearly understood, it is still useful to model it for the purposes of improving detection.

Previous studies of the adaptation effect separated detection and adaptation modeling [12, 13, 15–17], fitting the adaptation model to the estimated HRF obtained using the standard general linear model (GLM) [9]. This approach ignores the trial-to-trial variation. Work by Buxton *et al.* [4] introduced the biophysical balloon model for fMRI signals where the adaptation effect is captured through interactions among blood flow, blood volumes, and de-oxyhemoglobin content, instead of an explicit interaction between stimuli. Friston *et al.* [10] proposed a statistical model using the Volterra kernels to capture interaction between stimuli. The interaction can be efficiently estimated and statistically examined via the *F*-test. However, the physiological interpretation of the model parameters is challenging, since the model treats the stimuli symmetrically, effectively ignoring the causal nature of the adaptation effect.

In this work, we extend the basic GLM by incorporating a region-specific model of adaptation. In addition to the stimulus onset, our design matrix also depends on the ISIs between stimuli via a single-parameter exponential function. Specifically, this model captures the decrease in the magnitude of the hemodynamic response if the time interval to the preceding stimuli is short. In other words, we only model causal interactions among stimuli, in contrast to the bidirectional interaction model in [10]. By combining detection and adaptation modeling, the proposed method takes into account trial-to-trial variation. It is expected that the adaptation effect strengthens when more stimuli are presented prior to the current stimulus. We summarize this effect from multiple stimuli through a multiplicative model. This extension allows for a more flexible choice of an experimental paradigm in contrast to previous fMRI adaptation studies which were restricted to presentations of stimulus pairs [12, 13, 16, 17].

We jointly estimate the decay parameter of the exponential function for each region and the HRF for each location in the brain. The estimated parameter of the exponential function reflects the length of the adaptation period for the corresponding region, and the estimated HRF indicates the activation status of the corresponding location. Our experimental results demonstrate a significant improvement in detection sensitivity and confirm previously known adaptation phenomena in the sensory systems.

The univariate GLM [9] assumes that the fMRI signal * y_{n}* at location

$${\mathit{y}}_{n}=\mathbf{X}{\mathit{\beta}}_{n}+\mathbf{\Phi}{\mathit{\alpha}}_{n}+{\mathit{\epsilon}}_{n},$$

(1)

where **X** is the design matrix, constructed based on the experimental protocol. Columns of matrix **Φ** = [ϕ_{1}, , ϕ_{R}], often in the form of low-order polynomials, model the protocol-independent factors such as cardiac activity and breathing. The measurement noise **ε**_{n} can be modeled as white Gaussian noise or as colored noise with an auto-regressive (AR) structure [3]. Vectors **β**_{n} and **α**_{n} are the corresponding coefficients of the protocol-dependent and the protocol-independent signals, respectively. To determine the activation status at location *n*, one compares the estimated protocol-dependent coefficient _{n} to the null hypothesis or to the corresponding estimates for other experimental conditions.

Without loss of generality, we assume a single type of stimulus. The matrix form of GLM represents the convolution of experimental protocol and HRF:

$${y}_{n}(t)={\displaystyle {\sum}_{k=1}^{K}h(t-{s}_{k})\star {\beta}_{n}(t)+}{\displaystyle {\sum}_{r=1}^{R}{\alpha}_{\mathit{\text{rn}}}{\varphi}_{r}(t)+{\epsilon}_{n}(t)},$$

(2)

where ** s** = [

The above model fails to capture the fact that previous stimuli can decrease the hemodynamic response to subsequent stimuli if the recovery period is longer than the ISI presented. Therefore, we incorporate an adaptation model into the standard GLM by introducing a damping weight *w _{k}* for each stimulus. Due to the regionally varying neuronal and vascular architecture of the brain, we parameterize the weight

$${y}_{n}(t)={\displaystyle {\sum}_{k=1}^{K}[{w}_{k}(\mathit{s};{\theta}_{m})h(t-{s}_{k})]\star {\beta}_{n}(t)}+{\displaystyle {\sum}_{r=1}^{R}{\alpha}_{\mathit{\text{rn}}}{\varphi}_{r}(t)+{\epsilon}_{n}(t)}.$$

(3)

In the matrix notation this equation reads

$${\mathit{y}}_{n}=\tilde{\mathbf{X}}({\theta}_{m}){\mathit{\beta}}_{n}+\mathbf{\Phi}{\mathit{\alpha}}_{n}+{\mathit{\epsilon}}_{n}\forall n\in {V}_{m},$$

(4)

where *V _{m}* is the set of locations in region

We combine the adaptation effects for stimulus *k* using a multiplicative exponential model, ranging between zero and one:

$${w}_{k}(\mathit{s};{\theta}_{m})={\displaystyle {\prod}_{j=1}^{k-1}(1-{e}^{-{\theta}_{m}({s}_{k}-{s}_{j})}).}$$

(5)

The exponential decay parameter θ_{m} models the length of the adaptation period. A larger value of θ_{m} indicates a weaker adaptation effect, or a shorter period for the region to recover. The multiplicative nature of the model reflects the fact that multiple preceding stimuli can affect the amplitude of the response to a particular stimulus. In other words, the brain response is modeled as a Markov process of order K.

We estimate the unknown parameters {θ_{m}, {**β**_{n}, **α**_{n}}_{nεVm}} in Eq. (4) by minimizing the sum of squares of the residual errors for each region independently. For a given θ_{m}, the optimal linear parameters _{n} and _{n} can be found in a closed-form:

$$[{\widehat{\mathit{\beta}}}_{n},{\widehat{\mathit{\alpha}}}_{n}]={({\mathbf{H}}^{\mathrm{T}}({\theta}_{m})\mathbf{H}({\theta}_{m}))}^{-1}{\mathbf{H}}^{\mathrm{T}}({\theta}_{m}){\mathit{y}}_{n},$$

(6)

where **H**(θ_{m}) = [(θ_{m}) **Φ**]. Substituting Eq. (6) into the expression for the residual error, we obtain the optimal ${\widehat{\theta}}_{m}^{*}$:

$${\widehat{\theta}}_{m}^{*}=\underset{\theta}{\text{arg min}}\phantom{\rule{thinmathspace}{0ex}}{\displaystyle {\sum}_{n\in {V}_{m}}{\Vert (\mathbf{I}-{({\mathbf{H}}^{\mathrm{T}}(\theta )\mathbf{H}(\theta ))}^{-1}{\mathbf{H}}^{\mathrm{T}}(\theta )){\mathit{y}}_{n}\Vert}^{2}.}$$

(7)

With the proposed adaptation model, Eq. (7) is a nonlinear function of a single scalar parameter θ_{m}. We can simply exhaustively search for the parameter value within a specified range. We then obtain the optimal values ${\widehat{\mathit{\beta}}}_{n}^{*}\phantom{\rule{thinmathspace}{0ex}}\text{and}\phantom{\rule{thinmathspace}{0ex}}{\widehat{\mathit{\alpha}}}_{n}^{*}$ by substituting ${\widehat{\theta}}_{m}^{*}$ into Eq. (6).

Due to the nonlinearity of the model, the true value of θ is needed to compute the covariance of the estimate, $\text{Cov}({\widehat{\mathit{\beta}}}_{n}^{*})$. Since θ is not known in real experiments, we approximate $\text{Cov}({\widehat{\mathit{\beta}}}_{n}^{*})$ with $\text{Cov}({\widehat{\mathit{\beta}}}_{n}^{*};{\widehat{\theta}}_{m}^{*})$. Hence, the statistic ${{\widehat{\mathit{\beta}}}_{n}^{*}}^{\mathrm{T}}{\text{Cov}}^{-1}({\widehat{\mathit{\beta}}}_{n}^{*};{\widehat{\theta}}_{m}^{*}){\widehat{\mathit{\beta}}}_{n}^{*}/{N}_{{\mathit{\beta}}_{n}}$, where *N*_{βn} is the number of regression coefficients in **β**_{n}, does not follow a known probability distribution under the null hypothesis, in contrast to the *F*-distribution in the standard GLM analysis. This represents a challenge in testing significance similar to GLM with AR noise modeling [3]. The exact statistical test can be achieved with Markov Chain Monte Carlo simulation, but it is too computationally intensive for a standard analysis procedure. We will see in the next section that comparing the values of the statistic across locations provides insight into the adaptation effect. Developing efficient methods for assessing statistical significance of the statistic is clearly a direction for future research.

To summarize, by accounting for the adaptation effects, we obtain a more accurate estimate of the HRF which leads to more accurate detection results. The estimates of the adaptation parameter θ promise to provide an insight into the neuronal and vascular architecture across different brain regions. In the following, we refer to our approach as GLMA (GLM with adaptation).

Due to the lack of ground truth in real experiments, we first study GLMA’s sensitivity to noise and parameter settings using simulated data. We then compare GLMA to GLM using human fMRI data from a visual-auditory study.

In both simulations and analysis of human fMRI data, we constrain the detection to the cortex and define different brain regions based on the cortical folding patterns, obtained with the FreeSurfer software [6, 8], 35 regions per hemisphere. Moreover, since the adaptation weight *w _{k}* recovers exponentially with respect to ISI, we only consider the stimulus interactions within a 16 sec window.

Since our model is estimated for each region separately, it is sufficient to study the performance of the model for a range of values of θ using data in a single region. For each value of θ, we generated 100 activation time courses by convolving a two-gamma function [14] with a train of stimuli whose onset times were generated randomly (mean ISI=4.0 sec, std=3 sec). The adaptation effect was modeled according to Eq. (5). We also generated 100 time courses without activation. We then added i.i.d. Gaussian noise to these 200 time courses, and repeated the simulation procedure 50 times for each value of θ. We chose noise levels to be within the typical SNR range of real fMRI data.

We separately processed the data set using GLM and GLMA. In both cases, HRF was modeled with the two-gamma function used in the simulation. Thus, β is a scalar in this case. Based on the activation statistic (*)^{2}/Var (*), we evaluated the detection accuracy of both methods at two false positive rates, 5 × 10^{−4} and 5 × 10^{−2}, as shown in Fig. 1(a,b). As expected, a better SNR in the data allows for a higher detection rate over the range of θ we examined. When θ is larger than 0.5, there is essentially no adaptation effect present in the data. Hence, the performance of the two methods is almost identical. The adaptation effect strengthens as θ decreases. We can clearly see an up-to 80% higher detection rate achieved by modeling the effect.

The two panels on the left compare the true positive rates for GLM (solid) and GLMA (dashed) as a function of θ, for two false positive rates and four SNR values. The rightmost panel shows the estimates * of the adaptation parameter **...**

We also investigated the robustness of the estimation of the adaptation parameter θ. As illustrated in Fig. 1(c), the estimates * closely match the simulation. As θ decreases, the response to the subsequent stimuli is very small for the chosen mean ISI, and * starts to deviate from the true value θ for noisy data, i.e., SNR= −10 dB.

We illustrate the application of the proposed method in a visual-auditory fMRI study. In this experiment, three healthy right-handed subjects were presented with visual-auditory stimuli in a random rapid ER fMRI paradigm in three runs, with mean ISI of 1.5 sec (std=0.7 sec), 3.0 sec (std=1.3 sec), and 6.0 sec (std=2.1 sec), respectively. To keep subjects engaged throughout the experiment, they were asked to respond to a rare target stimulus by pressing a button. The fMRI data were acquired using a 3T Siemens Trio scanner (TR 1.15 sec, TE=30 msec, flip angle 90 degree, 3.1×3.1×4 mm^{3}). Structural T1-weighted MRI scans of the subjects were acquired with a 1.5T Siemens Avanto scanner. The anatomical images were processed with the FreeSurfer software [6, 8]. Individual subject functional scans were morphed through a spherical mapping into an atlas constructed with 40 subjects [7].

We applied GLM and GLMA to data combined from all three runs (Fig. 2 left) and data combined from the first quarter of each of the three runs (Fig. 2 right). Since the statistics across methods are not directly comparable, and the ground truth activation is not known, we select top 5% of vertices in a hemisphere with the highest statistics, and visually evaluate the results to assess the importance of modeling the adaptation effect. We emphasize that in order to develop valid detectors, further work is required in building statistical tests for assessing significance of the detected activations.

The thresholded statistical parametric maps for three subjects using GLM with full length data (first row) and with 1/4 of the data (third row), as well as using GLMA with full-length data (second row) and with 1/4 of the data (fourth row).

Fig. 2 shows that the two methods provide similar detection results in the auditory area. However, the visual region detected by GLM has a smaller spatial extent than the corresponding results using GLMA. Without modeling the adaptation effects, the activation statistics in the visual area are smaller than those in the auditory area. Therefore, many activations on the visual cortex will be missed if a single threshold is applied to the whole brain. Furthermore, compared to GLM, our detection results with shorter length data are more similar to the results using full-length data, indicating improved robustness.

For the selected auditory and visual areas, we present the time required for the attenuation coefficient *w _{k}* to recover to 90%, i.e.,

To further validate our method, we also applied GLM, with HRF modeled as a two-gamma function, separately to each of the three runs. Fig. 4 shows the average estimated HRFs for the 50 most active vertices in each of the three selected regions for Subject 2. We can clearly see that the lingual area exhibits the strongest adaptation effect (*t*_{90} = 5.8 sec in Fig. 3), indicated by substantial differences in response magnitude to stimuli presented with different mean ISI conditions. The difference in response magnitudes to the three conditions is smaller in the lateral-occipital area (*t*_{90} = 3.8 sec), reflecting a weaker adaptation. On the other hand, the estimated HRFs across different ISIs are almost identical in the superior-temporal area (*t*_{90} = 2.3 sec). That means the most active locations in the superior-temporal area can almost fully recover in about 1.5 sec. The results of the HRF analysis for separate conditions agrees with our estimates of the adaptation effects in Fig. 3.

We proposed and demonstrated a novel method for modeling the adaptation effects in fMRI. Experimental results indicate a significant improvement in detection sensitivity. The proposed method also quantifies the adaptation effects across brain regions, providing insight into neuronal and vascular organization of the brain. The current adaptation model focuses on the ISIs. We plan to extend it to include information about the amplitude of the response to previous stimuli, since the adaptation effect is expected to be more pronounced immediately after a strong response than after a weak response. Our framework can be readily modified to include different functional forms of the magnitude changes due to adaptation. We will explore alternative functions, such as the sigmoid function, in adaptation modeling in future work.

This work was supported in part by NIH R01 NS048279, NIH R01 EB006385, NIH NIBIB NAMIC U54-EB005149, NIH NCRR NAC P41-RR13218, NIH NCRR P41-RR14075 grants, NSF CAREER Award 0642971, Sigrid Juselius Foundation, Academy of Finland, and Finnish Cultural Foundation. Wanmei Ou is partially supported by the PHS training grant DA022759-03.

1. Birn RM, Saad ZS, Bandettini PA. Spatial heterogeneity of the nonlinear dynamics in the fMRI BOLD response. NeuroImage. 2001;14:817–826. [PubMed]

2. Boynton GM, Engel SA, Glover GH, Heeger DJ. Linear systems analysis of functional magnetic resonance imaging in human V1. J. NeuroSci. 1996;16:4207–4221. [PubMed]

3. Burock M, Dale A. Estimation and detection of event-related fMRI signals with temporally correlated noise: a statistically efficient and unbiased approach. Hum. Brain Mapp. 2000;11:249–260. [PubMed]

4. Buxton RB, Wong EC, Frank LR. Dynamics of blood flow and oxygenation changes during brain activation: the balloon model. Magn. Reson. Med. 1998;39:855–864. [PubMed]

5. Dale A, Buckner R. Selective averaging of rapidly presented individual trials using fMRI. Hum. Brain Mapp. 1997;5:329–340. [PubMed]

6. Dale AM, Fischl B, Sereno MI. Cortical surface-based analysis: I. segmentation and surface reconstruction. NeuroImage. 1999;9:179–194. [PubMed]

7. Desikan RS, Segonne F, Fischl B, Quinn BT, Dickerson BC, Blacker D, Buckner RL, Dale AM, Maguire RP, Hyman BT, Albert MS, Killiany RJ. An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. NeuroImage. 2006;31:968–980. [PubMed]

8. Fischl B, Sereno M, Dale AM. Cortical surface-based analysis: II. inflation, flattening, and a surface-based coordinate system. NeuroImage. 1999;9:195–207. [PubMed]

9. Friston KJ, Holmes AP, Worsley KJ, Poline J-B, Frith CD, Frackowiak R. Statistical parametric maps in functional imaging: a general linear approach. Hum. Brain Mapp. 1994;2:189–210.

10. Friston KJ, Mechelli A, Turner R, Price CJ. Nonlinear responses in fMRI: the Balloon model, Volterra kernels, and other hemodynamics. NeuroImage. 2000;12:466–477. [PubMed]

11. Heckman GM, Bouvier SE, Carr VA, Harley EM, Cardinal KS, Engel SA. Nonlinearities in rapid event-related fMRI explained by stimulus scaling. NeuroImage. 2007;34:651–660. [PMC free article] [PubMed]

12. Huettel S, McCarthy G. Evidence for a refractory period in the hemodynamic response to visual stimuli as measured by MRI. NeuroImage. 2000;11:547–553. [PubMed]

13. Huettel S, McCarthy G. Regional differences in the refractory period of the hemodynamic response: an event-related fMRI study. NeuroImage. 2001;14:967–976. [PubMed]

14. Jezzard P, Matthews PM, Smith SM. Functional MRI: an introduction to methods. Oxford university press; 2002.

15. Robson MD, Dorosz JL, Gore JC. Measurements of the temporal fMRI response of the human auditory cortex to trains of tones. NeuroImage. 1998;7:185–198. [PubMed]

16. Soon CS, Venkatraman V, Chee MW. Stimulus repetition and hemodynamic response refractoriness in event-related fMRI. Hum. Brain Mapp. 2003;20:1–12. [PubMed]

17. Zhang N, Zhu XH, Chen W. Investigating the source of BOLD nonlinearity in human visual cortex in response to paired visual stimuli. NeuroImage. 2008;43:204–212. [PMC free article] [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |