2.1.

Subjects

Five pediatric controls, Subjects 1–5 (two female and three male, $8.4\pm 2.3\text{years}$ old), and five children with hemiparetic CP, Subjects 6–10 (three female and two male, $8.4\pm 2.3\text{years}$ old), were included in this study. All the controls were right handed. Only children with CP that had a successful anatomical magnetic resonance imaging scan performed, indicating a single subcortical lesion in one of the two brain hemispheres, were included in this study. Three of the children with CP, Subjects 6–8, had a subcortical lesion in the left hemisphere (LH), which affected movement of their right hand. Subjects 9 and 10 had a subcortical lesion in the right hemisphere (RH), which affected movement of their left hand. Informed consent was obtained from all subjects and their legal guardians. Except for Subjects 4 and 6, it was possible to attain a second visit for fNIRS measurements within two months of the first. These studies were performed under the approval of The University of Texas Southwestern Medical Center at Dallas (UTSW) Institutional Review Board protocol (IRB No. :042007-064).

2.2.

Measurements

A continuous-wave (CW) fNIRS brain imager (CW-5, Techen Inc., Milford, Massachusetts) was used to map the HbO and Hb changes induced by motor cortex stimulation. The source-detector geometry is shown in Fig. 1 . Eight detectors were placed over each brain hemisphere to cover a relatively large area of the motor cortex. The geometrical arrangement of sources and detectors was attached onto the subjects’ heads by perforated Velcro straps. Eight laser sources emitted at $690\mathrm{nm}$ and eight at $830\mathrm{nm}$ , such that each optical fiber bundle delivered light of both wavelengths at each source location [circles in Fig. 1], simultaneously. Each source bundle had up to four detectors [squares in Fig. 1] within a $3\text{-}\mathrm{cm}$ distance, and each detector received signals from up to two source bundles. As a result, there were 28 possible source-detector channel combinations for each wavelength, as indicated by the straight lines connecting sources to detectors in Fig. 1. Activation in cortical areas within the probes’ field of view $(6\times 20\mathrm{cm})$ was monitored simultaneously by all source-detector pairs, because the CW-5 enables all laser sources to be on at the same time with distinct modulation frequencies ( $6.4–12.6\mathrm{kHz}$ , with an increment of $200\mathrm{Hz}$ ). The backreflected light detected was sampled at a rate of $100.16\mathrm{Hz}$ , which was later demodulated and downsampled to $20.03\mathrm{Hz}$ to reduce data-set sizes.

Fig. 1

(a) Placement geometry of the fNIRS probes over the motor cortex. The circles represent sources, and the squares represent detectors. (b) Overall instrumentation setup.

In addition to measuring the cortical hemodynamic response by optical means, the finger-tapping, respiration, and cardiac pulsation patterns were also measured simultaneously for each subject. The finger-tapping rate was measured using a custom-built, capacitance-based tapping board that produced a bilevel voltage output for each tap. The respiration and cardiac pulsation measurements were amplified with a Brownlee Model 410 amplifier (AutoMate Scientific, Berkeley, California), and recorded by the CW-5 system. The respiration and cardiac measurements were amplified by a factor of 200 and $30\mathrm{mV}∕\mathrm{mV}$ , respectively. Breathing patterns were measured using a respiration belt with a piezoelectric transducer (Sleepmate Technologies, Glen Burnie, Maryland). The belt was wrapped around each subject’s lower chest while keeping the transducer away from the heart in order to prevent cardiac pulsation noise from leaking into the respiration measurement. Cardiac pulsation was measured using a pulse oximeter (Nellcor Inc., Boulder, Colorado) attached to the index finger on the hand not in use for tapping. In addition, movement of the forearm flexor muscles controlling the fingers was measured by electromyography [(EMG), B & L Engineering, Santa Ana, California]. The EMG measurements were not used in our filtering procedures but instead served as quality control sentinels for our measurements, making sure that the specified hand was the only hand tapping, while the other was kept still. A diagram of the overall instrumentation setup is shown in Fig. 1, in which EMG (R) and EMG (L) stand for the EMG taken from the right forearm and left forearm, respectively. All nonoptical physiological measurements were fed into the CW-5 and digitized by an analog-to-digital converter, sampled at $100.16\mathrm{Hz}$ and later downsampled to $20.03\mathrm{Hz}$ , to reduce data-set sizes, as was done for the reflectance signals. Therefore, all signals had a common time base with the fNIRS measurements, which enabled the use of the reference physiological measurements in the subsequent filtering procedures described in Section 2.5.

