^{-7}) than between the BOLD and oxyhemoglobin or total hemoglobin measures (R=0.38, p=0.04|0.37, p=0.05, respectively). Similarly, we find that the correlation between the ASL measured blood flow and optically measured total and oxyhemoglobin is stronger (R=0.73, p=5×10

^{-6}and R=0.71, p=9×10

^{-6}, respectively) than the flow to deoxyhemoglobin spatial correlation (R=0.26, p=0.10).

## 1.

## Introduction

Diffuse optical imaging (DOI) is a spectroscopic method capable of noninvasively measuring changes in the concentrations of oxy- and deoxyhemoglobin in the human brain during functional activity.^{1, 2} As a tool for neurology, this technique offers several additional attributes that complement other existing imaging methods, such as functional MRI (fMRI).^{3} MRI and optical imaging complement each other by their spatial, temporal, and spectroscopic resolutions. In recent years, this has led to an increased interest in multimodality functional human imaging^{4} and image reconstruction techniques.^{5, 6, 7} However, as these multimodality methods are developed, careful studies must be performed to explore the underlying biophysical relationships between the different imaging modalities. Properly interpreting data together from both modalities requires an understanding of the relationships between the underlying physiological hemodynamic states and the physical measurements obtained from the instruments and, with that, experimental validation of the biophysical models supporting each modality.

In recent years, a number of studies have examined these relationships through investigations of the temporal correspondence between DOI and blood oxygen level-dependent (BOLD)-fMRI,
^{8, 9, 10, 11, 12} DOI and arterial spin labeling (ASL)-fMRI,^{12} or DOI and cerebral blood volume (CBV)-fMRI.^{11} These and many other collected works have attempted to reconcile theoretical models describing the relationship between these fMRI contrast mechanisms and the physiological deoxyhemoglobin, blood volume, and blood flow parameters through the measurements recorded by optical imaging. While it is clear that the differing vascular sensitivities of the modalities may play a role in precise details of this relationship, the conclusion of many of these studies has been that the temporal dynamics of measurements by fMRI and optical imaging are consistent with their respective theoretical biophysical underpinnings. However, although there have been many such investigations, which have explored the temporal relationships between these two modalities, little has been published exploring the spatial correspondence between fMRI and DOI in a quantitative way. Although a few qualitative reports suggest agreement between the spatial profiles of fMRI and optical methods,^{13} no detailed and quantitative comparisons have been published. This lack of quantitative comparisons might be attributed to the ill-posed problem of image reconstruction by optical imaging. Optical imaging is based on the topographic reconstruction of a discrete set of absorption change measurements between pairs of optical source-detector positions on the head. This reconstruction is generally underdetermined, with far more unknown parameters than actual measurements. Tomographic (three-dimensional) reconstruction is even more ill-posed, since the added depth degree of freedom adds more unknowns, as well as, an exponentially decaying sensitivity to absorption changes. Over the years, several groups have continued to make progress with better reconstruction algorithms and methodologies, which exploit anatomical
^{6, 7, 14, 15} or functional^{5} MRI information. By improving the spatial localization of reconstructed optical images, spatial priors may help to improve the quantitative accuracy of DOI. However, for these reconstructions, the method and amount of regularization used to constrain the solutions will have a quantitative effect on the resulting image. As a result, it has been difficult to quantitatively compare the spatial profiles of fMRI and DOI independently of reconstruction technique. Without this independent assessment, we cannot be confident in the use of fMRI information to guide the optical inverse problem.

In this work, we perform the first quantitative assessment of the spatial correlation between optical and fMR imaging. To avoid the optical inverse problem, rather than reconstructing images from the optical data using these regularization techniques, we directly tested the fMRI data as a solution to the optical inverse problem. Here, given the anatomical structure of the head from MRI, we use photon-migration theory and Monte Carlo techniques to calculate the Green’s function describing the optical sensitivity profile to absorption changes in the underlying structures of the head (i.e., the optical “forward” model) as described in Boas
^{16} The overlap of this optical sensitivity profile with the brain activation identified by fMRI provides a means of predicting the response amplitudes, which would be measured optically between each source and detector pair. This approach allows us to account for the decaying depth sensitivity and partial volume errors of the optical measurements and avoid the ill-posed inverse problem by using the optical forward equation in the forward direction. Rather than increasing the dimensionality of the optical data by image reconstruction through regularization, knowledge of the optical measurement model is used to reduce dimensionality of the fMRI data into the source-detector–based (measurement) space of the optical probe, thus allowing quantitative comparisons between DOI and fMRI. Performing comparisons of the spatially and temporally varying hemodynamic response across the probe layout allows us to explore the spatiotemporal relationships between optical and fMRI while avoiding the need for regularized inversion techniques.

In this experimental study, we performed DOI simultaneously with BOLD- and ASL-based fMRI during an event-related finger-tapping task in order to quantitatively test the spatial correlation between the optical and fMRI methods. This data was published in Huppert, ^{12} who investigated the temporal correlation between the region-of-interest averaged signals. Here, we provide a significant follow-up to that analysis to investigate the spatial correlation between these modalities.

## 1.1.

### Theory

## 1.1.1.

#### Photon migration

DOI of brain activation is generally based on the measurement of spatiotemporal changes in the absorption of light. Near-infrared light is introduced into the head and propagates through the dense scattering layers of the scalp and skull into the brain. Hemodynamic changes in oxy- and deoxyhemoglobin concentrations affect the absorption properties of the brain and result in a change in the intensity of light as it migrates back out of the head. Using multiple measurements taken between an array of light source and detector positions spaced several centimeters apart, DOI attempts to spatially resolve these changes in hemoglobin. However, the reconstruction of images requires knowledge of the spatial profile of the light propagation through the head, which determines the spatial sensitivity of these measurements. To derive the distribution of photons in a medium with a complex distribution of absorption and scattering properties [
${\mu}_{a}\left(r\right)$
and
${\mu}_{s}^{\prime}\left(r\right)$
, respectively], such as the human head, the photon density must be modeled by empirical means using computerized simulations. In Monte Carlo–based modeling, such as is performed in this study, photons are “launched” into the medium. The distribution of these photons is statistically modeled based on the probability of an absorption or scattering event at each region of space as described by
${\mu}_{a}\left(r\right)$
and
${\mu}_{s}^{\prime}\left(r\right)$
.^{16, 17, 18}

For brain activation, the changes in the optical properties are small, and thus a linear approximation is reasonably accurate for predicting the changes in optical measurements produced by localized changes in the optical properties. The forward model for such optical measurements is of the form

## Eq. 1

$$y\left(\lambda \right)=\mathsf{A}\left(\lambda \right)\bullet \delta x\left(\lambda \right),$$^{19}Because the propagation of light depends on the optical properties of the medium, each of these parameters ( $y$ , $\mathsf{A}$ , and $\delta x$ ) are dependent on the optical wavelength $\left(\lambda \right)$ .

## 1.1.2.

#### The forward propagation of fMRI

The forward model, shown in Eq. 1, describes the linear transformation of absorption changes, which occur at particular volume locations, to the measurement of these changes between the optodes of the DOI probe. Rather than inverting this equation to achieve an image of these changes as is typical of optical imaging, the forward model can be used to predict how the changes shown by fMRI images should look when measured optically. Thus, Eq. 1 can be used to “forward model” the optical measurements based on the changes measured at high spatial resolution by fMRI. This provides us with a direct way of testing the hypothesis that the optical and fMRI spatial profiles of the response amplitudes are consistent with one another.

However, the direct application of the optical forward operator
$\mathsf{A}$
to project the fMRI image requires information about the relationship between the fMRI signal and optical absorption changes. A change in the fMRI-BOLD signal is proportional to a change in reduced hemoglobin (HbR) and not necessarily proportional to changes in optical absorption at a single given wavelength since the optical absorption is a function of both HbR and oxygenated hemoglobin
$\left({\mathrm{HbO}}_{2}\right)$
changes. Thus, Eq. 1 cannot be directly applied without considering how light propagation differs for each optical wavelength. Rather then using the optical forward operator determined at a single wavelength, we must consider a multispectral forward equation.^{20, 21}

In order to derive a forward operator to directly project hemoglobin changes, we assume the absorption is dominated by hemoglobin. Optical absorption changes are linear combinations of the oxy- and deoxyhemoglobin changes, such that

## Eq. 2

$$\mathrm{\Delta}{\mu}_{a}\left(\lambda \right)={\epsilon}_{{\mathrm{HbO}}_{2}}\left(\lambda \right)\Delta \left[{\mathrm{HbO}}_{2}\right]+{\epsilon}_{HbR}\left(\lambda \right)\Delta \left[HbR\right],$$^{20}):

## 3.

## Eq. 3a

$$\left[\begin{array}{c}\hfill y\left({\lambda}_{1}\right)\hfill \\ \hfill y\left({\lambda}_{2}\right)\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill \mathsf{A}\left({\lambda}_{1}\right)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \mathsf{A}\left({\lambda}_{2}\right)\hfill \end{array}\right]\left[\begin{array}{c}\hfill \delta {\mu}_{a}\left({\lambda}_{1}\right)\hfill \\ \hfill \delta {\mu}_{a}\left({\lambda}_{2}\right)\hfill \end{array}\right]$$## Eq. 3b

$$=\mathsf{AE}\left[\begin{array}{c}\hfill \delta \left[{\mathrm{HbO}}_{2}\right]\hfill \\ \hfill \delta \left[HbR\right]\hfill \end{array}\right]$$## Eq. 3c

$$=\left[\begin{array}{cc}\hfill {\mathrm{\epsilon}}_{{\mathrm{HbO}}_{2}}\left({\lambda}_{1}\right)\mathsf{A}\left({\lambda}_{1}\right)\hfill & \hfill {\mathrm{\epsilon}}_{HbR}\left({\lambda}_{1}\right)\mathsf{A}\left({\lambda}_{1}\right)\hfill \\ \hfill {\mathrm{\epsilon}}_{{\mathrm{HbO}}_{2}}\left({\lambda}_{2}\right)\mathsf{A}\left({\lambda}_{2}\right)\hfill & \hfill {\mathrm{\epsilon}}_{HbR}\left({\lambda}_{2}\right)\mathsf{A}\left({\lambda}_{2}\right)\hfill \end{array}\right]\left[\begin{array}{c}\hfill \delta \left[{\mathrm{HbO}}_{2}\right]\hfill \\ \hfill \delta \left[HbR\right]\hfill \end{array}\right],$$## 3d.

