Front. Neurorobot. Frontiers in Neurorobotics Front. Neurorobot. 1662-5218 Frontiers Media S.A. 10.3389/fnbot.2020.00010 Neuroscience Original Research Reduction of Onset Delay in Functional Near-Infrared Spectroscopy: Prediction of HbO/HbR Signals Zafar Amad 1 2 Hong Keum-Shik 1 3 * 1School of Mechanical Engineering, Pusan National University, Busan, South Korea 2Department of Electrical Engineering, University of Wah, Wah Cantonment, Pakistan 3Department of Cogno-Mechatronics Engineering, Pusan National University, Busan, South Korea

Edited by: Ganesh R. Naik, Western Sydney University, Australia

Reviewed by: Michele Luigi Pierro, Vivonics, United States; Uma Shahani, Glasgow Caledonian University, United Kingdom

 *Correspondence: Keum-Shik Hong kshong@pusan.ac.kr 18 02 2020 2020 14 10 22 08 2019 30 01 2020 Copyright © 2020 Zafar and Hong. 2020 Zafar and Hong

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

An intrinsic problem when using hemodynamic responses for the brain-machine interface is the slow nature of the physiological process. In this paper, a novel method that estimates the oxyhemoglobin changes caused by neuronal activations is proposed and validated. In monitoring the time responses of blood-oxygen-level-dependent signals with functional near-infrared spectroscopy (fNIRS), the early trajectories of both oxy- and deoxy-hemoglobins in their phase space are scrutinized. Furthermore, to reduce the detection time, a prediction method based upon a kernel-based recursive least squares (KRLS) algorithm is implemented. In validating the proposed approach, the fNIRS signals of finger tapping tasks measured from the left motor cortex are examined. The results show that the KRLS algorithm using the Gaussian kernel yields the best fitting for both ΔHbO (i.e., 87.5%) and ΔHbR (i.e., 85.2%) at q = 15 steps ahead (i.e., 1.63 s ahead at a sampling frequency of 9.19 Hz). This concludes that a neuronal activation can be concluded in about 0.1 s with fNIRS using prediction, which enables an almost real-time practice if combined with EEG.

 hemodynamic response prediction tracking vector phase analysis brain-machine interface (BMI) functional near-infrared spectroscopy (fNIRS) NRF-2017 R1A4A1015627 NRF-2017R1A2A1A17069430 National Research Foundation of Korea10.13039/501100003725 Introduction

Similar to functional magnetic resonance imaging and electroencephalography (EEG), functional near-infrared spectroscopy (fNIRS) is a non-invasive neuroimaging technique that measures hemoglobin oxygenation changes in the brain (Kato et al., 1993; Villringer et al., 1993). fNIRS can measure the absolute as well as relative concentration changes of oxyhemoglobin (HbO/ΔHbO) and deoxyhemoglobin (HbR/ΔHbR) using multiple near-infrared lights within the range of 650~1,000 nm (Pellicer and Del Carmen Bravo, 2011; Boas et al., 2014; Nguyen et al., 2016). It offers several advantages, including acceptable temporal and spatial resolution (Hong and Naseer, 2016; Nguyen and Hong, 2016), portability, and low cost (Ferrari and Quaresima, 2012). With these advantages, fNIRS has successfully demonstrated its potential as a viable neuroimaging tool for applications to the health care industry (Hong and Yaqub, 2019), neurological disorders (Obrig, 2014; Ghafoor et al., 2019; Yang et al., 2019), psychiatric disorders (Ohi et al., 2017), behavioral and cognitive development (Watanabe et al., 2017; Yaqub et al., 2018), and brain-computer interfaces (BCIs) (Nicolas-Alonso and Gomez-Gil, 2012; Naseer and Hong, 2015; Schudlo and Chau, 2018; Shin and Im, 2018).

The measured fNIRS signals (i.e., ΔHbO, ΔHbR) can be categorized into three durations (Frostig et al., 1990; Ernst and Hennig, 1994): (i) the initial dip, which represents the early extraction of oxygen by the nearby active neurons causing the ΔHbO/ΔHbR to decrease/increase, (ii) the conventional hemodynamic response (HR) that is the large increase in cerebral blood flow (CBF) resulting in an increase/decrease in ΔHbO/ΔHbR, respectively, and (iii) the undershoot before going back to the rest state. The changes in ΔHbO/ΔHbR upon the functional stimulation can be translated into meaningful commands for BCI applications (Matthews et al., 2008). These converted signals can be further used to actuate external devices such as robotic arm/leg or wheelchairs for improving the quality of patient lives (Mcfarland and Wolpaw, 2010, 2011; Ortiz-Rosario and Adeli, 2013; Yazdani et al., 2018). In particular, fNIRS devices are portable and have shown great potential for BCI applications. The main limitation of fNIRS for BCI is its slow nature of the HR and the inherent delay from the onset of the neuronal activity (Jasdzewski et al., 2003; Cui et al., 2010; Ahn and Jun, 2017), which restricts its use for online BCI applications as well as hybridization with other rapid techniques such as EEG (Jiao et al., 2018; Li et al., 2018; Yang et al., 2018), magnetoencephalography, etc. Because of this limitation, various features in different temporal windows of 0–5, 2–7, 0–10, 0–15, 0–17, and 0–20 s were used in multi-class classification algorithms to classify HRs associated with the same or different brain regions for fNIRS-BCI applications (Power et al., 2011; Khan et al., 2014; Schudlo and Chau, 2014; Gateau et al., 2015; Khan and Hong, 2015; Hong et al., 2017; Shin et al., 2017; Liu et al., 2018; Yi et al., 2018). Thus far, the features frequently used from these windows include signal mean, signal slope, signal peak, skewness, kurtosis, variance, standard deviation, number and sum of peaks, root mean square, median, etc. (Hwang et al., 2016; Naseer et al., 2016; Liu and Hong, 2017; Hong et al., 2018b; Wibowo et al., 2018).

Another means of addressing this delay is to utilize the initial dip for fast fNIRS-BCI applications. The initial dip is an early change in oxygenation prior to any subsequent increase in CBF, which is spatially more specific to the site of neuronal activity (Vanzetta and Grinvald, 2008; Hong and Zafar, 2018). However, there is also a time lag in detecting the initial dip (Hong and Naseer, 2016). A previous study by Hong and Naseer (2016) showed that the initial dip could be detected using a vector phase diagram with a single threshold circle. The vector phase diagram is a computationally efficient method to detect both the initial dip and the HR by displaying the trajectories of ΔHbO and ΔHbR, as orthogonal components, in the ΔHbO-ΔHbR polar coordinates (Oka et al., 2015). It was further proposed to use q-step-ahead prediction scheme in combination with the vector phase diagram to reduce the time lag in detecting the initial dip. They showed that the initial dip could be detected in 0.9 s using the q-step-ahead prediction scheme, showing high potential for BCI applications. Later, Zafar and Hong (2017) attempted to find the features and temporal window size for classifying the initial dip duration in fNIRS signals of different mental tasks. They showed that the running temporal window size for fNIRS could be reduced from 5 to 2.5 s using initial dip features (i.e., signal mean and signal minimum) in the classification process. Li et al. (2017) also used the mean value of ΔHbO and ΔHbR signals in the 0–2 s window as an initial dip feature and achieved 85.5% classification accuracy for the classification of left- and right-hand movements. Similarly, Khan and Hong (2017) used signal minimum as an initial dip feature and achieved a classification accuracy of 75.6% in classifying four mental tasks in a reduced window size (i.e., 0–2 s).

The use of dual threshold circles in the vector phase diagram was proposed to improve the detection of both initial dip and the conventional HR (Zafar and Hong, 2018), see Figure 1. The threshold circles in the vector phase analysis helps to minimize the false detection of resting-state fluctuation and large fluctuations of ΔHbO and ΔHbR signals during the task period. The radius of the inner circle was set to the maximum HbO during the resting state, and the radius of the outer circle was set to the sum of the radius of the inner circle and 30% of the peak value of the main HR. The peak value of the HR was determined through the averaging over trials measured in the training phase. They showed that the use of dual threshold circles in the vector phase diagram resulted in an enhancement of the classification accuracies of two-finger tapping tasks. They also used the signal mean and the minimum signal value in 0–2.5 s time window to classify two-finger tapping tasks. However, windows of 0–2 s or 0–2.5 s are still too large for real-time BCI applications and hybridization of fNIRS with other rapid techniques such as EEG. Furthermore, the previously mentioned q-step-ahead prediction scheme by Hong and Naseer (2016) to reduce the delay was an offline analysis, and the validity of the predicted signals with multiple steps was not examined. Knowing the maximal prediction size of the q-step-ahead prediction method is important because the error of the predicted signals increases significantly with the increase of the number of step sizes. In addition, for real-time BCI applications, an online scheme is required to reduce the onset delay in fNIRS signals.