2.3.

Protocol

The fNIRS probe placement was done according to the measured coronal (ear-to-ear) and sagittal (forehead-to-back-of-head) distances. The center of the probe set was placed at half the distance of the aforementioned measurements, which was considered the estimated midpoint of the motor cortex. A PowerPoint animation on a computer screen guided subjects through the finger-tapping protocol. No set tapping frequency was required of the pediatric subjects. Thus, the tapping frequency was self-paced, which resulted in average finger-tapping rates in the $1–2\mathrm{Hz}$ range within the $15\text{-}\mathrm{s}$ tapping period. Tapping consisted of moving all fingers of one hand, while the wrist was held down by soft straps to keep it from moving. The subjects touched the tapping board while tapping and rested their hand on the board during the no-tapping intervals. The data acquisition protocol consisted of a $30\text{-}\mathrm{s}$ baseline (no tapping), immediately followed by a series of ten consecutive epochs of $15\mathrm{s}$ of tapping and $25\mathrm{s}$ of rest, and ended with a $20\text{-}\mathrm{s}$ baseline measurement, resulting in a total acquisition time of $450\mathrm{s}$ . Measurements were separately acquired for both left and right finger tapping. The subjects sat up straight with their head resting back in a quiet, dimly lit room. All subjects were video recorded during the measurements in order to verify, along with the aforementioned EMG measurements, that they performed the finger-tapping tasks properly.

2.4.

Visualization of fNIRS Signals

Processing and visualization of the time-series reflectance data acquired by the CW-5 system is usually performed by the open-source HOMER software implemented in MATLAB (Mathworks, Natick, Massachusetts).²⁰ HOMER filters the signal using low-pass and high-pass Butterworth filters to reduce noise from frequencies outside the known physiological range. Furthermore, the software applies principle component analysis (PCA) to reduce any global hemodynamic signals, detrends signals to reduce any long-term drifts, and averages the detected signal for individual source-detector pairs over user-selected time intervals of tapping and rest. Subsequently, using the reflectance measurements from all detector channels activation images are reconstructed by use of the Tikhonov perturbation solution to the photon-diffusion equation,²¹ which employs a regularized Moore–Penrose inversion scheme.²² The reconstructed, two-dimensional images $(21\times 20\text{pixel})$ represent maps of Hb and HbO changes on the brain cortex surface, within the detector’s field of view [rectangle in Fig. 1] as a result of finger tapping. As Hb dynamics are highly correlated with those of HbO²³ and have significantly lower amplitudes than the latter, making them susceptible to cross-talk from HbO²⁴ and to interference from physiological artifacts, we have focused on the analysis of HbO dynamics only.

2.5.

Filtering

In addition to detecting evoked hemodynamic changes, fNIRS is sensitive to cerebral hemodynamic fluctuations of systemic origin. Such systemic fluctuations can be caused by cardiac pulsation, respiration, and Mayer waves. For typical motor cortex activation protocols, the cortical hemodynamic response can be found in the $0.1–0.4\mathrm{Hz}$ frequency range, while physiological artifacts, such as cardiac pulsation, can be found between $0.8–2.0\mathrm{Hz}$ , respiration in the $0.1–0.3\mathrm{Hz}$ range, and Mayer waves at $\sim 0.1\mathrm{Hz}$ or lower.^{25, 26} Because there is a significant overlap between the frequency spectra of respiration and Mayer waves and of the hemodynamic response due to brain activity, and there is an overlap in frequency spectra of cardiac pulsation and the tapping frequency, bandpass filtering is not effective in removing such physiological artifacts. Previous attempts in the removal of these physiological artifacts from fNIRS signals have used component analysis^{27, 28, 29, 30} or adaptive filtering.^{30, 31, 32, 33, 34}

Two types of component analysis have previously been used to filter fNIRS signals, PCA,²⁷ and independent component analysis.^{28, 29, 30} Component analysis is advantageous because no additional measurements are used to track the physiological artifacts. The disadvantage of these algorithms is that they need training data to properly identify the physiological artifacts and that the physiological artifacts may contribute to the signal differently once the study begins, thus making it difficult to attribute each component to a specific artifact or to define the number of components needed to be removed.²⁷

