Effects of Hemodynamic Response Function Selection on Rat fMRI Statistical Analyses.

The selection of the appropriate hemodynamic response function (HRF) for signal modeling in functional magnetic resonance imaging (fMRI) is important. Although the use of the boxcar-shaped hemodynamic response function (BHRF) and canonical hemodynamic response (CHRF) has gained increasing popularity in rodent fMRI studies, whether the selected HRF affects the results of rodent fMRI has not been fully elucidated. Here we investigated the signal change and t-statistic sensitivities of BHRF, CHRF, and impulse response function (IRF). The effect of HRF selection on different tasks was analyzed by using data collected from two groups of rats receiving either 3 mA whisker pad or 3 mA forepaw electrical stimulations (n = 10 for each group). Under whisker pad stimulation with large blood-oxygen-level dependent (BOLD) signal change (4.31 ± 0.42%), BHRF significantly underestimated signal changes (P < 0.001) and t-statistics (P < 0.001) compared with CHRF or IRF. CHRF and IRF did not provide significantly different t-statistics (P > 0.05). Under forepaw stimulation with small BOLD signal change (1.71 ± 0.34%), different HRFs provided insignificantly different t-statistics (P > 0.05). Therefore, the selected HRF can influence data analysis in rodent fMRI experiments with large BOLD responses but not in those with small BOLD responses.

Introduction

Functional magnetic resonance imaging (fMRI) was originally introduced in 1990. Since then, it has been modified to enable investigations on different functional aspects of the brain. The most popular fMRI technique is blood-oxygen-level dependent (BOLD) contrast, which relies on local deoxyhemoglobin changes (Ogawa et al., 1990, 1992). Owing to its advantages of absent radiation burden and non-invasiveness, BOLD fMRI has become a pivotal method for understanding brain function and physiological conditions (Tsurugizawa et al., 2010; Rana et al., 2013; Wu et al., 2014; Nasrallah et al., 2015; Golestani et al., 2016). The applications of BOLD fMRI in animals such as rats, have recently received increased attention. The majority of rodent fMRI studies have been conducted by using electric stimulation to induce somatosensory stimulation and to estimate activations in the primary sensory cortex (Silva et al., 1999; Keilholz et al., 2004; Tenney et al., 2004; Weber et al., 2006; Masamoto et al., 2007; Pelled et al., 2007; Shih et al., 2009, 2011, 2013;Ramos-Cabrer et al., 2010; Just et al., 2013; Sanganahalli et al., 2013; Nasrallah et al., 2015). Previous works have emphasized the importance of rodent fMRI studies in elucidating crucial topics in neuroscience research. These topics include functional recovery (Pelled et al., 2007; Ramos-Cabrer et al., 2010), pain processing (Shih et al., 2009, 2011), and neurodegenerative diseases (Tenney et al., 2004; Sanganahalli et al., 2013).

The construction of the hemodynamic response function (HRF) of the signal response to an external stimulus is the essential step in the statistical analysis of fMRI data for identifying activation regions. The typical HRF used in rat fMRI studies is the boxcar-shape HRF (BHRF) (Keilholz et al., 2004; Pelled et al., 2007; Shih et al., 2009, 2011; Yang et al., 2013; Nasrallah et al., 2015), which is based on the standard on/off format of the external stimulation. Meanwhile, the use of the so-called canonical HRF (CHRF), a sophisticated HRF based the convolution of a BHRF with the sum of two gamma functions, has also been suggested by other groups (Kim et al., 2007; Yu et al., 2012). The advantage of BHRF is that BHRF can be used for the rapid assessment of brain activation to further refine the CHRF. Although HRF selection in human fMRI studies have been widely discussed (Aguirre et al., 1998; Handwerker et al., 2004; Shan et al., 2014), the body of research that compares the HRF selection in rodent fMRI studies remains insufficient (Chavarrias et al., 2016).

