Measurement of brain activity by near-infrared light

. We review our most recent results on near-IR studies of human brain activity, which have been evolving in two directions: detection of neuronal signals and measurements of functional hemodynamics. We discuss results obtained so far, describing in detail the techniques we developed for detecting neuronal activity, and present-ing results of a study that, as we believe, conﬁrms the feasibility of neuronal signal detection. We review our results on near-IR measurements of cerebral hemodynamics, which are performed simultaneously with functional magnetic resonance imaging (MRI) These results conﬁrm the cerebral origin of hemodynamic signals measured by optical techniques on the surface of the head. We also show how near-IR methods can be used to study the underlying physiology of functional MRI signals. © 2005SocietyofPhoto-OpticalInstrumentationEngineers. 10.1117/1.1854673]


Introduction
There is a series of outstanding questions regarding the detection of brain activity in measurements made on the surface of the head and about the origin of the observed changes in the optical parameters of tissues. Experiments on exposed cortex have shown that neuronal activity produces measurable changes in optical parameters. For example, in 1973 Cohen 1 showed that the amount of light scattered from a slice of brain tissue rapidly changes following stimulation. Also, reflectivity changes in the near IR have been measured on exposed cortex. 2,3 The question is whether changes can be observed through the skull. Only recently, noninvasive observations of fast changes in optical parameters were reported. 4 -7 In these pioneering experiments, near-IR light was brought to the external surface of the head using a source optical fiber and the reflected light was measured at a distance of about 3 to 4 cm using a collection optical fiber. In this configuration, the light must travel deep ͑1 to 1.5 cm͒ into the head to reach the detector fiber ͑see Sec. 2͒. We showed with experiments on phantoms and with simulations that changes in the optical properties of a small region at a depth of about 1.5 cm from the head surface can affect the amount of light collected. These experiments showed that, in principle, it is possible to measure changes in optical parameters occurring in a deep region. The question remains whether we have sufficient sensitivity to measure changes due to neuronal activity.
Another topical application of the optical techniques to brain studies is the near-IR spectroimaging ͑NIRS͒ of meta-bolic processes. The main advantage of NIRS over other techniques used in this area is its outstanding biochemical specificity, which is based on the classical optical spectroscopic methods. Recently, many studies have demonstrated that cerebral hemodynamic changes associated with functional brain activity can be assessed non-invasively not only in infants, but also in adult human subjects. Several types of brain activity have been accessed, including motor, [8][9][10][11] visual, 12,13 auditory, 14 and cognitive 15,16 activities. Apart from purely spectroscopic measurements aimed at local detection of functional cerebral hemodynamic changes, several groups have performed imaging-type studies aiming at localization of cerebral hemodynamics. 17,18 Recently, Franceschini et al. 11 presented a first demonstration of noninvasive real-time imaging of cerebral activation in humans. We then reported the first NIRS mapping of functional hemodynamics in human motor cortex with simultaneous functional magnetic resonance imaging. 18 In this paper, we review the information generated so far and we discuss the available evidence concerning the origin of the neuronal and hemodynamic effects observed using NIRS. We also show how near-IR methods can be used to study the underlying physiology of functional magnetic resonance imaging ͑MRI͒ signals. Since the field is still controversial and rapidly advancing, we mainly focus on our own work and a few other studies that we choose to illustrate our opinions.
For the determination of the absolute values of the optical parameters, there is a general consensus among researchers that it is necessary to measure not only the light attenuation, but also the time of flight of the photons in the tissue to accurately recover the absolute values of the optical parameters. This can be achieved using two alternative methods: the time-domain method and the frequency-domain method. There are substantial technical differences between the two methods, but from a conceptual point of view, they provide equivalent information, being related through the Fourier transform. In this paper, we discuss the frequency-domain approach, because this is the method we are employing in our laboratory.