To better track the physiological artifacts, additional physiological measurements can be simultaneously taken or the physiological artifacts and the hemodynamic response can be modeled, as in adaptive filtering. The least mean square (LMS) filter is a finite impulse response (FIR) filter in which the filter coefficients adjust in order to minimize the least-squared error in fitting the measured noise measurements to the noisy signal.^{30, 31, 32, 33} The extended Kalman filter, otherwise, can filter the fNIRS signals by statistically modeling the physiological artifacts and the hemodynamic response due to brain activity.³⁴ Unfortunately for the extended Kalman filter, a model of the hemodynamic response for children with CP is unknown and may vary, depending on the type of lesion the child has.

In order to overcome the disadvantages of individually filtering with component analysis or adaptive filtering, a combination of PCA and LMS filtering was used in this work. Using the additional reference physiological measurements, physiological artifacts that were within the same frequency range as the tapping frequency or the cortical hemodynamic response were filtered. More specifically, the fNIRS signals were first bandpass filtered between 0.01 and $2\mathrm{Hz}$ , using a Butterworth filter, to exclude frequencies of nonphysiological origin from the signal. Subsequently, baseline correction compensated for the baseline drift that could occur throughout the repeated cycles of activation and then PCA was applied to the fNIRS signals. In contrast to the baseline procedure implemented in HOMER, which uses only the first time point of the experiment as baseline reference, we corrected the baseline for each activation interval individually by taking the average of the $5\mathrm{s}$ of data before the beginning of each activation interval and then subtracting that value from the activation data in that interval. In the implementation of PCA, we empirically found that removing the two largest eigenvalues, which represented components in the frequency range of Mayer waves and cardiac pulsation physiological artifacts, resulted in the best compromise between the reducing global background signal while not substantially reducing the contrast of activation areas in the reconstructed images as was done by Zhang ²⁷ In our analysis, we saw that the percent of the first two eigenvalues to the sum of all eigenvalues represented most of the baseline signal $(93.04\mathrm{\%}\pm 3.88\mathrm{\%})$ . In addition, in some cases the second eigenvector contributed to larger variations in the baseline signal than the first eigenvector. It was thus decided that the removal of the first two eigenvectors was optimal for removing baseline noise without compromising signal quality. Finally, the signal was adaptively filtered for respiration first, and then cardiac pulsation,³⁵ using an LMS algorithm.^{31, 32, 33} The reference noise signals [Fig. 1] that were used in this adaptive filtering procedure were first low-pass filtered at $5\mathrm{Hz}$ to reduce the contribution of electronic noise. No difference was found in the order in which the adaptive filters were applied, but using PCA after adaptive filtering was found to greatly decrease the signal quality.

2.6.

Averaging

After filtering the signals as described in Section 2.5 and excluding the tapping-rest intervals that were not executed properly, as judged by video recording and EMG, data from all source-detector pairs were processed by HOMER. Because the sampling frequency was $20.03\mathrm{Hz}$ , each time point represents every $0.05\mathrm{s}$ . The properly executed $40\text{-}\mathrm{s}$ tapping-rest intervals were averaged at each time point (total of 800 time points) giving a $40\text{-}\mathrm{s}$ averaged temporal response for each source-detector pair. Activation images of the motor cortex were then reconstructed for every time point from the averaged temporal response of all source-detector pairs.¹⁷ Time-averaged HbO activation images were formed by averaging HbO activation images within a window that started $5\mathrm{s}$ after the tapping began to $5\mathrm{s}$ after the tapping had stopped ( $5–20\mathrm{s}$ relative to the beginning of each $15\text{-}\mathrm{s}$ tapping interval) because that was the time window during which most of the activation signal was observed.

2.7.

Identifying Areas of Activation

Activation areas were identified by a $k$ -means clustering algorithm³⁶ implemented in MATLAB. The algorithm was applied to identify areas of activation in two different images: the time-averaged HbO images and the two-dimensional spatiotemporal plots representing a concatenation of the $40\text{-}\mathrm{s}$ averaged temporal responses of each pixel.

