Investigation of human brain hemodynamics by simultaneous near-infrared spectroscopy and functional magnetic resonance imaging.

The aim of this study was to compare functional cerebral hemodynamic signals obtained simultaneously by near infrared spectroscopy ~ NIRS ! and by functional magnetic resonance imaging ~ fMRI ! . The contribution of superﬁcial layers ~ skin and skull ! to the NIRS signal was also assessed. Both methods were used to generate functional maps of the motor cortex area during a periodic sequence of stimulation by ﬁnger motion and rest. In all subjects we found a good collocation of the brain activity centers revealed by both methods. We also found a high temporal correlation between the BOLD signal ~ fMRI ! and the deoxy-hemoglobin concentration ~ NIRS ! in the subjects who exhibited low ﬂuctuations in superﬁcial head tissues. © 2001 American Association of Physicists in Medicine. @ DOI: 10.1118/1.1354627 #


I. INTRODUCTION
In recent years near-infrared spectroscopy ͑NIRS͒ has been proposed as a method to study brain hemodynamics, which is simple and inexpensive compared to such ''heavy-duty'' methods as functional magnetic resonance imaging ͑fMRI͒ and positron emission tomography ͑PET͒. Many studies about applying NIRS to monitor functional cerebral hemodynamics were reviewed in Ref. 1. However, the precise location of tissues contributing to the optical signals is still uncertain. This issue can be clarified using coregistration of brain hemodynamics by NIRS and fMRI, because the blood oxygenation level dependent ͑BOLD͒ signal employed in fMRI is proportional to the changes in deoxy-hemoglobin concentration. In several studies functional brain activity was monitored separately by NIRS and fMRI in the same subjects. 2,3 However, due to strong biological fluctuations in both NIRS and fMRI signals, such a separate monitoring makes it difficult to compare the signals directly. We only know one paper describing simultaneous NIRS and MRI recordings of cerebral blood oxygenation during functional stimulation. 4 Using a continuous-wave NIRS instrument with one source-detector pair, the authors demonstrated spatial congruence between the area on the surface of the head where NIRS deoxy-hemoglobin signal was maximal and the focus of the MRI signal in the motor cortex related to the finger-movement task. However, since Kleinschmidt et al. employed only one slice in their fMRI setup, the contribution of extracranial artifacts to the hemodynamic changes measured by NIRS remained unclear. Furthermore, their instrument with only one source-detector pair did not allow deter-mining differences between hemodynamic changes in the activated area and outside of it. The aim of our paper is to address these questions.
In our study we use motor stimulations to evoke functional cerebral hemodynamic changes, a multichannel frequency-domain NIRS instrument, and a set of five fMRI slices at different depths of the head from the skin to the gray matter.

