Fusion of fNIRS and fMRI data: identifying when and where hemodynamic signals are changing in human brains

In this study we implemented a new imaging method to fuse functional near infrared spectroscopy (fNIRS) measurements and functional magnetic resonance imaging (fMRI) data to reveal the spatiotemporal dynamics of the hemodynamic responses with high spatiotemporal resolution across the brain. We evaluated this method using multimodal data acquired from human right finger tapping tasks. And we found the proposed method is able to clearly identify from the linked components of fMRI and fNIRS where and when the hemodynamic signals are changing. In particular, the estimated associations between fNIRS and fMRI will be displayed as time varying spatial fMRI maps along with the fNIRS time courses. In addition, the joint components between fMRI and fNIRS are combined together to generate full spatiotemporal “snapshots” and movies, which provides an excellent way to examine the dynamic interplay between hemodynamic fNIRS and fMRI measurements.


INTRODUCTION
Functional near infrared spectroscopy (fNIRS) offers unsurpassed high temporal resolution and provides quantitative hemodynamic information for both oxyhemoglobin (HbO 2 ) and deoxyhemoglobin (HbR), which plays an important role in the in vivo study of cognitive processing in the human brain (Jobsis, 1977;Cope and Delpy, 1988;Hoshi, 2003;Singh et al., 2005;Huppert et al., 2009;Ye et al., 2009;Yuan et al., 2010a,b;Brunno et al., 2011;Egetemeir et al., 2011;Gagnon et al., 2012;Yuan, in press). Similar to its fNIRS counterpart, functional magnetic resonance imaging (fMRI) is also a non-invasive imaging method that measures the hemodynamic responses to even-related neural activity with excellent spatial resolution and low temporal resolution. To take advantages of the complementary information from these two imaging modalities, two broad methods for fNIRS and fMRI/MRI integration have been developed for different clinical cases: (1) "spatial constraint," in which spatial information from fMRI/MRI images are utilized to aid diffuse optical imaging based on fNIRS measurements (Carpenter et al., 2007;Ferradal et al., 2013); (2) "temporal correlation," where fMRI bold signals are processed to generate the correlation with HbO 2 and HbR concentration changes converted from fNIRS recordings (Cui et al., 2011;Tak et al., 2011;Gagnon et al., 2012). However, so far the linking between HbO 2 /HbR signals and fMRI spatial maps has not been extensively investigated, which represents one of the main challenges for fNIRS and fMRI fusion. Therefore, it is crucial for us to develop new imaging techniques to reveal the connections between these two measurements so that we are able to examine the dynamic interplay between space and time of hemodynamic responses with high spatiotemporal resolution.
Interestingly joint independent component analysis (jICA) method has been developed to compute the linked temporally independent event-related potential (ERP) components and spatial independent fMRI components, which enables inferences to be made using estimated associations between fMRI sources and ERP electromagnetic sources (Calhoun et al., 2006a,b;Sui et al., 2009). However, what has not been tried is a joint estimation of the temporal parts of fNIRS waveforms and the spatial maps of fMRI images. In this study, the jICA is extended to identify the spatiotemporal decompositions composed of fMRI spatial components indicating where the hemodynamic signals are changing and fNIRS components indicating when the hemodynamic signals are changing. The fNIRS-fMRI fusion method will involve calculating for given stimuli, the connection between the timelocked fNIRS waveforms and fMRI activation maps for all participants or different measurement sections from a single subject. It is anticipated that the results derived from jICA could be visualized by computing spatiotemporal "snapshot," which provides an effective way to examine the dynamic interplay between fNIRS and fMRI hemodynamic sources.

