**ctHbO**) and deoxyhemoglobin (ctHb) in three dimensions. Using this algorithm, the maximum changes in

_{2}**ctHbO**and ctHb were found to occur at

_{2}**0.29**±

**0.02**and

**0.66**±

**0.04**

**mm**beneath the surface of the cortex, respectively. Rytov tomographic reconstructions revealed maximal spatially localized increases and decreases in

**ctHbO**and ctHb of

_{2}**321**±

**53**and

**555**±

**96 nM**, respectively, with these maximum changes occurring at

**4**±

**0.2**

**s**poststimulus. The localized optical signals from the Rytov approximation were greater than those from modified Beer–Lambert, likely due in part to the inability of planar reflectance to account for partial volume effects.

## 1.

## Introduction

The study of cerebral hemodynamics is important for treating brain pathologies, neuroimaging, and fundamental neuroscience. Cerebral hemodynamics plays a critical role not only in cerebrovascular diseases such as stroke^{1} but also in neurodegenerative diseases such as Alzheimer’s disease^{2} and age-related neurodegeneration.^{3} Hemodynamic activity is also the basis for functional magnetic resonance imaging (fMRI), which does not measure neuronal activity directly but instead measures the hyperemic response that accompanies neuronal activity.^{4}^{,}^{5} Hemodynamics may play a direct role in information processing itself, in addition to its supportive role of delivering glucose and oxygen.^{6} Thus, the ability to quantitatively measure cerebral hemodynamics with high temporal and spatial resolution is important for fundamental neuroscience research and for developing new treatments and therapeutic targets for brain pathologies.

Intrinsic signal optical imaging (ISOI) is a powerful tool for imaging the hemodynamic responses to evoked activity of the brain in animal models.^{7}^{,}^{8} ISOI is a camera-based technique that relies on the fact that hemodynamic activity changes the light intensity remitted from an illuminated brain. It allows for high spatial ($\sim 10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$) and temporal ($<0.1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$) resolution over a wide field of view. It is low cost, noninvasive, and can be used in chronic imaging studies in animal models. However, despite its widespread use, ISOI only provides two-dimensional (2-D) planar images, and thus cannot localize the depth in the cortex at which hemodynamic changes occur.

The need for three-dimensional (3-D) determination of hemoglobin concentration and oxygenation state has been a major factor driving the growth of the field of diffuse optical tomography (DOT).^{9} Researchers have used the temporal profile of pulsed sources^{10}^{–}^{14} and large arrays of fibers^{15}^{–}^{19} to localize the depth of hemodynamic activity in humans. Optical tomography of cerebral hemodynamics of small animals has also been performed using optical fibers^{20}^{,}^{21} and fast scanning.^{22} However, fiber-based approaches tend to be constrained by the number of discrete (fiber optic) source and detector probes that are used, and by the vast dynamic range required for different source–detector separations.

For small-animal imaging experiments in which only the first $\sim 1-2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ of tissue needs to be probed, hyperspectal camera-based imaging may provide a simple approach to acquire depth-resolved images with high spatial and temporal resolution. This would provide ISOI with the ability to determine whether optical signals originate from vasculature in the parenchyma or at the surface of the cortex.^{23}^{,}^{24} Quantifying these signals would provide a more complete picture of the spatiotemporal features of the hemodynamic response and may help to determine the signaling mechanisms underlying the coupling between neuronal and vascular activity.^{23}^{,}^{25} While several studies have used confocal and multiphoton microscopy to look at the flow of red blood cells,^{26}^{,}^{27} the dilation of individual vessels,^{23} and the distribution of oxygen tension,^{28} these approaches are based on exogenous probes and are generally only viable for depths up to $\sim 500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. In addition, microscopy-based methods necessarily probe a limited number of microvascular structures in a small field of view.