Therefore, the central objective of this study is to systematically investigate whether fMRI activation detection can be affected by the selected HRF. To achieve this objective, three HRFs (Lu et al., 2005; Lindquist et al., 2009; Nasrallah et al., 2015) were employed to model the BOLD signal. The extent of brain activation by electric stimulation is task-dependent, with whisker pad stimulation projecting larger somatosensory regions than forepaw stimulation (Yu et al., 2010). Thus, to map different sensory processing in the brain cortex, we subjected two groups of rats to whisker pad and forepaw electric stimulations. The estimated BOLD signal changes and t-statistics among three HRFs were compared. Such a comparison may provide recommendations for future rat fMRI studies.

Materials and Methods

Functional Magnetic Resonance Imaging Experiments

A total of 20 male Sprague-Dawley rats (280–345 g) were used in this study. Laboratory animals were housed in plastic cages with soft bedding and were maintained on a 12-h light/dark cycle. Food and water were available ad libitum, and the room was temperature controlled. This study was carried out in accordance with the recommendations of National Institutes of Health guide for the care and use of Laboratory animals. The protocol was approved by the China Medical University.

All MRI experiments were conducted on a 7T animal MRI scanner (Bruker ClinScan 70/30, Germany) with a gradient strength of 630 mT/m. A volume coil and a surface coil were used for signal excitation and reception, respectively. All rats were initially anesthetized with 4% isoflurane (ISO), and then was reduced to 1–1.2% ISO during fMRI (Liu et al., 2004). Each rat was secured in a head holder with ear bars and a bite bar to prohibit head motion. The rats were placed on a heated water pad to maintain body temperature at ∼37°C while in the magnet.

Rats were subsequently divided into two groups. In the first group (n = 10), needle electrodes were inserted under the skin of the left whisker pad for mapping the primary somatosensory cortex barrel field (S1BF). In the second group (n = 10), needle electrodes were inserted under the skin of the left forepaw for mapping the primary somatosensory cortex forelimb region (S1FL). Electric stimulation was performed by a stimulator (Isolated Pulse Stimulator Model 2100, Washington, DC, United States) supplying 3 mA, 330 μs pulses repeated at 3 Hz to either the left whisker pad or the forepaw upon demand. The stimulation paradigm of the fMRI experiment consisted of a block design. The stimulation paradigm including an initial 75 s period of resting followed by five cycles alternating 15 s of electric stimulation with 75 s of resting was implemented, with a total duration of 525 s. BOLD MR images were simultaneously acquired during this period. The BOLD imaging parameters were field of view (FOV) = 30 mm × 30 mm, matrix size = 64 × 64, 7 coronal slices, thickness = 1 mm, no gap, repetition time (TR)/echo time (TE) = 1000 ms/25 ms, and single-shot gradient echo echo-planar imaging. Anatomical images were obtained by turbo-spin-echo with scanning parameters of TR = 2560 ms, TE = 38 ms, echo train length = 7, number of excitation = 1, FOV = 30 mm × 30 mm, matrix size = 320 × 320, and slice thickness = 1 mm.

Data Analysis

The data analysis for each animal was performed using first-level analyses in SPM8. The voxel-by-voxel statistical analysis of fMRI data was based on the general linear model (GLM) analysis. The dependent variable was the BOLD signal, and the first regressor was the HRF. Three different types of HRFs were employed in this study to test the influence of HRF selection on fMRI sensitivity. We first employed a block design stimulus function that consisted of alternating blocks of resting and active conditions. This block was designated as the BHRF. The second HRF was designed by convolving a BHRF with the sum of two gamma functions (from SPM8) and was designated as the CHRF (Friston, 2003; Lindquist et al., 2009). The third one was the impulse response function (IRF) that was fitted to a gamma-variate function [IRF(t) = ktbe-t/c] appropriate for cerebral blood volume (CBV) weighted fMRI signal under rat whisker stimulation, with k = 0.9, b = 0.64, and c = 4.42 (Lu et al., 2005). The second regressor was the intercept, a vector of ones. A high-pass filter of 1/128 Hz was used to detrend fMRI data (Tanabe et al., 2002).

Three different t-maps and magnitude estimate β maps were generated from the corresponding HRF. For different GLM models, the fractional signal change of each voxel was calculated using the same equation as follows:

