A cortical rat hemodynamic response function for improved detection of BOLD activation under common experimental conditions

NeuroImage 208 (2020) 116446

Contents lists available at ScienceDirect

NeuroImage journal homepage: www.elsevier.com/locate/neuroimage

A cortical rat hemodynamic response function for improved detection of BOLD activation under common experimental conditions Henriette Lambers, Martin Segeroth, Franziska Albers, Lydia Wachsmuth, Timo Mauritz van Alst, Cornelius Faber * Translational Research Imaging Center (TRIC), Department of Clinical Radiology, University Hospital Münster, Albert-Schweitzer-Campus 1, Münster, D-48149, Germany

A R T I C L E I N F O

A B S T R A C T

Keywords: BOLD fMRI Data analysis Rat animal model Hemodynamic response function Optogenetic stimulation Mechanical stimulation

For a reliable estimation of neuronal activation based on BOLD fMRI measurements an accurate model of the hemodynamic response is essential. Since a large part of basic neuroscience research is based on small animal data, it is necessary to characterize a hemodynamic response function (HRF) which is optimized for small animals. Therefore, we have determined and investigated the HRFs of rats obtained under a variety of experimental conditions in the primary somatosensory cortex. Measurements were performed on animals of different sex and strain, under different anesthetics, with and without ventilation and using different stimulation modalities. All modalities of stimulation used in this study induced neuronal activity in the primary somatosensory cortex or in subcortical regions. Since the HRFs of the BOLD responses in the primary somatosensory cortex showed a close concordance for the different conditions, we were able to determine a cortical rat HRF. This HRF is based on 143 BOLD measurements of 76 rats and can be used for statistical parametric mapping. It showed substantially faster progression than the human HRF, with a maximum after 2.8
0.8 s, and a following undershoot after 6.1
3.7 s. If the rat HRF was used statistical analysis of rat data showed a significantly improved detection performance in the somatosensory cortex in comparison to the commonly used HRF based on measurements in humans.

1. Introduction The calculation of brain activation maps from fMRI measurements requires a realistic model of the BOLD response. Since preclinical fMRI studies on small animals often use analysis tools that have been developed for the analysis of studies in humans, it is crucial to verify that the small animal BOLD response has been characterized correctly. Previous studies showed that the BOLD response of rats have substantially faster temporal dynamics than the response of humans (Silva et al., 2007; Zwart et al., 2005). Both studies identified a higher flow velocity in veins and capillaries of rats as possible cause for the different BOLD responses. Another factor raised by Silva et al. is that humans have larger blood vessels that cover longer distances, resulting in a higher contribution of large vessel effects to the signal. These results indicate that it is mandatory to use an analysis tool that has been optimized for rats for the analysis of rat BOLD data. Many of the commonly used analysis tools are based on the assumption that the BOLD response linearly follows the stimulation which is used to induce neuronal activation. Accordingly, the BOLD response corresponds to the convolution of neuronal activity with

the impulse response function of the vascular system. This function is named hemodynamic response function (HRF). The HRF is often characterized by two gamma functions and then named canonical HRF. Peng et al. (2019) recently showed that a canonical HRF is a good choice for the analysis of rat fMRI data. There are many methods for statistical analysis of fMRI data. A commonly used method is the general linear model (GLM) (Friston et al., 1994, 1995; Worsley and Friston, 1995). Briefly, time courses from each voxel are compared with a model of the expected signal during neuronal activity using a linear regression. There are several basis sets comprising canonical, gamma, Fourier, finite impulse response (FIR) (Henson and Friston, 2007), B-spline (Genovese, 2000), which can be used for GLM based analysis. Each set consists of different basis functions, which are used to form regressors for the analysis. The number of regressors is defined by the model order. The commonly used canonical basis set (Henson and Friston, 2007) consists of the HRF and its temporal and dispersion derivatives. Convolution of the HRF and derivatives with the stimulation paradigm provides the regressors of the basis set. However, the canonical basis set has been

* Corresponding author. Albert-Schweitzer-Campus 1, Building A16, 48149, Münster, Germany. E-mail address: [email protected] (C. Faber). https://doi.org/10.1016/j.neuroimage.2019.116446 Received 30 September 2019; Received in revised form 13 November 2019; Accepted 5 December 2019 Available online 14 December 2019 1053-8119/© 2019 Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

H. Lambers et al.

NeuroImage 208 (2020) 116446

2.1. BOLD measurements

considered as inflexible (Liu et al., 2017). In some cases a mismatch between HRF and the measured BOLD responses may result in both false negative (Handwerker et al., 2004) or false positive voxels (Lindquist et al., 2009) in the analysis. Liu et al. (2017) described the FIR basis set as one of the most flexible sets, and reported that the 7th to 9th order provide good balance between detection and characterization performance. The FIR set models the BOLD response with a number of contiguous boxcar functions. The number is given by the model order (Henson and Friston, 2007). However, the use of a flexible model may result in detection of physiologically implausible responses, leading to wrong interpretation of the data. One of the most commonly used software tools for analysis of fMRI data is Statistical Parametric Mapping (SPM). It provides a canonical HRF by default (Friston, 2017). This HRF is based on human datasets (human HRF), which were recorded in a study by Friston and colleagues (Friston et al., 1998; Henson and Friston, 2007). Due to the significant differences between the BOLD responses of humans and rats, application of the human HRF may not be appropriate for the analysis of small animal data. Furthermore, several studies indicate that experimental conditions like anesthesia (Masamoto et al., 2009; Schlegel et al., 2015; Schroeter et al., 2014; van Alst et al., 2019), as well as stimulation type (Albers et al., 2018) and stimulation duration (Huettel and McCarthy, 2000; Pfeuffer et al., 2003; Vazquez and Noll, 1998; Yes¸ilyurt et al., 2008) influence the BOLD response. Therefore, it is not clear whether a universal rat HRF is appropriate to analyze data obtained under different experimental conditions. The aim of this study was to improve the framework of GLM for statistical analysis of rodent data obtained under commonly used experimental conditions. For this purpose, we determined the HRFs of BOLD measurements which were partly obtained in three previous fMRI studies (Albers et al., 2018; Amirmohseni et al., 2016; van Alst et al., 2019). We derived a cortical rat HRF based on 143 BOLD measurements of the rat primary somatosensory cortex. Additionally, we developed a functional t-test and compared the time courses of HRFs across different rat strains and anesthetic regimens. Finally, the detection performance of the GLM framework after implementation of the cortical rat HRF was compared to the performance when using the human HRF or the FIR basis set.

2.1.1. Animal preparation The procedure of animal preparation for optogenetic stimulation has been described in detail previously (Albers et al., 2018; Schmid et al., 2017). fMRI measurements were performed on animals expressing the opsins C1V1TT in excitatory neurons (CamKIIa promoter) or Chronos pan-neuronally (Synapsin1 promoter) in the forelimb region of the sensory cortex (S1FL, 18 animals). Animal preparation for all fMRI measurements was performed under surgical depth of isoflurane anesthesia. For optogenetic stimulation a 200 μm optical fiber was implanted approximately 100 μm above the opsin expressing region and attached to the skull with UV glue. Bipolar needle electrodes were placed subcutaneously in one fore- or hindpaw for electrical paw stimulation. For mechanical paw stimulation, a von Frey filament (100 g) was positioned above one hindpaw. 39 animals were ventilated. These animals were intubated at the beginning of the preparation. The tubus was connected to a pneumatic valve, which was driven by a ventilator (MRI-1 Ventilator, CWE Inc., Ardmore, PA, USA). The endexpiratory CO2 was continuously measured with a CO2 analyzer (CapStar-100 CO2 Analyzer, CWE Inc., Ardmore, PA, USA). 2.1.2. Optogenetic, visual, electrical and mechanical stimulation Four different block paradigms with three different durations of stimulation were used: 10 s stimulation and 20 s rest, 5 s stimulation and 25 s rest and 1 s stimulation with 29 s or 14 s rest. For optogenetic stimulation, light was coupled into the implanted fiber according to the chosen paradigm. Pulses of 1 up to 11 ms duration with frequencies of 9 Hz, 12 Hz, 120 Hz and 150 Hz and an intensity of 64–162 mW/mm2 at the fiber tip were used. A blue laser at 488 nm wavelength delivered the light for the activation of Chronos, while C1V1TT was activated with a green laser at 552 nm wavelength (both OBIS, Coherent, Dieburg, Germany). Visual stimulation was also performed using the blue Laser (488 nm) with a pulse duration of 2 up to 11 ms and a stimulation frequency of 9 Hz. Electrical paw stimulation was performed with 1–2 ms pulse duration at frequencies of 5 Hz, 9 Hz or 12 Hz. A pulse strength of either 1–1.5 mA or 5 mA was used for innocuous or noxious electrical stimulation, respectively. Pulses were generated by a constant current stimulator (DS5, Digitimer, Welwyn Garden City, UK). For mechanical stimulation, a frequency of 1 Hz and a pulse duration of 0.5 s was used. The filament above the paw was pneumatically controlled with air pressure by a custom-built device (Amirmohseni et al., 2016). Pressures from 2 bar up to 4 bar were used. Since the pressure did not linearly reflect the force applied to the paw, the assignment whether a stimulation was noxious or innocuous was based on the occurrence of pain-related activations in respective BOLD maps as described by Amirmohseni et al. (2016). These maps were created in ImageJ using a t-test based custom-written script (Schneider et al., 2012) as described before (Schmid et al., 2017).