Many researchers have performed ISOI at multiple wavelengths and used the differing spectral features of oxy- and deoxyhemoglobin to produce 2-D maps of changes in oxy- and deoxyhemoglobin concentration. Practical constraints typically limit the ability of serial spectral scanning to simultaneously provide both spatial and spectroscopic information. Spectroscopic approaches generally sacrifice imaging performance by either collecting data at a single point or diffracting light passing through a slit in order to obtain a one-dimensional (1-D) “image.”^{29}^{–}^{31} Other methods have used a filter wheel or flashing LEDs to serially obtain images at multiple wavelengths.^{32}^{–}^{34} While effective, these methods are limited by the trade-off between the number of wavelengths measured and temporal resolution. Recent attempts to overcome these limitations include approaches using a combination of beam splitters and filters to image at four wavelengths simultaneously^{35} and methods using concurrent CCD-based imaging and fiber-based spectroscopy in a single experiment.^{36} Ideally, one would like to acquire high-resolution images that contain full spectral information at every location on the cortex, and use this data to reconstruct 3-D maps of the changes in tissue hemodynamics.

In this paper, we describe image mapping spectroscopy (IMS) to image intrinsic optical signals following sensory stimulation in the rat. IMS is a new method for hyperspectral imaging capable of forming planar 2-D images in multiple distinct wavelength bins.^{37}^{–}^{39} In this work, we employed an IMS camera operating at 5 Hz to image with 38 spectral bands from 484 to 652 nm in a single snapshot. Unlike other hyperspectral imaging methods, IMS has high optical throughput and does not require scanning or image reconstruction to produce a set of spectrally binned images. As a result, IMS can acquire hyperspectral images at the frame rate of the camera. Thus, IMS is ideal for high-resolution hyperspectral imaging of the small and transient changes in light intensity following sensory stimulation.

To take advantage of the IMS spectral content, we developed a hyperspectral optical tomography algorithm to generate 3-D images of poststimulus changes in oxy- and deoxyhemoglobin concentrations from the acquired reflectance data. Multiwavelength depth localization has recently been studied by several groups for subsurface fluorescence imaging,^{40}^{–}^{42} in which distortions in the emission profiles resulting from the absorption of the emitted light as it travels to the surface are used to estimate the depth of the fluorescent object. Our tomographic reconstruction algorithm, based on the Rytov approximation, takes advantage of the distortion of the reflectance spectra to localize the depth at which the concentrations of tissue absorbers (i.e., oxy- and deoxyhemoglobin) change. Previous studies have employed the Rytov approximation for tomographic reconstruction of data taken at multiple spatial points^{43}^{,}^{44} or with light of multiple spatial frequencies.^{45}^{,}^{46} However, to our knowledge, this work is the first to employ the Rytov algorithm for tomographic reconstruction of hyperspectral optical imaging data. A key feature of the algorithm is its speed, as each hemodynamic time course consists of many time points (72 in the present study). By performing the calculation in the spatial Fourier domain and directly solving for the concentrations of oxy- and deoxyhemoglobin, the algorithm takes $<15\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$ to compute 3-D images at each time point using the large ($\sim 2.5\times {10}^{6}$ data points) hyperspectral data sets. The combination of hyperspectral data collection methods and DOT algorithms potentially enhances ISOI by providing new information for localization and quantification of hemodynamic signals. To compare our results with established approaches and utilize the power of broadband spectral content, we have also calculated hemodynamic signal changes using the modified Beer–Lambert law, where the widely varying photon path length was accounted for using Monte Carlo modeling.

## 2.

## Methods

## 2.1.

### Animal Preparation

Rats were prepared for imaging using our standard protocol that has been described previously.^{47}^{,}^{48} Briefly, five male Sprague-Dawley rats were anesthetized with an intraperitoneal injection of Nembutal ($55\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mg}/\mathrm{kg}$ b.w.). Rats were also given an intramuscular injection of atropine ($0.05\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mg}/\mathrm{kg}$ b.w.) into the hind leg, and a subcutaneous injection of 5% glucose (3 mL). Supplemental injections of Nembutal ($27.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mg}/\mathrm{kg}$ b.w.) were given as necessary (approximately $1/\mathrm{h}$). Body temperature was monitored by a rectal probe and maintained at 37°C using a self-regulating heating blanket. Skin and muscle were removed from the area over the left primary somatosensory cortex, and a $6\times 8\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ imaging window was made by thinning the skull to a thickness of $\sim 150\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ using a dental drill. Petroleum jelly was built around the thinned skull region, filled with saline, and capped with a glass coverslip, keeping the thinned skull moist and translucent during the imaging procedure. All performed procedures were in compliance with National Institutes of Health guidelines and approved by the University of California Animal Care and Use Committee.

## 2.2.

### Imaging Protocol