Rather than using Eq. 1 to forward project the fMRI image, the multispectral forward equation given in Eq. 3 should be considered. The BOLD-fMRI signal is hypothesized to represent the spatio-temporal variation in $\delta \left[\mathrm{HbR}\right]$ . Equation 3 can thus be used to predict the optical signals $y\left({\lambda}_{690}\right)$ and $y\left({\lambda}_{830}\right)$ due to a change in deoxy-hemoglobin

## Eq. 4

$$\left[\begin{array}{c}\hfill \mathbf{\Delta}OD\left({\lambda}_{1}\right)\hfill \\ \hfill \mathbf{\Delta}OD\left({\lambda}_{2}\right)\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill {\mathbf{\epsilon}}_{HbO}\left({\lambda}_{1}\right)\mathsf{A}\left({\lambda}_{1}\right)\hfill & \hfill {\mathbf{\epsilon}}_{HbR}\left({\lambda}_{1}\right)\mathsf{A}\left({\lambda}_{1}\right)\hfill \\ \hfill {\mathbf{\epsilon}}_{HbO}\left({\lambda}_{2}\right)\mathsf{A}\left({\lambda}_{2}\right)\hfill & \hfill {\mathbf{\epsilon}}_{HbR}\left({\lambda}_{2}\right)\mathsf{A}\left({\lambda}_{2}\right)\hfill \end{array}\right]\times \left[\begin{array}{c}\hfill \mathbf{0}\hfill \\ \hfill \mathbf{\Delta}\mathbf{BOLD}\hfill \end{array}\right].$$^{22, 23}

## Eq. 5

$$\Delta {\left[HbR\right]}_{BOLD}=\frac{{\mathbf{\epsilon}}_{\mathrm{HbO}}\left({\lambda}_{2}\right)\bullet \mathrm{\Delta}OD\left({\lambda}_{1}\right)\u2215[L\bullet {l}_{DPF}\left({\lambda}_{1}\right)]-{\mathbf{\epsilon}}_{\mathrm{HbO}}\left({\lambda}_{1}\right)\bullet \mathrm{\Delta}OD\left({\lambda}_{2}\right)\u2215[L\bullet {l}_{DPF}\left({\lambda}_{2}\right)]}{{\mathbf{\epsilon}}_{HbR}\left({\lambda}_{1}\right)\bullet {\mathbf{\epsilon}}_{\mathrm{HbO}}\left({\lambda}_{2}\right)-{\mathbf{\epsilon}}_{\mathrm{HbO}}\left({\lambda}_{1}\right)\bullet {\mathbf{\epsilon}}_{HbR}\left({\lambda}_{2}\right)}.$$
In the modified Beer-Lambert equation [Eq. 5],
$L$
is the linear distance between the position of a source and detector and
${l}_{dpf}$
is the differential path-length factor, which is a unitless coefficient defined as the effective path length through the head divided by the source-detector seperation.^{22, 23, 24} Using Eqs. 4, 5, the BOLD signal is changed to a spectrally resolved optical absorption change, projected through the optical forward equation at each measured wavelength, and finally converted back to an estimated concentration change using the modified Beer-Lambert law. Thus, we obtain an estimate of the deoxyhemoglobin changes for each source-detector pair, predicted from the spatial profile of the BOLD signal. This overall procedure is similar to the method detailed by Strangman
^{24} for examining the errors in the estimate of the hemoglobin changes arising from the use of the modified Beer-Lambert equation. These errors have been found to be negligible for the wavelengths used in our study.^{24, 25, 26} Similarly, the fMRI-ASL signal depicts changes in blood flow, which we have previously noted has a strong temporal correlation with total hemoglobin changes
$\left(\right[\mathrm{HbT}]=[{\mathrm{HbO}}_{2}]+[\mathrm{HbR}\left]\right)$
.^{12} Again, because this represents absorption changes at multiple wavelengths and must be propagated through a multispectral forward model, we approximate the ASL signal as a spatiotemporal variation in
$\delta \left[\mathrm{HbT}\right]$
. We predict the optical signals measured for each source-detector pair using Eqs. 3, 4, 5, but now projecting both oxy- and deoxyhemoglobin changes assuming an oxygen saturation of 65%.^{24, 27}

## 2.

## Methods

The data used in this analysis has been previously reported in Huppert
^{12} for comparison of the temporal characteristics of optical and fMR imaging. The methods are briefly summarized here.

## 2.1.

### Subjects

In this study, we enrolled 11 healthy, right-handed subjects (8 male, 3 female). The subjects were 20 to 60 years old (mean $31\pm 10.5$ ). The Massachusetts General Hospital Institutional Review Board approved the study and all subjects gave written informed consent. For BOLD:DOI comparisons (study I), a total of six subjects were scanned. One subject failed to show any significant activation $(p<0.05)$ in either modality and was excluded from further data analysis. All of the remaining five subjects (four male, one female) showed significant localized activation patterns in both modalities. In the ASL-DOI experiments (study II), again a total of six subjects were scanned. One of these subjects had to be excluded because of a low signal-to-noise ratio in the DOI measurement, due to poor coupling of the probe to the head. The remaining five subjects (four male, one female) were included in the full analysis. One subject was repeated for both studies (subject B/F).

## 2.2.

### Protocol

Prior to placement in the bore of the MRI scanner, subjects were briefed on the details of the finger-tapping task. During data acquisition, they were instructed to watch a visual display from a laptop computer, which was projected from the back of the scanner room onto a screen in the bore of the magnet and was visible from a mirror mounted on the top of the MR head coil above the subjects’ eyes. On this display was shown the text “STOP” or “GO.” For the duration of the GO visual cue, subjects were instructed to tap their thumb and each of the fingers on their right hand sequentially at a self-paced rate (approximately 2 to $3\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$ ). After $2\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ , the GO text was changed to the STOP command; at which time, subjects were instructed to stop tapping and keep their hand relaxed. Aside from the finger-tapping, subjects were otherwise instructed to remain motionless for the duration of the scanning session.

During the scan, the interstimulus interval (ISI) between finger-tapping periods was pseudorandomly chosen and optimized to provide the uniform temporal coverage necessary for deconvolution with a 500-ms time step.^{28} The length of the ISI ranged from 4 to
$20\phantom{\rule{0.3em}{0ex}}\mathrm{s}$
with an average ISI period of
$12\phantom{\rule{0.3em}{0ex}}\mathrm{s}$
. The timing of the stimulus presentation was synchronized with the MR image acquisitions and generated with a custom-written MATLAB script (Mathworks, Sherborn, Massachusetts). Each run lasted
$6\phantom{\rule{0.3em}{0ex}}\mathrm{min}$
and consisted of between 27 and 32 stimulus periods. This was repeated 4 to 6 times for each subject during the course of 1 scan session (approximately
$90\phantom{\rule{0.3em}{0ex}}\mathrm{min}$
) for both the ASL and BOLD measurements. The experimental procedure, subject instructions, and stimulus parameters were identical for both the BOLD and ASL functional scans.

## 2.3.

### DOI Acquisition and Processing

In these experiments, we used a multichannel continuous-wave optical imager (CW4, TechEn Inc., Milford, Masschusetts) to obtain the measurements as previously described.^{29} The DOI imager has 18 lasers—9 lasers at
$690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
$\left(18\phantom{\rule{0.3em}{0ex}}\mathrm{mW}\right)$
and 9 at
$830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
$\left(7\phantom{\rule{0.3em}{0ex}}\mathrm{mW}\right)$
—and 16 detectors. Of these, only four source positions (two wavelengths each) and eight detectors were used here. Laser wavelengths were chosen to minimize cross talk between the two hemoglobin species at
$690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
and
$830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
(
$18\phantom{\rule{0.3em}{0ex}}\mathrm{mW}$
and
$7\phantom{\rule{0.3em}{0ex}}\mathrm{mW}$
, respectively).^{24, 25, 26}

For detectors, the imager employs 16 avalanche photodiodes [(APDs), Hamamatsu C5460-01), each of which is digitized at approximately
$40\phantom{\rule{0.3em}{0ex}}\mathrm{kHz}$
. A bandpass filter follows each APD module with a cut-on frequency of approximately
$500\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$
to reduce
$1\u2215f$
noise and the 60-Hz room light signal, and a cutoff frequency of approximately
$16\phantom{\rule{0.3em}{0ex}}\mathrm{kHz}$
to reduce the third harmonics of the square-wave signals. The signal is then passed to a programmable gain stage, which is used to match the signal levels with the acquisition level on the analog-to-digital converter within the computer. A digital demodulation and low-pass filter is used off-line to separate the individual source signals on each detector channel, which were frequency encoded as described in Franceschini
^{29} This low-pass filtering was done via a seven-pole infinite-impulse response filter Butterworth design with a 10-Hz passband.

The DOI probe is made of flexible plastic strips with plastic caps inserted in it to hold the ends of the 10-m source-detector fiber optic bundles. The probe consists of two rows of four detector fibers and one row of four source fibers arranged in a rectangular grid pattern and spaced
$2.9\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$
between nearest neighbor source-detector pairs as demonstrated in Fig. 1
. This separation distance was chosen to give enough depth penetration to provide sensitivity to the brain.^{30} The plastic probe was secured to the subject’s head via Velcro straps and foam padding and centered over the contralateral primary motor cortex (M1). The 10-m fibers were run through the magnet bore to the back of the scanner and through a port into the control room where they were connected to the DOI instrument. The CW4 instrument was operated from the MR control room.