Frequency-Domain Method to Measure Optical Properties of Turbid Media
In 1989 Patterson et al. showed that the optical parameters of a turbid medium can be obtained from time-resolved measurements of short light pulses propagating in the medium. 19 A fit of the intensity as a function of time measured at some distance from the source provides the values of both the absorption and the reduced scattering coefficients of the medium. During the same period, Gratton et al. developed a method of measurement using intensity modulated light sources. 20 In this method, the light intensity is modulated at high frequency ͑in the 100-MHz range͒ and the dc ͑mean intensity͒, ac ͑modulated amplitude͒, and phase of the detected light is measured at a distance from the source after the light has traveled through the turbid medium. Since the light pulse used in the first method can be considered as the sum of the harmonic signals, both methods are mathematically related by means of the Fourier transform. However, in practice frequency-domain methods have a higher signal-to-noise ratio ͑SNR͒ and are faster than time-domain methods. The frequency-domain method was adopted by many investigators including Chance, Patterson, and Wilson ͑see Refs. 21-28͒. The speed and precision of the frequency-domain method are crucial for the detection and measurement of neuronal activity in the human brain.
In the range of source-detector distances used in our experiments ͑1 to 5 cm͒, the physical model based on the diffusion approximation to the Boltzmann transport equation is adequate to describe the propagation of light in highly scattering materials, such as tissues. The diffusion approximation is based on using the photon density, which can be interpreted as a quantity proportional to the probability of finding a photon at a distance r from the source at a time t. It can be shown 29,30 that at large enough distances from the light source and at frequencies below 1 GHz, the photon density obeys the diffusion equation whereas the diffusion coefficient depends on the reduced scattering coefficient s Ј of the medium and on the absorption coefficient a .
In the case of a harmonically modulated light source, solutions of the diffusion equation represent harmonic photon density waves that have a constant ͑dc͒ component and the oscillatory component ͑ac͒ spreading from the light source. 31 Since solutions depend on the absorption and scattering coefficients of the medium, they can be used to measure solutions. In our paper, to measure the absorption and scattering coefficients of homogeneous media, after the analysis of several possibilities, we found it more convenient to measure the values of the phase and ac at several distances from the source. 32 We then used only the slope of the plots of phase, ln(rU dc ), and ln(rU ac ) as a function of the distance to obtain the values of the reduced scattering and absorption coefficients. The details of the method are described in Ref. 33. The method provides separation of the scattering from the absorption with negligible crosstalk.
To implement this idea an instrument was built that operates simultaneously at two different wavelengths, 758 and 830 nm ͑ISS, Champaign, Illinois͒. At these wavelengths, the major contributions to absorption in human tissue are due to deoxyhemoglobin and oxyhemoglobin, respectively. In this instrument, the light from laser diodes is rf modulated at 110 MHz. Heterodyning at the detectors, using well-established techniques, 34,35 converts this rf frequency to 5 kHz. An analog-to-digital ͑A-to-D͒ converter card in a computer slot simultaneously samples the signal of the four detectors at 160 kHz. A computer program processes the stream of digitized data in real time. Exactly eight points are digitized per period at the 5-kHz frequency. After a preset number of periods of the 5-kHz frequency have been acquired, an electronic multiplexer circuit turns on the next light source. The multiplexing of the light sources proceeds indefinitely. The fast Fourier transform ͑FFT͒ routine performs the calculation of the dc, ac, and phase of the signal for all four detectors. Since the A-to-D sampling is synchronous with the 5-kHz frequency, the phase of each detected signal has a well-defined relationship with the phase of the light source. The result of data acquisition and preprocessing is a stream of dc, ac, and phase values at a rate up to about 200 Hz/light source, since the data from the four channels are interlaced. In principle, data acquisition can be faster if one uses a larger value of the heterodyning frequency. The obtained values of ac and phase are then used to calculate the optical properties of the tissue. Theoretically, one can use either ac or dc. In practice, however, ac is preferable because unlike dc it is not sensitive to ambient light. Since absorption can be measured free of the scattering contribution, an absolute estimation of hemoglobin oxygen saturation and concentration is easily achieved.

Noninvasive Measurement of Brain Activity: Fast and Slow Effects
Several authors have shown that functional changes in the brain tissues associated with neuronal activity can be detected by optical methods based on measures of propagation of diffusive light. 5,36 -41 Two types of changes have been reported: fast effects ͑with a 50-to 300-ms latency͒ and slow effects ͑with a 2-to 7-s latency͒. These two effects probably have different origins, which can be distinguished because of their large difference in time scale. Fast effects are presumably directly related to changes in neuronal tissue. Slow effects have a hemodynamic origin.
Slow and fast effects are measured in different ways. Slow changes in the tissue parameters are observed using a sustained activation task, such as one involving repeated stimulation ͑or repeated movements͒. Such an approach is used in optical, functional magnetic resonance imaging 42 ͑fMRI͒ and positron emission tomography ͑PET͒ studies. Fast effects have been detected by time-locked averaging the changes elicited by many individual stimulations. Apart from studies employing optical methods, this approach is also used in electrically recorded potentials ͑ERPs͒ and magnetic encephalography ͑MEG͒ experiments.
Probably, fast and slow effects differ in terms of the optical parameters they affect. Slow effects can be observed using both steady state and time-resolved methods. 36,37,[39][40][41]43 The fast effects were first detected using the time-resolved methods, 5,39,43 although in later works they were also successfully studied using steady state techniques. 44 -46 The best method to measure fast changes is still being researched.
To a large extent, the slow effects can be attributed to changes in absorption that are mostly produced by oxygenation and hemodynamic changes that occur in active areas of the brain. 37,47 They may also be related to the effects observed with PET studies of blood flow. Although the fast effects have been attributed to changes in scattering of neuronal tissue, 5 their nature is still being investigated.

