Abstract
Energetic basis of neural activity provides a solid foundation for noninvasive neuroimaging with calibrated functional magnetic resonance imaging (fMRI). Calculating dynamic changes in cerebral oxidative energy utilization (CMRO2) is limited by uncertainties about whether or not the conventional blood oxygenation level-dependent (BOLD) model can be applied transiently using multimodal measurements of blood flow (CBF) and volume (CBV) that affect the BOLD signal. A prerequisite for dynamic calibrated fMRI is testing the linearity of multimodal signals within a temporal regimen, as assessed by signal strength (i.e., both intensity and width). If each hyperemic component (BOLD, CBV, CBF) is demonstrated to be linear with neural activity under various experimental conditions, then the respective transfer functions generated by deconvolution with neural activity should be time invariant and thus could potentially be used for calculating CMRO2 transients. Hyperemic components were investigated at 11.7 T in α-chloralose-anesthetized rats and combined with electrophysiological recordings of local field potential (LFP) and multiunit activity (MUA) from the cortex during forepaw stimulation, in which stimulus number and frequency were varied. Although relationships between neural activity and stimulus features ranged from linear to nonlinear, associations between hyperemic components and neural activity were linear. Specific to each hyperemic component, a universal transfer function (with LFP or MUA) yielded predictions in agreement with experimental measurements. The results identified a component of the BOLD signal that can be attributed to significant changes in CMRO2, even for temporal events separated by <200 ms.
Introduction
Continuous need for oxidative energy by neurons and astrocytes is critical for normal brain function (Hyder et al., 2006; Buzsáki et al., 2007). Because glucose oxidation provides the large majority of the demanded energy (Riera et al., 2008), it is desirable to map oxygen consumption (CMRO2) as close approximation for brain work (Hyder et al., 2002). Measurements of tissue oxygen tension (pO2) is used to infer local dynamic uncoupling between oxygen delivery and its subsequent utilization (Ances et al., 2001; Thompson et al., 2003; Caesar et al., 2008). Positron emission tomography (PET) can measure changes in CMRO2 (Vafaee and Gjedde, 2000; Ito et al., 2005). Although 13C magnetic resonance spectroscopy (MRS) has lower spatial resolution compared with magnetic resonance imaging (MRI) and PET, it can distinguish energetics of neurons and astrocytes (Rothman et al., 1999; Shulman et al., 2004).
Functional MRI (fMRI) maps brain activity noninvasively, but the blood oxygenation level-dependent (BOLD) signal is an indirect indicator of neural activity because it measures the hyperemic response, which has both hemodynamic and metabolic contributions (Ogawa et al., 1993). To extract CMRO2 from fMRI, the BOLD signal is calibrated by combining changes in blood flow (CBF) and volume (CBV) (Hoge and Pike, 2001; Kida and Hyder, 2006). Steady-state calibrated fMRI has been used to calculate CMRO2 changes in humans (Davis et al., 1998; Kastrup et al., 2002; Lu and van Zijl, 2005; Chiarelli et al., 2007) and rodents (Smith et al., 2002; Liu et al., 2004; Kida et al., 2006; Maandag et al., 2007; Shen et al., 2008). Similarly, optical imaging techniques, which measure changes in CBV (from oxyhemoglobin and deoxyhemoglobin) and CBF concurrently, can be used to estimate CMRO2 variations (Jones et al., 2001; Devor et al., 2003; Durduran et al., 2004; Sheth et al., 2004; Dunn et al., 2005). Although optical imaging has better temporal resolution than fMRI, mapping is limited to superficial regions because of poor light penetration in tissue and/or increased light scattering.
Because fMRI has good spatiotemporal resolution and allows ample brain coverage, it would be attractive to calculate CMRO2 transients with calibrated fMRI. However, a problem encountered with dynamic calibrated fMRI is uncertainty about whether the conventional BOLD model can be applied transiently. Furthermore, there is also the issue of lack of a validation for the predicted transient energetics. An advantage of steady-state calibrated fMRI is that calculated ΔCMRO2 can be validated by comparison with results from MRS or PET (Hyder et al., 2001; Zhu et al., 2005; Kudomi et al., 2007). To extract CMRO2 transients, we used electrophysiological and hyperemic components measured from the rat somatosensory cortex during forepaw stimulation. By varying stimulus number and frequency within a narrow temporal regimen, if linearity is experimentally demonstrated between hyperemic components and neural activity, then it may be possible to use calibrated fMRI in a dynamic manner, because then the respective transfer functions (of BOLD, CBV, CBF) would be linear and time invariant (Boynton et al., 1996; Friston et al., 1998; Ances et al., 2000).
Materials and Methods
Animal preparation.
All procedures were performed in accordance with protocols approved by the ethical committee of Yale University School of Medicine Institutional Animal Care and Use Committee and in agreement with the National Institutes of Health Guide for the Care and Use of Laboratory Animals. All experiments were conducted on artificially ventilated (1–2% halothane during surgery, plus 70% N2O/30% O2) adult male rats (n = 26; Sprague Dawley; 200–300 g; Charles River). Twelve rats were used for simultaneous electrophysiology and laser-Doppler flowmetry (LDF) experiments, and 14 were used for multimodal fMRI studies. Of the 12 rats used for combined electrophysiology and LDF measurements, we averaged data from six subjects from which histology data were obtained. Both BOLD and CBV experiments were performed on eight subjects, whereas six more rats were studied for CBV experiments only. Therefore, we compared averages of six and eight rats from the two respective groups. A femoral artery was cannulated for monitoring physiologic parameters (pCO2, pO2, pH, blood pressure) and a femoral vein was cannulated for infusion of MRI contrast agent for measuring relative changes in CBV. Physiological variables (pCO2, pO2, pH, blood pressure) were measured and kept within normal limits (34 ± 3 mmHg, 112 ± 18 mmHg, 7.35 ± 0.04, 112 ± 16 mmHg, respectively). An intraperitoneal catheter was used for administration of α-chloralose (∼40 mg · kg−1 · h−1) and d-tubocurarine chloride (∼0.3 mg · kg−1 · h−1). After all the surgical procedures, halothane was discontinued and anesthesia was maintained with α-chloralose. The animal core temperature was measured by rectal probe and controlled by a heated water blanket maintained at ∼37°C.
Stimulus presentation and paradigm.
Two subcutaneously placed copper needles were inserted into the contralateral forepaw (between the second and fourth digits) and all snout whiskers were shaved to avoid contaminating somatosensory signals. All stimulus presentation was controlled by a μ-1401 analog-to-digital converter unit (CED) running custom-written script to provide 0.3 ms duration pulses with 2 mA amplitude by an isolation unit (WPI). The interpulse intervals ranged from 167, 333, and 667 ms, which corresponded to stimulus frequencies of 6, 3, and 1.5 Hz, respectively, at which strong BOLD signals have been reported (Sanganahalli et al., 2008). The stimulus number was varied from 1 to 4. A 5 min resting period was allowed in-between experimental trials. The stimulus paradigm is shown in supplemental Figure 1 (available at www.jneurosci.org as supplemental material).
Electrophysiology and LDF experiments.
Rats (n = 12) were mounted on a stereotaxic frame (David Kopf Instruments) placed on a vibration-free table inside a Faraday cage. The scalp and the galea aponeurotica were removed and small burr holes were drilled for insertion of a dual-sensor device consisting of adjacent electrical and optical probes to simultaneously measure neural and CBF signals, respectively (Schridde et al., 2008). Neural signals were measured with a high impedance tungsten microelectrode (2–4 MΩ; FHC) and CBF was measured by an LDF probe (805 nm; Oxyflow 4000; Oxford Optronix). The tip of the microelectrode was placed ∼300 μm below the LDF probe to ensure overlapping sampling volumes from both measurements (Trübel et al., 2006). Recordings were localized to the middle cortical layers of the forelimb somatosensory cortex (4.4 mm lateral, 1.0 mm anterior to bregma, 0.9 ± 0.1 mm depth from cortical surface) and confirmed histologically (Schridde et al., 2008) to be compared with fMRI signals at the same location (supplemental Fig. 2, available at www.jneurosci.org as supplemental material). Electrical and optical signals were digitized with CED μ-1401 using Spike 2 software (CED) at 20 kHz and 50 Hz, respectively.
Local field potential (LFP) and multiunit activity (MUA) signals were obtained by splitting the electrical signal into low (<150 Hz) and high (0.4–10 kHz) frequency bands using a Butterworth filter (24 dB/oct attenuation). LDF signals were measured with respect to the prestimulus baseline. The intensity of the LDF signal (with 3 Hz stimulation) was calibrated to CBF data collected with arterial spin tagging (ASL) contrast MRI (3 Hz, 2 mA, 0.3 ms, 96 pulses) as previously shown by Kida et al. (2004). The high temporal resolution neural signals were integrated into 0.02 s running bins (equivalent to 50 Hz sampling) for comparison with the lower temporal resolution BOLD and CBV data (see below).
Multimodal fMRI experiments.
All fMRI data (n = 14) were obtained on a modified 11.7 T Bruker horizontal-bore spectrometer (Bruker) using a home-built 1H surface coil radiofrequency probe (1.4 cm diameter). Details of fMRI measurements for BOLD and CBV were discussed previously (Sanganahalli et al., 2008; Herman et al., 2009a). Briefly, BOLD signal was acquired with echoplanar imaging (EPI) with sequential sampling (Hyder et al., 1995) using gradient-echo contrast with repetition and echo times of 1000 and 13 ms, respectively. The field of view was 2.56 cm in both axes with in-plane matrix of 64 × 64 and slice thickness of 2 mm. In the same subjects after completion of the BOLD experiments, CBV signal was measured by the same EPI parameters but in the presence of an intravenous injection of iron oxide nanocolloid particles (Combidex; 15 mg/kg; AMAG), which has a very long half-life in blood circulation of rodents (Kida et al., 2000). The BOLD effect was removed from the CBV-weighted signal (Kida et al., 2004; Lu et al., 2007). Because with the MRI contrast agent on board stimulation-induced decrease in the EPI signal reflected an increase in volume, all CBV-weighted fMRI data were multiplied by −1 to reflect positive changes in CBV during functional hyperemia.
All fMRI data were subjected to a translational movement criterion using a center-of-mass analysis (Chahboune et al., 2007). If either center-of-mass value in a series deviated by >25% of a pixel, the entire dataset was discarded from additional analysis. Single-run data were used to create activation maps and time courses. Activation foci for fractional change in BOLD maps were obtained by applying Student's t test comparison of resting and stimulated data. Individual time courses of activated voxels were averaged across many subjects, and mean ± SD was used to provide the BOLD and CBV responses.
Convolution analysis.
Insights into a system can be gained by establishing connection between the input and the output signals by convolution theory. The input signal, i(t), when convolved with a transfer function, h(t), produces an output signal, r(t). In our experiments, we measured both i(t) and r(t) with high spatiotemporal resolution. LFP or MUA was used as the neural signal, whereas BOLD, CBV, and CBF each was used as an independent hyperemic response.
The transfer function, h(t), can be achieved by a deconvolution between r(t) and i(t). If the frequency domain representation of i(t), r(t), and h(t) are given by I(ω), R(ω), and H(ω), respectively, then it can be shown that
where t and ω are time and frequency domains. Inverse Fourier transform of H(ω) gives h(t).
The gamma variate function is widely used for transfer function modeling. We used a slightly different form of the original equation (Madsen, 1992), which is given by the following:
where y0 is the baseline shift, ym is the intensity of the peak, tm is the peak time, α is related to the rise and fall times, and the appearance time (t0) is defined as t = 0. The parameters in Equation 2 are independent from each other.
To find the best applicable transfer function, we used a least-square mean fitting (Gauss-Newton) method (Matlab), which had three steps (supplemental Fig. 3, available at www.jneurosci.org as supplemental material): (1) a transfer function was created with initial parameters, (2) a transfer function was convolved on the input, and (3) a residual signal was created by differencing the predicted and measured signals. If the predicted signal was significantly different from the measured signal, then the parameters of the transfer function were changed and the process was continued from the first step. The fitting process was completed usually within several hundred iterations and the last applied parameter set was used to define the transfer function. We found the residuals acceptable if all of their values were within the range of uncertainty of the corresponding measured signal, as represented by the SD about the mean.
Our goal was to characterize the somatosensory cortex behavior by a temporal universal transfer function so that the signals can be predicted for a wide range of stimulus conditions. Transfer functions were calculated separately for BOLD, CBV, and CBF signals using the measured LFP and MUA data. If the residuals were within the range of ±SD of each of the measured responses [e.g., as shown for BOLD and LFP in supplemental Fig. 3 (available at www.jneurosci.org as supplemental material) with two stimulus pulses] for the different stimulus conditions, then a time-invariant universal transfer function was obtained.
Estimation of CMRO2 changes.
Changes in CMRO2 were calculated using multimodal measurements of changes in BOLD signal (ΔS/S), blood flow (ΔCBF/CBF), and blood volume (ΔCBV/CBV) from the relationship originally described by Ogawa et al. (1993) as follows:
where A represents a BOLD calibration parameter (Hyder et al., 2001). Application of Equation 3 depended, however, on demonstrating linearity and time invariance (see Results). Because analyses using A values of 0.4–0.6 did not appreciably change CMRO2 predictions (Schridde et al., 2008), a value of 0.5 was used for all calculations.
Measured values of CBF, CBV, and BOLD signal along with their respective SDs (σ) were used to analytically describe the SD of CMRO2 data. For an analytical description of the error propagation (Ku, 1966), it is helpful to rewrite Equation 3 in the following manner:
where m = (1 + ΔCMRO2/CMRO2), f = (1 + ΔCBF/CBF), v = (1 + ΔCBV/CBV), and s = (1 + ΔS/S). If the measured signals are independent from each other, the error propagation is restricted only by the derivatives of σ for each individual component. However, in our case, the measured data (CBF, CBV, BOLD) could be highly correlated and therefore their effects on each other need to be taken into consideration. Thus, the following Pearson correlation coefficients were used in error propagation: (1) the correlation between CBF and the natural logarithm of the BOLD signal is given by the parameter p1, where p1 = ρ(f, ln(s)); (2) the correlation between the numerator and the denominator of the BOLD equation is given by the parameter p2, where
(3) the correlation between CBF and the BOLD signal adjusted derivative is given by the parameter p3, where
Using these coefficients of Equation 4, the SD of CMRO2 (σm) can be described by the following:
where
Given observed values of σf, σv, and σs, it can be shown from Equation 5 that σm is affected by measurement uncertainties of CBF much more than CBV or BOLD signal. Using steady-state measurements of these signals, we estimate that the contribution of σv is five to seven times lower than σf for estimating σm.
By a convolution analysis similar to that described above, a transfer function for CMRO2 was also obtained. However, parameterization of this transfer function was not a gamma variate. Instead, the CMRO2 impulse response function was better fitted by a constrained nonlinear half-order rational equation of order 2/2.5.
Results
Measured neural responses (LFP, MUA)
Electrophysiological recordings of LFP and MUA from the contralateral somatosensory cortex are shown in Figures 1 and 2, respectively. Forepaw stimulation, consisting of square pulses with 2 mA amplitude and 0.3 ms duration, was varied in stimulus number (1–4) and frequency (1.5–6 Hz). Trends of data averaged across six subjects (Figs. 1A, 2A) were similar to individual subject/trial data (Figs. 1B, 2B). Both LFP and MUA data were represented by running 20 ms bins, but because of the lower contrast-to-noise ratio of MUA data, the root mean square (RMS) approach was applied (Stark and Abeles, 2007). In each case, the evoked neural responses to multiple stimuli were normalized to the response of the first stimulus pulse.
Neural responses to varying transient stimuli. A, LFP recorded from the contralateral somatosensory cortex during forepaw stimulation (2 mA, 0.3 ms) by varying stimulus number (1–4) and frequency (1.5–6 Hz). Each trace is 4 s in duration. Responses to subsequent stimuli were normalized to the response of the first stimulus pulse. For stimulus frequency >3 Hz, the alternate responses to subsequent stimuli were significantly attenuated compared with the response of the first stimulus pulse. For stimulus frequency of 6 Hz, the responses were nearly identical for different numbers of stimulus pulses. Similar trends were observed with MUA data (Fig. 2). The short vertical bars on the time axis represent the stimuli, and the long vertical bars on the right side represent ±SD from the mean. B, Individual data with three-pulse stimulation.
Neural responses to varying transient stimuli. A, MUA recorded from the contralateral somatosensory cortex during forepaw stimulation (2 mA, 0.3 ms) by varying stimulus number (1–4) and frequency (1.5–6 Hz). Each trace is 4 s in duration. Responses to subsequent stimuli were normalized to the response of the first stimulus pulse. For stimulus frequency >3 Hz, the alternate responses to subsequent stimuli were significantly attenuated compared with the response of the first stimulus pulse. For stimulus frequency of 6 Hz, the responses were nearly identical for different numbers of stimulus pulses. Similar trends were observed with LFP data (Fig. 1). The short vertical bars on the time axis represent the stimuli, and the long vertical bars on the right side represent the ±SD from the mean. We applied a binned RMS approach (Stark and Abeles, 2007) to process the low signal-to-noise ratio MUA data. Similar results were generated when bins of different sizes (5, 50, 500 Hz) were used with the RMS approach. B, Individual data with three-pulse stimulation.
At the lowest stimulation frequency of 1.5 Hz, each evoked response for individual stimulus pulses was quite similar to the previous ones, with minimal response attenuation for additional stimulus pulses. In contrast, at the highest stimulation frequency of 6 Hz, there were small or negligible evoked responses detected for multiple stimulus pulses. For the intermediate stimulus frequency of 3 Hz, alternate evoked responses to subsequent stimuli were significantly attenuated compared with the response of the previous stimulus pulse (p < 0.05). Stimulus frequency-dependent variations in neural activity are consistent with previous observations using α-chloralose (Matsuura and Kanno, 2001; Sheth et al., 2004), whereas similar relationships have been observed with variation of stimulus duration using other anesthetics (Nemoto et al., 2004; Ureshi et al., 2004). In summary, at stimulus frequency of 6 Hz increasing the stimulus number changed the evoked neural response very little, whereas at the lower stimulus frequencies increasing the stimulus number augmented the evoked neural responses in a linear manner.
Measured hyperemic components (BOLD, CBV, CBF)
BOLD responses from the contralateral somatosensory cortex are shown in Figure 3. Data averaged across eight subjects (Fig. 3A) were quite similar to individual trial data (Fig. 3B). BOLD responses were measured with respect to the prestimulus baseline. Generally, both the intensity and width of each BOLD response increased with increasing stimulus number. For stimulus frequencies of 1.5 and 3 Hz, the BOLD response increased linearly with increasing number of stimulus pulses, whereas for the same number of stimulus pulses the BOLD response decreased with higher stimulus frequency. Enhancement of the BOLD response with increase in stimulus number is correlated with the increase in neural activity (Figs. 1, 2). For stimulus frequency of 6 Hz, the BOLD responses were nearly identical for different number of stimulus pulses, again similar to the neural data (Figs. 1, 2). The mean time-to-peak in BOLD responses was 3.9 ± 0.3 s. BOLD responses for four stimulus pulses was largest at stimulus frequency of 1.5 Hz (8.0 ± 1.3%) compared with 3 Hz (5.4 ± 3.2%) and 6 Hz (3.9 ± 3.2%). This is attributable to the stimulus frequency tuning responses under α-chloralose anesthesia (Sanganahalli et al., 2008). Early human studies found BOLD responses to be mostly linear with respect to variations in stimulus duration (Boynton et al., 1996), whereas recent studies have found some nonlinear tendencies (Pfeuffer et al., 2003; Soltysik et al., 2004). Chloralose-anesthetized animal studies investigating variation of BOLD response with stimulus frequency are in general agreement with our findings (Brinker et al., 1999; Van Camp et al., 2006).
Dynamics of BOLD signal with varying transient stimuli. A, Time courses of BOLD responses from the contralateral somatosensory cortex during forepaw stimulation (2 mA, 0.3 ms) by varying stimulus number (1–4) and frequency (1.5–6 Hz). Each trace is 40 s in duration. For stimulus frequencies of 1.5 and 3 Hz, the BOLD response increased linearly with higher number of stimulus pulses, whereas for the same number of stimulus pulses the BOLD response decreased with higher stimulus frequency. For stimulus frequency of 6 Hz, the BOLD responses were nearly identical for different numbers of stimulus pulses. The BOLD data were obtained at 11.7 T using gradient echo contrast with a temporal resolution of 1 s. The signals were represented as the fractional change from the prestimulus baseline. The short vertical bars on the time axis represent the stimuli, and the long vertical bars represent the ±SD from the mean. B, Individual data with three-pulse stimulation.
CBV-weighted fMRI, using long half-life and plasma-borne intravascular contrast agents, has become an attractive alternative to BOLD-weighted fMRI in animal models (Lu et al., 2007; Kida et al., 2007). We obtained the CBV data in the same subjects used to collect the BOLD data. Figure 4 shows the CBV responses from the contralateral somatosensory cortex. The CBV response was just within the threshold of detection for one stimulus pulse. Previously, it has been shown that the CBV contrast can be enhanced with a higher concentration of the contrast agent (Lu et al., 2007; Englot et al., 2008). However, higher contrast agent concentration would lower the signal-to-noise ratio. Similar to the BOLD results, both the intensity and width of the CBV response were augmented with increase in stimulus number. Averaged data from eight subjects (Fig. 4A) showed similar trends as individual subject data (Fig. 4B). The mean time-to-peak in CBV responses was 3.3 ± 0.7 s. CBV responses for four stimulus pulses were largest at stimulus frequencies of 1.5 Hz (10.9 ± 6.5%) and 3 Hz (9.5 ± 3.2%) compared with 6 Hz (3.5 ± 2.0%). For stimulus frequencies of 1.5 and 3 Hz, the CBV responses increased linearly with increasing number of stimulus pulses, whereas for the same number of stimulus pulses the CBV response decreased with higher stimulus frequency. For stimulus frequency of 6 Hz, the CBV responses were nearly identical for multiple numbers of stimulus pulses, as assessed by experimental error. Optical studies using combination of oxyhemoglobin and deoxyhemoglobin signals to estimate CBV changes demonstrate similar dependencies with variations in stimulus frequency and duration (Sheth et al., 2003; Nemoto et al., 2004).
Dynamics of CBV with varying transient stimuli. A, Time courses of CBV responses from the contralateral somatosensory cortex during forepaw stimulation (2 mA, 0.3 ms) by varying stimulus number (1–4) and frequency (1.5–6 Hz). Each trace is 40 s in duration. For stimulus frequencies of 1.5 and 3 Hz, the CBV response increased linearly with higher number of stimulus pulses, whereas for the same number of stimulus pulses the CBV response decreased with higher stimulus frequency. The CBV data were obtained at 11.7 T using gradient echo contrast in conjunction with Combidex as the contrast agent. The temporal resolution was 1 s. Note that, for single stimulus pulse, the CBV response was negligible. The signals were represented as the fractional change from the prestimulus baseline. The short vertical bars on the time axis represent the stimuli, and the long vertical bars represent the ±SD from the mean. B, Individual data with three-pulse stimulation.
Time courses of CBF responses from the contralateral somatosensory cortex are shown in Figure 5. Generally, the characteristics of CBF responses (intensity, width) were similar to the BOLD and CBV data, but with some caveats. For stimulus frequencies of 1.5 and 3 Hz, the CBF response increased linearly with higher number of stimulus pulses, whereas for the same number of stimulus pulses the CBF response decreased with higher stimulus frequency. For stimulus frequency of 6 Hz, the CBF responses were nearly identical for different number of stimulus pulses, as evaluated by experimental inaccuracies. Averaged CBF data from six subjects (Fig. 5A) showed similar trends as observed with individual subject/trial data (Fig. 5B). The mean time-to-peak in CBF responses was 3.2 ± 0.2 s. The CBF responses for four stimulus pulses were largest at stimulus frequency of 1.5 Hz (100 ± 34%) and 3 Hz (101 ± 29%) compared with 6 Hz (53 ± 22%). Our observation of CBF trends with variations in stimulus duration and frequency is in agreement with previous observations by several studies (Glover, 1999; Ances et al., 2000; Miller et al., 2001; Norup Nielsen and Lauritzen, 2001; Martindale et al., 2003; Pfeuffer et al., 2003). However, note that the nonlinearity of the CBF response (with stimulus number and frequency) was different from that observed with BOLD or CBV data (see below).
Dynamics of CBF with varying transient stimuli. A, Time courses of CBF responses from the contralateral somatosensory cortex during forepaw stimulation (2 mA, 0.3 ms) by varying stimulus number (1–4) and frequency (1.5–6 Hz). Each trace is 40 s in duration. For stimulus frequencies of 1.5 and 3 Hz, the CBF response increased linearly with higher number of stimulus pulses, whereas for the same number of stimulus pulses the CBF response decreased with higher stimulus frequency. For stimulus frequency of 6 Hz, the CBF responses were nearly identical for different numbers of stimulus pulses. The CBF data were obtained simultaneously with the electrical data using an LDF probe (805 nm). The signals were represented as the fractional change from the prestimulus baseline. The short vertical bars on the time axis represent the stimuli, and the long vertical bars represent the ±SD from the mean. B, Individual data with three-pulse stimulation.
Relationships between stimulus parameters and evoked responses
The idea of the “size” of a signal is crucial for establishing relationships between multimodal signals. Because the signal is a function of varying amplitude over time, a reasonable measurement of the signal strength (or energy) is the area under the curve of the signal. Because all multimodal signals had negligible portions below the prestimulation baseline, negative parts of an evoked signal did not contribute significantly to the signal strength assessment. We rationalized that signal strength would better capture the dynamic nuances of the evoked responses, given that in the case of nonelectrical components both signal intensity and width varied, whereas for neural activity both the number and intensity of evoked signals differed.
To demonstrate the effect of increasing stimulus number and frequency on neural activity, the data were normalized to responses obtained with four pulses for stimulation frequency of 1.5 Hz. Figure 6, A and B, shows linear relationships between stimulus number and neural activity at stimulation frequencies of 1.5 and 3 Hz (r2 > 0.97), whereas at stimulus frequency of 6 Hz the activity changed minimally with increasing number of stimulus pulses. The relationship between MUA and LFP was generally similar (Fig. 6C).
Relationships between multimodal signals and stimulus number or frequency. Relationships between stimulus number or frequency versus neural activity (A–C) and hyperemic components (D–F) are shown. Relationships between LFP (G–I) and MUA (J–L) versus hyperemic components are shown. At stimulus frequency of 6 Hz, the responses changed minimally with increasing number of stimulus pulses, whereas linear trends were observed for stimulus frequencies of 1.5 and 3 Hz.
To demonstrate the effect of increasing stimulus number and frequency on hyperemic components, the multimodal data were normalized (independently in each case) to responses obtained with four pulses for the stimulation frequency of 1.5 Hz. As with the neural data (see above), linear relationships were observed with the BOLD, CBV, CBF data versus stimulus number for stimulus frequencies of 1.5 and 3 Hz (Fig. 6D–F). At stimulus frequency of 6 Hz, the responses changed very little with increasing number of stimulus pulses.
However, in Figure 6D–F, note that the response attenuation shifting from low to high stimulation frequencies was significantly different for each hyperemic component. For example, with BOLD data there was nearly 60% attenuation in the linear trend when shifting from stimulus frequency of 1.5–3 Hz, whereas with CBF that attenuation was only ∼20%. However, with CBV there were nearly two 50% attenuation steps in the linear trends when shifting from stimulus frequency of 1.5–3 Hz versus 3–6 Hz. In contrast, the attenuation in the linear trends for BOLD and CBF were quite different when shifting from stimulus frequency of 1.5–3 Hz versus 3–6 Hz.
We plotted the normalized responses to quantify the relationship between hyperemic components (BOLD, CBV, CBF) and neural activity (Fig. 6G–L). Because the neural and hyperemic signals at stimulus frequency of 6 Hz showed very little change with increasing stimulus number, this analysis was omitted here. Relationship between hyperemic components and LFP (Fig. 6G–I) for stimulus frequencies of 1.5 and 3 Hz showed high linearity (r2 ≥ 0.94). However, analysis with data from either three versus four stimulus pulses did not reveal relationships that significantly favored power law and/or exponential dependencies (data not shown). Similar trends were observed with the MUA data (Fig. 6J–L). These results are generally consistent with previous findings that have used chloralose-anesthetized rats (Sheth et al., 2004; Ureshi et al., 2004; Van Camp et al., 2006).
Linearity and time invariance (universal transfer function)
We used convolution analysis to find a transfer function that connects the neural signal to each hyperemic component, the effectiveness of which was characterized by the residual signal as given by the difference between the measured and predicted signals. If temporal fluctuations of the residual signal were smaller than the uncertainty of the measured signal over time (i.e., SD of the experimental measurement), the convolution process would produce a universal impulse response function (supplemental Fig. 3, available at www.jneurosci.org as supplemental material). For each hyperemic component, the objective was to find a single (i.e., universal) transfer function that could be used to model the measured hyperemic component successfully for all stimulus parameters (i.e., residual “noise” in modeling was smaller than the measurement uncertainty). This was important because that would suggest linearity between the hyperemic component and neural activity and therefore the respective transfer function should be time invariant within the temporal regimen of the stimulus parameters.
For each hyperemic component, we parameterized the transfer function with a gamma variate (for convolution analysis results with LFP and MUA, see supplemental Tables 1 and 2, respectively, available at www.jneurosci.org as supplemental material). A detailed discussion of the transfer functions is beyond the scope here. But briefly, the transfer functions generated by convolution analysis with LFP and MUA were quite similar in shape and the time-to-peak for BOLD, CBV, and CBF impulses were ∼2.6, 2.1, and 1.8 s, respectively. The systematic difference between the time-to-peak of each measured response and each transfer function is consistent with previous observations (Ances et al., 2000; Buracas and Boynton, 2002; Silva et al., 2007). However, it should be noted that previous studies either did not include neural data or used integrated neural data for deconvolution (see Discussion). We used the strength of the neural data in our convolution analysis to specifically target the subtle nonlinearity of neural data over time and their impact on the impulse response function. Because the input signal is extremely short compared with the transfer function, the result is expected to be close to the theoretical impulse response behavior of the system.
To demonstrate linearity and time invariance of the transfer function, we conducted two systematic comparisons between predicted and measured signals. In one comparison, the RMS of the residual signal was compared with the SD of the measurement (for analysis with LFP and MUA, see supplemental Tables 3 and 4, respectively, available at www.jneurosci.org as supplemental material). The RMS of the residual signal was averaged across the entire data set and compared with the average of the measurement SD over the same duration. In each case, the comparison showed that the residual “noise” in modeling was smaller than the measurement uncertainty.
In another comparison, only the first 10 s of data (after stimulation onset), which approximately corresponded to the hyperemic phase, were compared (for analysis with LFP and MUA, see supplemental Fig. 4A,B, respectively, available at www.jneurosci.org as supplemental material). Given the lower temporal resolution of the BOLD and CBV data (one image per second), the CBF data were de-interpolated to a lower temporal resolution for this comparison. Because experimental conditions were varied by three stimulus pulses (2, 3, 4) and three stimulus frequencies (1.5, 3, 6 Hz) as well as the one stimulus pulse data, there were a total of 10 data sets per experimental criteria. Given the temporal resolution of 1 s, the data points shown in each part of supplemental Figure 4 (available at www.jneurosci.org as supplemental material) should be a total of 50. This graphical representation, similar to the RMS versus SD comparison above, shows that the predicted and measured signals were well correlated in each case, thereby suggesting linearity and time invariance of each transfer function generated by the convolution analysis of multimodal data.
CMRO2 transients (dynamic calibrated fMRI)
Detailed analysis of the multimodal data (Figs. 1⇑⇑⇑–5) displayed linear relationships between neural activity and each of the hyperemic components (Fig. 6), despite the fact that there were some nonlinearity trends between each signal and stimulus parameters. Furthermore, a comprehensive convolution analysis (supplemental Fig. 4, Tables 3, 4, available at www.jneurosci.org as supplemental material) revealed that the respective transfer functions (of BOLD, CBV, CBF) were time invariant for extremely brief stimuli, even down to events separated by <200 ms. Collectively, these results secured the possibility of using calibrated fMRI in a dynamic manner for calculation of CMRO2 transients.
Figure 7A shows the CMRO2 time courses obtained from group-averaged data. However, using BOLD, CBF, and CBV data from a single subject, we obtained similar CMRO2 time courses (supplemental Fig. 5A, available at www.jneurosci.org as supplemental material). For stimulus frequencies of 1.5 and 3 Hz, the CMRO2 response (intensity, width) increased linearly with higher number of stimulus pulses, whereas for the same number of stimulus pulses the CMRO2 response decreased with higher stimulus frequency. For stimulus frequency of 6 Hz, the CMRO2 responses were nearly identical for different number of stimulus pulses. The mean time-to-peak in CMRO2 responses was ∼3.5 s, which is slightly delayed compared with CBF. CMRO2 responses for four stimulus pulses were largest at stimulus frequency of 1.5 Hz (63 ± 38%) and 3 Hz (68 ± 39%) compared with 6 Hz (38 ± 27%).
Dynamics of CMRO2 with varying transient stimuli. A, Time courses of CMRO2 responses from the contralateral somatosensory cortex during forepaw stimulation (2 mA, 0.3 ms) by varying stimulus number (1–4) and frequency (1.5–6 Hz). Each trace is 40 s in duration. For stimulus frequencies of 1.5 and 3 Hz, the CMRO2 response increased linearly with higher number of stimulus pulses, whereas for the same number of stimulus pulses the CMRO2 response decreased with higher stimulus frequency. For stimulus frequency of 6 Hz, the CMRO2 responses were nearly identical for different numbers of stimulus pulses. Based on our observation of linearity of evoked responses (Fig. 6) and time invariance of the transfer function (supplemental Fig. 4, Tables 1, 2, available at www.jneurosci.org as supplemental material), the CMRO2 data were obtained by using Equation 3. The signals were represented as the fractional change from the prestimulus baseline. The short vertical bars on the time axis represent the stimuli, and the long vertical bars represent the ±SD from the mean. B, Relationships between CMRO2 versus stimulus number (top), CMRO2 versus LFP and MUA (middle) for 1.5 Hz (△) and 3.0 Hz (○), and CBF versus CBV and CMRO2 for all stimulus conditions (bottom).
Figure 7B shows the relationships between CMRO2 versus stimulus duration (top panel), neural activity (middle panel), and hemodynamic parameters (bottom panel). Existing methods cannot provide high temporal resolution CMRO2 data. Therefore, based on the thermodynamic principle of relating energy expended for the work done (Shulman et al., 2002), we used the neural activity as a pseudo-validation for the predicted CMRO2 transients. Furthermore, given that the BOLD signal is dependent on oxygen transport from blood to tissue (Hyder et al., 1998), we compared the CBF–CBV and CBF–CMRO2 couplings.
Linear relationships were observed for CMRO2 versus stimulus duration (Fig. 7B, top panel) and CMRO2 versus neural activity, with either LFP (Fig. 7B, middle left panel) or MUA (Fig. 7B, middle right panel). Because the multimodal signals at stimulus frequency of 6 Hz showed very little change with increasing stimulus number, this analysis was excluded here. Recent studies in primates have implied that LFP corresponds to the neural input and MUA corresponds to the neural output (Logothetis et al., 2001), which approximates presynaptic and postsynaptic signaling, respectively. As stated previously, LFP and MUA generally behaved similarly as a function of stimulus number within the range of stimulus frequencies tested (Fig. 6A,B), resulting in good LPF–MUA coupling (Fig. 6C). Therefore, the linear relationships depicted in Figures 6C and 7B, middle panel, suggest that the presynaptic and postsynaptic activities are approximately equally weighted for the rat somatosensory cortex, which is in good agreement with theoretical predictions (Attwell and Laughlin, 2001).
The impact of CBV on CMRO2 calculation (Eqs. 3–5) cannot be overlooked because of the dynamic relationship between changes in flow and volume [i.e., CBV = CBFφ, which is the so-called Grubb's law (Grubb et al., 1974)]. The Φ value calculated for the entire hyperemic phase was ∼0.15 (Fig. 7B, bottom left panel), which is in close agreement with recent animal studies (Jin and Kim, 2008; Shen et al., 2008). However, the value of Φ was allowed to vary over time for the dynamic CMRO2 calculation (Kida et al., 2007).
Our observation of CMRO2 behavior with stimulus duration (Fig. 7B, top panel) is similar to that observed with the CBF data (Fig. 6F), thereby demonstrating a tight CBF–CMRO2 dynamic coupling, ranging between 3:2 and 5:2 (Fig. 7B, bottom right panel). Although the characteristics of CBF dynamics appear to be quite similar to CMRO2 transients, careful examination of these two impulse response functions shows that one is fitted by a gamma variate (CBF), whereas the other one is not (CMRO2), as discussed below. The CBF–CMRO2 coupling of ∼2:1 is in agreement with predictions of oxygen transport models (Herman et al., 2006; Huppert et al., 2007) and results from steady-state calibrated fMRI studies (Hoge et al., 1999; Zhu et al., 2005; Chiarelli et al., 2007).
Using a similar convolution analysis as described previously (supplemental Fig. 3, available at www.jneurosci.org as supplemental material), transfer functions for CMRO2 can be obtained from LFP and MUA data. Supplemental Figure 5B (available at www.jneurosci.org as supplemental material) shows the comparison between CMRO2 measured (by using Eq. 3) and predicted by using the transfer functions. The transfer functions were parameterized by a constrained nonlinear half-order rational equation (for impulse response functions with LFP and MUA, see supplemental Tables 1 and 2, respectively, available at www.jneurosci.org as supplemental material). The good agreement between CMRO2 measured (by using Eq. 3) and predicted by the two independent transfer functions provides additional validation of dynamic calibrated fMRI.
These results of CMRO2 transients raise the issue of the early dip observed with recordings of tissue pO2, optical imaging of deoxyhemoglobin, and/or BOLD signal, all suggesting the possibility of an earlier CMRO2 increase (within the first 500–700 ms). The early tissue pO2 dip has been observed in anesthetized rats and cats (Ances et al., 2001; Thompson et al., 2003). Although optical imaging of deoxyhemoglobin has reported the early dip (Malonek and Grinvald, 1996), some inconsistencies remain about its visibility (Lindauer et al., 2001; Mayhew et al., 2001). The early dip in the BOLD signal has been noted primarily in some human studies (Menon et al., 1995; Hu et al., 1997), but it remains unclear whether this observation is specific to humans because no rat (Silva et al., 2000) or monkey (Logothetis et al., 2001) studies have reported the early dip, whereas some cat studies have reported the dip (Kim et al., 2000) and others have not (Jezzard et al., 1997). Furthermore, there is a possibility that the early BOLD signal dip may be sensitive to blood pCO2 (Harel et al., 2002). Despite these factors, the observation of the early dip remains a very important topic of research (Ances, 2004). However, its magnitude is near detection threshold of most techniques. Given the fact our CMRO2 calculation is based on combination of several signals (Eqs. 3–5), it is unlikely whether detection of the early BOLD signal dip (i.e., early CMRO2 increase) will be feasible with current sensitivity. Future fMRI sensitivity improvements may allow the early BOLD signal dip and thereby allow even earlier CMRO2 increases to be measured.
Discussion
Recent studies have reestablished the importance of energetics for brain function (Hyder et al., 2006; Buzsáki et al., 2007; Riera et al., 2008). Because MRI detects nearly the whole brain with decent spatiotemporal resolution, a desirable goal is to use calibrated fMRI to extract CMRO2 dynamics for event-related experiments, as commonly done with block-design paradigms (Hoge et al., 1999; Hyder et al., 2001). However, a limiting factor for estimating transient energetics from calibrated fMRI is whether the conventional BOLD model can be applied transiently. To do this, each hyperemic component has to be coupled to neural activity in such a way that all points of the hyperemic phase would carry the same calibration. If so, then the linear and time-invariant impulse response functions can be used in dynamic calibrated fMRI. The first step was to experimentally assess whether linearity exists between each hyperemic component and neural activity, as assessed by the strength of each signal type. The second step was to test the universality (or time invariance) of the respective transfer functions to accurately predict the measured responses (within uncertainties). Because of the lack of other techniques that can measure transient energetics, neural activity was used to validate the energetics based on a thermodynamic principle (energy demanded ≡ work done). Within limits of experimental error, we identified a component of the BOLD signal for transient CMRO2 changes.
An important consideration for interpretation of our results is practical limitations of the used methods. An obvious caution is consolidating measurements from a single probe (electrical, optical) with high-resolution fMRI data. Because of the type of microelectrodes we used, sampling from larger spatial regions could not be made simultaneously. Therefore, our electrical data were mainly confined to the middle cortical layers. However, preliminary electrophysiology work within our laboratory (data not shown) suggests that the current results may be extended to layers above, but not necessarily below. The multimodal fMRI studies could have been improved by ASL measurements of CBF dynamics, but that would come at a cost of perfusion sensitivity loss at higher temporal resolution (Kida et al., 2004). Because all MRI measurements (of BOLD, CBV, CBF) are made concurrently, the induced responses have to be reproducible (Sanganahalli et al., 2008). Although our results provide temporal transfer functions, the spatial transfer functions would need to consider spatially varying activity patterns.
Linearity of evoked responses
Relationships between neural activity and hyperemic components are an active research area (Friston et al., 1998; Logothetis et al., 2001). Although studies have emphasized either linear or nonlinear tendencies, a consistent finding is that amplitude/contrast variation has significant effects on nonlinearity trends because of their impact on response saturation and/or adaptation. Because nonlinear tendencies were observed between neural and imaging signals when stimulus amplitude was varied with different anesthetics (Devor et al., 2003; Sheth et al., 2004), we rationalized that temporal evolution of neural activity (and energetics) is likely to be impacted by temporal aspects of stimuli (i.e., duration, frequency) as suggested by other studies (Nemoto et al., 2004; Ureshi et al., 2004). Therefore, we used number of pulses and interval between them as stimulus variables.
Although different anesthetized states generate different magnitude of responses (Smith et al., 2002; Maandag et al., 2007; Masamoto et al., 2007), the coupling between changes in neural and imaging signals are quite well preserved across these baseline levels. In support of this, previous rat studies using similar stimulation paradigms but different anesthetics showed that neural activity is coupled with hemodynamic signals [e.g., when stimulus duration (Nemoto et al., 2004; Ureshi et al., 2004) or frequency (Matsuura and Kanno, 2001; Sheth et al., 2003) are varied]. Recently, in primates it has been shown that the neural and imaging signals are correlated, albeit with different magnitudes, in anesthetized and awake subjects (Goense and Logothetis, 2008). Because the transfer functions are dependent on the measured neural responses and each of the measured imaging signals, we expect that the predictions from this study could be applied to slight variations of baseline (Huttunen et al., 2008) with minor modifications of forepaw stimuli (Herman et al., 2009b). However, pending additional studies, caution should be applied for extending the results beyond the somatosensory cortex.
In our study, both LFP and MUA showed linearity with stimulus number at lower frequencies, but not at high frequency at which the alternate evoked potentials or commensurate unit activities were almost fully suppressed. As repetition time between stimulus pulses becomes shorter, neural refractory plays a major role in nonlinear behavior of evoked responses, specifically with high stimulus frequencies, as also observed by others (Sheth et al., 2004; Ureshi et al., 2004). It is likely attributable to synaptic depression and/or other mechanisms, which were not part of current study goals but are being investigated by others (Hellweg et al., 1977; Simons, 1985). Similar to observations of neural activity data, each hyperemic component showed linearity with stimulus number at lower frequencies. At high frequency, however, responses were nearly identical when stimulus number was varied as observed by others (Norup Nielsen and Lauritzen, 2001; Franceschini et al., 2008), which suggests that moment-to-moment neural activity variations can manifest onto each hyperemic component (Roy and Sherrington, 1890).
Convolution analysis
Many studies have applied convolution analysis with the BOLD response to report linearity (Boynton et al., 1996) or nonlinearity (Vazquez and Noll, 1998). Because of the lack of neural activity recordings, past studies used the stimulus itself as the input. Results from such studies, although insightful about convolution theory, should be interpreted with caution. Recent studies include neural activity to reveal specific trends of linearity or nonlinearity that agree with our observations (Matsuura and Kanno, 2001; Sheth et al., 2004). There are, however, far fewer studies that have applied convolution analysis with CBF and CBV (Ances et al., 2001; Martin et al., 2006; Silva et al., 2007). In general, the transfer functions generated for BOLD, CBV, and CBF in the current study are in good agreement with literature reports. Although a detailed discussion of each transfer function is beyond the scope here, general observations are that the CBF impulse peaks the earliest and the BOLD and CBF impulses are different in shape compared with the CBV impulse.
It should be emphasized that some previous studies have used an integration procedure between consecutive evoked neural signals (Ances et al., 2000; Mathiesen et al., 2000), whereas we used the signal strength. Interpolating between separate neural events creates the impression of a more weighted signal. For long stimulus durations, this process may have minimal impact on deconvolution because smoothing changes the signal envelope very little. In contrast, for brief stimuli this approach could generate nonlinearity from deconvolution because smoothing creates a signal that is quite different from the actual response. If the experimental goal is to include intricacies of the neural response from moment to moment, then integration would likely defeat that intent. In this study, treatment of each evoked neural response was a key feature used in both linearity assessments across different measurements and time invariance demonstration of transfer functions. Questions remain, however, whether these transfer functions could be successfully used for much longer trains of stimulus pulses (Martindale et al., 2005; Herman et al., 2009b).
Transient energetics and brain function
Agreement between CMRO2 and neural activity provides a pseudo-validation and suggests that CMRO2 may contribute to BOLD signal with consecutive events that are <200 ms apart. Mitochondrial energy metabolism is strongly linked to both neuronal and glial activities (Riera et al., 2008). During augmented workload, increased energy demand stimulates reduction of oxygen to water and generates reduced cofactors in the tricarboxylic acid cycle (i.e., NADH, FADH2). Therefore, the initial dip in tissue pO2 during activation may indicate the initial demand for oxygen (Ances et al., 2001; Thompson et al., 2003), but its subsequent overshoot induces difficulties in quantifying CMRO2 transients (Caesar et al., 2008). Because it is also well known that rapid Ca2+ influx stimulates mitochondrial electron transport chain (Szabadkai and Duchen, 2008), it is possible that neuronal and glial Ca2+ transients are related to CMRO2 dynamics. A recent study, using two-photon imaging of Ca2+ signals in ferret visual cortex (Schummers et al., 2008), demonstrated neuronal and glial time courses, which together mimics the CMRO2 transients demonstrated here (data not shown). However, more studies are needed to link ion flux to transient energetics in vivo.
Conclusion
Changes in neural activity are resultant of ion permeability variations occurring within the tens to hundreds millisecond level and are aided by a variety of active energy consuming processes. CMRO2 (or ATP synthesis) is a slower process involving mobilization of metabolic intermediates, enzymes, and barriers. Although CMRO2 transients appear slower than the actual neural activity, their coupling for events separated by <200 ms apart demonstrates the sensitivity of high field dynamic calibrated fMRI studies for event-related paradigms.
Footnotes
-
This work was supported by National Institutes of Health Grants R01 MH-067528, R01 DC-003710, and P30 NS-52519 (F.H.). We thank scientists and engineers at the Magnetic Resonance Research Center (http://mrrc.yale.edu) and the Core Center for Quantitative Neuroscience with Magnetic Resonance (http://qnmr.yale.edu). We also thank the anonymous reviewers for their helpful suggestions.
- Correspondence should be addressed to either Basavaraju G. Sanganahalli or Fahmeed Hyder, N143 TAC (Magnetic Resonance Research Center), 300 Cedar Street, Yale University, New Haven, CT 06520. basavaraju.ganganna{at}yale.edu or fahmeed.hyder{at}yale.edu