After collection, demodulation, and low-pass filtering at a 10-Hz bandwidth to separate individual source contributions, the data was further processed using a custom MATLAB data analysis program (HOMER).^{31} Signals were further low-pass filtered at
$0.8\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$
using a zero-phase forward and reverse digital filter to remove the heart signal. Changes in optical density for each source-detector pair were then high-pass filtered with a 1/30-Hz drift correction and finally converted to change in hemoglobin concentrations using the modified Beer-Lambert relationship^{22, 23, 32} with a differential path-length correction calculated from the Monte Carlo simulations (Sec.
2.5).

## 6.

## Eq. 6a

$$\Delta \left[HbR\right]=\frac{{\mathbf{\epsilon}}_{\mathrm{HbO}2}\left({\lambda}_{2}\right)\bullet \mathrm{\Delta}OD\left({\lambda}_{1}\right)\u2215[L\bullet {l}_{DPF}\left({\lambda}_{1}\right)]-{\mathbf{\epsilon}}_{\mathrm{HbO}2}\left({\lambda}_{1}\right)\bullet \mathrm{\Delta}OD\left({\lambda}_{2}\right)\u2215[L\bullet {l}_{DPF}\left({\lambda}_{2}\right)]}{{\mathbf{\epsilon}}_{HbR}\left({\lambda}_{1}\right)\bullet {\mathbf{\epsilon}}_{\mathrm{HbO}2}\left({\lambda}_{2}\right)-{\mathbf{\epsilon}}_{\mathrm{HbO}2}\left({\lambda}_{1}\right)\bullet {\mathbf{\epsilon}}_{HbR}\left({\lambda}_{2}\right)},$$## Eq. 6b

$$\Delta \left[{\mathrm{HbO}}_{2}\right]=\frac{{\mathbf{\epsilon}}_{HbR}\left({\lambda}_{1}\right)\bullet \mathrm{\Delta}OD\left({\lambda}_{2}\right)\u2215[L\bullet {l}_{DPF}\left({\lambda}_{2}\right)]-{\mathbf{\epsilon}}_{HbR}\left({\lambda}_{2}\right)\bullet \mathrm{\Delta}OD\left({\lambda}_{1}\right)\u2215[L\bullet {l}_{DPF}\left({\lambda}_{1}\right)]}{{\mathbf{\epsilon}}_{HbR}\left({\lambda}_{1}\right)\bullet {\mathbf{\epsilon}}_{\mathrm{HbO}2}\left({\lambda}_{2}\right)-{\mathbf{\epsilon}}_{\mathrm{HbO}2}\left({\lambda}_{1}\right)\bullet {\mathbf{\epsilon}}_{HbR}\left({\lambda}_{2}\right)}.$$From the concentration data, the individual subject hemodynamic responses were calculated using an ordinary least-squares (OLS) linear deconvolution and implemented within the HOMER program (e.g., Ref. 33). No nuisance variable polynomial trends were fit in the linear regression, since the low- and high-pass filtering steps had already removed these. In contrast to the generalized linear modeling (GLM) methods commonly used in fMRI analysis, which make use of assumed functional forms of the hemodynamic response [i.e., SPM (Frackowiak
^{34})], OLS deconvolution estimates the entire response function and thus does not assume any fixed canonic responses (i.e., Gaussian or gamma variants). Because both the fMRI and DOI time-series analyses are preformed in this way, this allows comparison of the response functions from both modalities without concerns of introducing errors that might have arisen from such assumed response functions. Epoch timing was synchronized to the MRI images so that subject response times were kept constant between the fMRI and DOI. Regions of time showing significant motion artifacts (as clearly evidenced by extremely large and sudden perturbations in the measurement time course) were rejected from the analysis. The impulse response functions were then averaged across all runs for each subject taken within the same scan session for each source-detector pair.

## 2.4.

### fMRI Acquisition and Processing

For study I, BOLD-fMRI measurements were preformed using a
$3\phantom{\rule{0.3em}{0ex}}\mathrm{T}$
Siemens Allegra MR scanner (Siemens Medical Systems, Erlangen, Germany). fMRI data were collected using a gradient echo-planar imaging (EPI) sequence
$[\mathrm{TR}\u2215\mathrm{TE}\u2215\alpha =500\phantom{\rule{0.3em}{0ex}}\mathrm{ms}\u221530\phantom{\rule{0.3em}{0ex}}\mathrm{ms}\u221590\phantom{\rule{0.3em}{0ex}}\mathrm{deg}]$
with five 6-mm oblique orientation slices (1-mm spacing) and 3.75-mm in-plane spatial resolution. Structural scans were performed using a T1-weighted magnetization-prepared rapid acquisition gradient echo (MPRAGE) sequence (
$1\times 1\times 1.33\text{-}\mathrm{mm}$
resolution,
$\mathrm{TR}\u2215\mathrm{TI}\u2215\mathrm{TE}\u2215\alpha =2530\phantom{\rule{0.3em}{0ex}}\mathrm{ms}\u22151100\phantom{\rule{0.3em}{0ex}}\mathrm{ms}\u22153.25\phantom{\rule{0.3em}{0ex}}\mathrm{ms}\u22157\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$
]. To calculate the BOLD-based hemodynamic response functions, the functional images were first motion-corrected^{35} and spatially smoothed with a 6-mm Gaussian kernel. The response functions were then calculated by an OLS linear deconvolution. A third-order polynomial was included to remove drift effects. For each subject, the effect and standard deviation maps were input to a mixed effects analysis to generate a map of
$t$
statistics.^{36} The mean of the response from 2 to
$7\phantom{\rule{0.3em}{0ex}}\mathrm{s}$
was used to calculate the effect size.

For study II, ASL-fMRI was carried out at
$3\phantom{\rule{0.3em}{0ex}}\mathrm{T}$
(same scanner as study I) using proximal inversion with control for off-resonance effects (PICORE) labeling geometry^{37} with Q2TIPS saturation^{38} to impose a controlled label duration. A postlabel delay of
$1400\phantom{\rule{0.3em}{0ex}}\mathrm{ms}$
and a label duration of
$700\phantom{\rule{0.3em}{0ex}}\mathrm{ms}$
were used, with repetition and echo times of
$2\phantom{\rule{0.3em}{0ex}}\mathrm{s}$
and
$20\phantom{\rule{0.3em}{0ex}}\mathrm{ms}$
, respectively
$[\alpha =90\phantom{\rule{0.3em}{0ex}}\mathrm{deg}]$
. The PICORE labeling scheme allowed the collection of BOLD signals using the control phase of the acquisition. EPI was used to image five 6-mm slices (1-mm spacing) with 3.75-mm in-plane spatial resolution. Structural scans were also performed with the same scan prescriptions as study I. Although the original image acquisition was at
$2\phantom{\rule{0.3em}{0ex}}\mathrm{s}$
, the stimulus had been jittered evenly on a 500-ms time step, which allowed for the response to be calculated at
$2\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$
, albeit with a lower signal-to-noise ratio. To estimate the flow response, the ASL-functional scans were first separated into the control and negatively labeled tag images. These traces were then independently interpolated to upsample the data points to
$2\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$
using a cubic spline model.^{39} Subtracting each negatively labeled tag scan from the immediately subsequent control scan generated the flow image series. This flow series was then deconvolved and the region of interest average calculated similar to its calculation for the BOLD response. The mean of the response from 1 to
$6\phantom{\rule{0.3em}{0ex}}\mathrm{s}$
was used to calculate the effect size. Figure 2
depicts the processing stream.

## 2.5.

### Monte Carlo Simulations

The anatomical MPRAGE MRI scans taken at the start of each experimental session were resampled to 1-mm (isotropic) resolution and then segmented into white and gray matter using anatomical reconstruction tools available as part of the FREESURFER package from the Massachusetts General Hospital (MGH).^{40} This process first consisted of intensity normalization, skull stripping, and white-gray matter differentiation as described by Dale
^{41} and implemented in the RECON-ALL program as part of the MGH FREESURFER package. All parameters were set to their default values. To further differentiate the anatomical images into skin, skull, and cerebral spinal fluid (CSF) layers, a similar additional intensity-based initialization followed by sorting of ambiguous voxels based on the iterative expectation-maximization algorithm was applied to the unassigned regions of the intensity normalized image that were not accounted for after the skull-stripping algorithm.^{42, 43} For this, a custom written routine was used, which required manual supervision of this segmentation process. After full segmentation, tissue assignments were manually inspected. This process was repeated individually for all the subjects.

To simulate photon migration through the segmented brains, the probe positions were first determined from the location of the vitamin E fiducial markers that had been placed on the DOI probe in the MPRAGE scans. After locating the probe positioning, Monte Carlo–based simulations of light propagation were performed as described in Boas
^{16} The optical properties of each of the segmented tissues that were used in these simulations are listed in Table 1
.^{24} In these simulations, the optical properties of gray and white matter were taken to be the composite properties for brain tissue, which is optically dominated by gray matter. For each optode position, the trajectories of
${10}^{8}$
photons were simulated at both the 690-nm and 830-nm wavelengths. Because the amount of detected light decays exponentially with distance, only the three-point Green’s functions for the nearest-neighbor source-detectors were calculated. These simulations were also used to determine the differential path-length corrections for each wavelength and source-detector pair and applied to the calculation of concentration changes in the modified Beer-Lambert law used in Eqs. 5, 6.

## Table 1

Optical absorption and scattering properties of the tissue types used in the Monte Carlo simulations. Monte Carlo simulations of light propagation were performed on the tissue segmented structural images for each subject using the tissue optical properties at both wavelengths recorded (830 and 690nm ). Optical values for the various tissue types were taken from Ref. 24.

830mn | 690nm | |||
---|---|---|---|---|

Tissue | μa (mm−1) | μs′ (mm−1) | μa (mm−1) | μs′ (mm−1) |

Skin | 0.0191 | 0.66 | 0.0159 | 0.80 |

Skull | 0.0136 | 0.86 | 0.0101 | 1.00 |

CSF | 0.0026 | 0.01 | 0.0004 | 0.01 |