Optical Detection of Neuronal Activity
At the early stage of our studies on the optical detection of fast neuronal effects, to investigate the ability of the optical method to detect fast changes in the brain, we performed experiments using modulated light with a wavelength of 715 nm, and a power of less than 1 mW. The source-detector distance was fixed 43 at 3 cm. The phase and the attenuation of the photon density wave were acquired using sampling rates varying between 12.5 and 50 Hz in different experiments. The measurements were obtained using a single-channel system ͑maps were obtained by repeating measurements at multiple locations͒. The pulsation of the vascular system, and in particular of large and medium arteries, which is associated with systolic activity of the heart, produces relatively large changes in the transmission of light through tissue. 48 The signal due to pulsation can be quite significant, which is exploited in studies of hemodynamic effects related to brain activity. However, when the interest is in detecting neuronal activity, pulsation may produce substantial artifacts that are difficult to eliminate with simple bandpass-filtering approaches. It is therefore important to separate the changes in intensity due to hemodynamics from the changes due to neuronal activity. An algorithm was developed for the estimation of the pulsation artifact 49 based on regression procedures in which the effect of systolic pulsation is estimated from the data. Summarizing, the algorithm extracts a mean shape of the pulsation by screening each trace separately for pulses. The period of each pulse is adjusted before averaging. This mean shape corresponds to the best estimate of each pulse and is used to remove each pulse from the data. This requires adjusting the period of the mean shape to the one of each pulse and scaling the shape by linear regression. Tests showed that, when used to estimate and subtract the effects of pulsation from the optical recordings, the procedure substantially reduces the impact of the pulsation artifact.
One of the major problems in the optical detection of neuronal signals is that signals measured using a single sourcedetector pair are most sensitive to changes occurring very close to the source or detector, while in the case of the brain, the region of interest is situated about 1.5 cm below the surface. Therefore, special measures should be taken to increase the sensitivity of the optical device to changes deeper in tissue and to reduce the effect of biological noise near the surface. In the following description, we assume that the detector is placed on the surface of the tissue at a distance d from the source. We define the region of sensitivity ͑ROS͒ of such a source-detector pair as that region of space that is traversed by the photons emitted by the source and collected by the detector. The ROS can be obtained by calculating for each volume element P the product of the probability of a photon emitted by the source arriving at P multiplied by the probability that the detector will measure a photon at P. This is known as the equation of the ''light bundle.'' The shape of the light bundle depends on the boundary conditions, tissue optical parameters, and in the frequency domain, on the frequency of light modulation. By considering the effect caused by a small particle at position P on the ac and phase, we can also define the ROS for the ac and phase signal. To this end, we denote by ac 0 ( 0 ) and ac P ( P ) the ac ͑and phase͒ measured without the particle and with the particle, respectively. The ac and phase ROS bundles are obtained by plotting the following fractions: 1Ϫac p /ac 0 and 1Ϫ p / 0 , respectively. Cross sections of the ROS of ''ac'' and ''phase'' bundles in the plane passing through the source-detector line perpendicularly to the surface are shown in Fig. 1. The bundles were computed using the theory based on the diffusion approximation. 30 Coordinate x corresponds to the direction along the surface, and coordinate y to the direction perpendicular to the surface. As Fig. 1͑a͒ shows, the ac bundle of the single source-detector pair has the highest sensitivity in the region just below the source and the detector. At the same time, the ROS is relatively weak midway between the source and the detector. The phase ROS ͓Fig. 1͑b͔͒ extends deeper into the tissue, but it is still quite sensitive in the regions close to the source and detector. Note that the phase ROS also strongly depends on whether the signal change arises from changes in the absorption or scattering ͓compare Figs. 1͑b͒ and 2͑b͔͒. We confirmed these properties of the ac and phase ROS experimentally. For this, we inserted a small absorbing object in the medium and measured the difference in the signal ͑dc, ac, and phase͒ caused by the introduction of the object ͑1-mm absorbing sphere͒ as we moved the object throughout the medium. The measured changes in the ac and phase signals as the functions of the object coordinates were in excellent agreement with the ROS shown in Figs. 1͑a͒, 1͑b͒, 2͑a͒, and 2͑b͒.
To reduce the sensitivity of the optical device to changes in the superficial regions under the source or detector fibers and to make it more sensitive to changes in deeper regions between the source and detector we developed the concept of the ''pi-sensor.'' In the following description, we use the term ''sensor'' to indicate a combination of light sources and detectors that can measure localized changes in the optical parameters of the medium. The basic idea of the pi-sensor is the simultaneous use of two or more source-detector pairs that produce ''light bundles'' in a ''cross'' configuration, as shown in the center panels of Fig. 1. The signal processing is designed to be sensitive only to changes that simultaneously occur in all source-detector pairs. A simple way to obtain such a ''correlator'' is to numerically cross-correlate the changes measured at different source-detector channels ͑the name of the ''pi-sensor'' was derived from the type of data processing employed to produce the sensor output͒. The ROS of the pisensor composed of two source-detector pairs for the absorbing inclusion is shown in Figs. 1͑c͒ and 1͑d͒ for the ac and phase signals, respectively. Note that in the ac ROS of the pi-sensor, the sensitivity is larger in the region under the cen-ter of the cross rather than immediately under the source or the light detector ͓compare Fig. 1͑c͒ with Fig. 1͑a͔͒. An even larger sensitivity to changes in the central region can be achieved by adding more source-detector pairs centered at a common point ͓see Figs. 1͑e͒ and 1͑f͒ and 2͑c͒ and 2͑d͒ for absorbing and scattering inclusions, respectively͔. The ac ROS of the pi-sensor for both absorbing and scattering inclusions has an average radius of about 3 mm when eight sourcedetector pairs are used with a semi-infinite geometry. The lower panel of Fig. 1 shows how the radius ͑at half maximum͒ of the ac ROS for the semi-infinite medium changes as the number of source-detector pairs is increased ͑the ROS has no spherical symmetry, the radius of the ROS corresponds to the average radius͒.
We used the pi-sensor together with a low-noise frequency-domain NIRS instrument and highly effective filtering algorithms to detect the fast neuronal signals from volunteers during motor stimulation and to determine the collo-cation with the slow hemodynamic changes. The protocol was approved by the Institutional Review Board of the University of Illinois at Urbana-Champaign. Written informed consent was obtained from all subjects prior to the measurements.
Optical fibers from four laser diodes ͑at either 758 or 830 nm͒ were bundled at each source location to achieve a higher light intensity, so that 16 sources of each wavelength were used. To distinguish the light from the four source locations, the laser diodes were multiplexed. The detector fibers had a 3-mm core diameter, the sampling rate was 96 Hz. The sensor had four sources and four detectors ͓photomultiplier tubes ͑PMTs͔͒ arranged equidistantly along the arc of a 3.2-cmdiam circle ͑16 light bundles in total͒ and was placed above the motor cortex ͑C3 position͒. The tapping frequency was set at 2.5 times the heart rate, which was measured for each subject at voluntary tapping conditions prior the experiment.
There was a sequence of alternating periods of 20 s tapping and 20 s rest for 5 min. In each subject, five measure- ments were carried out by tapping different fingers ͑index, middle, ring, fourth, and all fingers͒. To provide synchronization between the stimulus and the optical data we recorded the ''tapping signals'' ͑the S signals͒ synchronously with the optical data using the auxiliary input of the data acquisition system. These were signals corresponding to the contact between the metal caps put on the tapping finger tip and the thumb.
To reduce physiological noise from the arterial pulsatility the ac and phase data were filtered using an adaptive filter. 49 In addition, the data were detrended by a digital high-pass filter with the cutoff frequency equal to 2.2 times the mean heart rate and the cuton frequency about 2.4 times the mean heart rate. The cutoff frequency was selected to further reduce the effect of the hemodynamic pulsatility, but still to be low enough not to affect the fast neuronal signal at 2.5 times the mean heart rate.
The cross-correlation function ͑CCF͒ of the time lag between the optical data and the stimulus ͑S͒ signal was obtained as well as the CCF between the optical signals in different source-detector channels of the pi-sensor ͑pi-correlation͒. To calculate the CCF between two signals f (t) and g(t) we used the equation To determine how the level of noise depended on the number of data samples, the S signal was cross-correlated with the optical signals measured on a phantom block. To verify that the observed correlation was not due to an artifact, we compared the correlation for the data obtained during tapping ͑after filtering the noise͒ with the CCF for the data acquired on the phantom. The level of random noise was reduced as 1/N 1/2 , where N was the number of data points. In contrast, if the optical signal contained components, which were correlated with the stimulus, this genuine signal would maintain its amplitude independent of the number of data points. Using measurements on the phantom block we determined the noise reduction as a function of N. We also verified that this factor was the same for physiological measurements without stimulation.
Filtering the heartbeat was crucial for the detection of the fast signals. Figures 3͑a͒ and 3͑b͒ show the CCF for the simultaneously recorded stimulus and the unfiltered ac and phase signals, respectively. Comparing these CCFs with those for the filtered data ͓bold curves in Figs. 3͑c͒ and 3͑d͔͒, one can see that the curves corresponding to unfiltered data exhibit no correlation. The bold curves in Figs. 3͑c͒ and 3͑d͒ correspond to data acquired at the source-detector distance of 3.2 cm. Figures 3͑c͒ and 3͑d͒ also show the CCFs for the filtered signals acquired at a 1.3-cm source-detector distance ͑thin curves͒. One can see that these short-distance signals exhibit much lower correlation than the signals acquired at the larger distance.
Note that the period of the CCFs shown in Figs. 3͑c͒ and 3͑d͒ by the black curves is very close to the period of tapping presented using the metronome ͑about 3 s͒. This indicates that the correlation is not spurious but is due to the coherence between the optical signals and the tapping. Using the Fourier transform we verified that the mean tapping frequency corresponded to the one of the metronome. However, we found that the tapping signals recorded during different experiments ͑with the same or different subjects͒ were uncorrelated. To test whether the cross-correlation technique is sensitive only to the unique tapping signal corresponding to the given data set we obtained CCFs between the optical and stimulus signals recorded in different experiments ͑with the same subject͒. Such CCFs for the ''wrong'' stimulus signal are shown in Figs. 3͑e͒ and 3͑f͒ for ac and phase, correspondingly. One can see that the CCFs corresponding to the wrong tapping signals show significantly lower correlation than the CCFs corresponding to the tapping signal recorded simultaneously with the optical data ͓bold curves in Figs. 3͑c͒ and 3͑d͔͒. Figures 3͑g͒ and 3͑h͒ demonstrate the performance of the pi-detector. Fine curves in each figure show the CCFs of the S signal with the signals in the individual source-detector channels ͓ac in Fig. 3͑g͒ and phase in Fig. 3͑h͔͒. The bold curves show CCFs for signals from different channels. One can see that the pi-correlation of two single-channel signals shows the fast signal at the tapping frequency as clearly as the CCFs for single-channel signals and the S signal. The difference in the time lags of the peaks of different curves in Figs. 3͑g͒ and 3͑h͒ accounts for the fact that, as Figs. 1 and 2 show, the individual light channels have much broader sensitivity regions than the pi-detector. Neuronal responses in different regions of the motor cortex can exhibit different time lags with respect to the stimulus. Therefore, when the neuronal response is sampled over two large cortical areas, the overall signals may also have different time lags.
In principle, instead of CCF the block average can be used to detect signals correlated with the stimulus. However, it has lower SNR than the CCF. Since the neuronal response provides relatively small contribution to the overall optical signal, we chose CCF as a tool for neuronal signal detection. Figure 4 represents the back-projection map of the fast signal during tapping exercise. The map was obtained using the ac signals of a four-bundle pi-detector shown in the left panel, which was sequentially placed on different positions of the head in the motor cortex area. As Fig. 4 shows, the fast signal was highly localized and clearly showed individual spots of activity. On the contrary, we found no difference in the localization of the hemodynamic signal corresponding to different fingers.
Note how we concluded that the signals we observe were ''real'' and discuss the control experiments we performed.
1. If we do not remove the pulse signal, the crosscorrelation function between the raw signal and the tapping channel contains a relatively large contribution due to cardiac pulsation, but in some cases the neuronal signal can still be distinguished. However, after the pulse removal, the SNR increases dramatically ͓Figs. 3͑c͒ and 3͑d͔͒ and the neuronal activity can be easily detected.

