|Home | About | Journals | Submit | Contact Us | Français|
We describe the use of the three dimensional characteristics of the functional magnetic resonance imaging (fMRI) blood oxygenation level dependent (BOLD) and cerebral blood volume (CBV) MRI signal changes to refine a two dimensional optical imaging spectroscopy (OIS) algorithm. The cortical depth profiles of the BOLD and CBV changes following neural activation were used to parameterise a 5-layer heterogeneous tissue model used in the Monte Carlo simulations (MCS) of light transport through tissue in the OIS analysis algorithm. To transform the fMRI BOLD and CBV measurements into deoxy-haemoglobin (Hbr) profiles we inverted an MCS of extravascular MR signal attenuation under the assumption that the extra-/intravascular ratio is 2:1 at a magnetic field strength of 3T. The significant improvement in the quantitative accuracy of haemodynamic measurements using the new heterogeneous tissue model over the original homogeneous tissue model OIS algorithm was demonstrated on new concurrent OIS and fMRI data covering a range of stimulus durations.
Optical imaging spectroscopy (OIS) is a technique used in neuro-imaging for the quantification of changes in cerebral blood volume and blood oxygenation during neural activity, with high temporal and spatial resolution (Frostig et al., 1990; Malonek and Grinvald, 1996; Mayhew et al., 2000; Nemoto et al., 1997). Jezzard et al. (1994) combined OIS and functional magnetic resonance imaging (fMRI) techniques to investigate the blood oxygenation level dependent (BOLD) signal. They showed that NMR signal changes and oxygen saturation (Y) signal changes (measured by spectrophotometric reflectance) were in good temporal agreement for periods of small blood volume changes. They concluded that using concurrent OIS and fMRI studies led to a better understanding of the nature and time course of the BOLD signal change during neural activity. However, OIS techniques rely on many assumptions concerning the light absorption and scattering properties of the imaged tissue, while MRI offers a more fundamental measurement of MR signal attenuation. Hence, in this current work we take the opposite approach and use MRI measurements of BOLD and cerebral blood volume (CBV) changes, as well as CBV baseline measurements, to refine the optical imaging spectroscopic analysis.
In our previous work concerning optical imaging spectroscopy (Jones et al., 2001; Kennerley et al., 2005a; Martin et al., 2002; Mayhew et al., 2000), both the baseline total haemoglobin and oxygen saturation values (Hbt0 and Y0 respectively), and the relevant changes (ΔHbt and ΔY) were assumed to be spatially constant throughout the brain tissue. This homogeneous tissue model approach is standard in optical imaging experiments (Frostig et al., 1990; Malonek and Grinvald, 1996; Mayhew et al., 2000; Nemoto et al., 1997). However it is well known that in cortical brain tissue there is a heterogeneous distribution of vascular structures (Pawlik et al., 1981) suggesting a layered tissue model may be more appropriate. This is supported by various MRI measurements (Jin and Kim, 2008; Lu et al., 2004; Silva and Koretsky, 2002; Silva et al., 2000; Zhao et al., 2006; Zhao et al., 2004), in particular, BOLD and CBV changes (Kennerley et al., 2005a; Mandeville and Marota, 1999), and baseline blood volume (Kennerley et al., 2005a; Tropres et al., 2001) which have different spatial profiles down the cortical lamina (Figure 1). The novelty of the present work was to use fMRI data, from Kennerley et al. (2005a), to parameterise a heterogeneous tissue model for use in OIS.
We implemented the layered heterogeneous tissue model and measured light attenuation captured using interleaved illumination at four wavelengths in the visible region. Spectroscopic analysis produced a 2D spatial map and accompanying pixel by pixel based time series of the changes in blood volume (Hbt) and oxygenation (Y). The significant improvement in performance of the modified OIS algorithm in comparison with the original OIS algorithm assuming homogeneous tissue was demonstrated on new concurrent OIS and fMRI data covering a range of stimulus durations. We found that spectroscopic analysis using the proposed heterogeneous tissue model resulted in time series of changes in oxy- and deoxy-haemoglobin that when used as inputs to a biophysical model of MR signal attenuation, gave a significantly better prediction of concurrently measured BOLD signal changes.
In this work we use the BOLD signal and other relevant fMRI measurements to parameterise optical imaging spectroscopy algorithms. Preliminary work and results were reported in abstract form by Kennerley et al. (2005b; 2006) and Martindale et al. (2005).
OIS exploits the differences of absorption by chromophores as a function of wavelength, λ. In the simple case of light passing through a non-scattering medium in a laboratory cuvette with thickness L and having no reflections at the surface, the Beer-Lambert law states that the attenuation of light, A(λ), is linearly related to the concentrations, ci, of the absorbers in the medium by:
where εi(λ) is the specific extinction coefficient of the ith absorbent. However this linear relationship no longer holds for a turbid medium, such as brain tissue. The path length of a photon in such a scattering medium is likely to be longer than that in a non-scattering medium, and hence the probability of the photon being absorbed is increased. The relationship between light attenuation A(λ)and the absorbent concentration ci becomes nonlinear and dependent on the absorption properties of the medium.
A more tangible estimate of the path length often used in OIS algorithms is the differential path length, Ld(λ), which relates changes in attenuation, ΔA, to changes in the concentrations of absorbents, Δc, by:
The differential path length thus defined is related to the mean and partial path lengths, as explained below.
Changes in neuronal activity produce changes in the concentration of the chromophores oxy- and deoxy-haemoglobin (HbO2 and Hbr); the effects of changes in other chromophores (e.g. water, cytochrome & lipids) are generally assumed to be negligible in the visible wavelengths. Thus Equation 2 becomes:
We implemented an OIS algorithm to estimate changes in oxy- and deoxy-haemoglobin concentration from changes in remitted light attenuation measured at several wavelengths. At the core of our OIS algorithm is an MCS of light transport through brain tissue (Schmitt et al., 1994) which estimates the unknown path length, Ld, as a function of wavelength. From diffusion theory, the path length is dependent only on the absorption coefficient μa(λ)(εHbO2(λ)HbO2 + εHbr(λ)Hbr) and the reduced scattering coefficient of the tissue. It is well known that the total hemoglobin concentration (Hbt) and the blood oxygen saturation (Y) can be used interchangeably with HbO2 and Hbr. Therefore the MCS of light transport through tissue can be parameterised with values for baseline Hbt0 and Y0. Once the path length is estimated, the changes in Hbt and Y due to evoked changes in neural activity can be calculated using Equation 3.
Previously published optical imaging spectroscopy studies (Frostig et al., 1990; Jones et al., 2001; Kennerley et al., 2005a; Malonek and Grinvald, 1996; Martin et al., 2002; Mayhew et al., 2000; Nemoto et al., 1997) have adopted a homogeneous tissue approach in data analysis, i.e. assumed that both the baseline values Hbt0 and Y0 and the changes ΔHbt and ΔY were spatially constant throughout the brain tissue.
Here we formed a layered tissue model, with boundaries at 100, 300, 700, 1500 and 10000µm, for the MCS of light transport, where the partial path lengths (Hiraoka et al., 1993) were estimated for each layer (see supplementary material). The optical properties are encapsulated by the absorption coefficient, μa, scattering coefficient, μs and anisotropy, g, for each layer. The in vivo measurements of Johns et al. (2005) indicate that the reduced scattering coefficient is largely homogeneous over the superficial 2mm of tissue that the optical system is most sensitive to, hence we have assumed that the scattering property of the brain tissue is homogenous and unchanged during neural activation. Based on the work by van der Zee et al. (1992) the reduced scattering coefficient was set to . Due to the paucity of precise empirical information the anisotropy was set to g=0.85 for all layers, hence the model was homogeneous in scattering.
The tissue was therefore characterized by the parameters (μa1, …, μan; μ′s) where n is the number of layers. The memory requirements increased exponentially with the number of parameters, limiting us to five layers. The choice of layers boundaries was therefore a compromise between computer power and delineating the tissue based upon its physical and physiological significance (Hutsler et al., 2005; Kotter et al., 1998). The top three layers are each responsible for about 1/3 of the optical signal, with the first layer including the surface vessels and the third layer containing cortical lamina IV. The fourth layer was determined by the thickness of the cortex. The large depth of the deepest layer reflects the fact that using visible wavelengths, very little remitted illumination penetrates below 1500µm (Kennerley et al., 2005a), and is essentially a ‘catch all’ boundary condition.
The new MCS of light transport through tissue was parameterised with cortical depth profiles of i) baseline total haemoglobin, Hbt0, and ii) baseline oxygen saturation, Y0. The unknown mean path length, Ld(λ), used in the analysis is simply the weighted sum of the partial path lengths. The layer weight is determined using depth profiles of iii) normalised changes in total haemoglobin, and iv) changes in oxygen saturation. Below we detail how MRI data taken from Kennerley et al. (2005a) was used to fashion these profiles for use in optical imaging spectroscopy.
By measuring the transverse relaxation rate change induced by a highly paramagnetic contrast agent (MION) and replicating the methods of Tropres et al. (2001), a cortical depth profile of baseline blood volume fraction was measured (Figure 1b). Smoothing and partial volume effects blur the brain boundary around 0mm. The profile of blood volume fraction, νf, can be converted into baseline haemoglobin (Hbt0) concentration (Figure 1b) in brain tissue by
where Hct is the fraction of haemocrit (0.46), Rc/l is the ratio of cerebral tissue to large vessels (0.69) (Wyatt et al., 1990), [Hb]RBC the concentration of haemoglobin in red blood cells (340g/l) and MmHb the molecular mass of haemoglobin (64450g/M). Following conversion we found that in the superficial areas of cortex the baseline blood volume was 104µM decreasing to 60µM in the deeper layers (Kennerley et al., 2005a). The profile was quantised between the tissue model boundaries and normalised between 0 and 1 where 1 indicated the largest baseline volume (i.e. 104µM in the superficial layers).
Although BOLD fMRI signal changes primarily reflect oxygenation changes we currently have no MRI method to allow direct measurement of baseline oxygen saturation. Therefore, like our previous work, baseline oxygen saturation was assumed to be 50% (Gray and Steadman, 1964; Vovenko, 1999) and homogenous through the cortex. Implications of this are covered in the discussion section of this manuscript.
In our previous work (Kennerley et al., 2005a) we directly compared MRI measurements of CBV changes (plasma volume measured with MRI following contrast agent injection) and concurrent OIS measurements of Hbt changes. It was shown that the OIS measurements were predominantly from the initial 500µm of cortex. In these superficial cortical layers the magnitude of CBV and Hbt changes were found to be identical. Assuming a constant haematocrit level (Sheth et al., 2004), CBV and Hbt changes can be used interchangeably and the profile of the normalised haemoglobin changes is identical to the profile of measured CBV (plasma volume measured with MRI following contrast agent injection) changes:
From this the profile of the normalised changes in total haemoglobin can be obtained.
Kennerley et al. (2005a) showed that optical measurements of changes in Hbt were equal to CBV measurements (from MION-enhanced MRI) in the superficial cortex. Therefore, we normalised this profile to 1 in the superficial layers (0–700µm) increasing with depth as the CBV profile increased. In the final layer (1500–10000µm), outside the measured CBV profile, Hbt changes were assumed to be zero, as this depth is outside the cortex and less than 1% of visible light penetrates to this depth.
The relationship between oxygen saturation and deoxy-haemoglobin concentration is given by:
Thus the formation of a profile of changes in oxygen saturation for use in the heterogeneous tissue model requires an estimation of normalised changes in deoxy-haemoglobin Hbr/Hbr0. MRI cannot be used directly to provide such a measure. However, by inverting a Monte Carlo simulation look-up table (MCS-LuT) of MR signal attenuation (Martindale et al., 2008) it is possible to estimate this parameter using the measured profiles of the normalised BOLD and CBV changes (Figure 1c).
It is important to note that although simulations were run at baseline blood volume of 4%, Martindale et al. (2008) shows a simple relationship exists to scale results to any baseline blood volume. This is important here because the baseline blood volume fraction in rat cortex has been measured as 6% (Kennerley et al., 2005a). The use of the scaling relationship was justified by running simulations again for 100 vessel networks at different baseline blood volume fractions (data not shown).
The MCS-LuT of MR signal attenuation (Martindale et al., 2008) takes as an input measured changes in Hbr and Hbt and produces an estimate of the BOLD signal changes. We inverted the MCS-LuT and used the BOLD and CBV depth profiles to estimate the Hbr changes as a function of cortical depth, under the assumption that CBV and Hbt can be used interchangeably. For the purpose of the MCS-LuT we assumed an extra-vascular (EV) compartment with a mean vessel radius of 10µm. We consider only the EV compartment of the MR signal as we do not know the value of the intrinsic extra-/intravascular ratio (Martindale et al., 2008).
It has been shown that at 3T the EV compartment accounts for approximately 67% of the total BOLD signal changes (Lu and van Zijl, 2005) (this has been confirmed in studies using bi-polar crushing gradients, see Appendix A). Hence we take 67% of our measured total BOLD signal changes as an estimate of the EV space contribution. We used this estimation of EV BOLD signal changes along with profiles of blood volume through the cortex in an inverted MCS-LuT of MR signal attenuation, resulting in a profile of Hbr changes (Figure 1c). For illustrative purposes, the spatial distribution of the predicted Hbr profile was compared to a photomicrograph of a coronal histological section stained for the metabolic enzyme cytochrome oxidase. The predicted Hbr changes peak around layer IV (~500µm) of the somatosensory cortex. The section was taken from Jones et al. (2004) who used multi-channel electrodes to determine neuronal activity as a function of cortical depth following whisker stimulation. Interestingly the current source density analysis showed a peak in sink activity centred also on layer IV.
Using the predicted Hbr profile from the inverted MCS-LuT of MR signal attenuation along with the CBV profile, the saturation profile was calculated using Equation 6. Note that the profile of baseline blood volume fraction was also used to normalise the CBV changes in terms of micro-molar changes. It was sub-sampled between the boundaries and normalised between 0 and 1. The value 1 indicated the largest change in oxygen saturation.
It should be noted that this estimate relies on the assumptions of constant haematocrit level (Sheth et al., 2004) and the homogeneity of the baseline oxygen saturation.
Using the above we implemented the layered heterogeneous tissue model shown in Figure 2 and evaluated it on new data obtained using concurrent optical imaging and fMRI. The improvement in results and validation of the model are described below.
All aspects of these methods and their development were performed with UK Home Office approval under the Animals (Scientific Procedures) Act 1986. Female Hooded Lister rats weighing between 250g and 400g were kept in a 12-hour dark/light cycle environment at a temperature of 22° C with food and water ad libitum. Animals were anaesthetised with urethane. Rectal temperature was maintained at 37°C throughout surgical and experimental procedures using a homeothermic blanket (Harvard Apparatus). The left and right femoral veins and arteries were cannulated which allowed drug infusion and measurement of mean arterial blood pressure respectively. The skull overlying the right somatosensory cortex was thinned to translucency to allow direct illumination of the cortex and measurement of remitted light.
We followed Jezzard et al. (1994) and collected haemodynamic and BOLD changes following activation using concurrent OIS and fMRI. Experimental methods (including schematics of apparatus developed) for the concurrent acquisition of fMRI and OIS data have been previously described (Kennerley et al., 2005a). In brief, an RF surface coil (20mm diameter) was fixed to the animal’s cranium, ensuring that the thin window lay in the centre of the well. OIS image data capture within the bore of a magnet utilised a modified non-magnetic endoscope (Endoscan Ltd, London) with a 6mm diameter circular window. Non-magnetic electrodes were inserted into the left whisker pad in a posterior direction, ensuring the whole whisker pad was activated during electrical stimulation.
Once the animal was secured within the magnet-compatible holding capsule the endoscope was attached to the well. Inside the capsule, an electrically filtered and isolated heating blanket (Harvard Apparatus) and rectal probe, maintained the temperature of the animal at normal physiological values. The animal was artificially ventilated and blood pressure monitored throughout. 12 animals were used in total.
Once prepared, the animal holding capsule was placed into a 3 Tesla magnet (Magnex Scientific Ltd., MRRS console), fitted with a 10cm inner diameter, self-shielded and water-cooled gradient set (Magnex Scientific Ltd. 10 kHz/mm maximum strength per axis with 100µs ramps). Having identified suitable oblique slices, in the same plane as the OIS (Figure 3), high spatial resolution gradient echo scans were performed (256*256 pixels, FOV = 30mm, slice thickness = 2mm, TR/TE=1000/15ms, flip angle=90°, 2 averages).
Functional data were acquired from the oblique slice using a single shot MBEST Gradient Echo - Echo Planar Imaging (GE-EPI) sequence during electrical stimulation (raw data matrix = 32*64 zero padded to 64*64, data sampling interval = 5µs, FOV = 30mm, slice thickness = 2mm, TR/TE=1000/15ms, flip angle 90°, 2 averages). Read-out direction was left-right. A standard phase correction procedure was used to minimise Nyquist ghosting. The BOLD signal was calculated as fractional change normalised by the mean of the one minute preliminary baseline signal.
The cortex was illuminated with a white light source filtered through a four wavelength (575.4±6.9, 559.0±8.2, 495.0±15.3 and 587.0±4.5nm) switching galvanometer system (Lambda DG-4, Sutter Instruments Company). The spectrographic data were recorded at 32Hz with a 2D spatial resolution of 200µm*200µm. Measurements of changes in remitted light attenuation were passed into an OIS algorithm (Equation 3) to estimate changes in Hbr & HbO2 concentration with a temporal resolution of 125ms. The mean path length (Ld in Equation 3) was estimated using an MCS of light transport through a 5 layer heterogeneous brain tissue model. The heterogeneous tissue model was previously parameterised (see theoretical considerations) using the three dimensional characteristics of the functional magnetic resonance imaging (fMRI) blood oxygenation level dependent (BOLD) signal and cerebral blood volume (CBV) MRI signal from Kennerley et al. (2005a).
Haemodynamic changes in the somatosensory cortex were induced by electrical stimulation of the whiskers. The left whisker pad was electrically stimulated at an intensity of 1.2mA, frequency 5Hz for durations of 2, 4, 8 and 16s. Within a 780s trial and following an initial 60 second control period there were 28, 20, 13 and 7 stimuli respectively. The inter-stimulus intervals were set to allow the cortex to return to a baseline state before the next stimulus. A 5Hz stimulation frequency is known to result in the greatest magnitude of haemodynamic response in the anaesthetised preparation. All stimulus control was performed using a CED1401 which was stimulus locked to MR echo acquisition. Concurrent optical imaging spectroscopy data capture was triggered by hand and locked to stimulus output.
Data analysis was conducted using Matlab (The MathWorks Inc.). After GE- EPI to GE structural image registration, all fMRI data were statistically analysed using the general linear model (GLM) statistical parameter mapping (SPM) approach (Friston et al., 1991). The time series for each voxel was regressed against a design matrix of a representative haemodynamic response (gamma function), a ramp and a DC offset. Activation z-scores were calculated on a voxel by voxel basis. These scores were then mapped onto co-registered structural scans. The threshold for classification as an active region was five or more adjacent voxels, with a z-score for each voxel greater than four.
The 5 layer heterogeneous tissue model, parameterised with MR data from Kennerley et al. (2005a), was compared with the homogeneous model and tested against new concurrent fMRI and optical imaging data covering a range of stimulus durations.
Representative topographic section BOLD statistical parameter maps (superimposed onto 256*256 structural scans) for electrical stimulation (1.6mA, 5Hz) of the left whisker pad as a function of stimulus duration (2, 4, 8 and 16s) are shown in Figure 4a. Z-score maps show activation in the contra-lateral somatosensory cortex. As expected the area of activity increases for extended stimulus durations. Concurrent optical imaging functional maps of changes in Hbt, Hbr and HbO2 are compared to BOLD fMRI activation maps (across stimulus duration) in Figures 4b, c and d. For haemodynamic data, spatial maps from both the homogenous tissue model and the new 5 layer heterogeneous tissue model were investigated. The spatial maps of haemodynamic changes (Hbt, Hbr & HbO2) from both the homogenous tissue model and the new 5 layer heterogeneous tissue model were almost identical (see supplementary material) and only the magnitude of the saturation response changed (see below). The Hbt response was, in the most part, arterial and clearly shows the major arterial structure. This is in accord with the findings of Lee et al. (2001) using MRI techniques and Berwick et al. (2005) in a high resolution OIS study. In contrast the changes in Hbr occurred primarily in the parenchyma region corresponding to the barrel cortex.
The time series of fMRI and optical measurements at differing stimulus durations (mean across animals) are shown in Figure 5. Standard errors are given. Using a 2:1 EV/IV ratio (Lu and van Zijl, 2005) the changes in EV BOLD signal time series were estimated. For haemodynamic data, time series estimated from both the homogenous tissue model and the new 5 layer heterogeneous tissue model are shown (assuming 50% oxygen saturation). The Hbt changes obtained from the two tissue models were the same. This was important as Hbt and CBV were shown to be equal in the superficial layers of cortex supporting the assumption of constant haematocrit (see section on theoretical considerations). However, the important difference between the two models is that they give very different results for changes in Hbr. For all stimulus durations, the homogeneous model results in approximately a 5% decrease in Hbr, while the heterogeneous model results in a 10% decrease. This has important implications for the estimates of the BOLD signal which is highly dependent on Hbr changes. The Hbr response to electrical stimulation shows a bi-phasic response and although there is a substantial difference in the maximum decrease in Hbr, the inserts in Figure 5 show that the initial increase in Hbr (Mayhew et al., 2000) is un-affected by the introduction of a heterogeneous tissue model. Further observation reveals the Hbt response peaks before the Hbr reaches its lowest magnitude. Interestingly the response to 2 and 4 second stimuli, peak after stimulus offset indicating a minimum response function peaking just after 4s to any stimulus input. Both the 8 and 16s stimuli response also peak at this time. The response then plateaus until stimulus offset. At this point the response returns to baseline.
Optical data was taken from a region encompassing all aspects of the vasculature (arteries, parenchyma and veins) to make it comparable to the less resolute BOLD fMRI data.
To evaluate the performance of the heterogeneous spectroscopy algorithm we used the optical data across all stimulus durations, analysed using the heterogeneous and homogeneous algorithm in the MCS-LuT of MR signal attenuation, to predict the respective concurrently measured EV BOLD signal changes. We used a 6% blood volume fraction corresponding to the 104µM baseline blood volume used in the spectroscopy analysis. Predictions as a function of mean vessel radius (varied between realistic limits 3–20µm) were compared to EV BOLD; the residual was calculated (data not shown) and minimised to optimise the model fit. The optimised EV BOLD predictions from the haemodynamic data analysed with both the homo- and heterogeneous tissue models are shown in Figure 6. The haemodynamic time series from the heterogeneous spectroscopy algorithm predicts the measured EV BOLD response for all stimulus durations within experimental error. Optimal mean vessel radius for the stimulus durations 2, 4, 8 and 16s were 4.7, 4.8, 6.8 and 7.8µm respectively. These are physiologically plausible vessel radius values. Results reported by Pawlik et al. (1981) measured the mean microvascular radius of cat cerebral cortex using microscopy. This work, often used as a standard reference for microvascular morphology, shows the greatest distribution of vessels to have mean vessel radius of ~3µm with a full width half maximum (FWHM) of 1.5µm. Larger vessels (4–7µm) were also apparent and are likely to be smaller venules and arterioles. This suggests that our responses originated from a compartment containing a mix predominantly of capillaries but also of some larger venules. We hypothesise that although there is very little difference in mean vessel radius between the 2 and 16s data sets (~3µm) any difference could be attributed to larger venules being activated for a longer time period in the latter case. This is partly supported by our finding in Figure 4 showing the regions of activity increasing for increasing stimulus duration.
It is important to note that using the optical data as analysed with the old homogeneous model we do not predict the concurrently measured EV BOLD signal for any given vessel radius within the realistic limits set. For the 16s stimulus the F distribution between the residual of the homo- and heterogeneous tissue models gives an f-score of 2.33. This relates to a probability of p=0.9974 with similar results across the other durations. There is a clear improvement over the homogeneous analysis.
We used the 3D characteristics of fMRI to construct a heterogeneous 5 layer tissue model for use in the analysis of optical imaging spectroscopy data. The improvement over the usual homogeneous tissue model was explored. The spatial maps of haemodynamic changes were unaffected by the choice of tissue model. However, the time series of Hbr changes to neuronal activity were markedly different, although in both cases there was still evidence for the ‘deoxy-dip’ (a small increase in Hbr immediately following stimulus onset and lasting about 2s). We show that with a heterogeneous tissue model, spectroscopic measurements of blood oxygenation are significantly better predictors of a concurrently measured BOLD fMRI signal change to electrical stimulation of the rat whiskers. In order to parameterise the heterogeneous tissue model we made the following simplifying assumptions; i) we used a previous finding that normalised changes in CBV (plasma) equate to changes in Hbt (Kennerley et al., 2005a) suggesting a constant haematocrit level, ii) we used a baseline oxygen saturation of 50%, iii) as the model is only spatially heterogeneous with respect to the absorption coefficient (μa), we assumed a constant scattering coefficient (μs). The justification of these assumptions follows. In a final section we discuss the reasons and implications of finding that the initial increase in the Hbr response immediately following stimulation is unaffected by the introduction of a heterogeneous tissue model. Furthermore, the absence of a corresponding ‘initial dip’ in the BOLD data is discussed. Finer technical details of the extra/intravascular components of the measured BOLD signal and possible BOLD contamination of CBV measurements in Kennerley et al. (2005a) are covered in appendices A and B respectively.
In order to relate Hbt and CBV signal changes, we assumed a constant haematocrit level of 0.46. Thus the 3D characteristics of the fMRI CBV response can be used in the layered tissue OIS model. Although there are some reports (Keyeux et al., 1995; Pawlik et al., 1981) suggesting changes in Hct do occur, Bereczki et al. (1993) found a 7% hypercapnic challenge produced little alteration in microvascular haematocrit. In more recent work in rat barrel cortex following whisker stimulation Kleinfeld et al. (1998) and Sheth et al. (2004) found no evidence of significant changes. Hence we feel secure in assuming a constant haematocrit.
Throughout we assumed the mean baseline oxygen saturation (Y) was 50% for parenchyma regions. Berwick et al. (Berwick et al., 2005) showed that the optically measured Hbr response to neural activation depends importantly on the choice of baseline blood saturation (figure 7a). In comparison with results from human studies, where baseline oxygen saturation in the veins is measured to be 55% (Haacke et al., 1997), the mean value for rat parenchyma (50%) used in our work may seem low. Vovenko (1999) measured the mean oxygen tension PO2 for a region of parenchyma as 40mm Hg which using the oxyhaemoglobin dissociation curve of rat blood (Gray and Steadman, 1964) corresponds to a saturation of ~50%. Gray and Steadman (1964) showed the dissociation curve for human blood is steeper than that of rat blood and hence at a similar oxygen tension (40mm Hg) human blood oxygenation is 70%, which may in part explain the difference between measured saturation values in human and rat.
We investigated what differences in the model of predicted BOLD signal changes (figure 7b) would occur if the baseline oxygen saturation was varied between 30 to 90% in the spectroscopic analysis (figure 7a). It can be seen that the peak BOLD signal occurs with a baseline saturation of 50% with little change in a range of saturation values 30–70% (Figure 7c). Hence we feel a baseline oxygen saturation of 50% is appropriate.
The mean path length is dependent on the absorption coefficient μa(λ) (εHbO2(λ)HbO2 + εHbr(λ)Hbr) and the reduced scattering coefficient of the tissue. In this work we assumed that the scattering property of the brain tissue is homogenous and unchanged during neural activation. Johns et al. (2005) showed that the reduced scattering co-efficient is only constant for the first 2000µm of brain tissue and there is an approximate three times increase in the co-efficient, peaking at around 3500µm. In our MCS of light transport in tissue we assume a lower boundary of 10mm. We tested the robustness of our results by modifying our heterogeneous tissue model to include this effect. The scattering coefficient was held constant, (van der Zee, 1992), in the first four layers (0–1500µm) and increased to in the deepest layer (1500–10000µm). There was negligible effect on the results (data not shown). As previously stated, in the visible wavelengths very few remitted photons have plumbed these depths.
Malonek & Grinvald (1996) have shown that spectroscopic data analysed conventionally with a homogeneous tissue model has a biphasic response in the changes in deoxy-haemoglobin following stimulation. This was confirmed (Mayhew et al., 2000) using an improved spectroscopy algorithm (albeit still using a homogeneous tissue model) which took into account the differential path lengths for the different illumination wavelengths. These were obtained from a Monte Carlo simulation of light transport through tissue (Schmitt et al., 1994). This initial increase in Hbr is believed to be the cause of an observed ‘deoxy dip’ in time series of BOLD signal changes. Lindauer et al. (2001) argued that when the wavelength dependence of photon path length in the brain is taken into account, the initial rise in Hbr vanishes. However we have been unable to replicate their results.
In the current study, Monte Carlo simulations were used to estimate the differential path lengths under the assumption of a heterogeneous layered model (parameterised with 3D MRI data). Our Monte Carlo simulations of light transport through tissue for both a homogeneous and heterogeneous tissue model independently preserve the ‘deoxy dip’ (figure 5). The results from both models can be used to obtain a matrix that transforms results from the spectroscopic analysis from one to the other. In our case this transformation preserved the ‘deoxy dip’ while increasing the magnitude of the Hbr washout; it is possible that other transformations exist that do not preserve the dip. Other than this, we can offer no explanation for the finding of Lindauer et.al. (2001).
In summary we believe the ‘deoxy dip’ to be a reliable phenomenon. Since Ernst et al. (1994) and Hennig et al. (1994), numerous other fMRI studies with sufficient temporal resolution, have shown evidence for an initial dip in the time course of BOLD signal changes in response to neuronal activity (Hu et al., 1997; McIntosh et al., 1996; Menon et al., 1995; Yacoub and Hu, 1999). In this study, the TR was 2s which does not provide sufficient temporal resolution to reveal the ‘deoxy dip’, which itself only lasts 2 seconds (Figure 5). The predicted magnitude of the BOLD ‘deoxy dip’ from these haemodynamic measurements is shown in Figure 6. As can be seen, the magnitude of the modelled dip is less than 0.1%, which is below the statistical error in our BOLD data at field strength of 3T.
The cortical depth profiles of the BOLD and CBV changes following neural activation were used to parameterise a 5-layer heterogeneous tissue model used in the Monte Carlo simulations (MCS) of light transport through tissue. The heterogeneous optical imaging spectroscopy algorithm together with Monte Carlo simulations of MR signal attenuation gives consistent prediction of the concurrently measured BOLD signal. In contrast the homogeneous model under-estimates BOLD changes by approximately 40%. Using the three dimensional characteristics of the functional magnetic resonance imaging signal we have improved and refined our optical imaging spectroscopy analysis. This work will be exploited in our current research using MRI to parameterise a diffuse optical tomography algorithm in a project investigating the metabolic monitoring of cortical activity.
Total BOLD prediction. The Monte Carlo model of MR signal attenuation allows predictions of the total BOLD signal from the haemodynamics if one knows the intrinsic EV/IV ratio (Martindale et al., 2008). Obata et al. (2004) used a value of 1.4 for ε, measured at 1.5T. Here we compared the total BOLD prediction, assuming this ratio holds at 3T, to the measured BOLD signal. There is a slight under-estimation of the total BOLD signal, however the prediction was within experimental error of the measured data. Parameterisation of the intrinsic EV/IV ratio at 3T and TE=15ms may improve the prediction.
The effect of increased diffusion weighting (b, s/mm2) on the BOLD response. (A2a) BOLD signal peak magnitude as a function of diffusion weighting in response to 10% hypercapnia. As diffusion weighting increases the BOLD magnitude decreases. The response size plateaus for large diffusion weighting values where there is a 1/3 decrease in BOLD signal. The majority of IV effects have been ‘crushed’ or nulled. This corresponds to ‘capillary-like’ velocities. Responses to 16s and 2s 1.2mA, 5Hz electrical stimulation of the whisker pad are shown (A2b & A2c respectively) with and without diffusion weighting. Due to the small magnitude of the 2s response (1%), analysis is done in the frequency domain. All studies reveal a consistent 1/3 decrease in signal giving an EV:IV ratio of 2:1. Standard errors are shown.
We gratefully acknowledge support from the National Institutes of Health, Medical Research Council and the Wellcome Trust. We are grateful to Guerbet Pharmaceutical Technologies for providing the MRI contrast agent Sinerem. We thank the laboratory technical staff (Marion Simkins, Natalie Kennerley and Malcolm Benn) for their assistance. Additional thanks to Dr Nikos Papadakis for implementing the IVIM crushing sequence.
Our Monte Carlo simulations of MR signal attenuation result in estimates of BOLD signal changes in both the extravascular (EV) and intravascular (IV) compartments (Martindale et al., 2008). To make a comparison to fMRI measured total BOLD signal changes one must either know i) the intrinsic intra- and extravascular ratio, ε, for use in the model, or ii) directly know the percentage of the total signal changes which originated in the EV space. This value can then be applied to the current measured total BOLD data to give an estimation of the EV signal changes.
Using our MCS-LuT of MR signal attenuation we have been able to accurately predict the estimated EV BOLD signal changes from concurrently measured haemodynamic measurements (using the refined spectroscopy algorithm – see Figure 6). As discussed, our MCS-LuT has the advantage of separating the EV and IV spaces. The summation of the compartments (Equation 11 of Martindale et al., 2008) would allow direct comparison to the measured total BOLD fMRI signal. However, the analysis requires an estimation of the intrinsic intra- and extravascular ratio,ε; we follow Obata et al. (2004) who used a value of 1.4. The total BOLD prediction for the 16s duration stimulation data is shown in Figure A1. As shown there is good agreement to the measured BOLD signal changes (similar results for all stimulation durations – data not shown). However, ε was parameterized at 1.5T and a TE of 40ms. ε is highly dependent on both field strength and TE, and therefore should be parameterized at 3T and TE=15ms before any definitive conclusion can be made. This is not a straightforward task and is work for the future.
Comparing VASO and BOLD fMRI techniques Lu & van Zijl (2005) showed that the EV space accounts for 67% of the total BOLD signal changes. We applied this 2:1 EV to IV ratio to our measured total BOLD signal changes to estimate the changes in the EV compartment for direct comparison to MC simulations of MR signal attenuation.
To validate this assumption the IV effects in the BOLD signal can be eliminated using bipolar ‘crushing’ gradients (Boxerman et al., 1995; Zhong et al., 1998). It is beyond the remit of this manuscript to describe this method in detail and therefore the interested reader is referred to the work of Boxerman et al. (1995). In brief, to investigate suppression of IV effects we utilized an interleaved GE-EPI crushing/no-crushing pulse sequence. A series of Gradient Echo (GE) echo-planar images (EPI) were acquired with varying degrees of diffusion weighting (b ~ 0, 13, 27, 40, 53, 67, 133, 333 s.mm−2). We note that to allow for the bi-polar gradient the echo time was increased to 27ms. All other sequence parameters were the same as those described above. b values of 133/333 s/mm2 correspond to ‘crushed’ velocities of 2.88/1.82 mm/s. In a study on cat intra-cortical capillaries (mean length 188µm, mean diameter 5.7µm), the mean velocity was found to be 2.10mm/s (Pawlik et al., 1981) and so these b-values should eliminate capillary IV effects.
14 separate animals were used to investigate the effect of crushing on the BOLD signal response to electrical stimulation of the whisker pad (16s or 2s, 1.6mA, 5Hz). The paradigm for electrical stimulation of the whisker pad was the same as described above. 8 of these animals underwent a 2 minute 10% hypercapnia challenge to investigate the effect of gradient size on BOLD signal magnitude. For the hypercapnia experiments we acquired 300s (60s baseline, 120s on, 120s off) of GE-EPI (TE=27ms, TR=1s, Av=2) data between a diffusion weighting (bi-polar, x, y & z directed, δ=6ms) of 0 and b for each image.
The BOLD signal change in response to 10% hypercapnia was measured in a barrel cortex region as a function of diffusion weighting (b value). The relationship is shown in Figure A2a. The magnitude of the BOLD signal changes appears to plateau with large diffusion weighting (b value). This suggests the major intravascular effects are ‘crushed’ at the b values corresponding to the range of capillary flow velocities. Our data imply that the ratio of EV to IV effects in the BOLD signal at 3T is ~2:1. Previous work by Boxerman et al. (1995) at 1.5T showed a 1:2 ratio. As field strength increases the IV contribution is expected to decrease. This supports our resulting EV/IV ratio. The initial dip seen at low b (~10s/mm2) has been modelled by Zhong et al. (1998) where the signal is sensitive to the inflow effect. Their model suggests that a greater than 29% increase in inflow velocity will cause a dip in the fractional signal change. Our data support this possibility but require further investigation to confirm it.
We find the same ratio (2:1) between EV and IV signal contributions for long (Figure A2b) and short (Figure A2c) duration electrical stimuli. The magnitude of the BOLD signal in response to short duration (2s) electrical stimuli (ISI=25s) was very small but spectral analysis revealed the EV/IV signal ratio (~ 2:1 at the fundamental frequency 0.04Hz). It is important to note that in this parallel study the inclusion of a bi-polar pulse has increased TE compared with the MR data taken with concurrent optical imaging. There is also a possibility (in large vessels with homogeneous velocities for example) that full IV suppression (highest b-value) may not be possible with single direction flow suppression and some IV effects may still remain in the BOLD signal (Henkelman et al., 1992; Neil and Ackerman, 1992). However, using gradient echo at similar field strengths to that used in this current study, Duong et al. (2003) showed that the potential residual blood contribution was insignificant. However, as the 2:1 ratio was constant across studies and in full agreement with the work of Lu & van Zijl (2005) we feel justified in the application of this ratio to our current 3T data (where concurrent OIS measurements of blood volume and oxygenation were taken) to allow comparison with EV predictions from the MSC-LuT of MR signal attenuation.
Kennerley et al. (2005a) showed that MRI-measured changes in CBV, following electrical stimulation of the whisker pad, were commensurate with Hbt changes (measured with concurrent OIS) indicating a constant haematocrit level. This important finding allowed, in the current work, CBV and Hbt changes to be used interchangeably and thus we used an MR profile of CBV changes as a profile of normalised changes in Hbt for use in the heterogeneous tissue model. However, there may have been some BOLD contamination of the CBV-weighted MRI (Harel et al., 2003) generating a systematic error in the MION signal attenuation. The BOLD signal (a decrease in paramagnetism) will oppose the MION signal (an increase in paramagnetism). Mandeville et al. (2004) assumed that the relaxation rate due to the MION, R2*MION, is proportional to local CBV. However, due to competing BOLD effects R2*MION cannot be measured directly. Assuming additive relaxation rates during activation the measured relaxation rate is given by:
Measurements of CBV changes are systematically underestimated, due to the BOLD effect, with a percentage error given by:
Using data from Kennerley et al. (2005a) we estimate there to be a 7% error in measurements of fractional CBV changes. This systematic error was within the statistical error across subjects. We therefore conclude that the 10mg/kg dose of MION used in our experiments was sufficiently high enough that BOLD contamination of CBV-weighted MRI was minimised. Lu et al. (2004) who also concluded that BOLD contributions have little influence on the CBV-fMRI signal used a similar MION dose (12mg/kg) to that applied in this work.