2.1.

Participants

Fifteen subjects (9 males and 6 females, Caucasian, $\text{mean age}\pm \mathrm{SD}$: $27\pm 4.6\text{years}$) participated in our study. Only subjects that fulfilled all inclusion criteria for the use of TMS³² were recruited. Subjects 2, 4, and 15 were left-handed, the other twelve subjects were right-handed. The body-mass index was calculated for every subject, and hair root density, hair color. and hair thickness were assessed on a scale from 0 to 5 (0: no hairs, 5: dense, thick, dark hairs) upon visual inspection by the same experimenter (LS).

The experiments were approved by the ethical commission of ETH Zurich (2018-N-22), the Cantonal Ethics Commission of the Canton of Zurich (2018-01078), and were conducted in accordance with the Declaration of Helsinki. Written informed content was provided by all subjects.

2.2.

fNIRS Instrumentation

A custom-built fNIRS instrument called optoHIVE was used to detect cortical brain activity.^23,33 It is a lightweight, fiber-less system designed for wearable, high-quality measurements, and each optode includes a four-wavelength LED light source (774, 817, 865, and 892 nm) and a silicon photomultiplier for photodetection.³³ This system was preferred over commercial fNIRS instruments because, in addition to wearability, it offers the advantages of modular optode placement, a large number of short-distance channels, and high optical sensitivity.²³ Each optode module (containing a light source and detector) comprises a short-separation (SS) channel of 7.5 mm and measures with every other optode within 30 mm distance, as shown in Fig. 1(a). In this study, eight optode modules were fixed using a custom-built headgear made of silicone patches, three-dimensional-printed parts, and elastic strings [see Fig. 1(c)]. The optodes were symmetrically placed over the left and right primary motor cortices (M1), ventral premotor cortices, and dorsal premotor cortices, resulting in 16 long-separation (LS) channels of 30 mm and 8 SS (7.5 mm) channels. As always, two channels measured from the same brain location (i.e., their light-paths were overlapping), and the 16 channels were reduced to eight ROIs. A map of the optode configuration and the ROIs is shown in Fig. 1(b). Measurements performed by optoHIVE are controlled through a LabVIEW (Version 2015, National Instruments NI, Texas) interface, and data from each channel are collected at a sampling frequency of 8.98 Hz over an NI myRIO data acquisition device.³³

Fig. 1

Experimental setup. (a) Sensitivity maps of SS and long-separation measurements for two optodes. Two-dimensional sensitivity maps were obtained from Monte Carlo simulations using ValoMC.^29,34,35 (b) Arrangement of the optoHIVE optodes (red) and ROIs (green) according to the 10–20 system of EEG placement. Each ROI consists of two long-separation channels that probed the same brain regions. (c) Subjects were seated with arms resting comfortably on cushioned armrests and the hands grasping a foam handlebar. If desired, a support pillow was added to additionally support the arms. OptoHIVE was placed over the left and right motor areas to record fNIRS signals. Different biosignals were concomitantly acquired (not used in this work).

2.3.

Study Protocol

2.4.

Data Processing

Data processing was performed in MATLAB (R2017a, Mathworks Inc.). Motion artifacts were removed using spline interpolation.³⁷ Raw optical intensities were converted to concentration changes of ${\mathrm{O}}_2\mathrm{Hb}$ ($[{\mathrm{O}}_2\mathrm{Hb}]$) and HHb ([HHb]) using the modified Beer–Lambert law. The absorption coefficients were adopted from Moaveni,^38,39 and the differential pathlength factors for the four wavelengths (6.2, 6.2, 5.9, 5.5) were from Cope.⁴⁰ For the removal of drift and cardiac pulsation of $[{\mathrm{O}}_2\mathrm{Hb}]$ and [HHb], different methods were considered according to Pinti et al.,⁴¹ but the best results were achieved following their suggested optimal filter (finite impulse response, order 1000) at cutoff frequencies of 0.015 and 0.35 Hz. Building on our previous work,²⁷ multichannel SCR based on non-negative least squares (${\mathrm{GLM}}^{\mathrm{multiSS}}$) was applied to best reduce the influence of physiological changes and, thus, to separate the hemodynamic changes from the extracerebral tissue layer. With ${\mathrm{GLM}}^{\mathrm{multiSS}}$, all short-channel distances are included as regressors in the general linear model (GLM) with the precautional measure of allowing for only positive estimates (non-negative least squares regression). Short channel signal quality was verified following the approach of Perdue et al.,⁴² which is based on the signal content (heart rate) rather than purely on the SNR.²³ It was found that 90% of the short channels were of good quality during our measurements when considering a signal quality threshold of 12 dB.²⁷ Amplitudes of MW oscillations were obtained from [${\mathrm{O}}_2\mathrm{Hb}$] by normalizing the band-power (0.07 to 0.14 Hz) with its pulse band-power (0.6 to 2 Hz), and the median value of all long-separation channels was extracted.²⁷