If we perform a test experiment without tapping, there
is no cross-correlation signal and no neuronal activity could be detected. 3. If we measure optical signals in several locations in the region of the motor cortex, there is a definite spatial dependence of the strength of the cross-correlation signal ͑Fig. 4͒. 4. If we analyze light bundles that are quite superficial ͑source-detector separation distance 1.3 cm͒ the crosscorrelation signal is either absent ͑in the phase signal͒ or much reduced ͑by more than a factor of 10͒ with respect to deeper bundles ͓Figs. 3͑c͒ and 3͑d͔͒. 5. One possible artifact can result from poor statistics. We analyzed a set of optical data using the S signal from a different experiment ͑similar in tapping frequency͒. We found that there was no correlation in this case ͓Figs. 3͑e͒ and 3͑f͔͒. 6. The measurements were repeated four times on different days and with different tapping frequencies. In all cases, we found very similar results with good SNRs.

Simultaneous fMRI and NIR Mapping of Brain Activity
Recently NIRS has been proposed as a method to study brain hemodynamics. It is simple and inexpensive compared to ''heavy-duty'' methods such as fMRI and PET. Many studies applying NIRS to monitor functional cerebral hemodynamics have been reported 8 -16,39,40 ͑a review can be found in Ref. 41͒.
However, the precise location of tissues contributing to the optical signals was uncertain. This issue could be clarified using coregistration of brain hemodynamics by NIRS and by fMRI, because the blood-oxygenation-level-dependent ͑BOLD͒ signal employed in fMRI is believed to be proportional to changes in the deoxyhemoglobin concentration. We performed a study with the aim of comparing functional cerebral hemodynamic signals measured simultaneously by the optical method and by fMRI. The contribution of superficial layers ͑skin and skull͒ to the optical signals was also assessed. Both optical and fMRI methods were used to generate maps of functional hemodynamics in the motor cortex area during a periodic sequence of stimulation by finger motion. For NIR measurements we use a two-wavelength ͑758 and 830 nm͒ frequency-domain ͑110-MHz modulation frequency͒ Oximeter ͑ISS, Champaign, Illinois͒, which had 16 laser diodes ͑8 per wavelength͒ and two PMT detectors. At a wavelength of 758-nm light absorption by the deoxyhemoglobin ͑HHb͒ substantially exceeded absorption by the oxyhemoglobin (O 2 Hb), while at 830 nm the O 2 Hb absorption was dominant. The laser diodes operated in a sequential multiplexing mode with 10 ms of ''on'' time per diode. Light emitted by these laser diodes was guided to the tissue through 10-m-long, multimode, silica optical fibers. Two 10-m-long glass fiber bundles collected the scattered light and brought it to the detectors. The paired ͑758 and 830 nm wavelength͒ source fibers were attached to the probe at eight positions. Together with two detectors, they provided 10 biwavelength sourcedetector channels with a source-detector distance of 3 cm. The probe covered an area of 9ϫ6 cm 2 . The geometric layout of the optical probe and its position on the head are shown in Fig. 5. In this figure, the numbers 1 through 8 correspond to the light sources and the letters A and B to the detectors. Using this notation, 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. 50 Changes in the concentration of O 2 Hb and HHb were calculated from the changes in the light intensity at two wavelengths using the modified Lambert-Beer law. 41 Three multimodality radiological markers were 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.
MRI was performed using a 1.5-T whole-body MR scanner ͑Sigma, General Electric Medical Systems, Milwaukee, Wisconsin͒ 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 echo-planar images were acquired using a data matrix of 64ϫ64 complex points, temporal resolution ͑TR͒ϭ640 ms, echo time ͑TE͒ϭ40 ms, field of view ͑FOV͒ϭ240 mm, slice thicknessϭ7 mm, no interslice The middle slice ␥ was set between the skull and the brain surface at the C3 position. This slice mostly contained cerebrospinal fluid ͑CSF͒, dura mater, arachnoidal tissue, and pia mater, but could also reach the cortical tissue depending on the subject's anatomy. Two deeper slices, ␣ and ␤, mostly contained brain tissue, while two outer slices, ␦ and ⑀, included the skull, the skin and the markers. Figure 7͑a͒ represents the anatomical image corresponding to the deepest slice ␥ shown in Fig. 6. The numbers 1 through 8 show the source locations and the letters A and B indicate the detectors. 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 preexercise epoch, 10 20-s stimulation epochs, separated by 10 20-s control epochs, and a 50-s afterexercise epoch. During stimulation epochs, subjects performed light palm squeezing with the right hand following the rhythm presented by the sound from the vacuum pump of the MRI scanner. 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. The 750 images acquired during each run resulted in an 8-min scan time. The synchronization of the exercise sequence with the fMRI and optical records was provided by a computer program generating the commands for the subject. Using head landmarks for the C3 position, in most cases 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. 7͑a͔͒.
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. Usually, the largest size of the activated area was in the middle slice ␥ ͓see Fig. 7͑b͔͒. Figure 7͑b͒ shows a typical correlation map of the functional BOLD fMRI signal for the slice ␥ ͓the corresponding anatomical image is shown in Fig.  7͑a͔͒. The highest intensity white color in Fig. 7͑b͒ corre-sponds to the voxels with an activation-response correlation z score larger than 9.0. In the next slice deeper in the brain ͑slice ␤ in Fig. 6͒, which was mostly cortical tissue, the activation area was almost as large as in the middle slice ␥ shown in Fig. 7. The deepest slice ␣ typically showed a significantly smaller activation area. In two upper slices ␦ and ⑀, the coefficient of correlation between BOLD signal and the paradigm function was typically less then 0.1 and the corresponding z score was less then 2.0.
The arrow in Fig. 7͑a͒ points to the cluster of voxels ͑shown as a dark color͒ exhibiting the highest correlation ͑z score higher than 10͒ between the ͓HHb͔ signal from light channel A1 and the negative of the voxel BOLD signal. One can see that this cluster is situated very close to the middle of the area between the source 5 and detector A. This indicates that it is the tissue corresponding to this voxel cluster that contributes to the time course of ͓HHb͔ measured in the light channel A5. Figure 8͑a͒ shows a BOLD signal from the slice ␥ with r ϭ0.7 and a z score of 12.5. Figures 8͑b͒ and 8͑c͒ show the  ͓HHb͔ and the ͓O 2 Hb͔ signals, respectively, from the light channel A5 ͑see Fig. 6͒ acquired simultaneously with the BOLD signal shown in Fig. 8͑a͒. One can see that during the activation periods indicated by the gray bars, the BOLD signal shows the increase, while the ͓HHb͔ shows the decrease, and during the rest periods the BOLD signal decreases, while the ͓HHb͔ increases. Due to fluctuations in the upper tissue layers corresponding to slices ␦ and ⑀, in some cases the ͓O 2 Hb͔ and ͓HHb͔ signals from the activated area were noisier than the signals shown in Fig. 8. However, in all subjects the folding average analysis ͑see Fig. 9͒ revealed the same type of ͓O 2 Hb͔ and ͓HHb͔ behavior for light channels situated above the major activated area. Figure 9 shows a map of folding average ͓HHb͔ and ͓O 2 Hb͔ traces acquired simultaneously with the fMRI maps presented in Fig. 7. One can see that the light channels situated above the activated area ͑channels A1, A5, B1, and B5) exhibit a prominent decrease of ͓HHb͔ during the stimulation, which was concurrent with a significant increase of the oxyhemoglobin concentration. No significant decrease in the folding average ͓HHb͔ traces concurrent with the significant ͓O 2 Hb͔ increase was observed during stimulations in the light channels outside the activated area.
The decrease in ͓HHb͔ concurrent with the increase in ͓O 2 Hb͔ agrees with the results of a previous simultaneous fMRI-optical study, 51 BOLD signal theory, and with basic knowledge of brain physiology. 52,53 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 proves that NIR monitors the functional brain hemodynamics. Additional support for this conclusion is the temporal correlation between the intracranial BOLD signals and the respective NIR signals. We observed such a strong correlation in several 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 result indicated that the task-related hemodynamic changes measured by NIRS are not due to artifacts, but have an intracranial origin.

