# 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

**DOI**: 10.5812/ans.82585

### 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, CO_{2} 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

### 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:

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

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 mm^{3}.

#### 3.5. FMRI Task Design

For fMR imaging in block design, CO_{2} gas was administered using a face-mask under physiological monitoring. The FMRI data were registered by administration of CO_{2} 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 CO_{2} 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.

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.

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.

^{a}

Standard | TD | DD | GAM | FIR | IL | |
---|---|---|---|---|---|---|

Time to peak | 7.5 | 7.5 | 7 | 5.5 | 7 | 8.3 |

Height | 1 | 1.008 | 1.178 | 0.7912 | 0.9992 | 0.8912 |

Width | 5.5 | 5.5 | 5 | 5.5 | 6 | 5.8 |

MSE | - | 0.052 | 0.105 | 0.146 | .055 | 0.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.

^{a}The 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.

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.

IL Model | FIR Model | TD Model | |
---|---|---|---|

Height | 0.076577 | 0.071911 | 0.121912 |

Time to peak | 4.184738 | 12.79348 | 4.195652 |

Width | 3.934783 | 2.413043 | 3.043478 |

MSE | 0.143363 | 0.23641 | 0.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.

IL Model | FIR Model | TD Model | |
---|---|---|---|

Height | 0.089164 | 0.099497 | 0.139391 |

Time to peak | 4.21001 | 14.23241 | 4.265332 |

Width | 3.334783 | 3.251252 | 3.122598 |

MSE | 0.091099 | 0.186892 | 0.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.

IL Model | FIR Model | TD Model | |
---|---|---|---|

Height | 0.093541 | 0.075221 | 0.116532 |

Time to peak | 4.205297 | 12.442018 | 4.152324 |

Width | 3.652218 | 3.102328 | 3.25138 |

MSE | 0.121870 | 0.220289 | 0.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.

### Acknowledgements

### Footnotes

### 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].

## LEAVE A COMMENT HERE: