2.1.1.

Concurrent resting state study

A concurrent fNIRS and fMRI resting state study was conducted on six healthy volunteers (four male and two female, aged 28.0 ± 4 years). A MRI-compatible fNIRS probe was placed on the right side of the prefrontal area of each subject (roughly between Fp1 and F7 according to the 10–20 system) to collect fNIRS data while he/she was scanned by a Siemens TIM Trio 3T scanner (Siemens Medical Systems, Malvern, Pennsylvania). The probe was designed to collect data from short and long source-detector pairs (1.0, 3.0 cm) at prefrontal cortex, where a well-known resting state network exists [Fig. 1a]. The total scan time was about 6 min and 30 s with 1.5 s repetition time (TR). The acquisition rate of fNIRS was 12.5 Hz. The fNIRS probe had three detector fibers and two source fibers with source-detector distances at 1 or 3 cm [Fig. 1a].

Fig. 1

Schematic diagrams of the optical probes and their placement in the (a) resting state and (b) finger tapping studies.

2.1.2.

Concurrent finger tapping study

A block design finger tapping task was performed by a healthy subject (male, aged 30 years). The task included six periods of right-hand tapping alternated with rest periods (20 s tapping and 20 s rest). The MRI-compatible fNIRS helmet was worn by the subject during the experiment. The helmet consisted of four detectors and 16 sources arranged in a configuration that covered both hemispheres [Fig. 1b]. The acquisition rate of fNIRS was 6.25 Hz, while the TR was 2.02 s.

2.2.1.

Resting state data analysis

In the case of resting state studies, four steps were required: 1. raw fNIRS time courses (690 and 830 nm data) were converted into two time courses representing Δ[HbO] and Δ[Hb] according to the differential path length factor method²¹ using Matlab (The Mathworks, Natick, Massachusetts); 2. the time courses of Δ[HbO] from C1 (due to its high SNR for all subjects) and its temporally shifted versions (61: −7.2 to +7.2 s with steps of 0.24 s) were then anti-aliased and down-sampled to the fMRI acquisition frequency of 0.67 Hz (1/1.5 s), which had 260 data points; 3. each down-sampled time course was used as a regressor in fMRI data analysis, a total of 61 analyses were performed; and 4. the resulting 61 z-statistic maps were combined and the maximum z value was chosen for each voxel. The final maximum z-statistic map would include all of the voxels that had high correlation with fNIRS at certain time points.

2.2.2.

Finger tapping experiment and data analysis

In order to study the contralateral and ipsilateral brain activation under the finger tapping task, the NIRS probe was designed to have the maximal coverage of motor cortex areas on both sides [Fig. 1b]. Both NIRS and fMRI data were simultaneously collected. A similar RIPTiDe procedure was applied on the finger tapping data with the following modifications: 1. The time course of Δ[HbO] and its time-shifted versions (21 time shifts: −2.4 to +2.4 s with steps of 0.24 s) from each source-detector pair were used as regressors in the analysis. In total, 21 analyses were performed for each source-detector pair; 2. In order to separate the neuronal motor activation from the physiological signal in fNIRS and BOLD fMRI data, the block diagram of the finger tapping periods was convolved with the hemodynamic delay function. The convolved timecourse and its derivative were utilized as the first two regressors in the BOLD fMRI analysis, followed by the delayed fNIRS regressor. The derivative is included to compensate for variations in the hemodynamic response and slice timing. In this case, we were interested in the localization of the source of the physiological signal that overlays the neuronal signal in fNIRS. By using the task timecourse as well as the data from fNIRS as the regressors, the convolved finger tapping timecourse and its derivative were treated as confounds to regress out any task-related portion of the fNIRS regressor, leaving the corresponding z-statistic maps of the fNIRS regressor to show only the significant response to the physiological signal. Figure 2 shows three regressors for one particular channel (C3) directly over the activated area as an example: 1. the finger tapping task (boxcar function) convolved with the hemodynamic function (double gamma function offered by the FSL package), 2. the temporal derivative, and 3. the resampled time traces of Δ[HbO] for C3.

Fig. 2

Three regressors used in the analysis of the fMRI data obtained during a finger tapping task: 1. block diagram; 2. derivative of the block diagram; 3. resampled temporal trace of the Δ[HbO] from C3. Shaded areas indicate the tapping periods.

3.1.1.

Temporal traces

Figure 3 shows the resampled time traces of Δ[HbO] and Δ[Hb] measured at C1 and B1 on the foreheads of six subjects in the resting state studies. The resampling was required so that the fNIRS data could be used as regressors in the GLM analysis of the fMRI data. The anti-aliasing filter used in resampling applies a low pass filter (low pass = 0.3 Hz due to TR = 1.5 s).¹⁸ Therefore, the heart beat signal is not visible in Fig. 3. Prominent features of the data are: 1. the Δ[HbO] signal is much bigger than that of Δ[Hb] as seen in Figs. 3a, 3b, which is probably due to the high oxygenation of the blood throughout and 2. the phase of Δ[HbO] is not always 180 deg (anticorrelated) from that of Δ[Hb]. In fNIRS studies, a 180 deg phase difference between Δ[HbO] and Δ[Hb] obtained by fNIRS is commonly observed in the activated brain areas. However, the relationship between Δ[HbO] and Δ[Hb] in fNIRS resting state studies is not fully understood;²² and 3. the Δ[HbO] measured from B1 (1 cm separation) in the same subject is highly correlated [as seen in Figs. 3c, 3d] with that from C1 (3 cm separation). The same is true for the Δ[Hb]. The average correlation coefficients are 0.82 and 0.80, respectively. This indicates that a significant proportion of the fNIRS signal recorded by both source-detector pairs (which are collinear but with different distances) is the same, which in turn means that either the source of most of the signal is within the range of both probes or the source of B1, which is at a shallow region and source of C1, which is at a deep region are highly correlated.