NIR Study of the Roles of Deoxyhemoglobin and Blood Volume Changes in the fMRI BOLD Signal
The results described above qualitatively confirmed that the BOLD fMRI signal is closely related to local changes in  ͓HHb͔. However, recently some theoretical and experimental arguments were put forward that not only ͓HHb͔, but also the cerebral blood volume changes can be important for the quantitative BOLD signal description. 53 The experimental argument was limited by the fMRI measurements of the blood flow changes, which were not directly and immediately linked to the volume changes. To better understand the roles of deoxyhemoglobin and blood volume changes in the fMRI BOLD signal, we decided to use frequency-domain NIRs simultaneously with BOLD fMRI. Unlike the previously described study on simultaneous fMRI and NIR mapping of brain activity, in this study we used the multidistant method described in Sec. 2. We selected this method because it enables us to measure not only relative changes in ͓HHb͔ and ͓tHb͔ϭ͓HHb͔ϩ͓O 2 Hb͔ ͑the total hemoglobin concentration, which is proportional to the blood volume V), but also to assess the background values ͓HHb͔ 0 and ͓tHb͔ 0 . These values are necessary because, according to theory, the BOLD signal is given 53

by a linear combination of the values Pϭ1
Ϫ⌬͓HHb͔/͓HHb͔ 0 and Rϭ⌬V/V 0 Ϫ1. During this study we found another advantage of the multidistance method, namely, that it was much less sensitive to fluctuations in superficial tissues than the method using only one source-detector pair for measuring hemoglobin concentration changes. We related this feature to the fact that the mean square spatial fit of the multisource data efficiently filters out uncorrelated fluctuations occurring near particular sources.
In this study, we used a small ͑1.5ϫ3.5-cm͒ optical sensor that had two detectors and 16 paired sources at eight locations ͑758 and 830 nm at each location͒. The distances between sources and detector ranged between 2 and 4 cm. The optical sensor was centered at the measured C3 position of the left hemisphere. The fMRI and NIRS instrumentation and the activation paradigm were the same as in the study already described: each experiment included 10 10-s stimulation epochs, which were separated by 10 17-s control epochs. During stimulation epochs subjects performed light palm squeezing with the right hand. The total length of one data record was about 5 min. The measurements were performed on four healthy right-handed male hairless volunteers, 18 to 57 years old. Informed consents were obtained from all subjects. To increase the SNR in the measured changes ͓HHb͔ 0 -͓HHb͔ and ͓tHb͔ 0 -͓tHb͔ ͑denoted as ⌬͓HbO 2 ͔ and ⌬͓HHb͔, respec-tively͒ we performed time-locked averaging of the traces related to the 10 stimulation-relaxation cycles. The BOLD signal in the activated voxels of echo planar image ͑EPIs͒ was also subjected to time-locked averaging. In order to assess the contribution of the deoxyhemoglobin and blood volume changes to the BOLD signal, we modeled the fMRI time series as the linear combination of Pϭ1Ϫ⌬͓HHb͔/͓HHb͔ 0 and Rϭ⌬V/V 0 Ϫ1: This equation is equivalent to the model of the BOLD signal proposed in Ref. 53. It represents the physical fact that the BOLD signal can depend on both cerebral blood volume and oxygenation, and a mathematical fact that small changes in these values should contribute to the BOLD signal approxi-mately as their linear combination. The details on the derivation of this equation can be found in Ref. 53. We computed the coefficients A P and A R using the bivariate regression analysis procedure, in which we considered P(t) and R(t) as the regressors for the BOLD signal. We also computed the correlation coefficients between the BOLD signal and the hemodynamic variables P(t) and R(t). In this study, we found that temporal variations in R were much less correlated with the BOLD signal than variations in P, whose time course was very close to the time course of the BOLD signal ͑see Table  1͒. The values of the correlation coefficient between P and BOLD signal in all measurements were very close to unity, indicating a high covariance between the changes in ͓HHb͔ and the negative of the BOLD signal. The values of the correlation coefficient between R and the BOLD signal were significantly smaller, and for one subject they were even negative. The ratio of regression coefficients A P and A R was always positive, indicating that an increase in the blood volume correlated with an increase in the BOLD contrast. However, the values of the r 2 statistic for the regression only increased slightly when the blood volume changes were included in the regression model ͑see Table 1͒. From this we concluded that the contribution of blood volume changes to the BOLD signal was very small compared to the contribution of the deoxyhemoglobin changes. This conclusion was also supported by the fact that the changes in R were significantly smaller than changes in P ͑see Table 1͒.