The algorithm employed three cluster means that were initialized at the maximum pixel value, minimum pixel value, and zero. The algorithm then placed each image pixel into one of the clusters, as judged by the closest Euclidean distance [Eq. 1] of the value of each pixel (of coordinates $x,y$ ) to the mean of the three clusters. However, if the minimum distance was equal between two clusters, then the pixel was placed into the cluster with the smaller mean

The algorithm recalculated the mean of each cluster after processing the whole image and repeated the process until the means of each cluster no longer changed, in which a maximum of 11 iterations were necessary for the time-averaged HbO images. The cluster that initially had a mean equal to the maximum pixel value in the image represented pixel values related to activation [Fig. 2 ]. The cluster that initially had a mean equal to zero represented pixel values related to the baseline. Lastly, the cluster that initially had a mean equal to the minimum pixel value in the image represented deactivation.

Fig. 2

Pictorial descriptions of the spatial and temporal image metrics: (a) Measurement of distance from center for both hemispheres. (b) Duration was computed as the average temporal length of activation (white areas) per pixel. The pixels were concatenated from the LH to the RH. (c) The area difference metric measured as the difference in activation area between the middle area (MA) and the contralateral hemisphere to the tapping hand (i.e., RH for left finger tapping or LH for right finger tapping). (d) Time-to-peak measurement in which the arrow represents the time that the signal took to reach its peak from the beginning of the tapping interval.

Similarly, the temporal responses of each pixel were clustered to identify activation, baseline, and deactivation. The $40\text{-}\mathrm{s}$ averaged temporal responses of each pixel were concatenated from the LH to RH to give a two-dimensional plot, in which the $y$ -axis represented the pixels and the $x$ -axis represented time. Then the previously discussed $k$ -means clustering algorithm³⁶ was applied to the two-dimensional plot, in which a maximum of 20 iterations were needed to identify the three clusters. The pixels with activation in the two-dimensional plot [Fig. 2] showed the time and location of areas of activation during finger tapping. Subsequently, a number of spatial and temporal metrics were derived from the resulting time-averaged and two-dimensional activation images, as described here below.

2.8.

Metrics Derived from the Reconstructed Images

We identified spatial metrics (distance from center and area difference) from the time-averaged HbO images and temporal metrics (duration and time-to-peak) from the two-dimensional plots. The identified spatial and temporal metrics presented differences between controls and children with CP. We have also identified similar temporal patterns spatially over the motor cortex, as will be shown in Section 3, which has great promise for differentiating between the two subject groups.

2.8.4.

Similarity

Image areas, which were spatially distinct from the primary activation region, were observed to have demonstrated temporal profiles of HbO change similar to those of the primary activation region. Because the amplitudes in HbO change of these secondary regions were lower in magnitude, they were not always assigned to the cluster of pixels with activation, as identified by the clustering algorithm discussed above.³⁶ To identify these areas, a similarity concept was defined that quantified aspects of resemblance between the temporal profile of the maximum activation area and other regions in the motor cortex.

Only pixels that had a value greater than or equal to zero from the time-averaged HbO activation images were compared in the similarity analysis. Before checking the similarity between the temporal profiles of each pixel, these were first smoothed to remove high-frequency instrumentation noise using a third-order Savitzky–Golay FIR filter.^{38, 39} After smoothing, each pixel was normalized to the peak activation value during the $40\text{-}\mathrm{s}$ tapping-rest interval because the similarity analysis focuses on the shape of activation curves and not their magnitude.

Subsequently, two metrics were defined for describing each pixel’s temporal activation profile. These were the angle defined by the tangent $\left(d\Theta \right)$ and the change of signal amplitude $\left(ds\right)$ at every time point, as shown in Fig. 3 .⁴⁰ A cost function was then defined to quantify the degree of similarity for these two metrics between the seed pixel (i.e., the peak activation pixel) and other pixels in the image. Though previous works have used the Euclidean distance to find the similarity between two time series,^{41, 42, 43} a dynamic time-warping (DTW)–based cost function was selected to measure similarity to accommodate for small changes in phase between the seed pixel and activation occurring elsewhere in the image.^{44, 45} A brief overview of how DTW was adapted to the fNIRS image analysis is explained next.

Fig. 3

