Optimal hemodynamic response model for functional near-infrared spectroscopy

Functional near-infrared spectroscopy (fNIRS) is an emerging non-invasive brain imaging technique and measures brain activities by means of near-infrared light of 650–950 nm wavelengths. The cortical hemodynamic response (HR) differs in attributes at different brain regions and on repetition of trials, even if the experimental paradigm is kept exactly the same. Therefore, an HR model that can estimate such variations in the response is the objective of this research. The canonical hemodynamic response function (cHRF) is modeled by two Gamma functions with six unknown parameters (four of them to model the shape and other two to scale and baseline respectively). The HRF model is supposed to be a linear combination of HRF, baseline, and physiological noises (amplitudes and frequencies of physiological noises are supposed to be unknown). An objective function is developed as a square of the residuals with constraints on 12 free parameters. The formulated problem is solved by using an iterative optimization algorithm to estimate the unknown parameters in the model. Inter-subject variations in HRF and physiological noises have been estimated for better cortical functional maps. The accuracy of the algorithm has been verified using 10 real and 15 simulated data sets. Ten healthy subjects participated in the experiment and their HRF for finger-tapping tasks have been estimated and analyzed. The statistical significance of the estimated activity strength parameters has been verified by employing statistical analysis (i.e., t-value > t_critical and p-value < 0.05).

Introduction

Functional near-infrared spectroscopy (fNIRS) is a non-invasive and an emerging neuro-imaging technique (Santosa et al., 2013, 2014). The brain functional information is decoded through the interpretation of the variation in the optical properties of near-infrared (NIR) light (Naseer et al., 2014). NIRS monitors regional cerebral blood flow (rCBF) variations through the absorption changes of the NIR light at wavelengths between 650–950 nm. The oxy-hemoglobin (HbO) and deoxy-hemoglobin (HbR) are two major chromospheres in the blood which absorb NIR light (Cope and Delpy, 1988). The concentration of HbO and HbR varies in the capillary blood during the rest and task sessions (Hu et al., 2013). Thus, brain functional information can be revealed by the estimation of HbO and HbR. fNIRS, with the ability to estimate both chromospheres, is a potential brain imaging modality (Kamran and Hong, 2014). Functional-magnetic resonance imaging (fMRI) (Cohen et al., 2014; Zhou et al., 2014) and electro-encephalography (EEG) are most frequently used cortical imaging modalities in past. In comparison, the huge size of fMRI and its paramagnetic constraints and low spatial resolution of EEG (Soekadar et al., 2014) enhance the potential position of fNIRS with good temporal resolution applicable to brain-computer interface (BCI) applications (Hu et al., 2010). Its spatial resolution is affected by the low penetration depth but lies in between fMRI and EEG (Boudriay et al., 2014). In addition to these attributes, its portability, safety and low cost features makes it on top position for rehabilitation and BCI applications (Khan et al., 2014).

The detection of neuronal activation in a particular cortical area refers to find out a specific waveform in the hemodynamic response (HR) (Ciftçi et al., 2008). In hemodynamic related neuro-imaging modality, likewise fNIRS, the existence of such waveform is indicated as a statistical comparison to a specific time series shape, known as canonical hemodynamic response function (cHRF) (Abdelnour and Huppert, 2009). The cHRF has a key role in the analysis of fNIRS time series as its shape varies between different brain regions, repetition of trials, and among subjects as well (Hong and Nugyen, 2014). The difference in the dynamic shape of HRF during event-related motor and visual paradigms revealed that the peak times of HbO, HbR, and total hemoglobin (HbT) for visual paradigm are approximately equal unlike for motor paradigm (Jasdzewski et al., 2003). Additionally, it is found that the wavelength dependent differential path length factor (DPF) and age can also affect the characteristics of HR (Duncan et al., 1996). A mismatch between these features could result as a decrease in the detection performance (Ciftçi et al., 2008).