In each imaging trial, sensory stimulation consisted of five deflections at 5 Hz to the right C2 whisker. For each deflection, a copper wire displaced the whisker 0.8 mm along the rostral–caudal axis at a distance of 5 mm from the C2 follicle for a duration of 100 ms. Images were acquired with an exposure time of 50 ms at a rate of 5 Hz beginning 1.5 s before stimulus onset and lasting 15 s. For each rat, two blocks of 64 trials each were acquired with an intertrial interval of 15 s plus a random delay between 0 and 5 s. The imaging window was evenly illuminated by broadband light from a tungsten-halogen quartz lamp.

## 2.3.

### Hyperspectral Imaging

We have described our IMS system previously.^{37}^{–}^{39} Briefly, light remitted from the cortex is focused on a custom-made redirecting mirror composed of many long strips. Each strip has a slightly different 2-D tilt angle that reflects light in a different direction. A prism subsequently disperses the light from each strip before it arrives at a large-format CCD array. Each pixel on the CCD array corresponds to the light within a narrow wavelength band remitted from a point on the cortex. Therefore, images at the various wavelength bands are obtained by a simple mapping procedure, eliminating the need to reconstruct images as is done with other hyperspectral snapshot techniques such as computed tomography imaging spectrometry^{49} and coded aperture snapshot spectral imaging.^{50} For this paper, we acquired intrinsic signal images simultaneously at 38 separate wavelengths between 484 and 652 nm over a grid of $350\times 350\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{pixels}$ using the image mapping spectrometer system. Our field of view was approximately $6\times 6\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ such that each pixel corresponded to a $17\times 17\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ area of the cortex. The width of each wavelength band $\mathrm{\Delta}\lambda (\lambda )$ varied continuously from $\mathrm{\Delta}\lambda =3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ at $\lambda =484\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ to $\mathrm{\Delta}\lambda =8\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ at $\lambda =652\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$. A representative hyperspectral snapshot of a rat cortex acquired with a single 50-ms exposure appears in Fig. 1. The key advantage provided by the IMS is that it simultaneously acquires images at all 38 wavelengths, enabling brain imaging at multiple frames per second, a critical feature for resolving the rapid hemodynamic signatures of neurovascular coupling in response to stimulus.

## 2.4.

### Data Analysis

## 2.4.1.

#### Intrinsic signal

For each block of 64 trials, we averaged the imaged reflectance at each time point over the 64 trials. For each time point, the averaged data set consisted of a data cube consisting of 38 $350\times 350\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{pixel}$ images, with each image corresponding to one wavelength bin. We averaged the three exposures immediately preceding stimulation onset to calculate the baseline intensity. The intrinsic signal for each time point was then calculated for each pixel and wavelength separately as $\phi (t)=\mathrm{log}[R(t)/{R}_{0}]$, where $R(t)$ and ${R}_{0}$ are the measured light intensity at time $t$ and the baseline intensity, respectively. We related the intrinsic signal $\phi (t)$ to the changes in oxy- and deoxyhemoglobin in the cortex resulting from sensory stimulation using both the modified Beer–Lambert law and the Rytov approximation.

## 2.4.2.

#### Modified Beer–Lambert law

In order to generate 2-D images of the changes in hemoglobin concentrations, we used the modified Beer–Lambert law, which appears as

## Eq. (1)

$$\phi (t)=-[{\u03f5}_{{\mathrm{HbO}}_{2}}(\lambda )\delta {\mathrm{ctHbO}}_{2}(t)+{\u03f5}_{\mathrm{Hb}}(\lambda )\delta \mathrm{ctHb}(t)]L(\lambda ).$$^{51}$\delta {\mathrm{ctHbO}}_{2}(t)$ and $\delta \mathrm{ctHb}(t)$ are the changes in tissue concentration of those chromophores relative to their baseline values, and $L(\lambda )$ is the mean photon path length in the cortex as a function of wavelength. Monte Carlo simulations were used to determine the mean path lengths traveled by photons of differing wavelengths. In our Monte Carlo simulations, we began by using tissue oxy- and deoxyhemoglobin concentrations of ${\mathrm{ctHbO}}_{2}=60\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{M}$ and $\mathrm{ctHb}=40\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{M}$, respectively, tissue-scattering coefficient of ${\mu}_{\mathrm{s}}=10\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, and tissue-scattering anisotropy of $g=0.9$, similar to values we have measured previously in this animal model using quantitative spatial frequency domain imaging.

