1.1.

Theory

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\mathrm{Hz}$ ). After $2\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.²⁸ The length of the ISI ranged from 4 to $20\mathrm{s}$ with an average ISI period of $12\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\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\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.²⁹ The DOI imager has 18 lasers—9 lasers at $690\mathrm{nm}$ $\left(18\mathrm{mW}\right)$ and 9 at $830\mathrm{nm}$ $\left(7\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\mathrm{nm}$ and $830\mathrm{nm}$ ( $18\mathrm{mW}$ and $7\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\mathrm{kHz}$ . A bandpass filter follows each APD module with a cut-on frequency of approximately $500\mathrm{Hz}$ to reduce $1∕f$ noise and the 60-Hz room light signal, and a cutoff frequency of approximately $16\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 ²⁹ 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\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.³⁰ 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.

Fig. 1

Here we show the probe design and placement for DOI measurements. This probe consisted of eight detector positions (most medial and lateral rows) and four source positions (middle row). The source-detector separation was $2.9\mathrm{cm}$ . This probe was positioned approximately over the subject’s contralateral primary motor cortex as the subject lay in the MRI scanner. The probe also contains vitamin E rings placed on the caps of the optodes, which allowed registration of the optode array in the MR structural image.

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).³¹ Signals were further low-pass filtered at $0.8\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).

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 ³⁴)], 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\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}∕\mathrm{TE}∕\alpha =500\mathrm{ms}∕30\mathrm{ms}∕90\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}∕\mathrm{TI}∕\mathrm{TE}∕\alpha =2530\mathrm{ms}∕1100\mathrm{ms}∕3.25\mathrm{ms}∕7\mathrm{deg}$ ]. To calculate the BOLD-based hemodynamic response functions, the functional images were first motion-corrected³⁵ 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.³⁶ The mean of the response from 2 to $7\mathrm{s}$ was used to calculate the effect size.