where Sbase and Sact are the signals at baseline and activation, respectively. β1 and β2 are the estimated parameters for the two regressors from the GLM (Miao et al., 2014). The resulting t- value map and signal-change map were used for the following analysis. Voxels with t-values higher than the threshold of 4.8 (corresponding to P = 10-6) and only groups of at least four activated pixels (Keilholz et al., 2006; Weber et al., 2006) were regarded as significantly activated.

Voxels with t-values greater than 4.8 from three HRF analyses were considered as activated voxels and used for regions of interest analysis. Averaged t-values and signal changes within S1BF and S1FL were calculated by averaging the t-values and signal changes of its constituent activated voxels as defined above, respectively. Differences in t-values and estimated BOLD signal changes in the primary somatosensory cortex (S1) among the three HRFs were tested through one-way analysis of variance (ANOVA) tests with repeated measures. If the effect was observed in the ANOVA test results, post hoc Tukey’s honest significant difference test was employed. Task-based BOLD data were analyzed using in-house Matlab (The MathWorks, Natick, MA, United States) scripts. Data were expressed as mean ± standard error.

Results

Robust fMRI activations in the contralateral side of the brain were detected in all rats under whisker pad or forepaw stimulation. fMRI signal time curves from the S1 of two representative rats under whisker pad stimulation or forepaw stimulation are shown in Figure 1A,B, respectively. The fMRI signals gradually increased and then gradually decreased to the baseline.

FIGURE 1

Functional magnetic resonance imaging (fMRI) signal time curves from the primary somatosensory cortex of two representative rats under (A) whisker pad stimulation or (B) forepaw stimulation. The shaded regions indicate the 15-s stimulation period.

Figure 2 shows the group-level activation maps of 10 rats under electric stimulation. The maps shown in this figure were obtained through GLM with CHRF. Consistent with previous reports, robust activations were detected in the S1BF, secondary somatosensory cortex (S2), and primary somatosensory cortex upper lip region (S1ULp) of rats under whisker pad stimulation (Yu et al., 2010). Significant activations in the S1BF associated with forepaw stimulation were found. In line with previous studies (Silva et al., 1999; Shih et al., 2013, 2014), activation in the S2 was not easily detectable in animals under forepaw stimulation. Notably, the volume of the active region in rats under forepaw stimulation was smaller than that in rats under whisker pad stimulation as indicated by the decreased cortical somatosensory representation in the rat brain. The presented activations under whisker pad and forepaw stimulation were similar when employing GLM with the BHRF or IRF was employed.

FIGURE 2

Group-level result of brain regions activated by mysticial pad electrical stimulation and forepaw electrical stimulation. Analyses were performed using a one-sample Student’s t-test with a cluster size of four voxels.

Blood-oxygen-level dependent signal changes and t-values quantified through GLM analysis are plotted in Figure 3A,B, respectively. The estimated BOLD percentage changes in the S1BF were 4.02 ± 0.38, 4.31 ± 0.42, and 3.22 ± 0.34% when the GLM used CHRF, IRF, and BHRF, respectively. One-way ANOVA with repeated measures showed that HRF had a significant effect on BOLD signal change estimation (P < 0.001), where the estimated BOLD signal change in S1BF for GLM with BHRF was significantly lower than those with CHRF (P < 0.01) and IRF (P < 0.01).

FIGURE 3

Comparisons of (A) estimated BOLD signal changes and (B) t-statistics between two hemodynamic response functions (HRFs) (∗∗P < 0.01, ∗P < 0.05).

The selection of HRFs also had significant effects on the estimated BOLD signal changes in response to forepaw stimulation (P < 0.01). The estimated BOLD signal change in the S1FL for GLM with IRF (1.71 ± 0.34%) was significantly higher than that with BHRF (1.44 ± 0.31%, P < 0.05). The comparisons between CHRF and IRF and between CHRF and BHRF were not significantly different (both P > 0.05).

The analysis of the influence of HRF selections yielded a similar pattern for the quantification of t-values. For whisker pad stimulation, the statistical power of the t-values derived from GLM with CHRF or IRF significantly improved relative to those of the t-values derived from BHRF (Figure 3B, P < 0.001), suggesting that the use of CHRF or IRF improved the statistical significance of voxels. The comparison between CHRF and IRF showed insignificant differences (P > 0.05). No significant differences in the quantified t-values were detected for the data obtained under forepaw stimulation (P = 0.13) when the GLM used any HRF.

