Cortical Signal Analysis and Advances in Functional Near-Infrared Spectroscopy Signal: A Review.

Functional near-infrared spectroscopy (fNIRS) is a non-invasive neuroimaging modality that measures the concentration changes of oxy-hemoglobin (HbO) and de-oxy hemoglobin (HbR) at the same time. It is an emerging cortical imaging modality with a good temporal resolution that is acceptable for brain-computer interface applications. Researchers have developed several methods in last two decades to extract the neuronal activation related waveform from the observed fNIRS time series. But still there is no standard method for analysis of fNIRS data. This article presents a brief review of existing methodologies to model and analyze the activation signal. The purpose of this review article is to give a general overview of variety of existing methodologies to extract useful information from measured fNIRS data including pre-processing steps, effects of differential path length factor (DPF), variations and attributes of hemodynamic response function (HRF), extraction of evoked response, removal of physiological noises, instrumentation, and environmental noises and resting/activation state functional connectivity. Finally, the challenges in the analysis of fNIRS signal are summarized.

Introduction

Near-infrared spectroscopy (NIRS) is an emerging non-invasive brain-imaging methodology that utilizes near-infrared (NIR) light of 650–900 nm to determine cerebral oxygenation, blood flow, and the metabolic status of a localized region of the brain (Saager and Berger, 2008; Yamada et al., 2009; Khan et al., 2014; Molavi et al., 2014; Santosa et al., 2014; Naseer and Hong, 2015). Activation in a particular part of the brain causes an increase in the regional cerebral blood flow (rCBF) (Zhang et al., 2011b; Umeyama and Yamada, 2013; Kopton and Kenning, 2014). The rate of rCBF increase exceeds that of the regional cerebral oxygen metabolic rate (rCMRO2), which is the major cause of de-oxy hemoglobin (HbR) decrease in venous blood (Sitaram et al., 2009). Thus, cortical activation causes an increase in total hemoglobin (HbT) and oxy-hemoglobin (HbO), with a corresponding decrease in HbR. The absorption of NIR light changes with changes in the concentration of HbO and HbR during activation and rest periods. The attenuation of NIR light due to the absorption change reflects, according to the modified-Beer Lambert law (MBLL), the concentrations of HbO and HbR. Among neuro-imaging modalities, functional near-infrared spectroscopy (fNIRS)'s simplicity, portability, low cost, good temporal resolution (suitable for real-time imaging), and high signal-to-noise ratio make it a favorable option (Hu et al., 2013; Chang et al., 2014; Herff et al., 2014). fNIRS also has been considered as a potential multi-modality imaging methodology (Yunjie and Blaise, 2012). One disadvantage of fNIRS, however, is its low penetration depth. Details on the pros and cones of fNIRS can be found in Gervain et al. (2011), Barati et al. (2013), and Tak and Ye (2013).

The increase in HbR at particular area of brain is an indicator of the neuronal activity in nearby area. The detection of the neuronal activation in a particular cortical area is nothing but extraction of a specific waveform in the hemodynamic response (HR) (Ciftçi et al., 2008). In past, the canonical hemodynamic response function (cHRF) is frequently be used as desired impulse response in hemodynamic signal. Of course, it could vary in its shape, time to peak, relaxation time, and full width half maximum (FWHM). These variation in the characteristics of HRF are observed in different brain areas, among subjects and on repetition of trials. Such variations in the attributes of cHRF measured by fNIRS, has been observed in HbO concentration changes (Hong and Nugyen, 2014). The major cause of this phenomenon is the brain's continuous consciously/unconsciously processing for several tasks at the same instant of time. Even if the subject is instructed to relax and sit comfortably during an experiment, the brain consciously or unconsciously processes for many past, present and future events. Several studies have reported such findings during the analysis of fNIRS data. The difference in the dynamical 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). An additional source of such variations in hemodynamic signal measured through fNIRS, could also be as a result of certain artifacts (Yamada et al., 2009; Umeyama and Yamada, 2013). The artifacts could be related to instrumentation noise, not proper fixation of NIRS optodes, and motion of subjects likewise body tilt, breathing hold, and head nodding etc. (Yamada et al., 2009; Robertson et al., 2010; Umeyama and Yamada, 2013). Among these considerations, there is another factor that can affect the shape of HRF. This factor is known as differential path length factor (DPF). It determines the additional distance traveled by light photon due to scattered behavior of brain tissues. It is found that the wavelength dependent 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).