A GLM^27,43 was applied on the time course of the long-separation channel (i.e., with and without SCR). As an evaluation metric, $t$-values were obtained. The $t$-values give an indication on the signal strength of a fitted hemodynamic response curve in relation to the residuals. The used GLM consisted of a modeled hemodynamic response time course, obtained from the convolution of the boxcar function and the canonical hemodynamic response,⁴³ its time and dispersion derivatives, and a constant offset, which were fitted into the fNIRS data of each recording channel. The time and dispersion derivatives⁴⁴ were included to correct for deviations of the onset and the shape of the hemodynamic response, respectively. The $t$-values were stored in a vector with 1920 entries (15 subjects $\times 16$ channels $\times 2$ hands $\times 2$ runs $\times 2$ sessions) for [${\mathrm{O}}_2\mathrm{Hb}$] and [HHb]. When ROI analysis was performed, the average of the $t$-values of the corresponding channels was used. GLM analysis was applied on non-regressed (NR) long-separation measurements and on SCR data.

2.5.

Statistical Analysis

Statistical analysis was performed in R (Version 3.6.3, RStudio Inc.).⁴⁵ To find a threshold to distinguish between active and inactive channels, GLM was applied on the baseline data (i.e., random task onsets during rest condition without systemic brain activity), and $t$-values were extracted. From the obtained $t$-values distributed around 0, the threshold, below which the probability is $>95\%$ that the brain is in rest condition, was extracted. Consequently, the 5% significance level to indicate if a hemodynamic response, representing brain activity above chance, was present, was found to be $t\ge 30$ for SCR data and $t\ge 22$ for unregressed data.

Reproducibility between sessions was assessed using linear correlation analysis and intraclass correlation coefficient (ICC) analysis applied to the $t$-values provided by the GLM. Linear correlation between sessions 1 and 2 on a group level was calculated based on Pearson correlation coefficients applied on the $t$-values of the M1. For this purpose, the $t$-values for left M1 (right hand task) and right M1 (left hand task) were extracted, averaged for the two runs per session, and correlated between sessions for the 15 subjects. Test-retest reliability was determined using ICC based on an absolute agreement, two-way random effects model with repeated measures.^17,46 Single (ICC(2,1)) and average (ICC(2,k)) measures and their 95% confidence intervals were calculated as suggested by Li et al.¹⁷ The ICC gives an indication of the reliability of measurements by comparing the variability of different tests of the same individuals with the total variation across all ratings and all individuals. A high ICC (close to 1) indicates low intrasubject variability relative to the intersubject variability, whereas a low ICC (close to 0) means that values from the same group are not similar.⁴⁶ Thresholds for interpreting ICCs vary in literature; we used the definition according to Li et al.:¹⁷ poor ($\mathrm{ICC}<0.40$), fair ($0.40\le \mathrm{ICC}<0.60$), good ($0.60\le 0.75$), and excellent ($0.75\le \mathrm{ICC}<1.00$). To estimate the change in $t$-values across measurement sessions, the mean absolute scaled error (MAE%) between sessions 2 and 1 was calculated and normalized with respect to the range of observed values.