METHODS
In terms of jICA, the joint spatial and temporal independences of fNRIS and fMRI are assumed to satisfy the following generative model for the data (Calhoun et al., 2006a,b), in which X fNIRS is the group data from the chromophore concentration change of HbO 2 or HbR for n subjects/n sections of a single subject and X fNIRS = FIGURE 2 | jICA decomposition of HbO 2 and fMRI joint data for right finger tapping tasks: three components were found to be significantly correlated with finger tapping tasks [Independent components (ICs) 2, 4, and 5] while two were well-correlated with the physiology and body movement noise (ICs 1 and 3). Each of these components is shown in a separate panel in the figure: IC2 (A), IC4 (B), IC5 (C), IC1 (D) and IC3 (E).
The fMRI maps are thresholded at |Z | > 1.5 for display purposes.
We employ the infomax ICA method for jICA of Equation 3, which utilizes a gradient ascent iteration algorithm to maximize the entropy of the output of a single layer neural network (Bell and Sejnowski, 1995). The resulting updated equation for the algorithm to calculate the shared unmixing matrix W (i.e., the inversion of A), the fused independent fNIRS sources u fNIRS and fMRI sources u fMRI is as follows, is the non-linearity in the neural network. The initial value for W, W(0) is a matrix composed of random vectors. The jICA method generally doesn't provide us the details on how the components between fNIRS and fMRI interact with one another. To achieve this, the spatiotemporal "snapshots" of the most significant components are generated in two ways (Calhoun et al., 2006a,b). First we need to calculate the linear combination of the fMRI components weighted by their joint fNIRS time courses for a specific point in time. If the n spatial (fMRI) and temporal (fNIRS) components are given by S = s 1 s 2 . . . s n , and T = t 1 t 2 . . . t n , where t i is a T × 1 vector containing the fNIRS time points and s i is a V × 1 vector represents the V brain voxels, the fMRI movie is computed as M fMRI = |T| × S T . It is noted the absolute value is utilized here because the joint components are fused using a single parameter. And a change in