II. INSTRUMENTATION AND EXPERIMENTAL PROTOCOL
For NIRS measurements we use a two-wavelength ͑758 and 830 nm͒ frequency-domain ͑110 MHz modulation fre-quency͒ Oximeter ͑ISS, Champaign, IL͒, which has 16 laser diodes ͑eight per each wavelength͒ and two photomultiplier tube detectors. At a wavelength of 758 nm light absorption by the deoxy-hemoglobin ͑HHb͒ substantially exceeds absorption by the oxy-hemoglobin (O 2 Hb), while at 830 nm the O 2 Hb absorption prevails over the HHb absorption. The laser diodes operate in a sequential multiplexing mode with 10 ms ''on'' time per each diode. Light emitted by these laser diodes is guided to the tissue through 10-meter long multimode silica optical fibers. Two 10-m long glass fiber bundles collect the scattered light and bring it to the detectors. The paired ͑758 and 830 nm wavelength͒ source fibers are attached to the probe at eight positions. Together with two detectors, they provide 10 bi-wavelength source-detector channels with a source-detector distance of 3 cm. The probe covers an area 9ϫ6 cm 2 . The geometric layout of the optical probe and its position on the head are shown in Fig. 1. The numbers 1 through 8 correspond to the light sources and the letters A and B to the detectors. Using these notations, in this paper we denote the source-detector light channels as A1, A2 and so on through B8. The probe is centered at the measured C3 position according to the International 10-20 System. 5 Three multimodality radiological markers ͑IZI Medical Products Corp, Baltimore, MD͒ are embedded into the optical probe to facilitate correct orientation of the MRI slices with respect to the probe and to enable recovery of the probe orientation for data analysis.
Magnetic resonance imaging was performed using a 1.5 Tesla whole body MR scanner ͑Signa, General Electric Medical Systems, Milwaukee, WI͒ equipped with echospeed gradients and a standard circularly polarized birdcage headcoil. Sagittal T1-weighted localizer scans were used to determine the correct plane for the functional scans. Gradientecho echo-planar images were acquired using a data matrix of 64ϫ64 complex points, TRϭ640 ms, TEϭ40 ms, FOV ϭ240 mm, slice thicknessϭ7 mm, no interslice gap, receiver bandwidth 62.5 kHz, and tip angle 90 degrees. Figure 2 shows the arrangement of a typical set of five slices denoted as ␣, ␤, ␥, ␦, and ⑀. The slices are parallel to the plane of three radiological markers on the optical probe. The middle slice ␥ was set between the skull and the brain surface at the C3 position. This slice was mostly filled with the cerebrospinal fluid ͑CSF͒, dura mater, arachnoidal tissue, and pia mater, but also could touch the cortical tissue depending on the subject's anatomy. Two deeper slices, ␣ and ␤, mostly contained the brain tissue, while two outer slices, ␦ and ⑀, included the skull, the skin and the markers ͑see The studies were performed in six healthy right-handed male volunteers, 18 to 37 years old. Informed consent was obtained from all subjects. Each exercise run consisted of a 30-s pre-exercise epoch, ten 20-s stimulation epochs separated by ten 20-s control epochs, and a 50-s after-exercise epoch. Thus, the stimulation period was 40 s, which corresponds to 0.025 Hz in the frequency domain. During stimulation epochs subjects performed light palm squeezing with the right hand with the frequency about 1.25 Hz following the rhythm presented by the sound from the vacuum pump of the MRI scanner. This rhythm was selected to maintain the same amount of finger movements ͑about 50͒ per stimulation period. The commands to begin and to stop the exercise were presented via the speaker installed in the scanner room and by a light, which was on during the exercise epoch and off during the control epoch. 750 images acquired during each run resulted in an 8-min scan time. The synchronization of the exercise sequence with the MRI and NIRS recording was provided by a computer program generating the commands for the subject and the scanner operator based on the preset command timing.

III. SIGNAL ANALYSIS
To convert optical intensity data into hemodynamic concentration changes, we used a model of light transport in strongly scattering medium based on the modified Lambert-Beer law. 1 We assumed that the changes in the light intensity occurred mainly due to the absorption changes rather than due to the scattering ones and that the absorption variations in the tissue were caused only by changes in HHb and O 2 Hb concentrations. The modified Lambert-Beer law gives the following equation to compute the change in the absorption coefficient ⌬ a at a given wavelength from the initial intensity I 0 and the current intensity I : where d is the distance between the source and the detector fibers, and is the differential pathlength factor ͑DPF͒. We use DPF values for the adult head of 6.2 at 758 nm and 5.9 at 830 nm ͑Ref. 6͒. The DPF dependence on the wavelength reflects a weak wavelength dependence of the tissue scattering. To avoid noise from the ambient light, we use the amplitude ͑AC͒ of the intensity-modulated light as the source data for I . Changes of ͓O 2 Hb͔ and ͓HHb͔ corresponding to the changes in the absorption coefficients ⌬ a 1,2 of the tissue at the wavelengths 1 and 2 can be obtained using the
Equations ͑1͒, ͑2͒, and ͑3͒ are derived assuming homogeneity of the scattering medium, and therefore they are not expected to provide accurate quantitative intracranial hemodynamics assessment. We do not aim at quantifying the absolute hemodynamic changes in the motor cortex, but rather we aim at the analysis of temporal characteristics, changes in ͓HHb͔ relative to changes in ͓O 2 Hb͔, and relative changes in different regions of the cortex. For this reason we quantify ͓HHb͔ and ͓O 2 Hb͔ changes by arbitrary units, which correspond to M under the hypothetical assumption of homogeneity of tissue optical properties.
The time series obtained from the optical data using Eqs. ͑2͒ and ͑3͒ were bandpass filtered to exclude fluctuations of the time scale faster than 2 s and to eliminate the changes of the time periods significantly ͑approximately three times͒ longer than the stimulation/relaxation period.
Assuming synchronization between stimulation and hemoglobin response, we separate the functional response from background fluctuations by means of a folding average ͑time-locked average 4 ͒, i.e., the data corresponding to the same point on the stimulation/relaxation cycle is averaged over a number of cycles. The beginning of the time-locked period was chosen in the middle of the relaxation epoch.
The analysis of the fMRI data was performed using MEDx software ͑Sensor Systems, Inc͒. Prior to the statistical analysis, the data was subjected to the motion detection, spatial and temporal filtering, and global intensity normalization. Data with significant motion artifacts were discarded and measurements repeated. The resulting statistical images show the maps of coefficient r measuring correlation between voxel BOLD intensity and a reference function. As the reference function we used either the boxcar paradigm function with the ''ON'' and ''OFF'' durations equal to the stimulation and relaxation phases of the exercise, or the ͓HHb͔ signal obtained by NIRS taken with the opposite sign.