Concept of vector phase diagram with dual threshold circles (Hong and Zafar, 2018; Zafar and Hong, 2018).

In this study, the use of a kernel recursive least squares algorithm (KRLS) is proposed for the q-step-ahead prediction of fNIRS signals. Three most commonly used kernels (i.e., Gaussian, polynomial, and sigmoid) are tested to compare the errors in the predicted fNIRS signals. Then, the effectiveness of the proposed prediction scheme was evaluated using fNIRS signals of finger tapping tasks measured from the left motor cortex of eleven subjects. The results of the proposed scheme were compared with those of the commonly used recursive least squares (RLS) algorithm. This paper further presents the applicability of the q-step-ahead prediction scheme to reduce the time lag in detecting the initial dips in fNIRS signals.

 Methods Brain Activity Model and Kernel Recursive Least Square

In this paper, a brain activity is modeled in a linear form using the autoregressive moving average with exogenous signals (ARMAX) model as follows.

yi(k) =∑n = 1noaniyi(k−n)+∑m=1mobmi(m)u(k−m)            +∑p=1pocpi(p)wi(k−p)+co⋅1+εi(k)

where i represents the channel number; y is the measured ΔHbO/ΔHbR; u is the desired hemodynamic response function (dHRF); w is the physiological noise; ε is the zero-mean Gaussian noise; an, bm, cp, and co are unknown coefficients that are recursively estimated; and no, mo, and po are the orders of the system, input, and exogenous signals, respectively. For fNIRS, the exogenous signal w consists of specifically three sinusoidal signals representing cardiac, Mayer, and respiration related physiological noises (Abdelnour and Huppert, 2009; Nguyen H.-D. et al., 2018). Also, the exogenous signals can be dropped out in the estimation process (i.e., po = 0) if the prediction/tracking of the measured signal is focused. Nevertheless, the fNIRS signals were low- and high-pass filtered to minimize the effect of the physiological noises before the estimation process. Equation (1) can be written in a simplified vector form as follows.

yi(k)=φT(k)θi(k)+ei(k) φT(k)=[y(k−1)⋯y(k−no)                  u(k−1)⋯u(k−mo) w(k−1)⋯w(k−po)   1] θi(k)=[a1i⋯anoib1i⋯bmoic1i⋯cpo1co]

where φ(k) ∈ ℜ(n + m + p + 1) × 1 is the regression vector and superscript T stands for the transpose operator. Figure 2 shows the estimation/prediction scheme.

Online estimation/prediction scheme.

In this study, dHRF [i.e., u(k)] was generated by convolving the canonical HRF (cHRF), denoted by h(k), with a stimulation period, s(k), as follows.