Additionally, NIRS signals include physiological noises associated with heart beats, respiration rhythms and low-frequency fluctuations. A special algorithm that not only to suppress physiological signals present in optical signals, measured through fNIRS, but also other unwanted (activation not related to experimental paradigm) signals due to continuous brain processing, therefore is required. In fNIRS signal analysis, most of studies have been reported in relation with reduction of physiological and instrumentation noises or to extract neuronal activation related waveform. But recent research in this field has also been turned toward the analysis of functional connectivity of brain regions during resting states (Lu et al., 2010; Hu et al., 2013). Additionally, does this resting state connection stays during task periods of a particular region (Zhang et al., 2011b; Hu et al., 2012, 2013)? Until now fNIRS has a limitation that NIRS optodes cannot cover full skull at once to study/analyze complete functional connectivity (Lu et al., 2010). Figure 1 summarizes different subfields in the area of fNIRS signal analysis for development of a standard methodology.

Figure 1

Stages of fNIRS signal analysis.

Brain optical signal

Data acquisition and pre-processing

In recent years, different types of NIRS imaging systems have been developed, which can be grouped into continuous wave (CW), frequency domain (FD), and time-resolved (TR) categories. CW-fNIRS instrument measures the concentration changes of HbO, HbR, and total hemoglobin with assumption that scattering remains constant, while FD-fNIRS and TD-fNIRS detects the absolute concentrations of HbO and HbR. TD-fNIRS is based on the principle of time of flight measurement and the most expensive instrument. The central differentiating element among these instruments is the estimation of the path length traveled by the photon due to scattering. The CW system, the least expensive, provides relative-change information in the forms of the concentrations of HbO and HbR, and was the version utilized most frequently. Figure 2 shows the geometry of fNIRS signal acquisition.

Figure 2

The geometry of fNIRS signal acquisition.