Fig. 3

Resampled temporal trace of the Δ[HbO] (solid black line) and Δ[Hb] (dotted black line) obtained during resting state from six subjects calculated from (a) C1 and (b) B1. To demonstrate the similarities between time traces of Δ[HbO] and Δ[Hb] measured in C1 and B1 for all subjects, (c) resampled temporal traces of the Δ[HbO] in B1 (gray) and C1 (black) were plotted for every subject, while in (d) resampled temporal traces of the Δ[Hb] in B1 and C1 were plotted.

3.1.2.

Spatial distribution

In order to better locate the source of the fNIRS signal (in C1 and B1), we conducted voxelwise correlation between the fNIRS signals and the BOLD fMRI using RIPTiDe.¹⁹ Because Δ[HbO] has a much higher signal to noise ratio (SNR), and is demonstrated to be highly correlated with cerebral blood flow,²³ the results of using Δ[HbO] were the focus of the discussion. Figure 4 depicts the maximum z-statistic results of six subjects and the average of all subjects using Δ[HbO] from C1 [Fig. 4a] and Δ[HbO] from B1 [Fig. 4b]. The coronal (in the middle panel) and axial (in the right panel) views of Figs. 4a, 4b were chosen to show the prefrontal area underneath the probe (the markers of the probe is visible in both graphs). In the result corresponding to C1 [Fig. 4a], the voxels underneath the probe were highly correlated with the fNIRS signal as expected. However, the activation patterns are global and symmetrical throughout the brain (the probe was on the right side only). The highest correlation was present in the voxels resembling the superior sagittal sinus, the inferior sagittal sinus, and straight sinus as displayed in the left panels of Fig. 4a. Detailed discussion and graphs on the similarities of the correlation patterns of the subjects and their venous system (mainly on the surface of the cortex) were presented in our previous work,¹⁸ in which the same resting state data was analyzed. In short, a 3D-rendered image of all of the activated voxels was similar to the phase contrast image (MRI technique to map the blood vessel only) and the passing of the activated voxels throughout the brain followed the same passage of the cerebral blood flow and matched the time of blood circulation. These findings indicate that part of the Δ[HbO] signals obtained by fNIRS probes (C1) might reflect the changes in the vein and venules on the surface of the brain. The fact that the brain area underneath the fNIRS probe and its contralateral equivalent were prominently visible in Fig. 4, implies that the neuronal signal (e.g., resting state network) may also be present in the fNIRS data. The results corresponding to B1 [Fig. 4b] are, to a large extent, similar to the ones shown in Fig. 4a, as the symmetrical venous system. This is not surprising regarding the similarities between Δ[HbO] from C1 and B1 [seen in Fig. 3c]. However, for most of the subjects, smaller regions of activation underneath the probe were observed for B1 in Fig. 4b, which is even more obvious in the averaged results, shown at the bottom of Fig. 4. To further quantitatively assess the difference in the activated areas (obtained from B1 and C1), we defined a ROI (on the averaged head of six subjects) underneath the probe that covers a big area from the surface of the cortex to 2.5 to 3.0 cm deep into the brain [Fig. 5a]. All of the activated voxels (z > 2.3) observed by B1 and C1 within this ROI were counted for each subject. The ratios (number of voxels from B1/C1) were calculated and showed in Fig. 5b. From Fig. 5b, we can see that only subject 4 has the ratio bigger than 1 (1.11 to be exact), indicating that a slightly bigger area was observed for B1. For the rest of the five subjects, the ratios are smaller than 1 (ranging from 0.10 to 0.88), which confirmed that the smaller regions of activation underneath the probe were observed for B1. This indicates that most likely the sensitivity of B1 is lower in the cortex than that of C1, but it may nevertheless reach the surface veins.

Fig. 4

Maximum z-statistic maps for six subjects in resting state studies were obtained by employing the RIPTiDe for Δ[HbO] from (a) C1 and (b) B1. The z-statistic maps were registered onto the averaged anatomical head images of the subjects. The averaged z-statistic maps for all subjects are shown in the bottom. The markers of the fNIRS probe can be seen in the green circles.

Fig. 5

(a) 3D rendered average head with ROI (gray blob) chosen underneath the NIRS probe. (b) Ratios between activated voxels (within the ROI) observed by B1 and that by C1 for each subject.

3.2.1.

Temporal correlation

