Creating a combined HRF tmap

In our particular setting, the analysis of the EEG-fMRI data consists in an event-related fMRI paradigm, where the event consists in the occurrence of epileptic events, such as interictal spikes, measured with the EEG in the scanner. Compared to cognitive fMRI studies, the paradigm is not completely controlled, as epileptic events occur unpredictably. After the identification of those events by an expert electroencephalographer, the usual general linear model approach is used to analyze fMRI data.

In a typical linear model analysis, a single HRF consisting of a positive gamma function followed by a smaller negative gamma function (the undershoot) is used. A commonly used HRF is the Glover HRF (Glover 1999). Its first gamma function peaks at 5.4 seconds, with a full width at half maximum (FWHM) of 5.2 seconds. The second gamma function (the undershoot) peaks at 10.8 seconds, has a FWHM of 7.35 seconds and an amplitude of 35% of the first gamma function. This model has widely been shown to be successful at detecting fMRI activations and has been used extensively. However, the possibility of varying local responses and inter-subject variability may mean that an analysis using only one HRF may not be optimal.


Our group examined patient specific HRFs (Kang et al 2003) . Patients were analyzed using the Glover HRF. Once activations were obtained, a patient specific local HRF was calculated. Data from this patient were then reanalyzed with this HRF. Activation in the original region predictably increased, but new regions of activation were also uncovered. However, this method has a fundamental drawback. If the initial regions of activation are defined by an analysis using the Glover HRF, the HRFs for those regions will tend to be similar to the Glover HRF. If a region has an HRF response that does not well match the Glover HRF, it is less likely to be significantly activated or included in the voxels used to calculate the patient specific HRF. To overcome this limitation, we developed a multi-HRF approach.


In this approach, we conduct 4 separate analysis of each data set (Bagshaw et al 2004). Each analysis uses a different HRF, consisting of a single positive gamma function peaking at 3, 5, 7, or 9 seconds after event onset, with a FWHM of 5.2 seconds. We do not model an undershoot in this analysis because the undershoot is a more variable component of the HRF, sometimes absent and variable in amplitude and time to peak. A separate analysis is conducted using each of the four models. Using fMRIstat, this means fmrilm should be run four times for each scan for each subject. In SPM, the analysis should be run four times, each time using a different model HRF.

The four models used in our multi HRF method.

The four models used in our multi HRF method.

fMRIstat instructions

The use of multiple HRFs simply involves running multiple independent analyses. fMRIstat analyses should be conducted as usual for each HRF. For example, in our analysis we run fMRIstat for peak3, peak5, peak7, and peak9.

Combining HRFs

While it is possible to simply look at the ‘best’ tmap from each subject, defining which tmap is best can be problematic. For example, is a higher t value or a larger number of voxels more important? Also, these separate analyses do not address the issue of within subject variability. Combining the separate t-maps into a single map partially addresses this problem.

The combined method is quite simple. A new image file is created using the t-maps for the individual models as the template. Each voxel is independently examined, and the value of the tmap with the highest absolute value is written into that voxel. We use 5 contiguous voxels with a t value of over 3.1 as our threshold, which corresponds to p < 0.05 corrected for the multiple analysis in a scan of 64x64x25 voxels of 5mm x 5mm x 5mm.

Our investigations have shown this method to be quite useful. In our epilepsy studies, we do not always see activation even if we record a large number of events. This is probably an issue of the signal to noise ratio of the BOLD response to spikes. However, the use of the multi HRF method increased the percentage of patients with responses from 45% to 62.5% (Bagshaw et al 2004). In addition, several areas are better modeled by the multi HRFs than the standard Glover HRF. While the Glover HRF is generally quite good at detecting activations, the multi HRFs analysis was often better for deactivations, which tended to peak later (but were not an undershoot of an earlier response).

Examples of tmaps for activation. In this case, peak3 and peak5 seem to be the best tmaps. The combined map has the best elements of all the four peak maps.

 

Note the activation on the right side of this image is strongest peak7, while the left sided activation is strongest in peak5. This is a good example of an activation that is optimized by the use of more than one model HRF.

Examples of tmaps for activation. In this case, peak3 and peak5 seem to be the best tmaps. The combined map has the best elements of all the four peak maps.

  Note the activation on the right side of this image is strongest peak7, while the left sided activation is strongest in peak5. This is a good example of an activation that is optimized by the use of more than one model HRF.
 

Deactivation maps that show a widespread decrease in deactivation that is not present in peak3 and modeled poorly by peak5. This is an excellent example of how a single HRF model (such as peak5) can miss large widespread deactivations.

Deactivation maps that show a widespread decrease in deactivation that is not present in peak3 and modeled poorly by peak5. This is an excellent example of how a single HRF model (such as peak5)can miss large widespread deactivations.

Back to top