u(k)=∑l=0k-1h(l)s(k-l), s(k)={1,  if k∈task,0,  if k∈rest,

where task and rest represent the task period and the rest period, respectively (task = 10 s and rest = 20 s in this study). cHRF was generated as a linear combination of three gamma functions by the following equation (Shan et al., 2014).

h(k)=∑j=13Ajkαj-1βjαje-βjkΓ(αj),

where j represents the number of gamma functions, Aj is the amplitude, αj and βj tune the shape and the scale, respectively, and k is the time step (in this work, A1 = −1.5, A2 = 7, A3 = −2, α1 = 1.5, α2 = 6, α3 = 16, and β1 = β2 = β3 = 1 were used). The unknown coefficients in Equation (2) are estimated and updated using the KRLS based on the optimization of the cost function given by

minθ(k)JKRLS=∑k=1NλN-k|y(k)-κ(φ(k), ·)Tθ(k-1)|2                    +RλN‖θ(k-1)‖H2, Φ(k)=[κ(φ(1), ·)  κ(φ(2), ·)                    ⋯   κ(φ(k), ·)]T

where κ represents the Mercer kernel, Φ is the kernel matrix of all k input data points, R is a positive number known as the regularization parameter, H represents the reproducing kernel Hilbert space (RKHS) associated with the Mercer kernel, and λ (0.98 in this study) is the forgetting factor. The performances of the following three most commonly used kernels in improving the prediction of the fNIRS signals are tested (Muller et al., 2001):

(i) Gaussian kernel

κ(φ,φ′)=exp(-||φ-φ′||22σ2)

where σ is a scaling factor, and φ′represents the new upcoming data.

(ii) Polynomial kernel

κ(φ,φ′)=(φTφ′+c)p

where c is a non-negative constant, and p is the order of the polynomial kernel.

(iii) Sigmoid kernel

κ(φ,φ′)=tanh(s(φTφ′)+t)

where s and t are suitable non-negative constants.

The basic idea is to map input data points to a high dimensional feature space (i.e., RKHS). This process allows the transformation of linear inner products into RKHS by simply changing their inner product into kernels (Schölkopf and Smola, 2002; Liu et al., 2010). The transformed feature space is then solved using the linear algorithm. The advantage of kernel-based algorithms is that they have a unique global solution that can be derived by solving a convex optimization problem (Chen et al., 2014). Furthermore, if data show a non-linear relationship, linear regression techniques cannot model them adequately. The kernel method can address this issue by moving to another feature space that is more likely to correspond to a linear model. However, the kernel method suffers from the overfitting problem because the Hilbert space induces high dimensionality of data. To address the issue of overfitting, the solution is penalized by limiting it to the L2 norm, as shown in Equation (8) (Evgeniou et al., 2000; Pillonetto et al., 2014), which is solved and updated as follows (Liu et al., 2010).

θ(k)=Φ(k)[Rλ+Φ(k)TΦ(k)]-1y(k),θ(k)=Φ(k)a(k),  a(k)=Q(k)y(k), Q(1)=[Rλ+κ(φ(1),φ(1))]-1,  a(1)=Q(1)y(1), K(k)=K(k−1)Tκ(φ(k), ·)=[κ(φ(1), φ(k)),     ⋯,                  κ(φ(k−1),φ(k))]T, z(k)=Q(k-1)K(k), δ(k)=Rλ+κ(φ(k),φ(k))-zT(k)K(k), Q(k)=δ-1(k)[Q(k-1)δ(k)+z(k)zT(k)-z(k)-z(k)1], e(k)=y(k)-KT(k)a(k-1), a(k)=[a(k-1)-z(k)δ-1(k)e(k)δ-1(k)e(k)].

As the kernel matrix grows linearly with the number of observations, the computational complexity of KRLS increases. The complexity is reduced by using the approximate linear dependency (ALD) criterion (Engel et al., 2004). The KRLS-ALD algorithm has been implemented using the MatlabTM kernel adaptive filtering toolbox (Van Vaerenbergh and Santamaría, 2013; Van Vaerenbergh, 2017).

Accordingly, using Equation (2), the estimated brain activity model can be represented as

y^i(k)=φT(k)θ^i(k)+ei(k).

For q-step-ahead prediction, Equation (21) can be written as follows.

y^i(k+q)=φT(k+q)θ^i(k)+ei(k).

The performance of the algorithm was tested using the percentage fitting (%FIT) criterion as follows (Pillonetto et al., 2014).

FIT=100(1-∑k=1N(yi(k)-y^i(k))2∑k=1N(yi(k)-mean(yi))2).

The performance criterion in (23) quantifies how much of the variance of y is captured by the q-step-ahead predicted signal (Pillonetto et al., 2014). Furthermore, %FIT criteria measure how accurately the q-step-ahead predicted signals are estimated.

 Experimental Data

Previously published experimental data (Zafar and Hong, 2018) of right-hand (thumb and little) finger tapping sessions from 11 subjects were used for validating the proposed q-step-ahead prediction scheme. Brain signals generated by the finger-tapping were acquired from the left motor cortex using the frequency domain fNIRS system (ISS Imagent, ISS Inc.) at a sampling rate of 9.19 Hz. The electrode placement and the corresponding emitter-detector distances are shown in Figure 3A. A total of 36 channels were formed using emitter-detector pairs.

Emitter-detector placement and experimental paradigm for the right-hand finger tapping task: (A) Emitter-detector placement and their distances, (B) experimental paradigm.

The experimental paradigm is shown in Figure 3B. The experimental paradigm consists of two sessions of finger tapping tasks. A session is composed of six trials of 30 s. Each trial includes a 10 s activity task followed by a 20 s rest. During the task period, the subjects were instructed to tap their right-hand finger as fast as they could without paying attention to the number of taps. The raw data (ΔHbO and ΔHbR) obtained from the ISS Imagent data acquisition and analysis software (ISS-Boxy) were pre-processed to remove physiological noises related to respiration, cardiac, and low-frequency drift signals. Fourth-order Butterworth low- and high-pass filters with cutoff frequencies of 0.15 and 0.01 Hz, respectively, were used to minimize the respiration, cardiac, and low-frequency drift signals from the obtained fNIRS signals.

 Detection of Initial Dip

The initial dip will be detected through the vector phase analysis with dual threshold circles (Yoshino and Kato, 2012; Hong and Naseer, 2016; Zafar and Hong, 2018), see Figure 1. Vector phase analysis is a polar coordinate plane method defined by ΔHbO and ΔHbR as orthogonal vector components. Two other vector components, cerebral oxygen exchange (ΔCOE) and cerebral blood volume (ΔCBV), are obtained by rotating the vector coordinate system by 45° counterclockwise using the following equations (Yoshino et al., 2013; Khan et al., 2018).

ΔCBV=12(ΔHbO+ΔHbR), ΔCOE=12(ΔHbR-ΔHbO).

The magnitude and the phase of a vector p = (ΔHbO, ΔHbR) in the phase plane can be calculated as follows.

|p|=ΔHbO2+ΔHbR2, ∠p=tan-1( ΔHbR ΔHbO)=tan-1( ΔCOE ΔCBV)+45o.

The degree of oxygen exchange is defined by the ratio of ΔCOE and ΔCBV. Therefore, the oxygen exchange in the blood vessel is represented by the change in ΔCOE. Using the abovementioned four indices, eight phases are defined on the vector phase diagram, see Figure 1. Phases 1–5 (i.e., Phase 1: 0 < ΔHbR < ΔHbO, ΔCOE < 0 < ΔCBV; Phase 2: 0 < ΔHbO < ΔHbR, 0 < ΔCOE < ΔCBV; Phase 3: ΔHbO < 0 < ΔHbR, 0 < ΔCBV < ΔCOE; Phase 4: ΔHbO < 0 < ΔHbR, ΔCBV < 0 < ΔCOE; Phase 5: ΔHbO < ΔHbR < 0, ΔCBV < 0 < ΔCOE) are defined as the initial dip phases because they reflect an increase in either ΔHbR or ΔCOE, whereas Phases 6 to 8 (i.e., Phase 6: ΔHbR < ΔHbO < 0, ΔCBV < ΔCOE < 0; Phase 7: ΔHbR < 0 < ΔHbO, ΔCOE < ΔCBV < 0; Phase 8: ΔHbR < 0 < ΔHbO, ΔCOE < 0 < ΔCBV) are defined as HR phases. If there are no threshold circles in the vector diagram, the resting-state fluctuation and large fluctuations of ΔHbO and ΔHbR signals during the task period with ΔCOE > 0 can easily be interpreted as an initial dip. Threshold circles (i.e., dual threshold circles) incorporated in the vector phase analysis help in minimizing the detection of false dips (Hong and Naseer, 2016; Zafar and Hong, 2018). The radius of the first (inner) threshold circle in Figure 1 was determined during the resting state period as follows (Hong and Naseer, 2016).

r1=max(ΔHbOresting2+ΔHbRresting2).

The single (inner) threshold circle can help to separate the resting-state fluctuation from the initial dip and task-related HR. However, a large fluctuation of ΔHbO and ΔHbR above the threshold circle can still falsely be interpreted as an initial dip. Therefore, a second (outer) threshold circle as an upper bound is drawn on the vector diagram to separate large ΔHbO and ΔHbR fluctuations from the initial dip. The radius for the second threshold circle was determined using the following equation (Zafar and Hong, 2018).

r2=r1+0.3(p1+SD),

where p1 and SD are the peak value and standard deviation of the averaged ΔHbO trial over several trials from the most active channel, where the most active channel means the channel that shows the largest difference between the maximum ΔHbO values in the resting state and the HR of the first trial during the training stage. The initial dips are detected if the trajectory lies in any phase from Phase 3 to Phase 5 and remains within the two threshold circles within first 2–4 s of the task period, and it moves to either Phases 7 or 8, after 2–4 s. The first (i.e., inner) threshold circle is used to detect the time instance of the occurrence of an initial dip and the HR. Any trajectory going outside the secondary threshold circle is considered as a false dip or noise.

 Results

The data during the resting-state and the first session were used in the training stage, whereas the second session data were used to test the proposed method. The parameters of the Gaussian (i.e., σ), polynomial (i.e., c and p), and sigmoid (i.e., s and t) kernels were determined iteratively, and the value with the maximum %FIT for each kernel was selected for further analysis. Using the data of 792 channels [i.e., 11 subjects × (36 HbO + 36 HbR)], the values of parameters σ, c, p, s, and t for the Gaussian, polynomial, and sigmoid kernels were found to be 1 in the training stage. For the regularization parameter (R), several different values were tested through trial and error, and R = 10−8 was found to achieve the best fitting (i.e., %FIT) of the predicted signals. R > 10−8 did not affect the %FIT of the predicted signals, but lower values decreased %FIT. Figure 4 shows the fitting of the one-step-ahead predicted ΔHbO and ΔHbR signals on top of the measured signals (ΔHbO, ΔHbR) using the RLS method and the Gaussian-kernel RLS method for both active (i.e., Ch. 18) and non-active (i.e., Ch. 3) channels of Subject 1, respectively.

Fitting of ΔHbO and ΔHbR signals for active (Ch. 18) and non-active (Ch. 3) channels at q = 1 for Sub. 1: (A,B,E,F) were obtained by using the RLS method; (C,D,G,H) were obtained by using the Gaussian-kernel RLS method.

Tables 1–4 reports the %FIT of ΔHbO and ΔHbR for individual subjects using RLS and KRLS. The statistical significance of the %FIT was verified using two-sample t-tests. Signal information in the predicted signals (i.e., %FIT) significantly decreases (p < 0.05) as the step size increases. Table 5 shows a comparison of the averaged %FITs of RLS and KRLS with the Gaussian, polynomial, and sigmoid kernels for different q-step-ahead predicted fNIRS (ΔHbO, ΔHbR) signals.

Averaged %FIT for ΔHbO and ΔHbR over all channels using RLS after training and testing.

Sub. %FIT
q = 1 (0.1 s) q = 5 (0.54 s) q = 10 (1.08 s) q = 15 (1.63 s) q = 20 (2.17 s)
ΔHbO ΔHbR ΔHbO ΔHbR ΔHbO ΔHbR ΔHbO ΔHbR ΔHbO ΔHbR
1 69.8 ± 8.6 68.3 ± 6.3 66.6 ± 9.3 64.9 ± 6.9 62.6 ± 10.0 60.8 ± 7.5 59.2 ± 10.6 57.2 ± 8.1 56.6 ± 10.9 54.2 ± 8.5
2 65.5 ± 9.3 63.7 ± 7.4 61.6 ± 10.1 59.7 ± 8.1 57.1 ± 11.2 54.9 ± 8.9 52.9 ± 12.1 50.5 ± 9.6 49.6 ± 12.8 46.9 ± 10.2
3 67.0 ± 6.6 65.7 ± 5.4 63.2 ± 7.2 61.6 ± 5.9 58.7 ± 7.8 56.9 ± 6.5 54.8 ± 8.4 52.8 ± 7.1 51.5 ± 8.9 49.4 ± 7.5
4 69.4 ± 9.8 68.7 ± 6.5 66.0 ± 10.4 65.3 ± 7.0 61.9 ± 11.2 61.3 ± 7.5 58.2 ± 11.8 57.7 ± 8.4 55.0 ± 12.3 54.6 ± 8.4
5 65.4 ± 7.3 64.4 ± 6.2 61.6 ± 7.8 60.6 ± 6.6 57.1 ± 8.4 55.9 ± 7.2 53.1 ± 8.9 51.7 ± 7.8 49.7 ± 9.3 48.2 ± 8.2
6 73.0 ± 9.4 69.4 ± 7.5 70.1 ± 10.4 66.1 ± 8.3 66.5 ± 11.6 62.1 ± 9.2 63.3 ± 12.6 58.6 ± 10.1 60.5 ± 13.4 55.6 ± 10.8
7 76.1 ± 6.1 75.3 ± 4.3 73.2 ± 6.4 72.5 ± 4.7 69.8 ± 6.8 69.2 ± 5.2 66.7 ± 7.3 66.2 ± 5.6 64.1 ± 7.7 63.6 ± 6.0
8 73.6 ± 8.4 67.9 ± 6.1 70.5 ± 8.9 64.4 ± 6.6 66.9 ± 9.6 60.2 ± 7.4 63.7 ± 10.2 56.5 ± 8.1 61.1 ± 10.6 53.5 ± 8.6
9 68.8 ± 6.9 64.9 ± 4.5 65.4 ± 7.4 61.2 ± 4.8 61.5 ± 8.1 56.7 ± 5.2 58.0 ± 8.5 52.8 ± 5.4 55.1 ± 8.9 49.6 ± 5.7
10 76.1 ± 9.5 71.3 ± 7.1 73.2 ± 10.2 67.9 ± 7.7 69.7 ± 11.0 63.9 ± 8.5 66.4 ± 11.6 60.3 ± 9.3 63.5 ± 11.9 57.2 ± 9.9
11 70.1 ± 8.1 66.8 ± 6.0 67.0 ± 8.8 63.2 ± 6.5 63.2 ± 9.5 58.9 ± 7.1 59.9 ± 10.1 55.2 ± 7.6 57.2 ± 10.4 52.2 ± 8.1
Mean ± SD 70.4 ± 8.2 67.9 ± 6.1 67.1 ± 8.8 64.3 ± 6.6 63.2 ± 9.6 60.1 ± 7.3 59.7 ± 10.2 56.3 ± 7.9 56.7 ± 10.9 53.2 ± 8.4

Averaged %FIT for ΔHbO and ΔHbR over all channels using KRLS with the Gaussian kernel after training and testing.

Sub. %FIT
q = 1 (0.1 s) q = 5 (0.54 s) q = 10 (1.08 s) q = 15 (1.63 s) q = 20 (2.17 s)
ΔHbO ΔHbR ΔHbO ΔHbR ΔHbO ΔHbR ΔHbO ΔHbR ΔHbO ΔHbR
1 94.4 ± 0.8 93.6 ± 1.2 92.7 ± 1.6 91.5 ± 1.8 90.4 ± 2.7 89.0 ± 2.8 88.1 ± 3.7 86.6 ± 3.5 85.9 ± 4.6 84.4 ± 4.2
2 94.6 ± 1.0 92.2 ± 4.1 92.0 ± 2.2 89.4 ± 4.2 88.7 ± 4.1 85.5 ± 4.6 85.5 ± 5.9 81.9 ± 5.6 82.7 ± 7.5 78.8 ± 6.7
3 94.3 ± 0.7 91.2 ± 4.0 91.9 ± 1.7 88.4 ± 3.8 88.5 ± 3.1 84.9 ± 4.1 85.3 ± 4.5 81.7 ± 4.6 82.4 ± 5.7 78.9 ± 5.3
4 94.7 ± 1.0 92.8 ± 1.8 92.9 ± 1.3 90.9 ± 2.2 90.4 ± 2.7 88.3 ± 3.0 87.8 ± 4.1 85.9 ± 3.8 85.4 ± 5.4 83.6 ± 4.5
5 94.4 ± 2.1 92.9 ± 3.2 92.3 ± 2.1 90.5 ± 3.1 89.5 ± 2.6 86.9 ± 3.3 86.7 ± 3.6 83.6 ± 4.0 84.3 ± 4.6 80.6 ± 4.6
6 94.8 ± 0.9 93.1 ± 2.1 93.3 ± 1.6 91.1 ± 2.6 91.2 ± 2.6 88.5 ± 3.3 88.9 ± 3.8 85.9 ± 4.1 86.6 ± 5.0 83.6 ± 4.9
7 95.1 ± 0.9 94.9 ± 0.8 93.2 ± 1.4 93.1 ± 1.2 90.7 ± 2.5 90.8 ± 2.0 88.3 ± 4.0 88.6 ± 2.7 86.1 ± 4.6 86.5 ± 3.5
8 94.6 ± 0.9 92.2 ± 2.4 92.2 ± 1.6 89.9 ± 2.8 89.5 ± 2.5 87.1 ± 3.6 87.1 ± 3.4 84.0 ± 5.4 84.9 ± 4.2 82.1 ± 5.0
9 94.3 ± 0.6 93.9 ± 0.5 92.4 ± 1.4 92.0 ± 1.1 89.8 ± 2.4 89.1 ± 2.2 87.3 ± 3.3 86.3 ± 3.4 85.0 ± 4.2 83.6 ± 4.5
10 95.4 ± 1.0 94.4 ± 1.4 93.7 ± 1.4 91.8 ± 1.4 91.5 ± 2.4 88.7 ± 2.2 89.3 ± 3.4 85.8 ± 3.5 87.5 ± 4.3 83.3 ± 4.7
11 94.4 ± 1.1 93.6 ± 1.1 92.8 ± 1.6 91.7 ± 1.3 90.7 ± 2.4 89.1 ± 2.1 88.5 ± 3.4 86.5 ± 3.1 86.6 ± 4.2 84.1 ± 4.1
Mean ± SD 94.6 ± 1.0 93.2 ± 2.1 92.7 ± 1.6 90.9 ± 2.3 90.1 ± 2.7 88.0 ± 3.0 87.5 ± 3.9 85.2 ± 4.0 85.2 ± 4.9 82.7 ± 4.7

Averaged %FIT for ΔHbO and ΔHbR over all channels using KRLS with the polynomial kernel after training and testing.

Sub. %FIT
q = 1 (0.1 s) q = 5(0.54 s) q = 10 (1.08 s) q = 15 (1.63 s) q = 20 (2.17 s)
ΔHbO ΔHbR ΔHbO ΔHbR ΔHbO ΔHbR ΔHbO ΔHbR ΔHbO ΔHbR
1 93.5 ± 2.4 92.6 ± 2.4 91.6 ± 3.4 90.6 ± 2.6 89.1 ± 5.1 88.1 ± 3.2 86.5 ± 6.6 85.7 ± 3.8 84.1 ± 7.9 83.6 ± 4.3
2 94.3 ± 0.9 87.5 ± 14.1 91.7 ± 2.1 84.7 ± 13.7 88.4 ± 4.0 81.0 ± 13.4 85.2 ± 5.7 77.5 ± 13.4 82.4 ± 7.4 74.4 ± 13.7
3 94.0 ± 1.2 83.1 ± 19.7 91.4 ± 1.9 80.4 ± 19.7 88.1 ± 3.2 77.1 ± 20.0 84.8 ± 4.6 73.8 ± 20.4 81.9 ± 5.8 71.1 ± 21.0
4 94.5 ± 1.1 89.5 ± 7.3 92.6 ± 1.5 87.6 ± 7.3 89.9 ± 2.7 85.1 ± 7.4 87.3 ± 4.0 82.5 ± 7.7 84.8 ± 5.2 80.3 ± 8.1
5 94.2 ± 1.8 90.5 ± 8.0 91.9 ± 2.2 87.8 ± 7.6 89.1 ± 2.3 84.2 ± 7.3 85.7 ± 4.5 80.8 ± 7.2 83.1 ± 4.9 78.0 ± 7.3
6 94.7 ± 0.8 91.8 ± 4.9 93.3 ± 1.5 89.7 ± 5.2 91.2 ± 2.5 87.1 ± 5.8 89.0 ± 3.7 84.4 ± 6.6 86.7 ± 4.8 82.1 ± 7.5
7 95.1 ± 0.9 94.1 ± 2.9 93.1 ± 1.4 92.2 ± 2.8 90.6 ± 2.4 89.9 ± 3.1 88.2 ± 3.5 87.6 ± 3.4 86.0 ± 4.5 85.5 ± 3.9
8 94.5 ± 0.9 91.2 ± 4.3 92.1 ± 1.6 88.9 ± 4.5 89.3 ± 2.5 86.0 ± 4.9 86.9 ± 3.4 83.2 ± 5.5 84.8 ± 4.2 80.8 ± 6.1
9 94.1 ± 0.9 93.6 ± 0.8 92.2 ± 1.5 91.5 ± 1.3 89.5 ± 2.4 88.5 ± 2.9 86.9 ± 3.4 85.5 ± 3.6 84.6 ± 4.3 82.8 ± 4.7
10 95.3 ± 0.9 93.8 ± 2.4 93.6 ± 1.4 91.2 ± 2.2 91.4 ± 2.4 88.1 ± 2.8 89.4 ± 3.4 85.2 ± 3.8 87.4 ± 4.8 82.6 ± 4.9
11 94.4 ± 0.8 93.3 ± 1.6 92.8 ± 1.4 91.3 ± 1.8 90.6 ± 2.3 88.7 ± 2.5 88.4 ± 3.1 86.1 ± 3.4 86.4 ± 3.8 83.8 ± 4.3
Mean ± SD 94.4 ± 1.1 91.0 ± 6.2 92.4 ± 1.8 88.7 ± 6.2 89.7 ± 2.9 85.8 ± 6.7 87.1 ± 4.2 82.9 ± 7.2 84.7 ± 5.2 80.5 ± 7.8

Averaged %FIT for ΔHbO and ΔHbR over all channels using KRLS with the sigmoid kernel after training and testing.

Sub. %FIT
q = 1 (0.1 s) q = 5 (0.54 s) q = 10 (1.08 s) q = 15 (1.63 s) q = 20 (2.17 s)
ΔHbO ΔHbR ΔHbO ΔHbR ΔHbO ΔHbR ΔHbO ΔHbR ΔHbO ΔHbR
1 86.3 ± 13.2 66.1 ± 24.5 84.1 ± 12.9 64.2 ± 24.2 81.5 ± 12.7 61.8 ± 23.9 78.9 ± 12.6 59.5 ± 23.7 76.7 ± 12.6 57.4 ± 23.6
2 87.5 ± 9.7 63.6 ± 26.2 84.5 ± 9.2 60.9 ± 25.1 80.8 ± 9.1 57.4 ± 23.8 77.3 ± 9.2 54.1 ± 22.8 74.2 ± 9.7 51.2 ± 22.1
3 84.0 ± 13.2 53.8 ± 27.3 81.0 ± 12.8 51.5 ± 26.6 77.2 ± 12.4 48.7 ± 25.8 73.7 ± 12.3 45.8 ± 25.1 70.7 ± 12.3 43.4 ± 24.5
4 86.8 ± 13.6 61.6 ± 27.1 84.6 ± 13.5 59.7 ± 26.8 81.7 ± 13.4 57.2 ± 26.6 78.9 ± 13.3 54.7 ± 26.4 76.3 ± 13.4 52.5 ± 26.3
5 89.1 ± 6.2 70.2 ± 19.1 86.4 ± 6.1 67.8 ± 18.6 82.9 ± 6.2 64.6 ± 18.0 79.7 ± 6.6 61.7 ± 27.5 76.8 ± 7.3 59.1 ± 17.3
6 87.9 ± 10.7 73.9 ± 22.9 86.2 ± 10.5 72 ± 22.6 83.8 ± 10.4 69.4 ± 22.2 81.4 ± 10.4 66.9 ± 21.9 79.1 ± 10.5 64.7 ± 21.7
7 94.0 ± 3.2 84.7 ± 16.7 91.8 ± 3.4 82.8 ± 16.4 89.2 ± 3.8 80.4 ± 16.1 86.6 ± 4.4 78.1 ± 15.8 84.2 ± 5.1 76.0 ± 15.7
8 84.3 ± 16.5 66.8 ± 20.1 82.1 ± 16.1 64.8 ± 19.7 79.4 ± 15.6 62.2 ± 19.2 77.1 ± 15.3 59.7 ± 18.8 74.9 ± 15.0 57.4 ± 18.5
9 91.1 ± 5.8 82.3 ± 13.4 88.7 ± 5.5 79.8 ± 13.0 85.8 ± 5.2 76.6 ± 12.6 83.1 ± 5.2 73.7 ± 12.3 80.8 ± 5.4 71.2 ± 12.2
10 90.3 ± 6.4 75.4 ± 18.4 88.4 ± 6.3 73.1 ± 17.8 86.1 ± 6.3 70.2 ± 17.2 83.8 ± 6.5 67.5 ± 17.0 81.7 ± 7.0 65.1 ± 16.9
11 91.3 ± 4.5 76.4 ± 22.7 89.4 ± 4.5 74.0 ± 22.2 86.9 ± 4.6 70.9 ± 21.5 84.5 ± 4.9 68.1 ± 21.1 82.4 ± 5.2 65.5 ± 20.9
Mean ± SD 88.4 ± 9.3 70.4 ± 21.6 86.1 ± 9.2 68.2 ± 21.2 83.2 ± 9.1 65.4 ± 20.6 80.4 ± 9.2 62.7 ± 21.1 77.9 ± 9.4 60.3 ± 19.9

Comparison of the averaged %Fits of RLS and KRLS (Gaussian, polynomial, sigmoid) for different q-step-ahead predicted ΔHbO and ΔHbR signals after training and testing.

q-step %Fit of fNIRS signals
ΔHbO ΔHbR
RLS KRLS with Gaussian KRLS with polynomial KRLS with sigmoid RLS KRLS with Gaussian KRLS with polynomial KRLS with sigmoid
1 (0.1 s) 70.4 ± 8.3 94.6 ± 1.0 94.4 ± 1.1 88.4 ± 9.3 67.9 ± 6.1 93.2 ± 2.1 91.0 ± 6.2 70.4 ± 21.6
5 (0.54 s) 67.1 ± 8.8 92.7 ± 1.6 92.4 ± 1.8 86.1 ± 9.2 64.3 ± 6.6 90.9 ± 2.3 88.7 ± 6.2 68.2 ± 21.2
10 (1.08 s) 63.2 ± 9.6 90.1 ± 2.7 89.7 ± 2.9 83.2 ± 9.1 60.1 ± 7.3 88.0 ± 3.0 85.8 ± 6.7 65.4 ± 20.6
15 (1.63 s) 59.7 ± 10.2 87.5 ± 3.9 87.1 ± 4.2 80.4 ± 9.2 56.3 ± 7.9 85.2 ± 4.0 82.9 ± 7.2 62.7 ± 21.1
20 (2.17 s) 56.7 ± 10.9 85.2 ± 4.9 84.7 ± 5.2 77.9 ± 9.4 53.2 ± 8.4 82.7 ± 4.7 80.5 ± 7.8 60.3 ± 19.9

A number of previous studies reported that the peak of the initial dip occurred at approximately 1.9–2.5 s (Hu and Yacoub, 2012; Zafar and Hong, 2017). Therefore, q = 15 (i.e., 1.63 s since the sampling frequency was 9.19 Hz in this study) was selected for further analysis. Table 1 shows that KRLS with the Gaussian kernel yielded the best fitting (p < 0.05) for the 1.63 s ahead predicted ΔHbO (i.e., 87.5%) and ΔHbR (i.e., 85.2%) signals as compared to all other methods. Therefore, the Gaussian-kernel RLS was further used for reducing the delay in detecting initial dips in the fNIRS signals. Figure 5 shows the 15-step-ahead predicted ΔHbO and ΔHbR signals of active channels for different subjects. It can be clearly seen that the predicted signals are well-ahead (i.e., blue dotted lines) and perfectly tracking the measured signals (solid red line).

Measured and 1.63 s ahead predicted signals (q = 15) of HbO (left) and HbR (right) with the Gaussian-kernel RLS algorithm: (A) Sub. 1 (Ch. 21), (B) Sub. 8 (Ch. 18), and (C) Sub. 10 (Ch. 30).

A comparison of vector-phase trajectories using measured and 1.63 s ahead predicted fNIRS signals for Subject 3 is shown in Figure 6. Table 6 shows the times of initial dip detection using 15-step-ahead predicted ΔHbO and ΔHbR signals for active channels of all subjects.

Comparison between measured and predicted (1.63 s) signals (Sub. 3, Ch. 21).

Time of initial dip detection using 1.63 s (q = 15) ahead prediction.

Sub. Channel Time of initial dip detection (sec)
Trial 1 Trial 2 Trial 3 Trial 4 Trial 5 Trial 6
1 21 0.65 0.11 0.11 0.11 Not detected Not detected
2 18 0.11 0.11 0.21 0.65 Not detected 0.65
3 21 0.11 0.11 Not detected 0.11 Not detected 0.11
4 29 0.11 0.11 0.11 0.32 0.11 0.21
5 33 0.11 0.32 0.21 Not detected 0.21 0.11
6 29 0.11 0.43 Not detected Not detected Not detected 0.11
7 17 0.11 Not detected 0.21 0.21 0.11 Not detected
8 21 Not detected 0.11 Not detected 0.32 0.11 Not detected
9 17 Not detected 0.21 0.43 0.21 0.21 Not detected
10 29 0.11 0.32 0.65 0.65 0.11 0.11
11 33 0.11 Not detected 0.11 0.11 Not detected 0.11
Discussions

The newly emerging neuroimaging modality (i.e., fNIRS) has a disadvantage of inherent onset delay from neuronal activation, which limits its application for rapid BCIs. To overcome this, the use of a kernel method for q-step-ahead prediction of fNIRS signals was proposed for the first time. The novelty of this study lies in using an online prediction scheme to reduce this onset delay for online applications.

A previous study by Hong and Naseer (2016) could reduce the delay in detecting initial dip in fNIRS signals to approximately 0.9 s using an offline q-step-ahead ARMAX model-based prediction scheme. Our results (Tables 1–5) reveal that the fitting accuracy of the q-step-ahead predicted signals decreased significantly (p < 0.05) with the increase of prediction step sizes. Therefore, the selection of a proper step size is very crucial to ensure that the predicted signals contain the maximum information of the measured signals.

In this study, a linear combination of three gamma functions was used (i.e., dHRF, see Figure 2). Most early studies used only s two-gamma-function dHRF to analyze the fNIRS time-series (Abdelnour and Huppert, 2009; Ye et al., 2009; Hu et al., 2010). A key drawback in using two gamma functions is that the initial dip duration is neglected in the estimation/prediction process. This limitation was overcome by using three gamma functions, which provides an extra degree of freedom by including the initial dip in the dHRF model for better estimation/prediction of the fNIRS signal.

The KRLS algorithm improves the fitting of the predicted signals as compared to the RLS algorithm by moving from the input space to the transformed feature space, i.e., a high dimensional space (see Tables 1–5). The non-linear relationship in the data cannot be adequately modeled by using linear regression techniques. The advantage of moving to a higher dimensional space is that there is a high probability that the data corresponds to a linear model, and it can be solved using the linear algorithms (Liu et al., 2010). Regarding the kernels, the Gaussian kernel yielded the best fitting of the predicted ΔHbO (i.e., 87.5%) and ΔHbR (i.e., 85.2%) signals at q = 15 step-ahead. The polynomial kernel also yielded good results for the ΔHbO signals, but the fitting slightly decreased for the ΔHbR signals. In contrast, the fitting of both ΔHbO and ΔHbR signals significantly decreased for the sigmoid kernel. Furthermore, the fitting of the predicted ΔHbR signals was lower than that of the predicted ΔHbO.

Early studies reported that the peak of initial dip occurred around 1.9–2.5 s (Malonek and Grinvald, 1996; Yacoub and Hu, 2001; Yacoub et al., 2001; Hu and Yacoub, 2012; Zafar and Hong, 2017). From this viewpoint, 1.63 s ahead prediction was selected in this study for an early detection of initial dips. Nevertheless, the peak of an initial dip depends on various factors, such as the type of task performed, the duration of the task period, and the brain area under investigation. The trajectories (Figure 6) for both measured and predicted signals were almost the same, showing that the predicted signals were well-tracking the measured signals. However, if the fitting of the predicted signal is not adequate, the trajectory can lead to a wrong decision regarding the detection of initial dip or HR. With 1.63 s ahead prediction (Table 6), the initial dips were detected in minimum 0.11 s (maximum 0.65 s), which is much lower than that of Hong and Naseer (2016) (i.e., 0.9 s). Furthermore, the initial dip phenomenon did not occur in some trials. In the literature, this issue has been discussed considering several issues. One interesting report is that it is due to the use of caffeine before the experiment (Behzadi and Liu, 2006; Hong and Zafar, 2018). In addition, the detection time of the initial dip varies among trials and subjects (Hu et al., 2013).

Finally, this study demonstrated a step moving toward the development of a real-time BCI and a brain monitoring system using fNIRS. The significance of this study lies in the fast detection of activity-related responses in fNIRS signals. Even if an initial dip is not present, the inherent onset delay in the conventional HR can be reduced using the proposed q-step-head prediction scheme. Moreover, the use of q-step-head prediction with improved fitting can help in the hybridization of fNIRS with other rapid modalities such as EEG. Nevertheless, further research is still required to improve the fitting of the predicted fNIRS signals with an accuracy more than 90% using advanced signal processing (Ghafoor et al., 2017; Chen et al., 2018; Hong et al., 2018a) and adaptive algorithms (Iqbal et al., 2018; Nguyen Q. C. et al., 2018). In the future, other types of kernels should also be investigated for further improvement of the predicted fNIRS signals. The limitations of this study are as follows: (i) the order of the system (an) and the input (bn) was set as 1 to ensure low computational complexity. Therefore, the optimal order of the system and the input for the prediction of fNIRS signals should be investigated further. (ii) Exogenous signals were excluded from the estimation/prediction process. These signals should be considered for further improvement of the predicted fNIRS signals in the future.

 Conclusion

In this study, the q-step-ahead prediction scheme based on KRLS was used to reduce the onset delay from the neuronal activation in fNIRS signals. fNIRS signals of right-hand finger tapping task acquired from the left motor cortex were used to evaluate the performance of the prediction scheme. The results show that the Gaussian kernel yields the best fitting for both ΔHbO (i.e., 87.5%) and ΔHbR (i.e., 85.2%) signals at q = 15 step ahead prediction (i.e., 1.63 s with the sampling frequency of 9.19 Hz). The application of the scheme was found to reduce the delay in detecting the initial dip. The improvement in the fitting of 1.63 s ahead predicted fNIRS signals enabled the detection of initial dip in 0.1 s. The reduction in the onset delay is a significant improvement in the development of real-time BCI applications using fNIRS.

 Data Availability Statement

The datasets analyzed in this article are not publicly available. Requests to access the datasets should be directed to kshong@pusan.ac.kr.

 Ethics Statement

The studies involving human participants were reviewed and approved by Pusan National University Institutional Review Board. The patients/participants provided their written informed consent to participate in this study.

 Author Contributions

AZ carried out the data processing and wrote the first draft of the manuscript. K-SH suggested the theoretical aspects of the current study, corrected the manuscript, and supervised the entire process from the beginning.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 References Abdelnour A. F. Huppert T. (2009). Real-time imaging of human brain function by near-infrared spectroscopy using an adaptive general linear model. Neuroimage 46, 133–143. 10.1016/j.neuroimage.2009.01.03319457389 Ahn S. Jun S. C. (2017). Multi-modal integration of EEG-fNIRS for brain-computer interfaces- current limitations and future directions. Front. Hum. Neurosci. 11:503. 10.3389/fnhum.2017.0050329093673 Behzadi Y. Liu T. T. (2006). Caffeine reduces the initial dip in the visual BOLD response at 3 T. Neuroimage 32, 9–15. 10.1016/j.neuroimage.2006.03.00516635577 Boas D. A. Elwell C. E. Ferrari M. Taga G. (2014). Twenty years of functional near-infrared spectroscopy: introduction for the special issue. Neuroimage 85(Pt 1), 1–5. 10.1016/j.neuroimage.2013.11.03324321364 Chen H. T. Jiang B. Lu N. Y. (2018). A multi-mode incipient sensor fault detection and diagnosis method for electrical traction systems. Int. J. Control Autom. Syst. 16, 1783–1793. 10.1007/s12555-017-0533-0 Chen T. Andersen M. S. Ljung L. Chiuso A. Pillonetto G. (2014). System identification via sparse multiple kernel-based regularization using sequential convex optimization techniques. IEEE Trans. Autom. Control 59, 2933–2945. 10.1109/TAC.2014.2351851 Cui X. Bray S. Reiss A. L. (2010). Speeded near infrared spectroscopy (NIRS) response detection. PLoS ONE 5:e15474. 10.1371/journal.pone.001547421085607 Engel Y. Mannor S. Meir R. (2004). The kernel recursive least-squares algorithm. IEEE Trans. Signal Process 52, 2275–2285. 10.1109/TSP.2004.830985 Ernst T. Hennig J. (1994). Observation of a fast-response in functional MR. Magn. Reson. Med. 32, 146–149. 10.1002/mrm.19103201228084231 Evgeniou T. Pontil M. Poggio T. (2000). Regularization networks and support vector machines. Adv. Comput. Math. 13:1. 10.1023/A:1018946025316 Ferrari M. Quaresima V. (2012). A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application. Neuroimage 63, 921–935. 10.1016/j.neuroimage.2012.03.04922510258 Frostig R. D. Lieke E. E. Tso D. Y. Grinvald A. (1990). Cortical functional architecture and local coupling between neuronal-activity and the microcirculation revealed by invivo high-resolution optical imaging of intrinsic signals. Proc. Natl. Acad. Sci. U.S.A. 87, 6082–6086. 10.1073/pnas.87.16.6082 Gateau T. Durantin G. Lancelot F. Scannella S. Dehais F. (2015). Real-time state estimation in a flight simulator using fNIRS. PLoS ONE 10:e0121279. 10.1371/journal.pone.012127925816347 Ghafoor U. Kim S. Hong K.-S. (2017). Selectivity and longevity of peripheral-nerve and machine interfaces: a review. Front. Neurorobotics 11:59. 10.3389/fnbot.2017.0005929163122 Ghafoor U. Lee J.-H. Hong K.-S. Park S.-S. Kim J. Yoo H.-R. (2019). Effects of acupuncture therapy on MCI patients using functional near-infrared spectroscopy. Front. Aging Neurosci. 11:237. 10.3389/fnagi.2019.0023731543811 Hong K.-S. Aziz N. Ghafoor U. (2018a). Motor-commands decoding using peripheral nerve signals: a review. J. Neural Eng. 15:031004. 10.1088/1741-2552/aab38329498358 Hong K.-S. Bhutta M. R. Liu X. L. Shin Y.-I. (2017). Classification of somatosensory cortex activities using fNIRS. Behav. Brain Res. 333, 225–234. 10.1016/j.bbr.2017.06.03428668280 Hong K.-S. Khan M. J. Hong M. J. (2018b). Feature extraction and classification methods for hybrid fNIRS-EEG brain-computer interfaces. Front. Hum. Neurosci. 12:246. 10.3389/fnhum.2018.0024630002623 Hong K.-S. Naseer N. (2016). Reduction of delay in detecting initial dips from functional near-infrared spectroscopy signals using vector-based phase analysis. Int. J. Neural Syst. 26:1650012. 10.1142/S012906571650012X26971785 Hong K.-S. Yaqub M. A. (2019). Application of functional near-infrared spectroscopy in the health industry: a review. J. Innov. Opt. Heal. Sci. 12:1930012. 10.1142/S179354581930012X Hong K.-S. Zafar A. (2018). Existence of initial dip for BCI: an illusion or reality. Front. Neurorobotics 12:69. 10.3389/fnbot.2018.00069 Hu X.-S. Hong K.-S. Ge S. S. (2013). Reduction of trial-to-trial variability in functional near-infrared spectroscopy signals by accounting for resting- state functional connectivity. J. Biomed. Opt. 18:017003. 10.1117/1.JBO.18.1.01700323291618 Hu X.-S. Hong K.-S. Ge S. S. Jeong M.-Y. (2010). Kalman estimator- and general liner model-based on-line brain activation mapping by near-infrared spectroscopy. Biomed. Eng. Online 9:82. 10.1186/1475-925X-9-8221138595 Hu X. P. Yacoub E. (2012). The story of the initial dip in fMRI. Neuroimage 62, 1103–1108. 10.1016/j.neuroimage.2012.03.00522426348 Hwang H. J. Choi H. Kim J. Y. Chang W. D. Kim D. W. Kim K. . (2016). Toward more intuitive brain-computer interfacing: classification of binary covert intentions using functional near-infrared spectroscopy. J. Biomed. Opt. 21:091303. 10.1117/1.JBO.21.9.09130327050535 Iqbal M. Rehan M. Hong K.-S. (2018). Robust adaptive synchronization of ring configured uncertain chaotic FitzHugh-Nagumo neurons under direction-dependent coupling. Front. Neurorobotics 12:6. 10.3389/fnbot.2018.0000629535622 Jasdzewski G. Strangman G. Wagner J. Kwong K. K. Poldrack R. A. Boas D. A. (2003). Differences in the hemodynamic response to event-related motor and visual paradigms as measured by near-infrared spectroscopy. Neuroimage 20, 479–488. 10.1016/S1053-8119(03)00311-214527608 Jiao Y. Zhang Y. Wang Y. Wang B. Jin J. Wang X. (2018). A novel multilayer correlation maximization model for improving CCA-based frequency recognition in SSVEP brain-computer interface. Int. J. Neural Syst. 28:1750039. 10.1142/S012906571750039328982285 Kato T. Kamei A. Takashima S. Ozaki T. (1993). Human visual cortical function during photic-stimulation monitoring by means of near-infrared spectroscopy. J. Cereb. Blood Flow Metab. 13, 516–520. 10.1038/jcbfm.1993.668478409 Khan M. J. Ghafoor U. Hong K.-S. (2018). Early detection of hemodynamic responses using EEG: a hybrid EEG-fNIRS study. Front. Hum. Neurosci. 12:479. 10.3389/fnhum.2018.0047930555313 Khan M. J. Hong K.-S. (2015). Passive BCI based on drowsiness detection: an fNIRS study. Biomed. Opt. Express 6, 4063–4078. 10.1364/BOE.6.00406326504654 Khan M. J. Hong K.-S. (2017). Hybrid EEG-fNIRS-based eight-command decoding for BCI: application to quadcopter control. Front. Neurorobotics 11:6. 10.3389/fnbot.2017.0000628261084 Khan M. J. Hong M. J. Hong K.-S. (2014). Decoding of four movement directions using hybrid NIRS-EEG brain-computer interface. Front. Hum. Neurosci. 8:244. 10.3389/fnhum.2014.0024424808844 Li R. H. Potter T. Huang W. T. Zhang Y. C. (2017). Enhancing performance of a hybrid EEG-fNIRS system using channel selection and early temporal features. Front. Hum. Neurosci. 11:462. 10.3389/fnhum.2017.0046228966581 Li W. Li M. Zhou H. Chen G. Jin J. Duan F. (2018). A dual stimuli approach combined with convolutional neural network to improve information transfer rate of event-related potential-based brain-computer interface. Int. J. Neural Syst. 28:1850034. 10.1142/S012906571850034X30185104 Liu W. Principe J. C. Haykin S. (2010). Kernel Adaptive Filtering: A Comprehensive Introduction. Hoboken, NJ: John Wiley & Sons, Inc. 10.1002/9780470608593 Liu X. Kim C.-S. Hong K.-S. (2018). An fNIRS-based investigation of visual merchandising displays for fashion stores. PLoS ONE 13:e0208843. 10.1371/journal.pone.020884330533055 Liu X. L. Hong K.-S. (2017). Detection of primary RGB colors projected on a screen using fNIRS. J. Innov. Opt. Heal. Sci. 10:1750006. 10.1142/S1793545817500067 Malonek D. Grinvald A. (1996). Interactions between electrical activity and cortical microcirculation revealed by imaging spectroscopy: implications for functional brain mapping. Science 272, 551–554. 10.1126/science.272.5261.5518614805 Matthews F. Pearlmutter B. A. Ward T. E. Soraghan C. Markham C. (2008). Hemodynamics for brain-computer interfaces. IEEE Signal Process. Mag. 25, 87–94. 10.1109/MSP.2008.4408445 Mcfarland D. J. Wolpaw J. R. (2010). Brain-computer interfaces for the operation of robotic and prosthetic devices. Adv. Comput. 79, 169–187. 10.1016/S0065-2458(10)79004-5 Mcfarland D. J. Wolpaw J. R. (2011). Brain-computer interfaces for communication and control. Commun. Acm 54, 60–66. 10.1145/1941487.194150621984822 Muller K.-L. Mika S. Ratsch S. Tsuda K. Scholkopf B. (2001). An introduction to kernel-based learning algorithms. IEEE Trans. Neural Netw. 12, 181–201. 10.1109/72.914517 Naseer N. Hong K.-S. (2015). fNIRS-based brain-computer interfaces: a review. Front. Hum. Neurosci. 9:3. 10.3389/fnhum.2015.00003 Naseer N. Noori F. M. Qureshi N. K. Hong K.-S. (2016). Determining optimal feature-combination for LDA classification of functional near-infrared spectroscopy signals in brain-computer interface application. Front. Hum. Neurosci. 10:237. 10.3389/fnhum.2016.0023727252637 Nguyen H.-D. Hong K.-S. (2016). Bundled-optode implementation for 3D imaging in functional near-infrared spectroscopy. Biomed. Opt. Express 7:3491–3507. 10.1364/BOE.7.00349127699115 Nguyen H.-D. Hong K.-S. Shin Y.-I. (2016). Bundled optode method in functional near-infrared spectroscopy. PLoS ONE 11:e0165146. 10.1371/journal.pone.016514627788178 Nguyen H.-D. Yoo S.-H. Bhutta M. R. Hong K.-S. (2018). Adaptive filtering of physiological noises in fNIRS data. Biomed. Eng. Online 17:180. 10.1186/s12938-018-0613-230514303 Nguyen Q. C. Piao M. Hong K.-S. (2018). Multivariable adaptive control of the rewinding process of a roll-to-roll system governed by hyperbolic partial differential equations. Int. J. Control Autom. Syst. 16, 2177–2186. 10.1007/s12555-017-0205-0 Nicolas-Alonso L. F. Gomez-Gil J. (2012). Brain computer interfaces, a review. Sensors 12, 1211–1279. 10.3390/s12020121122438708 Obrig H. (2014). NIRS in clinical neurology – a ‘promising’ tool. Neuroimage 85, 535–546. 10.1016/j.neuroimage.2013.03.04523558099 Ohi K. Shimada T. Kihara H. Yasuyama T. Sawai K. Matsuda Y. . (2017). Impact of familial loading on prefrontal activation in major psychiatric disorders:a near-infrared spectroscopy (NIRS) study. Sci. Rep. 7:44628. 10.1038/srep4426828295013 Oka N. Yoshino K. Yamamoto K. Takahashi H. Li S. Sugimachi T. . (2015). Greater activity in the frontal cortex on left curves: a vector-based fNIRS study of left and right curve driving. PLoS ONE 10:e0127594. 10.1371/journal.pone.012759425993263 Ortiz-Rosario A. Adeli H. (2013). Brain-computer interface technologies: from signal to action. Rev. Neurosci. 24, 537–552. 10.1515/revneuro-2013-003224077619 Pellicer A. Del Carmen Bravo M. (2011). Near-infrared spectroscopy: a methodology-focused review. Semin. Fetal Neonatal Med. 16, 42–49. 10.1016/j.siny.2010.05.00320580625 Pillonetto G. Dinuzzo F. Chen T. De Nicolao G. Ljung L. (2014). Kernel methods in system identification, machine learning and function estimation: a survey. Automatica 50, 657–682. 10.1016/j.automatica.2014.01.001 Power S. D. Kushki A. Chau T. (2011). Towards a system-paced near-infrared spectroscopy brain-computer interface: differentiating prefrontal activity due to mental arithmetic and mental singing from the no-control state. J. Neural Eng. 8:066004. 10.1088/1741-2560/8/6/06600421975364 Schölkopf B. Smola A. J. (2002). Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond. London: The MIT Press. Schudlo L. C. Chau T. (2014). Dynamic topographical pattern classification of multichannel prefrontal NIRS signals: II. Online differentiation of mental arithmetic and rest. J. Neural Eng. 11:016003. 10.1088/1741-2560/11/1/01600324311057 Schudlo L. C. Chau T. (2018). Development of a ternary near-infrared spectroscopy brain-computer interface: online classification of verbal fluency task, stroop task and rest. Int. J. Neural Syst. 28:1750052. 10.1142/S012906571750052629281922 Shan Z. Y. Wright M. J. Thompson P. M. Mcmahon K. L. Blokland G. G. de Zubicaray G. I. . (2014). Modeling of the hemodynamic responses in block design fMRI studies. J. Cereb. Blood Flow Metab. 34, 316–324. 10.1038/jcbfm.2013.20024252847 Shin J. Im C. H. (2018). Performance prediction for a near-infrared spectroscopy-brain-computer interface using resting-state functional connectivity of the prefrontal cortex. Int. J. Neural Syst. 28:1850023. 10.1142/S012906571850023529914312 Shin J. Kwon J. Choi J. Im C.-H. (2017). Performance enhancement of a brain-computer interface using high-density multi-distance NIRS. Sci. Rep. 7:16545. 10.1038/s41598-017-16639-029185494 Van Vaerenbergh S. (2017). Kernel Adaptive Filtering Toolbox. Matlab. Available online at: https://www.mathworks.com/matlabcentral/fileexchange/46747-kernel-adaptive-filtering-toolbox (accessed September 12, 2018). Van Vaerenbergh S. Santamaría I. (2013). A comparative study of kernel adaptive filtering algorithms, in IEEE Digital Signal Processing and Signal Processing Education Meeting (DSP/SPE) (Napa, CA), 181–186. 10.1109/DSP-SPE.2013.6642587 Vanzetta I. Grinvald A. (2008). Coupling between neuronal activity and microcirculation: implications for functional brain imaging. HFSP J. 2, 79–98. 10.2976/1.288961819404475 Villringer A. Planck J. Hock C. Schleinkofer L. Dirnagl U. (1993). Near-infrared spectroscopy (NIRS) - a new tool to study hemodynamic-changes during activation of brain-function in human adults. Neurosci. Lett. 154, 101–104. 10.1016/0304-3940(93)90181-J8361619 Watanabe H. Shitara Y. Aoki Y. Inoue T. Tsuchida S. Takahashi N. . (2017). Hemoglobin phase of oxygenation and deoxygenation in early brain development measured using fNIRS. Proc. Natl. Acad. Sci. U.S.A. 114, E1737–E1744. 10.1073/pnas.161686611428196885 Wibowo S. A. Lee H. Kim E. K. Kim S. (2018). Collaborative learning based on convolutional features and correlation filter for visual tracking. Int. J. Control Autom. Syst. 16, 335–349. 10.1007/s12555-017-0062-x Yacoub E. Hu X. (2001). Detection of the early decrease in fMRI signal in the motor area. Magn. Reson. Med. 45, 184–190. 10.1002/1522-2594(200102)45:2<184::AID-MRM1024>3.0.CO;2-C11180423 Yacoub E. Shmuel A. Pfeuffer J. Van De Moortele P. F. Adriany G. Ugurbil K. . (2001). Investigation of the initial dip in fMRI at 7 Tesla. NMR Biomed. 14, 408–412. 10.1002/nbm.71511746933 Yang C. Han X. Wang Y. Saab R. Gao S. Gao X. (2018). A dynamic window recognition algorithm for SSVEP-based brain-computer interfaces using a spatio-temporal equalizer. Int. J. Neural Syst. 28:1850028. 10.1142/S012906571850028430105920 Yang D. Hong K.-S. Yoo S.-H. Kim C.-S. (2019). Evaluation of neural degeneration biomarkers in the prefrontal cortex for early identification of patients with mild cognitive impairment: an fNIRS study. Front. Hum. Neurosci. 13:317. 10.3389/fnhum.2019.0031731551741 Yaqub M. A. Woo S.-W. Hong K.-S. (2018). Effects of HD-tDCS on resting-state functional connectivity in the prefrontal cortex: an fNIRS study. Complexity 2018:1613402. 10.1155/2018/1613402 Yazdani M. Salarieh H. Foumani M. S. (2018). Bio-inspired decentralized architecture for walking of a 5-link biped robot with compliant knee joints. Int. J. Control Autom. Syst. 16, 2935–2947. 10.1007/s12555-017-0578-0 Ye J. C. Tak S. Jang K. E. Jung J. W. Jang J. D. (2009). NIRS-SPM: statistical parametric mapping for near-infrared spectroscopy. Neuroimage 44, 428–447. 10.1016/j.neuroimage.2008.08.03618848897 Yi G. Mao J. X. Wang Y. N. Guo S. Y. Miao Z. Q. (2018). Adaptive tracking control of nonholonomic mobile manipulators using recurrent neural networks. Int. J. Control Autom. Syst. 16, 1390–1403. 10.1007/s12555-017-0309-6 Yoshino K. Kato T. (2012). Vector-based phase classification of initial dips during word listening using near-infrared spectroscopy. Neuroreport 23, 947–951. 10.1097/WNR.0b013e328359833b22989928 Yoshino K. Oka N. Yamamoto K. Takahashi H. Kato T. (2013). Correlation of prefrontal cortical activation with changing vehicle speeds in actual driving: a vector-based functional near-infrared spectroscopy study. Front. Hum. Neurosci. 7:895. 10.3389/fnhum.2013.0089524399953 Zafar A. Hong K.-S. (2017). Detection and classification of three-class initial dips from prefrontal cortex. Biomed. Opt. Express 8, 367–383. 10.1364/BOE.8.00036728101424 Zafar A. Hong K.-S. (2018). Neuronal activation detection using vector phase analysis with dual threshold circles: a functional near-infrared spectroscopy study. Int. J. Neural Syst. 28:1850031. 10.1142/S012906571850031430045647

Funding. This work was supported by the National Research Foundation (NRF) of Korea under the auspices of the Ministry of Science and ICT, Republic of Korea (Grant No. NRF-2017R1A2A1A17069430).