Brain (Gray and White) | 0.0186 | 1.11 | 0.0178 | 1.25 |

## 2.6.

### fMRI Projection

In order to test the spatial correspondence of fMRI activation with the DOI measurements, the fMRI hemodynamic responses were multiplied by the optical forward matrix $\left(\mathsf{A}\right)$ to project the changes into source-detector space as described in Sec. 2.5. This was done by first registering and resampling the forward matrices from the 1-mm isotropic resolution of the anatomical images and Monte Carlo simulations to the lower spatial resolution of the functional scans. The forward matrix is then multiplied onto the data matrix holding the fMRI hemodynamic time courses from the BOLD or ASL scans. In order to reduce the noise contribution from the fMRI, the activation $t$ maps from the fMRI were used to first mask the fMRI data $(p<0.1)$ prior to projection. Although, only voxels under this criterion were used, this high choice of $p$ value should be considered quite lenient. The optical probe extended beyond the spatial coverage of the fMRI scan, which spanned roughly only $3.5\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ axially [5 slices $\times $ (6mm $+$ 1-mm gap $=7\left)\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\right]$ . In order to avoid the introduction of partial volume errors, the projection and analysis were limited to only those source-detector pairs for which the forward operator was spanned by the functional MRI sample volume (i.e., the superior row of source and detector pairs).

In Fig. 3 , we illustrate the concept of the fMRI projection by showing the overlay of the three-point Green’s function (shown as contour lines in logarithmic scale) and the anatomical MRI volume. The functional activation region (shown only above half-maximum amplitude) is inset, showing peak activation approximately $21\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ below this particular source-detector pair. The forward equation sums these changes over the volume using the weights of the sensitivity matrix (Green’s function) to predict the measurements between each source-detector pair.

Because the BOLD signal is proportional to deoxyhemoglobin (but this proportionality is a matter of debate),^{12, 44, 45} the resulting deoxyhemoglobin predictions need to be normalized. A single scaling factor, to normalize to unity the maximum change across the entire probe, was applied for each subject.

Although we report the results from using a spectrally defined forward operator [Eqs. 3, 4], we also examined the results using the Green’s functions at both 690 and $830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ . Due to the coarse spatial resolution of the fMRI relative to any wavelength-specific variations noted in the Green’s functions, the results were nearly identical with all three choices of forward operators. For consistency, we report only the spectral operator.

## 2.7.

### Statistical Analysis

To compare the optical measurements and their predicted values from the fMRI, we used linear regression analysis to report correlation values for the spatial or temporal amplitude profiles. Pearson coefficients (i.e.,
$R$
values) describe the zeroth-order correlation between two vectors of values. The associated
$p$
value for this test describes the probability that random chance alone could produce a correlation equal or better than this
$R$
value.^{46}

Because Pearson coefficients are nonlinear, to compare the correlations between hemodynamic parameters, we use a $z$ transformation of the $R$ value. The $z$ transformation is defined as

Statistical differences between two different comparisons (i.e., correlations) can be tested by a paired $t$ test of their $z$ transformed $R$ values.^{46}

## 2.8.

### Depth Sensitivity of Optical Imaging

The differential path-length factors (DPF) for each source-detector pairing were calculated from the Monte Carlo data from each subject. The DPF for a given source-detector pair was calculated as the actual mean path length traveled by all photons leaving the source and reaching the detector, divided by the linear source-detector separation.^{22, 23, 24}

In addition to allowing us to forward model the fMRI signal into the optical imaging space, the determination of the three-point Green’s functions in this work allows us to comment on the depth sensitivity of the optical measurements. We calculated the mean activation depth of the BOLD and ASL responses for each subject by computing the first moment of distance to the nearest point on the head’s surface. The mean depth of the optical contrast was calculated in a similar way by weighting by the appropriate Green’s functions $\left(\mathsf{A}\right)$ . That is,

## 8.

## Eq. 8a

$${\u27e8depth\u27e9}_{fMRI}=\frac{\sum _{\mathit{r}}\parallel {\mathit{r}}_{surface}-\mathit{r}\parallel \bullet fMR{I}_{\mathit{r}}}{\sum _{\mathit{r}}fMR{I}_{\mathit{r}}}.$$## Eq. 8b

$${\u27e8depth\u27e9}_{DOI}=\frac{\sum _{i=\left\{src\text{-}det\right\}}\sum _{\mathit{r}}\parallel {\mathit{r}}_{surface}-\mathit{r}\parallel \bullet {A}_{\mathit{r},i}\bullet fMR{I}_{\mathit{r}}}{\sum _{i=\left\{src\text{-}det\right\}}\sum _{\mathit{r}}{A}_{\mathit{r},i}\bullet fMR{I}_{\mathit{r}}}.$$## 3.

## Results

## 3.1.

### Study I: BOLD/DOI

In study I, we performed simultaneous optical and fMRI-BOLD recordings. Individual subject measurements were gathered and processed as described in the previous sections.

In Fig. 4 , we show the results from the projection of the BOLD signal through the optical forward equation, as described in Eq. 3, for each of the five subjects. The time course for each (nearest-neighbor) source-detector pair is shown for both the DOI HbR recording and the projected BOLD trace. Both of these time courses have been normalized to their respective maximal time points across all channels, such that the relative amplitudes between each pair are comparable. Both time courses are shown as positive changes for direct comparison.

In Fig. 4, we note that there is generally good spatial and temporal agreement between the BOLD and optical signals. In several of the subjects, we see the characteristics of poststimulus undershoots in both the BOLD and HbR data. Among the channels that did not significantly $(p>0.05)$ respond to the stimulus, we still see some degree of correlation in the time courses, suggesting that the noise in these channels has a physiological source common to both modalities.

In order to further investigate the spatial correlation between the optical and BOLD measurements, we calculated the maximum response amplitude for each source-detector pair. This is presented as a bar graph in Fig. 5 . We also calculated the mean response over the time period from 2 to $7\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ for the HbR and BOLD responses, yielding results similar to those found using the maximum of the response (data not shown). The comparison of response amplitudes presented in Fig. 5 shows more clearly the spatial correspondence of the BOLD and optical signals. The distribution of response amplitudes across the probe is similar in both the measured optical data and predictions from the BOLD signal. In addition, we also performed a similar analysis to compare the BOLD and optically measured total hemoglobin (HbT) and ${\mathrm{HbO}}_{2}$ responses. The results of the spatial correlation of BOLD with the optically measured hemoglobin parameters are given in Table 2 .

## Table 2

Spatial correlation (study I) results show the spatial correlation of BOLD to the optically measured HbR, HbO2 , and HbT parameters. The BOLD signal was projected through the HbR forward model (shown in Fig. 2) and then the mean amplitude for each channel was calculated (shown in Fig. 3). These amplitudes were then correlated between the two modalities. Displayed here are the correlation (Pearson’s) coefficients for each comparison for the five subjects used in study I. The probability ( p value) of each correlation is shown in brackets for each comparison.

Subject | HbR:BOLD | HbO2:BOLD | HbT:BOLD |
---|---|---|---|

A | 0.95 $(1\times {10}^{-3})$ | −0.94 $(2\times {10}^{-3})$ | −0.85 $(1\times {10}^{-2})$ |

B | 0.80 $(2\times {10}^{-2})$ | 0.59 $(1\times {10}^{-1})$ | 0.55 $(2\times {10}^{-1})$ |

C | 0.64 $(1\times {10}^{-1})$ | 0.72 $(6\times {10}^{-2})$ | 0.82 $(2\times {10}^{-2})$ |

D | 0.84 $(1\times {10}^{-2})$ | 0.72 $(7\times {10}^{-2})$ | 0.60 $(2\times {10}^{-1})$ |

E | 0.62 $(1\times {10}^{-1})$ | 0.40 $(3\times {10}^{-1})$ | 0.51 $(2\times {10}^{-1})$ |

All | 0.68 $(6\times {10}^{-6})$ | 0.38 $(4\times {10}^{-2})$ | 0.37 $(5\times {10}^{-2})$ |

The spatial correlation between the BOLD and optically measured HbR amplitudes was strongly significant ( $R=0.68$ , $p=6\times {10}^{-6}$ ) across all five subjects. Additionally, the BOLD/HbR spatial correlation was better than that observed between $\mathrm{BOLD}\u2215{\mathrm{HbO}}_{2}$ $(R=0.38)$ or BOLD/HbT $(R=0.37)$ . A paired $t$ test on these values allows us to test the probability that the correlation between BOLD/HbR is statistically better than that between $\mathrm{BOLD}\u2215{\mathrm{HbO}}_{2}$ $(p=0.12)$ and BOLD/HbT $(p=0.12)$ ( $Z$ -transformed $t$ test).

## 3.2.

### Study II: ASL:DOI

In the second study, we repeated the same motor task experiment and simultaneously recorded ASL-based fMRI and optical signals. These measurements were analyzed similarly as in study I. In Fig. 6 , we present in optical measurement space, the time courses for the fMRI-measured blood flow, which has been projected from the ASL-fMRI data by the simulated optical forward matrix $\mathsf{A}$ [Eq. 3]. The amplitudes have been scaled to match the single maximum amplitude of the fMRI and optical time courses, allowing for a comparison of the relative amplitudes across space and time. Here we see a good correspondence between those source-detector pairs predicted to have the largest response from the fMRI measurements and the optical data.

In Fig. 7 , we present the maximum response amplitudes across all source-detector pairs to clarify the spatial correlation between the fMRI and optical measurements. The calculated correlations are presented in Table 3 .

## Table 3

Spatial correlation (study II) results show the spatial correlation of ASL to the optically measured HbR, HbO2 , and HbT parameters. The ASL signal was projected through the HbO2 forward model (shown in Fig. 4) and then the mean amplitude for each channel was calculated (shown in Fig. 5). These amplitudes were then correlated between the two modalities. Displayed here are the correlation (Pearson’s) coefficients for each comparison for the five subjects used in study II. The probability ( p value) of each correlation is shown in brackets for each comparison.