(a) Measurements made for the change in angle defined by the tangent $\left(d\Theta \right)$ and the change in amplitude $\left(ds\right)$ for every time point, (b) illustration of a DTW path, and (c) histogram of the dynamic warping costs for all control subjects.

Suppose there are two time series, $S$ and $P$ , of length $N$ and $M$ , respectively, where $S=s_1$ , $s_2,\dots$ , $s_N$ and $P=p_1$ , $p_2,\dots ,p_M$ . An $N\times M$ matrix is then formed to align the two sequences in which the entry $(n,m)$ consists of the distance, $d(s_n,p_m)$ , of the two points $s_n$ and $p_m$ [i.e., ${(s_n-p_m)}^2$ ]. A warping path⁴⁵ is then found in which a contiguous set of matrix elements defines a mapping between $S$ and $P$ [Fig. 3] and that the warping cost, $\mathrm{DTW}(S,P)$ , is minimized such that

Eq. 2

$$\mathrm{DTW}(S,P)=\min \left[\sqrt{\left(\sum _{k=1}^K{w}_k\right)}\right],$$

where $K$ is the number of warping elements and $w_k$ is the distance found at ${(n,m)}_k$ . The warping path is subject to specific constraints such that the warping starts and finishes in diagonally opposite corners, steps are only allowed to adjacent cells (including diagonally adjacent cells), and each $w_k$ must be monotonically spaced in time.⁴⁵ Furthermore, the time warping is globally constrained using the Sakoe–Chiba band,⁴⁶ where $r=10$ (i.e., $0.5\mathrm{s}$ ). The purpose of this constraint is to prevent a relatively small section of one time series to be mapped to a relatively large section of another. In this case, a maximum of $0.5\mathrm{s}$ of one time series can be mapped to a single point of another (i.e., a shift of $0.5\mathrm{s}$ can be detected).

In this study, DTW was not directly applied to the pixel temporal responses, but to the two metrics that defined the shape of the temporal responses, $d\Theta$ and $ds$ . $N\times M$ matrices were formed to find the warping costs for both $d\Theta$ and $ds$ . Instead of looking for the minimum warping cost for $d\Theta$ and $ds$ separately, the warping cost was redefined such that $w_k$ in Eq. 2 is represented by

Eq. 3

$$w_k={(d{\Theta }_n-d{\Theta }_m)}^2+{(d{s}_n-d{s}_m)}^2.$$

Subsequently, a threshold needed to be defined to identify the range of DTW cost corresponding to activation image pixels that were considered similar to the seed pixel. Because the dynamic warping cost is directly comparable between controls and subjects with CP, the dynamic warping costs of the controls, for both right and left finger tapping, were pooled together in a single histogram that presented the frequency of each dynamic warping cost for all control subjects [e.g., Fig. 3]. Because there were no clear divisions found in the histogram, a $k$ -means algorithm was used to divide the histogram into five levels, in which the clusters were initialized similar to that of Liu and Yu.⁴⁷ Though the $k$ -means algorithm could divide the histogram into more levels for a more stringent similarity definition, dividing it into five levels was empirically found to be sufficient for finding similar time series, as will be shown in Section 3. A flowchart of the similarity algorithm is shown in Fig. 4 .

Fig. 4

Block diagram of the similarity algorithm.

2.9.

Statistical Analysis

Two-sampled t-tests were performed, using SAS 9.1 (SAS Institute Inc., Cary, North Carolina) to see if there was a significant difference $(p<0.05)$ of metric means between control subjects and ones with CP. In cases where the variances between the two subject groups were found to be significantly different $(p<0.05)$ , the Mann–Whitney test was used instead. The null hypothesis was defined as a zero difference in the mean between the two groups. For the cases where a statistically significant difference was found between the two subject groups, a post hoc power analysis was performed to validate the statistical power $\left(1\text{-}\beta \right)$ of the t-test, which indicates the probability of not having a false significant difference for the given subject population size.⁴⁸

3.1.

Distance from Center