The voxel-wise comparisons of t-values between GLM with CHRF and BHRF are displayed in Figure 4. Compared with BHRF, GLM with CHRF significantly improved the activation maps mainly in S1BF under whisker pad stimulation. This result indicated that the CHRF is more appropriate for stimulations with large BOLD signal change. By contrast, BHRF did not increase sensitivity relative to CHRF. No difference was detected under forepaw stimulation in either direction. The comparisons between GLM with IRF and BHRF showed similar patterns.

FIGURE 4

Brain regions showing significant differences in t-value quantifications between two HRFs. Relative to that with BHR, GLM with CHRF significantly improved the t-statistics for S1BF under whisker pad stimulation. No difference was detected in the opposite direction or under forepaw stimulation in either direction. Analyses were performed using paired Student’s t-test with a cluster size of four voxels.

Discussion

In this study, we found that the choice of the HRF is crucial in the computation of activations in rat fMRI studies, especially in studies involving stimulations with large BOLD signal changes, such as whisker pad stimulation. We tested three types of HRFs: CHRF, IRF, and BHRF. BHRF is simple and popular among rat fMRI studies (Keilholz et al., 2004; Pelled et al., 2007; Nasrallah et al., 2015). However, its statistical estimations tend to vary with the type of responses. Experimental fMRI data have been used to demonstrate that BHRF affects the results of statistical analyses by underestimating BOLD magnitude changes and t-values.

Our results demonstrate that BHRF significantly underestimates signal changes and t-statistics. Given that the ground truth is unknown, CHRF or IRF may also be overestimating signal changes and t-statistics. In addition, the effect of HRF selection may be dependent on the brain area. To address this issue, we simulated the fMRI signal time curves (Shan et al., 2014) with ground truth HRFs (Supplementary Figure 1) to test the performance of the different HRF models. Numerical simulations showed that for the fMRI signal time curve simulated from the BHRF with BOLD signal changes of more than 1.7%, GLM with CHRF produced t-values that were significantly lower than those produced by GLM with BHRF (Supplementary Figure 2A, P < 0.05). The same scenario occurred when the fMRI signal time curve was simulated with CHRF but using BHRF in the GLM (Supplementary Figure 2B, P < 0.05). The potential explanation to this phenomenon is that when the BOLD signal change is small and the signal-to-noise ratio is low, everything is buried with noise and no difference could be detected among the HRFs, thereby reducing the relative advantage of CHRF. When the BOLD signal change is larger and accompanied with a clear peak instead of a plateau, CHRF is preferred over BHRF due to its existing peak. These simulation results may provide further evidence for the possible underestimation of rodent fMRI results by BHRF, particularly when the BOLD signal change is large and has a clear peak. Moreover, the effect of HRF selection may be independent of the brain area.

Early work done by Chavarrias et al. (2016) showed that a simple approach using a boxcar response provides better model fitting results than complex approaches. This conclusion is not congruent with our present findings. We found that CHRF and IRF improve the statistical power, especially for stimulations with large BOLD signal changes. This deviation may be attributed to the low temporal resolution of 3 s employed in the aforementioned study compared with the 1 s that we used in functional volume acquisition. When temporal resolution is low, the intrinsic hemodynamic response could be blurred and could also cloud the true response, yielding biased estimations (Kim et al., 1997). Additional comprehensive experiments with different temporal resolution settings may help to further parse out this issue.