A linear mixed effects model with restricted maximum likelihood estimation (lmer in $R$) was applied on the $t$-values to investigate the statistical significance of factors that could affect the estimation of brain activation. Two mixed effects models were established, one for [${\mathrm{O}}_2\mathrm{Hb}$] and one for [HHb]. The dependent variable consisted of the 1920 $t$-values (channelwise). The fixed effects were selected as (1) the interaction between hand and channel hand*channel, (2) run, (3) raw signal strength, and (4) MW amplitude. As random variables, an intercept for subject and a nested random effect of channel per subject were considered to allow levels of $t$-value to vary across channels and subjects. Square-root transformation was applied to ensure that the model residuals were normally distributed. The goodness of fit was verified from a normal distribution and homoscedasticity of the model residuals. After model fitting, the estimates were backtransformed, and their effect was investigated by multiple comparisons with Tukey contrasts.^47,48

2.6.

Classification

A classifier was trained and tested on the obtained fNIRS data to link the aspect of reproducibility with a potential application such as a BCI. Therefore, pseudo-online classification (i.e., continuous samplewise processing of data) was performed using a support vector machine (SVM) with $L1$-norm regularization⁴⁹ to classify between right and left hand grasping. First, [${\mathrm{O}}_2\mathrm{Hb}$] and [HHb] were bandpass-filtered with a forward filter between the cutoff frequencies 0.005 Hz (Chebyshev type II, order 2) and 0.35 Hz (Butterworth, order 4). Second, adaptive filtering based on non-negative weight estimation was used for the samplewise regression (i.e., SCR) of the SS channels from the long-separation channels. Third, feature extraction was performed by extracting three feature types from [${\mathrm{O}}_2\mathrm{Hb}$] and [HHb]: amplitude, slope, and correlation-based signal improvement.⁵⁰ Fourth, an $L1$-norm SVM was used for joint training of the classifier and feature selection.⁵¹ Its regularization parameter was found with fivefold cross-validation on the training data. The classifier was trained on the first run and tested on the second run of the same session. The intrasession accuracy was calculated over the 30 trials per testing run for sessions 1 and 2, separately. The entire 16 s trial window was used for classification, and a decision was obtained at the end of each trial. Although other window sizes or features could be considered, we selected the entire duration to keep the processing pipeline simple as similar results would be expected with a shorter window and delay until a decision is made was not deemed critical during this study. Three classifier scenarios were trained to investigate the effect of the selected dataset (i.e., SS signals only, long-separation only, or SCR signals) on classification accuracy.

3.1.

Spatial Activation Patterns