^{52}

## 2.4.3.

#### Rytov approximation

While the modified Beer–Lambert law has been used extensively to produce images of changes in hemoglobin concentration, Eq. (1) assumes that the changes occur uniformly throughout the brain. The Rytov approximation is a generalization of the modified Beer–Lambert law that allows one to solve for spatially varying changes in chromophore concentrations.^{53} By taking advantage of the differing penetration depths of visible and near-infrared light in the cortex, we can implement the Rytov approximation to generate 3-D images of the changes in Hb and ${\mathrm{HbO}}_{2}$ concentrations.

The propagation of multiple-scattered incoherent light in tissue is governed by the radiative transport equation (RTE)

## Eq. (2)

$$\widehat{\mathbf{s}}\xb7\nabla I(\mathbf{r},\widehat{\mathbf{s}})+{\mu}_{\mathrm{t}}I(\mathbf{r},\widehat{\mathbf{s}})={\mu}_{\mathrm{s}}\int {\mathrm{d}}^{2}sA(\widehat{\mathbf{s}},{\widehat{\mathbf{s}}}^{\prime})I(\mathbf{r},\widehat{\mathbf{s}})+S(\mathbf{r},\widehat{\mathbf{s}}),$$## Eq. (3)

$$\phi (\lambda ,{\mathbf{\rho}}_{\mathrm{d}},{\widehat{\mathbf{s}}}_{\mathrm{d}})=-\int {\mathrm{d}}^{3}r\text{\hspace{0.17em}}{\mathrm{d}}^{2}sI(\lambda ,z,\widehat{\mathbf{s}})G(\lambda ,\mathbf{r},\mathbf{s},{\mathbf{r}}_{\mathrm{d}},{\widehat{\mathbf{s}}}_{\mathrm{d}})\delta {\mu}_{\mathrm{a}}(\lambda ,\mathbf{r}),$$## Eq. (4)

$$\phi (\lambda ,{\mathbf{\rho}}_{\mathrm{d}},{\widehat{\mathbf{s}}}_{\mathrm{d}})=I(\lambda ,z=0,\widehat{\mathbf{s}})\mathrm{log}\left[\frac{R(\lambda ,{\mathbf{\rho}}_{\mathrm{d}},{\widehat{\mathbf{s}}}_{\mathrm{d}})}{{R}_{0}(\lambda ,{\mathbf{\rho}}_{\mathrm{d}},{\widehat{\mathbf{s}}}_{\mathrm{d}})}\right]$$## Eq. (5)

$$\delta {\mu}_{\mathrm{a}}(\lambda ,\mathbf{r})={\u03f5}_{\mathrm{Hb}}(\lambda )\delta \mathrm{ctHb}(\mathbf{r})+{\u03f5}_{{\mathrm{HbO}}_{2}}(\lambda )\delta {\mathrm{ctHbO}}_{2}(\mathbf{r}).$$^{54}

^{,}

^{55}Using this correction, we can rewrite Eq. (4) as

## Eq. (6)

$$\phi (\lambda ,{\mathbf{\rho}}_{\mathrm{d}})=\int {\mathrm{d}}^{3}r[I(\lambda ,z)G(\lambda ,\mathbf{r},{\mathbf{r}}_{\mathrm{d}})-\frac{\ell {*}^{2}}{3}\nabla I(\lambda ,z)\xb7\nabla G(\lambda ,\mathbf{r},{\mathbf{r}}_{\mathrm{d}})]c\delta {\mu}_{\mathrm{a}}(\lambda ,\mathbf{r}).$$## Eq. (7)

$$[D(\lambda ){\nabla}^{2}+c{\mu}_{\mathrm{a}}(\lambda )]G(\lambda ,\mathbf{r},{\mathbf{r}}^{\prime})=\delta (\mathbf{r}-{\mathbf{r}}^{\prime}),$$Equation (7) is the most easily inverted in the Fourier domain where it is block diagonal. Taking a 2-D Fourier transform with respect to the detector position ${\mathbf{\rho}}_{\mathrm{d}}$ yields a 1-D integral for each value of the Fourier coefficient $\mathbf{q}$

## Eq. (8)