BEHAVIOR TASKS AND fNIRS-fMRI RECORDINGS
The fNIRS tests are implemented with a block design for a right finger tapping task. The experiment is performed with a 24chnnel fNIRS system (Oxymon MKIII, Artinis), which has 8 sources, 8 detectors, and 24 channels. In this system, two continuous wave lights at wavelengths 781 and 856 nm are emitted at each source fiber. In the case of block design for right finger tapping tasks, the onset time for the first trigger was at 42 s, then followed by a 21 s period of activation alternated with a 30 s period of rest. This was repeated 10 times for the subject. As such, the total recording time was 552 s. During the task period, subject was instructed to perform a finger flexion and extension action repeatedly. Data segmentation, which is also known as epoching in signal processing, is utilized to chop up the continuous fNIRS data into small time periods. The general way to do this is to extract segments surrounding the event codes from the experiments, e.g., from −35 s prior to the event onset until 35 s after the event code in this study. The original photon density datasets could be downloaded from (http://bisp.kaist.ac.kr/NIRS-SPM/Sample_data). The converted HbO 2 and HbR measurements (the unitless differential path length factor DPF = 4; sampling rate: 9.75 Hz) are filtered, segmented and plotted in Figures 1A,B, respectively. The configurations of 24 channels located on the scalp are also provided in Figure 1C. We found channels 8-12 are very unique because their locations are close to the left motor cortex so that they are able to reveal the hemodynamic responses from both the right finger tapping stimuli and other events. As a result, signal averaging is implemented for the datasets from these five channels to generate the representative HbO 2 / HbR measurements with reduced noise. In particular, five-section (runs 2, 4, 6, 8, and 10) measurements are selected from the representative ten-run fNIRS recordings for further fusion analysis since they show the best signal noise ratio during the stimulus processing.
For fMRI recordings, the fiber length was 10 m to connect the optodes in the MR scanner to the NIRS instrument in the MR control room. A 3.0 T MRI system (ISOL, Republic of Korea) was used to measure the BOLD response. During the blocked task paradigm, the echo planar imaging (EPI) sequence was used with TR/TE = 3000/35 ms, flipangle = 80 • , 35 slices, and 4 mm slice thickness (Ye et al., 2009). The fMRI data could be downloaded partly from the same website for the same subject or requested from the KAIST lab. Based on the datasets provided, we generated five mean MRI images which are considered as five-section fMRI measurements (Calhoun et al., 2001;Ye et al., 2009). Then we will combine the five-section averaged fNIRS recordings in temporal domain and five mean fMRI images in space to extract the joint independent components for hemodynamic fusion. It should be noted here the mean fMRI images are not generated from the raw fMRI images at each TR. As such, the fMRI image for each section may not accurately match its correlated fNIRS measurements in terms of recording times, which may have some influence on the final results. However, the influences should not be that significant since we use mean data from each section for fusion. The five-section averaged fMRI image and HbO 2 / HbR results are presented in Figure 1D for references.

RESULTS OF fNIRS-fMRI FUSION
The computed joint components are provided in Figures 2, 3, in which the spatial components and regions of brain activity from the fMRI maps are plotted on the right while the correlated temporal joint components along with the mean HbO 2 / HbR signals are given on the left. We can see from Figures 2, 3 that the time courses of temporal HbO 2 / HbR components correspond well to different peaks on the mean HbO 2 / HbR Frontiers in Human Neuroscience www.frontiersin.org October 2013 | Volume 7 | Article 676 | 5

FIGURE 5 | HbO 2 components and fMRI "snapshots": on the left of each window is shown a linear combination of the fMRI maps that are well-correlated with physiology and body movement noise (composite ICs 1 and 3), weighted by the HbO 2 part of the components at a specific point in time.
On the right of each window is shown the estimated HbO 2 components that are correlated physiology and body movement noise. The time courses for IC1 (in blue) and IC3 (in green) are also plotted on the right of each window. Such a display provides a dynamic way to visualize how the noise affects the brain activity at different time signals. Interestingly, when put together the fMRI maps and fNIRS signals, we found the five spatial joint components correlate very well with the five temporal joint components. For example, the first negative peak of HbO 2 in Figure 2A basically identifies the physiology noise such as visual activity and eye blinks. Also visible for this waveform of the second temporal joint component in Figure 2A is a late positive peak. During the positive peak, the hemodynamic activity is mainly found in the left motor cortex. Then the joint components from HbO 2 / HbR and fMRI are combined together to generate full spatiotemporal movies to more clearly display the dynamic interplay between the fNIRS and fMRI measurements. In particular, two movies are created to show the spatiotemporal dynamics of the fMRI and HbO 2 fusion, in which Movie 1 displays the hemodynamic responses of right finger tapping events (composite spatiotemporal joint components 2, 4, and 5) while Movie 2 basically reveals brain activity from physiology and body movement(composite joint components 1 and 3). The spatiotemporal hemodynamic changes derived from fMRI and HbR fusion are also recorded by Movies 3, 4 for right finger tapping tasks (composite joint component 1 and 3) and physiology and body movement noise (composite joint component 2, 4, and 5), respectively. Figure 4 shows seven "snapshots" snipped from Movie 1 produced from fused HbO 2 and fMRI data for finger tapping tasks. On the left of the "snapshots" a linear sum of the fMRI maps weighted by their respective HbO 2 time courses are provided while the mean HbO 2 waveforms combined the estimated HbO 2 components are plotted on the right of each "snapshot." Figure 5 plots five "snapshots" captured from fused the fMRI and HbO 2 data, which basically reveals the spatiotemporal responses of the physiology and body movement noise. Likewise, the "snapshots" from Movie 3, 4 produced from fused fMRI and HbR components are given in Figures 6, 7 for finger tapping tasks and physiology and body movement noise, respectively.

DISCUSSION
We show, for the first time, a spatiotemporal reconstruction of the human brain's responses with high temporal and spatial resolution by fusing together the fNIRS signals and fMRI images. The jICA is a very key and unique technique for connecting hemodynamic fNIRS temporal signal and hemodynamic fMRI spatial information according to the joint constraints of temporal and spatial independences. Such direct correlation further provides a very valuable tool to localize non-invasively those high resolution structures that underlines temporally well-resolved fNIRS responses. Meanwhile, deep brain structures, which are not readily detectable with scalp fNIRS alone, can now be identified by the fusion technique when aided by the correlated spatial components from fMRI. It is observed from Figure 2 that the first and third HbO 2 -fMRI joint components basically identifies body movement and negative eye blink activities in pre-frontal lobe while the second, fourth and fifth HbO 2 -fMRI joint components are fused with hemodynamic activity in left primary motor cortex (PMC), supplementary motor area (SMA), and motor association cortex such as parietal cortex. In particular, the fifth joint component basically identifies the neural activity in the SMA and parietal cortex. These observations are further validated by Movies 1, 2, which show the dynamic interplay between space and time of HbO 2 -fMRI hemodynamic responses. It is noted Movie 1 displays the spatiotemporal changes of neural activation patterns for finger tapping tasks. Further, we can see from Figure 3 that the second, fourth, and fifth HbR-fMRI joint components basically reveal body movement and other physiology noise while the first and third HbR-fMRI joint components identify the patterns of brain activity in the left motor cortex and motor association cortex. In particular, the first HbR-fMRI joint component reveals brain activity in SMA and parietal cortex. These spatiotemporal findings from fused HbR-fMRI data are further processed to generate Movies 3, 4, in which we can see the hemodynamic changes for a finger tapping task with high temporal and spatial resolution, and we can also see clearly how physiology and body movement noise corresponds to the neural activity. Importantly, we can observe from Figures 2, 3 as well as Movies 1, 3 that strong brain activity mainly occurs in the left PMC, which validates stimuli in the right finger tapping task will yield the activations in the left PMC. Alternatively, we can see from Figure 4 the estimated HbO 2 time courses (i.e., linear combinations of the significant HbO 2 components) at specific positions in the cerebral cortex. The first "snapshot" in Figure 4A is associated with physiology and body movement noise while the second to fourth "snapshots" show the HbO 2 concentration change during onset and stimulus processing of right finger tapping, which are validated by the spatial maps of fMRI image on the left. We can further observe from Figure 4 that the fifth "snapshot" captures the maximized spatiotemporal hemodynamic responses while the sixth and seventh ones show the decreased trends of neural activity. So if we examine the different "snapshots" in Figure 4, the spatiotemporal dynamics of the right finger tapping are revealed. In addition, the negative and positive peaks of the estimated HbO 2 time courses in Figure 5 directly correspond to the physiology and body movement noise though partly they capture the weak and negative hemodynamic responses of right finger tapping in the case of strong body movement.
The HbR-fMRI "snapshots" in Figure 6 identify SMA and parietal cortex with early peak hemodynamic responses. Later responses occurring in left PMC show peaks of the HbR responses during stimulus processing. In particular, "snapshots" 2-7 display the hemodynamic changes in detail with high spatiotemporal resolution for finger tapping task. Figure 7 provides the portions of estimated HbR time courses at specific positions of fMRI images, in which we observe the peaks basically corresponding to the activities of physiology noise, eye blinks, and body movement.
In summary, we have used a combination of hemodynamic fNIRS and fMRI data to visualize the neural activation patterns with high spatiotemporal resolution. The fusion techniques have shown the potential by using the joint hemodynamic data to identify unique neural information that cannot be revealed in either technique alone. In particular, we implemented a joint decomposition of fNIRS and fMRI data, which was linked or fused by a common mixing parameter. The present method does not involve the solution of inverse problem for diffuse optical imaging or involve the use of the threshold for the fMRI data. However, we do need to decompose the combined data into specific components by using jICA, which are composed of fNIRS and fMRI portions. Separating the data into joint components will provide us a useful way to examine component specific differences among different patient groups or different stimulus tasks.
We have to point out the present investigation has significant limitations though we are encouraged that the results we have found seem meaningful and valid. For example, we only validate data from a single subject even without accurate measurements from fMRI. In future we would like to examine our method to incorporate different subject groups and complex neural stimuli as well as multiple time points acquired from multiple channels of fNIRS systems.