The most commonly used model for cHRF is composed of two gamma functions to characterize the shape and undershoot, respectively. It has been frequently be implemented in the analysis of fMRI temporal data. The performance accuracy of detection is improved by modifying the basis set, incorporating temporal derivative (TD) and dispersion derivation (DD) along with blood oxygen level dependent (BOLD) in the design matrix (Friston et al., 1998). Thus, modeled HR is represented as a linear combination of three waveforms. The characteristics of BOLD response are similar to cHRF used in NIRS data analysis. But fNIRS signal has additional challenge of temporal correlation present in the optical signal caused by physiological noises (Hu et al., 2010). The feature values for basis set were imposed constraints (time-to-peak, number of positive and negative peaks, time to- and magnitude of undershoot) to improve the extraction of the specific wave-pattern that formulate the dynamic shape of cHRF (Ciftçi et al., 2008). Finally, a general linear model (GLM) framework is being utilized to tune the unknown parameters in the model using Bayesian approach (Ciftçi et al., 2008). A new public statistical toolbox (NIRS-SPM freely available at http://bispl.weebly.com/nirs-spm.html#/) for the analysis of NIRS data was introduced in 2009, incorporating GLM based estimation of the cortical activity (Ye et al., 2009). NIRS-SPM is an extension of statistical parameter mapping toolbox for fMRI, thus it uses the GLM approach with basis set (Friston et al., 1998) to map the neuronal activities on brain templates. A detailed comparison of modeling techniques for HRF in fMRI regarding assumptions in the models, the complexity in their design and interpretation shows that it is difficult to accurately recover true task-evoked changes (Lindquist et al., 2009). In their study, the gradient approach has been utilized for the estimation of the free parameters that define the shapes of different HRF models. The fNIRS time series is contaminated with physiological noises. Thus, addition of physiological signals in the design matrix could improve the detection of task-related HR and its application to BCI (Abdelnour and Huppert, 2009). The parameters of cHRF were assumed to be fixed in Abdelnour and Huppert (2009) and activity strength parameters have been estimated using Kalman filters. The conventional averaging techniques have been used most frequently in past but its major drawback is the number of trials necessary to derive the stable HRF (Scarpa et al., 2010). Scarpa et al. (2013) proposed the methodology of near/closed channels to remove the physiological noises with fix parameters to extract cHRF. Several studies in past have presented the idea of combining HR model and adaptive signal processing algorithms to recursively tune the model parameters (Kamran and Hong, 2013; Santosa et al., 2013; Hong and Nugyen, 2014). Thus, an optimal HR model is still a topic of interest for many researchers in fNIRS area.

In this paper, an optimal HR model has been proposed for the analysis of fNIRS time series. The measured HR is modeled as a linear combination of evoked-HR, the physiological noises (cardiac pulsation, respiratory beat and low frequency Mayer waves) and base-line correction. The evoked-HR is the convolution of cHRF and the experimental paradigm. The cHRF has been modeled as a linear combination of two gamma functions (Lindquist et al., 2009). Six parameters in the cHRF model have been supposed as free parameters (delay of response relative to onset, delay of undershoot, dispersion of response, dispersion of undershoot, baseline, and a scaling factor) (Friston et al., 1994). The selection of optimal parameters in cHRF is the crucial step due to the variability of HR in different brain regions, repetition of trials and among subjects as well. In addition, the variation in the frequency and amplitudes of the physiological noises is a common phenomenon in the optical signal (Abdelnour and Huppert, 2009). Thus, these parameters in physiological noises are also supposed as free. The optimal parameters whose subsequent results best fit to the measured HR are found through an iterative optimization process. Initial values for these free parameters have been used from existing literature (Friston et al., 1994; Abdelnour and Huppert, 2009). Finally, the brain-activation model is formulated as an optimization problem with 12 free parameters. The formulated pre-optimized model is passed to an iterative simplex method with initial parameter vector. The simplex algorithm (Spendley et al., 1962) and its modified version (Nelder and Mead, 1965) has frequently been used in past for many signal processing and engineering-design optimization applications (Luersen and Riche, 2004). Fifteen simulated data sets have been generated with known parameters to verify the correctness of the proposed algorithm. Simulated data sets have been generated through method described in Prince et al. (2003). A low error in the estimation of free parameters shows great potential of the proposed algorithm in this field. Ten healthy participants have been examined for motor related, typically box-car based rest-task-rest experiment. Finally, brain functional maps have been shown to localize the cortical activity.

Materials and Methods

Data Acquisition and Pre-Processing

Three different types of optical neuro-imaging systems are available commercially, namely, continuous wave (CW), frequency domain (FD), and time-resolved spectroscopy (Hu et al., 2012; Schudlo et al., 2013). CW is the least expensive and most frequently used approach for BCI applications. It provides the relative change in the concentration of HbO and HbR. The CW-NIRS imaging system (DYNOT: Dynamic Near-infrared Optical Tomography; NIRx Medical Technologies, Brooklyn, NY) was used in this study with two wavelengths of NIR light (760 and 830 nm). It has 32 optodes which can be configured as emitters or detectors according to the experimental requirement with data acquisition frequency of 1.81 Hz. The optodes were placed on the left motor cortex at 16 different emitter–detector–pair locations. The source–detector separation was approximately 3 cm. The optode configuration has been shown in Figure 1A.

FIGURE 1

Figure 1. Source/detector locations and distribution of channels (A) and experimental paradigm (B).

The optical density variations measured through NIRS imaging system is converted into relative concentration changes of HbO and HbR using modified-Beer Lambert law (MBLL). According to MBLL, the optical densities at two different wavelengths can be solved by simple algebra to estimate the relative concentration changes of HbO and HbR with assumption of constant scattering (Power et al., 2011; Kamran and Hong, 2013)

$\begin{array}{ll}\Delta Hb{O}^i(k)=\frac{({\epsilon }_{HbR}^{{\lambda }_1}\frac{\Delta O{D}^{{\lambda }_2}(k)}{DP{F}^{{\lambda }_2}})-({\epsilon }_{HbR}^{\lambda 2}\frac{\Delta O{D}^{{\lambda }_1}(k)}{DP{F}^{{\lambda }_1}})}{l^i({\epsilon }_{HbR}^{{\lambda }_1}{\epsilon }_{HbO}^{{\lambda }_2}-{\epsilon }_{HbR}^{{\lambda }_2}{\epsilon }_{HbO}^{{\lambda }_1})},\hfill & \left(1\right)\hfill \end{array}$

$\begin{array}{ll}\Delta Hb{R}^i(k)=\frac{({\epsilon }_{HbO}^{{\lambda }_2}\frac{\Delta O{D}^{{\lambda }_1}(k)}{DP{F}^{{\lambda }_1}})-({\epsilon }_{HbO}^{{\lambda }_1}\frac{\Delta O{D}^{{\lambda }_2}(k)}{DP{F}^{{\lambda }_2}})}{l^i({\epsilon }_{HbR}^{{\lambda }_1}{\epsilon }_{HbO}^{{\lambda }_2}-{\epsilon }_{HbR}^{{\lambda }_2}{\epsilon }_{HbO}^{{\lambda }_1})},\hfill & \left(2\right)\hfill \end{array}$

where ΔHbOi (k) and ΔHbRi (k) are relative concentration changes of HbO and HbR, respectively, k is the discrete time, i represents the ith-channel of emitter-detector pair, λ1 and λ2 represent 760 and 830 nm wavelengths, ελ1HbO, ελ1HbR, ελ2HbO, and ελ2HbR indicates the extinction coefficients (refers to the measure of absorption of light) of HbO and HbR at two different wavelengths, respectively, ΔODλj (k) is the optical density variation at kth-sample time at particular wavelength (j = 1, 2), li is the source-detector separation and DPFλj is the DPF at particular wavelength (j = 1, 2). The extinction coefficients corresponding to 760 nm are 1.486 (for HbO) and 3.843 (for HbR) and those corresponding to 830 nm are 2.231 (for HbO) and 1.791 (for HbR) (Kamran and Hong, 2014).

Experimental Setup and Paradigm

Ten right-handed healthy subjects (age: 28 ± 7 years) participated in this study. None of the subject has the neuronal-disorder history before the experimentation. The written consent of each participant was collected before experimentation. The experiment was in accordance with the latest version of the Declaration of Helsinki. The subjects were completely introduced about the experiment and instrument before the start of the experimentation. The subjects were advised to avoid the head motion as much as possible. The load of NIRS optode fibers were supported through the hanger available with the instrument. The experiment includes a typical box-car rest-task-rest session. The experiment includes an initial rest of 10 s followed by a task session of 10 s and 30 s of a rest at the end. The subjects were instructed to tap their right index finger during task. A monitor screen was placed in front of the subject at a distance of approximately 110 cm. It remained blank during rest sessions and showed “finger tapping” during task session. The experimental paradigm has been shown in Figure 1B.

Simulated Subjects Data

The proposed algorithm is based upon an iterative optimization algorithm. Thus, it is necessary to verify the algorithm through simulated data sets with known values of free parameters. The simulated data is generated and supposed to be combination of HRF, three physiological signals, baseline term and random Gaussian noise. Fifteen different simulated data sets were generated using different values of free parameters. The data set has been generated with the methods described in existing literature (Prince et al., 2003; Abdelnour and Huppert, 2009)

$\begin{array}{ll}HRF(k)=h{(k)}^{\ast }u(k),\hfill & \left(3\right)\hfill \end{array}$

$\begin{array}{ll}\begin{array}{l}h(k)=\left[\frac{k^{{\alpha }_1-1}{\beta }_1^{{\alpha }_1}{e}^{-{\beta }_{1}k}}{\Gamma ({\alpha }_1)}-\frac{k^{{\alpha }_2-1}{\beta }_2^{{\alpha }_2}{e}^{-{\beta }_{2}k}}{6\Gamma ({\alpha }_2)}\right]\\ \end{array}\hfill & \left(4\right)\hfill \end{array}$

where u is the experimental paradigm, h is the cHRF, α₁ is the delay of the response, α₂ is the delay of the undershoot, β₁ is the dispersion of the response, β₂ is the dispersion of the undershoot and Γ represents the Gamma distribution. The physiological signals in simulated data have been generated through the linear combination of three sinusoids (Abdelnour and Huppert, 2009; Kamran and Hong, 2014). The specific values of free parameters used for all 15 data sets have been listed in Table 1.

TABLE 1

Table 1. The results of 15 simulated data-sets: the actual values of parameters (A) and estimated ones through the proposed algorithm (E).

Linear Brain Model and Parameter Optimization

GLM is a statistical linear model to decompose the output into predefined regressors. The existence of a particular regressor of interest depends upon the intensity of the activity strength parameter (Santosa et al., 2013). The positive value of the activity strength parameter with increasing t-value shows the significant existence of the particular regressor in the measured waveform (Hu et al., 2010). In this study, the measured HbO concentration change is supposed to be the linear combination of evoked-HRF, physiological noises and base-line correction

$\begin{array}{ll}\begin{array}{l}{y}_{HbO}^i(k)=a_o+a_{1}HRF(k)+a_c\sin (2\pi f_{c}k)+a_r\sin (2\pi f_{r}k)\\ +a_m\sin (2\pi f_{m}k)+{\epsilon }^i(k),\hfill \end{array}\hfill & \left(5\right)\hfill \end{array}$

where yⁱ_HbO is the measured HbO time series at ith-channel, a_o is the baseline, a₁ is the activity strength parameter, a_c, a_r, a_m, f_c, f_r, f_m are the amplitudes and the frequencies of the cardiac, respiratory and Mayer wave respectively and εⁱ (k) is the zero mean Gaussian noise at kth-sample time. Let us define a cost function J as sum of squares of residuals

$\begin{array}{ll}\begin{array}{l}J={\displaystyle \sum _{k=1}^N\{}{y}_{HbO}^{i}i(k)-(a_o+a_{1}HRF(k)+a_c\sin (2\pi f_{c}k)+\\ a_r\sin (2\pi f_{r}k)+a_m\sin (2\pi f_{m}k)){\}}^2.\hfill \end{array}\hfill & \left(6\right)\hfill \end{array}$

The above cost function can be formulated in an optimization environment with constraints

$\begin{array}{ll}\begin{array}{l}\min J({\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,a_o,a_1,a_c,a_m,a_r,f_c,f_r,f_m)s.t\hfill \\ C_1:2\le {\alpha }_1\le 10,C_7:0\le a_c\le 2,\hfill \\ C_2:6\le {\alpha }_2\le 20,C_8:0\le a_r\le 2,\hfill \\ C_3:0.5\le {\beta }_1\le 2,C_9:0\le a_m\le 2,\hfill \\ C_4:0\le {\beta }_2\le 1.5,C_{10}:0.5\le f_c\le 1.5,\hfill \\ C_5:0\le a_o\le 20,C_{11}:0.2\le f_r\le 0.3,\hfill \\ C_6:0\le a_1\le 15.C_{12}:0.09\le f_m\le \mathrm{0.1.}\hfill \end{array}\hfill & \left(7\right)\hfill \end{array}$

The optimal values of free parameters (α*₁,α*₂,β*₁,β*₂,a*_o,a*₁,a*_c,a*_m,a*_r,f*_c,f*_r,f*_m) are estimated by improved version of simplex method [later named as Nelder–Mead simplex method (NMSM)]. The iteration of NMSM can be performed by three steps, namely, ordering, centroid and transformation. The simplex of size a is defined at initial point (Haftka et al., 1990)

$\begin{array}{ll}{x}_j=x_o+p{e}_j+{\displaystyle \sum _{\begin{array}{l}k=1\hfill \\ k\ne j\hfill \end{array}}^{n}q{e}_k};j=1,2,\dots ,n,\hfill & \left(8\right)\hfill \end{array}$

$\begin{array}{ll}p=\frac{a}{n\sqrt{2}}(\sqrt{n+1}+n-1)\text{ \& }q=\frac{a}{n\sqrt{2}}(\sqrt{n+1}-1).\hfill & \left(9\right)\hfill \end{array}$

where x_j (j = 1, 2,.., n) represent the vertices, x_o is the initial guess, n = 12 in this study and represents the number of free parameters, e_j represents the unit vector in the direction of jth vertex. The next step is to order the function in increasing order at all vertices and it is easy to sort as

$\begin{array}{ll}J\text{(}{x}_l\text{) < }J\text{(}{x}_s\text{) < }J\text{(}{x}_h\text{)}\text{.}\hfill & \left(10\right)\hfill \end{array}$

where x_l, x_h, and x_s, are the vertices with minimum value, maximum value and second highest value of the cost function, respectively. The next step is to discard the highest value by defining the centroid

$\begin{array}{ll}\overline{x}=\frac{1}{n}{\displaystyle \sum _{\begin{array}{l}i=0\hfill \\ i\ne h\hfill \end{array}}^n{x}_i},\hfill & \left(11\right)\hfill \end{array}$

where x is the centroid. The replacement of upper bound vertex of the cost function is done by reflection, expansion, contraction, and shrinkage (Lagarias et al., 1998). The mathematical equations for all these steps are given below

$\begin{array}{ll}\text{Reflection}:x_r=\overline{x}+{\delta }_1(x_h-\overline{x}),\hfill & \left(12\right)\hfill \end{array}$

$\begin{array}{ll}\text{Expansion}:x_e=\overline{x}+{\delta }_2(x_r-\overline{x}),\hfill & \left(13\right)\hfill \end{array}$

$\begin{array}{ll}\text{Contraction}:x_c=\overline{x}+{\delta }_3(x_h-\overline{x}),\hfill & \left(14\right)\hfill \end{array}$

$\begin{array}{ll}\text{Shrinkage}:x_e=\overline{x}+{\delta }_4(x_l-x_i);i=0,1,\dots ,n,\hfill & \left(15\right)\hfill \end{array}$

where δ₁, δ₂, δ₃, and δ₄ are coefficients of reflection, expansion, contraction, and shrinkage, respectively. The typical values of these coefficients have been chosen as 1, 2, 0.5, and 0.5, respectively (Lagarias et al., 1998; Luersen and Riche, 2004). The schematic of the algorithm is shown in Figure 2. The updated value at any step of iteration shall be replaced with bounded value, if it crosses bound at any step. The detail about the algorithm can be found in Haftka et al. (1990), Lagarias et al. (1998), and Luersen and Riche (2004).

FIGURE 2

Figure 2. Schematic of Nelder–Mead simplex method.

Functional Brain Maps and Statistical Significance

The estimation of the cortical activation and its localization is a challenging task in the analysis of fNIRS data series. Previous studies showed that localization of the cortical activation could be statistically estimated by fitting the estimated HRF to a pre-defined HRF (Hu et al., 2010; Kamran and Hong, 2014; Santosa et al., 2014). Let the optimal brain activation model is

$\begin{array}{ll}\begin{array}{l}{y}_{HbO}^i(k)=a_0^{\ast }+a_1^{\ast }HR{F}^{\ast }(k)+a_c^{\ast }\sin (2\pi f_c^{\ast }k)+a_r^{\ast }\sin (2\pi f_r^{\ast }k)\\ +a_m^{\ast }\sin (2\pi f_m^{\ast }k)+{\epsilon }^i(k),\hfill \end{array}\hfill & \left(16\right)\hfill \end{array}$

The estimated optimal value of the activity strength parameter, a₁, related to HRF indicates the activation of the particular brain region with proper statistics (Hu et al., 2010). The basic idea is to test whether the estimated value of the activity strength parameter is greater or less than a target value zero with statistically significance (t-value > t_critical and p-value < 0.05). Thus, it is equivalent of testing a null hypothesis H_o with proper statistics i.e.,

$\begin{array}{ll}{H}_o:a_1^{\ast }=0\hfill & \left(17\right)\hfill \end{array}$

$\begin{array}{ll}{t}_{value}=\frac{a_1^{\ast }-0}{SE(a_1^{\ast })}.\hfill & \left(18\right)\hfill \end{array}$

where SE is the standard error of the estimated coefficient.

Results

In this study an online recursive optimization algorithm is proposed for cortical activation detection to display discrete brain maps. The algorithm is verified through 15 synthetic data sets and implemented to real data sets of 10 healthy subjects. The cHRF (Figure 3, top right panel) is modeled as two Gamma functions (Figure 3, top left panel). All the simulated data sets have been displayed in Figure 3 (bottom panel). It is obvious to note that different width, height and undershoot have been considered for verifications. Table 1 summarized the values of free parameters and their estimate through proposed algorithm. The comparison of HRF with actual parameter values and estimated parameter values has been shown in Figure 4. The results of the estimated-evoked-HR of 10 subjects have been presented in Figure 5. The significance of results has been verified using t-test. The t-maps of the cortical activations have been presented in Figure 6. Figures 7, 8 display a comparison of the estimated parameters related to most active channel in each subject of real data set and simulated data set, respectively.

FIGURE 3

Figure 3. Hemodynamic response function generations: two Gamma functions for generation of cHRF (top left), the standard cHRF (top right) and different simulated HRF (bottom).

FIGURE 4

Figure 4. HRF using actual values of free parameter (solid) and HRF using estimated values of free parameter (circular blue).

FIGURE 5

Figure 5. Results of estimated HRF related to most active channel corresponding to all subjects.

FIGURE 6

Figure 6. t-Maps of each subject and all channels.

FIGURE 7

Figure 7. Variations in the estimated parameters in real data sets.

FIGURE 8

Figure 8. Variations in the estimated parameters in simulated data sets.

Discussion

The non-invasive neuro-imaging techniques have a favorable position due to an increasing demand of BCI applications in the rehabilitation and medical diagnostics. There are several studies reported for the estimation of HRF in fMRI with numerical optimization techniques (Lindquist et al., 2009; Shah et al., 2014). But in the case of fNIRS, it constitutes an additional challenge of the physiological noise in the optical signal. Recently, several studies have been reported to analyze fNIRS time series using existing/new and/or modified versions of existing HRF models (Abdelnour and Huppert, 2009; Hu et al., 2010; Kamran and Hong, 2013, 2014; Santosa et al., 2013; Scarpa et al., 2013; Hong and Nugyen, 2014). The approaches vary in their implementation from simple estimation algorithms to more complex adaptive algorithms (Kamran and Hong, 2013) and blind signal processing (Santosa et al., 2013).

Hu et al. (2010) decomposed measured HRF into predefined regressors (evoked-HRF, base line correction and three others were included to design a set of high pass filter). Santosa et al. (2013) implemented the independent component analysis (ICA) framework to extract the statistically significance of a known wave pattern in the observed fNIRS data. Kamran and Hong (2013) explored the idea of adaptive signal processing to tune the variations in the measured HRF using parameter varying methodology. Later Kamran and Hong (2014) proposed to decompose HbO signal using ARMAX model for better cortical estimation as compared to existing ones. Abdelnour and Huppert (2009) proposed the adaptive framework to tune the HR with pre-built HRF in the model. Scarpa et al. (2013) emphases to remove the physiological noises by incorporating a near-detector (<0.7 cm from source) and to model the remaining signal as linear combination of pre-set evoked-HR and base line in Bayesian framework for adaptive tuning.

The DPF is a wavelength dependent factor and also varies with age of the subject causing variations in the hemodynamic signal (Duncan et al., 1996). Jasdzewski et al. (2003) reported that difference found in the characteristics of HRF in multiple brain regions. Their results suggest that the characteristics of HRF in different brain regions show variations. It is also observed in their study that some of the features like initial dip could be found in certain brain regions but not all. Hong and Nugyen (2014) developed a state-space model for different brain regions using adaptive signal processing framework. Their results revealed that there exist a significant difference between responses of different brain regions. It is a well-known fact that the hemodynamic signal has inter-subject variability as well as inter-trial variations. Hu et al. (2013) analyzed the reduction of trial-to-trial variations by analyzing correlation in the observed signal of different channels.

Thus existing literature suggest that the HR varies in its shape and characteristics not only in different brain areas, but it differs corresponding to the different mental task complexity, repetition, inter-subject and inter-trials as well. Some fNIRS-based BCI studies suggest that learning can improve the response, that is, less effort is required to repeat the same mental task. Thus, it is very important to model the HR in an adaptive framework together with a setup in which the parameters of HRF could be optimized as per real-time information in the measured response. Therefore, a recursive optimization algorithm have been presented in this study to model the variations in the HR. In contrast to existing fNIRS data analysis models, the proposed model has the capability to track time-varying characteristics (if exist) of HRF within same experiment as well. The estimation of pre-defined parameters of simulated data is shown in Table 1. It is obvious to observe that the proposed model estimated these parameters with a significant accuracy. After the validation of the algorithm using simulated data set, the methodology is applied to real data set of 10 subjects. The optimized values of HRF model parameters for 10 subjects have been listed in Table 2. The t-maps of each subject and all channels have been presented in Figure 6. It is evident from Figures 5, 6 that inter-subject difference exist in HRF parameters. Generally, inter-subject variability is due to the individual's differences in anatomical factors likewise skull and cerebrospinal fluid (CSF) structure, vessels distributions and the ratios of the arteries and veins. Thus, some of the subjects have more activation as compared to others.

TABLE 2

Table 2. The values of free parameter estimated through proposed algorithm in most active channel of each subject.

Conclusion

An optimal HR model has been proposed that can extract the shape and scale of the HRF in addition to the amplitudes and the frequencies of the physiological sinusoids. Twelve parameters in the HR model have been supposed free with bounded constraints. The problem is formulated as an optimization problem and solved through an iterative optimization framework. The algorithm is first verified through different simulated data sets with known values of free parameters. A low error in estimation shows the accuracy of the proposed methodology. Furthermore, the algorithm is implemented to real data sets of 10 healthy participants. The parameters in HRF while repeating same trials are found to be different. Thus, it shall be beneficial for fNIRS data analysis, as the proposed model can track the characteristics of changes in HRF and physiological noises.

Conflict of Interest Statement

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.

Acknowledgments

This work was supported by a 2-Year Research Grant of Pusan National University.

References

Abdelnour, A. F., and Huppert, T. (2009). Real-time imaging of human brain function by near-infrared spectroscopy using an adaptive general linear model. Neuroimage 46, 133–143. doi: 10.1016/j.neuroimage.2009.01.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Boudriay, Y., Feltane, A., and Besio, W. (2014). Significant improvement in one-dimensional cursor control using Laplacian electroencephalography over electroencephalography. J. Neural Eng. 11:035014. doi: 10.1088/1741-2560/11/3/035014

PubMed Abstract | CrossRef Full Text | Google Scholar

Ciftçi, K., Sankur, B., Kahya, Y. P., and Akin, A. (2008). Constraining the general linear model for sensible hemodynamic response function waveforms. Med. Biol. Eng. Comput. 46, 779–787. doi: 10.1007/s11517-008-0347-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, O., Koppel, M., Malach, R., and Friedman, D. (2014). Controlling an avatar by thought using real-time fMRI. J. Neural Eng. 11:035006. doi: 10.1088/1741-2560/11/3/035006

PubMed Abstract | CrossRef Full Text | Google Scholar

Cope, M., and Delpy, D. T. (1988). System for long-term measurement of cerebral blood and tissue oxygenation on newborn-infants by near-infrared trans-illumination. Med. Biol. Eng. Comput. 26, 289–294. doi: 10.1007/BF02447083

PubMed Abstract | CrossRef Full Text | Google Scholar

Duncan, A., Meek, J. H., Clemence, M., Elwell, C. E., Fallon, P., Tyszczuk, L., et al. (1996). Measurement of cranial optical path length as a function of age using phase resolved near infrared spectroscopy. Pediatr. Res. 39, 889–894.

PubMed Abstract | Google Scholar

Friston, K. J., Fletcher, P., Joseph, O., Holmes, A., Rugg, M. D., and Tunner, R. (1998). Event-related fMRI: characterizing differential responces NeuroImage 7, 30–40.

PubMed Abstract | Google Scholar

Friston, K. J., Holmes, A. P., Worsley, J.-P., Frith, C. D., and Frackowiak, R. S. J. (1994). Statistical parameter maps in functional imaging: a general linear model approach. Hum. Brain Mapp. 2, 189–210. doi: 10.1002/hbm.460020402

CrossRef Full Text | Google Scholar

Haftka, R. T., Gurdal, Z., and Kamat, M. P. (1990). Element of Structural Optimization. Dordrecht: Kluwer Academic Publisher. doi: 10.1007/978-94-015-7862-2

CrossRef Full Text

Hong, K.-S., and Nugyen, H.-D. (2014). State-space models of impulse hemodynamic responses over motor, somatosensory, and visual cortices. Biomed. Opt. Express 5, 1778–1798. doi: 10.1364/BOE.5.001778

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, X. S., Hong, K.-S., and Ge, S. S. (2012). fNIRS-based online deception decoding. J. Neural Eng. 9:026012. doi: 10.1088/1741-2560/9/2/026012

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, X. S., Hong, K.-S., and 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. doi: 10.1117/1.JBO.18.1.017003

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, X. S., Hong, K.-S., Ge, S. S., and Jeong, M. Y. (2010). Kalman estimator-and general linear model-based on-line brain activation mapping by near-infrared spectroscopy. BioMed. Eng. Online 9:82. doi: 10.1186/1475-925x-9-82

PubMed Abstract | CrossRef Full Text | Google Scholar

Jasdzewski, G., Starngman, G., Wanger, J., Kwong, K. K., Poldrack, R. A., and 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. doi: 10.1016/S1053-8119(03)00311-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Kamran, M. A., and Hong, K.-S. (2013). Linear parameter-varying model and adaptive filtering technique for detecting neuronal activities: an fNIRS study. J. Neural Eng. 10:056002. doi: 10.1088/1741-2560/10/5/056002

PubMed Abstract | CrossRef Full Text | Google Scholar

Kamran, M. A., and Hong, K.-S. (2014). Reduction of physiological effects in fNIRS waveforms for efficient brain-state decoding. Neurosci. Lett. 580, 130–136. doi: 10.1016/j.neulet.2014.07.058

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, M. J., Hong, M. J., and Hong, K.-S. (2014). Decoding of four movement directions using hybrid NIRS-EEG brain-computer interface. Front. Hum. Neurosci. 8:244. doi: 10.3389/fnhum.2014.00244

CrossRef Full Text | Google Scholar

Lagarias, J. C., Reeds, J. A., Wright, M. H., and Wright, P. E. (1998). Convergence properties of Nelder-Mead simplex method in low dimension. SIAM J. Optim. 9, 112–147. doi: 10.1137/S1052623496303470

CrossRef Full Text | Google Scholar

Lindquist, A., Meng, L. J., Atlas, L. Y., and Wager, T. D. (2009). Modeling the hemodynamic response function in fMRI: efficiency, bias and mis-modeling. NeuroImage 45, 187–198. doi: 10.1016/j.neuroimage.2008.10.065

CrossRef Full Text | Google Scholar

Luersen, M. A., and Riche, R. L. (2004). Globalized Nelder–Mead method for engineering optimization. Comput. Struct. 82, 2251–2260. doi: 10.1016/j.compstruc.2004.03.072

PubMed Abstract | CrossRef Full Text | Google Scholar

Naseer, N., Hong, M. J., and Hong, K.-S. (2014). Online binary decision decoding using functional near-infrared spectroscopy for the development of brain-computer interface. Exp. Brain Res. 232, 555–564. doi: 10.1007/s00221-013-3764-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Nelder, J. A., and Mead, R. (1965). A simplex method for functional minimization. J. Comput. 7, 308–313. doi: 10.1093/comjnl/7.4.308

PubMed Abstract | CrossRef Full Text | Google Scholar

Power, S. D., Kushki, A., and 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. doi: 10.1088/1741-2560/8/6/066004

PubMed Abstract | CrossRef Full Text | Google Scholar

Prince, S., Kolehmainen, V., Kaipio, J. P., Franceschini, M. A., Boas, D., and Arridge, S. R. (2003). Time series estimation of biological factors in optical diffusion tomography. Phys. Med. Biol. 48, 1491–1504. doi: 10.1088/0031-9155/48/11/301

PubMed Abstract | CrossRef Full Text | Google Scholar

Santosa, H., Hong, M. J., and Hong, K.-S. (2014). Lateralization of music processing with noises in the auditory cortex: an fNIRS study. Front. Behav. Neurosci. 8:418. doi: 10.3389/fnbeh.2014.00418

PubMed Abstract | CrossRef Full Text | Google Scholar

Santosa, H., Hong, M. J., Kim, S.-P., and Hong, K.-S. (2013). Noise reduction in functional near-infrared spectroscopy signals by independent component analysis. Rev. Sci. Instrum. 84:073106. doi: 10.1063/1.4812785

PubMed Abstract | CrossRef Full Text | Google Scholar

Scarpa, F., Brigadoi, S., Cutini, S., Scatturin, P., Zorzi, M., Dell'Acqua, R., et al. (2013). A reference-channel based methodology to improve estimation of event-related hemodynamic response from fNIRS measurements. NeuroImage 72, 106–119. doi: 10.1016/j.neuroimage.2013.01.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Scarpa, F., Cutini, S., Scatturin, P., Dell'Acqua, R., and Sparacino, G. (2010). Bayesian filtering of human brain hemodynamic activity elicited by visual short-term maintenance recorded through functional near-infrared spectroscopy (fNIRS). Opt. Express 18, 26550–26568. doi: 10.1364/OE.18.026550

PubMed Abstract | CrossRef Full Text | Google Scholar

Schudlo, L. C., Power, S. D., and Chau, T. (2013). Dynamic topographical pattern classification of multichannel prefrontal NIRS signals. J. Neural Eng. 10:046018. doi: 10.1088/1741-2560/10/4/046018

PubMed Abstract | CrossRef Full Text | Google Scholar

Shah, Z. Y., Wright, M. J., Thompson, P. M., McMahon, K. L., Blokland, G. G., de Zubicaray, G. I., et al. (2014). Modeling of hemodynamic response function in block design fMRI studies. J. Cereb. Blood Flow Metab. 34, 316–324. doi: 10.1038/jcbfm.2013.200

PubMed Abstract | CrossRef Full Text | Google Scholar

Soekadar, S. R., Witkowski, M., Cossio, E. G., Birbaumer, N., and Cohen, L. G. (2014). Learned EEG-based brain self-regulation of motor-related oscillations during application of trans-cranial electric brain stimulation: feasibility limitations. Front. Behav. Neurosci. 8:93. doi: 10.3389/fnbeh.2014.00093

PubMed Abstract | CrossRef Full Text | Google Scholar

Spendley, W., Hext, G. R., and Himsworth, F. R. (1962). Sequential application of simplex designs in optimization and evolutionary operation. Technometrics 4, 441–468. doi: 10.1080/00401706.1962.10490033

CrossRef Full Text | Google Scholar

Ye, J. C., Tak, S., Jang, K. E., Jung, J., and Jang, J. (2009). NIRS-SPM: statistical parametric mapping for near-infrared spectroscopy. Neuroimage 44, 428–447. doi: 10.1016/j.neuroimage.2008.08.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Li, S., Dunn, J., Li, H., Qin, W., Zhu, M., et al. (2014). The neural correlates of risk propensity in males and females using resting-state fMRI. Front. Behav. Neurosci. 8:2. doi: 10.3389/fnbeh.2014.00002

CrossRef Full Text | Google Scholar