2. Materials and methods Animal experiments were carried out according to the German Tierschutzgesetz and were approved by the Landesamt für Natur, Umwelt und Verbraucherschutz of Nordrhein-Westfalen, Germany (approval IDs 87-51.04.2010.A274, 87-51.04.2011.A065, 84-02.04.2014.A347 and 84-02.04.2015.A427). We analyzed 281 datasets from experiments in 117 rats with a body weight between 155 g and 465 g. Firstly, the BOLD responses of 251 fMRI measurements from 23 different groups were extracted and examined. Groups were defined by the experimental conditions, comprising rat strain (Fischer or SpragueDawley (SD)), sex, anesthesia (medetomidine or medetomidine with 0.7 % isoflurane), ventilation (ventilated or not ventilated), pain (innocuous or noxious), stimulation duration (1 s, 5 s, or 10 s stimulation), frequency (5 Hz, 9 Hz, 12 Hz, 120 Hz, 150 Hz) and, stimulation mode (optogenetic stimulation of cortex, visual stimulation, electrical or mechanical paw stimulation), activated region (cortex or subcortical regions) and acquisition mode (2D echo planar imaging (EPI) or 1D line scanning (LS)) (Fig. 1). Suppl. Table 1 lists the studies from which the datasets of the different groups were taken and summarizes their experimental conditions. Secondly, statistical analysis was performed on 30 datasets for testing the detection performance of the GLM framework after implementation of different HRFs.

2.1.3. fMRI measurements All MRI experiments were performed on a 9.4 T Bruker Biospec 94/20 small animal scanner with a 0.7 T/m gradient system (Bruker BioSpin GmbH, Ettlingen, Germany), using a 10-mm or 20-mm surface coil or a four-element phased-array head coil. Firstly, an anatomical image was acquired using a 2D fast spin echo sequence (RARE), TR/TEEff: 2000/50 ms, RARE factor 8, Matrix: 256x256, FOV: 28  26 mm2 or 30  30 mm2, slice thickness: 1.2 mm, 8–14 contiguous slices. Subsequently, a localized shim of the area of interest was performed using Mapshim (Bruker). 2D BOLD measurements were carried out for 10 min using a single-shot gradient EPI sequence (TR/TE: 1000/18 ms, Matrix: 80x80, same geometry as anatomy, flip angle: 60 or 65 , 20 repetitions of the 30 s paradigm). In 2

H. Lambers et al.

NeuroImage 208 (2020) 116446

Fig. 1. BOLD measurements from 23 different groups were examined. Groups were defined by different experimental conditions: MRI sequence, strain, sex, anesthesia, ventilation, stimulation duration and frequency, pain, mode of stimulation and by the activated region. The groups and their experimental conditions are listed in Suppl. Table 1.

and consequently SNR was too low. The visual pathway included the superior colliculus, dorsal lateral geniculate and lateral posterior nucleus (Pawela et al., 2008). Within these ROIs a voxelwise Mann-Whitney U test (U test) determined whether the mean signal during stimulation differed significantly (p < 0.05) from the mean signal during the rest period (Fig. 2B). The custom written MATLAB script is available as a supplement (Suppl. A). To take into account the delayed onset of the BOLD response (due to the slow vascular response), for statistical testing both start and end of stimulation were shifted by 2 s. p-values were Bonferroni corrected by the number of voxels in the ROI. If the mean signal during stimulation was larger, the signal was considered positive. For calculation of the BOLD time courses, the signal was summed over all voxels showing a significant and positive BOLD response, and subsequently averaged over all stimulation cycles. A baseline correction was performed by setting the mean signal of the pre-stimulation period (5 s prior to stimulation) to zero. Only measurements with more than four voxels that showed a significant and positive signal change were used in this study. The approach to use individually identified voxels for each measurement might introduce bias due to different vascular structures in individual animals. We therefore repeated the analysis using fixed ROIs. In agreement to human studies by Boynton et al. (1996), no differences were observed between time courses obtained with the different analysis methods (Suppl. Fig. 1).

addition to the 2D measurements, 1D BOLD measurements were performed on 9 animals using a linescanning (LS) sequence (Albers et al., 2018) (TR/TE: 50/18 ms, Matrix: 1x128, FOV: 2,1  10 mm2, slice thickness: 1.2 mm, 1 slice, flip angle: 13 , duration: 32 min). LS employs a FLASH sequence without phase encoding resulting in 1D lines instead of images. This allows measurements with a significantly improved temporal resolution. For LS measurements, the FOV was positioned over the center of the S1Fl region and the direction of frequency encoding was perpendicular to the surface of the cortex. LS measurements were performed during electrical paw or optogenetic stimulation. Optogenetic LS measurements were only performed on animals that expressed C1V1TT. fMRI measurements were performed under medetomidine sedation, which was induced by a subcutaneous bolus injection of 0.04 mg/kg and followed by a continuous infusion of 0.05 mg/kg/h. In the next 10 min, isoflurane was stepwise reduced to zero. The fMRI measurements were started at least 30 min after the medetomidine bolus injection. On 7 ventilated Fischer rats electrical stimulation was also performed with a combination of isoflurane (0.7 % in 1 l/min 80/20 air/O2) and medetomidine (same dose as above). During measurements, animals were placed on a temperaturecontrolled bed and fixated with bite and ear bars. Rectal temperature was held at 36.4
1.3  C. Respiration rate of non-ventilated animals was held at 67
22 bpm. For ventilated animals ventilation settings were adjusted to hold the end-tidal CO2 at 2.8
0.4 %, thereby avoiding hypercapnia, which is known to reduce the BOLD response (Cohen et al., 2002; Jones et al., 2005).

2.2.2. Extraction of BOLD responses from LS measurements LS data were analyzed as previously described using a customized MATLAB script (Albers et al., 2018). Briefly, data of the complete first stimulation cycle and the first 20 acquisitions from the following cycles were discarded, because in these acquisitions magnetization had not reached the steady-state. After performing a discrete Fourier transform along the direction of frequency encoding, data points were assigned to the cortex according to the distance (150–1300 μm) from the cortical surface. BOLD time courses were determined by summing up the signal of these data points and subsequent averaging over all stimulation cycles.

2.2. Determination of the rat HRFs 2.2.1. Extraction of BOLD responses from EPI measurements fMRI datasets from the EPI measurements were realigned using SPM 8 (Statistical Parametric Mapping, Functional Imaging Laboratory, Wellcome Trust Centre for Neuroimaging, London, UK). The first five scans of each measurement were discarded to avoid pre-steady-state artefacts. Afterwards, the positive BOLD response was extracted using a customwritten MATLAB script (Release 2014b; The MathWorks, Inc., Natick, MA, USA). With exception of visual stimulation all modalities of stimulation induced neuronal activity in the primary sensory cortex (S1): Sensory stimulation in the regions of fore or hind limb (S1Fl or S1Hl), optogenetic stimulation of the cortex in S1Fl. Therefore, the entire S1 region on the activated side of the brain was marked as region of interest (ROI). For visual stimulation only subcortical regions belonging to the visual pathway were analyzed (marked as ROI in Fig. 2A), because the position of the surface coil was not optimized to cover the visual cortex,

2.2.3. Fit of the extracted BOLD response According to the linear model it was assumed that the BOLD response can be predicted by the convolution of neuronal activity and HRF. Furthermore, it was assumed that the neuronal activity has a boxcar shape following the stimulation paradigm. Therefore, the convolution of the stimulation paradigm and an HRF was chosen as predictor for the BOLD response. This predicted response was fitted to the time course of the positive BOLD response (Fig. 2C). Shan et al. (2014) have shown that 3

H. Lambers et al.

NeuroImage 208 (2020) 116446