The transportation of light photon through tissue is a complex process. When photon of different wavelengths incident on tissues, the characteristics of the detected light depends upon the combination of scattering, absorption, and reflection (Cope and Delpy, 1988). It can be modeled using the modified Beer-Lambert law given as where I(λ), I(λ) are the incident and the detected light respectively, μ(λ) is the absorption, d is the source detector separation, DPF(λ) is the DPF, and G(λ) is the geometry-dependent parameter. The first step to express chromophore changes from optical signal is to find the optical density (OD) defined as (Cope and Delpy, 1988; Duncan et al., 1996; Kamran and Hong, 2014) or Considering two chromophores, i.e., HbO and HbR and assuming the phenomenon of scattering to be constant, where λ is the wavelength of the incident light and and are the extinction coefficients of HbO and HbR, respectively. By considering two different wavelength of light, the above equation could be rearranged as follows (Ye et al., 2009; Kamran and Hong, 2013; Santosa et al., 2013) where ΔHbR(k) and ΔHbR(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 two different wavelengths, , , , and indicates the extinction coefficients of HbO and HbR at two different wavelengths respectively, is the optical density variation at kth-sample time at particular wavelength (j = 1, 2), d is the source-detector separation and is the DPF at particular wavelength (j = 1, 2).

Hemodynamic response analysis

Effects of differential path length factor

The scattering behavior of human brain tissue to NIR light entails that DPF is required to correct the reading observed through fNIRS. Initially it was in practice to use DPF value between 3 and 6. Later, time of flight methodology or intensity modulated spectroscopy were used to estimate the values of DPF. Duncan et al. (1995) found the DPF values for adult heads for four different wavelengths in a population of 100 (50 males and 50 females) subjects. It was found that the value of DPF are 6.51 ± 1.13, 6.53 ± 0.99, 6.26 ± 0.88, and 5.86±0.98 for 690, 744, 807, and 832 nm, respectively. The magnitude of the DPF determines the magnitude of the calculated concentration changes (Kohl et al., 1998). Therefore, DPF has an important role for any instrument claiming accurate measurement of chromophores changes. Kohl et al. (1998) has used the key idea that DPF is proportional to the rate of change of absorbance with respect to absorption. Thus, they found wavelength dependent DPF by the ratio that rate of change of absorbance to the absorption spectrum of arterial blood. Duncan et al. (1996) analyzed 283 subjects (age between 1 day and 50 years) and developed equations to express DPF values as a function of age at different wavelengths. Their results summarizes the DPF values for four different wavelengths as under; But CW-NIRS systems could be used for different wavelengths for their equipment, thus a general equation was the requirement that could be used for any wavelength and for any age. Schroeter et al. (2003) analyzed fNIRS data from 14 young (23.9 ± 3.1 years old) and 14 elderly (65.1 ± 3.1) subjects and suggested that DPF is not only effected by age but also with different brain regions. They concluded that the hemodynamic response can be decreased by age in the frontal association cortex during functional activation and proposed to calculate effect size to analyze age-related effects in fNIRS studies. Scholkmann and Wolf (2013) realized this need and fitted a polynomial of degree three on the available data set for the values of DPF. They represented DPF as a function of wavelength and age as under The values of the unknown parameters were found by using Levenberg-Marquardt algorithm (LMA) and least absolute residuals (LAR). There results suggested that The above equation is a generalized form of DPF correction depending upon age and wavelength. This equation is advantageous to use for any researcher because any one can easily evaluate the DPF value at any wavelength and age. The published articles presenting the values of DPF for different age and wavelengths have been summarized in Table 1.

Table 1

The published references for DPF values.

References	Cortex	Subjects	Age	Wavelengths	S-D separation
van der Zee et al., 1992	Temporal frontal	10 adults, 10 Infants	22~54 Y	761 nm	2.5 cm
Duncan et al., 1995	Left forehead	100 adults, 35 Infants	21~ 9 Y and Infants	690, 744, 807, 832 nm	>4 cm
Duncan et al., 1996	Left forehead	283	1 day ~50 Y	690, 744, 807, 832 nm	4.5 cm
Cooper et al., 1996	Temporal frontal	19	23~38 weeks	730 and 830 nm	4.9 cm
Kohl et al., 1998	Right occipital	10	23~40 Y	700–1000 nm	3 cm
Scholkmann and Wolf, 2013	NA	NA	Applicable to all	Applicable to any	NA

Variations in HRF pattern

The neural activation indication, measured through fNIRS, may be confounded with individual anatomical or systemic physiological sources of variance (Heinzel et al., 2013). 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. Barati et al. (2013) observed that the variability in the stimulus condition for HbO was revealed in the slope, amplitude, and timing of the peak response. Jasdzewski et al. (2003) analyzed the difference in the dynamical shape of HRF during event-related motor and visual paradigms. Their results 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, their results have been analyzed for different values of source-detector separation. But if the source-detector separation is greater than 3 cm than the results are not much affected by DPF values (Duncan et al., 1996). Power et al. (2012) analyzed two very important questions: (1) is it possible to distinguish activation task from baseline or from other tasks? And if so (2) are the spatiotemporal characteristics of the response consistent across sessions? Their results concluded that the mental arithmetic tasks can be distinguished from base line but the characteristics of the response changes from session to session. Hong and Nugyen (2014) analyzed 19 subjects to conclude variations in the impulse responses at three different brain regions, somatosensory cortex (SC), motor cortex (MC), and visual cortex (VC). Their findings suggest that the activation- and the undershoot-peak of the HbO in MC are higher than those in SC and VC. Additionally, the time-to-peaks of the HbO in three brain regions are almost the same (about 6.76 76 ± 0.2 s) and the time to undershoot peak in VC is the largest among three.

Constrained basis set

The detection of cortical activation related waveform from neuroimaging discrete signal is nothing but a search for a consistent and specific wave pattern (Koray et al., 2008). This is equaling to fit the measured signal to a known waveform up to certain accuracy. A mismatch of such fitting might lead to misleading results. The cHRF attributes includes magnitude of initial dip, time to the first peak, time to the undershoot, magnitude of the undershoot etc. In literature, cHRF consisting of two gamma functions have been used most frequently. In which, first peak is to tackle the main response and second peak is for undershoot after the response. Likewise fNIRS, functional magnetic resonance imaging (fMRI) modeling requires flexible HRF modeling, with the HRF being allowed to vary spatially, on repetition of trials and between subjects (Woolrich et al., 2004). Thus, a Bayesian constrained frame work is described in Woolrich et al. (2004) to best choose the HRF in the measured data. Koray et al. (2008) proposed constraint GLM model parameters such as main response (time to first peak) must be within 3–8 s, not more than one positive peak, not more than two dips, initial dip magnitude must be lesser than quarter of the magnitude of onset, an undershoot after 2–8 s of time to peak and magnitude of post stimulus undershoot must be lesser than half of the magnitude of onset. A 3D volume of parameter values of the canonical basis set then supposed as prior distribution in the Bayesian analysis. Finally, Gibbs sampling in this volume is used to find out the parameter of interest. This method is advantageous as it tries to constraint the basis set with possible reduction in solution space. The solution space is reduced on the basis of HRF physical properties reported in fMRI data in past.

Extraction of evoked-response

Jobsis (1977) was the first to present an idea that there is a possibility to detect changes of cortical oxygen using NIR light. Later, Cope and Delpy (1988) designed an NIR system with four different wavelengths (778, 813, 867, and 904 nm). It is well-known fact that the neuronal activity generates an early de-oxygenation in a particular area of the brain from where the activity is started. The characteristics of HRF include an early rise after 1–2 s of stimulation and reaches to peak around 5–6 s. Finally, it starts to drop down and reaches a base line level after having a slight undershoot. The total duration of HRF for impulse stimulation is around 26–30 s. Friston et al. (1994) introduced the statistical parameter mapping (SPM) software for fMRI signal analysis, modeling the oxygen dependent signal as linear combination of two Gamma functions.

This two Gamma function model is most frequently used to account the first peak and final undershoot of the oxygen dependent waveform. Figure 3 displays two Gamma functions (on left side of Figure 3) and final cHRF by employing their linear combinations (right side of Figure 3). The standard values to generate the shape of cHRF are given in SPM (Friston et al., 1994, 1998). The mathematical form to generate this type of HRF is described below where u is the experimental paradigm, h represents the cHRF, α1 is the delay of response, α2 is the delay of the undershoot, β1 is the dispersion of the response, β2 is the dispersion of the undershoot and Γ represents the Gamma distribution. Boynton et al. (1996) presented the idea to model neuronal related fNIRS waveform by employing only one Gamma function with two free parameters.

Figure 3

Canonical hemodynamic response function (cHRF).

In contrast with fMRI, fNIRS observed signal is contaminated with several physiological noises. In most of studies, such signals are filtered out in pre-processing steps with known cutoff frequencies. But, Prince et al. (2003) modeled such signals as a linear combination of different sinusoids. In past several studies have modeled HRF using the model presented in Equations (13) and (14) in GLM based analysis framework. Jasdzewski et al. (2003) analyzed impulse response attributes in the form of a linear model. Later, Koh et al. (2007) introduced a functional optical signal analysis (fOSA) software based upon the GLM methodology. Plichta et al. (2006, 2007) presented the functional brain maps of visual cortex by employing GLM methodology with ordinary least square estimation (OLSE) to extract the values of activity strength parameters. Taga et al. (2007) analyzed the effects of source-detector separation to extract neuronal activity related fNIRS signal. Koray et al. (2008) estimated HRF by fitting constrained parameters of cHRF in Bayesian framework. NIRS-SPM (Ye et al., 2009) is the extension of SPM (Friston et al., 1994, 1998), frequently been used in fMRI analysis. This software package has employed the concept of GLM with known regressors to decompose the measured NIRS time series. The brain signal model could be represented mathematically in the form of known n-regressors, The GLM methodology has been employed quite frequently while analyzing the fMRI time series. For this purpose, a basis set including predicted HRF (pHRF) and a base line correction is been used. Furthermore, the temporal and dispersion derivatives have been added to tackle the temporal and spatial effects in HRF (Friston et al., 1994, 1998). Likewise fMRI, fNIRS instrument monitors the concentration changes of HbO/HbR, thus an identical regression vector with fMRI is used in optical signal analysis. But fNIRS time series has additional challenges of existence of physiological signals in the measured waveform. Thus, Abdelnour and Huppert (2009) described the basis set as linear combination of pHRF, a base line correction and three sinusoids for physiological signals. Hu et al. (2010), supposed regression vector to be a combination of five components; pHRF, a baseline correction and remaining three forms a set of high pass filter with cut-off frequency 0.0006 Hz. Zhang et al. (2011b, 2012) introduced the use of recursive algorithms for better extraction of neuronal related concentration changes in observed fNIRS data. Aqil et al. (2012a) supposed fNIRS signal in a standard GLM framework with estimation of activity strength parameters using recursive algorithm. They modeled the fNIRS time series as a linear combination of pHRF, first derivative of pHRF (temporal derivative), second order derivative of pHRF (dispersion derivative) and a base line correction. Kamran and Hong (2013) has modified the method and analyzed the measured optical data in the form of linear parameter varying approach with a recursive technique that can estimate activity strength parameters in a Lagrangian framework. Scarpa et al. (2013) presented the idea to consider a reference channel with a source-detector separation of <0.7 cm. This reference channel contains only the physiological signals and useful component from nearby channel could be extracted by subtracting the data measured through reference channel. Later, Kamran and Hong (2014) modeled cortical signal in the form of auto-regressive moving average with exogenous signal (ARMAX). The physiological signal have been considered as a known amplitude and frequencies and variation in HRF is modeled by employing ARMA model. Hong and Nugyen (2014) has used the same basis set as described in Aqil et al. (2012b) to model impulse response as a state-space model.

Independent component analysis (ICA) is a powerful blind signal processing technique (Katura et al., 2008). It can extract independent components from measured discrete series by using statistical concepts. Morren et al. (2004) analyzed and detected fast neuronal signal with a source-detector separation of 3 cm using ICA technique. Zhang et al. (2013) has employed ICA methodology to explore the existence of particular wave form, modeled as Gamma variants. Similarly, Santosa et al. (2013) has used ICA to extract the pHRF from regression vector including pHRF, a baseline correction and physiological noises. NIRS data analysis is performed in medical field as well for detection of different brain diseases. Machado et al. (2011) have used the GLM methodology to estimate the existence of hemodynamic responses to epileptic activity. As a general comparison, the most of above methods could be analyzed on the basis of computational cost. For example, ICA needs more computational cost as compared to ordinary/recursive least squares estimation algorithm for unknown parameter estimation. Initial dip is an important and most crucial attribute and indicator of neuronal activity. It gives the idea of a particular location which is involved in originating the said activity. Thus, using Gamma functions as basis set could be more advantageous if the factor of initial dip is added and analyzed. Table 2 summarizes the list of studies describing methodologies to extract the neuronal activity related wave pattern from fNIRS signal.

Table 2

Signal processing methodologies for extraction of evoked-hemodynamic response.

References	Methodological details
Jobsis, 1977	Possibility to detect changes of cortical oxygen using NIR light.
Cope and Delpy, 1988	Design of NIR system with four wavelengths (778, 813, 867, and 904 nm) with applying modified Beer-Lambert law for data conversion.
Friston et al., 1994	Statistical parameter mapping software for fMRI but later used for fNIRS data analysis with modifications.
Boynton et al., 1996	HRF model with one Gamma function with two free parameters.
Prince et al., 2003	Biological signals modeled as sum of sinusoids.
Jasdzewski et al., 2003	Impulse response, initial dip, and time to peak analysis in fNIRS signal.
Koh et al., 2007	A software functional optical signal analysis (FOSA) was introduced based on GLM methodology.
Plichta et al., 2006, 2007	GLM methodology with ordinary least square estimation to generate functional maps of visual cortex.
Taga et al., 2007	Analysis of effect of source-detector separation to fNIRS hemodynamic response.
Koray et al., 2008	Estimation of constrained HRF parameters in Bayesian frame work.
Abdelnour and Huppert, 2009	GLM based methodology with Kalman filter to estimate handedness.
Ye et al., 2009	GLM based NIRS-SPM software package for analysis of fNIRS data.
Hu et al., 2010	Brain functional maps by using GLM and Kalman filtering.
Zhang et al., 2011b	Recursive least squares (RLS)-empirical mode decomposition for noise reduction.
Zhang et al., 2012	RLS estimation with forgetting factor to remove physiological noise.
Aqil et al., 2012a	GLM and RLSE for estimation of brain functional maps.
Aqil et al., 2012b	Generation of cHRF using state-space approach.
Scarpa et al., 2013	Reference channel based methodology for estimation of evoked-response
Santosa et al., 2013	ICA methodology to estimate pre-defined cortical activation signal.
Kamran and Hong, 2014	Linear parameter varying model and adaptive filtering to estimate HRF and functional maps of brain.
Barati et al., 2013	Principle component analysis to continuous fNIRS data (using spline method).
Kamran and Hong, 2014	Auto-regressive moving average model with exogenous signal (ARMAX) model for cortical activation estimation.
Hong and Nugyen, 2014	State-space model for impulse response using fNIRS.

HRF model using state-space model

fNIRS measured time series is a discrete data series that could be converted into state-space model for further analysis. Aqil et al. (2012b) summarizes the state space model of optical time series using standard subspace-based approach, but the numeric values of final matrices were not been displayed. The brain model used in their work is similar with Aqil et al. (2012a). Later, Hong and Nugyen (2014) converted fNIRS cortical signal model (same as Aqil et al., 2012a) into state-space model by using standard subspace based approach. They have summarized the mathematical derivation as well and finally they used built in Matlab function to extract final state-space model and matrices of order six. They have also summarized the numeric values of finalized matrices for different brain regions as well. Kamran and Hong (2013) presented the idea that linear parameter varying model could be beneficial to tackle the time varying characteristics of the human brain signal. In their work, the measured optical data is modeled as a state-space model whose matrices are dependent upon time varying parameter. But a final state-space model has not been reported in their work. Modeling using state space methods could be beneficial over other estimation methodologies with recursive algorithms because the analysis in state space model would be much easier if a model is developed that can cater the attributes of HRF as variable parameters.

Physiological noises

fNIRS data analysis incorporates an additional challenge of temporal correlation presented in the data due to the physiological signals. The physiological noises includes cardiac beat, respiration rhythm and low frequency fluctuations known as Mayer waves. In most of studies, the physiological signals are pre-filtered by using standard signal filtering techniques. Prince et al. (2003) has presented the idea to model the biological signals as a set of sinusoids. Zhang et al. (2007) has proposed to nullify the effects of global interferences by using multi-separation probe configuration (placing a detector close to the source) and adaptive filtering. Abdelnour and Huppert (2009) included physiological signals as known regressors in their regression set. Hu et al. (2010) have added a set of high pass filter with cut-off frequency 0.0006 Hz to tackle the physiological signals. Zhang et al. (2011b) analyzed multi-distant source-detector separation further by decomposing short distance source-detector measurement into intrinsic mode functions (IMFs). An estimate of global interference is derived by analyzing weight coefficients of IMFs. Cooper et al. (2012) has analyzed the simultaneous fNIRS and fMRI recordings to reduce the physiological effects. They calculated the variance of the residual error in a GLM of the base line fMRI signal and the observed variance is reduced by incorporating NIRS signal in the model. Zhang et al. (2012) presented to remove physiological effects from simulated fNIRS data set of near and far detectors using recursive algorithm. Yamada et al. (2012) proposed that functional and systemic responses could be separated on the basis of a negative and positive linear relationship between HbO and HbR changes of the functional and the systemic signals. Later, Kirilina et al. (2012) included an additional predictor to account for systemic changes in the skin to analyze time course, localization and physiological origin of task related superficial signals in fNIRS measured data. They found that skin blood volume depends upon the cortical state. Additionally they found that origin of the task related systemic signals in fNIRS are co-localized with veins draining the scalp. Frederick et al. (2012) generated regressors for systemic blood flow and oxygenation fluctuation effects by applying a voxel-specific time delay to concurrently acquired fNIRS-fMRI time series. Kamran and Hong (2014) added three sinusoids with known frequencies and amplitudes as exogenous signal in their ARMAX frame work. It is very important to find out the frequencies and amplitudes of existing sinusoid in measured data in addition with its origin. Since, the frequency of experimental paradigm coincident with physiological noises could result in the generation of harmonics. Most of the previous studies were using a generic idea regarding a fixed frequency and amplitude of known sinusoid to cater the physiological signals in fNIRS data. A recently published Kamran et al. (2015) is more advantageous estimation algorithm because it allows the user to estimate the frequency and amplitude of physiological signals automatically from measured data instead of using a fix pattern. The studies related to removal of physiological noises from fNIRS signal have been summarized in Table 3 describing cortical area, number of subjects, nature of mental task, methodology used and source-detector separation.

Table 3

Physiological noise estimation and reduction in fNIRS measured signal.

References	Cortex	Subjects	Task	Methodology	S-D separation
Prince et al., 2003	motor	5	Hand tapping/Rest	Stochastic model with extended Kalman filter	>2 cm
Zhang et al., 2007	Five layer slab model	Simulated data	Block design paradigm	Multi-separation probe configuration and Monte Carlo	1.5~4.5 cm
Abdelnour and Huppert, 2009	Motor	3	Finger tapping	GLM with Kalman estimator	3.1 cm
Zhang et al., 2011b	Five layer slab model	Simulated data	Block design paradigm	Multi-distance approach with empirical mode decomposition	0.1 cm
Zhang et al., 2012	Five layer slab model	Simulated data	Block design paradigm	Recursive least square estimation (RLSE) filtering	0.15~4.5 cm
Yamada et al., 2012	Primary motor	7	Finger tapping	Negative and positive correlation	1~4 cm
Frederick et al., 2012	Right frontal lobe	6	Resting state	Voxel-specific time delay	0.1 and 3 cm
Kamran and Hong, 2014	Motor	6	Finger tapping	ARMAX	~ 2.5 cm
Erdogan et al., 2014	Pre-frontal	18	Mental arithmetic	Extended superficial signal regression method.	2.5 cm
Kirilina et al., 2013	Frontal lobe	15	German words recognition	Time-domain fNIRS with wavelet coherence analysis.	3 cm
Bauernfeind et al., 2014	Motor cortex	12	Cue-based right hand (RH) and both feet (FE) motor execution	A common reference method, ICA and transfer function models	3 cm
Zhang et al., 2013	Five layer slab model	Simulated data	Epoch block	Multi-distance probe configuration and ICA	0.5 and 4.5 cm
Barker et al., 2013	Pre-frontal	22	Resting state	Regression analysis using GLM	–
Tong et al., 2013	Middle hand and left Big toe	7	Resting state	Group ICA	1.5 cm
Santosa et al., 2013	Pre-frontal	8	Arithmetic task	ICA with pre-defined regressors	2.2, 2.5, and 4.3 cm
Scarpa et al., 2013	Motor	10	Key pressing with left or right index finger	Reference channel based noise removal	1.5 and 3 cm
Kirilina et al., 2012	Pre-frontal	15	Semantic Continuous performance task	Concurrent fNIRS and fMRI with Bio-signals	3 cm
Cooper et al., 2012	Frontal and temporal lobe	7	Resting state	Variance of residues in GLM for concurrent fNIRS and fMRI	1 and~ 3cm
Hu et al., 2010	Motor	5	Finger tapping	Kalman filters and GLM	2 cm
Scarpa et al., 2010	Parieto-occipital	13	Visual graphics	Bayesian filtering	3 cm
Katura et al., 2008	Sensrimotor	30	Finger tapping	ICA	3 cm
Saager and Berger, 2008	Left pre-frontal	21	Resting state	Multi-detector CW-fNIRS	3.3 cm
Zhang et al., 2005	Sensrimotor	10	Finger movement task	Eigen vector based spatial filtering	> 3 cm
Cui et al., 2010	Motor	10	Finger tapping with head motion	Maximization of negative correlation
Haeussinger et al., 2011	Frontal	24	Working memory	Identification of channels with major extra-cranial signal contributions

Resting state functional connectivity

Human brain generates continuous low frequency fluctuations during resting state. These low frequency fluctuations could be used as an informative source to understand the mechanism of different brain regions. The reason that these fluctuations are useful because they are correlated with different brain regions. These findings opened up a new way of thinking to explore a topic named as resting state functional connectivity (RSFC). There is no consensus available at the moment to specifically designating a range of frequencies for resting state low frequency fluctuations of hemodynamic measured waveform (Fox and Raichle, 2007; Lu et al., 2010). fMRI studies have reported a high level of inter-hemispheric correlations in different brain regions (Damoiseaux et al., 2006; De Luca et al., 2006). White et al. (2009) have reported their analysis regarding RSFC in motor in visual brain cortices. They calculated correlations using Pearson correlation coefficient. Their results suggest inter-hemispheric correlations exist in both motor and visual networks. Lu et al. (2010) analyzed the RSFC maps of the sensory motor and the auditory cortices using seed-based correlation analysis and data driven cluster analysis during resting state and motor-localizer task sessions. Their results suggested RSFCs were detected both within the ipsilateral and between the bilateral sensorimotor seed-regions. Additionally, it was found that significant correlation exist within the ipsilateral and between the bilateral temporal auditory cortices but not between temporal auditory areas. Zhang et al. (2011a) raised an issue of reliable RSFC maps. They analyzed test-reset reliability at three different scales; maps cluster and channel wise at individual and group levels. Their finding suggest that one should be very careful when interpreting the individual channel wise RSFC. But individual level and group level RSFC has excellent map-/cluster wise reliability. The trial-to-trial variability (TTV) in fNIRS signal exist even if the experimental procedure is kept constant. Hu et al. (2013) suggested to reduce TTV using RSFC information. They concluded that low frequency fluctuations are significant source of TTV and TTV decreases after removing the effects of bilateral connectivity. Since, fNIRS optodes cannot cover whole head surface which makes it difficult to analyze the RSFCs between all different brain regions (Lu et al., 2010). One possible way of covering whole skull is to increase the number of optodes. But this shall affect the temporal resolution of the equipment.

Environmental and instrumental effects/artifacts/noises

fNIRS measured data series includes hemodynamic signal related to certain artifacts. These artifacts could be related to biological processes or from outer sources. Instrumental noise is one of major external sources whose affect could be reduced/nullify by proper calibration of the instrument. The other artifacts could be as a result of not good contact between skull and NIRS optodes. The uncoupling of optodes and skull is the source of fluctuations in the detected intensity, leading to uncorrected results. Thus, it must be insured that optodes and skull have a proper contact at correct angle. Removal of unnecessary hair at contact could also improve the quality of the signal. Another artifact is due the motion (tilt in the body, slight head moment, and breath holding) of the subject. This motion cause changes in the blood flow which is a major reason of fluctuations in the measured hemodynamic response. A crude way to remove motion artifacts is to average a certain number of experiments that its effects could be nullified. But for real-time BCI applications, it is necessary to estimate the actual contributions of the motion artifacts in single-trial. Yamada et al. (2009) presented a theoretical analysis of optical signal using Monte Carlo simulation. They proposed a multi-distant probe arrangement can reduce/eliminate artifacts in fNIRS measured data. Later, Robertson et al. (2010) experimented in a co-located channel configuration to analyze the known motion artifacts around three axis. They found motion related hemodynamic signal is detectable at co-located channels but not at unique channel. Cui et al. (2010) studied the effect of artifacts on fNIRS data induced by the head motion and they found that oxy- and deoxy-Hb are generally negatively correlated, head motion causes the correlation to become more positive. They proposed a correlation based signal improvement method to maximize the negative correlation between oxy- and deoxy-Hb signals. Haeussinger et al. (2011) developed a method to identify channels with major extra-cranial signal contributions and subtracted the average of these channels from all channels to obtain improved fNIRS signals. Thus, it would be concluded that moment in different direction caused changes in the absorption of light at different brain regions.

Statistical significance and functional maps

fNIRS signal are highly corrupted by several measurement noises and physiological interferences. Therefore, a careful statistical analysis is required to extract neuronal-activity related signal from observed optical data (Tak and Ye, 2013; Kamran and Hong, 2014). The existence of a particular response HRF(k) in the measured data is found by a statistical analysis known as t-test (Hu et al., 2010; Kamran and Hong, 2013; Santosa et al., 2013, 2014; Hong and Nugyen, 2014). 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 > tcritical and p < 0.05). Thus, it is equivalent of testing a null hypothesis Ho with proper statistics i.e.,