A detailed analysis of the neuronal response to finger tapping in fMRI and fNIRS can be found in our previous work.¹⁷ In this study, we will primarily focus on the physiological signal detected by fNIRS as opposed to the neuronal activity. Figure 6 shows the time traces of Δ[HbO] and Δ[Hb] for each channel surrounding four detectors. First, task activation was detected in C3 and C7, which were located over the motor cortex contralateral to the tapping hand. On the ipsilateral side, task-related deactivations were observed in A3 and A7 with negative z-values (z = −4.4 and −3.7 obtained from the correlation with task regressor). The reason why the finger tapping task evoked “activation” and “deactivation” on the contralateral and ipsilateral side of the brain is unclear and more subjects are needed to confirm the phenomenon. However, even in these task-related channels (C3, C7, A3, and A7), the physiological signal was nevertheless visible on top of the activation/deactivation patterns. In the other channels where task activation was not obviously visible, the physiological signal became the dominant signal. To further understand the internal relationship between the physiological signal detected by all of the channels, a high pass filter was applied to all 20 channels that are shown in Fig. 6, with low cutoff frequency at 0.025 Hz to filter out the neuronal effects of the finger tapping task (1/40 s = 0.025 Hz). The resulting time traces of Δ[HbO] and Δ[Hb] are shown in Figs. 7a, 7c, respectively. The correlation coefficients among these channels are illustrated as tables in Figs. 7b, 7d. Figure 7a shows that the temporal courses of Δ[HbO] from A1 to D7 (except A3 and A7) are highly correlated with each other despite their distributed (spatial) locations. This is further demonstrated by the correlation coefficient table in Fig. 7b in which the average correlation coefficient (except A3 and A7) is 0.69. This global correlation implies a strong underlying physiological fluctuation that affects most channels. The temporal traces of Δ[HbO] in A3 and A7 are correlated to each other with a correlation with coefficient as high as 0.93. Figure 7c displays the timecourses of Δ[Hb] from A1 to D7. Unlike Δ[HbO], in which the global correlation was observed, the location of the highly correlated channels of the Δ[Hb] are more (spatially) clustered. For example, high correlation was observed in the following groups: 1. A3 and A7; 2. B4, B5, B6, B7, and C1; and 3. C7 and D7. The average correlation coefficient is only 0.37.

Fig. 6

Resampled temporal traces of Δ[HbO] (solid line) and Δ[Hb] (dotted line) obtained in (a) finger tapping task by detector A and its surrounding five sources (indicated by the numbers next to the temporal traces), (b) detector B, (c) detector C, and (d) detector D. Shaded areas indicate the tapping periods.

Fig. 7

High pass filtered resampled temporal traces of (a) Δ[HbO] and (c) Δ[Hb] from all the channels in the finger tapping task. The correlation coefficient table of temporal traces of (b) Δ[HbO] and (d) Δ[Hb]. Shaded areas indicate the tapping periods.

3.2.2.

Spatial location

In order to localize the global signal observed in almost every channel [Fig. 7a], we applied the modified RIPTiDe for each channel as described in Sec. 2. The depiction of the Δ[HbO] from channel C3 (in Fig. 2) was chosen to demonstrate that when used as a regressor in conjunction with the task regressor, the portion highly correlated with the finger tapping signal was essentially removed from the fNIRS regressor, allowing separate detection of the effect of the task and the physiological signal superimposed on the task response (removing the correlated portion of the two regressors). Figures 8a, 8c, 8e show the corresponding z-statistic results of channel C3 in the brain areas activated by the task in Fig. 8a and the brain areas highly correlated with the physiological signal of fNIRS in Fig. 8c and their axial slices in Fig. 8e. The corresponding averaged z-statistic results from all 20 channels were calculated and overlayed onto the structure brain of the subject and rendered into 3D images depicted in Figs. 8b, 8d. The axial slices of the averaged results are also illustrated in Fig. 8f.

Fig. 8

z-statistic maps of the finger tapping task block diagram registered on (a) the subject's anatomical head and (b) its 3D rendered image with marked fNIRS probe. The maximum z-statistic maps obtained by RIPTiDe corresponding to the physiological fluctuation of C3 registered (c) on the subject's anatomical head and (e) its axial images. (d) 3D rendered images of the averaged maximum z-statistic maps of 20 channels (corresponding to the physiological fluctuation) and (f) their axial images.

The activation of the finger tapping task was clearly observed in both Figs. 8a, 8b as expected. More interestingly, the voxels correlated with the global physiological fluctuations, shown in Figs. 8c, 8d, 8e, 8f, are symmetrically spreading throughout the brain as visible in Fig. 4. Since we did not perform the whole head fMRI scan for the subject during the finger tapping task, only the middle portion of the brain is depicted in Figs. 8c, 8d, 8e, 8f. Nevertheless, the results show that areas that are highly correlated with the physiological signal seem to overlay the cerebrovascular system, especially the veins and venules as discussed in Sec. 3.1.2 of the spatial distribution of the resting state. This is further supported by the dynamic pattern of physiological activation generated by time-shifted NIRS regressors that follows the course of the circulating blood (data not shown).