Fig. 2. Methods for extracting and fitting the BOLD response and for determining the HRFs. (A) The entire S1 region on the activated side of the brain or subcortical regions belonging to the visual pathway were marked as ROI. (B) Within this ROI the signals from rest and stimulation periods were compared in a voxelwise manner using a U test. If signals differed significantly and showed a positive signal change, the voxel was taken as activated. The time courses of the BOLD responses were calculated by averaging across all activated voxels and all stimulation cycles. (C) The convolution of stimulation paradigm and canonical HRF was fitted to the time course of the positive BOLD response. This delivered the parameters (A; b; p1 ; p2 ; V) of the canonical HRF of the BOLD response. (D) The HRFs of the measurements were normalized and averaged per group. (E) The parameters of the HRFs were determined by fitting the canonical HRF to their time courses. (F) It was evaluated if there was a statistically significant difference between pairs of HRFs by using a functional t-test: First, the t-value T of each time point τ was determined and the maximum t-value was considered. Subsequently, for each permutation the maximum t-value was determined. The p-value is given by the probability for the maximum t-value of the real distribution being smaller than the maximum t-value of the permutations.

4

H. Lambers et al.

NeuroImage 208 (2020) 116446

2.4. Statistical assessment of the different HRFs

an HRF consisting of two gamma functions is the best choice for modeling of the BOLD response in block design studies. Therefore, we used an HRF consisting of two superimposed gamma functions. Such an HRF is named canonical HRF. One gamma function models the maximum, the second the undershoot. Here, the canonical HRF was defined as in SPM (Friston, 2017), expanded by a parameter which allows variable amplitudes:
A ⋅ ebt ⋅

bp1 bp2 ⋅ tp1 1  ⋅ tp2 1 Γðp1 Þ V ⋅ Γðp2 Þ

Statistical testing of the HRFs for the 23 individual experimental conditions was performed using two different approaches. 2.4.1. Functional t-test for the comparison of HRFs Time courses of the normalized HRFs for the different experimental conditions were compared pairwise. To test for a statistically significant difference between the HRFs of two groups, consisting of a number of n1 and n2 HRFs, a functional t-test as described by Ramsay et al. (2009) was used. Fig. 2F illustrates the method of the functional t-test: Firstly, the t-value T of each time point τ was determined.

 (1)

In analogy to the HRF implemented in SPM, it was assumed that both gamma functions of the canonical HRF have the same dispersion parameter, which describes the width of the gamma functions. Furthermore, the HRF was characterized by the position parameters and (representing width and peak positions of the gamma functions) and the ratio parameter (the ratio of the amplitudes of peak and undershoot). The parameters (A; b; p1 ; p2 ; V) were adjusted to fit the predicted response to the BOLD response. Using a custom written MATLAB script, the fit was optimized with the method of least squares to minimize the deviation between fit and measured BOLD response. This optimization procedure was repeated for 40,000 starting points to find the global minimum of the fitting of EPI measurements. Starting points were generated within given limits of the parameters using the function multistart. Chosen limits are given in Suppl. Table 2. Due to the high temporal resolution of LS measurements (50 ms) the optimization was repeated for 2000 starting points for these data. To avoid unphysiological results, the fit had to fulfill two linear inequalities: Firstly, the time tp , at which the HRF was at maximum, had to be smaller than the time tu , at which the undershoot reached its minimum (tp  tu ). Secondly, the difference between tu and tp had to be less than 6 s (tu  tp  6 s). The latter inequality was chosen, because the undershoot of the HRFs did not always occur directly after the peak, if the difference was greater than 6 s. For measurements with 1 s or 5 s stimulation, only the first 20 s of the extracted time courses were fitted, because the BOLD response had already subsided within this time frame. For EPI measurements the fit of the time course and the time course itself have been normalized to the maximum of the time course. The mean squared error (MSE) during the peak of the BOLD response was determined for these normalized time courses. If the MSE was smaller than 0.035 and no parameter reached its limits, the fit was considered successful. For LS measurements no MSE was determined. These fits were considered successful if no parameter reached its limits. In cases where the amplitude A reached the upper limit given by the fitting procedure, time courses of the BOLD responses were normalized and fitted again. This normalization was deemed acceptable since the amplitude of the fit was not relevant for the following evaluation.

TðτÞ ¼ jHRFs1 ðτÞ  HRFs2 ðτÞj qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 var½HRFs1 ðτÞ þ n12 var½HRFs2 ðτÞ n1

(2)

Subsequently, the maximum t-value max (T (τ)) was considered. To decide if this t-value represented a significant difference, the HRFs were permutated randomly as often as possible and for each permutation the maximum t-value was determined. Subsequently, the probability was calculated that the maximum t-value of the original distribution was smaller than the maximum t-values of the permutations. This probability was expressed as p-value of the functional t-test. A total of 33 functional ttests between the experimental groups were made. For this purpose, a MATLAB script was written which is available as supplement (Suppl. B). The resulting p-values were Bonferroni corrected. 2.4.2. Statistical analysis using condition-specific and non-specific HRFs Here, we specifically tested whether analysis of data obtained at experimental conditions that had been identified to show significantly different HRFs requires the application of a condition-specific HRF in order to achieve optimal detection performance. For this purpose, the condition-specific HRFs of the different groups were implemented in SPM 12 (Friston, 2017). Subsequently, the detection performance of SPM using these different HRFs was compared. Comparisons were made for analysis of 30 EPI datasets using SPM 12. These datasets had not previously been used for determination of HRFs, and were acquired in experiments with either 5 s (n ¼ 10) or 10 s (n ¼ 10) innocuous electrical paw stimulation of non-ventilated female Fischer rats, or from 5 s stimulation of ventilated animals (n ¼ 10). Prior to analysis, the datasets were realigned and smoothed using a 0.5 mm Gaussian kernel using SPM. Analysis was performed by using the 1st order canonical basis set with the convolution of the HRF and the stimulation paradigm as regressor. BOLD maps were calculated using a t-test (p < 0.05, family wise error correction). Each dataset was first analyzed using its condition-specific HRF. The resulting maximum t-values and cluster sizes, which indicated the number of activated voxels, were considered as gold standard of reference. Subsequently, each dataset was analyzed using HRFs, specific for different experimental conditions. The resulting cluster sizes and maximum t-values were normalized to the gold standard values. The obtained relative cluster sizes and maximum tvalues allowed for direct comparison between detection performance of the condition-specific HRF and non-specific HRFs. A less than 25 % difference to the gold standard was arbitrarily defined as equally well performing.

2.3. Generation and characterization of condition-specific HRFs HRFs were normalized and averaged across all measurements of single groups (Fig. 2D). The resulting time courses were used to derive the parameters of the canonical HRFs which were specified for the experimental conditions of the examined groups (condition-specific HRFs). In SPM the canonical HRF is represented by equation (1) without the scaling factor A. Therefore, the parameters (b; p1 ; p2 ; V) were determined as outlined in Fig. 2E. Firstly, the HRFs were normalized to the maximum of the human HRF. Secondly, the canonical HRF (equation (1) without A) was fitted to the HRF time courses using MATLAB in analogy to the procedure described in 2.2.3. The fit was optimized using the method of least squared errors, the optimization problem was repeated 40,000 times and the parameter limits were set as b ¼ [0.1, 5.0], p1 ¼ [1, 50], V ¼ [1, 15] and p2 ¼ [5, 90]. For an implementation in SPM 12 (Friston, 2017) the resulting parameters had to be inserted in a subprogram which can be found in the supplement (Suppl. C).

2.5. Generation and characterization of the cortical rat HRF HRFs of all measurements of all groups, except for those groups that had been found to differ significantly in the functional t-test analysis were normalized and averaged. For implementation of the cortical rat HRF in SPM the parameters (b; p1 ; p2 ; V) of the corresponding canonical HRF have been characterized using the method described in section 2.3.

5

H. Lambers et al.

NeuroImage 208 (2020) 116446

Parameters for the rat HRFs were obtained by fitting the canonical HRF convolved with the stimulation paradigm to the measured response for each experimental condition. Fig. 3 shows four examples of the extracted and fitted responses. Most of the time courses showed the typical shape of a BOLD response. The amplitude of the extracted time courses was 2.0
0.8 %. The signal began to increase 1–2 s after the start of stimulation. The response to short stimulations decreased immediately after reaching its maximum. In contrast, BOLD responses to 10 s stimulation remained on a plateau after reaching the maximum and decreased after the end of the stimulation. In total, 78 % of the extracted responses were fitted successfully.

2.6. Detection performance of the rat HRF Detection performance of analysis using the cortical rat HRF was assessed using the 10 datasets of 5 s electrical paw stimulation of nonventilated rats and procedures to calculate relative cluster sizes and tvalues as described in 2.4.2. Firstly, detection performance of the 1st order canonical basis set using the rat HRF was compared with the performance when using the human HRF. Secondly, the detection performance of the 3rd order canonical basis set using the rat HRF was compared with the performance when using the human HRF. For this purpose, the HRF and its temporal and dispersion derivatives were used as basis functions and were convolved with the stimulation paradigm, and subsequently used as regressors. Finally, the datasets were also analyzed with the 7th order FIR basis set, which consisted of a set of seven boxcar functions. The width of each box equaled one seventh of 14 s, which represented the length of the BOLD response at 5 s stimulation. BOLD maps of analysis with the 3rd order canonical basis set or the FIR set were calculated using an F-test (p < 0.05, family wise error correction). The relative cluster sizes and F-values were calculated in analogy to the procedure for the first order.