Subject | HbR:ASL | HbO2:ASL | HbT:ASL |
---|---|---|---|

F | 0.41 $(4\times {10}^{-1})$ | 0.93 $(6\times {10}^{-3})$ | 0.95 $(3\times {10}^{-3})$ |

G | −0.44 $(3\times {10}^{-1})$ | 0.44 $(3\times {10}^{-1})$ | 0.50 $(3\times {10}^{-1})$ |

H | −0.78 $(7\times {10}^{-2})$ | 0.92 $(1\times {10}^{-2})$ | 0.92 $(1\times {10}^{-2})$ |

I | −0.93 $(6\times {10}^{-3})$ | 0.87 $(2\times {10}^{-2})$ | 0.89 $(2\times {10}^{-2})$ |

J | −0.74 $(6\times {10}^{-1})$ | 0.59 $(2\times {10}^{-1})$ | 0.53 $(2\times {10}^{-1})$ |

All | 0.26 $(1\times {10}^{-1})$ | 0.73 $(3\times {10}^{-6})$ | 0.75 $(1\times {10}^{-6})$ |

The spatial correlation between the ASL and optically measured ${\mathrm{HbO}}_{2}$ or HbT amplitudes is strongly significant ( $R=0.73$ , $p=3\times {10}^{-6}$ and $R=0.75$ , $p=1\times {10}^{-6}$ , respectively) across all five subjects. The ASL/HbT and $\mathrm{ASL}\u2215{\mathrm{HbO}}_{2}$ spatial correlations are better than that observed between ASL/HbR ( $R=0.26$ , $p=0.10$ ). These correlations between $\mathrm{ASL}\u2215{\mathrm{HbO}}_{2}$ and ASL/HbT were determined to be nearly significantly ( $p=0.08$ and $p=0.08$ ) better than the correlation between ASL/HbR ( $Z$ -transform $t$ -test).

Since the pulsed ASL scan used in this study allows the extraction of BOLD time courses from the interweaved control images, we were able to repeat the analysis to examine the BOLD/HbR correlation from this set of subjects. Although not shown, these results are similar to those from study I, which used a dedicated BOLD scan. We found that the spatial correlation between BOLD and HbR was 0.43 $(p=0.07)$ . This result is slightly weaker than that reported in study I, due to the lower contrast-to-noise ratio of the BOLD signal inherent in the $p\mathrm{ASL}$ scan sequence $(\mathrm{TE}=20\phantom{\rule{0.3em}{0ex}}\mathrm{ms})$ . However, similar to the results presented in Table 2, the correlation between BOLD and HbR is higher than that between $\mathrm{BOLD}\u2215{\mathrm{HbO}}_{2}$ $(p=0.01)$ and BOLD/HbT $(p=0.02)$ ( $Z$ -transform $t$ -test).

## 3.3.

### Comparison of BOLD versus ASL-based fMRI

Our results from Secs. 3.1, 3.2 suggest that the spatial profiles of HbR and ${\mathrm{HbO}}_{2}$ are distinct, since these measurements correlated better to the BOLD and ASL activations, respectively. To further investigate this result, we examined the spatial comparison between BOLD and ASL directly. Because the pulsed ASL sequence used in study II interweaves unlabeled, control images among the inversion labeled flow images, this allowed the extraction of both BOLD and flow-weighted time courses. Statistical maps of the activations from both of these time-courses were calculated.

As reviewed by Detre and Wang,^{47} ASL and BOLD imaging have different spatiotemporal characteristics. These two fMRI methods, in particular BOLD imaging, depend on details of the vascular anatomy, which in turn affect the measurements and the resulting images. It was therefore not surprising that we observed a slight disparity between the two activation maps from these two fMRI modalities; this is consistent with similar observations by Silva, ^{48} Luh, ^{49} and Yang
^{50, 51} Shown in Fig. 8
are the ASL and BOLD activation maps for one such subject (F).

## Table 4

Spatial comparison of ASL and BOLD activation profiles. For each of these five subjects, we calculated the center of mass of the BOLD and ASL activations. This was done by calculating the positional moment weighted by the normalized p -value map as described in the text. The second radial moments were also calculated to describe the area of the activation.

ASL | BOLD | ||||
---|---|---|---|---|---|

Subject | Activation Center (First Moment) | Activation Area (Second Moment) | Activation Center (First Moment) | Activation Area (Second Moment) | Displacement |

F | [127, 56, 142] | $7351\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{2}$ | [132, 54, 139] | $6343\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{2}$ | $6.5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ |

G | [128, 52, 132] | $8067\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{2}$ | [130, 47, 136] | $6113\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{2}$ | $6.0\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ |

H | [128, 67, 140] | $8395\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{2}$ | [128, 65, 141] | $6656\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{2}$ | $2.6\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ |

I | [129, 52, 137] | $8373\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{2}$ | [132, 48, 144] | $5881\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{2}$ | $9.1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ |

J | [125, 35, 152] | $5913\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{2}$ | [128, 35, 142] | $6865\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{2}$ | $10.7\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ |

All | — | — | — | — | 7.0 $(\pm 3.1)\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ |

For each of these five subjects we calculated the “center of mass” of the BOLD and ASL activations (shown in Table 4 ). This was done by calculating the positional moment weighted by the normalized $p$ -value map (i.e., ${\sum}_{\mathit{r}}\mathit{r}\bullet {p}_{\mathit{r}})$ . The second radial moments were also calculated to describe the area of the activation $\left[{\sum}_{\mathit{r}}\right(\mathit{r}-\u27e8\mathit{r}\u27e9{)}^{2}\bullet {p}_{\mathit{r}}]$ . As shown in Table 4, we observed a displacement between the BOLD and ASL centers of activation (mean $7.0\pm 3.1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ ).

## 3.4.

### Temporal Correlation between DOI and fMRI

In our previous work,^{12} we reported a strong temporal correlation between the region-of-interest averaged time courses of the optically measured HbR and BOLD signals and between the ASL and HbT signals. However, in that analysis, we did not take into account the effects of DOI sensitivity or partial volume effects across the probe. Previously, these region-of-interest averages had been calculated from the inclusion of all source-detector pairs or fMRI voxels that met a
$p<0.05$
criterion. In this current approach, by using the DOI forward matrix to project a prediction of each source-detector measurement from the BOLD signal, we have the ability to better address these effects. The projection of the fMRI onto the optical probe using the optical forward model accounts for these partial volume effects across the probe. One aspect of this current analysis, however, is that the resulting time courses have a smaller contrast-to-noise ratio because the projection preferentially weights the fMRI response from a small number of superficial voxels.

The mean temporal response from all source-detector pairs from the top half of the DOI probe for all subjects was used to calculate the response for each modality. The resultant average time courses are shown in Fig. 9
. Consistent with the results from our previous work,^{12} the results presented in Table 5
show statistically strong temporal correlations between the BOLD and optically measured HbR responses (
$R=0.93$
,
$p=1\times {10}^{-13}$
) and again between the ASL and HbT responses (
$R=0.91$
,
$p=3\times {10}^{-12}$
). The BOLD/HbR correlation tends to be higher than that observed between the
$\mathrm{BOLD}\u2215{\mathrm{HbO}}_{2}$
or BOLD/HbT. This trend, however, is not significant (
$p=0.23$
and
$p=0.20$
, respectively) (
$Z$
-transform
$t$
-test). Similarly, the ASL and HbT correlation is stronger than that between ASL and HbR. This result was nearly significant in a test of the paired
$Z$
-transformed Pearson’s coefficients
$(p=0.07)$
. The correlation between ASL and HbT is also higher than that between
$\mathrm{ASL}\text{-}{\mathrm{HbO}}_{2}$
, although, this trend is again not significant
$(p=0.14)$
(
$Z$
-transform
$t$
-test).

## Table 5

Temporal correlation results show the temporal correlation of BOLD and ASL to the optically measured HbR, HbO2 , and HbT parameters. The time courses from the projected fMRI and DOI were averaged over all source-detector channels and then correlated. Displayed are the correlation (Pearson’s) coefficients for each comparison for the five subjects used in study II. The probability ( p value) of each correlation is shown in brackets for each comparison.

Subject | HbR:BOLD | HbO2:BOLD | Hbr:BOLD |
---|---|---|---|

A | 0.61 $(4\times {10}^{-4})$ | 0.39 $(3\times {10}^{-2})$ | 0.34 $(6\times {10}^{-2})$ |

B | 0.83 $(4\times {10}^{-8})$ | 0.83 $(2\times {10}^{-8})$ | 0.83 $(2\times {10}^{-8})$ |

C | 0.68 $(3\times {10}^{-5})$ | 0.78 $(3\times {10}^{-7})$ | 0.82 $(2\times {10}^{-8})$ |

D | 0.65 $(1\times {10}^{-4})$ | 0.60 $(5\times {10}^{-4})$ | 0.48 $(7\times {10}^{-3})$ |

E | −0.08 $(6\times {10}^{-1})$ | 0.02 $(9\times {10}^{-1})$ | 0.07 $(3\times {10}^{-1})$ |

All | 0.93 $(1\times {10}^{-13})$ | 0.89 $(4\times {10}^{-11})$ | 0.88 $(1\times {10}^{-10})$ |

HbR:ASL | ${\mathrm{HbO}}_{2}:\mathrm{ASL}$ | HbT:ASL | |

F | 0.11 $(2\times {10}^{-1})$ | 0.69 $(6\times {10}^{-19})$ | 0.70 $(4\times {10}^{-20})$ |

G | 0.42 $(7\times {10}^{-7})$ | 0.29 $(8\times {10}^{-4})$ | 0.37 $(2\times {10}^{-5})$ |

H | −0.49 $(3\times {10}^{-9})$ | 0.51 $(1\times {10}^{-9})$ | 0.58 $(5\times {10}^{-13})$ |

I | 0.07 $(4\times {10}^{-1})$ | 0.77 $(7\times {10}^{-27})$ | 0.62 $(1\times {10}^{-14})$ |

