# Hemodynamic Response Function Modeling to Determine the Areas with High Blood Supply in Block-Design fMRI Experiments

AUTHORS

Seyedeh Mahboobe Seyed Abbasi 1 , Mohammad Ali Oghabian 2 , Seyed Salman Zakariaee 3 , * , Abbas Rahimiforoushani 1 , **

1 Department of Epidemiology and Biostatistics, Faculty of Public Health, Tehran University of Medical Sciences, Tehran, Iran

2 Department of Medical Physics and Biomedical Engineering, Faculty of Medicine, Tehran University of Medical Sciences, Tehran, Iran

3 Department of Medical Physics, Faculty of Paramedical Sciences, Ilam University of Medical Sciences, Ilam, Iran

Corresponding Authors:

How to Cite: Seyed Abbasi S M , Oghabian M A, Zakariaee S S, Rahimiforoushani A . Hemodynamic Response Function Modeling to Determine the Areas with High Blood Supply in Block-Design fMRI Experiments, Arch Neurosci. Online ahead of Print ; 6(Brain Mapping):e82585. doi: 10.5812/ans.82585.

ARTICLE INFORMATION

Archives of Neuroscience: 6 (Brain Mapping); e82585
Published Online: March 18, 2019
Article Type: Research Article
Received: July 24, 2018
Revised: December 2, 2018
Accepted: February 27, 2019

# Crossmark

#### CHEKING

##### Abstract

Background: Determining the activated brain areas due to different activities is one of the most common targets in functional magnetic resonance imaging (fMRI) data analysis, which could be carried out by hemodynamic response function (HRF) evaluation. The HR functions reflect changes of cerebral blood flow (CBF) in response to neural activity.

Objectives: In this study, five models of HRF estimation were evaluated based on a simulated dataset. Models with higher accuracy were used to determine HRF parameters of the block-design fMRI data.

Methods: The fMRI data were acquired in a 3 Tesla scanner. For block-design fMR imaging, CO2 gas was administered using a face-mask under physiological monitoring. Three patients with brain tumors were scanned. The fMRI data analysis was performed using the SPM 12 MATLAB toolbox. Akaike’s information criterion (AIC), Schwarz’ Bayesian (SBC), and mean square error (MSE) criteria were used to select the best HRF estimation model.

Results: In simulation studies, the estimated HRFs by the canonical HRF plus its temporal derivative (TD), finite impulse response (FIR), and inverse logistic (IL) models were almost equal to the standard HRF. Mean square error, AIC, and SBC indices were ignorable for TD, FIR, and IL models (MSE/AIC/SBC magnitudes for TD, FIR, and IL models were 0.052/-1235.1/-1223.9, 0.055/-1206.4/-1194.9, and 0.068/-1091.5/-1049.2, respectively), which indicates that these models could accurately estimate HRF in block design fMRI studies.

Conclusions: The HRF models could non-invasively evaluate the change of MR signal intensity under cerebrovascular reactivity (CVR) conditions and they might be helpful to investigate changes in human cerebral blood flow.

## Keywords

Functional Magnetic Resonance Imaging Hemodynamic Response Function Inverse Logistic Model Finite Impulse Response Model Canonical HRF Plus Its Temporal Derivative Model