3.2. Comparison of the determined HRFs To assess whether experimental conditions have an impact on the HRF of rats, all HRFs were normalized and experimental conditions were compared pairwise using the functional t-test. The threshold (p ¼ 0.050) for significant comparisons was corrected for multiple comparisons using Bonferroni correction, resulting in a threshold of p ¼ 0.003. As shown in Table 1 almost no significant differences were found between different stimulation modalities (electrical paw, mechanical paw, optogenetic stimulation), between male and female rats, between Fisher and SD rats, between ventilated and non-ventilated animals, between anesthesia with isoflurane versus medetomidine, between innocuous and noxious stimulation, for stimulation of 5 s versus 10 s, and finally for different stimulation frequencies (5, 9, 12 Hz, electrical paw stimulation; 9, 12, 120, 150 Hz, optogenetic stimulation). Only three out of 33 comparisons showed a p-value smaller than 0.003: Significant differences in the temporal dynamics of peak and undershoot were observed between the HRFs of subcortical regions and the S1Fl (visual versus electrical paw stimulation, p < 0.001, Fig. 4A), and between the HRFs of electrical paw versus optogenetic stimulation using the opsin Chronos (p ¼ 0.002, Fig. 4E). The HRFs of 1 s versus 5 s stimulation duration showed significant differences in the time to peak (p ¼ 0.001, Fig. 4C). Differences in the time to peak, however non-significant, were also seen between the HRFs of electrical paw versus optogenetic stimulation using the opsin C1V1 (p ¼ 0.018, Fig. 4K). Finally, non-significant differences in the temporal dynamics of the undershoot were seen between the HRFs of 5 s versus 10 s stimulation duration (p ¼ 0.011, Fig. 4G), and between the

2.7. Data and code availability statement All data and codes used in this study are available upon direct request. Additionally, MATLAB codes for the U test to find activated voxels, the functional t-test and HRF implementation are available as supplements to this publication: Suppl. A contains the codes for the U test (voxelwise_u_test). Suppl. B contains the codes for the functional t-test (f_t_test1, f_t_test2 and f_t_test_T_max). Suppl. C contains the code for the implementation of the rat HRF in SPM (spm_my_defaults). Data and code sharing comply with the requirements of the funding body of this study and with the institutional ethics approval. 3. Results 3.1. Fit of the extracted BOLD response To determine the HRFs of rats under different experimental conditions, the BOLD responses of 251 fMRI measurements were extracted.

Fig. 3. Representative BOLD time courses for different stimulation lengths and sequences. Extracted (black) and fitted (green) BOLD time courses for EPI measurements of group 1 (1 s stimulation, A), group 16 (5 s stimulation, B), group 20 (10 s stimulation, C) and for one linescanning measurement of group 22 (D). 6

H. Lambers et al.

NeuroImage 208 (2020) 116446

Table 1 Pairwise comparisons of the calculated HRFs. The table summarizes all group comparisons, conditions which have been compared and the p-values, resulting from the functional t-test. p-values below 0.003 remained significant after Bonferroni correction, and are indicated by an asterisk. Comparisons for which relative cluster sizes and t-values have been calculated are printed in bold and shown in Fig. 4. Groups are indicated by their group number. n gives the number of measurements per group. difference

group a

group b

compared conditions

p value

anesthesia

3 (n ¼ 8)

4 (n ¼ 8)

medetomidine

medetomidine þ 0.7 % isoflurane

0.185

7 (n ¼ 5) 6 (n¼10) 1 (n ¼ 11) 1 (n¼11)

3 (n ¼ 8) 2 (n¼11) 7 (n ¼ 5) 3 (n¼8)

10 s

5s

0.239

10 s

5s

0.011

1s

10 s

0.033

1s

5s

* 0.001

14 (n ¼ 7) 13 (n ¼ 7) 12 (n ¼ 7) 13 (n ¼ 7) 12 (n ¼ 7) 12 (n ¼ 7) 7 (n ¼ 5) 8 (n ¼ 11) 7 (n ¼ 5)

15 (n ¼ 6) 15 (n ¼ 6) 15 (n ¼ 6) 14 (n ¼ 7) 14 (n ¼ 7) 13 (n ¼ 7) 8 (n ¼ 11) 5 (n ¼ 9) 5 (n ¼ 9)

120 Hz

150 Hz

0.811

12 Hz

150 Hz

0.894

9 Hz

150 Hz

0.492

12 Hz

120 Hz

0.763

9 Hz

120 Hz

0.517

9 Hz

12 Hz

0.182

9 Hz

12 Hz

0.713

12 Hz

5 Hz

0.077

9 Hz

5 Hz

0.115

sex

16 (n ¼ 13)

19 (n ¼ 9)

male

female

0.172

mode

2 (n ¼ 11) 3 (n¼8) 22 (n ¼ 8) 3 (n¼8) 11 (n ¼ 6) 17 (n ¼ 7) 18 (n ¼ 7)

10 (n ¼ 9) 11 (n¼6) 23 (n ¼ 6) 12 (n¼7) 12 (n ¼ 7) 20 (n ¼ 10) 21 (n ¼ 6)

el. paw stim.

optog. stim. (C1V1TT) optog. stim. (C1V1TT) optog. stim. (C1V1TT) optog. stim. (Chronos) optog. stim. (Chronos) mech. paw stim. mech. paw stim.

0.351

17 (n ¼ 7) 20 (n ¼ 10)

18 (n ¼ 7) 21 (n ¼ 6)

3 (n¼8)

duration

frequency

pain

region/ mode

sequence

el. paw stim. el. paw stim. el. paw stim. optog. stim. (C1V1TT) el. paw. stim el. paw. stim

Table 1 (continued ) difference

group a

group b

2 (n ¼ 11) 10 (n ¼ 9)

22 (n ¼ 8) 23 (n ¼ 6)

compared conditions

p value

EPI

LS

0.798

strain

2 (n ¼ 11)

19 (n ¼ 9)

Fischer

SD

0.088

ventilation

2 (n ¼ 11) 6 (n¼10) 10 (n ¼ 9)

3 (n ¼ 8) 7 (n¼5) 11 (n ¼ 6)

not ventilated

ventilated

0.246

not ventilated not ventilated

ventilated

0.012

ventilated

0.306

HRFs of ventilated versus non-ventilated animals (p ¼ 0.012, Fig. 4I). To assess the practical relevance of the differences between the individual HRFs, an independent analysis was performed. The six pairs of experimental groups described before, were subjected to an evaluation whether condition-specific HRFs were required for an optimal detection performance. These six pairs involved eight different experimental groups, for each of which a condition-specific HRF was implemented in the GLM framework of SPM 12 (Friston, 2017), as summarized in Suppl. Table 3. Datasets were first analyzed using its condition-specific HRF, and the resulting cluster size and maximum t-value were defined as gold standard for detection performance. Subsequently, each dataset was analyzed using a “wrong” condition-specific HRFs, defined by the mutual partner in the six pairs. Median values for the resulting cluster sizes and maximum t-values, normalized to the gold standard values, were below a threshold of 0.75 for the three pairs that were significantly different in the functional t-test analysis (Fig. 4B, D and F). For one of the three remaining pairs, the median values of the cluster sizes and maximum t-values, normalized to the gold standard values, were at this threshold of 0.75 (Fig. 4L). For the other two pairs the median values were over the threshold (Fig. 4H and J). Therefore, we conclude that the result of the statistical analysis agrees with the result of the functional t-test and that most experimental parameters assessed here do not significantly alter the HRF. Yet, two highly relevant exceptions were observed. The BOLD response in subcortical regions deviates from the response in cortex and cannot be modeled with the same HRF. Furthermore, the rat HRF changed with stimulation duration. Short stimuli (1 s) resulted in a significantly different HRF, as compared to longer stimulation durations.

0.018 0.014 * 0.002 0.226

3.3. A cortical rat HRF

0.663