Discussion
Numerous studies show that neuronal activity can be detected optically. Basically, there are two lines of evidence, one based on optical studies of exposed cortex of animals, 2,3 and the other based on measurements through the skull. 5,6,[43][44][45][46] While the measurements on exposed cortex directly visualize areas of neuronal activation, the spatial resolution that can be achieved by imaging though the skull is insufficient to delineate the columnar areas of cortical activity. It is also likely that the optical contrast mechanism could be different for the two types of approaches due to the different ways light interrogates the tissue. In the direct exposed cortex measurement, the contrast results from the changes of reflectivity, due to either changes in absorption or scattering at the brain surface. 2,3 In the through-the-skull measurements, we must make some hypothesis about the path that the photons follow through the bone and brain and back out. Most of the controversy about the noninvasive optical measurements arises from the particular model used to describe light propagation through the head. While we have quite an accurate theory to describe light propagation in multiple scattering media, the problem of applying this theory to the intact head is that we need to know the details of local anatomy and of the optical parameters of each tissue element. Probably, the most successful approach has been that of considering only average tissue properties, and segmenting different tissues with a coarse resolution, of the order of a centimeter. 54 The theoretical justification for this approach should be sought in the measured values of the mean-free path for photons in the different tissues. This parameter strongly depends on the nature of the tissue. Under the hypothesis that we are in the multiple scattering regime, we can use the value of the reduced scattering coefficient as an indicator of the mean-free path and of the length over which we can average the optical properties. 30 For tissues in the head, this value is typically of the order of millimeters. This consideration has led to the hypothesis that we can divide various tissues into three major segments: the superficial skin layer, the skull, and the brain surface.
In a study on the visual cortex activation, 6 the observed fast changes of the optical properties suggested that this signal occurs simultaneously with, or before, the surface electrical potential ͑evoked potential or ERP͒. This feature makes the fast optical signal an ideal marker for indexing the occurrence of neuronal activity. The fast optical signal also appears to be very localized, enabling discrimination between different parts of the same neuro-anatomical area of the cortex ͑as in the case of Brodmann's area 17 ͒. These data indicate 7 that such localization is consistent with that obtained using fMRI. This may be useful in two ways. First, it may help us to study the relative timing of neuronal activity in different brain areas. This may have profound theoretical implications, since we expect most brain functions ͑such as vision, memory, attention, language, movement control, etc.͒ to depend on the dynamic interactions between different brain areas. Second, the possibility of acquiring data with high resolution both in the temporal and the spatial domains may help to integrate various noninvasive methods for studying brain function, and in particular electrophysiological and hemodynamic techniques.
Today the optical techniques for measuring functional cerebral hemodynamic changes in humans [8][9][10][11][12][13][14][15][16][17][18] are established much better than the methods to measure neuronal activity. However, even this class of potential applications is still in a developing stage. There are two different types of measurements that one can perform on the brain. The first attempts to measure the absolute optical properties of the brain tissue. 54 A second approach is concerned only with the measurement of the variations in the optical parameters. 10,11 Of course, for the measurements of the variations in optical parameters, the knowledge of the absolute value is less critical. Based on perturbation approach, it can be shown that the variations of optical parameters can be recovered much more accurately than the absolute values. For most of the brain functional studies, this is the approach that is normally followed. How-ever, for clinical applications it is of great interest to accurately determine the absolute values of the tissue optical parameters.
In spite of such an immature stage of the NIRS technology, it already yields valuable results. In particular, as we demonstrated in Sec. 5, the remarkable biochemical specificity of NIRS enabled us to investigate the physiological aspects of the functional magnetic resonance imaging of the brain. The main advantages of fMRI in brain mapping studies are high spatial resolution and the absence of any penetration limits. 42,55 The problem, however, is that the BOLD effect is a complex biophysical phenomenon, so that the BOLD signal may depend upon multiple physiological parameters. 53 At the University of Illinois at Urbana-Champaign ͑UIUC͒ we have an exciting opportunity to integrate the optical method with fMRI, which is the long-term goal of our research. Currently, we are very active in the development of a multimodal near-IR-fMRI instrument. Together with the optically detected neuronal hemodynamic changes, the combination of fMRI and NIRS can provide quantification that is difficult to achieve with fMRI alone.

