|Home | About | Journals | Submit | Contact Us | Français|
Functional near infrared (fNIR) imaging was used to identify spatiotemporal relations between spatially distinct cortical regions activated during various hand and arm motion protocols. Imaging was performed over a field of view (FOV, 12 x 8.4 cm) including the secondary motor, primary sensorimotor, and the posterior parietal cortices over a single brain hemisphere. This is a more extended FOV than typically used in current fNIR studies. Three subjects performed four motor tasks that induced activation over this extended FOV. The tasks included card flipping (pronation and supination) that, to our knowledge, has not been performed in previous functional magnetic resonance imaging (fMRI) or fNIR studies. An earlier rise and a longer duration of the hemodynamic activation response were found in tasks requiring increased physical or mental effort. Additionally, analysis of activation images by cluster component analysis (CCA) demonstrated that cortical regions can be grouped into clusters, which can be adjacent or distant from each other, that have similar temporal activation patterns depending on whether the performed motor task is guided by visual or tactile feedback. These analyses highlight the future potential of fNIR imaging to tackle clinically relevant questions regarding the spatiotemporal relations between different sensorimotor cortex regions, e.g. ones involved in the rehabilitation response to motor impairments.
Functional near infrared (fNIR) spectroscopy detects changes in light absorption and scattering in tissue caused by changes in concentration of oxyhemoglobin (HbO) and deoxyhemoglobin (Hb) secondary to neuronal activity, occurring through neurovascular coupling mechanisms that are still under investigation . With the use of multiple source-detector pairs over a specified field of view (FOV), images can be reconstructed by mapping the detected signals to the Green’s solution of the linearized perturbation diffusion model for tissue-light interactions [2,3]. This method is referred to as fNIR imaging. Though only being able to image hemodynamic activity in the cortex, its portability and lower cost compared to established radiological imaging modalities have motivated its use for the monitoring of a wide range of neurological conditions [4–10]. Of increasing interest in recent years has been the use of fNIR imaging to monitor cortical plasticity in response to motor skill rehabilitation [5,8] as an alternative technology to functional magnetic resonance imaging (fMRI) [11–14]. Though fNIR imaging has lower spatial resolution compared to fMRI, it can detect hemodynamic change patterns with higher sensitivity and temporal resolution [15,16], and enables researchers to perform longitudinal studies more easily. Additionally, fMRI requires subjects to be immobilized, since it is sensitive to motion artifacts [17–19], whereas fNIR imaging is relatively robust to motion artifacts allowing subjects to perform a larger range of motions [17,18] and thus offers more ways to probe sensorimotor cortex function.
These promising features of fNIR spectroscopy technology are tempered by the limited number of source-detector pairs that current brain imaging systems can accommodate. This limitation constrains users to either measure sparse activation data over a large area over the cortex that result in no images , or in low resolution images [21–23], or perform high resolution imaging studies over a limited FOV size . More specifically, in a previous study a large FOV of 13 x 13 cm was used to image the primary motor cortex (M1), primary sensory cortex (S1), premotor cortex (PMC), and supplementary motor area (SMA) of both hemispheres during the learning of a rotary arm movement . The minimum source-detector distance was 3 cm, and the spatial resolution of this study was limited due to the sparse source-detector geometry. Other studies with similar source-detector separations, but with smaller FOVs (6 x 6 cm, 3.75 x 2.5 cm, and 5.6 x 8.4 cm) focusing on M1 and S1, have been reported to image cortical activation due to finger tapping, forceful pinching, sequential tapping, and palm squeezing protocols [21,25,26]. In an attempt to further improve spatial resolution, a study limiting its FOV to 3.75 x 3.00 cm used bifurcated fibers separated by 0.75 cm to image cortical activity of finger tapping and vibrotactile stimulation of individual fingers . Though in most activation studies only the function of the S1 and M1 cortices are probed by finger tapping [1,15,25–28], palm squeezing , and vibrotactile stimulation , other surrounding cortical regions contributing to the planning and execution of motions are left outside of the chosen FOV.
In this work we used an array of bifurcated source-detector fiber bundles to image cortical activation over a FOV, spanning the PMC/SMA, M1, S1, and posterior parietal cortex (PPC) of a single brain hemisphere. Cortical activity was caused by finger tapping, tactile stimulation, squeezing of a stress ball (palm squeezing), sequential finger tapping, and card flipping (supination and pronation). This work presents, to our knowledge for the first time, fNIR images over an extended FOV (12 x 8.4 cm) in a single brain hemisphere for a number of different hand and arm motion activation protocols. The unique spatiotemporal relation of the PMC/SMA, M1, S1, and PPC activation patterns for each type of movement is demonstrated. We also demonstrate that all activation patterns can be measured without motion artifacts and reproducibly. In addition, application of a cluster component analysis (CCA) algorithm  to these fNIR signals shows that depending on the motor task different cortical regions share similar temporal activation patterns. Though these results are limited to one side of the brain due to the limited number of source-detector fibers available, they showcase the rich information that can be obtained from spatiotemporal activation patterns when a larger cortical area is imaged. Importantly, these results provide a paradigm for exploring the function of the secondary motor, M1, S1, and PPC simultaneously, which may be a useful tool for future use in rehabilitation response monitoring and other clinical applications of fNIR imaging.
Three right handed subjects were included in this study (one female and two male, 25 ± 3 years old). For each subject one repeat session was performed 3 to 4 weeks from the first visit to test the reproducibility of cortical activation in the same FOV. These studies were performed under the approval of The University of Texas at Arlington (UTA) Institutional Review Board protocol (IRB No. 2011-0193).
A continuous wave fNIR spectroscopy brain imager (DYNOT, NIRx Medical Technologies, LLC., Glen Head, New York) was used to map the HbO and Hb changes due to cortical activity induced by block design hand and arm motion protocols performed sequentially, as described below. In this brain imaging system, two wavelengths (760 nm and 830 nm) were simultaneously provided by two laser diodes, whose light was coupled sequentially into a maximum of 32 different fiber bundles that delivered the light to user-selected positions on the subject’s scalp. Each fiber bundle was given an equal interval of illumination time, the duration of which was determined by the number of source-detector pairs that the system had to cycle through (time-multiplexing). Though only one source was on at any one time, the intensity from each source was sine-wave modulated at a frequency in the 4 – 11 kHz range as a means of tagging its identity. Each of the fiber bundles touching the scalp’s surface was bifurcated, thus allowing each to act as both source and detector. Since the weight of the fibers, in the absence of any external support, caused discomfort to the subjects a stand was constructed to support part of the weight of each fiber bundle. With this stand in place only the length needed to adjust the position and angle of the fibers on the head surface was left to be supported by the subject’s head [Fig. 1(a) ]. The probe stand consisted of two wooden planks, a fixed metal column, and an adjustable metal column. One of the wooden planks was placed at the base to support the balance of the entire stand assembly, and the other plank was placed at the top of the metal column to support the weight of the fiber bundles. The fixed metal column was a metal mailbox stand and the adjustable metal column was used to change the height of the wooden fiber bundle support plank as some subjects were taller than others. Additionally, a wearable probe holder was made to maintain the fiber bundles in a fixed position onto the scalp during fNIR measurements. The probe holder was created in layers that were glued together. The outer layer was made of vinyl rubber, the middle was Buna rubber, and the lower layer was self-adhesive foam. The vinyl rubber was used for its rigidity and the Buna rubber for its flexibility so that their combination kept the probes pointing straight onto the scalp with the entire assembly following the head’s curvature. The foam layer that was glued onto one side of the composite rubber layers was in contact with the head’s surface so as to keep the subjects comfortable. Finally, a plastic adjustable headband, originally intended to hold a welder’s mask, was used to keep the probe holder in place onto the subject’s head. The geometrical positioning of the probes fit a 12 x 8.4 cm FOV, in which the probes were placed 2.1 cm apart laterally and 2.4 cm apart in the anterior-posterior direction [Fig. 1(b)]. This geometry enabled measurements from source-detector distances of 2.1 cm, 2.4 cm, and 3.2 cm, totaling in 178 source-detector pairs. Since the sampling rate decreased with an increasing number of source-detector combinations due to the increased time multiplexing, the fNIR spectroscopy signals could be sampled at 1.81 Hz for the chosen source-detector arrangement. Two additional measurements were also taken at the corners of the FOV, with a source-detector separation of 1.1 cm [encircled in Fig. 1(b)]. These additional measurements were used to measure the respiration and Mayer wave hemodynamics, in order to subsequently filter these out of the activation signals [31–34].
The fNIR spectroscopy probe placement was done according to the measured coronal (ear-to-ear) and sagittal (nasion-to-inion) distances. The midpoint of the right edge of the probe set was placed at the intersection of the aforementioned measurements, while parallel to the sagittal line. This way the entire assembly would cover the PMC/SMA, M1, S1 and PPC areas of one hemisphere of the brain. Since all tasks were performed with the dominant hand (right hand) the probe assembly was placed over the brain hemisphere contralateral (left hemisphere) to the hand and arm performing the movements indicated by the activation protocols. The subjects sat up straight with their head resting back in a quiet, dimly lit room. A PowerPoint animation on a computer screen guided subjects through the protocols. The data acquisition paradigm consisted of a 30 s baseline (no stimulation), immediately followed by a series of ten consecutive epochs of 15 s of stimulation and 25 s of rest, and ended with a 20 s baseline measurement, resulting in a total acquisition time of 450 s. Data was acquired for each subject for five consecutive periods lasting 450 s each, with one of five different activation protocols being executed each time. Four protocols requiring movement were finger tapping, palm squeezing, sequential finger tapping, and card flipping (supination and pronation). Finger tapping consisted of the subject repeatedly tapping with all fingers and the palm squeezing protocol consisted of the subject repeatedly forming a fist by squeezing a stress ball. For sequential tapping, the index finger, middle finger, ring finger, and pinky were correspondingly numbered from 1 to 4. The PowerPoint animation presented the numbers in a random order, as determined by a random number generator in MATLAB (MathWorks, Natick, Massachusetts), during the stimulation portion of the protocol. The last activation protocol was a card flipping task in which the subject was asked to flip individual cards from one deck to a separate deck of cards. The two decks were placed apart by a distance equal to the subject’s shoulder width and subjects took a card from the left deck and flipped it as they placed it onto the right deck. Importantly, during card flipping subjects were looking at the computer screen, not at the deck of cards, and therefore had to rely on hand sensory input to perform this task. Additionally, a protocol in which the bristle of a toothbrush was rubbed across the tips of the index, middle, and ring fingers was also performed without any concurrent hand movements to locate S1. No set stimulation frequency was required of the subjects. Thus, all motions were self-paced and were in the 1-2 Hz range for each protocol.
In addition to detecting evoked hemodynamic changes, fNIR spectroscopy is sensitive to cerebral hemodynamic fluctuations of systemic origin. Such systemic fluctuations can be caused by cardiac pulsation, respiration, and Mayer waves [15,31,35]. For typical motor activation protocols, the cortical hemodynamic response can be found in the 0.01–0.4 Hz frequency range, while physiological artifacts, such as cardiac pulsation, can be found between 0.8 and 2.0 Hz, respiration in the 0.1–0.3 Hz range, and Mayer waves at ~0.1 Hz or lower [15,35,36]. Because there is a significant overlap between the frequency spectra of respiration and Mayer waves and of the hemodynamic response due to brain activity, band-pass filtering is not effective in removing such physiological artifacts. Recent studies have used component analysis [33,35,37,38], adaptive filtering [31–34], or a combination [5,33] to remove physiological artifacts.
This study used a combination of band-pass filtering, adaptive filtering, and component analysis to filter the fNIR signals . The first step in the signal filtering was to band-pass filter fNIR signals from all source-detector separations, including two short-distance reference hemodynamics measurements [encircled in Fig. 1(b)], in the 0.01 and 0.4 Hz range as this frequency window was found to contain the power spectrum of activation signals as well as that of respiration and Mayer waves. Cardiac pulsation, which was in the 1 Hz, range was filtered out of this window though some signal aliasing is possible due to the low sampling rate in this study (1.81 Hz). Subsequently, a combination of an adaptive least mean square (LMS) filter and principal component analysis (PCA) was used to filter the band-pass filtered fNIR signals . The physiological hemodynamics reference measurement at the top left FOV was used as the adaptive filter reference for the top half of the optodes shown in Fig. 1(b) and the hemodynamics reference measurement at the bottom right FOV was used for the bottom half of this optode arrangement.
Though the fiber bundles for fNIR measurements were placed on the head with good stability, hair still played a significant role in obstructing these fiber bundles from obtaining good optical contact with the scalp. Including fNIR data from fibers with poor optical contact into the analysis resulted in significant additional noise propagating into the resulting reconstructed images. Therefore a method to remove low signal-to-noise ratio (SNR) pairs was needed. In sparse source-detector arrangements data from such detector fibers can be removed after visual inspection of their time-series signals. Since a large number of source-detector pairs was involved in this work a quantitative automated method, somewhat akin to the jackknife method , was developed to identify and remove, or ‘prune’ the noisy detector channels, as described in more detail here below.
The flowchart of the algorithm used to identify the source-detector pairs with low SNR after filtering is shown in Fig. 2(a) . Firstly, SNR values for all source-detector pairs were determined by Eq. (1):
where SNRs,d was the SNR value for a specific source-detector pair, Ps,d was the power within the 0.01 Hz to 0.4 Hz frequency range for a specific source-detector pair after being filtered, and Pd was the power within the same frequency range for a dark measurement (no light source on) with the specific detector. As SNR was dependent on detected light intensity, which in turn depended on source-detector separation, comparisons of SNR were first performed within groups of the same source-detector separation. Three groups were thus formed corresponding to 2.1 cm, 2.4 cm, and 3.2 cm source-detector separations. In each group SNR values were ordered from least to greatest and an iterative procedure was applied, the first step of which was to remove the smallest SNR value of the group creating a pruned group. Subsequently, a comparison was made between the mean and standard deviation of this pruned group and the corresponding values in the original, non-pruned group. These pair-wise comparisons between groups were performed by using the T-test and F-test, respectively. If a significant difference was not found (p > 0.05) for either a change in the mean or standard deviation of the pruned group, the lowest SNR value was discarded and the process was repeated for the detector signal with the next smallest SNR value in that group. The first SNR value that, when removed, caused a statistically significant change in either the mean or the standard deviation of the remaining group members was chosen as the noise threshold for that group. Once the SNR threshold values for each of the three source-detector separation groups were determined, the global SNR threshold was set as the minimum of the three group threshold values. Detectors at all three source-detector separations below this global threshold were excluded from further analysis. Selection of a minimum global threshold was empirically found to reject detectors with SNR values significantly lower than the mean SNR of all source-detector pairs in the FOV. The average SNR for this study was 27.33 ± 18.12 dB, while the SNR threshold was found to be 5.32 ± 0.64 dB and the number of source-detector pairs removed was 13.17 ± 7.35 out of a total 178. After pruning, the remaining data was further processed for signal visualization and image reconstruction.
Reconstruction and visualization of fNIR activation images resulting from the acquired reflectance data that passed the SNR selection criteria outlined above was performed by the open-source HomER software implemented in MATLAB . In this software activation images were reconstructed by use of the Tikhonov perturbation solution to the photon diffusion equation [41,42], which employed a regularized Moore–Penrose inversion scheme [41,42]. The reconstructed, two-dimensional images (21 x 21 pixels) represented maps of Hb and HbO changes on the cortical surface, within the detector’s FOV [Fig. 1(b)]. Furthermore, 40-s rest-activation intervals were removed from the analysis of all channels if motion artifacts were detected within a rest-activation interval. Motion artifacts were determined by spikes that were more than three standard deviations above the mean of the rest-activation intervals . If a 40-s rest-activation interval was removed, the two adjacent rest-activation intervals were concatenated together to form again a continuous time signal for each pixel. Images were then reconstructed for every 0.55 s time interval and time-averaged images of activation were created [Fig. 2(b), panel (1)] by averaging all images in the interval between 5 s and 20 s after the beginning of each activation protocol, as these times presented activation that was visually discernible over baseline values for all protocols presented in this work.
After the fNIR images were reconstructed, regions of activation were determined by a T-test in an approach similar to that used in fMRI image analysis . Areas of activation were identified by comparing, for the time-series of each pixel, the mean and standard deviation of the changes in HbO and Hb during the activation period (5 s to 20 s) relative to changes in baseline values acquired just prior to activation (−10 s to 0 s). After considering Bonferroni’s correction for multiple comparisons , pixels in reconstructed images were considered to have significant activation when p < 0.0001, and color-coded T-value maps were plotted [Fig. 2(b), panel (2)]. The resulting T-value maps showed that the activation areas expanded further than what was seen from the time-averaged activation images.
Once activation regions were identified it was subsequently investigated if there were cases where spatially distinct sub-regions had similar time-series activation profiles for the different hand and arm motions performed in each activation protocol. To that end, the time-series pattern of fNIR signals with significant activation (p < 0.0001) were decomposed into a signal harmonic subspace and then clustered into source-detector pairs with similar harmonic content using a cluster component analysis (CCA) algorithm . The details of how this algorithm was applied to this study are described in the Appendix. CCA grouped the source-detector pairs with similar temporal patterns into clusters that were displayed as connected circles having a different color for each cluster [Fig. 2(b), panel (3)]. In that panel the grey circles indicate the location of bifurcated source-detector fibers with no significant activation as determined by the above mentioned T-test.
It should be noted that due to the diffuse light scattering occurring for fNIR spectroscopy signals a detector channel can see activation changes from multiple spatial locations, which can result in some fiber locations being assigned to more than one cluster. Nevertheless, it was found that CCA showed little overlap between different source-detector pair clusters. For example, for a palm squeezing task in Subject 2 there was overlap for only a single source/detector fiber location [Fig. 2(b), panel (3), right edge]. An example of the difference in time-series profiles between different source-detector clusters is shown in Fig. (2), panel (4). Finally, it is important to mention that, in conjunction with CCA, it is critical to use the fNIR signal filtering procedure described in Section 2.4 above to avoid clustering errors due to activation signal contamination from large amplitude physiological hemodynamics.
After creating the activation images and determining the regions of activation, spatial and temporal metrics were used to quantify the repeatability of the location and temporal pattern of the hemodynamic response for each protocol. These metrics were applied to those pixels in the activation images which were considered significantly active by the T-test (p < 0.0001) where the spatial metric, the center of mass of the activation area, was applied to the time-averaged images and the temporal metrics, time-to-peak (TtP) and duration of activation, were applied to the time-series images. The center of mass for the activation area was determined separately for the x- and y-directions, and was defined by Eq. (2), where CoM was the center of mass, wx,y was the change in HbO for a pixel, and d is the distance of a pixel from the left edge of the FOV for the x-direction or the distance from the FOV bottom for the y-direction.
Additionally, in analyzing time-series images the TtP was defined as the difference between the beginning of the activation interval and the time point of maximum HbO change . To determine activation duration the time-series of each pixel was first placed into a matrix where the rows were the pixels and the columns were pixel values at each 0.55 s time point. Then a k-means clustering algorithm  was applied to divide the time-series pixels into activation, baseline, and deactivation. The duration was then defined as the total amount of time a pixel was found to have activation .
All statistical tests were performed, using SAS 9.1 (SAS Institute Inc., Cary, North Carolina) to see if there was a significant difference (p < 0.05) between visits or between groups. One-sampled T-tests were performed to see if there was a difference in spatial and temporal metric means between visits. The null hypothesis was defined as a zero difference in the means between the two visits. After source-detector pairs were clustered with CCA, the temporal metrics, TtP and duration, were compared between clusters for sequential finger tapping and card flipping using a two-sampled T-test. Furthermore, the frequency that a given cluster was assigned to each source-detector pair was tested across motion tasks using the Fisher’s exact test. The analysis was first applied for all tasks as a group, and then applied as a post-hoc analysis to compare each possible pair of tasks.
For each of the five activation protocols performed by each subject the detector signal intensities for the properly executed (no motion artifacts) 40-s stimulation-rest intervals of each source-detector pair retained after data pruning were time-averaged for each 0.55 s time point (a total of 72 time points). Activation images over the measured FOV were then reconstructed for every time point from the averaged temporal response of all kept source-detector pairs . The resulting images were used to determine spatial locations and temporal patterns of the hemodynamic response due to each protocol, and their repeatability between visits. Both HbO and Hb data were analyzed. However, it was found that HbO and Hb activation patterns were qualitatively similar, except that Hb changes were negative and delayed by 1-2 s, consistent with prior literature . Therefore, for brevity, the results presented in this work are primarily for HbO, though some Hb results are also presented to demonstrate this similarity.
In Fig. 3 , time-averaged images for the 5 s to 20 s interval are shown for each task, for both visits of each subject. From these images the locations of activation were visually identified for each protocol, and were seen to be repeatable for each subject where the center of mass of the activation area within a cortical region was not significantly different between visits in either the x-direction (0.47 ± 0.88 cm, p = 0.39) or the y-direction (0.64 ± 0.94 cm, p = 0.31). The probe placement accuracy between repeat sessions was ~3 mm in the x- and y-directions, as estimated from sagittal and coronal measurements done for the probe placement procedures described in the Methods section. No concurrent MRI measurements were performed to mark fiber bundle placement with respect to the cortical surface for each subject, thus identification of the PMC/SMA, M1, S1, and PPC regions was based on inference from previously reported fMRI and PET studies involving similar types of hand and arm motions [16,22,24,46–56]. For example, the placement of the dividing line separating S1 from M1 in Fig. 3 was made under the assumption that these cortical regions were each ~2 cm in width in the anterior-posterior direction, as MRI images have shown for the adult human cortex [15,24,57]. Also, as was found in previous fMRI studies, finger tapping elicits activation in the M1 and S1/PPC cortical regions [16,46–48], and sensory stimulation by the bristle end of a toothbrush results in activation of the S1/PPC regions [24,46,49–51].
Activation in the S1 and M1 regions for finger tapping and in the S1/PPC and M1 regions for sensory stimulation was also seen with fNIR imaging (Fig. 3, rows 1 and 2), but the activation of the PPC during finger tapping was not seen. Nevertheless, T-value maps of the fNIR data did show that there was significant activation (p < 0.0001) in the PPC (Fig. 5 below, column 1). The palm squeezing, sequential tapping and card flipping protocols increased the amount of effort compared to finger tapping either physically or mentally. The squeezing of a stress ball (palm squeeze) increased the physical effort and induced activation in the M1, S1, and PPC regions of the cortex (Fig. 3, row 3), in agreement with previous fMRI studies . Next a sequential tapping protocol increased the mental effort required of the subjects in executing a physical task. Images for this protocol (Fig. 3, row 4) show additional activation in the PMC/SMA cortex. Finally, the same cortical regions were found to be activated due to card flipping (Fig. 3, row 5) as for sequential finger tapping. It is important to iterate that sensory feedback was important for the card flipping protocol due to the subject’s need to feel for the card’s location by hand without looking at the deck, since the subject’s visual attention was directed towards the computer screen. The results for both the sequential tapping and card flipping protocols are consistent with the notion that when sensory feedback is needed in order to execute a motion, there is activity in the PMC/SMA cortical areas [22,53–56].
In addition to identifying the spatial location of activation areas, fNIR imaging also enables investigation of the temporal dynamics of each cortical region due to the different activation protocols. Figure 4 shows the evolution of activation patterns for each protocol, in time-averaged 5 s blocks for the first 30 s for Subject 1. Similar patterns were found for the other two subjects (time-series images not shown) and were repeatable between the two imaging sessions for all subjects. In order to make some quantitative comparisons between temporal activation patterns, the TtP and duration of activation elicited by each protocol in each spatially distinct activation region was quantified as described in the Methods Section. The TtP (1.17 ± 2.38 s) and duration (1.36 ± 2.97 s) difference between visits for each protocol of each subject were not found to be significant (p > 0.05). The corresponding results for TtP are shown in Table 1 and for duration in Table 2 .
From observing Table 1 it is seen that cortical regions involved with higher amount of physical or mental effort for any given activation protocol have shorter TtP values, i.e. they reach an activation peak more quickly. For example, the M1 column of Table 1 shows that finger tapping, the simplest motor protocol, had a TtP of 18.13 ± 2.35 s, whereas palm squeezing (14.21 ± 1.84 s) and card flipping (9.55 ± 2.47 s) had earlier TtPs since they used more arm and hand muscles. Though there is a strong trend, the small number of subjects (n = 3) prevented any statistical analysis to conclude a significant change due to the low statistical power (β = 0.32). However, by grouping those source-detector pairs that CCA determined to have similar temporal patterns (Section 2.6) during the sequential tapping and card flipping tasks, significant differences were found in the TtP between motion tasks, as shown in Section 3.2.2 below. Furthermore, the more refined movements of sequential tapping also had an early TtP (12.39 ± 3.55 s) matching with previous fNIR results . In addition, the S1/PPC column of Table 1 shows a much shorter TtP value for card flipping due to the need for sensory information from the hand to find the decks of cards. Comparing TtP values across the rows of Table 1, it is seen that sequential tapping peaks relatively early in M1, but card flipping peaks early both in M1 and S1/PPC. These results are consistent with the idea that both protocols require more effort compared to the other protocols used in this work, but card flipping is more demanding in terms of sensory feedback. Finally, as a point of validation, the constant sensory stimulation from a brush resulted in a hemodynamic response with a strong early rise and peaked at 14.21 ± 1.84 s, which agreed with previous fMRI studies using vibrotactile stimulation of the finger tips [46,58].
Further investigation shows that increased effort also increases the duration of the hemodynamic response, as shown in Table 2. Palm squeezing, the most physically demanding task, had a very large duration in comparison to the other tasks, as seen in the columns of M1 (35.13 ± 3.21 s) and S1/PPC (33.07 ± 3.01 s). The palm squeezing hemodynamic response did decay before the measured 40 s (data not shown), though the figure only shows up to 30 s due to space constraints. Moreover across the rows, sequential tapping and card flipping both had shorter durations in M1 (15.43 ± 2.46 s and 16.23 ± 2.87 s, respectively) than those found in PMC/SMA (22.67 ± 3.57 s and 23.18 ± 2.84 s, respectively) and S1/PPC (22.13 ± 3.24 s and 23.17 ± 2.67 s, respectively). For the sequential tapping and card flipping tasks, by using CCA to group source-detector pairs with similar temporal patterns (Section 2.6) into a single cluster, the difference in activation duration of the cluster containing the PMC/SMA and S1/PPC regions versus duration of the M1 cluster can be shown to be statistically significant (Section 3.2.2). The large amount of continuous external sensory information and feedback needed for these latter tasks explains the extended activation duration for the PMC/SMA and S1/PPC regions. Though the increased effort shortened the TtP in M1 (Table 1), the amount of physical strain was not to the extent of that in palm squeezing and this may be the reason why the duration of hemodynamic response for sequential tapping and card flipping in M1 was not as prolonged.
Though TtP and duration were found to be useful metrics for identifying similarities and differences between time-series activation patterns, these were not necessarily indicative of the overall temporal pattern similarity between different cortical regions. It was hypothesized that if similar time-series activation patterns existed between spatially distinct cortical regions, irrespective of activation amplitude, this pattern similarity could indicate a possible functional relation between these regions.
In this analysis the source-detector pairs with statistically significant activation over the background hemodynamic fluctuations were first identified using a T-test , and subsequently fNIR signals were clustered into regions with respect to their similarity in temporal patterns by use of the CCA algorithm , as described in the Methods Section. The results for finger tapping, sensory stimulation, and palm squeeze T-value maps and clustered regions are shown in Fig. 5 for Subject 1. Similar results were found for the other subjects, but are not shown here for brevity. The HbO T-tests indicated activation region sizes that appear larger than the activation images shown in Figs. 3 and and4.4. Moreover, in the cluster images (Fig. 5) S1 and PPC were grouped into one cluster and the M1 region into another for finger tapping and palm squeezing, suggesting communication between the cortical regions, which is consistent with previous fMRI and electrophysiological studies [15,16,46,48,59]. Additionally, the Hb T-maps presented smaller areas of activation in the same locations as that found for HbO, which is consistent with results of previous studies [15,60].
The cluster analysis was then focused on the sequential tapping and card flipping protocols since these protocols stimulated multiple cortical regions within the FOV due to the complexity of the tasks involved. Comparison, of the cluster map results for these two tasks showed differences in the temporal relation between the involved cortical regions (Fig. 6 ). After CCA, both temporal and spatial group analysis were performed to prove the statistical significance of these observed differences. The temporal group analyses proved the statistical significance of the trends seen in Tables 1 and and2,2, above. Specifically, analysis on the card flipping data was performed by combining CCA-selected source-detector pairs with similar patterns in the M1 region to compute a cluster-average TtP (Table 1, Section 3.1) and this value was compared with that of corresponding clusters created from CCA of the finger tapping and palm squeezing tasks. The TtP for card flipping was found to be significantly different (p = 0.019) with respect to these other two motion tasks with high statistical power (β = 0.98). Similarly, if the clusters with similar patterns for sequential finger tapping were merged, the duration (Table 2, Section 3.1) of activation in M1 was found to be significantly different from the S1/PPC and PMC/SMA cluster duration (p = 0.04). The spatial group analyses showed that source-detector patterns with activation were arranged in different spatial patterns depending on the motion task performed. Specifically, the difference in the frequency that a given cluster was assigned to each source-detector pair was tested across motion tasks with Fisher’s exact test and was found to be significantly different between all tasks (p < 0.0001). Also, the post-hoc analysis found that each task was significantly different from all the others (p < 0.05), except between palm squeezing and card flipping (p = 0.56). Moreover, group T-maps created by the T-test of the concatenated time-series pixels of each subject for the sequential tapping and card flipping tasks showed that the S1/PPC, M1, and PMC/SMA regions had significant activation consistently for all subjects (p < 0.0001) .
When comparing T-maps with corresponding HbO activation images (e.g. Figure 4 with Fig. 5, or Fig. 3 with Fig. 6) it was found that, even after applying Bonferroni’s correction (p < 0.0001), T-maps typically indicated larger areas of activation that engulfed the areas seen in HbO images. This apparent activation pattern difference occurred because T-maps included areas of low amplitude activation that were statistically significant, as deduced from univariate analysis of each pixel’s amplitude fluctuations with respect to baseline values. In contrast, in HbO images the color assigned by the HomER analysis software for a pixel with a given activation amplitude depended on its relative magnitude with respect to the image maximum and minimum values. In that case the low amplitude activation pixels were assigned colors very close to those of background and were therefore not visually discernible, thus making the apparent size of activation areas smaller. Furthermore, due to the better ability of T-maps to indicate lower amplitude activation areas, there were few instances where activation in functionally distinct regions was only discernible in the T-maps. An example can be seen for Subject 1 doing finger tapping, where no activation was discernible in the PPC (FOV bottom) for HbO images (Fig. 4 – top row) but activation was seen in the PPC for T-maps (Fig. 5 – top left panel), consistent with prior fMRI studies [16,46–48]. This difference was also observed in the other subjects (images not shown for brevity).
Fascinatingly, in Fig. 7 it is seen that the secondary motor cortices (PMC and SMA) clustered together with the S1 and PPC areas for sequential tapping, which depended on visual cues to execute the protocol. In contrast, for the card flipping task, which depended more heavily on sensory information from the hand and arm to identify the location of the decks of cards, the S1, PPC, and M1 regions clustered together. Though previous fMRI studies have used sensory based protocols to determine the relation of the PPC to the motor and sensory systems [46,48,60], there are no studies to our knowledge, which compare the temporal hemodynamic profiles of these separate cortical regions. On the other hand electrophysiological studies in monkeys have shown that the PMC/SMA region is active before the movement, between movements, and during the movement if there is an external cue [49,53,54,56,62]. The fact that the spatiotemporal clustering results were consistent with prior literature suggests that similarity in hemodynamic activation patterns, though slower and secondary to neuronal activation, could be indicative of the functional relations between these different cortical regions.
Resting state functional connectivity has also been used to determine the functional relations between cortical regions by looking at their temporal correlations [63,64]. Since this is done at a resting state, a cortical region may be found to be well correlated with several regions of the brain, though some of these connections may not be functionally relevant for a specified task. Previous fNIR studies, including this study (data not shown), found broad cortical connectivity between motor and sensory cortical regions using resting state and functional connectivity analysis. On the other hand, the use of CCA for a stimulation paradigm gives functional relations between cortical regions during a specified task, thus giving a more indicative task-specific functional connection between cortical regions.
This work demonstrated the feasibility of performing fNIR measurements over an extended FOV (12 x 8.4 cm) spanning the primary sensorimotor and secondary motor cortex and identified spatiotemporal interrelations between different activation regions for various hand and arm movements. In addition to activation by the popular finger tapping and sensory stimulation protocols [46–48], palm squeezing, sequential tapping, and card flipping (supination and pronation) protocols were also demonstrated for this extended FOV. The variety of protocols enabled detecting concurrent activation from several cortical regions within the measured FOV including the SMA, PMC, M1, S1 and PPC. Our results of the areas activated and temporal activation characteristics of the hemodynamic response for each protocol matched well with previous fMRI and electrophysiological studies [49,56,62], though to our knowledge, card flipping has not been done before in previous fMRI or fNIR studies. Additionally, a CCA algorithm was used to cluster cortical regions with similar temporal activation patterns , indicating functional relations between cortical regions for a specified task. Analysis of time-series fNIR data for the different protocols performed showed earlier peaking of activation as well as longer activation duration for the tasks requiring higher physical effort (palm squeezing), or higher mental effort (sequential tapping and card flipping). As can be seen from Table 1, with increasing physical demand or sensory feedback needed to perform a task, there was an earlier TtP found in M1, S1 and PPC with finger tapping being the simplest motor task. On the other hand, sequential tapping did not follow this pattern in S1 and PPC since there was no increase of sensory stimuli from the arm or hand, but there was an increase of sensory information from the visual cortex in order to identify the number with the proper finger. Additional analysis in Table 2 presented an increase in duration in M1, S1, and PPC if there was large increase in physical demand (palm squeezing), or in the PMC/SMA, S1, and PPC if there was a large increase in sensory feedback (sequential tapping and card flipping).
A clustering analysis method for comparing the similarity of time-series activation data irrespective of amplitude was also used to compare patterns between spatially distinct activation regions . It was found that when the subject relied more heavily on visual guidance to perform a protocol, temporal patterns between the PPC/S1 and the secondary motor regions (PMC/SMA) of the cortex had similar temporal patterns. Furthermore, when sensory information from the hand directed the movement of flipping playing cards, the PPC/S1 had similar temporal profiles with the M1 region of the cortex. With the large FOV and the diverse protocols used in this study, there is noticeable importance in imaging the PPC and secondary motor cortical regions (PMC/SMA) due to the sensory feedback they provide for the execution of more complex tasks.
The findings of this work indicate that there is rich information to be obtained from the spatiotemporal relations between different cortical activation regions when performing fNIR imaging over an extended area encompassing the sensorimotor, secondary motor cortex, and PPC. Unfortunately, the limited number of source-detector channels in our currently available fNIR imaging system limited functional imaging to only one brain hemisphere for this FOV size. The commercial availability of future fNIR imaging systems with substantially higher source-detector channels will enable performing extended FOV measurements bilaterally in the brain. Such measurements will be of great use for many functional imaging studies, including ones underlying the motivation for the present work, namely the longitudinal monitoring of cortical plasticity during rehabilitation of patients with motor deficits [5,7,9].
The authors would like to thank Dr. Hanli Liu for making the fNIR DYNOT imager available for this work, Dr. Fenghua Tian for his assistance in data collection, and the subjects who volunteered their time to participate in this study.
This Section describes how the CCA algorithm, originally used in the fMRI field , was adapted in this work to perform the analysis of fNIR imaging data. In the fNIR spectroscopy field, blood concentration and oxygenation changes secondary to neuronal activation are reconstructed into maps of absorbance change after applying a modified Beer-Lambert law . More specifically, maps for the change in optical density (ODλ) for every wavelength (λ) detected by a given number of source-detector pairs (SD) at any one time point are reconstructed by inverting Eq. (A1). In that equation yλ is an SD x 1 array of the ODλ for each source-detector, Aλ is the joint probability distribution for each source-detector pairs , and xλ is the reconstructed ODλ image (N pixels) for one time point. A ODλ image is then reconstructed by taking the regularized inverse of Aλ as shown in Eq. (A2), where AT is the transpose of A, α is the regularization parameter (0.01), smax is the largest eigenvalue of E[AAT], and I is the identity matrix :
However, because the target of this work was to identify and cluster fNIR signals with similar temporal patterns it was necessary to consider data for all time points at once. Therefore time-series source-detector signals and their corresponding time-series images were expressed in a form described by Eq. (A3), where Xλ (time x N) was the time-series for each pixel from the reconstructed ODλ images and Yλ (SD x time) were the ODλ time signals for each source-detector pair:
Before the fNIR signals could be decomposed into a signal harmonic subspace, a model matrix was created in which each column consisted of a full time-series (450 s) of sine waves whose frequencies were multiples of the fundamental frequency. A fundamental frequency of 0.02 Hz was assigned as that was the frequency equal to the inverse of the duration of a task epoch. Each column of the model matrix did not extend the full 450-s time-series if some of the 40-s activation-rest intervals were removed due to motion artifacts (see Section 2.6). All harmonics (f) of the fundamental frequency up to half the sampling rate (Nyquist frequency) were input into the model matrix H (time x f). From Eq. (A3), each ODλ was then decomposed into harmonic components using a general linear model (YλT = Hλθλ + v), where θλ (f x SD) was the harmonic image matrix consisting of the weight coefficients (arbitrary units of ODλ) to each frequency harmonic and v (time x SD) was the noise estimated form the residuals of the linear model. This noise consisted of instrumentation dark noise and the residual signal from the physiological hemodynamics that were not removed by the PCA and adaptive filtering procedures applied to the measured fNIR signals. This noise was determined to be effectively white. The above considerations resulted in Eq. (A4):
Beginning with the reconstructed ODλ images for the whole time series, a least squares estimate of θ () was computed [Eq.(A5)], which consisted of the weight coefficients to each frequency harmonic. The resulting estimation error was then found by Eq. (A6):
By computing the covariance matrix of the least squares estimate [Eq. (A7)] and the estimation error of the harmonic image matrix [Eq. (A8)], the covariance matrix of the noiseless harmonic image matrix θ was found [Eq. (A9)]. By computing the eigenvalues of the covariance matrix of θ [Eq. (A9)], using Schur decomposition , and eliminating all components with negative amplitude, the harmonic subspace was reduced to include only those harmonics that spanned the ODλ time-series signals. The reduced harmonic subspace for those ODλ signals was determined by Eq. (A10) where was the reduced eigenvector matrix that included only the eigenvectors corresponding to positive eigenvalues Λλ [from Eq. (A9)]:
Since the harmonic subspace was in arbitrary units of ODλ, the conversion to change in HbO and Hb, Eqs. (A11)-(A12), for each matrix entry was performed by applying the modified Beer-Lambert Law at the two wavelengths . In these equations εHbOλ was the extinction coefficient of HbO at wavelength λ, εHbλ was the extinction coefficient of Hb at wavelength λ, and L was the pathlength of the light through the tissue (calculated as source-detector separation times a differential pathlength factor of six ):
Finally, only those source-detector pairs with significant HbO and Hb activation were fed into the CCA algorithm. Significant activation was identified by comparing, for the time-series of each fNIR signal, the mean and standard deviation of the changes in Hb and HbO during the activation period (5 s to 20 s) relative to changes in baseline values acquired just prior to activation (−10 s to 0 s). After applying Bonferroni’s correction for multiple comparisons  source-detector pairs were considered to have significant activation when p < 0.0001. The CCA algorithm  then clustered fNIR signals into a self-determined number of groups with similarly weighted harmonic components and therefore similar time-series shapes irrespective of activation magnitude.