To derive a cortical rat HRF all determined HRFs were normalized and averaged, except for those three conditions that had been found to differ significantly in the functional t-test analysis. We thus obtained a cortical rat HRF based on 143 BOLD measurements in the primary somatosensory cortex of 76 animals. Fig. 5A illustrates that this cortical rat HRF deviates substantially from the human HRF. The maximum already occurred after 2.8
0.8 s and the undershoot after 6.1
3.7 s, compared to 5 s and 15 s for the human HRF (Friston, 2017). Also the full width at half maximum (FWHM) of the rat HRF (2.3
1.0 s) was much smaller than the FWHM of the human HRF (5.3 s) (Friston, 2017). Errors were calculated based on the standard deviation over the fitted HRFs. The differences between the cortical rat HRF and the human HRF are significantly larger than the differences between the condition-specific HRFs (Fig. 5B). To determine the parameters of the canonical HRF which are needed for the use in SPM, the cortical rat HRF was fitted with the canonical HRF (Suppl. Fig. 2). The cortical rat HRF was characterized by the dispersion parameter b ¼ 2.0 the peak parameters p1 ¼ 7.4 and p2 ¼ 8.9 and the ratio parameter V ¼ 1.5.

0.431

innocuous

noxious

0.860

innocuous

noxious

0.344

9 (n¼14)

el. paw stim./ S1Fl

* <0.001

11 (n ¼ 6)

9 (n ¼ 14)

optog. stim. (C1V1TT)/S1Fl

12 (n ¼ 7)

9 (n ¼ 14)

optog. stim. (Chronos)/ S1Fl

visual stim./ subcortical regions visual stim./ subcortical regions visual stim./ subcortical regions

EPI

LS

0.043

0.040

0.717

7

H. Lambers et al.

NeuroImage 208 (2020) 116446

Fig. 4. Selected comparisons of six HRF pairs. Averaged HRFs which have been compared using the functional t-test are shown (A, C, E, G, I and K) for six pairs consisting of HRFs from eight different experimental groups (group number indicated by first number in the legend). Significant differences were found for HRFs from subcortical regions versus S1FL (A), for HRFs for 1 s versus 5 s stimulation (C), and for HRFs for sensory versus optogenetic stimulation using the opsin Chronos (E). Comparisons of 5 s versus 10 s stimulation (G), of ventilation versus no ventilation (I) and sensory versus optogenetic stimulation (C1V1TT, K) are also shown. Headings indicate the differences between the compared groups. The shaded area of the HRFs indicates their confidence interval. n gives the number of measurements per group. Boxplots (B, D, F, H, J, and L) show relative cluster sizes (cs) and t-values (t), which resulted from the analysis of 10 measurements each. For determination of the relative values, datasets were first analyzed with the condition-specific HRF for one of the compared groups. The resulting cluster size and maximum t-value were taken as gold standard of the detection performance. Secondly, the same datasets were analyzed with the condition-specific HRF of the other group. The resulting cluster sizes and maximum t-values were normalized to the gold standard values. For relative cluster sizes and t-values of at least 0.75 (red line) differences in detection performance were assumed negligible. In the boxplots, the central mark represents the median, boxes include the 25th and 75th percentiles, whiskers extend to the last data points within the 1.5*interquartile, and outliers were marked as dots.

either the human or the cortical rat HRF. Finally, datasets were analyzed using the 7th order FIR set. Analysis with the 1st order canonical set using the human HRF failed to detect activation for 60 % of the datasets. All other regressor sets were able to reliably detect activation. However, the clusters of the activation maps differed in size and the magnitude of the t- or F-values (Fig. 6). For analysis with the 1st and 3rd order of the canonical model using the human HRF and with the 7th order FIR, strong differences to the analysis with the optimized HRF were detected. Analysis using the 1st order

3.4. Assessment of detection performance of the cortical rat HRF To test the detection performance of the cortical rat HRF, statistical analysis of 10 datasets (all not ventilated and 5s electrical paw stimulation) was performed with different sets of regressors. Datasets were first analyzed with the 1st and 3rd order canonical set using a conditionspecific HRF, and the resulting cluster size and maximum t- or F-values were defined as gold standard of detection performance. Secondly, datasets were analyzed with the 1st and 3rd order canonical set, using

Fig. 5. Cortical rat HRF versus human HRF. (A) The cortical rat HRF (green) and its standard deviation (green shading) obtained from 143 individual datasets versus the human HRF (black), which is implemented in SPM (Friston, 2017) by default. The cortical rat HRF proceeds significantly faster than the human HRF. (B) The differences between the cortical rat HRF (dark green) and the human HRF (black) are significantly larger than the differences between the condition-specific HRFs (thin lines). Also the differences between the HRFs of the conditions that are included the cortical rat HRF (green) and the HRFs of the conditions that were excluded from the cortical rat HRF (blue) are substantially smaller than differences to the human HRF. 8

H. Lambers et al.

NeuroImage 208 (2020) 116446

HRFs. We found that the latter differences do not notably affect the statistical analysis of BOLD measurements of rats, while exchanging HRFs of humans and rats has a substantial influence on the analysis. As shown in Fig. 6, analysis using the cortical rat HRF resulted in a substantially better detection performance than the analysis using the human HRF or the 7th order FIR basis set.

canonical model using the human HRF revealed low relative cluster sizes (median ¼ 0) and t-values (median ¼ 0). Also for analysis with the 3rd order model using the human HRF, the median of the relative F-values (0.60) is below the threshold of 0.75. Analysis with the 7th order FIR showed a high coincidence with the cluster size of the analysis using the condition-specific HRF. However, the relative F-values were substantially lower for the FIR analysis (median ¼ 0.56). In contrast, analysis using the rat HRF showed the best detection performance with relative cluster sizes and t- or F-values higher than 0.80.

4.2. Methodological aspects of the HRF calculation The aim of this study was to improve the framework of the GLM for statistical analysis of rodent data acquired under commonly used experimental conditions. Of course, it was not possible to examine all possible conditions that are commonly applied. Our study compared different settings of ten experimental parameters and investigated measurements from 23 groups with different experimental conditions. These cover a substantial range of conditions found in typical fMRI experiments in rat. Yet, it has to be acknowledged that the use of only two rat strains, only two anesthetic regimens, and only two brain regions set limitations for generalization of our findings. A second limitation being the small group size (n ¼ 5–14), which is typical for rodent studies, was overcome by development of the functional t-test. This statistical procedure was able to provide robust and significant results despite the small group sizes.

4. Discussion In this study, BOLD responses of rats obtained under a variety of experimental conditions were investigated. 33 pairs of HRFs from 23 different experimental groups were compared, to identify the impact of individual experimental parameters. Only three parameters, as revealed by the pairwise comparisons, showed significant differences. Our data therefore suggest that the HRF in the rat brain is largely independent of commonly used experimental parameters. Significant differences have been seen between the HRFs of 1 s and 5s stimulation, possibly indicating a deviation from linearity of the BOLD response. Secondly, significant differences were detected between the HRFs for sensory stimulation and optogenetic stimulation using the pan-neuronally expressed opsin Chronos. Finally, significant differences were found for the HRFs of visual activation of subcortical regions and activation of the cortex during electrical paw stimulation. Based on the close concordance of the time courses of the BOLD response under all other tested experimental conditions, a cortical rat HRF for the somatosensory cortex could be derived. The application of the cortical rat HRF substantially improved the detection performance compared to analysis using the standard basis set implemented by default in SPM.

4.3. HRFs for different brain regions It was shown previously for both humans (Handwerker et al., 2004) and rodents (Amirmohseni et al., 2016; Schlegel et al., 2015; Schroeter et al., 2014) that the BOLD response varies across brain regions. We found significant differences between the HRFs in somatosensory cortex and subcortical visual pathways, suggesting that our cortical rat HRF is only valid for analysis of cortical regions. This notion is supported by a study reporting differences between the BOLD responses in the visual cortex and subcortical regions upon visual stimulation (Pawela et al., 2008). In contrast, a study by Bailey et al. (2013) found no significant differences in the impulse response functions of the primary visual cortex and superior colliculus for visual stimulation with a frequency of 1–5 Hz, but reported inherent neurovascular uncoupling for a stimulation frequency of 10 Hz. Similarly, Devonshire et al. (2012), found linear neurovascular coupling in the cortex and non-linear coupling in non-cortical sites. Considering that our cortical rat HRF was only determined on the basis of BOLD responses in the somatosensory cortex, the question arises whether our HRF is also valid for other cortical brain regions. It is

4.1. Rat HRFs differ significantly from the human HRF The cortical rat HRF showed a significantly faster progression than the human HRF from Friston (2017). Peak position und FWHM were consistent with the impulse response functions, which have been characterized in other studies (Bailey et al., 2013; Silva et al., 2007; Zwart et al., 2005). Silva et al. and Zwart et al. attributed this difference to a higher flow velocity in rodent capillary and venous vasculature, and a higher contribution of large vessel effects to the signal. As shown in Fig. 5 B, the differences between cortical rat HRF and the human HRF are by far larger than the differences between the different condition-specific rat