A. BOLD signal
In all subjects the analysis of the BOLD signal revealed an area under the optical probe where the signal was highly correlated with the paradigm boxcar function. It was an area in the primary motor cortex with the center close to the central sulcus. The size of this area varied depending on the slice. Figures 3͑d͒, 3͑e͒, and 3͑f͒ show the typical correlation coefficient maps of the functional signals for the slices ␣, ␤, and ␥ shown in Figs. 3͑a͒, 3͑b͒, and 3͑c͒ correspondingly. White color of the highest intensity corresponds to the voxels with r greater than 0.5 ͑correlation z score larger than 9.0͒. Correlation z score is the voxel correlation coefficient divided by the standard deviation of the distribution of correlation coefficients over the whole image. It represents a significance of correlation in the given voxel relative to the voxels whose correlation is within standard deviation from the average one. Usually, the largest size of the activated area was in the middle slice ␥ ͓see Fig. 3͑f͔͒. In the next slice deeper in the brain ␤, which was mostly cortical tissue, the activation area was almost as large as in the middle slice ͓Fig. 3͑e͔͒. The deepest slice ␣ ͓see Fig. 3͑d͔͒ typically showed a significantly smaller activation area. Using head landmarks for the C3 position, the center of the optical probe was placed very close to the central sulcus so that light channels A1, A5, B1, and B5 were above the activated area ͑see Fig. 3͒. In two subjects the probe center was displaced from the central sulcus by 1.5-2 cm toward the frontal lobe, and the major activated area appeared under the channels B5, B6, B7, and B8. In two upper slices ␦ and ⑀ the coefficient of correlation between BOLD signal and the paradigm function was typically less than 0.1 and the corresponding z score was less than 2.0. Figure 4͑a͒ shows a BOLD signal from the slice ␤ ͓see Fig. 3͑e͔͒ with rϭ0.7, z score of 12.5. Figures 4͑b͒ and 4͑c͒ show, respectively, the ͓HHb͔ and the ͓O 2 Hb͔ signals from the light channel A5 ͑Fig. 3͒ acquired simultaneously with the BOLD signal shown in Fig.  4͑a͒. One can see that the directions of changes in BOLD and ͓HHb͔ signals are almost perfectly opposite. In three subjects the task-related changes were superimposed with the significant hemodynamic fluctuations in the upper tissue layers corresponding to slices ␦ and ⑀. However, in all subjects the folding average analysis revealed the same type of ͓O 2 Hb͔ and ͓HHb͔ behavior for light channels situated above the major activated area. Figure 5 shows a map of folding average ͓HHb͔ and ͓O 2 Hb͔ traces acquired simultaneously with the fMRI maps presented in Fig. 3. For the light channels situated above the activated area ͑channels A1 and A5͒, the characteristic feature was a significant decrease of ͓HHb͔ during the stimulation, which was concurrent with a significant increase of the oxy-hemoglobin concentration. Typically rapid changes in ͓HHb͔ and ͓O 2 Hb͔ began 2-3 s after the stimulation onset and continued during next 7-15 s. The maximal task-related change in ͓O 2 Hb͔ occurred from 0.5 to 4 s before the maximal change in ͓HHb͔. Then ͓HHb͔ fluctuated near its low level for the rest of the stimulation epoch and the beginning of the resting epoch. After achieving its maximum, ͓O 2 Hb͔ either fluctuated at its high level or exhibited a slight decrease ͑see channels A1, A5, and B1 in Fig. 5͒. A rapid recovery toward the baseline level begins 4-6 s after the onset of the rest epoch in both ͓HHb͔ and ͓O 2 Hb͔. Such a behavior of ͓HHb͔ and ͓O 2 Hb͔ in the activated area was qualitatively the same in all six subjects. Table I Fig. 3͑e͔͒ with rϭ0.7, z score of 12.5; ͑b͒ and ͑c͒ ͓HHb͔ and ͓O 2 Hb͔ signals, respectively, from the light channel A5 acquired simultaneously with the BOLD signal shown in ͑a͒.  error for each value of ͓HHb͔ or ͓O 2 Hb͔ folding-average trace is calculated by averaging over the 10 exercise periods. Then the average over the period is taken. One can see that for all subjects the task-related change exceeds the average standard error at least by a factor of 10. The values of the ratio of ͓HHb͔ and ͓O 2 Hb͔ changes presented in the last column show that ͓HHb͔ changes were always smaller than the corresponding ͓O 2 Hb͔ changes. No significant decrease in the folding average ͓HHb͔ traces concurrent with the significant ͓O 2 Hb͔ increase were observed during stimulations in the light channels outside the activated area. Usually ͓HHb͔ in such channels ͑see channels A2, A3, A4, B8, B7, and B6 in Fig. 3 and the corresponding panels in Fig. 5͒ fluctuated without correlation with the paradigm function. In some cases ͓O 2 Hb͔ increased ͑see channel B8 in Fig. 5͒, but without significant ͓HHb͔ change. In two subjects both ͓HHb͔ and ͓O 2 Hb͔ decreased in some channels during the stimulation indicating a decrease in blood volume in these nonactivated regions.