Conclusion
Although changes in optical parameters ͑NIR reflectivity and scattering͒ have been measured in the exposed brain, the detection of these changes through the skull remained controversial. It was commonly accepted that neuronal signals obtained by currently available techniques are buried in the strong hemodynamic fluctuations and, in general, in instrument noise. Several researchers have concluded that studies over a population are necessary to extract the signal from the background noise. Our goal was to enhance the quality of the measurements to the level that a measurement on a single subject and on a relatively short time scale would have sufficient SNR to be detectable. We believe that we have now demonstrated that this goal can be achieved and we have understood some of the reasons why earlier studies were affected by large background noise. Note that, if we can routinely obtain the quality of the signal shown in Fig. 3, the optical technique could make a tremendous leap into the methods we use routinely for studying the functioning of the brain.
A remarkable biochemical specificity of the optical method enables researches to understand physiological mechanisms of complex metabolic processes, such as functional cerebral hemodynamics. Using NIRS we compared spatial distribution of hemodynamic changes with those imaged by fMRI during motor activation. A good collocation of the fMRI and NIRS signals confirmed the intracranial origin of the optical response to stimulation. We have also shown that, at 1.5 T, the change in deoxyhemoglobin concentration is the major factor determining the time course of the BOLD signal in the human motor cortex during motor stimulation. The influence of changes in regional cerebral blood flow ͑rCBV͒ is qualitatively in agreement with recently established theory, 53 but it is much smaller than that of changes in the deoxyhemoglobin concentration. The increase in cerebral blood oxygenation during functional activation is due to an increase in the rCBF velocity, and occurs without a significant swelling of the blood vessels.
Further integration of NIRS and MRI techniques promises a significant mutual enhancement of both by combining a remarkable spatial resolution of MRI with outstanding temporal and spectral resolutions of NIRS.