Functional mapping of the human brain with near-infrared spectroscopy in the frequency-domain Angelo Sassaroli*a, Yunjie Tonga, Francesco Fabbria, Blaise Frederickb, Perry Renshawb, and Sergio Fantinia a Department of Biomedical Engineering, Bioengineering Center,Tufts University, 4 Colby Street, Medford, MA 02155 b McLean Hospital, 115 Mill St. Belmont, MA 02478 ABSTRACT Comparison of the spatial and temporal information retrieved from near-infrared phase and average intensity (DC) data reveals that these data types can play a complementary role in the study of the temporal and spatial features of the optical response associated with brain activation during a finger-tapping test. The optical data have been collected with a frequency-domain tissue imager at two wavelengths (690 and 830 nm) and have been analyzed using standard filtering and folding-average procedures. DC and phase data show different temporal and spatial features. A plausible explanation of the different behavior of DC and phase data has been attempted by using Monte Carlo simulations. Keywords: Near-infrared spectroscopy, frequency-domain, brain imaging, Monte Carlo. INTRODUCTION Significant research efforts have been recently devoted to functional imaging of the brain with noninvasive optical methods in the near-infrared [1]. In the frequency-domain, near infrared spectroscopy (NIRS) employs intensitymodulated illumination and phase-sensitive detection and provides average intensity (DC), amplitude (AC), and phase data. Theoretical [2] and tissue-phantom studies [3] indicate that DC and phase data are sensitive to different regions within a turbid medium, so that they can play a complementary role in monitoring optical changes in tissues. However, there has not yet been a clear demonstration of the effectiveness of combining DC and phase data for functional imaging of the brain. One problem of using phase data is that it is typically associated with a lower signal-to-noise ratio with respect to DC data. In this work, we have collected DC and phase data on the head of a human subject during a finger-tapping test by using a helmet devised for the spatial mapping of the changes in the concentrations of deoxy-hemoglobin ([Hb]) and oxy-hemoglobin ([HbO]) during brain activation. In particular, we studied the temporal traces of DC and phase measured by the source-detector channels most sensitive to the brain activation. The spatial and temporal behavior of DC and phase data revealed different features that might suggest a way to use both data types for better localizing the perturbation (i.e. the hemodynamic change) associated with brain activity and its evolution with time. MATERIALS AND METHODS The experiments were carried out with a frequency-domain tissue imager (Imagent, ISS, Inc. Champaign, IL) comprising sixteen laser sources at 690 nm, sixteen laser sources at 830 nm and four optical detectors (photomultipliers tubes, Hamamatsu Photonics R928). The average output power of the lasers was about 1 mW. For brain mapping, the lasers were multiplexed typically at a frequency of 6.25 Hz so that each laser was on for 10 ms and two laser sources were on at the same time. The laser diodes were coupled to optical fibers (0.04 cm in core diameter and 10 m long), that deliver the light to the head of the subject. The optical signal was collected at the head of the subject by means of four optical fiber bundles (3 mm in diameter and 10 m long) that deliver the collected light to the detectors. A special helmet was designed to hold the source and detector fibers and to fit comfortably on the head of each subject. Figure 1 shows the locations of the optical fibers on the head. The distance between each source and its closest detector is fixed at 3.5 cm. P * angelo.sassaroli@tufts.edu; phone 1 617 627 4321 371 23,31 22,30 24,32 A B 21,29 17,25 20,28 19,27 6,14 detector 18,26 7,15 8,16 source 5,13 D 4,12 C 3,11 1,9 2,10 Fig. 1 View from above of the source and detector fibers arrangement for brain mapping. Sources 1&17, 2&24, 3&23 etc. at 690 nm are simultaneously on, respectively, and a full cycle through all 690-nm light sources takes 80 ms. The second turn-on cycle is completed for the sources at 830 nm (9&25, 10&32, etc.). The distance between each source and the closest detector is fixed at 3.5 cm. The finger tapping protocol consists of the following: a) Two minutes of baseline acquisition. In this phase the patient is asked to lie comfortably on a chair, while staying motionless and relaxing. b) A series of 7-8 activation/rest periods, each made of 20 s of activation (finger tapping) followed by 20 s of rest. The data were analyzed by using standard algorithms to filter out unwanted frequencies prior to applying folding average procedures [4]. RESULTS An example of raw DC and phase temporal traces collected on a human subject at 830 nm and their fast Fourier transform (FFT) is shown in Fig. 2. As we can see from Fig. 2, the DC data has strong frequency components at the heart rate and at the respiratory frequency (peaks at 1.3 Hz and 0.3 Hz in the FFT graph, respectively), which on the contrary are absent in the phase data. We investigated the possible origin of the different effects of these time-varying contributions on the DC and phase data by using a three-layer, frequency-domain Monte Carlo (MC) simulation [5] in a plane geometry. We chose the thickness of the first two layers (3 mm and 5 mm, respectively) to be representative of the adult scalp and skull, respectively, and the third layer (the brain) was considered to be infinitely extended. The optical properties of the layers (0.1, 0.1, and 0.1 cm-1, respectively, for the absorption coefficient; 8, 10, and 15 cm-1, respectively, for the reduced scattering coefficient) were chosen on the basis of values reported in the literature. Amplitude (AC), DC and phase data were calculated from large statistics of detected photons (60,000) for diffuse reflectance at 3.5 cm from the injection point. We modulated the absorption coefficient of one or more layers with a sinusoidal wave having amplitude of about 1% of the baseline value and a frequency of 1 Hz in order to simulate the perturbation due to the arterial pulsation at the heartbeat frequency. In the third layer we always introduced a modulation of the absorption coefficient (to represent cortical activation) having an amplitude of 1% of the baseline value and a frequency of 0.025 Hz (1/40 s-1). For the case of spatially uniform (i.e. occurring in all three layers) absorption perturbations due to the arterial pulsation, the results of the Monte Carlo simulation (not shown) indicate that both DC and phase data have spectral components at the heart rate. The fact that this is typically not the case for experimental data (see Fig. 2) suggests a nonuniform distribution of the absorption changes associated with the arterial pulsation. We have considered a non-uniform case, where the absorption modulation at the heartbeat frequency (1 Hz) only occurs on the top layer (scalp), while the second layer (skull) has a constant absorption coefficient (0.1 cm-1), and the third layer (brain) has an absorption coefficient modulated by the tapping/rest sequences (0.025 Hz). The results of the Monte Carlo simulation under these conditions are shown in Fig. 3. As we can see (Fig. 3, bottom panels) the phase data are not sensitive to the 1-Hz absorption changes occurring in the upper layer, which instead are clearly visible in the DC data. The agreement between the results of this simulation and the experimental data suggests that the optical perturbation at the heartbeat frequency is mostly originated in the superficial, extra cerebral layers of skin and scalp. 372 Temporal traces FFT 4 0.15 3.5 FFT (arbitrary units) 0.2 Δ DC/DC 0.1 0.05 0 -0.05 -0.1 -0.15 3 2.5 2 1.5 1 0.5 0 -0.2 50 100 150 200 250 300 0 0.5 1 1.5 2 frequency (Hz) 2.5 3 0 0.5 1 1.5 2 frequency (Hz) 2.5 3 time (sec) 0.25 4 FFT (arbitrary units) 0.2 Δ Phase (rad.) 0.15 0.1 0.05 0 -0.05 -0.1 -0.15 3.5 3 2.5 2 1.5 1 0.5 -0.2 0 -0.25 50 100 150 200 250 300 time (sec) Fig. 2 Experimental traces of ΔDC/DC = (DC(t)-DC0)/DC0 (top left) and ΔPhase = Phase(t)-Phase0 (bottom left) as a function of time, and their FFT (DC: top right, Phase: bottom right, respectively). DC0 and Phase0 are average values calculated during the baseline acquisition. The vertical lines indicate the start of the 20-s finger tapping task. The wavelength and the source-detector distance are 830 nm and 3.5 cm, respectively. The temporal traces are normalized to the integral of the squared values. Temporal traces FFT 2 1.95E-05 1.8 1.94E-05 FFT (arbitrary units) 1.6 -2 DC (cm ) 1.93E-05 1.92E-05 1.91E-05 1.90E-05 1.4 1.2 1 0.8 0.6 0.4 1.89E-05 0.2 0 1.88E-05 0 10 20 time (sec) 30 0 40 0.2 0.4 0.6 0.8 1 1.2 1.4 1 1.2 1.4 frequency (Hz) 43.8 2 1.8 1.6 FFT (arbitrary units) Phase (deg) 43.75 43.7 43.65 43.6 1.4 1.2 1 0.8 0.6 0.4 43.55 0.2 43.5 0 0 10 20 time (sec) 30 40 0 0.2 0.4 0.6 0.8 frequency (Hz) Fig. 3 Temporal traces of the DC (top left) and phase (bottom left) calculated by the MC code; the right panel shows the FFT of the temporal traces for the DC (top) and phase (bottom). The MC results refer to a three-layer medium with the following properties: s1 = 3 mm, μ’s1 = 8 cm-1, s2 = 5 mm, μ’s2 = 10 cm-1, s3 → ∞, μ’s3 = 15 cm-1. The absorption coefficients of the three layers are: μa1 (t) = [0.1+0.001sin(2πf1t)] cm-1, μa2(t) = 0.1 cm-1, μa3 (t) = [0.1+0.001sin(2πf2t)] cm-1, where f1 = 1 Hz and f2 = 0.025 Hz. The source-detector distance is 3.5 cm, and the relative and absolute refractive indices are n =1.4. 373 B A B D3 C D A C 4 detector source detector source (a) (b) Fig. 4 Temporal traces of ΔDC/DC (top) and Δphase (bottom) collected at the source-detector location D-3 (panel (a)) and D-4 (panel (b)). The shaded area indicates the finger-tapping period. source source 3 cm a) μa1 b) detector μa2 3 mm 3 cm μa1 detector μa2 μa3 3 mm absorption perturbation r = 2 mm 7 mm μa3 absorption perturbation r = 2 mm 7 mm Fig. 5 Geometry used for the MC simulations. The medium includes two outer layers having thicknesses of s1 = 3 mm and s2 = 7 mm, respectively, and absorption coefficients of μa1 = 0.07 cm-1 and μa2 = 0.1 cm-1, respectively. The radius of the outer spherical boundary is 8 cm and the absorption coefficient of the inner region is μa3 = 0.15 cm-1. Spherical perturbations having a radius of 0.2 cm are added in the inner region (on the plane defined by the center of the spherical geometry, the source, and the detector) at a distance of 1.21 cm from the outer boundary. In (a) a single sphere is initially placed symmetrically with respect to source and detector, and in the following steps pairs of identical spheres are added on each side of the existing spheres until the detected signal does not change noticeably. In (b) one sphere is initially placed to the right of the detector, and then one sphere at a time is added to the left side of the previous sphere(s) until the detected signal saturates. The absorption coefficient of the spheres is μas = 0.3 cm-1. The MC simulation was run for a uniform reduced scattering coefficient μ’s = 7 cm-1. A tapping-task experiment on a different subject was carried out with a data acquisition setting where 10 consecutive multiplexing cycles were averaged for each AC, DC and phase data point. A non-recursive time-domain filter implemented in Matlab was used on the experimental data to cut off the lowest (< 0.0125 Hz) and the highest (> 0.1 Hz) frequencies in the raw data. A standard folding average procedure enhanced the change in DC and phase data due to the finger tapping task. Spatial and temporal mapping of the activated area was carried out by a simple back- 374 projection algorithm [4]. Figures 4(a) and 4(b) show the activation traces collected at 830 nm from two different sourcedetector channels (D-3 and D-4, respectively, see Fig. 1), after the filtering and folding average procedures. We notice that the DC and phase temporal traces relative to these two source-detector pairs share some common features, but also display distinct differences. In particular, despite the common decrease of both data types during the activation, we observe that the onset of response in the phase is delayed with respect to the DC in channel D-3 (Fig. 4(a)), while the onset of decrease is more synchronous in channel D-4 (Fig. 4(b)). The decrease in both the DC and phase data during brain activation indicates an increase in the absorption coefficient. To explain the different relative behavior of the phase and DC at channels D-3 and D-4, we resorted again to MC simulations. In this case, we implemented a MC code for a spherical geometry, which can accommodate an arbitrary number of layers and spherical absorbing defects. The reduced scattering coefficient was assumed uniform (μ’s = 7 cm-1) across the various layers. A schematic representation of the MC geometry is shown in Fig. 5. The spherical defects are added according to two different schemes as shown in Fig. 5. In a) initially one sphere is placed half way between source and detector; and then two identical spheres are added in each step to the left and to the right of the first sphere until the detected DC, AC and phase saturate. In b) one sphere is initially placed to the right of the detector, in a region where it is barely detectable; then one sphere at a time is added to the left side of the previous one until the detected DC, AC and phase saturate. Consistent with experimental data, we fixed a threshold of detect ability for DC and phase of 0.5% of the baseline value and 0.1º, respectively. The results of these MC simulations are shown in Fig. 6. We notice that for case a) one spherical defect is enough for the DC to reach the detect ability threshold, while three spherical defects are necessary for the phase. For case b) two spherical defects are necessary for the DC signal and three spherical defects for the phase. In other words, when the perturbation starts from a location sandwiched by the source and the detector, the phase signal should lag behind the DC signal more than in the case when the perturbation starts closer to the detector (or source) location. Similar effects are also found when the perturbations are out of the plane defined by source, and detector (results not shown), and this might explain the experimental results shown in figure 4, where the activation most probably started in the region D3 and spread out over region D4. The different information content of the phase and DC data is also evident in the spatial maps of the optical changes. Figure 7 shows two maps of ΔDC/DC and Δphase at times t = 8 s and t = 21 s, respectively, after the onset of stimulation. At early times, the phase response appears more localized (to channel D-4) with respect to the DC response over channels D-3 and D-4. This result is consistent with the delayed phase response in channel D-3, as reported in Fig. 4(a). 0 -0.05 Δ DC/DC -1.00E-02 Δ Phase (deg) 0.00E+00 -5.00E-03 CW -1.50E-02 threshold -2.00E-02 -2.50E-02 -0.1 -0.15 PH (fi=0) threshold -0.2 -0.25 -3.00E-02 -0.3 -3.50E-02 0 0 3 6 9 12 15 18 21 3 6 24 a) 9 12 15 18 21 24 number of spheres number of spheres 0.00E+00 0 -5.00E-03 -0.05 Δ phase (deg) Δ DC/DC -1.00E-02 cw -1.50E-02 threshould -2.00E-02 -0.1 -0.15 -2.50E-02 -3.00E-02 ph threshould -0.2 -0.25 -3.50E-02 0 2 4 6 8 number of spheres 10 12 -0.3 0 2 4 6 8 number of spheres 10 12 b) Fig. 6 The results of the MC simulation for the case a) and b) of figure 5 are shown in the top and bottom row respectively. On the left (right) panel ΔDC/DC (Δphase) is plotted against the number of spherical defects added in the inner region. 375 t = 8 sec -0.015 t =21 sec -0.006 -0.038 t = 8 sec -0.14º -0.0043 t =21 sec 0.05º -0.24 º 0.03 º Fig. 7 Spatial maps of ΔDC/DC (top row) and Δphase (bottom row) at times t = 8 sec (left) and t = 21 sec (right) after the onset of stimulation by right hand tapping. The imaged area corresponds to the left brain hemisphere (detector channels C and D, as shown in Fig. 1). CONCLUSION In this study, we have investigated the possibility of combining phase and DC data in near-infrared, frequency-domain imaging of the brain for the characterization of the spatial and temporal features of the activation-induced hemodynamic changes. The temporal traces and their Fourier transform showed that the phase and DC data have different features. For example, the arterial pulsation and the modulation due to respiration are typically seen only in the DC data. On the basis of Monte Carlo simulations, we found that this result is consistent with a dominant contribution of the arterial pulsation to the optical data that originates from superficial tissue layers (skin, scalp) rather than from deeper brain tissue. We found that the DC and phase responses to brain activation also present different spatial/temporal features. At some locations, the phase response lags the DC response, while in other locations the DC and phase responses are synchronous. Monte Carlo simulations suggest that the spatial expansion of the absorption perturbation associated with brain activation, and the specific location of the activated area with respect to the illumination and collection points may account for this latter result. Further studies will be focused on the development of an algorithm to use both DC and phase data for the localization and temporal characterization of the hemodynamic response to brain activation. 376 ACKNOWLEDGEMENTS We thank Maria Angela Franceschini for her contribution in the collection of the data shown in Fig. 2. We acknowledge support from the National Institutes of Health (Grant No. DA14178) and the National Science Foundation (Award No. BES-93840). REFERENCES 1. 2. 3. 4. 5. See, for example, the Special Section on Optical Imaging in Psychophysiology 40(4), 487-571 (2003). M. Firbank, E. Okada, D. T. Delpy, “A theoretical study of the signal contribution of regions of the adult head to near-infrared spectroscopy studies of visual evoked responses,” Neuroimage 8 69-78 (1998) S. Fantini, M. A. Franceschini, S. A. Walker, J. S. Maier, and E. Gratton, “Photon path distributions in turbid media: Applications for imaging,” Proc. SPIE 2389, 340-349 (1995). M. A. Franceschini, V. Toronov, M. E. Filiaci, E. Gratton, and S. Fantini, “On-line optical imaging of the human brain with 160-ms temporal resolution,” Optics Express 6 49-57 (2000). M. Testorf, U. Osterberg, B. Pogue, and K. Paulsen, “Sampling of time- and frequency-domain signals in Monte Carlo simulations of photon migration,” Appl. Opt. 38, 236-245 (1999). 377