Visual inspection of activation patterns in children with CP suggested these subjects had time-averaged activation areas closer to the centerline of the motor cortex compared to controls. The distance from center metric was therefore defined (Section 2.8.1) to measure this distance. Figures 5 and 5 are bar graphs showing the distances separately for the RH (right horizontal columns) and LH (left horizontal columns) for control subjects (Subj1–Subj5), ones with right hemiparetic CP (Subj6–Subj8), and ones with left hemiparetic CP (Subj9–Subj10). In Figs. 5 and 5, it is seen that controls consistently had activation further from the centerline in the ipsilateral hemisphere than in the contralateral hemisphere relative to the tapping hand. In contrast, the ipsilateral activation for children with CP was closer to the centerline, and in a few cases, it was even closer than the contralateral activation for the same subject.

Fig. 5

Bar graphs presenting the distance from center measurements from the time-averaged HbO activation images for (a) left finger tapping and (b) right finger tapping. The distances are presented separately for the RH (right horizontal columns) and LH (left horizontal columns) for control subjects (Subj1–Subj5, LH–black, RH–white), ones with right hemiparetic CP (Subj6–Subj8, LH–brick, RH–dots) and ones with left hemiparetic CP (Subj9–Subj10, LH–trellis pattern, RH–vertical lines). $\mathbf{x}$ indicates the exclusion of the measurement since no activation was identified in that hemisphere.

A high statistical significance was found $(p=0.0008$ , $1–\beta =0.98$ , Table 1 ) for the mean value difference in ipsilateral distance from center between controls $(7.29\mathrm{cm}\pm 1.98\mathrm{cm})$ and children with CP $(3.14\mathrm{cm}\pm 1.46\mathrm{cm})$ , but no significant difference $(p=0.54)$ was observed in the contralateral distance from center between controls $(3.80\mathrm{cm}\pm 1.48\mathrm{cm})$ and ones with CP $(3.30\mathrm{cm}\pm 2.06\mathrm{cm})$ .

Table 1

Efficacy of the derived metrics for differentiating between controls and subjects with CP. Comparisons are made between the right or left finger tapping (RFT/LFT) for controls versus tapping with the affected hand for subjects with CP.

3.2.

Area Difference

Although the distance from center metric gave information about differences in activation location between controls and children with CP, this metric did not capture differences in the spatial extent of the corresponding activation areas. To quantify the latter, the area-difference metric was defined as described in Section 2.8.2.

In Fig. 6 , the results of the area difference metric show that for controls the area of activation in the contralateral hemisphere of tapping was always greater than the activation in the middle area of the motor cortex by at least $2{\mathrm{cm}}^2$ . Three of the five subjects with CP (Subjects 6, 9, and 10) showed greater activation area in the middle portion of the motor cortex when tapping with one of their two hands. Interestingly, Subject 8, a right hemiparetic subject, presented an area difference of $<1{\mathrm{cm}}^2$ when right finger tapping (i.e., even though he was using his affected hand). Subject 7 had similar results as the controls for this metric. Although in some cases, subjects with CP presented larger activation in the middle region of the motor cortex than in the contralateral hemisphere, for the current population size there was no significant difference $(p=0.068)$ in area difference metrics between controls $(9.04{\mathrm{cm}}^2\pm 5.96{\mathrm{cm}}^2)$ and children with CP $(3.40{\mathrm{cm}}^2\pm 7.01{\mathrm{cm}}^2)$ . Also, there was no statistically significant difference in the mean values of the area difference metric when subjects with CP were tapping with their affected $(0.53{\mathrm{cm}}^2\pm 5.81{\mathrm{cm}}^2)$ or unaffected hand ( $6.27{\mathrm{cm}}^2\pm 7.51{\mathrm{cm}}^2$ , $p=0.2135$ ).

Fig. 6

Activation area difference for control subjects (Subj1–Subj5, left tapping–black, right tapping–white), ones with right hemiparetic CP (Subj6–Subj8, left tapping–brick, right tapping–dotted), and ones with left hemiparetic CP (Subj9–Subj10, left tapping–trellis pattern, right tapping–horizontal lines).

3.3.

Duration and Time-to-Peak

The ratio of duration over time-to-peak was identified as a potentially useful metric for differentiating between the two subject groups. With one exception, the controls had a duration over time-to-peak ratio of $⩾1$ for both right and left finger tapping (Fig. 7 ), i.e., duration was greater than or equal to time-to-peak. In contrast, duration was found to be less than the time-to-peak when children with CP were tapping with their affected hand (right, dotted column for Subjects 6–8 and left, trellis patterned column for Subjects 9 and 10). In addition, when children with CP were tapping with their unaffected hand, this ratio often assumed values that were substantially higher or lower than unity for some individuals, though no consistent pattern was observed.