Our present results provide some interesting insights into the HRF selection in rat fMRI studies. Our data suggest that when the BOLD signal change is large, such as that under whisker pad stimulation, CHRF and IRF are appropriate candidates for modeling the BOLD response, even though IRF is derived from CBV-fMRI, which may provide a response different from that provided by BOLD. The CHRF is derived from the sum of two gamma functions, whereas the IRF is derived from a single gamma function. The results of these HRFs are similar and comparable. This finding is in agreement with that of a previous human study showing that two gamma functions are neither better nor worse than a single gamma function (Handwerker et al., 2004). In CHRF, the second gamma function is included to model the post-stimulus undershoot. We carefully inspected our data and found that the post-stimulus undershoot is not observed in the data for every rat. Thus, inter-subject variation may restrict the statistical power of the approach and may imply that as long as a peak exists instead of a plateau in HRF, the BOLD response can be correctly modeled. It should also be noted that the parameters of the two gamma functions used in this study were empirically derived from SPM. The application of the default setting from SPM to rodent studies may not be appropriate since the parameters in SPM were originally designed for human studies. However, determining the hemodynamic parameters for animal fMRI studies is non-trivial (Silva et al., 2007) and may not be generally performed across studies. In addition, other factors such as the anesthesia regime and the different targeted brain regions may contribute to the variations in hemodynamic parameters. Therefore, the use of the default setting of the two gamma functions from SPM is simple and convenient, thus gaining increasing popularity in rodent fMRI studies (Just et al., 2013; Niranjan et al., 2016). Although the parameters from CHRF were not optimized in this study, the advantage of CHRF over BHRF in improving the statistical power was demonstrated in this work. Additional research similar to the current one but with different hemodynamic parameters can be an important area for future work.

In the field of rodent fMRI studies, Student’s t-test (Tenney et al., 2004; Weber et al., 2006; Masamoto et al., 2007; Sanganahalli et al., 2013; Poplawsky and Kim, 2014) and cross-correlation (CC) analysis (Keilholz et al., 2004; Pelled et al., 2007; Shih et al., 2009, 2011; Yang et al., 2013; Nasrallah et al., 2015) are popular statistical strategies for localizing brain regions activated by a task. The principle of Student’s t-test is to compare the data between “baseline (off)” and “stimulus (on)” phases, thus providing high t-scores for large differences with small standard deviations, and low t-scores for small differences with large standard deviations. Notably, the comparison between “on” and “off” corresponds to BHRF and may imply that similar to BHRF, the Student’s t-test may underestimate the t-value when the BOLD response is large. CC analysis takes a HRF of expected neural responses and correlates it with the MRI signal variations of each voxel. Correlation coefficients are calculated and converted to t-values (Hinkle et al., 2003) to generate the activation map. In this regard, the CHRF or IRF can be considered as complementary HRF when employing CC analysis to assess functional activities.

The results in the present work should be interpreted in the context of several limitations. First, in the present study, each rat was subjected to either whisker pad or forepaw electrical stimulation. Therefore, we were unable to make within-rat comparisons. The HRF effect on whisker pad or forepaw stimulations may be affected by physiological differences across rats. The duration of ISO anesthesia is known to influence the functional connectivity in rats (Magnuson et al., 2014). Nevertheless, the time-dependent effects of ISO on electric stimulation fMRI studies remain unclear. As a result, we used separate groups to maintain the same duration of anesthesia in our experiments. The different periods of anesthetization should not be a major concern in the experimental design. Second, electrical stimulation parameters are often dependent on the type of anesthetic (Huttunen et al., 2008; Liu et al., 2013) and sensory system used (Melzer et al., 2006; Just et al., 2010). In this study, the stimulation parameters were 3 Hz and 330 μs electrical pulses. These parameters were first identified to induce robust BOLD response to rat forepaw somatosensory stimulus under alpha-chloralose (Brinker et al., 1999; Silva et al., 2007). Thus, optimal stimulus parameters must be employed to clarify the effects of HRF selection on the analysis of fMRI data obtained through whisker pad stimulation under ISO anesthesia.

Conclusion

We demonstrated that rat fMRI results could be influenced by HRF selection, especially for stimulations with large BOLD response. BHRF is a simple and straightforward HRF but may underestimate the magnitude of BOLD response and the t-values of statistical tests. Sophisticated HRFs, such as CHRF and IRF, provide robust estimation. Our results suggest that CHRF and IRF could serve as complementary HRFs in the analysis of rat fMRI data.

Author Contributions

S-LP conceptualized the study. S-LP and C-MC designed the experiments and prepared the manuscript. S-CC, C-YH, and W-CS performed the experiments and conducted data analysis. C-TS, C-WH, and S-CC participated in the data collection.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.