$$\tilde{\phi}(\lambda ,\mathbf{q})=\int \mathrm{d}z\kappa (\lambda ,\mathbf{q},z)\delta {\tilde{\mu}}_{\mathrm{a}}(\lambda ,\mathbf{q},z).$$^{55}In the semi-infinite geometry, we have

## Eq. (9)

$$\kappa (\lambda ,\mathbf{q},z)={\left(\frac{\ell}{D}\right)}^{2}[1-\frac{\ell {*}^{2}}{3}Q(\lambda ,q)Q(\lambda ,0)]\left(\frac{\mathrm{exp}\{[Q(\lambda ,q)+Q(\lambda ,0)]z\}}{[Q(\lambda ,q)\ell +1][Q(\lambda ,0)\ell +1]}\right),$$^{56}and $Q(\lambda ,q)={(c{\mu}_{\mathrm{a}}/D+{q}^{2})}^{1/2}$. The changes in ctHb and ${\mathrm{ctHbO}}_{2}$, relative to baseline, are then calculated according to

## Eq. (10)

$$\delta \mathrm{ctHb}(\mathbf{\rho},z)=\int \frac{{\mathrm{d}}^{2}q}{{(2\pi )}^{2}}\text{\hspace{0.17em}}\mathrm{exp}(-i\mathbf{q}\xb7\mathbf{\rho})\sum _{m,n}{\u03f5}_{\mathrm{Hb}}({\lambda}_{m})\kappa *({\lambda}_{m},\mathbf{q},z){M}_{m,n}^{-1}(\mathbf{q})\tilde{\phi}({\lambda}_{n},\mathbf{q}),$$## Eq. (11)

$$\delta {\mathrm{ctHbO}}_{2}(\mathbf{\rho},z)=\int \frac{{\mathrm{d}}^{2}q}{{(2\pi )}^{2}}\text{\hspace{0.17em}}\mathrm{exp}(-i\mathbf{q}\xb7\mathbf{\rho})\sum _{m,n}{\u03f5}_{{\mathrm{HbO}}_{2}}({\lambda}_{m})\kappa *({\lambda}_{m},\mathbf{q},z){M}_{m,n}^{-1}(\mathbf{q})\tilde{\phi}({\lambda}_{n},\mathbf{q}).$$## Eq. (12)

$${M}_{m,n}(\mathbf{q})={\int}_{0}^{\infty}\mathrm{d}z[{\u03f5}_{\mathrm{Hb}}({\lambda}_{m}){\u03f5}_{\mathrm{Hb}}({\lambda}_{n})+{\u03f5}_{{\mathrm{HbO}}_{2}}({\lambda}_{m}){\u03f5}_{{\mathrm{HbO}}_{2}}({\lambda}_{n})]\kappa ({\lambda}_{m},\mathbf{q},z)\kappa *({\lambda}_{n},\mathbf{q},z),$$The framework of this Rytov-approximation-based algorithm has previously been shown to provide accurate estimates of depth of buried fluorescent inclusions in small animals.^{46} In addition, we have recently performed a validation of the hyperspectral tomography algorithm on data generated using Monte Carlo simulations^{57} with a 2-mm thick absorbing inclusion buried 0.5 mm beneath the surface. For this case (Fig. 2), the multispectral tomography algorithm was able to accurately reconstruct the center depth of the inclusion to within 23% and the inclusion thickness to within 18%. A more detailed discussion of tomography validation studies, employing both tissue-simulating phantoms and computational models of the forward problem, will be the subject of a future report (currently in preparation).

## 3.

## Results and Discussion

Figure 3 demonstrates our ability to image the changes in light intensity due to functional activation with high spatial and temporal resolution simultaneously at many wavelengths through a thinned skull in a representative rat. Figure 3(a) shows images of changes in light remitted from the cortex at 529 nm (isosbestic point equally sensitive to oxy- and deoxyhemoglobin) and 630 nm (about $9\times $ more sensitive to deoxyhemoglobin than oxyhemoglobin).

Color maps of the change in light intensity are superimposed over black and white images at 576 nm, which clearly show the underlying vasculature. As can be seen by looking at the time dependence of the signal over the area of peak activation shown in Fig. 3(b), the time points chosen in Fig. 1(a) correspond to the initial decrease in the 630 nm signal, the time at which the isosbestic wavelength (529 nm) has its largest signal, and the increase of signal at 630 nm. Complete spectral information for these three time points is shown in Fig. 3(c), and the intrinsic signal is represented in Fig. 3(d) in false color as a function of both wavelength and time. These results are in agreement with our previous ISOI imaging at a single wavelength.^{48}

