There is currently no standard method of topographic NIRS data analysis for brain mapping. In NIRS detection of hemodynamic responses, the light attenuation measured by the equipment needs to be converted to HbO and HbR concentration changes via the modified Beer-Lambert law (MBLL) [10]. Hence, a differential path length factor (DPF, intra- and inter-subject varying) in the MBLL has to be assumed to account for the increase of the path length between a source and a detector [11].

The classical approach in topographic NIRS data analysis is a paired *t*-test to determine if a concentration change between two states (for instance, "rest" vs. "task") is statistically significant. Many researchers nowadays use this approach [12–14], because it is simple and, thus, can provide a quick assessment to the task. One of the most popular tools in this regard is a Matlab-based program known as HomER ([15]
http://www.nmr.mgh.harvard.edu/PMI/).

However, there are limitations to the classical *t*-test. First, a maximum activation period needs to be predefined. The remaining temporal information not included in the defined activation period therefore is ignored, leading to underestimation of brain activation. Another problem is the DPF assumption. Since the DPF is intra- and inter-subject variant and impossible to be measured for every measurement location with commonly obtainable continuous wave NIRS equipment, use of a constant DPF leads to biased estimation of concentration changes [16].

To overcome such problems of the classical *t*-test, a number of research groups have used various general linear model (GLM)-based methods for analysis of NIRS data [16–20]. The GLM-based methods were initially developed for fMRI-based functional brain mapping [21]. The GLM is a statistical linear concept that explains measured data in the form of a linear combination of several explanatory variables plus an error term. The explanatory variables, modelled according to the time course, separately account for the brain-activity-evoked signals and noises. Therefore, the estimation of brain activity is reduced to the problem of estimating the relevant coefficients with proper statistics.

The GLM-based methods negate the need for user-defined rest and task periods, because the response is modelled according to the entire time course. The temporal information over the entire time course, thus, is examined. On the other hand, the GLM-based methods investigate the temporal variation pattern of the signal, and estimate the coefficients with statistics at different measurement locations, separately. Therefore, these methods are robust in cases where an assumed constant DPF is used.

GLM-based methods, however, cannot offer on-line analysis. This constrains their use in applications where real-time information or feedback is required. Real-time brain imaging data analysis in comparison with off-line methods, significantly, improves the information acquisition rate and the feedback speed. Furthermore, it may potentially benefit the development of brain-computer interfaces (BCI).

The critical task in modifying a GLM-based method on-line is the recursive estimation of GLM coefficients. Previous researchers have used different methods to achieve recursive estimation, including recursive least square [22], Cholesky-decomposition-based recursive least square [23], and Kalman filtering [24, 25]. All these approaches can effectively and recursively estimate GLM coefficients.

It is not sufficient to draw an updated brain activation map only by estimating GLM coefficients on-line. It is known that NIRS time series contain noises from different sources. Very-low-frequency noises caused by optodes shifts or slow cardiac/vascular artifacts [26], for example, might lead to biased estimation. The temporal correlation caused by physiological (cardiac, respiratory, blood pressure) noises might lead to an inflated *t*-value and, thereby, overestimation of brain activation. Furthermore, for functional brain mapping applications, a statistical test is very important, since it will provide significance verification of the derived estimation. All these issues with regard to on-line versions of GLM-based methods need to be addressed.

On-line estimation of GLM coefficients is generally discussed in [24]. In [25], the whole framework provided on-line estimation of GLM coefficients, and a relevant statistical test analysed the fMRI data. However, this work did not consider temporal correlation in fMRI data, which might also exist in NIRS data. In [26], a framework for NIRS-based BCI applications was developed that can estimate GLM coefficients without statistical information in classifying different hand tasks in real-time. In [27], the feasibility of estimating GLM coefficients on-line using a Kalman filter by studying a set of fMRI data was examined.

In the present study, we develop an on-line Kalman-estimator- and GLM-based NIRS data processing framework for brain activation mapping. We aim to answer two questions. (i) Is it possible to covert the off-line GLM-based method to an on-line version for NIRS-based brain activation mapping? The proposed method can provide updated brain activation maps on-line as well as track task-related brain-active areas from an early stage while the experiment is running. (ii) Can the proposed method prevent noises that distort the estimation while data is sequentially incorporated?

### NIRS measurement system and experimental procedure

Five right-handed healthy volunteers (all male, aged 24 to 34 years) participated in the experiment. None of the participants had a history of any neurological disorder. All of the participants provided written informed consent. The experiment was conducted in accordance with the latest Declaration of Helsinki. The data were acquired with a continuous-wave NIRS imaging system (DYNOT: DYnamic Near-infrared Optical Tomography) obtained from NIRx Medical Technologies, Brooklyn, NY, at a sampling rate of 1.81 Hz. The system emits laser lights of different wavelengths (760 nm and 830 nm) from each source. Figure 1 shows the channel distribution and measurement location. The distance between different optodes is 2 mm.

In the experiment, the subjects were asked to perform a finger-tapping task. The experiment consisted of a 42 sec preparation period and 10 sessions. Each session included a 21 sec finger-tapping period and a 30 sec rest period. Accordingly, the total duration of the experiment was 552 sec.

### Analysis framework of NIRS data

A schematic summarization of our method is as follows. (i) The relative concentration changes of two blood chromospheres, HbO and HbR, are calculated via MBLL [10, 27]; (ii) A linear model is built according to the experimental procedure to fit the relative HbO concentration change. The model describes both the signals corresponding to the brain activity and the noise; (iii) The model coefficients are recursively estimated using the Kalman estimator; (iv) At every time step, the brain-activity-related coefficient is selected, and then passed to the *t*-statistical test to determine if its value is statistically greater than zero (a value greater than zero indicates brain activity): In this way, a probability brain activation map is drawn. Figure 2 is a schematic flow chart of the framework.

### NIRS measurement model

In NIRS measurement, the optical density variation (Δ*OD*) can be expressed as a linear combination of hemoglobin concentration changes (Δ*C*
_{
HbO
} and Δ*C*
_{
HbR
} ) multiplied by proper coefficients. Their relationship is described in MBLL terms as

\Delta O{D}^{\lambda ,i}={\alpha}^{\lambda}\Delta C{L}^{\lambda ,i}DP{F}^{\lambda}

(1)

where {\alpha}^{\lambda}\Delta C={\alpha}_{HbO}^{\lambda}\Delta {C}_{HbO}+{\alpha}_{HbR}^{\lambda}\Delta {C}_{HbR}, *λ* is the wavelength of the laser source, *i* indicates channel number, *α*
_{
HbO
} [μM^{-1}mm^{-1}] and *α*
_{
HbR
} [μM^{-1}mm^{-1}] are the extinction coefficients of the HbO and HbR, *L* is the distance between the source and the detector, and DPF is the differential path length factor. In the present study, the optical density variation was derived by dividing the light intensity measured at each time step by the light intensity measured at the first time step.

### General linear model (GLM)

In NIRS-based studies, both HbO and HbR concentration changes can reflect changes in the rCBF. However, it has been suggested that HbO is a more sensitive indicator of such changes [28]. Therefore, only the HbO concentration change data was considered in the present study.

The GLM design process is described in detail in [29]; we provide only a brief description here. A design matrix *H* including a set of explanatory variables is predefined in order to model the observed NIRS time series. Five explanatory variables are considered. The first variable models the HbO concentration changes (the brain activity signals) using a stimulus vector convolved with the basis function (BF, a double-gamma model; [30]). The second one models the baseline level, and the remaining variables represent a set of high-pass filters (discrete cosine transform, DCT) [30] with a cut-off frequency of 0.0006 Hz.

At time *k*, *y*^{i} (*k*), the measured NIRS data of channel *i*, is predicted by *H*(*k*), specifically by multiplying a coefficient vector plus an error term *ε*^{i} (*k*). The model can then be expressed as

{y}^{i}(k)=H(k){\beta}^{i}(k)+{\epsilon}^{i}(k)

(2)

where *β*^{i} (*k*) is the coefficient vector quantifying the magnitude of the explanatory variables. In vector *β*^{i} (*k*), we are interested in the component *β*^{i}
_{1}(*k*), which reflects the magnitude of the task-evoked brain response: By statistically determining if it is greater than zero, the existence of brain activity at the area covered by channel *i* can be confirmed (it is greater than zero) or ruled out (it is less than zero).

Several physiological processes are known to produce temporal correlation in NIRS data, which might lead to inflated *t*-values, and thus underestimation of brain activity. One way to deal with this problem is to make the model fit an AR (*p*) model (an autoregressive model of the order *p*). This leads to a decomposition of the error term *ε* into a systematic and a model conform error part. After this, the AR transformation coefficient is applied to both sides of the regression equation

{y}^{i}(k)\text{-}\rho {y}^{i}(k\text{-}1)=H(k)\text{-}\rho H(k\text{-}1){\beta}^{i}(k)+{u}^{i}(k)

(3)

where *ρ* is the estimated autocorrelation coefficient in an AR(1) process, and *u*(*k*) = *ρε*(*k-1*) + *ε*(*k*). By redefining each transformed variable, that is, *y*^{i}
***(*k*) = *y*^{i} (*k*) - *ρy*^{i} (*k-* 1), *H**(*k*) = *H*(*k*) - *ρH*(*k-* 1), one can simplify Equation (3) to

{y}^{i*}(k)={H}^{*}(k){\beta}^{i}(k)+{u}^{i}(k)

(4)

It is worth noting that we make an assumption, |*β*(*k*) - *β*(*k-* 1)| < *ζ*, for the AR(1) model used on-line in the current study, where *ζ* > 0 is an arbitrary small number. This assumption compromises the model's robustness. However, as the result shows, the temporal correlation can be effectively reduced. In the present study, *β*^{i} (*k*) was updated with a Kalman estimator.

### Kalman estimator

The Kalman filtering method is a recursive tracking scheme that estimates the state of a process using an updated regularized linear inversion routine [31]. After decades of development, the Kalman filtering is very mature. Due to its remarkable estimation performance, the Kalman filtering is widely used in many areas [32–34] including neuroscience [24, 25, 35, 36]. In the present study, the Kalman filter was used as a model coefficients estimator. The model coefficients from all of the 24 channels were updated in parallel. For a given channel, the state vector, transition equation and observation equation can be described in the form

X(k)={[\begin{array}{cccc}{\beta}_{1}(k)& {\beta}_{2}(k)& ...& {\beta}_{L}(k)\end{array}]}^{T}

(5)

X(k)=A\widehat{X}(k-1)+w(k)

(6)

where *L* is the number of explanatory variables. The state is assumed to follow a random walk with zero drift over time: Thus, *A* equals the identity matrix, and the process noise *w*(*k*) ~ *N*(0, *Q*), *y*(*k*) is the measured data, *H*(*k*) is the vector of explanatory variables, and the observation noise *v*(*k*) ~ *N*(0, *R*). The filter performs state estimation by the iterative process

{\widehat{X}}^{-}(k)=A\widehat{X}(k-1)

(8)

{P}^{-}(k)=AP(k-1){A}^{T}+Q

(9)

K(k)={P}^{-}(k){H}^{T}(k){E}^{-1}(k)

(10)

\widehat{X}(k)={\widehat{X}}^{-}(k)+K(k)\Delta y(k)

(11)

P(k)=(I-K(k)H(k)){P}^{-}(k),

(12)

where *E*(*k*) = *H*(*k*)*P*^{-}(*k*)*H*^{T} +*R*, \Delta y(k)=y(k)-H(k){X}^{-}(k), *K*(*k*) is the Kalman gain, and *P*(*k*) is the updated error covariance matrix. In this notation, the superscript (-) refers to the intermediate state and covariance predictions provided by the state update model, which are then modified by the measured data to produce the next state value. In the present study, the state vector was initialized to zero. The a priori estimates of the process and observation noise covariances (*Q* and *R* respectively) were (1%/sec)^{2} and (0.5 μM/sec)^{2}, according to a restricted maximum likelihood (ReML) estimation and an empirical-experimental performance check based on a set of training data. We collected the training data during 3 sessions of finger tapping for each subject. We estimated the *Q* and *R* values in two steps: They were estimated separately from each of the subjects by ReML, averaged, and then adjusted according to the performance in practically estimation based on the training data.

###
*t*-Statistics

The estimated model coefficient vector *β* was used to calculate a relevant *t*-value for a one-tailed *t*-test to test the null hypothesis *c*^{T} *β* = 0 [29]. In the present study, the *t*-statistics of channel *i* at time step *k* were obtained using

{t}^{i}(k)=\frac{{c}^{T}{\beta}^{i}(k)}{\sqrt{{\widehat{\sigma}}^{i2}(k){c}^{T}{[{\displaystyle \sum _{\text{1}}^{k}{H}_{k}^{T}{H}_{k}}]}^{-\text{1}}c}},

(13)

where *c* is a vector of contrast for selecting the coefficient of interest [29], and {\widehat{\sigma}}^{i2} is the residual sum-of-squares divided by the appropriate degrees of freedom, and is given by

{\widehat{\sigma}}^{i2}(k)=\frac{1}{k-L}{\displaystyle \sum _{1}^{k}{[{y}_{k}^{i}-{H}_{k}{\beta}_{k}^{i}]}^{2}},

(14)

where *L* is the number of regressors. Therefore, the null hypothesis *c*^{T}
*β*
_{
i
} (*k*) = 0 was assessed by comparing *t*^{i} (*k*) with a *t*-distribution with *k-L* degrees of freedom. By setting proper *p*-values with Bonferroni correction, a statistical activation map of the detected area could be displayed.