Fig. 6. Detection performance of five different sets of regressors. Detection performances for the 1st order canonical model using the human or rat HRF (left), for the 3rd order canonical model using the human or rat HRF, and the 7th order FIR (right) are shown relative to analysis using the conditionspecific HRF. Top: BOLD maps of one exemplary measurement. Bottom: Boxplots of relative cluster sizes (green), t-values and F-values (blue), which resulted from the analysis of 10 measurements each. Relative values were normalized to the results from analysis using the condition-specific HRF. The red line shows the threshold of 0.75. In the boxplots the central mark represents the median, boxes include the 25th and 75th percentiles, whiskers extend to last data points within the 1.5*interquartile, and outliers are indicated as dots.

9

H. Lambers et al.

NeuroImage 208 (2020) 116446

neuronal activation (Brinker et al., 1999; Herman et al., 2009; Logothetis et al., 2001; Ogawa et al., 2000), suggesting a nonlinearity between stimulation and neuronal activity. There are models which include nonlinearities of the BOLD response (Friston et al., 1998, 2000; Wager et al., 2005). However, these models require more parameters, which can lead to overfitting and unphysiological results. Our cortical rat HRF, which includes the effects of neuronal adaptation, therefore appears to be a feasible alternative providing a good approximation for analyzing rodent fMRI data acquired at common experimental conditions.

commonly assumed that the cellular structure is very similar throughout the cortex. Recently, it was shown that the density of penetrating arterioles is nearly uniform across different cortical regions (Adams et al., 2018). For the human cortex, temporal parameters of the BOLD response upon audiovisual stimulation were found to remain fairly stable across the neocortex, despite small but significant variations (Taylor et al., 2018). We, therefore, assume that our rat HRF provides a better approximation than the human HRF for analysis of the BOLD response in the entire rat cortex. 4.4. Linearity of the BOLD response

4.5. Impact of anesthetic regimen and pain The procedures for extracting the HRF implemented in this manuscript are based on the assumption that the hemodynamic BOLD response to neuronal activation can be described by a linear model, and that a linear relationship exists between stimulation paradigm and the induced neuronal activation. Accordingly, the HRF should be independent of the stimulation duration. However, we have observed a significant difference between the HRFs for 1 s and 5 s stimulation, possibly indicating a deviation from the linearity of the BOLD response. One possible cause for this deviation may be neuronal adaptation. Local field potential measurements have previously shown that neuronal adaptation occurred after less than 5 s for a stimulation frequency of 9 Hz (Schmid et al., 2016). Thus, the shape of the time course of neuronal activation may deviate from the box car shape of the stimulation paradigm for longer stimulation durations, possibly explaining the difference between 1 s and 5 s stimulation. Consequently, our cortical rat HRF inherently includes effects of neuronal adaptation which occur for longer stimulation durations. We think that our cortical rat HRF therefore provides an appropriate approximation for most paradigms using stimulation durations of 5 s or longer. Additionally, an HRF for measurements with 1 s stimulation duration has been determined (1-s HRF, Suppl. Table 3) and can be implemented in SPM (Suppl. C). The 1-s HRF is most likely the best choice for analysis of measurements with short stimulation durations or for analysis of measurements for which the actual shape of the neuronal response is known. Assuming linearity of the BOLD response, the HRF should be independent of stimulation frequency. However, previous studies have reported a strong dependence of the hemodynamic response on stimulation frequency (Bailey et al., 2013; Kim et al., 2010; Schmid et al., 2016; Zhao et al., 2008). Zhao et al. (2008) and Schmid et al. (2016) both reported that the amplitude of the BOLD response was the highest at 9 Hz electrical paw stimulation and small at either lower (1 Hz and 3 Hz) or higher frequencies (12–18 Hz) under medetomidine anesthesia. These differences may also be due to neuronal adaptation. Kim et al. (2010) reported that the differences in evoked field potentials were prominent between low frequencies (1–3 Hz) and higher frequencies (10–12 Hz) for electrical paw stimulation. We therefore analyzed the shape of the HRF only for frequencies between 5 Hz and 12 Hz for paw stimulation and found no substantial differences in this narrow range. Optogenetic stimulation, which directly drives neuronal spiking, resulted in robust BOLD responses for stimulation frequencies up to 150 Hz. No frequency-dependent effects were observed, suggesting that the vascular response is largely independent of stimulation frequency. Our study did not investigate whether the length of individual pulses had an influence on the linearity of the BOLD response. A study by Ngai et al. (1995) has shown that extending the pulse length by a factor of ten strongly influenced pial arteriolar diameter and blood flow. Therefore, it has to be assumed that the BOLD response is also dependent on the pulse length. The observed deviations from linearity of the BOLD response are consistent with several human studies (Huettel and McCarthy, 2000; Miller et al., 2001; Pfeuffer et al., 2003; Vazquez and Noll, 1998; Wager et al., 2005; Yes¸ilyurt et al., 2008), in which such deviations have partly been attributed to neuronal adaptation. However, several rodent studies showed that the BOLD response changed approximately linearly with

The shape of the BOLD response is also influenced by other factors. A substantial impact of the anesthetic regimen on neurovascular coupling, brain state and neuronal responses has been recognized (Masamoto and Kanno, 2012). For mice it was shown that temporal characteristics of the BOLD responses depend on anesthesia (Schlegel et al., 2015; Schroeter et al., 2014). In rats a complex impact of isoflurane on both neurovascular and neuronal response has been shown. This includes a dose dependent effect of the amplitude of the cerebral blood flow (Masamoto et al., 2009) and a reduction of the BOLD amplitude (van Alst et al., 2019). However, according to our analysis the temporal shape of the HRF showed no difference after the addition of low doses of isoflurane to medetomidine. Furthermore, according to previous reports pain influences the shape of the BOLD response. In different rat models of pain a biphasic shape of the BOLD response has been observed (Amirmohseni et al., 2016). In contrast, our analysis of the cortical HRF showed no difference between innocuous and noxious stimulation. The impact of the complex phenomenon of pain on the neurovascular response in different brain regions remains a matter of further investigation and it is questionable whether the BOLD response can still be described by the linear model under painful conditions. 4.6. Sex related differences of the HRF For humans a gender related difference of the BOLD amplitude has been reported (Kaufmann et al., 2001; Levin et al., 1998). A possible cause of this difference is a higher global cerebral blood flow in women (Gur and Gur, 1990). For rats, a sex related difference of the response of cerebral blood flow has been reported for hypercapnia, but not for normal physiological conditions (Ances et al., 2001). However, only differences in amplitude were observed and no differences in the temporal evolution of the BOLD response were detected. For the method used to determine the HRF in this study, only the temporal evolution of the BOLD response is relevant, whereas amplitude is not. Therefore, our result of no significant difference between the HRFs of female and male rats is consistent with the literature. However, it is possible that the BOLD maps of male and female animals cannot be compared quantitatively. 4.7. Impact of stimulation pathways Our study provides only limited insight into the impact of neuronal processing on the HRF. In a study by Liu et al. (2015), optogenetic stimulation of the central thalamus resulted in negative signal changes in cortex upon 10 Hz stimulation, while stimulation with 40 Hz and 100 Hz resulted in a positive BOLD signal. We did not observe an impact of the stimulation frequency on the HRF, nor did any of the sensory inputs result in different HRFs. Optogenetic stimulation using the opsin C1V1TT, expressed only in excitatory neurons, also did not result in an altered HRF. However, for optogenetic stimulation using the opsin Chronos expressed in all neurons, including inhibitory interneurons, a significantly different HRF was observed. This difference supports the recently published notion that inhibitory interneurons are crucial in shaping the BOLD response (Aksenov et al., 2019). Further, our study only considered positive signal changes of the BOLD responses and did not assess whether 10

H. Lambers et al.

NeuroImage 208 (2020) 116446

negative signal changes can also be modeled by the cortical rat HRF. Such analysis, however, may require alternative models, since the mechanism of negative BOLD is a matter of debate (Schridde et al., 2008; Shih et al., 2009; Shmuel et al., 2002).