Representative images of the changes in ctHb and ${\mathrm{ctHbO}}_{2}$ along with traces of the region of interest calculated using the modified Beer–Lambert law appear in Figs. 4(a) and 4(b). Changes in concentration (Table 1, averaged over all rats) were determined using the modified Beer–Lambert law, where optical path lengths were determined using Monte Carlo simulations of photon transport. The peak increase in ${\mathrm{ctHbO}}_{2}$ ($\delta {\mathrm{ctHbO}}_{2}=248\pm 80\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nM}$) was $\sim 3\times $ greater than the ctHb decrease ($\delta \mathrm{ctHb}=-84\pm 20\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nM}$). The peak changes in concentration did not occur at the same time; the peak change in total hemoglobin arose first (${t}_{\mathrm{HbT}}=2.6\pm 0.4\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$), followed by the peak in ${\mathrm{ctHbO}}_{2}$ (${t}_{{\mathrm{HbO}}_{2}}=3.1\pm 0.3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$), and the maximum decrease in ctHb (${t}_{\mathrm{Hb}}=4.3\pm 0.3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$).

## Table 1

Peak times and maximum changes in the concentrations of oxyhemoglobin (ctHbO2) and deoxyhemoglobin (ctHb) using the modified Beer–Lambert law and the Rytov approximation.

Modified Beer–Lambert law | ||
---|---|---|

ctHbO2 increase | ctHb decrease | |

Time to peak (s) | $3.1\pm 0.3$ | $4.3\pm 0.3$ |

Magnitude (nM) | $248\pm 80$ | $-84\pm 19$ |

Rytov approximation | ||

Time to peak (s) | $4\pm 0.2$ | $4\pm 0.2$ |

Magnitude (nM) | $321\pm 53$ | $-555\pm 96$ |

We varied the number of wavelengths used to fit for ctHb and ${\mathrm{ctHbO}}_{2}$ to determine the effect on the recovered changes in hemoglobin concentration. We found that using fewer wavelengths led to additional noise in the recovered hemoglobin values. For example, Figs. 4(b) and 4(c) shows the hemoglobin time courses for a single rat using all 38 wavelengths [Fig. 4(b)], and only using two wavelengths [Fig. 4(c)]: an isosbetic wavelength (529 nm) and a deoxyhemoglobin sensitive wavelength (630 nm). When reducing the number of wavelengths from 38 to 2, the standard error of hemoglobin concentration changes in the regions of interest increased from 6 to 20 nM for ${\mathrm{ctHbO}}_{2}$ and from 3.5 to 5.6 nM for ctHb. However, despite the increase in noise, removing wavelengths did not significantly alter the timing or peak magnitudes of the hemodynamic time courses.

We also conducted Monte Carlo simulations over a wide range of parameters to determine the effect of path length on the calculated values for changes in ctHb and ${\mathrm{ctHbO}}_{2}$. Figure 5 illustrates the variability in the magnitude of hemoglobin changes calculated when different baseline optical properties are assumed. The variability of the initial ctHb increase, the ${\mathrm{ctHbO}}_{2}$ increase, and the subsequent ctHb decrease are reported in Table 2. Increasing the baseline value of total hemoglobin or scattering leads to an increase in the magnitude of all parameters. This occurs because an increase in absorption or scattering decreases the mean photon path length. Since the data remain fixed, this leads to an increase in the calculated hemoglobin concentrations (see Sec. 2). All changes due to varying the wavelength dependence of scattering, described by $A{\lambda}^{-b}$ where the $A$ and $b$ parameters refer to scattering amplitude and power, respectively,^{46}^{,}^{52} were less than 15% despite changing scattering power, $b$, over a large range ($b=0\u20132.5$).

## Table 2

Magnitude of changes in hemoglobin concentrations for the differing optical properties used in the Monte Carlo simulations.

Parameter varied in Monte Carlo simulation | ||||||||
---|---|---|---|---|---|---|---|---|

HbT (ctHbO2+ctHb) | StO2 (ctHbO2/HbT) | μs (Aλ−b) | b | |||||