J | 0.27 $(2\times {10}^{-3})$ | 0.32 $(2\times {10}^{-4})$ | 0.25 $(4\times {10}^{-3})$ |

All | 0.59 $(4\times {10}^{-4})$ | 0.80 $(9\times {10}^{-8})$ | 0.91 $(3\times {10}^{-12})$ |

## 3.5.

### Depth Sensitivity of Optical Imaging

We calculated the mean depth of contrast for the fMRI and DOI measurements as defined by Eq. 8. These values are reported in Table 6 . We found that the contrast depth measured by the fMRI was $20\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ $(\pm 5\phantom{\rule{0.3em}{0ex}}\mathrm{mm})$ for the BOLD and $17\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ $(\pm 5\phantom{\rule{0.3em}{0ex}}\mathrm{mm})$ for the ASL. The mean depth of the optical contrast was shallower at only $14\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ $(\pm 4\phantom{\rule{0.3em}{0ex}}\mathrm{mm})$ for both the ASL- and BOLD-derived signals. This depth was approximately equal to the cortical depth (mean distance from the scalp to the cortex), which was $13\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ $(\pm 2\phantom{\rule{0.3em}{0ex}}\mathrm{mm})$ . The thicknesses of the superficial layers (skin, skull, and CSF) are also shown in Table 6. The skull thickness and CSF layer varied considerably over the nine different anatomical volumes (B and F are the same subject). In some cases, the mean optical contrast was shallower than the cortical depth. This is due to the lower resolution of the functional MR images $\left(3.75\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\right)$ compared to the anatomical volume $\left(1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\right)$ .

## Table 6

Contrast depth and sensitivity comparison between the depth of functional contrast from the surface of the head for the fMRI and optical measurements. These were calculated as described in the text in Eqs. 8a, 8b. This reflects the mean depth of the signal giving rise to each measurement. For comparison, we also show the cortical depth, defined as the mean distance of the cortex from the surface of the head. The thickness of the three superficial layers is given in millimeters. The mean DPF for each subject is shown with the standard deviation of the value over all source-detector pairs.

Cortical Depth | Functional Depth | Photon Migration | |||||
---|---|---|---|---|---|---|---|

Subject | Total [mm] | [Skin, Skull, CSF] thickness [mm] | fMRI Contrast [mm] | Optical Contrast [mm] | Differential Path-Length Factor | ||

A | 12.2 | [2.6, 3.9, 5.7] | BOLD | 24.9 | 13.1 | $830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.46 (0.38) |

$690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.48 (0.38) | ||||||

B | 14.4 | [3.1, 4.6, 6.7] | BOLD | 14.9 | 13.1 | $830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.45 (0.38) |

$690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.48 (0.41) | ||||||

C | 12.6 | [2.7, 2.7, 7.2] | BOLD | 22.4 | 21.0 | $830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.63 (0.25) |

$690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.73 (0.23) | ||||||

D | 11.7 | [2.5, 3.6, 5.5] | BOLD | 20.7 | 7.6 | $830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.00 (0.30) |

$690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.11 (0.36) | ||||||

E | 12.8 | [2.7, 5.9, 4.2] | BOLD | 26.8 | 14.6 | $830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 4.79 (0.36) |

$690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 4.84 (0.36) | ||||||

F | 14.4 | [3.1, 4.6, 6.7] | BOLD | 14.2 | 14.6 | $830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.12 (0.33) |

ASL | 13.8 | 13.1 | $690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.15 (0.32) | |||

G | 16.9 | [3.1, 6.7, 7.0] | BOLD | 19.5 | 19.3 | $830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.53 (0.90) |

ASL | 17.4 | 19.4 | $690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.59 (0.91) | |||

H | 13.6 | [2.6, 3.1, 7.9] | BOLD | 21.6 | 16.5 | $830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.59 (0.23) |

ASL | 21.3 | 17.5 | $690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.61 (0.24) | |||

I | 11.8 | [1.8, 4.2, 5.9] | BOLD | 12.2 | 9.9 | $830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.93 (1.11) |

ASL | 11.4 | 10.3 | $690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.93 (1.09) | |||

J | 10.8 | [2.5, 3.2, 5.2] | BOLD | 22.9 | 10.8 | $830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 4.99 (0.47) |

ASL | 22.9 | 11.3 | $690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 4.96 (0.47) | |||

AII | 13.0 | [2.6, 4.2, 6.1] | BOLD | 20.0 | 14.0 | $830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.35 (0.35) |

ASL | 17.4 | 14.3 | $690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ | 5.39 (0.36) |

Related to the mean sensitivity depth, in Table 6, we also show the mean differential path-length factor for each subject, averaged over all source-detector pairs. Consistent with previously published theoretical findings from anatomical head geometries,^{24, 30} we find that the DPF was approximately 5.4 and 5.3
$(\pm 0.4)$
for the 690- and 830-nm wavelengths. For subjects G and I, there was a fairly large variation in the DPF between optode positions, although the ratio of the DPFs for each wavelength showed less variation. From visual inspection, it appeared that the position of the probe on these two subjects was slightly more posterior than on the other subjects. A possible explanation of the higher variation in the DPF for these two subjects could be the variation in skull thickness, increasing toward the back of the head.

## 4.

## Discussion

In this paper, we investigated the spatial relationships between optical and fMRI measures of the hemodynamic response to motor activity. By using Monte Carlo simulations of light propagation, we were able to calculate the sensitivity profiles for each optical source-detector pair on the DOI probe for each subject. Rather than performing image reconstruction of the optical data, we used these sensitivity profiles to forward project the fMRI data enabling us to predict the measurements between each optical channel. This method allows quantitative spatial comparison of optical and fMRI data without using the regularized inversion techniques standard to optical imaging.

## 4.1.

### Spatial Correlation between fMRI and DOI

In the first part of this paper, we compared diffuse optical measurements to simultaneously acquired BOLD-fMRI. The BOLD signal measures localized changes in deoxyhemoglobin (e.g., Buxton, ^{44} Obata
^{45}). However, because MRI derives its signal from the hydrogen nuclei of water molecules, the BOLD signal reflects changes in the resonance lifetime of water rather than a direct measure of deoxyhemoglobin. The ferrous heme group in deoxyhemoglobin has a paramagnetic electron structure. This creates localized magnetic perturbations in and around blood vessels, which give rise to the shortening of the
${t}_{2}$
relaxation time of water molecules as they migrate through these regions. The change in relaxation times
$\left(\mathrm{\Delta}{R}_{2}\right)$
, in turn, gives rise to changes in signal intensity, which underlies the BOLD response. Because it derives from hydrogen nuclei, the BOLD signal contains contributions from both water molecules inside the blood vessels, which can interact directly with hemoglobin, and from water outside the blood vessels, which interact with the magnetic perturbations extending from the vessels into the parenchyma tissue.^{44, 45} The BOLD signal is thus an indirect measurement of HbR, depending on details of the fMRI acquisition sequence, baseline properties of the brain, and vascular sensitivities.^{52} Spatially, the BOLD signal should reflect the underlying vascular anatomy. It is generally believed that functional changes may be biased to larger venous structures.^{47, 53} This observation is the result of two factors. First, deoxyhemoglobin changes are believed to be larger in the veins than in the arterial and capillary compartments. Second, for gradient echo imaging, as used here, the extravascular signal is most sensitive to larger vessels, selectively biasing the BOLD signal to the draining pial veins.^{52} On the other hand, optical measures of HbR will similarly be biased to the venous areas because this is where these changes are the largest. Thus, one would expect that both the optically measured HbR and the BOLD signal would be biased distal to the center of neuronal activity.

The fMRI-ASL contrast is thought to arise mainly from blood flow changes in the arterial and capillary vessels, leading to a signal that is better localized to the cortex, as opposed to the larger vessel structures.^{47, 48} This result, along with the expectation of venous-weighted BOLD, helps explain the slight spatial discrepancies that have been noted between the BOLD and ASL modalities.
^{48, 49, 50, 51}

In this experiment, as was expected based on the differential vascular compartments giving rise to each of the measurements, we found spatial agreement between HbR and BOLD and between the HbT and ASL profiles. In study I, we found that the amplitude of the BOLD signal was slightly better correlated to the optically measured HbR signal than to either the HbT or ${\mathrm{HbO}}_{2}$ signals. In study II, when the complementary ASL and DOI comparisons were made, we found better correlation between the optically measured $\mathrm{HbT}\u2215{\mathrm{HbO}}_{2}$ and ASL signals than between HbR and ASL.

## 4.2.

### Correlation of “Noise” between DOI and fMRI

Although the aim of this study was to compare the functional hemodynamic responses measured with DOI and fMRI, we observed that even some of the noise in the data was present in both the optical and fMRI time courses. This noise, which we consider to be the spurious early or anomalous undershoots in the time courses of the measurements, not typical of the traditional canonical responses of hemodynamic changes, can be seen in both data sets. While, conventionally, many of these features would be considered noise in the response and do not meet statistical significance in effects analysis, the clear presence in both modalities indicates that they most likely have true physiological underpinnings. While perhaps not a component of the functional response, these artifacts are indeed physiological and probably due to blood pressure or respiratory effects. Subject F from study II is an example of this feature showing both early dips in flow and volume and negative going response curves, which are present in both the measured optical and ASL data.

The correlation of these time courses between modalities and varying across regions of the probe indicates that this noise is of physiological origin. Such physiological fluctuations have been noted in previous literature in both fMRI and DOI studies (i.e., Toronov, ^{54} Bhattacharyya and Lowe,^{55} Wise
^{56}). This suggests that accounting for physiological fluctuations may therefore improve future analysis of fMRI and/or optical data, because this is clear evidence of nonwhitened noise contributions to the signals. Because most popular statistical mapping methods used in fMRI or optical analysis often make assumptions of random, white, or signal noise, the presence of such strong contributions suggests more care should be taken to investigate these assumptions. Approaches to modeling these systemic signals have recently been described for analysis of optical data and fMRI and similar methods might be extended to combined multimodality analysis in the future.
^{57, 58, 59, 60}