Fig. 7

Bar graph presenting the ratio of duration over time-to-peak metrics for control subjects (Subj1–Subj5, left tapping–black, right tapping–white), ones with right hemiparetic CP (Subj6–Subj8, left tapping–brick, right tapping–dotted), and ones with left hemiparetic CP (Subj9–Subj10, left tapping–trellis pattern, right tapping–horizontal lines).

No significant difference $(p=0.99)$ was observed between the mean values of duration over time-to-peak ratios between controls that were left finger tapping $(1.19\pm 0.09)$ or right finger tapping $(1.19\pm 0.26)$ . Subsequently, these results were pooled into a single normal subject group and a two-sampled t-test verified the high statistical significance ( $p=0.0002$ , $1–\beta =0.99$ , Table 1) in the difference between the means of the duration over time-to-peak ratios between controls $(1.19\pm 0.18)$ and children with CP tapping with their affected hand $(0.60\pm 0.25)$ .

3.4.

Similarity

By using the above-described similarity algorithm (Section 2.8.4), pixel areas were identified having similar temporal profiles of HbO concentration change with respect to the primary area of activation. It was found that these similar areas did not necessarily coincide with the areas having the highest activation in the time-averaged images, but could have an HbO concentration change in the baseline range. Furthermore, subjects with CP demonstrated a greater diversity of temporal patterns, as quantified by the range of DTW penalty term values, within the areas of activation than control subjects. To illustrate these points, a sample activation image [Fig. 8 ], an image processed by an activation cluster algorithm [Fig. 8], an image processed by the similarity algorithm [Fig. 8], and the color-coordinated normalized, average temporal profiles of the areas in the similarity image [Fig. 8] are shown for a control subject. For purposes of comparison, the same types of images are shown for a subject with CP [Figs. 8, 8, 8, 8]. It is seen that the area of activation identified by the clustering algorithm [Fig. 8] does not coincide with the area of pixels with similar activation [Fig. 8]. In contrast, children with CP had areas of activation [Fig. 8] with dissimilar activation profiles, except for pixels immediately neighboring the seed-pixel location [Fig. 8]. The temporal profiles of activation are also shown [Figs. 8 and 8], in order to provide a visual example of the degree of similarity between the seed pixel (black curves) and the most similar cluster of pixels (level 1). Because the dynamic warping costs of the control subjects were pooled together, not all subjects individually present all five levels of similarity as seen in Figs. 8 and 8.

Fig. 8

Images shown are for a control subject, Subject 2, in which the (a) normalized image of the average change in HbO concentration, (b) the thresholded activation image, (c) the similarity image with color-coordinated levels, and (d) the corresponding average temporal responses for the color-coordinated similarity levels are shown. Not all five levels of similarity are seen for this control subject. The images at the bottom are for a subject with CP, Subject 7, in which (e–h) correspond to (a–d) for the control subject.

The area difference of the most similar pixels to the primary activation area (level 1) in the similarity images was investigated as a metric for differentiating between controls and children with CP. In Fig. 9 , it is seen that controls typically had a greater area difference than children with CP. This difference between the mean values of area difference of controls $(12.14{\mathrm{cm}}^2\pm 12.85{\mathrm{cm}}^2)$ and children with CP $(-0.29{\mathrm{cm}}^2\pm 4.66{\mathrm{cm}}^2)$ was indeed found to be significant ( $p=0.0383$ , $1–\beta =0.77$ , Table 1). The results for area difference demonstrate that the similarity concept can highlight differences between controls and subjects with CP much better than the corresponding metric derived from applying a threshold algorithm to the same subject’s activation images. Additionally, there was no discernible correlation between area difference and whether children with CP were tapping with their affected $(1.82{\mathrm{cm}}^2\pm 5.09{\mathrm{cm}}^2)$ or unaffected hand ( $1.25{\mathrm{cm}}^2\pm 4.14{\mathrm{cm}}^2$ , $p=0.3259$ ).

Fig. 9

Area difference in similarity images for control subjects (Subj1–Subj5, left tapping–black, right tapping–white), ones with right hemiparetic CP (Subj6–Subj8, left tapping–brick, right tapping–dotted), and ones with left hemiparetic CP (Subj9–Subj10, left tapping–trellis pattern, right tapping–horizontal lines).