50 μM | 20 μM | 50% | 100% | 5 mm−1 | 20 mm−1 | 0 | 2.5 | |

${\mathrm{ctHbO}}_{2}$ increase (nM) | 142 | 450 | 237 | 247 | 225 | 283 | 241 | 259 |

ctHb decrease (nM) | $-52$ | $-143$ | $-84$ | $-67$ | $-71$ | $-103$ | $-83$ | $-86$ |

Future work will include a quantitative comparison between the tomography model used in this paper and a Monte Carlo-based tomography model that uses photon path length information stored in the simulations to provide accurate quantification of absorber depth and concentration. Future work will also include a detailed computational investigation into the influence of tissue background absorption and scattering properties on the depth sensitivity of the optical tomography algorithm. This future study will also include a systematic investigation into how the background optical properties and the set of wavelengths employed affects the ability of the algorithm to fully separate oxygenated hemoglobin absorption from deoxygenated hemoglobin absorption. This future work will include a quantitative comparison between the tomography model used in this paper and a Monte Carlo-based tomography model that uses photon pathlength information stored in the simulations to provide accurate quantification of absorber depth and concentration.

Figure 6 shows representative depth-resolved snapshots of the changes in oxy- and deoxyhemoglobin at 4 s poststimulus for one rat. These images were calculated according to the Rytov approximation method described in Sec. 2.4.2. Eleven wavelengths from 586 to 652 nm were employed in the reconstructions. Slices of the 3-D images of ctHb and ${\mathrm{ctHbO}}_{2}$ are shown at six depths from $100\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ to 1.1 mm beneath the surface of the cortex. Depth profiles of the changes in oxy- and deoxyhemoglobin, as a function of depth beneath the surface of the cortex, are shown in Fig. 7(a) for the same rat as in Fig. 6. Figure 7(b) shows time courses of these changes in oxy- and deoxyhemoglobin, at a depth of 0.3 mm beneath the surface of the cortex, averaged over four rats.

The time courses shown in Fig. 7(b) illustrate an initial decrease in the concentration of oxyhemoglobin and an initial increase in the concentration of deoxyhemoglobin over the first 2 s poststimulus, followed by an increase in the oxyhemoglobin concentration and a decrease in the deoxyhemoglobin concentration over the subsequent 5 s. Using the Rytov tomography algorithm, the maximum changes in ${\mathrm{ctHbO}}_{2}$ and ctHb were found to occur at depths of $0.29\pm 0.02$ and $0.66\pm 0.04\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ beneath the surface of the cortex, respectively. The peak times and magnitudes of concentration changes for the 3-D image reconstructions are summarized in Table 1. For the Rytov tomography algorithm, the point of peak concentration change (Table 1) was defined as the peak time point at the depth ($z=0.3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ beneath the surface of the cortex) at which the peak change in oxyhemoglobin concentration occurred. The magnitudes of the peak changes in oxy- and deoxyhemoglobin concentration shown in Table 1 are larger for the Rytov tomography algorithm than the Beer–Lambert method; this finding is expected because the Rytov tomography algorithm localizes the changes in 3-D space.

The results shown in Figs. 6 and 7 are indicative of the depletion of oxygen by the cortex in response to the stimulus, followed by the reperfusion of oxygenated blood into the region. Figures 6 and 7 also suggest that the relative magnitudes of the changes in oxy- and deoxyhemoglobin may vary with depth. This effect could possibly be explained by differences in depth of diving arterioles and venuoles, with the extraction of oxygen occurring deeper in the cortex and the reperfusion of oxygenated blood occurring closer to the surface. However, this effect may also represent a limitation of the ability of the measurement technique and/or mathematical model to fully resolve the depths of two different absorbers (oxy- and deoxyhemoglobin), whose spectral features are close together in the interrogated wavelength range.