Furthermore, although optical imaging is very susceptible to superficial, systemic noise effects, such as cardiac or respiration signals, the presence of these correlated artifacts in the fMRI data suggests this noise is not entirely arising from the surface. Although, this superficial, systemic physiology clearly has an impact on the optical data, this alone cannot explain the correlated temporal noise between the fMRI and DOI data, because the fMRI signals are arising from mostly cortical tissues.

## 4.3.

### Generality of Results

Although this experimental design yielded evidence for the correlation of these two methodologies, we must be cautious in drawing conclusions as to the generality of these results to other experimental situations. For both DOI and fMRI, the observed hemodynamic parameters, such as blood volume,
${\mathrm{HbO}}_{2}$
, or HbR, are spatially and temporally filtered versions of their true underlying changes in the brain.^{61, 62} In general, these filters are not expected to be identical between the DOI and fMRI. As a result, the measured results may differ because the filter functions may differ. Because these filters will depend on a number of factors, including the acquisition scheme employed, it is possible that other echo, scan parameters, or even brain regions will yield slightly differing results. At different magnetic field strengths, the extra- or intravascular weightings of the BOLD signal will vary, giving rise to images that are differentially biased to macrovascular structure.^{52, 63} Although this data clearly suggest that there is correspondence between the fMRI and DOI modalities, further investigation of the explicit details of these results should be pursued.

In this paper, we have shown that the ASL signal correlates with the spatial profile of blood volume changes in the brain. This observation is the result of the primarily arteriolar contributions to these two measurements. The correspondence of volume and flow changes will not be strictly obeyed and further investigation will be needed. In particular, the correspondence of ASL with optical measurements of blood flow, such as diffuse correlation spectroscopy,^{64, 65} should be explored in future work.

## 4.4.

### Implications toward Multimodality Integration

Understanding the spatial and temporal correspondence between fMRI and optical methods is an important prerequisite for the multimodality integration of these two techniques and a motivation of the current study. Incorporation of proper measurement models into the processing of multimodality data will enable the ability to exploit analysis schemes that combine the mutual information of all modalities appropriately. In the future, such approaches, designed to maximize the covariance between the two sets of measurements, might greatly improve our discrimination of functional contrast from measurement noise. Placed in a Bayesian formulation, calculating the hemodynamic response that maximizes the *a posteriori* estimate of the joint set of measurements, as opposed to treating the two modalities as independent measurements, could provide a more accurate estimate of the significance of the functional response.^{7, 66} Thus, analyzing these modalities together as measurements derived from a common physiological state variable, rather than as independent variables could allow us to better cancel instrument-associated noise factors. These combined analysis schemes are currently being pursued.

In addition, the spatial correspondence between optical and fMRI techniques may also have an impact on multimodality schemes for optical image reconstruction methods. Although DOI offers several advantages for studying the temporal dynamics of the hemodynamic response, as a tool for imaging or brain mapping, it has a number of limitations as an independent imaging tool. Because DOI is based on the measurement of back-reflected, diffuse light, optical imaging is restricted to superficial depths of 0.5 to
$2\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$
into the brain.^{67} In addition, continuous-wave imaging systems, as used in this study, have virtually no ability to discern activation depth and are exponentially more sensitive to superficial activation. As a result, the varying activation depth across the region of interest creates partial-volume errors. This was the specific motivation for not using image reconstructions in this study.

With the addition of the high-spatial resolution anatomical and functional information that MRI can provide, the optical path-length and partial-volume corrections can be included to improve the quantitative accuracy of DOI (e.g., Intes, ^{5} Boas and Dale^{6}). With functional MRI to provide prior information for the location of activation, image reconstruction schemes can make use of this to constrain solutions.^{7} Because this paper has directly shown that the fMRI data is indeed a possible solution to the inversion problem, we now have more confidence that these functionally constrained methods are valid. Although we did show that the spatial profiles of both modalities are consistent with one another, we found that the optical HbR more closely matched the BOLD data, whereas the
$\mathrm{HbT}\u2215{\mathrm{HbO}}_{2}$
more closely matched the ASL data. This suggests that prior spatial constraints used with these inversion methods should be used in a hemoglobin species-specific manner. While BOLD could provide a reliable spatial constraint on HbR, it does not correlate as well to the arterial-weighted
${\mathrm{HbO}}_{2}$
or HbT changes. Likewise, ASL could provide a spatial prior to HbT, but might provide a less reliable constraint for HbR.

## 5.

## Conclusions

The synergistic information content provided by the DOI and fMRI methodologies opens up the possibility of future multimodality studies. It is clear that, although both methods measure aspects of the hemodynamic response, they do so using very different means. As a result, measurement sensitivity, data processing, and the very different physics of acquisition all can play important roles in interpretation of these results. Vital to the correct interpretation of such multimodality studies, proper control experiments must be first performed. By using the approach of projecting the fMRI signal through the calculated DOI sensitivity profile, we believe we can more quantitatively examine the spatial and temporal correlation between these measurements. In the first part of this study, the common measurement of HbR between the BOLD-fMRI and DOI recordings provided one such control. Our results here show a statistical temporal $(p=1\times {10}^{-13})$ and spatial $(p=6\times {10}^{-6})$ correlation between these measurements. In addition, in study II, we showed that optically measured ${\mathrm{HbO}}_{2}$ and HbT were correlated with the ASL response both temporally ( $p=9\times {10}^{-8}\mid p=3\times {10}^{-12}$ , respectively) and spatially $(p=3\times {10}^{-6}\mid p=1\times {10}^{-6})$ .

These results confirm the hypothesis that the optical and fMRI measurements arise from common hemodynamic physiological parameters. Furthermore, they suggest the measurements of ${\mathrm{HbO}}_{2}$ , HbT, and ASL probe different vascular territories from HbR and the BOLD signal. These results have implications toward the future of the multimodality integration of DOI and fMRI.

## Acknowledgments

We would like to thank Christiana Andre, Danny Joseph, Elizabeth Alf, and Monica Allen for helping in organizing and performing the experimental sessions. T. J. H. is funded by the Howard Hughes Medical Institute predoctorial fellowship program. This work was supported by the National Institutes of Health (Grants No. R01-EB002482, RO1-EB001954, and P41-RR14075).

## References

**,” Neurosci. Lett., 14 (1–2), 101 –104 (1993). 0304-3940 Google Scholar**

*Near infrared spectroscopy (DOI): A new tool to study hemodynamic changes during activation of brain function in human adults***,” J. Appl. Physiol., 75 (5), 1842 –1846 (1993). 8750-7587 Google Scholar**

*Dynamic multichannel near-infrared optical imaging of human brain activity***,” Adv. Exp. Med. Biol., 454 103 –113 (1998). 0065-2598 Google Scholar**

*BOLD MRI vs. NIR spectrophotometry. Will the best technique come forward?***,” Neuroimage, 25 701 –707 (2005). 1053-8119 Google Scholar**

*Simultaneous recording of task-induced changes in blood oxygenation, volume, and flow using diffuse optical imaging and arterial spin-labeling MRI***,” Phys. Med. Biol., 49 (12), N155 –N163 (2004). https://doi.org/10.1088/0031-9155/49/12/N01 0031-9155 Google Scholar**

*Diffuse optical tomography with physiological and spatial**a priori*constraints**,” Appl. Opt., 44 (10), 1957 –1968 (2005). https://doi.org/10.1364/AO.44.001957 0003-6935 Google Scholar**

*Simulation study of magnetic resonance imaging-guided cortically constrained diffuse optical tomography of human brain function***,” Phys. Med. Biol., 50 (12), 2837 –2858 (2005). https://doi.org/10.1088/0031-9155/50/12/008 0031-9155 Google Scholar**

*Diffuse optical tomography with**a priori*anatomical information**,” Phys. Med. Biol., 48 (4), 417 –427 (2003). https://doi.org/10.1088/0031-9155/48/4/301 0031-9155 Google Scholar**

*Correlation between near-infrared spectroscopy and magnetic resonance imaging of rat brain oxygenation modulation***,” Phys. Med. Biol., 48 1391 –1403 (2003). https://doi.org/10.1088/0031-9155/48/10/311 0031-9155 Google Scholar**

*Temporal comparison of functional brain imaging with diffuse optical tomography and fMRI during rat forepaw stimulation***,” Neuroimage, 19 (4), 1521 –1531 (2003). 1053-8119 Google Scholar**

*The roles of changes in deoxyhemoglobin concentration and regional cerebral blood volume in the fMRI BOLD signal***,” Magn. Reson. Med., 54 354 –365 (2005). https://doi.org/10.1002/mrm.20511 0740-3194 Google Scholar**

*Concurrent fMRI and optical measures for the investigation of the hemodynamic response function***,” Neuroimage, 29 (2), 368 –382 (2006). 1053-8119 Google Scholar**

*A temporal comparison of BOLD, ASL, and DOI hemodynamic responses to motor stimuli in adult humans***,” Med. Phys., 28 (4), 521 –527 (2001). https://doi.org/10.1118/1.1354627 0094-2405 Google Scholar**

*Investigation of human brain hemodynamics by simultaneous near-infrared spectroscopy and functional magnetic resonance imaging***,” IEEE Comput. Sci. Eng., 2 63 –77 (1995). https://doi.org/10.1109/99.476370 1070-9924 Google Scholar**

*MRI-guided optical tomography: Prospectives and computation for a new imaging method***,” Opt. Lett., 23 1716 –1718 (1998). 0146-9592 Google Scholar**

*High-resolution near-infrared tomographic imaging of the rat cranium by use of**a priori*magnetic resonance imaging structural information**,” Opt. Express, 10 159 –170 (2002). 1094-4087 Google Scholar**

*Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head***,” Comput. Methods Programs Biomed., 47 131 –146 (1995). https://doi.org/10.1016/0169-2607(95)01640-F 0169-2607 Google Scholar**

*MCML—Monte Carlo modeling of light transport in multi-layered tissues***,” Optical-Thermal Response of Laser-Irradiated Tissue, Plenum, New York (1995). Google Scholar**