C. Temporal correlation between BOLD and NIRS signals
Significant changes in the folding average traces corresponding to the activated area demonstrate a correlation between the hemodynamic changes and the paradigm function.
To determine the location of the tissue contributing to such a correlation we performed a correlation analysis of the BOLD signals using the inverse NIRS ͓HHb͔ signals as the reference function. ͓͑HHb͔ signal was taken with the opposite sign because an increase in BOLD signal should correspond to a decrease in ͓HHb͔.͒ In Figs. 3͑a͒, 3͑b͒, and 3͑c͒ the arrows point at the clusters of voxels where z score for the correlation coefficient between BOLD and ͓HHb͔ signals in channel A5 exceeded 10.0. One can see a good collocation of these clusters with the channel A5. Such a direct temporal correlation between optical and BOLD signals in cerebral and near-cerebral tissues was detected in three subjects.

D. Influence of background hemodynamic fluctuations
In other subjects the correlation coefficient between the intracranial BOLD signal and optically measured hemodynamic signals did not exceed 0.4, and the correlation images did not contain clusters of z scores larger than 5, although the folding average traces exhibited features similar to the ones of channels A1, A5, and B5 in Fig. 5 ͑see Fig. 6͒. We believe that this low temporal correlation between the intracranial BOLD and optical signals was due to strong background hemodynamic fluctuations in optical signals contributed by the superficial tissues. Such a case is illustrated in Figs. 6 and 7. Figure 6 ͑a͒ shows task-related ͓HHb͔ changes superimposed with the strong fluctuations of the period that is approximately twice as long as the task period. Because the period of fluctuations was close to the exercise period, it was difficult to effectively eliminate the corresponding fluctuations by high-pass filtering without eliminating task-related changes. Although the amplitude of the background long-period changes is so strong that the task-related changes are almost imperceptible in the time trace in Fig. 6͑a͒, the folding average reveals the task-related ͓HHb͔ change ͓Fig. 6͑b͔͒. Figure 7͑a͒ shows the power spectra of BOLD signals from the activated area in slice ␥ ͑thick line͒ and the top slice ⑀ ͑thin line͒. The thick line has a high peak at the frequency corresponding to the paradigm function period ͑0.025 Hz͒, while the thin line peaks at the frequency of the background fluctuations and its harmonics. Figure 7͑b͒ shows the spectrum of the ͓HHb͔ signal from the light channel situated above the activated area. One can see that although the ͓HHb͔ spectrum has significant power at the exercise frequency, there is a much higher peak at the frequency of the highest peak in the superficial BOLD signal ͓a thin line in Fig. 7͑a͔͒.