On a group level, spatially specific patterns of brain activity were detected for the left and right grasping tasks, shown in Fig. 2. The strongest hemodynamic response for the right hand task was observed over the contralateral (left) M1 and vice versa for the left hand task. In particular, ROI 1 showed the strongest magnitude for the right hand task with an average change of $0.26/-0.12\mu \mathrm{M}$ for [${\mathrm{O}}_2\mathrm{Hb}$] and [HHb], whereas for the left hand task, ROI 5 changed the strongest with $0.21/-0.09\mu \mathrm{M}$ ([$[{\mathrm{O}}_2\mathrm{Hb}]/[\mathrm{HHb}]$). Weaker activation was observed at the adjacent frontal brain areas (ROI 2+3 for the right hand and ROI 6+7 for the left hand) and the ipsilateral M1 (ROI 5 for the right hand and ROI 1 for the left hand). The remaining ROI (ROI 4+6+7+8 for the right hand task and ROI 2+3+4+8 for the left hand task) did not exhibit significant activation. The hemodynamic responses for the NR signals had slightly larger magnitudes mainly for [${\mathrm{O}}_2\mathrm{Hb}$] than for the SCR signals with $0.33/-0.13\mu \mathrm{M}$ versus $0.26/0.12\mu \mathrm{M}$ for the right hand task and $0.27/-0.10\mu \mathrm{M}$ versus $0.26/0.09\mu \mathrm{M}$ for the left hand task. This observation was accompanied by a smaller standard deviation of the signal. In comparison with the NR signal, the SCR signal did not exhibit an oscillatory response after task onset.

Fig. 2

Group average of [${\mathrm{O}}_2\mathrm{Hb}$] and [HHb]. The hemodynamic responses of the 15 subjects were averaged ($\text{mean}\pm \mathrm{SD}$) for each ROI for (a) the right and (b) left hand grasping tasks. The gray bars indicate the task period when the grasping task was conducted. The colored lines represent relative concentration changes of ${\mathrm{O}}_2\mathrm{Hb}$ (red) and HHb (blue) when SCR was applied, with the shaded areas indicating the standard deviation (SD). The black (dotted) lines represent the mean (and SD) of the unregressed (NR) signal. The numbers in the upper left corners indicate the ROI. Units are in $\mu \mathrm{M}$.

On a single-subject level (see A), task-evoked brain activation became visible for nine out of 15 subjects, whereas six subjects showed only minimal or no activation (subjects 5, 6, 7, 8, 14, and 15). These qualitative results can be put into relation with the individual $t$-values (visible in Fig. 3 or 5), which were the lowest for the six subjects with weak activation. More specifically, the $t$-values for subjects with weak hemodynamic responses were in the range of 14 to 24 after SCR, which is below the threshold of 30 (as determined from baseline measurements). The lowest $t$-value of the other subjects was 31. Subsequently, the nine subjects with distinct spatial activation and $t$-values $>30$ are denoted “strong responders,” and the other six subjects with $t$-values $\le 30$ are “weak responders.”

Fig. 3

Correlation between test and retest sessions. Reproducibility of $t$-values between the two sessions was assessed for [${\mathrm{O}}_2\mathrm{Hb}$] and [HHb] of the right and left hand grasping tasks. $t$-values give an indication of the quality of the measured hemodynamic response. For each task, the correlation plots for the ROIs over the contralateral M1 (i.e., left M1 for right hand grasping and right M1 for left hand grasping) are shown. (a) Results for NR $t$-values. (b) Results for regressed (after SCR). Bright and dark gray areas indicate that the $t$-values exceed baseline noise (i.e., 22 for no regression and 30 for SCR) for one or both sessions, respectively. Pearson’s correlation coefficient, its confidence bounds, and $p$-values between days are displayed in the upper left corner of each scatter plot. Data points are labeled with the subject number, and the red dashed line indicates the confidence bounds.

3.2.

Reproducibility

In the correlation plots in Fig. 3, the agreement between the test (session 1) and retest (session 2) sessions among all subjects is presented, taking into account $t$-values from single subjects. The test-retest agreement changed depending on the chromophore (${\mathrm{O}}_2\mathrm{Hb}$ or HHb, respectively) and signal processing step. The lowest agreement with correlations of 0.71/0.51 (left/right hand) was observed for [${\mathrm{O}}_2\mathrm{Hb}$] when no regression was applied and increased after SCR to 0.81/0.54. Correlation for [HHb] was only slightly higher for SCR over NR with 0.81/0.70 to 0.81/0.64. A distinct difference between right and left hand grasping was observed, with the left hand task always scoring more than 0.1 points less in correlation coefficients. The MAE% in Table 1 supports this observation, as the MAE% was generally lower for the right hand grasping task than the left hand grasping task (20.1% versus 13.4% for SCR ${\mathrm{O}}_2\mathrm{Hb}$).

Table 1

ICC and MAE% single and average ICCs are given for ROIs above contralateral M1 for right (ROI 1) and left (ROI 5) hand grasping.17 Results for NR and SCR data are shown. The t-values of the runs per session were averaged. In brackets, the 95% confidence intervals are given. MAE% indicates the mean absolute scaled error relative to the range of the data.

Test-retest reliability based on ICC is presented in Table 1. The contralateral fNIRS channel above M1 was investigated for each task, i.e., left M1 for the right hand task and right M1 for the left hand task. The ICC values ranged between 0.58 and 0.81 for the single measurements and 0.85 and 0.94 for the averaged metric. It was found that the ICCs depend on the chromophore ([${\mathrm{O}}_2\mathrm{Hb}$] and [HHb]) and processing steps (NR and SCR). There was a trend that ICCs were higher after SCR in comparison with NR, as the ${\mathrm{ICC}}_{\text{single}}$ values were fair–good for NR and increased to fair–excellent for SCR. The ${\mathrm{ICC}}_{\text{average}}$ was excellent for all investigated combinations. The ICCs were higher for the right hand than the left hand grasping task.

3.3.

Linear Mixed Effects Model

The influence of different variables on $t$-value estimation was determined by two linear-mixed effects models applied on [${\mathrm{O}}_2\mathrm{Hb}$] and [HHb].

Analysis of variance of [${\mathrm{O}}_2\mathrm{Hb}$] revealed a strong significant main effect of the hand*channel interaction ($F_{(\mathrm{7,827})}=21.78$, $p<.001$). Also, signal strength ($F_{(\mathrm{1,235})}=19.88$, $p<.001$), run ($F_{(\mathrm{3,843})}=3.04$, $p<.05$), and MW amplitude ($F_{(\mathrm{1,925})}=5.80$, $p<.05$) were found to have a significant effect on $t$-values. For the right hand grasping (${\mathrm{O}}_2\mathrm{Hb}$), the highest activation was found in left M1, and the lowest activity was detected in right M1. The left hand condition (${\mathrm{O}}_2\mathrm{Hb}$) had the highest activity in right M1 and the lowest in left M1.

The results for HHb were similar to ${\mathrm{O}}_2\mathrm{Hb}$. A strong significant main effect of hand*channel interaction was found ($F_{(\mathrm{7,827})}=15.52$, $p<.001$). Signal strength ($F_{(\mathrm{1,239})}=13.47$, $p<.001$), run ($F_{(\mathrm{3,842})}=6.23$, $p<.001$) and MW amplitude ($F_{(\mathrm{1,921})}=17.68$, $p<.001$) had a highly significant effect on $t$-values. The highest activity for the right hand (HHb) was found in left M1, and the lowest activity was detected in right M1. For the left hand grasping (HHb), the highest activity was in right M1 and the lowest in left M1.

3.4.

Confounding Factors

The effect of signal strength and MW amplitude on $t$-values is visualized in Fig. 4. In Fig. 4(a), the raw signal magnitudes are plotted against $t$-values, with the maximal values for all runs averaged per subject. For better visibility, log-transformation was applied on the signal strength. A trend can be observed that smaller signal strength correlates with smaller $t$-values. In particular, four out of six weak responders (subjects 6, 7, 14, and 15) had low signal strength. An outlier was subject 5, which had low $t$-values but high signal strength. In Fig. 4(b), MW amplitude is compared with $t$-values. When ignoring the subjects with low signal strength (subjects 6, 7, 14, and 15), a negative correlation between $t$-values and MW amplitude is observed.

Fig. 4

Influence of signal strength and MW amplitude on brain activity estimates. (a) $t$-values are plotted against signal magnitude for NR and SCR signals. (b) $t$-values are plotted against the median MW amplitude for NR and SCR signals. The highest $t$-values of the two runs were averaged per session. The labels correspond to the subject number. Pearson’s correlation coefficient, confidence bounds, and $p$-values are displayed in the upper left corner of each scatter plot.

Figure 5 further presents influencing factors in relation to $t$-values on a single-subject level. It again becomes visible that the subjects with high $t$-values often have high signal strength and low Mayer-wave amplitudes. High $t$-values are closely linked to high signal strength, which is dependent on hair characteristics: subjects with dense and dark hair tend to have lower signal strength than those with blonde and thin hair. A correlation between classification accuracy and $t$-values is observed, with subjects having high $t$-values also scoring high classification accuracy.

Fig. 5

Subject-specific signal assessment. Data are sorted according to SCR $t$-values from lowest to highest subject-specific values (first row). The maximal $t$-value of the four runs per subject were averaged. The red lines for $t$-values indicate the threshold of 30 and 22 for the separation of subjects into strong and weak activators, which were found based on a 5% level of significance above chance. Signal strength is given in voltage, and the red line defines the threshold of 0.06 V, which is required for measurements with sufficient quality. In the last row, classification accuracy per subject averaged over the two testing runs is shown, with the red line at 70% indicating the 1% level of significance above chance. SCR, short-channel regression; NR, no regression.

3.5.

Classification

Classification accuracies were in good agreement with the results from the previous analyses: strong responders had high accuracies and weak responders had low accuracies. In the bottom line of Fig. 5, the single-subject classification accuracy graphically presents this trend as a dependency of $t$-values. There were seven subjects with accuracies >95% (subjects 1, 2, 3, 4, 9, 12, and 13). When separating the subjects into strong and weak responders, a significant difference between the two groups is observed [95% versus 69%, Fig. 6(a)]. No significant difference was observed between the classification accuracies of the test and retest sessions [Fig. 6(b)]. When feeding different input data to the classifier, there was no significant advantage of using SCR over NR data in terms of classification accuracy. When only SS data were used as input for the classifier, as a way to validate the absence of brain activity in the SS measurements, the accuracy was close to 50%, which corresponds to chance level and confirms the assumption that no brain activity is present in SS channels.

Fig. 6

Classification accuracies. (a) Strong versus weak activators (SCR signals). (b) First versus second sessions (SCR signals). (c) Different input signals for the classifier. In particular, signals from SS and long-separation NR channels, as well as SCR signals, were used as input for the classifier pipeline. The red line indicates the 1% significance level for the classifier.⁵² **: $p<0.001$ using a two-tailed, paired $t$-test.

6.

Appendices

6.1.

A Appendix

In Fig. 7, correlation between the test and retest sessions is shown on the group level with the average $t$-values over all subjects per ROI. This plot confirms the findings from Plichta et al.⁶⁰ that fNIRS are highly reproducible across sessions for group-averaged optode locations.

Fig. 7

Group-level correlation of the eight ROIs between test and retest sessions. Reproducibility of $t$-values between the two sessions was assessed for [${\mathrm{O}}_2\mathrm{Hb}$] and [HHb] of the right and left hand grasping tasks. $t$-values give an indication of the quality of the measured hemodynamic response. For each task, the correlation plots for the eight ROIs averaged over all subjects are shown. (a) Results for NR $t$-values. (b) Results for regressed (after SCR). Gray areas indicate the $t$-value range of statistical significance (i.e., 22 for no regression and 30 for SCR). Pearson’s correlation coefficient, its confidence bounds, and $p$-values between days are displayed in the upper left corner of each scatter plot. Data points are labeled with the ROI number, and the red dashed line indicates the confidence bounds.

In Figs. 8Fig. 9–10, the block averages of each subjects are shown for the right and left hand grasping tasks. Individual patterns for individual subjects become apparent. For example, there were subjects (e.g., S2, S3, and S4) exhibiting strong and easily visible hemodynamic responses, others (e.g., S1, S5, S15) having strong task-evoked systemic activity, and some (e.g., S7, S8, S14) showing no to minimal activation.

Fig. 8

Block average of [${\mathrm{O}}_2\mathrm{Hb}$] and [HHb] for subjects 1 to 5. The hemodynamic responses of the four runs per subject were averaged ($\text{mean}\pm \mathrm{SD}$) for each ROI. The spatial patterns for the right and left hand grasping tasks are shown. The gray bars indicate the task period when grasping with either the left or right hands was conducted. Units are in $\mu \mathrm{M}$.

Fig. 9

Block average of [${\mathrm{O}}_2\mathrm{Hb}$] and [HHb] for subjects 6 to 10. The hemodynamic responses of the four runs per subject were averaged ($\text{mean}\pm \mathrm{SD}$) for each ROI. The spatial patterns for the right and left hand grasping tasks are shown. The gray bars indicate the task period when grasping with either the left or the right hand was conducted. Units are in $\mu \mathrm{M}$.

Fig. 10

Block average of [${\mathrm{O}}_2\mathrm{Hb}$] and [HHb] for subjects 11 to 15. The hemodynamic responses of the four runs per subject were averaged ($\text{mean}\pm \mathrm{SD}$) for each ROI. The spatial patterns for the right and left hand grasping tasks are shown. The gray bars indicate the task period when grasping with either the left or the right hand was conducted. Units are in $\mu \mathrm{M}$.