*Monte Carlo modeling of light transport in tissues***,” Inverse Probl., 15 14 –93 (1999). https://doi.org/10.1088/0266-5611/15/2/022 0266-5611 Google Scholar**

*Optical tomography in medical imaging***,” Opt. Lett., 29 (3), 256 –258 (2004). https://doi.org/10.1364/OL.29.000256 0146-9592 Google Scholar**

*Reconstructing chromosphere concentration images directly by continuous-wave diffuse optical tomography***,” Opt. Lett., 30 (15), 1968 –1970 (2005). https://doi.org/10.1364/OL.30.001968 0146-9592 Google Scholar**

*Spectral priors improve near-infrared diffuse tomography more than spatial priors***,” Med. Biol. Eng. Comput., 26 289 –294 (1988). https://doi.org/10.1007/BF02447083 0140-0118 Google Scholar**

*System for long-term measurement of cerebral blood flow and tissue oxygenation on newborn infants by infrared transillumination***,” Phys. Med. Biol., 33 1433 –1442 (1988). https://doi.org/10.1088/0031-9155/33/12/008 0031-9155 Google Scholar**

*Estimation of optical pathlength through tissue from direct time of flight measurement***,” Neuroimage, 18 (4), 865 –879 (2003). 1053-8119 Google Scholar**

*Factors affecting the accuracy of near-infrared spectroscopy concentration calculations for the focal changes in oxygenation parameters***,” Med. Phys., 28 1108 –1114 (2001). https://doi.org/10.1118/1.1373401 0094-2405 Google Scholar**

*Wavelength dependence of the precision of noninvasive optical measurement of oxy-, deoxy-, and total-hemoglobin concentration***,” Neuroimage, 21 1554 –1562 (2004). 1053-8119 Google Scholar**

*Practicality of wavelength selection to improve signal-to-noise ratio in near-infrared spectroscopy***,” Phys. Med. Biol., 46 2227 –2237 (2001). https://doi.org/10.1088/0031-9155/46/8/313 0031-9155 Google Scholar**

*In vivo*optical characterization of human tissues from 610 to $1010\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ by time-resolved reflectance spectroscopy**,” Hum. Brain Mapp, 8 109 –114 (1999). 1065-9471 Google Scholar**

*Optimal experimental design for event-related fMRI***,” Psychophysiology, 42 (16), 3063 –3072 (2003). 0048-5772 Google Scholar**

*Hemodynamic evoked response of the sensorimotor cortex measured non-invasively with near-infrared optical imaging***,” Appl. Opt., 42 (16), 2915 –2922 (2003). 0003-6935 Google Scholar**

*Near-infrared light propagation in an adult head model. II. Effect of superficial tissue thickness on the sensitivity of the near-infrared spectroscopy signal***,” Neuroimage, 23 S275 –S288 (2004). 1053-8119 Google Scholar**

*Diffuse optical imaging of brain activation: Approaches to optimizing image sensitivity, resolution and accuracy***,” Neuroimage, 21 1690 –1700 (2004). 1053-8119 Google Scholar**

*A comparison of methods for characterizing the event-related BOLD time series in rapid fMRI***,” Magn. Reson. Med., 42 1014 –1018 (1999). https://doi.org/10.1002/(SICI)1522-2594(199912)42:6<1014::AID-MRM4>3.0.CO;2-F 0740-3194 Google Scholar**

*Real-time 3D image registration for functional MRI***,” Neuroimage, 2 (3), 173 –181 (1995). 1053-8119 Google Scholar**

*Analysis of fmri time-series revisited—again***,” NMR Biomed., 4–5 (10), 237 –249 (1997). 0952-3480 Google Scholar**

*Implementation of quantitative perfusion imaging techniques for functional brain mapping using pulsed arterial spin labeling***,” Magn. Reson. Med., 6 (41), 1246 –1254 (1999). https://doi.org/10.1002/(SICI)1522-2594(199906)41:6<1246::AID-MRM22>3.0.CO;2-N 0740-3194 Google Scholar**

*QUIPSS II with thin-slice ti1 periodic saturation: A method for improving accuracy of quantitative perfusion imaging using pulsed arterial spin labeling***,” Neuroimage, 9 179 –194 (1999). 1053-8119 Google Scholar**

*Cortical surface-based analysis I: Segmentation and Surface Reconstruction***,” Med. Image Anal, 1 (2), 109 –127 (1996). 1361-8415 Google Scholar**

*Segmentation of brain tissue from magnetic resonance images***,” Neuroimage, 9 (5), 461 –476 (1999). 1053-8119 Google Scholar**

*Brain segmentation and the generation of cortical surfaces***,” Magn. Reson. Med., 39 (6), 855 –864 (1998). 0740-3194 Google Scholar**

*Dynamics of blood flow and oxygenation changes during brain activation: the balloon model***,” Neuroimage, 21 144 –153 (2004). 1053-8119 Google Scholar**

*Discrepancies between BOLD and flow dynamics in primary and supplementary motor areas: Applications of the balloon model to the interpretation of BOLD transients***,” Probability and Statistics for Engineering and the Sciences, 474 –522 4th ed.Wadsworth Inc., Belmont, California (1995). Google Scholar**

*Simple linear regression and correlation***,” Clin. Neurophysiol., 113 (5), 621 –634 (2002). 1388-2457 Google Scholar**

*Technical aspects and utility of fMRI using BOLD and ASL”***,” J. Cereb. Blood Flow Metab., 19 (8), 871 –879 (1999). https://doi.org/10.1097/00004647-199908000-00006 0271-678X Google Scholar**

*Simultaneous blood oxygenation level-dependent and cerebral blood flow functional magnetic resonance imaging during forepaw stimulation in the rat***,” Magn. Reson. Med., 44 137 –143 (2000). https://doi.org/10.1002/1522-2594(200007)44:1<137::AID-MRM20>3.3.CO;2-I 0740-3194 Google Scholar**

*Comparison of simultaneously measured perfusion and BOLD signal increases during brain activation with T-1-based tissue identification***,” Neuroimage, 12 287 –297 (2000). 1053-8119 Google Scholar**

*A CBF based event-related brain activation paradigm: Characterization of impulse-response function and comparison to BOLD***,” Magn. Reson. Med., 52 1407 –1417 (2004). https://doi.org/10.1002/mrm.20302 0740-3194 Google Scholar**

*Simultaneous MRI acquisition of blood volume, blood flow, and blood oxygenation information during brain activation***,” Magn. Reson. Med., 34 4 –10 (1995). 0740-3194 Google Scholar**

*The intravascular contributions of fMRI signal change: Monte Carlo modeling and diffusion-weighted studies**in vivo***,” Magn. Reson. Med., 45 233 –246 (2001). https://doi.org/10.1002/1522-2594(200102)45:2<233::AID-MRM1032>3.0.CO;2-W 0740-3194 Google Scholar**

*Quantitative differentiation between BOLD models in fMRI***,” Med. Phys., 27 (4), 801 –815 (2000). https://doi.org/10.1118/1.598943 0094-2405 Google Scholar**

*Near-infrared study of fluctuations in cerebral hemodynamics during rest and motor stimulation: Temporal analysis and spatial mapping***,” Magn. Reson. Imaging, 22 (1), 9 –13 (2004). https://doi.org/10.1016/j.mri.2003.08.003 0730-725X Google Scholar**

*Cardiac-induced physiologic noise in tissue is a direct observation of cardiac-induced fluctuations***,” Neuroimage, 21 (4), 1652 –1664 (2004). 1053-8119 Google Scholar**

*Resting fluctuations in arterial carbon dioxide induce significant low frequency variations in BOLD signal***,” Magn. Reson. Med., 35 (1), 107 –113 (1996). 0740-3194 Google Scholar**

*Reduction of physiological fluctuations in fMRI using digital filters***,” Hum. Brain Mapp, 6 (3), 160 –188 (1998). 1065-9471 Google Scholar**

*Analysis of fMRI data by blind separation into independent spatial components***,” Neuroimage, 30 (1), 88 –101 (2006). 1053-8119 Google Scholar**

*Dynamic physiological modeling for functional difuse optical tomography***,” J. Biomed. Opt., 10 011014 (2005). https://doi.org/10.1117/1.1852552 1083-3668 Google Scholar**

*Eigenvector-based spatial filtering for reduction of physiological interference in diffuse optical imaging***,” NMR Biomed., 7 (1–2), 45 –53 (1994). 0952-3480 Google Scholar**

*Brain or vein–oxygenation or flow? On signal physiology in functional MRI of human brain activation***,” Phys. Med. Biol., 47 (18), N249 –N257 (2002). https://doi.org/10.1088/0031-9155/47/18/402 0031-9155 Google Scholar**

*A haemodynamic model for the physiological interpretation of**in vivo*measurements of the concentration and oxygen saturation of haemoglobin**,” Magn. Reson. Med., 38 (2), 296 –302 (1997). 0740-3194 Google Scholar**

*Experimental determination of the BOLD field strength dependence in vessels and tissue***,” J. Biomed. Opt., 10 044002 (2005). https://doi.org/10.1117/1.2007987 1083-3668 Google Scholar**

*Noninvasive detection of functional brain activity with near-infrared diffusing-wave spectroscopy***,” Opt. Lett., 29 (15), 1766 –1768 (2004). https://doi.org/10.1364/OL.29.001766 0146-9592 Google Scholar**

*Diffuse optical measurement of blood flow, blood oxygenation, and metabolism in a human brain during sensorimotor cortex activation***,” IEEE Trans. Med. Imaging, 16 (6), 750 –761 (1997). https://doi.org/10.1109/42.650872 0278-0062 Google Scholar**

*Multimodality Bayesian algorithm for image reconstruction in positron emission tomography: A tissue composition model***,” Appl. Opt., 42 (16), 2881 –2887 (2003). 0003-6935 Google Scholar**

*Monte Carlo prediction of near-infrared light propagation in realistic adult and neonatal head models*