Copyright © 2019, Archives of Neuroscience. This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial 4.0 International License (http://creativecommons.org/licenses/by-nc/4.0/) which permits copy and redistribute the material just in noncommercial usages, provided the original work is properly cited.

### 1. Background

In functional magnetic resonance imaging (fMRI) studies, the physical parameters of the brain change in response to a given task. Regional increase in blood flow rate is the parameter, which could be activated with a specific stimulation. Cerebrovascular reactivity (CVR) is the change of cerebral blood flow in response to a vasoactive stimulus (1, 2)

Determining activated brain areas due to different activities is one of the most common targets in fMRI data analysis, which could be carried out by hemodynamic response function (HRF) evaluation (1, 3-6).

Neurovascular connections can be non-invasively studied based on the HRF information. Thus, HRF modeling plays an important role in fMRI data analysis and evaluation of the cerebrovascular diseases (4, 7, 8). The HR functions reflect changes of cerebral blood flow (CBF) in response to neural activity (1, 3, 4). This function is different between subjects and has a significant variation between brain regions. Therefore, correct HRF modeling could enhance the accuracy of fMRI studies (1). Intensity, onset latency, and duration of the underlying brain metabolic activity were evaluated based on HRF shape characteristics (including height, delay, and duration) (5). Hemodynamic response function parameters, including time to peak (T) and height (H) provide information about brain area activation, and full-width at half-max (W) of a stimulated HRF determines the time of activity (5).

The fMRI task design has a considerable effect on the efficiency and detection power of the study. Previous studies indicated that event-related designs have high estimation efficiency with poor detection power. However, good detection power at the cost of low estimation efficiency is achieved using block designs (9). Block designs are mostly used in BOLD fMRI studies. Therefore, selection of the best model to estimate HR functions in block design studies could have a critical role in patient data analysis.

### 2. Objectives

In this study, five models including gamma (GAM) (10), inverse logistic (IL) (11), the canonical HRF plus its temporal and dispersion derivatives (DD), canonical HRF plus its temporal derivative (TD) (12, 13), and finite impulse response (FIR) models (14, 15) were evaluated to estimate HR function of a simulated dataset. Models with higher accuracy were used to determine HRF parameters of patient data.

### 3. Methods

#### 3.1. Hemodynamic Response Function Models

The relationship between blood oxygenation level dependent (BOLD) response and stimulus is presented using the Equation 1:

Equation 1.$Y=(s×h)(t)$

Where y (t) is the signal at time t and is the convolution between the hemodynamic response (h(t)) and stimulus (s(t)) functions. In these studies, a fixed canonical shape is assumed for h(t) (10).

The HRF was estimated based on five models, including gamma (GAM), inverse logistic (IL), canonical HRF plus its temporal and DD, canonical HRF plus its TD, and finite impulse response (FIR) models. The Lindquist approach was used to determine the HRF parameters (5).

#### 3.2. HRF Model Selection

Akaike’s information criterion (AIC), Schwarz’ Bayesian (SBC), and mean square error (MSE) criteria were used to select the best HRF estimation model. The AIC, SBC, and MSE indices are described in Equations 2 - 4

Equation 2.
Equation 3.
Equation 4.$MSE=∑i=1nei2/(n-2)$

Where n and p are sample size and number of parameters, respectively. The model with the least amounts of AIC, SBC, and MSE had a better performance for HRF estimations.

#### 3.3. Data Simulation

The simulated BOLD signal was produced based on real data. Changes have been regularly quantified in terms of commencement and period of neural activity. The stimulus function was convoluted with SPM’s canonical HRF to produce the simulated signal. It was assumed that repetition time (TR) and duration of stimulation were 3 and 120 seconds, respectively. The onset stimuli were presented in 60 and 240 seconds. The dataset consisted of the BOLD signal plus white noise. Furthermore, a random interpersonal variance with a standard deviation equal to one-third of interpersonal variance was added to each time series. On the simulation study, five models were fitted on the simulated BOLD signal. The simulation studies were repeated one hundred times. The HRF features, including T, H, and W indices, were extracted for each fitting model.

#### 3.4. FMRI Data Acquisition

The fMRI data were acquired in a 3 Tesla Siemens (Munich, Germany) MAGNETOM TIM Trio scanner using a 12-channel Matrix head coil (Siemens Medical Systems, Erlangen, Germany). Functional MR images with BOLD contrast were acquired using an echo-planar imaging (EPI) sequence. The applied scanning parameters were: TR = 3 seconds, TE = 30 msec, flip angle = 90°, matrix size = 64 × 64, and voxel size = 3 × 3 × 3 mm3.

#### 3.5. FMRI Task Design

For fMR imaging in block design, CO2 gas was administered using a face-mask under physiological monitoring. The FMRI data were registered by administration of CO2 for 120 seconds with 60-second rest intervals. The HRF in the cerebral areas with higher blood supply was determined based on changes in blood CO2 level.

#### 3.6. Data Processing

The FMRI data were processed using the MATLAB software (version 8.6, The MathWorks TM, Natick, Massachusetts, United States). The simulation procedures were repeated 100 times and the best HRF estimation models were selected. Then, these models were implemented on the patient data. The fMRI data analysis was performed using the SPM 12 MATLAB toolbox (Wellcome Trust Center for Neuroimaging).

#### 3.7. Patient Data

Three patients with brain tumors (two males and one female; mean age, 28.33 years; age range, 21 to 34 years) were imaged. The study was approved by the local committee for medical research ethics. Informed consent was obtained from the patients prior to the study.

### 4. Results

#### 4.1. Data Simulation

The simulated signal is depicted in Figure 1A. For time series estimation, the efficacies of the models were evaluated. The standard and estimated time series by five models are shown in Figure 1B. As it could be seen from the Figure 1, the estimated signals by IL, FIR, and TD models are in acceptable consistency with the standard time series.

Figure 1. A, the simulated signal in block-design method, TR = 3 seconds, onset = 60 and 240 seconds, duration = 120 seconds, block length = 420 seconds. B, The estimated signals by five mentioned models. The standard time series is also depicted. Abbreviations: DD, the canonical HRF plus its temporal and dispersion derivatives, FIR, finite impulse response; GAM, gamma; IL, inverse logistic; TD, the canonical HRF plus its temporal derivative models.

For determining the best HRF estimation model, an HRF was produced as the standard function and it was estimated by five HRF estimation models.

The standard HRF is illustrated in Figure 2A. Figure 2B shows the estimated HRFs by five models. In Figure 2, the standard HRF is also depicted. The estimated HRFs by IL, FIR, and TD models were almost equal to the standard HRF.

Figure 2. A, illustration of the standard HRF. B, The standard and estimated HRFs by five models. Abbreviations: DD, the canonical HRF plus its temporal and dispersion derivatives, FIR, finite impulse response; GAM, gamma; IL, inverse logistic; TD, the canonical HRF plus its temporal derivative models.

A scatter plot was used to evaluate the model efficiency for signal estimation. For the best model with acceptable efficiency to estimate the simulated signal, spots would be more closely located around the bisector of the first quadrant. The results of IL, FIR, and TD models are almost placed around the bisector of the first quadrant. The simulation procedures were repeated 100 times and the average values of H, T, and W parameters were registered. The characteristics of the estimated HRFs (including H, T, and W) are listed in Table 1. For each HRF estimation model, MSE, AIC, and SBC indices were also listed in the Table 1. The standard values for H, T, and W indices were 1, 7.5, and 5.5, respectively. From the Table 1 and the visual evaluation, it could be concluded that IL, FIR, and TD models present better HRF estimation than the other models.

Table 1. The Characteristics of the Estimated HRFs by Five Modelsa
StandardTDDDGAMFIRIL
Time to peak7.57.575.578.3
Height11.0081.1780.79120.99920.8912
Width5.55.555.565.8
MSE-0.0520.1050.146.0550.068
AIC--1235.1-883.5-175.8-1206.4-1091.5
SBC--1223.9-871.5-165.7-1194.9-1049.2

Abbreviations: AIC, Akaike’s information criterion; DD, the canonical HRF plus its temporal and dispersion derivatives, FIR, finite impulse response; GAM, gamma; IL, inverse logistic; MSE, mean square error; SBC, Schwarz’Bayesian criteria; TD, the canonical HRF plus its temporal derivative models.

aThe characteristics of the standard signal are also listed. For each HRF estimation model, MSE, AIC, and SBC indices are listed.

#### 4.2. Patient Data Analysis

The efficiency magnitudes of five HRF estimation models were evaluated based on a simulated time series. The best HRF estimation models were selected to analyze the patient data. Three patients with brain tumors were imaged in block-design fMRI paradigms. The selected models were used to estimate HR functions of the patient data.

Figure 3 shows cerebral regions with higher blood supply in three patients with brain tumors. The subjects’ data were analyzed with P value = 0.001, voxel threshold = 20, onset stimulus = 60 and 240, and duration of stimulation = 120 seconds. For the patients, the areas with higher blood supply were placed at the frontal, temporal, and occipital lobes, respectively. Figure 4 shows the estimated signals using models at the cerebral area with higher blood supply. The estimated HRFs by IL, FIR, and TD models are illustrated in Figure 5.

Figure 3. The standard views of the brain (transverse, sagittal, and coronal views) for three patients with brain tumors. For the patients, the areas with higher blood supply are placed at the frontal (A), temporal (B), and occipital (C) lobes, respectively.
Figure 4. The block-design fMRI signals of the patients and the results of the fitted IL, FIR, and TD models on the signal. For each patient, the block-design fMRI signal and the related fitted models were depicted in separate sub figures. Abbreviations: FIR, finite impulse response; IL, inverse logistic; TD, the canonical HRF plus its temporal derivative models.
Figure 5. The HRFs estimated using the IL, FIR, and TD models. For each patient, the HR functions estimated by the models were depicted in separate sub figures. Abbreviations: FIR, finite impulse response; IL, inverse logistic; TD, the canonical HRF plus its temporal derivative models.

The average values of H, T, and W indices at the regions of interest (ROIs) are listed in Tables 2 - 4, respectively. The MSE, AIC, and SBC indices are also listed in Tables 2 - 4 for each model.

Table 2. The characteristics of HRF for a 32-year-old woman with brain tumor
IL ModelFIR ModelTD Model
Height0.0765770.0719110.121912
Time to peak4.18473812.793484.195652
Width3.9347832.4130433.043478
MSE0.1433630.236410.135698
AIC-266.9885-196.4844-273.88888
SBC-258.1636-187.6595-265.0638

Abbreviations: AIC, Akaike’s information criterion; FIR, finite impulse response; IL, inverse logistic; MSE, mean square error; SBC, Schwarz’Bayesian criteria; TD, the canonical HRF plus its temporal derivative models.

Table 3. The characteristics of HRF for a 34-year-old man with brain tumor
IL ModelFIR ModelTD Model
Height0.0891640.0994970.139391
Time to peak4.2100114.232414.265332
Width3.3347833.2512523.122598
MSE0.0910990.1868920.088393
AIC-330.416-229.8149-334.6386
SBC-325.8137-220.9900-321.5914

Abbreviations: AIC, Akaike’s information criterion; FIR, finite impulse response; IL, inverse logistic; MSE, mean square error; SBC, Schwarz’Bayesian criteria; TD, the canonical HRF plus its temporal derivative models.

Table 4. The characteristics of HRF for a 21-year-old man with brain tumor
IL ModelFIR ModelTD Model
Height0.0935410.0752210.116532
Time to peak4.20529712.4420184.152324
Width3.6522183.1023283.25138
MSE0.1218700.2202890.108532
AIC-289.6759-206.7976-305.9025
SBC-280.8510-197.9726-297.0776

Abbreviations: AIC, Akaike’s information criterion; FIR, finite impulse response; IL, inverse logistic; MSE, mean square error; SBC, Schwarz’Bayesian criteria; TD, the canonical HRF plus its temporal derivative models.

### 5. Discussion

The mentioned models were fitted on the simulated signal in block-design studies. Figure 2A shows that all of the models, except gamma, can fairly estimate simulated BOLD signal. The standard and estimated HRFs by each model are shown in Figure 2B. Figure 2 visually shows that TD, FIR, and IL models estimate the simulated HRF. The estimated HRF by the DD model has a higher height than the simulated HRF and the estimated HRF using gamma model has a different time to peak than the simulated HRF. Also, the estimated HRF using the IL model has a different time to peak and height yet these differences could be ignored. For simulation studies, which were repeated 100 times, average magnitudes of T, W, and H indices were registered. For the canonical HRF plus its temporal derivative model, the characteristics of the estimated HRF were closer to the standard HRF. The MSE, AIC, and SBC indices were ignorable for TD, FIR, and IL models, which indicates that these models could accurately estimate HRF in a block-design fMRI study. From the model selection criteria and the characteristics of the estimated HRFs, it could be concluded that TD is the best model for HRF estimation.

In Lindquist et al.’s studies, seven models, including canonical HRF, canonical HRF plus TD, canonical plus TD and DD, FIR, sFIR, two-parameter gamma distribution, and total three IL models were evaluated to estimate HRF. The patients wore a heat arm band scaled in event-related method (warm, mild painful, fairly painful and pain threshold) and then they were imaged. The results showed that the IL model estimates HRF better than TD and sFIR models. The TD model yields the weakest estimation of HRF (5, 11).

In this study three selected models were evaluated on real data. Figure 5 shows the signals of cerebral areas with higher blood supply and the results of the fitting models. As it could be seen from the Figure 5, the models have been acceptably fitted on the given signal, and TD model has the best fitting trend line. A square region of interest (3 × 3 pixels in size) was extracted to evaluate the cerebral areas with higher blood supply and their averages of T, W, and H magnitudes. The MSE, AIC, and SBC indices were smaller for the TD model. Unlike Lindquist et al.’s study, the TD model has a better HRF estimation for block-design fMRI data. In Lindquist et al.’s study, the IL model provides a better HRF estimation for event-related data (5). In Shan et al.’s study HR function (in block-design method) was modeled using five models, including the canonical HRF used in SPM8, canonical HRF, the sum of two and three gamma functions, and the sum of three inverse logit functions used by Lindquist et al. It was found that in the event-related method, four IL models had a higher efficiency for HRF estimation. However, the results of real data (in block-design method) showed that two-gamma with six parameters had the best HRF estimation (16).

#### 5.1. Conclusions

In the present study, HRF modeling of the cerebral area with higher blood supply was investigated. For a simulated dataset, HRF modeling was evaluated using five models. Based on the results, TD, FIR, and IL models were selected for HRF modeling in a block-design method. These models were used to analyze the block-design fMRI data for patients with brain tumors. The results showed that the TD model yields better HRF estimation to evaluate the changes in human cerebral blood flow. These models could be helpful to investigate the change of signal intensity under CVR conditions, non-invasively.

### References

• 1.

Aguirre GK, Zarahn E, D'Esposito M. The variability of human, BOLD hemodynamic responses. Neuroimage. 1998;8(4):360-9. doi: 10.1006/nimg.1998.0369. [PubMed: 9811554].

• 2.

Logothetis NK, Pauls J, Augath M, Trinath T, Oeltermann A. Neurophysiological investigation of the basis of the fMRI signal. Nature. 2001;412(6843):150-7. doi: 10.1038/35084005. [PubMed: 11449264].

• 3.

Bellgowan PS, Saad ZS, Bandettini PA. Understanding neural system dynamics through task modulation and measurement of functional MRI amplitude, latency, and width. Proc Natl Acad Sci U S A. 2003;100(3):1415-9. doi: 10.1073/pnas.0337747100. [PubMed: 12552093]. [PubMed Central: PMC298787].

• 4.

Handwerker DA, Gonzalez-Castillo J, D'Esposito M, Bandettini PA. The continuing challenge of understanding and modeling hemodynamic variation in fMRI. Neuroimage. 2012;62(2):1017-23. doi: 10.1016/j.neuroimage.2012.02.015. [PubMed: 22366081]. [PubMed Central: PMC4180210].

• 5.

Lindquist MA, Meng Loh J, Atlas LY, Wager TD. Modeling the hemodynamic response function in fMRI: efficiency, bias and mis-modeling. Neuroimage. 2009;45(1 Suppl):S187-98. doi: 10.1016/j.neuroimage.2008.10.065. [PubMed: 19084070]. [PubMed Central: PMC3318970].

• 6.

Muthukumaraswamy SD, Edden RA, Jones DK, Swettenham JB, Singh KD. Resting GABA concentration predicts peak gamma frequency and fMRI amplitude in response to visual stimulation in humans. Proc Natl Acad Sci U S A. 2009;106(20):8356-61. doi: 10.1073/pnas.0900728106. [PubMed: 19416820]. [PubMed Central: PMC2688873].

• 7.

D'Esposito M, Deouell LY, Gazzaley A. Alterations in the BOLD fMRI signal with ageing and disease: A challenge for neuroimaging. Nat Rev Neurosci. 2003;4(11):863-72. doi: 10.1038/nrn1246. [PubMed: 14595398].

• 8.

Iadecola C. Neurovascular regulation in the normal brain and in Alzheimer's disease. Nat Rev Neurosci. 2004;5(5):347-60. doi: 10.1038/nrn1387. [PubMed: 15100718].

• 9.

Maus B, van Breukelen GJ, Goebel R, Berger MP. Optimal design for nonlinear estimation of the hemodynamic response function. Hum Brain Mapp. 2012;33(6):1253-67. doi: 10.1002/hbm.21289. [PubMed: 21567658].

• 10.

Worsley KJ, Friston KJ. Analysis of fMRI time-series revisited--again. Neuroimage. 1995;2(3):173-81. doi: 10.1006/nimg.1995.1023. [PubMed: 9343600].

• 11.

Lindquist MA, Waugh C, Wager TD. Modeling state-related fMRI activity using change-point theory. Neuroimage. 2007;35(3):1125-41. doi: 10.1016/j.neuroimage.2007.01.004. [PubMed: 17360198].

• 12.

Friston KJ. Imaging neuroscience: Principles or maps? Proc Natl Acad Sci U S A. 1998;95(3):796-802. doi: 10.1073/pnas.95.3.796. [PubMed: 9448243]. [PubMed Central: PMC33800].

• 13.

Friston KJ, Glaser DE, Henson RN, Kiebel S, Phillips C, Ashburner J. Classical and Bayesian inference in neuroimaging: Applications. Neuroimage. 2002;16(2):484-512. doi: 10.1006/nimg.2002.1091. [PubMed: 12030833].

• 14.

Glover GH. Deconvolution of impulse response in event-related BOLD fMRI. Neuroimage. 1999;9(4):416-29. doi: 10.1006/nimg.1998.0419. [PubMed: 10191170].

• 15.

Goutte C, Nielsen FA, Hansen LK. Modeling the haemodynamic response in fMRI using smooth FIR filters. IEEE Trans Med Imaging. 2000;19(12):1188-201. doi: 10.1109/42.897811. [PubMed: 11212367].

• 16.

Shan ZY, Wright MJ, Thompson PM, McMahon KL, Blokland GG, de Zubicaray GI, et al. Modeling of the hemodynamic responses in block design fMRI studies. J Cereb Blood Flow Metab. 2014;34(2):316-24. doi: 10.1038/jcbfm.2013.200. [PubMed: 24252847]. [PubMed Central: PMC3915209].