It is important to note that the peak magnitudes of concentration change were higher in the 3-D tomographic images (using the Rytov approximation) than for the modified Beer–Lambert law analysis. This is not surprising since the modified Beer–Lambert law results represent a depth-averaged hemodynamic change, while the tomographic results localize this information. When the reconstructed 3-D maps of concentration changes were averaged over depths from 0.1 to 1.3 mm beneath the surface of the cortex to create 2-D maps of changes in oxy- and deoxyhemoglobin, the magnitude of the change in oxyhemoglobin concentration decreased to $154\pm 29\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nM}$, but the change in deoxyhemoglobin concentration remained essentially the same as in the 3-D map ($555\pm 99\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nM}$). This result is related to the inability of the Rytov tomography algorithm to localize the changes in deoxyhemoglobin concentration as precisely as the changes in oxyhemoglobin concentration (Figs. 6 and 7). This discrepancy could be due in part to the fact that the absorption features of oxy- and deoxyhemoglobin were not completely decoupled from each other due to the size of the wavelength bands detected by the image mapping spectrometer. Future studies will quantitatively examine this issue by convolving the oxy- and deoxyhemoglobin absorption spectra with the bandwidth profile $\mathrm{\Delta}\lambda (\lambda )$ of the detector.

## 4.

## Conclusion

We have employed a parallel hyperspectral imaging technique, IMS, to record intrinsic optical signals in the rat somatosensory cortex over 38 wavelength bands (from 484 to 652 nm) at 5 Hz following mechanical movement of a single whisker. This approach was combined with a hyperspectral optical tomography algorithm to quantify and localize time-dependent concentration changes in oxygenated and deoxygenated hemoglobin that are characteristic of the response to the evoked stimulus. Our image reconstruction algorithm uses the Rytov approximation and takes advantage of the fact that the measured spectral changes in reflectance have a strong depth-dependent distortion. Data were also analyzed using a modified Beer–Lambert law in which the photon pathlength was accounted for using Monte Carlo modeling.

The Rytov approach revealed that the maximum measured changes in ${\mathrm{ctHbO}}_{2}$ and ctHb were primarily confined to layers II–IV of the cortex. However, the depth distribution of the Rytov-calculated ctHb change had a broader tail than that of the ${\mathrm{ctHbO}}_{2}$ change. Rytov tomographic reconstructions revealed spatially localized increases and decreases in ${\mathrm{ctHbO}}_{2}$ and ctHb of $321\pm 53$ and $555\pm 96\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nM}$, respectively, at $4\pm 0.2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$ poststimulus. The modified Beer–Lambert model revealed a $248\pm 80\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nM}$ increase in ${\mathrm{ctHbO}}_{2}$ at $3.1\pm 0.3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$ poststimulus and an $84\pm 19\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nM}$ decrease in ctHb at $4.3\pm 0.3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$ poststimulus.

Localized optical BOLD signals for both ${\mathrm{ctHbO}}_{2}$ and ctHb were greater using the Rytov approximation than those estimated using the modified Beer–Lambert approach. One contribution to these differences in signal magnitude likely comes from the inability of planar reflectance methods to account for partial volume effects. The additional discrepancy between the ctHb values calculated with the Rytov and Beer–Lambert approaches is likely an artifact of the inability of the algorithm to completely separate the signal contributions from oxy- and deoxyhemoglobin. This is due, in part, to the fact that we did not rigorously characterize the effect of the wavelength dependence of detector bandwidth $\mathrm{\Delta}\lambda (\lambda )$. Future studies will investigate potential improvements to the Rytov model as well as imposing further optical path length constraints through the use of structured light projection patterns.^{45}^{,}^{46} Despite these limitations, this study provides, for the first time, a general framework for localizing and characterizing hemodynamic signals in the cortex with $\sim \mathrm{mm}$ resolution over a wide, scalable field of view. One application of great potential biomedical importance is to employ multispectral optical tomography in combination with spatial frequency domain and laser speckle contrast imaging methods in order to separate the hemodynamic effects of blood flow from metabolism in the brain. These studies are currently in progress in our laboratory, with emphasis on characterization of the brain’s response to traumatic perturbations such as cardiac arrest and stroke in preclinical animal models.

## Acknowledgments

Support for this work was provided by the NIH Laser Microbeam and Medical Program (LAMMP, P41-EB015890), NIH R21 NS078634, R21 EB014440, R21EB009186, and R01CA124319, and NINDS NS 066001 and NS 055832. We acknowledge the use of Virtual Photonics software resources ( http://virtualphotonics.org/) supported by LAMMP. Soren D. Konecky was supported by a postdoctoral fellowship from the Hewitt Foundation for Medical Research. Robert H. Wilson is also supported by a postdoctoral fellowship from the Hewitt Foundation for Medical Research. Beckman Laser Institute programmatic support from the Beckman Foundation is also acknowledged.