compared to sensory stimulation. Neuroimage 164, 144–154. https://doi.org/ 10.1016/j.neuroimage.2016.12.059. Amirmohseni, S., Segelcke, D., Reichl, S., Wachsmuth, L., G€ orlich, D., Faber, C., PogatzkiZahn, E., 2016. Characterization of incisional and inflammatory pain in rats using functional tools of MRI. Neuroimage 127, 110–122. https://doi.org/10.1016/ j.neuroimage.2015.11.052. Ances, B.M., Greenberg, J.H., Detre, J.A., 2001. Sex differences in the cerebral blood flow response after brief hypercapnia in the rat. Neurosci. Lett. 304 (1–2), 57–60. Bailey, C.J., Sanganahalli, B.G., Herman, P., Blumenfeld, H., Gjedde, A., Hyder, F., 2013. Analysis of time and space invariance of BOLD responses in the rat visual system. Cerebr. Cortex 23 (1), 210–222. https://doi.org/10.1093/cercor/bhs008. New York, N.Y. : 1991. Boynton, G.M., Engel, S.A., Glover, G.H., Heeger, D.J., 1996. Linear systems analysis of functional magnetic resonance imaging in human V1. J. Neurosci. 16 (13), 4207–4221. https://doi.org/10.1523/JNEUROSCI.16-13-04207.1996. Brinker, G., Bock, C., Busch, E., Krep, H., Hossmann, K.-A., Hoehn-Berlage, M., 1999. Simultaneous recording of evoked potentials and T*2-weighted MR images during somatosensory stimulation of rat. Magn. Reson. Med. 41 (3), 469–473. https:// doi.org/10.1002/(SICI)1522-2594 (199903)41:3<469:AID-MRM7>3.0.CO;2-9. Cohen, E.R., Ugurbil, K., Kim, S.-G., 2002. Effect of basal conditions on the magnitude and dynamics of the blood oxygenation level-dependent fMRI response. J. Cereb. Blood Flow Metab. 22 (9), 1042–1053. Devonshire, I.M., Papadakis, N.G., Port, M., Berwick, J., Kennerley, A.J., Mayhew, J.E.W., Overton, P.G., 2012. Neurovascular coupling is brain region-dependent. Neuroimage 59 (3), 1997–2006. https://doi.org/10.1016/j.neuroimage.2011.09.050. Friston, K.J., 2017. SPM12. Wellcome Trust Centre for Neuroimaging. Friston, K.J., Holmes, A.P., Worsley, K.J., Poline, J.-P., Frith, C.D., Frackowiak, R.S.J., 1995. Statistical parametric maps in functional imaging: a general linear approach. Hum. Brain Mapp. (2), 189–210. Friston, K.J., Jezzard, P., Turner, R., 1994. Analysis of functional MRI time-series. Hum. Brain Mapp. 1 (2), 153–171. https://doi.org/10.1002/hbm.460010207. Friston, K.J., Josephs, O., Rees, G., Turner, R., 1998. Nonlinear event-related responses in fMRI. Magn. Reson. Med. 39 (1), 41–52. https://doi.org/10.1002/mrm.1910390109. Friston, K.J., Mechelli, A., Turner, R., Price, C.J., 2000. Nonlinear responses in fMRI: the Balloon model, Volterra kernels, and other hemodynamics. Neuroimage 12 (4), 466–477. https://doi.org/10.1006/nimg.2000.0630. Genovese, C.R., 2000. A bayesian time-course model for functional magnetic resonance imaging data. J. Am. Stat. Assoc. 95 (451), 691–703. https://doi.org/10.1080/ 01621459.2000.10474253. Gur, R.E., Gur, R.C., 1990. Gender differences in regional cerebral blood flow. Schizophr. Bull. 16 (2), 247–254. Handwerker, D.A., Ollinger, J.M., D’Esposito, M., 2004. Variation of BOLD hemodynamic responses across subjects and brain regions and their effects on statistical analyses. Neuroimage 21 (4), 1639–1651. https://doi.org/10.1016/ j.neuroimage.2003.11.029. Henson, R., Friston, K.J., 2007. Convolution models for fMRI. In: Friston, K.J. (Ed.), Statistical Parametric Mapping: the Analysis of Functional Brain Images, first ed. Elsevier/Academic Press, Amsterdam, Boston. Herman, P., Sanganahalli, B.G., Blumenfeld, H., Hyder, F., 2009. Cerebral oxygen demand for short-lived and steady-state events. J. Neurochem. 109 (Suppl. 1), 73–79. https:// doi.org/10.1111/j.1471-4159.2009.05844.x. Huettel, S.A., McCarthy, G., 2000. Evidence for a refractory period in the hemodynamic response to visual stimuli as measured by MRI. Neuroimage 11 (5 Pt 1), 547–553. https://doi.org/10.1006/nimg.2000.0553. Jones, M., Berwick, J., Hewson-Stoate, N., Gias, C., Mayhew, J., 2005. The effect of hypercapnia on the neural and hemodynamic responses to somatosensory stimulation. Neuroimage 27 (3), 609–623. https://doi.org/10.1016/ j.neuroimage.2005.04.036. Kaufmann, C., Elbel, G.-K., G€ ossl, C., Pütz, Benno, Auer, Dorothee P., 2001. Frequency dependence and gender effects in visual cortical regions involved in temporal frequency dependent pattern processing. Hum. Brain Mapp. 14 (1), 28–38. Kim, T., Masamoto, K., Fukuda, M., Vazquez, A., Kim, S.-G., 2010. Frequency-dependent neural activity, CBF, and BOLD fMRI to somatosensory stimuli in isofluraneanesthetized rats. Neuroimage 52 (1), 224–233. https://doi.org/10.1016/ j.neuroimage.2010.03.064. Levin, J.M., Ross, M.H., Mendelson, J.H., Mello, N.K., Cohen, Bruce M., Renshaw, Perry F., 1998. Sex differences in blood-oxygenation-level-dependent functional MRI with primary visual stimulation. Am. J. Psychiatry 155 (3), 434–436. Lindquist, M.A., Meng Loh, J., Atlas, L.Y., Wager, T.D., 2009. Modeling the hemodynamic response function in fMRI: efficiency, bias and mis-modeling. Neuroimage 45 (Suppl. l), S187–S198. https://doi.org/10.1016/j.neuroimage.2008.10.065, 1. Liu, J., Duffy, B.A., Bernal-Casas, D., Fang, Z., Lee, J.H., 2017. Comparison of fMRI analysis methods for heterogeneous BOLD responses in block design studies. Neuroimage 147, 390–408. https://doi.org/10.1016/j.neuroimage.2016.12.045. Liu, J., Lee, H.J., Weitz, A.J., Fang, Z., Lin, P., Choy, M., Fisher, R., Pinskiy, V., Tolpygo, A., Mitra, P., Schiff, N., Lee, J.H., 2015. Frequency-selective control of cortical and subcortical networks by central thalamus. eLife 4. https://doi.org/ 10.7554/eLife.09215.001. Logothetis, N.K., Pauls, J., Augath, M., Trinath, T., Oeltermann, A., 2001. Neurophysiological investigation of the basis of the fMRI signal. Nature 412 (6843), 150–157. https://doi.org/10.1038/35084005. Masamoto, K., Fukuda, M., Vazquez, A., Kim, S.-G., 2009. Dose-dependent effect of isoflurane on neurovascular coupling in rat cerebral cortex. Eur. J. Neurosci. 30 (2), 242–250. https://doi.org/10.1111/j.1460-9568.2009.06812.x.

5. Conclusions We have determined a cortical rat HRF based on 143 BOLD measurements in the primary somatosensory cortex of 76 animals, which proceeds significantly faster than a human HRF. GLM analysis of rat fMRI data with the appropriate rat HRF implemented in SPM showed better detection performance as compared to using the human HRF or the FIR approach. Strictly, our cortical rat HRF is valid only for the range of experimental conditions tested, which were male or female Fischer or SD rats, mechanical or electrical paw stimulation (both innocuous and noxious) or optogenetic stimulation, using frequencies from 5 Hz up to 12 Hz and durations of 5 s or 10 s, under medetomidine anesthesia or medetomidine supplemented with 0.7 % isoflurane, using ventilation or spontaneous breathing animals. Yet, the small differences between the different condition-specific HRFs as compared to the human HRF, suggest that our cortical rat HRF is a good choice for analysis of rat fMRI data. If experimental conditions differ from those investigated here, caution is warranted, and results should be verified by either deriving a conditions-specific HRF or using the FIR approach. This is in particular valid for subcortical regions and stimulation with very short pulses, and may also apply for other anesthetic regimens. Author contributions HL and CF designed the study. TvA, FA, LW and HL conducted fMRI measurements. MS developed the functional t-test. HL implemented the voxelwise u-test and the fitting procedure and analyzed the data. All authors discussed data and edited the manuscript. HL and CF wrote the manuscript. Ethics statement Animal experiments were carried out according to the German Tierschutzgesetz and were approved by the Landesamt für Natur, Umwelt und Verbraucherschutz of Nordrhein-Westfalen, Germany (approval IDs 87-51.04.2010.A274, 87-51.04.2011.A065, 84-02.04.2014.A347 and 84-02.04.2015.A427). Declaration of competing interest None. Acknowledgements We thank Zakhar Kabluchko for pointing out the functional t-test to us, and Saeedeh Amirmohseni and Florian Schmid for providing their data. This work was funded by the DFG (Fa474/5-1 and Fa4747/6-1). Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.neuroimage.2019.116446. References Adams, M.D., Winder, A.T., Blinder, P., Drew, P.J., 2018. The pial vasculature of the mouse develops according to a sensory-independent program. Sci. Rep. 8 (1), 9860. https://doi.org/10.1038/s41598-018-27910-3. Aksenov, D.P., Li, L., Miller, M.J., Wyrwicz, A.M., 2019. Role of the inhibitory system in shaping the BOLD fMRI response. Neuroimage 201. Albers, F., Schmid, F., Wachsmuth, L., Faber, C., 2018. Line scanning fMRI reveals earlier onset of optogenetically evoked BOLD response in rat somatosensory cortex as 11