Null hypothesis H : β1 = 0 (16)

Alternative hypothesis: β1 ≠ 0.

Finally, t-value could be evaluated as under where SE is the standard error of the estimated coefficient.

However, in practical situations while analyzing fNIRS data, multiple comparison problems are often required to be addressed. Thus, there is chance that such analysis shall result in increase of inference error. Therefore, in addition to analyze the multi-data gathered through fNIRS modality, it is also required to put necessary checks that each individual data is analyzed with more care and strong level of evidence so that to reduce inference error in proceeding steps. For example, in Plichta et al. (2006), the Bonferroni correction and the Dubey/Armitage–Parmar alpha boundary were used for statistical inference of activated channels to estimate the statistical significance of fNIRS response during task periods. A detailed review on statistical analysis of fNIRS data could be found in Ye et al. (2009) and Tak and Ye (2013).

Challenges

fNIRS signal is not consistent among subjects, repeated trails and repetition of experiment even if the conditions are assumed to be similar. Therefore, the optical signal model constitutes an additional challenge to the researchers working in this area. Since, fNIRS time series is a combination of several physiological signals (Hu et al., 2010). These physiological signals are cardiac beat (~ 1 Hz), respiratory rhythm (~ 0.2–0.3 Hz) and low frequency fluctuations (< 0.1 Hz) among others. The first step is to estimate the frequency of particular signal present in the observed time series. In addition to this, it is also required to estimate the amplitudes of the signal present in the data. The next target is to find out the attributes of the HRF. There are several important characteristics, e.g., initial dip, FWHM, time to peak, peak height, post stimulus undershoot etc. Thus, a non-linear model is required to incorporate such attributes. For instant, two Gamma function model is most attractive because each Gamma function represents each peak (actual response and post dip). Furthermore, an iterative non-linear optimization algorithm is needed to estimate the free parameters in the model with significant accuracy. Finally, more precise statistical support is required to state that the estimation in the model is statistically significant. Another important factor is the design of experimental paradigm. The frequency in experimental paradigm must be different from physiological signals to avoid harmonics. Brain functionality is complex and coupled non-linear system. The analysis of different brain region's coupling is also a fundamental step toward improvement of fNIRS analysis. Considering the fact that fNIRS optodes cannot cover full skull, the number of optodes could be increased but it shall reduce the temporal resolution as well. Thus, an optimal source-detector separation is needed to establish for significant and maximum surface analysis.

