Study of local cerebral hemodynamics by frequency-domain near-infrared spectroscopy and correlation with simultaneously acquired functional magnetic resonance imaging

: The aim of our study was to explore the possibility of detecting hemodynamic changes in the brain using the phase of the intensity modulated optical signal. To obtain optical signals with the highest possible signal-to-noise ratio, we performed a series of simultaneous NIRS-fMRI measurements, with subsequent correlation of the time courses of both measurements. The cognitive paradigm used arithmetic calculations, with optical signals acquired with sensors placed on the forehead. Measurements were done on seven healthy subjects. In five subjects we demonstrated correlation between the hemodynamic signals obtained using NIRS and BOLD fMRI. In four subjects correlation was found for the hemodynamic signal obtained using the phase of the intensity modulated signal. ____________________________________________________________________________________


Introduction
The method of frequency-domain near-infrared spectroscopy (NIRS) of biological tissues has a number of advantages over continuous-wave NIRS [1][2][3]. Using numerical simulations of light propagation in the human head Firbank et al. [4] showed that the spatial profile of the sensitivity of phase to changes in optical properties spreads deeper into the brain tissue compared to the sensitivity profile of the amplitude. Also, frequency-domain NIRS allows better separation of absorption and scattering properties of tissue [1], which is important for imaging heterogeneous tissues [5]. However, it is usually difficult to obtain a sufficiently high signal-to-noise ratio in the phase signal to perform dynamic measurements. Although many results on measuring brain hemodynamics using continuous-wave NIRS have been published [6][7][8][9][10], no measurement of hemodynamic changes using the phase of frequency-domain NIRS signal have been reported to date. The goal of the study presented here was to explore the possibility of detecting hemodynamic changes in the brain using frequency-domain NIRS.
To obtain phase signals at the highest possible signal-to-noise ratio, we performed a series of simultaneous NIRS-fMRI measurements, with subsequent correlation of the time courses of both measurements. The cognitive paradigm used arithmetic calculations, with optical signals acquired with sensors placed on the forehead. The arithmetic calculation paradigm was selected because it is known to cause activation in the regions of the frontal cortex situated under the areas of the head surface that are usually free from hair [6,7]. To measure cerebral hemodynamics independently from NIRS, simultaneous fMRI of the brain was employed. In our previous studies of functional brain hemodynamics by simultaneous NIRS and fMRI we observed distortions of functional signals measured by NIRS, which were presumably caused by the hemodynamic fluctuations in extracranial tissues [10]. To account for the hemodynamic changes in different tissue layers we designed an optical sensor having a set of short and long source-detector distances. The short-distance signals were used to correct for the contributions of superficial and systemic fluctuations in the long-distance optical signals. Following our previous study of functional hemodynamics in motor cortex by ( simultaneous NIRS and fMRI [10], we used correlation analysis to reveal correlation in the time course of NIRS signals with the cerebral hemodynamic changes.
In this study, we found that, not only can functional hemodynamic changes in a large volume of the brain tissue be detected using the phase of intensity modulated light, but also that the phase can be sensitive to relatively weak local spontaneous cerebral fluctuations.

Instrumentation and experimental setup
For NIRS measurements we used a two-wavelength (758 and 830 nm) frequency-domain (110 MHz modulation frequency) Oximeter (ISS, Champaign, IL), which had sixteen laser diodes (eight for 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 is much higher than the HHb absorption. The laser diodes operated in a sequential multiplexing mode with 40 ms "on" time for each diode. The whole acquisition cycle for 16 light sources was 640 ms. Light emitted by the laser diodes was guided to the tissue through 10m long multi-mode silica optical fibers (CUDA, FL). Two 10m long glass fiber bundles collected the scattered light and guided it to the detectors. The paired (758 and 830 nm wavelength) source fibers were attached to the sensor at 8 positions. For this study we designed the optical sensor shown in Fig. 1. The sensor had two detectors, marked as A and B in the figure, and sixteen paired sources at eight locations marked by numbers 1 through 8 (758 and 830 nm at each location). The distances between sources at locations 3-8 and detector A ranged between 2.5 and 3.6 cm. These light channels were used for measuring hemodynamic changes in the brain. The distances between sources at locations 3-8 and the detector B, as well as the distances between sources at locations 1 and 2 and the detector A, were 0.8 cm or 1.1 cm (see Fig. 1). These light channels provided information on the hemodynamic fluctuations in superficial tissue near the main detector A and near the main sources at locations 3-8. This information was used for the correction of the signals measured at large source-detector separations with respect to the contribution of superficial and systemic fluctuations (see Section 3). To ensure that the light intensities from all sources were within the dynamic range of the data acquisition card, the signal from detector A and the sources at locations 7 and 8 were attenuated using neutral density filters (Kodak). The optical sensor was attached to the left or right side of the subject's forehead above the eyebrow above the sinuses.
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 head coil. Sagittal T 1 -weighted localizer scans were used to determine the correct plane for the functional scans. Gradient-echo echoplanar images were acquired using a data matrix of 64 x 64 complex points, TR=1280 ms, TE = 40 ms, FOV = 240 mm, slice thickness = 7 mm, no inter-slice gap, receiver bandwidth 62.5 kHz, and flip angle 60 o . Multi-modality radiological markers (IZI Medical Products Corp, Baltimore, MD) were embedded into the optical sensor to facilitate correct orientation of the MRI slices with respect to the sensor and to enable recovery of the sensor orientation for data analysis.

Data analysis
The major absorbers of near-infrared light in the human tissues are oxy-and deoxyhemoglobin. Changes in the oxy-and deoxyhemoglobin concentrations (∆[O 2 Hb] and ∆ [HHb], respectively) can be related to the changes of the tissue absorption coefficient [ HHb are the extinction coefficients of oxy-and deoxyhemoglobin [11].
On the other hand, λ µ a ∆ can be related to the changes in the recorded near-infrared signal using solutions of the diffusion equation describing the transport of light in the highly scattering medium. The simplest solution of the frequency-domain diffusion equation is [12] 2 where c is the speed of light in the medium, U DC , U AC and φ are the intensity, modulation amplitude (AC) and phase of the photon density wave respectively, r is the distance from the source to the point of measurement, µ' s is the reduced scattering coefficient, and is the diffusion coefficient where the value of the parameter α may vary between 0.2 and 1 depending on the experimental conditions [13][14][15]. S 0 and A are the strength and the modulation amplitude of the light source, respectively. The value of the dimensionless parameter x=ω/cµ a is small at the modulation frequency ω of our instrument (110 MHz). The solution shown in Eqns. (2) is valid for an infinite homogeneous scattering medium and gives a good approximation to the exact solution for a homogeneous semiinfinite medium in the case of a large source-detector distance (r»D) [12]. Although the optical properties of the head are not homogeneous, Eqns. (2) may be used for qualitative estimation of the time course of changes occurring in the brain provided that there are no significant fluctuations in the extracranial tissue or if such fluctuations can be removed from the signal by means of some correction procedure. Given the small size of x and keeping in mind that for human tissues µ a «µ' s , one can see from Eqns. (2,b) and (2, . Also we can assume that relative variations in the scattering due to changes in the oxy-and deoxyhemoglobin concentrations are small compared to the relative variations in absorption. There are two reasons for that. First, due to the blood flow increase during functional activation no significant increase in the total blood volume occurs but rather exchange between the relative concentrations of HHb and HbO 2 [16,17]. Second, moderate changes in tissue blood volume cause stronger relative changes in absorption than in scattering at light wavelengths of 758 nm and 830 nm [18]. This is because at 758 nm and 830 nm the absorption is caused almost exclusively by hemoglobin, while the scattering is caused by the whole tissue, not only by the blood cells. Neglecting variation in scattering, change in the recorded signal can be related to small absorption change as: where AC σ , DC σ , and φ σ are the differential path length factors (DPF) for the signal modulation amplitude, intensity and phase, respectively [19].
In their numerical simulation study of the sensitivity of the frequency-domain near infrared signal to hemodynamic changes in different layers of the human skull Firbank et al. [4] showed that the signal is most sensitive to the changes in the superficial tissue adjacent to the source and detector but the sensitivity decreases rapidly with depth. In our previous NIRS studies of functional hemodynamics in the brain we found possible signal distortion by hemodynamic fluctuations in extracranial tissues [10]. These facts prompted correcting the long-distance optical signal with respect to fluctuations contributed by superficial tissue layers. For this we fitted every long-distance phase and AC signal from light channels between source locations 3-8 and detector A with a linear combination of all short-distant The average signals were calculated from the six corrected long-distance AC or phase signals and then used to calculate deoxyhemoglobin changes. The motivation for this procedure was that the subtraction of the multiple regression fit from a long-distance signal should eliminate components in its time course that are correlated with the short-distance signals, and the averaging of the corrected long-distant signals should emphasize the signal variations common for all long-distance signals. We expected this procedure to reduce both the contribution of superficial fluctuations and systemic changes, i.e. changes which were correlated in all layers of the head. One should note that the signal correction procedure using multiple regression fit does not distort the signals since, if there is no correlation, the regression coefficients are equal to zero. Additionally, the data were linearly detrended and low-pass filtered with a cut-off frequency of 0.2 Hz. This low-pass filtering excluded fluctuations of the time scale significantly shorter than the activation period, particularly fluctuations due to respiration [9]. It is known that the blood oxygen level dependent (BOLD) fMRI signal (the change in the intensity of the EPI images) is due to the changes in the deoxyhemoglobin concentration [10]. Namely, the increase of the BOLD signal corresponds to a decrease in [HHb]. The -∆[HHb] time series obtained from the AC and phase records with original acquisition period of 640 ms were locked to the fMRI acquisition period (1280 ms) and then used as an indicator function (a predictor) in the correlation analysis of the fMRI data.
For the analysis of fMRI data we used the MEDx 3.4 image processing package (Sensor Systems, Inc.). Every measurement produced 500 sets of six axial 2D EPI slices (64×64 matrix). with a repetition time of 1280 ms. The time series corresponding to the EPI image intensity (the BOLD signal) at each voxel were cross-correlated with the predictor (-∆[HHb] time series) to calculate the correlation coefficient. Values of correlation coefficients were then transformed to a z-score (which is a Gaussian distribution with zero mean and unit variance [20]). Then the test of statistical significance of the z-values was performed using an uncorrected p-value of 0.005. The Bonferroni correction was applied to reduce the probability of false significance [21]. The resulting critical z-value corresponding to the corrected probability threshold was used as a threshold in the final correlation maps. To obtain fMRI correlation maps showing activation due to the task a smoothed boxcar predictor function instead of -∆[HHb] time series was used. The maximum of this predictor function coincided with the middle of the activation periods, and the minimum corresponded to the middle of the relaxation periods. Correlation coefficients were calculated using the autocorrelation correction and assuming a hemodynamic response delay of 5 s. The assessment of statistical significance and thresholding was performed in the same way as when using the -∆[HHb] predictor.

Measurement protocol
The pairs of numbers separated by the "+" symbol were presented to the subjects using the fiber-optic goggles connected to a computer outside the magnet room. Every subject performed six 40s long calculation trials separated by the 20s periods of fixation. During fixation periods only the "+" symbol was presented. The beginning of the fMRI recording was triggered by the program presenting the stimulus. The beginning of the fMRI recording and every activation period were marked in the NIRS data to provide synchronization between fMRI and optical data records.
Seven healthy right-handed 20-38 year old subjects were screened. The optical sensor was attached to the left or right side of the subject's forehead above the eyebrow and the sinuses. The detector line was oriented parallel to the eyebrow. After the activation analysis measurements were repeated with those subjects who exhibited activation under the non-hairy area of the forehead with the optical sensor attached to the place where the activation occurred during the first measurement. The above protocol was approved by the Institutional Review Board (protocol number 01075.)

Results
Analyzing task-related brain activation from fMRI data we found that in all subjects the calculation stimulus caused most pronounced activation (in terms of z-score) in dorsal lateral areas of the frontal lobe, which were under the hairy skin and away from the optical sensor location. In some subjects activation was also detected in the lower part of the superior and middle frontal gyri under the non-hairy skin. However, during the second measurement the frontal activation was reproduced at exactly the same location only in two subjects. In one subject activation under the optical sensor was detected during the first measurement. Statistically significant correlation in the cortical tissue near the location of the optical sensor was detected using both intensity and phase data in all three subjects who showed activation in the area of NIRS measurement. Fig. 2 displays thresholded statistical images (z-score>4.8) showing statistically significant temporal correlation of the BOLD fMRI signal with the -∆[HHb] signal obtained from the intensity data ( Fig.2 (a)) and from the phase data ( Fig.2 (b)). Fig. 2 (c) shows the activated areas in the brain obtained by correlation of the fMRI signals with the smoothed boxcar predictor. The arrows in all subfigures point to the location of the center of the optical sensor. In all three figures the threshold z-score value is the same and corresponds to the critical value for the activation map (Fig. 2 (c)), which had the highest critical value of the three (determined from the statistical significance analysis). In Figs. 2 (a) and (b) one can see that almost the same group of voxels that shows correlation with -∆[HHb] signal obtained from the AC data also exhibits high correlation with the signal calculated from the phase (located within the circled area near the sensor location). Note also that these voxels belong to the circled activated area in Fig. 2 (c). Apart from the circled right frontal area, where the sensor was positioned, Fig. 2(c) also shows activation in the left prefrontal area and in the occipital area in the back of the brain. Although these areas are more activated than the one located near the optical sensor (see Fig. 2 (c)), they exhibit lower correlation with both phase-and AC-generated -∆[HHb] signals than the circled areas in Figs.  2 (a) and (b). This indicates that those may be the local cerebral hemodynamic changes including non-task-related fluctuations that contribute to the changes in the optical signals. In two subjects statistically significant correlation between BOLD and -∆[HHb] signals were detected in brain areas that showed no correlation between the BOLD signal and the smoothed boxcar predictor. In one of them correlation was found only for hemodynamic signal obtained from the AC data; in the other -in the signals calculated from both the AC and phase data. Fig. 4 shows correlation maps for the latter subject. As in Fig.2 the threshold z-score value is the same in all three figures and corresponds to the critical value for the activation map (Fig. 4 (c)). In Figs. 4 (a) and (b) one can see that cortical areas near the optical sensor location exhibit statistically significant correlation between the BOLD and -∆[HHb] signals, although the areas of correlation are much narrower than those in Figs. 2 (a) and (b) at close threshold z-score value (4.8 for Fig.2  Since the critical values of the z-score corresponding to the statistical significance for the maps showing correlation between BOLD and -∆[HHb] signals were slightly lower than those for the activation maps, the actual areas of statistically significant correlation between BOLD and -∆[HHb] signals were slightly wider than it is shown in Figs. 2 (a) and (b) and in Figs. 4 (a) and (b), but only due to the voxels exhibiting lower correlation z-score than those which are shown. Correlation was also detected in the slices adjacent to the slices shown in Figs. 2 and 4. One should note that these figures were obtained using corrected optical signals, so that Figs. 2 (a) and (b) and Figs. 4 (a) and (b) reveal correlation in specific local fluctuations. The correlation of -∆[HHb] signals with BOLD signals in some groups of voxels situated away from the sensor location could be due to correlation of local hemodynamic changes in different parts of the brain.

Discussion
We detected statistically significant correlation between BOLD and deoxyhemoglobin signals in five out of seven measured subjects. These five cases include all three cases when the cortical tissue adjacent to the optical sensor location was activated, so that a relatively large volume of the brain tissue exhibited strong hemodynamic fluctuations synchronous with the stimulus presentation. In two other cases activated brain areas were found some distance from the sensor location, so that only spontaneous fluctuations occurred near the optical sensor. In four cases, including one case without activation under the optical sensor, correlation was found in the hemodynamic signals obtained using both AC and phase data. This is the first time demonstration of the fact that the phase signal-to-noise ratio is high enough to detect cortical hemodynamic changes.
Calculating ∆[HHb] signal from the optical AC and phase data we assumed that hemodynamic fluctuations do not cause significant changes in the scattering properties of tissue. The collocation of the areas of correlation obtained using phase and AC data favors this assumption. On the face of it, this result seems to support the idea of quantitative measurement of cerebral hemodynamic changes using Equations (3) and (4), provided that the DPF values were measured. However, this idea does not seem to be practical because of the inhomogeneity of the head tissue. Although the fact that we found statistically significant correlation between -∆[HHb] and BOLD signals indicates that the DPF method may give correct time course of hemodynamic changes, the calibration of these changes is still an open question.
To obtain cerebral hemodynamic signals we subjected the original data to a correction procedure, which uses signals from the short-distance light channels (see Section 3). The fact that the resulting statistical images display highest correlation z-score in the voxels close to the location of the optical sensor indicates that our method is efficient for the exclusion of superficial and systemic hemodynamic fluctuations.
The fact that statistically significant correlation was detected even for spontaneous hemodynamic fluctuations indicates that NIRS can be sensitive even to relatively weak local fluctuations in the human brain. However, this is only a preliminary result, which should be confirmed by larger number of observations.

7.Conclusion
In five subjects we demonstrated correlation in hemodynamic signals obtained from the NIRS and BOLD fMRI data, which were simultaneously recorded during brain activation. In four subjects correlation was found for the hemodynamic signal obtained using the phase of intensity modulated NIRS signal. This is the first-time demonstration of the possibility of using the phase of the intensity-modulated optical signal for the detection of hemodynamic changes in the human brain. Our results show that the DPF method correctly reproduces the time course of the cerebral hemodynamic signals. The effect of superficial and systemic fluctuations on the overall optical signal can be reduced by a data correction procedure, which uses signals acquired at small source-detector distances (1 cm and less).