H. Lambers et al.

NeuroImage 208 (2020) 116446 Schroeter, A., Schlegel, F., Seuwen, A., Grandjean, J., Rudin, M., 2014. Specificity of stimulus-evoked fMRI responses in the mouse: the influence of systemic physiological changes associated with innocuous stimulation under four different anesthetics. Neuroimage 94, 372–384. https://doi.org/10.1016/j.neuroimage.2014.01.046. Shan, Z.Y., Wright, M.J., Thompson, P.M., McMahon, K.L., Blokland, G.G.A.M., de Zubicaray, G.I., Martin, N.G., Vinkhuyzen, A.A.E., Reutens, D.C., 2014. Modeling of the hemodynamic responses in block design fMRI studies. J. Cereb. Blood Flow Metab. : Off. J. Int. Soc. Cereb. Blood Flow Metab. 34 (2), 316–324. https://doi.org/ 10.1038/jcbfm.2013.200. Shih, Y.-Y.I., Chen, C.-C.V., Shyu, B.-C., Lin, Z.-J., Chiang, Y.-C., Jaw, F.-S., Chen, Y.-Y., Chang, C., 2009. A new scenario for negative functional magnetic resonance imaging signals: endogenous neurotransmission. J. Neurosci. : Off. J. Soc. Neurosci. 29 (10), 3036–3044. https://doi.org/10.1523/JNEUROSCI.3447-08.2009. Shmuel, A., Yacoub, E., Pfeuffer, J., van de Moortele, P.-F., Adriany, G., Hu, X., Ugurbil, K., 2002. Sustained negative BOLD, blood flow and oxygen consumption response and its coupling to the positive response in the human brain. Neuron 36 (6), 1195–1210. https://doi.org/10.1016/S0896-6273(02)01061-9. Silva, A.C., Koretsky, A.P., Duyn, J.H., 2007. Functional MRI impulse response for BOLD and CBV contrast in rat somatosensory cortex. Magn. Reson. Med. 57 (6), 1110–1118. https://doi.org/10.1002/mrm.21246. Taylor, A.J., Kim, J.H., Ress, D., 2018. Characterization of the hemodynamic response function across the majority of human cerebral cortex. Neuroimage 173, 322–331. https://doi.org/10.1016/j.neuroimage.2018.02.061. van Alst, T.M., Wachsmuth, L., Datunashvili, M., Albers, F., Just, N., Budde, T., Faber, C., 2019. Anesthesia differentially modulates neuronal and vascular contributions to the BOLD signal. Neuroimage 195, 89–103. https://doi.org/10.1016/ j.neuroimage.2019.03.057. Vazquez, A.L., Noll, D.C., 1998. Nonlinear aspects of the BOLD response in functional MRI. Neuroimage 7 (2), 108–118. https://doi.org/10.1006/nimg.1997.0316. Wager, T.D., Vazquez, A., Hernandez, L., Noll, D.C., 2005. Accounting for nonlinear BOLD effects in fMRI: parameter estimates and a model for prediction in rapid event-related studies. Neuroimage 25 (1), 206–218. https://doi.org/10.1016/ j.neuroimage.2004.11.008. Worsley, K.J., Friston, K.J., 1995. Analysis of fMRI time-series revisited - again. Neuroimage (2), 173–181. Yes¸ilyurt, B., U gurbil, K., Uluda g, K., 2008. Dynamics and nonlinearities of the BOLD response at very short stimulus durations. Magn. Reson. Imaging 26 (7), 853–862. https://doi.org/10.1016/j.mri.2008.01.008. Zhao, F., Zhao, T., Zhou, L., Wu, Q., Hu, X., 2008. BOLD study of stimulation-induced neural activity and resting-state connectivity in medetomidine-sedated rat. Neuroimage 39 (1), 248–260. https://doi.org/10.1016/j.neuroimage.2007.07.063. de Zwart, J.A., Silva, A.C., van Gelderen, P., Kellman, P., Fukunaga, M., Chu, R., Koretsky, A.P., Frank, J.A., Duyn, J.H., 2005. Temporal dynamics of the BOLD fMRI impulse response. Neuroimage 24 (3), 667–677. https://doi.org/10.1016/ j.neuroimage.2004.09.013.

Masamoto, K., Kanno, I., 2012. Anesthesia and the quantitative evaluation of neurovascular coupling. J. Cereb. Blood Flow Metab. : Off. J. Int. Soc. Cereb. Blood Flow Metab. 32 (7), 1233–1247. https://doi.org/10.1038/jcbfm.2012.50. Miller, K.L., Luh, W.-M., Liu, T.T., Martinez, A., Obata, T., Wong, E.C., Frank, L.R., Buxton, R.B., 2001. Nonlinear temporal dynamics of the cerebral blood flow response. Hum. Brain Mapp. 13 (1), 1–12. https://doi.org/10.1002/hbm.1020. Ngai, A.C., Meno, J.R., Winn, H.R., 1995. Simultaneous measurements of pial arteriolar diameter and laser-Doppler flow during somatosensory stimulation. J. Cereb. Blood Flow Metab. (15), 124–127. Ogawa, S., Lee, T.-M., Stepnoski, R., Chen, W., Zhu, X.-H., Ugurbil, K., 2000. An approach to probe some neural systems interaction by functional MRI at neural time scale down to milliseconds. Proc. Natl. Acad. Sci. 97 (20), 11026–11031. https://doi.org/ 10.1073/pnas.97.20.11026. Pawela, C.P., Hudetz, A.G., Ward, B.D., Schulte, M.L., Li, R., Kao, D.S., Mauck, M.C., Cho, Y.R., Neitz, J., Hyde, J.S., 2008. Modeling of region-specific fMRI BOLD neurovascular response functions in rat brain reveals residual differences that correlate with the differences in regional evoked potentials. Neuroimage 41 (2), 525–534. https://doi.org/10.1016/j.neuroimage.2008.02.022. Peng, S.-L., Chen, C.-M., Huang, C.-Y., Shih, C.-T., Huang, C.-W., Chiu, S.-C., Shen, W.-C., 2019. Effects of hemodynamic response function selection on rat fMRI statistical analyses. Front. Neurosci. 13, 400. https://doi.org/10.3389/fnins.2019.00400. Pfeuffer, J., McCullough, J.C., van de Moortele, P.-F., Ugurbil, K., Hu, X., 2003. Spatial dependence of the nonlinear BOLD response at short stimulus duration. Neuroimage 18 (4), 990–1000. https://doi.org/10.1016/S1053-8119(03)00035-1. Ramsay, J.O., Hooker, G., Spencer, G., 2009. Functional Data Analysis with R and MATLAB, first ed. Springer, Dordrecht, Heidelberg, London, New York. Schlegel, F., Schroeter, A., Rudin, M., 2015. The hemodynamic response to somatosensory stimulation in mice depends on the anesthetic used: implications on analysis of mouse fMRI data. Neuroimage 116, 40–49. https://doi.org/10.1016/ j.neuroimage.2015.05.013. Schmid, F., Wachsmuth, L., Albers, F., Schwalm, M., Stroh, A., Faber, C., 2017. True and apparent optogenetic BOLD fMRI signals. Magn. Reson. Med. 77 (1), 126–136. https://doi.org/10.1002/mrm.26095. Schmid, F., Wachsmuth, L., Schwalm, M., Prouvot, P.-H., Jubal, E.R., Fois, C., Pramanik, G., Zimmer, C., Faber, C., Stroh, A., 2016. Assessing sensory versus optogenetic network activation by combining (o)fMRI with optical Ca2þ recordings. J. Cereb. Blood Flow Metab. : Off. J. Int. Soc. Cereb. Blood Flow Metab. 36 (11), 1885–1900. https://doi.org/10.1177/0271678X15619428. Schneider, C.A., Rasband, W.S., Eliceiri, K.W., 2012. NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9 (7), 671–675. https://doi.org/10.1038/nmeth.2089. Schridde, U., Khubchandani, M., Motelow, J.E., Sanganahalli, B.G., Hyder, F., Blumenfeld, H., 2008. Negative BOLD with large increases in neuronal activity. Cerebr. Cortex 18 (8), 1814–1827. https://doi.org/10.1093/cercor/bhm208. New York, N.Y. : 1991.

12