Conclusion

Brain engineering is a multi-disciplinary field with a focus to extract useful information from cortical signal observed by neuroimaging equipment. In this article recent advancements in the analysis of the optical signal observed through fNIRS are summarized. It is important for new researcher to understand the importance of pre-processing steps, effects of DPF and other factors during analysis. The recent conclusion for such factor is also presented. Additionally, different methodologies that have been developed in past to extract the neuronal activation related waveform (pre-processing steps, effects of DPF, variations and attributes of hemodynamic response function (HRF), extraction of evoked response, removal of physiological noises, instrumentation, and environmental noises, and resting/activation state functional connectivity), are summarized. Since systemic, instrumentation, and environmental noises effect the measured signal and the analysis. Therefore, reduction/removal of such noises must be performed carefully. Special consideration must be given for the selection of experimental paradigm to avoid physiological harmonics. Several methodologies that have been reported in past decade for noise removal, have also been summarized here. It is well-known by fMRI studies that different brain regions have connections during resting and task periods. Thus, it is very important to analyze such connections of brain using fNIRS as well. A brief review of RSFC is also added in this article.

Author contributions

MK did the literature review and wrote the paper. MN participated in literature review and revising the article. MJ has suggested the theoretical aspects of the current study and supervised all the process from the beginning. All the authors have approved the final manuscript.

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.