V. DISCUSSION
The decrease of ͓HHb͔ concurrent with the increase in ͓O 2 Hb͔ agrees with the results of the previous simultaneous fMRI-NIRS study, 4 BOLD signal theory and with the basic knowledge of brain physiology. 9 Therefore, the collocation of the light channels exhibiting such a behavior in the folding average ͓HHb͔ and ͓O 2 Hb͔ traces with the activated area revealed by fMRI in all subjects confirms that NIRS monitors the functional brain hemodynamics. Additional support for this conclusion is the temporal correlation between the intracranial BOLD signals and the respective NIRS signals. We observed such a correlation in 50% of subjects in the deep slices ␣, ␤, and ␥, which included brain tissue, dura mater, pia mater, and arachnoidal tissue. The area of correlation was close to the central sulcus and collocated with the corresponding light channels. This indicates that the taskrelated hemodynamic changes measured by NIRS are not due to artifacts, but have intracranial origin.
Although our results on the task-related hemodynamic changes are in qualitative agreement with the ones of Kleinschmidt et al., 4 there is a quantitative discrepancy in the relative magnitudes of ͓HHb͔ and ͓O 2 Hb͔ changes. Kleinschmidt et al. 4 reported that in their measurements the ͓HHb͔ changes were larger than the ͓O 2 Hb͔ ones, we always observed greater changes in ͓O 2 Hb͔. We believe that the reason for this discrepancy is a different set of wavelengths of light sources in our instrument. Kleinschmidt et al. 4 used a four-wavelength instrument ͑NIRO-500͒ with 775, 825, 850, and 904 nm light sources. All these wavelengths are longer than the local maximum in the ͓HHb͔ extinction coefficient at 758 nm, while in our instrument one of two wavelengths was 758 nm. This allows us to perform more accurate measurement of ͓HHb͔ changes. Another source of errors in measurements using NIRO-500 could be a high contribution of water absorption to the signal at 904 nm.
A comparison of simultaneously acquired BOLD and NIRS signals shows that in some subjects the hemodynamic fluctuations in superficial tissues at frequencies close to the exercise frequency make the detection of functional brain signals difficult and may introduce additional errors in measurements. It is a known problem in the NIRS method that the changes in superficial layers contribute much more to the overall signal change than the deep layers. 10 According to Firbank et al., 10 a 1 M change in blood concentration in the brain will give rise only to approximately 0.0005 OD change in attenuation at 3-cm source-detector distance ͓Fig. 9͑a͒ in Ref. 10, no wavelength specified͔. In our measurements the order of magnitude of task-related attenuation changes at both 830 nm and 758 nm was about 0.005 OD. This implies that if the estimation by Firbank et al. 10 is correct, the actual hemoglobin concentration change should be significantly larger than values predicted by the homogeneous model ͑Table I͒. We present changes in ͓HHb͔ and ͓O 2 Hb͔ as the ones measured in ''arbitrary units'' because we believe that the layered nonhomogeneous structure of the head tissues should significantly affect the actual scale of changes calculated using the assumption of homogeneity. To achieve the quantitative accuracy in NIRS of brain hemodynamic changes one should perform the three-dimensional reconstruction of the optical properties of different tissue layers. This should also help to solve the problem of the contamination of the brain signals by superficial fluctuations. However, the three-dimensional reconstruction can only be achieved using multiple pairs of light sources and detectors separated by a variety of different distances. It cannot be applied to our data, which was acquired using the arrangement of equidistant source-detector pairs.
Another source of errors in quantitative hemodynamic measurements could be the assumption of the constancy of the scattering properties. This problem can be solved using the frequency-domain NIRS, which allows assessment of both absorption and scattering changes.

VI. CONCLUSION
Analyzing near-infrared and BOLD signals acquired simultaneously under the motor stimulation conditions, we found a good collocation between the light channels with significant task-related folding average hemodynamic changes and the functionally activated area in the motor cortex in six human subjects. A direct temporal correlation between NIRS and the intracranial BOLD signals was demonstrated in 50% of subjects. These results show the intracranial origin of the NIRS signals obtained under the periodical stimulation conditions. The lack of temporal correlation between the optical and fMRI signals in three subjects was due to the contamination of the optical signal by hemodynamic fluctuations in the superficial tissue layers. The problem of the signal contamination by superficial fluctuations and the spatial nonhomogeneity of the head tissues urge applying methods of three-dimensional reconstruction of the tissue optical properties.

ACKNOWLEDGMENTS
This work was supported by National Institutes of Health ͑NIH͒ CA57032 and RR10966. a͒