For study II, ASL-fMRI was carried out at $3\mathrm{T}$ (same scanner as study I) using proximal inversion with control for off-resonance effects (PICORE) labeling geometry³⁷ with Q2TIPS saturation³⁸ to impose a controlled label duration. A postlabel delay of $1400\mathrm{ms}$ and a label duration of $700\mathrm{ms}$ were used, with repetition and echo times of $2\mathrm{s}$ and $20\mathrm{ms}$ , respectively $[\alpha =90\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\mathrm{s}$ , the stimulus had been jittered evenly on a 500-ms time step, which allowed for the response to be calculated at $2\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\mathrm{Hz}$ using a cubic spline model.³⁹ 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\mathrm{s}$ was used to calculate the effect size. Figure 2 depicts the processing stream.

Fig. 2

Functional and structural data was processed through several stages before comparison between DOI and fMRI. During the experimental session, anatomical and functional MRI scans were preformed along with the simulations optical imaging. The anatomical images were registered and segmented into a five-layered model (skin, skull, CSF, and gray and white matter). This was used to perform Monte Carlo simulations to determine the light propagation through the head and estimate path-length factors. These forward models were then used to project the functional MRI data into the same source-detector measurement space as the optical data.

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).⁴⁰ This process first consisted of intensity normalization, skull stripping, and white-gray matter differentiation as described by Dale ⁴¹ 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 ¹⁶ The optical properties of each of the segmented tissues that were used in these simulations are listed in Table 1 .²⁴ 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.

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\mathrm{cm}$ axially [5 slices $\times$ (6mm $+$ 1-mm gap $=7\left)\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\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.

Fig. 3

Monte Carlo simulations were used to generate the optical sensitivity profiles for each source and detector pair. A representative profile is shown here overlaying the anatomical volume (subject B). This image is shown in logarithmic scale and contours are presented for each order of magnitude. The inset at the top right shows the location of this coronal slice relative to the optical probe. To show the location of the activation, the $t$ map for the BOLD signal (masked at half maximum) is shown overlaying the anatomical volume. The forward modeling of the BOLD signal involves the multiplication of the BOLD activation and the sensitivity matrix $\left(\mathsf{A}\right)$ . This is equivalent to the sensitivity-weighted sum over the volume of the BOLD contrast changes.

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\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.⁴⁶

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

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,

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.

Fig. 4

Following the processing stream outlined in Fig. 2, BOLD data was projected through the simulated optical sensitivity (forward) matrices to produce predicted spatiotemporal profiles of the response across the optical probe. Shown above, correspondence was seen between the response predicted from the BOLD-fMRI data and the measured optical HbR response. The DOI and BOLD time courses above have been normalized to the single maximum response to allow comparisons of the amplitudes across the probe. The HbR time courses have been inverted for this 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\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 .

Fig. 5

In order to better investigate the spatial correlation between the time course shown in Fig. 4, the mean response amplitude of the period of 3 to 9 seconds was calculated for each source-detector pair. Shown as bar graphs above, the amplitude variation across the probe was coincident between the two modalities.

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.

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}∕{\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}∕{\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.

Fig. 6

Similar to the data presented in Fig. 4, ASL data was projected through the simulated optical sensitivity (forward) matrices to produce predicted spatiotemporal profiles of the response across the optical probe. Shown above, correspondence was seen between the response predicted from the ASL-fMRI data and the measured optical HbT response. The DOI and ASL time courses above have been normalized to the single maximum response to allow comparisons of the amplitudes across the probe.

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 .

Fig. 7

The maximum response amplitude was calculated from the optically measured HbT and ASL data for each source-detector pair. Shown as bar graphs above, the amplitude variation across the probe was coincident between the two modalities.

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.

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}∕{\mathrm{HbO}}_2$ spatial correlations are better than that observed between ASL/HbR ( $R=0.26$ , $p=0.10$ ). These correlations between $\mathrm{ASL}∕{\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\mathrm{ms})$ . However, similar to the results presented in Table 2, the correlation between BOLD and HbR is higher than that between $\mathrm{BOLD}∕{\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,⁴⁷ 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, ⁴⁸ Luh, ⁴⁹ and Yang ^{50, 51} Shown in Fig. 8 are the ASL and BOLD activation maps for one such subject (F).

Fig. 8

Although both fMRI methods measure hemodynamic responses, the BOLD and ASL activations were observed to locate to slightly displaced positions in the brain. Shown above is a single axial slice from one of the subjects (F) with both the BOLD and ASL effects maps overlaid. Here the $t$ statistics maps are shown with a half-maximum mask. In both images, crosses mark the positions of the centroids of activation in this slice. The average displacement between the BOLD and ASL maps was $7.0\mathrm{mm}\pm 3.1\mathrm{mm}$ as shown in Table 4. The approximate location of the optical probe is given as an intensity projection and indicated in dotted circles overlaying these images.

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.

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}-⟨\mathit{r}⟩{)}^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\mathrm{mm}$ ).

3.4.

Temporal Correlation between DOI and fMRI

In our previous work,¹² 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,¹² 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}∕{\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.

Fig. 9

Using the fMRI following projection by the simulated optical forward matrices, region-of-interest average time courses were calculated from the average of all source-detector pairs. The relative amplitudes have not been rescaled following averaging such that response amplitudes within each comparison are comparable.

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\mathrm{mm}$ $(\pm 5\mathrm{mm})$ for the BOLD and $17\mathrm{mm}$ $(\pm 5\mathrm{mm})$ for the ASL. The mean depth of the optical contrast was shallower at only $14\mathrm{mm}$ $(\pm 4\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\mathrm{mm}$ $(\pm 2\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\mathrm{mm}\right)$ compared to the anatomical volume $\left(1\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.

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.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, ⁴⁴ Obata ⁴⁵). 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.⁵² 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.⁵² 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}∕{\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, ⁵⁴ Bhattacharyya and Lowe,⁵⁵ Wise ⁵⁶). 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\mathrm{cm}$ into the brain.⁶⁷ 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, ⁵ Boas and Dale⁶). With functional MRI to provide prior information for the location of activation, image reconstruction schemes can make use of this to constrain solutions.⁷ 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}∕{\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.