Although children with CP had smaller areas of similarity, their most similar areas (level 1) were spatially closer to their seed pixel than in controls (Fig. 10 ). To present this difference the maximum distance among the most similar areas (level 1) from the seed pixel was found, in which the difference between the controls $(7.43\mathrm{cm}\pm 6.00\mathrm{cm})$ and children with CP $(2.18\mathrm{cm}\pm 1.17\mathrm{cm})$ was significant ( $p=0.0223$ , $1–\beta =0.73$ , Table 1). The larger distance in controls can be explained by some similarity areas being seen in the ipsilateral hemisphere and the larger areas of similarity neighboring the seed pixel, whereas in children with CP this is not the case. Furthermore, there is no significant difference $(p=0.7427)$ when the children with CP were tapping with their affected $(2.05\mathrm{cm}\pm 1.53\mathrm{cm})$ or unaffected hand $(2.31\mathrm{cm}\pm 0.84\mathrm{cm})$ .

Fig. 10

Maximum distance of the pixel with the highest similarity (level 1) to the seed pixel in similarity images for control subjects (Subj1–Subj5, left tapping–black, right tapping–white), ones with right hemiparetic CP (Subj6–Subj8, left tapping–brick, right tapping–dotted), and ones with left hemiparetic CP (Subj9–Subj10, left tapping–trellis pattern, right tapping–horizontal lines).

3.5.

Validation Analysis to Differentiate between Controls and Subjects with CP

Though four of the five presented metrics gave statistically significant differences between controls and children with CP (distance from center, duration/time-to-peak, area difference of similarity images, and maximum distance from seed), there were still false positives or false negatives found. Table 1 shows the sensitivity, specificity, positive predictive value (PPV), and negative predictive value (NPV) for each of these metrics. Subjects were considered to be in the CP group if their metrics were beyond one standard deviation of the corresponding mean value of the same metric for controls.

From Table 1, it is clear that each fNIRS metric performs differently in each category. For example, the specificity, sensitivity, PPV, and NPV values for area difference in activation images are low, which supports the t-test conclusion $(p\text{-}\text{value}=0.068)$ that this metric may be less useful for differentiating between the two subject populations. On the other hand, distance from center and distance from seed clearly tell if a subject is a control ( $\mathrm{PPV}=80\mathrm{\%}$ , $\mathrm{NPV}=100\mathrm{\%}$ ), while duration/time-to-peak clearly tells if a subject is in the CP group ( $\mathrm{PPV}=100\mathrm{\%}$ , $\mathrm{NPV}=80\mathrm{\%}$ ). However, the area difference of similarity does not clearly show if a subject is a control or with CP, but it does have a higher specificity (80%) than the area difference in activation images (60%). Because no single metric has 100% success, an approach for distinguishing between the two subject populations where multiple metrics are quantified is warranted.

3.6.

Repeatability

The exact anatomical placement of probes and the quality of optical contact with the scalp are common causes of concern for the repeatability of fNIRS measurements. To reduce such concerns in this study, measurements were performed twice for four subjects in each group (Table 2 ). As an example, the reflectance time series for the source-detector pair with maximum signal, which was the same for both visits, are presented in Figs 11 and 11 for a control subject (Subject 2) and for a subject with CP (Subject 8) tapping with their right hand. It is seen that these time-series signals, as well as the reconstructed activation and activation-threshold images [Figs. 11 and 11] have similar patterns for both visits. Furthermore, spatial metrics (distance from center and area difference) were quantified for both visits and a paired t-test was performed on these metrics for each subject group. For both groups, there was no significant difference in these metrics between the two visits (Table 2, X is shown when no activation was found in the ipsilateral hemisphere). The repeatability findings were similar for all other metrics (not shown).

Fig. 11

Repeated fNIRS measurements for a control (Subject 2) and a subject with CP (Subject 8) performing right finger tapping: (a,b) show the time-series reflectance patterns for the same source-detector pair, while (c,d) show the activation images (top) and their corresponding activation-threshold images (bottom) for each visit.

Table 2

The ipsilateral distance from center and area difference metrics measured in two separate visits for four subjects in each group.