Differential Path-Length Factor's Effect on the Characterization of Brain's Hemodynamic Response Function: A Functional Near-Infrared Study

Introduction

Optical spectroscopy is an emerging neuro-imaging modality that can indicate cortical functionality with good temporal resolution relative to the other modalities (e.g., Mannan et al., 2016a; Meszlényi et al., 2017). The optical signal observed through functional near-infrared spectroscopy (fNIRS) reflects the interaction of light with matter, and will differ according to the properties of matter (Maikala, 2010). fNIRS utilizes near-infrared (NIR) light in the spectral range of 650–900 nm. Thus, absorption and scattering are among the fundamental characteristics/properties that produce signal/NIR light attenuation (van der Zee et al., 1990; Kohl et al., 1998). It is well-known that biological tissue is a medium that highly scatters NIR light (Duncan et al., 1995). The absorption of NIR light is mainly caused by two chromophores in the blood, namely, oxy-hemoglobin (HbO) and deoxy-hemoglobin (HbR) (Schroeter et al., 2004). Continuous wave near-infrared spectroscopy (CW-NIRS) measures the concentration changes of HbO (ΔHbo) and HbR (ΔHbR), respectively, at each time step, making it at attractive option for brain-computer interface (BCI) applications, among others (Scholkmann and Wolf, 2013; Talukdar et al., 2013; Naseer and Hong, 2015; Mannan et al., 2016b; Cavazza et al., 2017; García-Prieto et al., 2017).

The fNIRS system consists mainly of a source-emitting NIR light of dual or more wavelengths of particular frequencies, along with detector(s) to measure the intensity of attenuated light. NIR light of multiple wavelengths (e.g., 760 and 830 nm) is thrown onto the surface of the scalp, with a physical connection between the source and scalp. This incident NIR light passes through the scalp, skull, cerebral spinal fluid (CSF), gray matter and white matter, a portion of it being received at a nearby photo-detector. It is assumed that the NIR light emitted by the source and detected by a near-by detector follows a banana-shaped path (Kamran et al., 2016). Thus, the actual path traveled by NIR light is much longer than the source-detector separation on the surface of the scalp (Delpy et al., 1988; Hiraoka et al., 1993; Chatterjee et al., 2015; Piao et al., 2015).

It is necessary to accurately estimate/measure the actual distance traveled by NIR light (Kohl et al., 1998). A parameter, namely the differential path-length factor (DPF), has been used in previous studies to cover the extra distance traveled by NIR light (Kamran et al., 2016). This is a wavelength-dependent scaling factor that indicates how many times farther (than the actual source/detector separation) the detected NIR light has traveled (Scholkmann and Wolf, 2013). It has been common practice to use a DPF value between three and six (Delpy et al., 1988; Duncan et al., 1995; Scholkmann and Wolf, 2013). The value of DPF is generally determined through measurement of the mean time of flight (TOF) of a pico-second pulse of NIR light traveling through biological tissue (van der Zee et al., 1990, 1992; Ferrari et al., 1992). Another possibility for DPF measurement is via frequency-domain NIRS (FD-fNIRS) (Lackowicz and Berndt, 1990; Duncan et al., 1993, 1995). Duncan et al. (1995) determined the values of DPF at four different wavelengths (i.e., 690, 744, 807, 832 nm) for adult arm, head, leg, and new-born-baby head, respectively. Later, Duncan et al. (1996) experimented with 283 subjects (age range: 1 day ~ 50 years) and developed four different equations for the age dependency of DPF at particular wavelengths. The equations are as follows:

$\begin{array}{rcll}DP{F}^{690}& =& 5.38+0.049^{*}(A^{0.877}),& (1)\end{array}$

$\begin{array}{rcll}DP{F}^{744}& =& 5.11+0.106^{*}(A^{0.723}),& (2)\end{array}$

$\begin{array}{rcll}DP{F}^{807}& =& 4.99+0.067^{*}(A^{0.814}),& (3)\end{array}$

$\begin{array}{rcll}DP{F}^{832}& =& 4.67+0.062^{*}(A^{0.819}),& (4)\end{array}$

Later, Kohl et al. (1998) concluded that the DPF is the ratio of the attenuation change to the change in absorption coefficient, which is to say, ($DPF=\frac{\partial A}{\partial {\mu }_a}{r}^{-1}$), where A is the attenuation, [μ_a] is the absorption coefficient, and r is the source-detector separation. Schroeter et al. (2003) analyzed the age dependency of the hemodynamic response function on the DPF. Their results suggest that aging decreases the hemodynamic response and that the decrease might be caused by the changing of tissue properties with age. In their subsequent work, Schroeter et al. (2004) pointed out that the DPF can vary intra-individually between several pixels and inter-individually between subjects of the same age. Their observation was based on previous studies (Essenpreis et al., 1993; Zhao et al., 2002). Since none of the commercially available optical imaging systems measures DPF for each pixel, such studies are impossible (Schroeter et al., 2004). The multiple layers of the human scalp and skull can affect the DPF, thereby leading to systematic errors (Talukdar et al., 2013). Recently, Scholkmann and Wolf (2013) modeled the DPF as a function of age and wavelength, deriving the following equation (Equation 5) for calculation of DPF at any wavelength and age A:

$\begin{array}{rcl}DPF(\lambda ,A)& =& 223.3+0.05624A^{0.8493}-5.723^{*}10^{-7}{\lambda }^3\end{array}$

$\begin{array}{rc}+0.001245{\lambda }^2-0.9025\lambda .& (5)\end{array}$

It is obvious from the previous studies that the DPF can affect optically measured signals and incur, thereby, error and inaccurate results (Jasdzewski et al., 2003). Therefore, it is very important to analyze the effects of the DPF on NIRS waveforms. In the current study, the authors analyzed the effect of DPF on optically measured signals. In the first step, simulated data sets were generated using a method described in previous studies (Prince et al., 2003; Kamran et al., 2015). The simulated HbO signal was separated into its constituent optical densities under different constraints, and the accuracy of the results were confirmed before further processing. The results were observed for variations of DPF${}^{{\lambda }_1}$ and DPF${}^{{\lambda }_2}$ and compared with the actual signals. Later, NIRS experiments were performed on six healthy subjects. The problem was formulated as an iterative optimization problem for estimation of the best possible fit of HRF and physiological noises with free parameters. Comparative results were derived, and a detailed analysis was performed for different values of DPF${}^{{\lambda }_1}$ and DPF${}^{{\lambda }_2}$. The results serve to comprehensively summarize the effect of DPF variations on HRF.

Theory

Derivation of Concentration Changes of Chromophores

The intensity of attenuated NIR light received by a detector depends on the reflection, absorption and scattering properties of the tissues that light has passed through (Delpy et al., 1988; Duncan et al., 1996; Maikala, 2010; Kamran et al., 2016). This phenomenon can be modeled according to the modified beer-Lambert law (MBLL) as

$\begin{array}{rc}OD=ln(\frac{I_{in}(\lambda )}{I_o(\lambda )})& (6)\end{array}$

Or

$\begin{array}{rc}OD=\Delta {\mu }_a(\lambda )dDPF(\lambda )+G(\lambda )& (7)\end{array}$

where I_in(λ) and I₀(λ) are the incident and attenuated detected light, respectively, μ_a is the absorption coefficient, d is the source-detector separation, and G(λ) is the geometrical parameter for the light's scattering properties. The absorption coefficient is the sum of the products of the extinction coefficients and the concentrations of the different chromophores present in the medium under experimentation (Scholkmann and Wolf, 2013; Kamran et al., 2016). The above equation can be modified, by scattering to be considered as a constant (Maikala, 2010). Finally, we obtain

$\begin{array}{rc}\Delta O{D}^{{\lambda }_i}=({\epsilon }_{HbO}^{{\lambda }_i}\Delta HbO+{\epsilon }_{HbR}^{{\lambda }_i}\Delta HbR)dDPF(\lambda )& (8)\end{array}$

where [λ_i] is the wavelength of incident NIR light, and ${\epsilon }_{HbO}^{{\lambda }_i}$ and ${\epsilon }_{HbR}^{{\lambda }_i}$ are the extinction coefficients of HbO and HbR, respectively. Assuming two wavelengths λ₁ and λ₂, we can rewrite the above equation as

$\begin{array}{rc}\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})}& (9)\end{array}$

and

$\begin{array}{rc}\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})}& (10)\end{array}$

where ΔHbOⁱ(k) and ΔHbRⁱ(k) are the relative concentration changes of HbO and HbR, respectively, k is the step time, i represents the ith-channel of the source-detector set, λ₁ and λ₂ are NIR light of 760 and 830 nm wavelength, respectively, ${\epsilon }_{HbO}^{{\lambda }_1}$, ${\epsilon }_{HbR}^{{\lambda }_1}$, ${\epsilon }_{HbO}^{{\lambda }_2}$, and ${\epsilon }_{HbR}^{{\lambda }_2}$ represent the extinction coefficients of HbO and HbR at two different wavelengths, respectively, $\Delta {OD}^{{\lambda }_j}(k)$ is the optical density variation at the kth-sample time and a particular wavelength (j = 1, 2), lⁱ is the source-detector separation, and ${DPF}^{{\lambda }_j}$ is the differential path-length factor at a particular wavelength (j = 1, 2).

Hemodynamic Response Model

Since the fNIRS-measured cortical signal resembles one measured through functional magnetic resonance imaging (fMRI), Friston et al. (1994) introduced statistical parameter mapping (SPM) software for analysis of fMRI data. The general characteristics of the hemodynamic response (HR) for stimulus are well-established in the literature. There is early de-oxygenation at 1–3 s, followed by a positive peak at around 5–6 s, a drop-down to the baseline, and, finally, a settling down at 25–30 s (Friston et al., 1998; Koray et al., 2008). The canonical shape of this model can be represented as a linear combination of two Gamma functions (Friston et al., 1994; Kamran et al., 2015),

$\begin{array}{rc}HRF(k)=h{(k)}^{*}u(k),& (11)\end{array}$

$\begin{array}{rc}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]& (12)\end{array}$

where u he experimental procedure, h is the canonical hemodynamic response function (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.

Optical Signal Model

In the case of fNIRS, the observed processed signal is a linear combination of cortical activity relevant to a particular experiment, physiological signals and certain unknown signals termed as “Gaussain noise.” The signal related to cortical activity is modeled in Equation (11). Physiological noises are respiratory rhythm, heartbeat, and low-frequency Mayer waves (Prince et al., 2003; Abdelnour and Huppert, 2009; Kamran et al., 2015). These physiological noises are periodic in nature for healthy human beings. Therefore, their corresponding signals can be added as sinusoidal signals (Prince et al., 2003; Kamran et al., 2015). Accordingly, the optical signal can be modeled mathematically as

$\begin{array}{rcl}{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)\end{array}$

$\begin{array}{rc}+a_{m}sin(2\pi f_{m}k)+{\epsilon }^i(k)& (13)\end{array}$

where $y_{HbO}^i$ is the observed HbO time-series at the ith-channel, a₀ is the baseline correction, a₁ is the activity-strength parameter, a_c, a_r, a_m, f_c, f_r, f_m are the amplitudes and frequencies of the cardiac, respiratory and Mayer waves, and εⁱ(k) is the zero mean Gaussian noise at the kth-sample time.

Experimental Procedure and Paradigm

Six healthy subjects participated in the experiments. Their average age was 28 years with a standard deviation of 7 years. Each subject was asked about his/her medical history, and it was determined that none had had any neuronal disorder. All of the subjects had 6 × 6 eyesight with/without glasses. A written consent attesting to his willingness to participate in the experiment was signed by each subject. The experiment was conducted in accordance with the latest version of the Declaration of Helsinki. The study was approved by ethical Institutional Review Board, Pusan National University, South Korea. As it is standard good practice, all of the subjects were informed of the experimental details for the best possible results. They were advised to remain calm while sitting on a comfortable seat, so as to avoid artifacts and noise. The experiment included an initial rest period of 15 s followed by four task/rest trials. Each trial included a 15–s task session followed by a 15–s rest session. During the task sessions, different arithmetic operations with answers were shown; the subject was instructed to tap his right index finger if the answer was correct, and to keep the finger still if the answer was incorrect. It was observed that all of the subjects properly tapped their finger only on display of a correct answer. The arithmetic operations were shown on a screen positioned 100 cm away from the subject. The data were acquired at a sampling rate of 1.81 Hz. Later, it was resampled at a rate of 100 Hz for further processing.

The algorithm was verified on simulated data sets of different combinations of stimulus patterns and task durations (10, 20, and 30 s). The simulated data sets were generated using Equations (11–13) according to the method described in Kamran et al. (2015).

Problem Formulation

There are three commercially available NIRS systems: time domain (TD), frequency domain (FD), and continuous wave (CW) (Talukdar et al., 2013). In continuous-wave NIRS (CW-NIRS), usually, two different wavelengths of NIR light are thrown onto the scalp, and their corresponding changes in optical density are received by detector/detectors. In that way, the optical densities corresponding to wavelengths λ₁ and λ₂ are recorded (Ye et al., 2009). These optical densities $\Delta {OD}^{{\lambda }_1}$ and $\Delta {OD}^{{\lambda }_2}$ are used to obtain the concentration changes of HbO and HbR by Equations (9, 10). The problem with respect to simulated data sets is formulated as a cost function j₁,

$\begin{array}{rc}{J}_1={\displaystyle \sum _{k=1}^N}{\left\{y_{HbO}^i(k)-0.0375\Delta O{D}^{{\lambda }_2}(k)+0.01644\Delta O{D}^{{\lambda }_1}(k)\right\}}^2.& (14)\end{array}$

The Equation (14) is solved for optical densities $\Delta {OD}^{{\lambda }_1}$ and $\Delta {OD}^{{\lambda }_2}$] using an optimization algorithm to minimize J₁ under certain constraints. The above equation has infinitely many solutions, which are divided into four categories. Correspondingly, four differently constrained problems are solved for optical densities: (1): $\Delta {OD}^{{\lambda }_1},\Delta {OD}^{{\lambda }_2}>0$, (2): $\Delta {OD}^{{\lambda }_1}>,\Delta {OD}^{{\lambda }_2}<0$, (3): $\Delta {OD}^{{\lambda }_1}<0,\Delta {OD}^{{\lambda }_2}>0$, and (4): $\Delta {OD}^{{\lambda }_1}<0,\Delta {OD}^{{\lambda }_2}<0$. The solutions for Equation (14) are fed into the equation

$\begin{array}{rc}{y}_{HbO}^i(k)=0.2170\frac{\Delta O{D}^{{\lambda }_2}(k)}{DP{F}^{{\lambda }_2}}-0.1015\frac{\Delta O{D}^{{\lambda }_1}(k)}{DP{F}^{{\lambda }_1}},& (15)\end{array}$

and the results are compared with the actual data to verify the correctness of the solutions for actual values of ${DPF}^{{\lambda }_1}$ and ${DPF}^{{\lambda }_2}$. Later, different values of ${DPF}^{{\lambda }_1}$ and ${DPF}^{{\lambda }_2}$ are used to analyze the effects of these parameters on HRF attributes. In the case of real data sets, most of the previous studies have utilized the outputs of Equations (9, 10) for further processing and modeling of HRF signals (Jasdzewski et al., 2003; Abdelnour and Huppert, 2009; Kamran and Hong, 2013, 2014; Kamran et al., 2015). In the present study, a proposed alternative means of optimization for estimation of the HRF directly from given optical densities of real data sets was tested. For this purpose, a cost function J₂,

$\begin{array}{l}{J}_2={\displaystyle \sum _{k=1}^N\{(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))-0.2170\frac{\Delta O{D}^{{\lambda }_1}(k)}{DP{F}^{{\lambda }_2}}-0.1015\frac{\Delta O{D}^{{\lambda }_2}(k)}{DP{F}^{{\lambda }_1}}{\}}^2,\end{array}$

was defined, which can be formulated in an optimization problem according to the following constraints:

$\begin{array}{rc}\begin{array}{c}min{J}_2({\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\\ C_1:2\le {\alpha }_1\le 10,C_7:0\le a_c\le 2,\\ C_2:6\le {\alpha }_2\le 20,C_8:0\le a_r\le 2,\\ C_3:0.5\le {\beta }_1\le 2,C_9:0\le a_m\le 2,\\ C_4:0\le {\beta }_2\le 1.5,C_{10}:0.5\le f_c\le 1.5,\\ C_5:0\le a_o\le 20,C_{11}:0.2\le f_r\le 0.3,\\ C_6:0\le a_1\le 15.C_{12}:0.09\le f_m\le 0.1.\end{array}& (16)\end{array}$

The optimal values of the free parameters (${\alpha }_1^{*},{\alpha }_2^{*},{\beta }_1^{*},{\beta }_1^{*},{\alpha }_0^{*},{\alpha }_1^{*},{\alpha }_c^{*},{\alpha }_m^{*},{\alpha }_r^{*},f_c^{*},f_r^{*},f_m^{*},$) for specific values of $DP{F}^{{\lambda }_1}$ and $DP{F}^{{\lambda }_2}$ are estimated using an improved version of the simplex method [later named the Nelder-Mead simplex method (NMSM)]. The iteration of NMSM can be performed in three steps, namely, ordering, centroid, and transformation. The details on this algorithm are available in Haftka et al. (1990), Lagarias et al. (1998), Luersen and Riche (2004), and Kamran et al. (2015).

Results

A schematic of the algorithm is shown in Figure 1. It displays the general idea of the study. Each time different set of values of DPF were used to analyze the output of neuro-activation optimization model and results were compared to comprehensively analyze the effects. The details of the experimental paradigm and source-detector separation and localization are shown in Figure 2. Sixteen set of source-detector pairs were utilized to measure motor cortex. The channel 3 is located at C3 location in 10–20 system. Figure 3 plots the simulated HRFs and their corresponding cHRFs for different task sessions. Figures 4–7 plot the effects of the variation of $DP{F}^{{\lambda }_1}$ on the HRFs corresponding to different stimuli (St₁-St₅). Similarly, Figures 8–11 show the effects of variations of $DP{F}^{{\lambda }_2}$ on HRFs corresponding to different stimuli (St₁-St₅). The x-axis in Figures 4–11 shows the sample time as mentioned before that the data is resampled at 100 Hz. Thus, each step of sample time corresponds to 0.01 s. It is obvious to analyze from Figures 4–11 that DPF has strong relation with the shape of HRF. A change in DPF value can affect the peak value of the main response as well as the depth of post-stimulus undershoot. In some cases DPF has no effect on main response (Figures 6, 10) and in some cases it has nothing to do with post-stimulus undershoot (Figures 7, 11). The reason behind this is that positive output is not possible in case: 3 (See section Discussion) and negative output is not possible in case: 4 (See section Discussion). Figure 12 summarizes the results for the six healthy subjects with respect to the effects of $DP{F}^{{\lambda }_1}$ change and DPF1 in the label represents the $DP{F}^{{\lambda }_1}$. Figure 13 plots the effects of $DP{F}^{{\lambda }_2}$ variation on real data sets and DPF2 in the label represents the $DP{F}^{{\lambda }_2}$. Figure 14 projects the effects of changes in the peak values of HRF as $DP{F}^{{\lambda }_1}$ (top plot) and $DP{F}^{{\lambda }_2}$ (bottom plot) are varied. The concentration change of HbO is calculated by using $DP{F}^{{\lambda }_1}$ and $DP{F}^{{\lambda }_2}$ from Equation (9). Figure 14 top plot displays the effects of variation of $DP{F}^{{\lambda }_1}$ keeping $DP{F}^{{\lambda }_2}$ at fixed value. The value of $DP{F}^{{\lambda }_2}$ is displayed by dark blue vertical line. Similarly, the bottom plot shows the value of $DP{F}^{{\lambda }_1}$ by vertical blue line while displaying the variation of $DP{F}^{{\lambda }_2}$ for all subjects.

FIGURE 1

Figure 1. Schematic of algorithm.

FIGURE 2

Figure 2. (A) Experimental paradigm, (B) source-detector localization and separation.

FIGURE 3

Figure 3. Simulated data with initial rest (10 s) different stimulus (St₁-St₅) and task durations (10, 20, and 30 s) followed by 30 s rest, 10 s activation (top left), 20 s activation (top right), 30 s activation (bottom left), and canonical hemodynamic response to all stimuli (bottom right).

FIGURE 4

Figure 4. Results for both positive cases (Case: 1) with variation in value of $DP{F}^{{\lambda }_1}$ (3–8) under stimulation St₁-St₅ (top left, top middle, top right, bottom left, and bottom middle).

Discussion

fNIRS is an emerging neuro-imaging methodology that can be used in medical as well as BCI applications (Kamran et al., 2016). The observed time-series of optical data includes unnecessary (unrelated to cortical activity) information. This information is noise that renders extraction of the neuron-related information on a particular task more challenging. Noise is of different types: physiological noise, motion artifacts, systematic noise, and others. There also exist other factors operative in the process of converting fNIRS' observed optical densities, those can affect signal shape. The DPF is the important parameter in the conversion step, as it can affect signal attributes (Jasdzewski et al., 2003). The DPF is included to cover the actual distance traveled by light in the signal acquisition process (Maikala, 2010; Kamran et al., 2016). Several studies have reported the effect of the DPF on fNIRS signals (Lackowicz and Berndt, 1990; Duncan et al., 1993, 1995; Kohl et al., 1998; Jasdzewski et al., 2003; Schroeter et al., 2003, 2004; Talukdar et al., 2013). Hiraoka et al. (1993) reported that DPF has only valid for a homogenous medium but real human brain and head tissues have different optical properties. Thus, accurate quantification of NIRS signal is possible with information regarding nature of light transport through in-homogenous medium. Chatterjee et al. (2015) analyzed a Monte-Carlo based computational model within a single layer of tissue like human brain. Their finding concluded that optical path changes by changing the source-detector separation. The brain is highly scatter medium of NIR light. Thus, variation in path-length traveled by light photons is also possible. DPF depends upon the extra path traveled by light photons. Schroeter et al. (2004) commented on a very important point that DPF can vary intra-individually between several pixels and none of commercially available optical imaging system measures DPF for each pixel. Therefore, it is necessary to analyze its effects on the signal. In the present study, the effects of DPF change were analyzed for both simulated data sets and the healthy subjects' real data sets. The simulated data sets were generated using the method described in Kamran et al. (2015). These data sets were subjected to Equation (14) to solve for optical densities under constraints. Four different cases were resolved with Equation (14), and the results were compared with actual data using Equation (15). Later, different values of $DP{F}^{{\lambda }_1}$ and $DP{F}^{{\lambda }_2}$ were used to analyze the effects.

Case 1: Both of the optical densities were constrained to be positive. In this case, as $DP{F}^{{\lambda }_2}$ was changed from the lower value of 3 to the higher value of 8, a translation was observed along the y-axis in the upward direction (Figure 4), but as $DP{F}^{{\lambda }_1}$ was increased more, the translation decreased. Additionally, the post-stimuli undershoot decreased as $DP{F}^{{\lambda }_2}$ increased. Similarly, a translation along the y-axis in the downward direction was observed for the variation of $DP{F}^{{\lambda }_2}$ (Figure 8). The change in $DP{F}^{{\lambda }_2}$ did not have much affect on the post-stimulus undershoot. Case 2: Both of the optical densities were constrained to be negative. This case showed a translation along the y-axis in the downward direction as $DP{F}^{{\lambda }_1}$ as varied. The translation step became small with proximity to the actual HRF (Figure 5). There is no effects on the post-stimulus undershoot. A translation also was observed along the y-axis in the upward direction as $DP{F}^{{\lambda }_2}$ as changed from lower to higher values (Figure 9). The post-stimulus undershoot decreased with increasing $DP{F}^{{\lambda }_2}$, and the translation became too small as value $DP{F}^{{\lambda }_2}$ is increased. Case 3: The optical density related to wavelength λ₁ was constrained to be positive, and the one related to wavelength λ₂ as negative. In this case, there was no effect on the positive part of the HRF (Figures 6, 10) for variation in either $DP{F}^{{\lambda }_1}$ or $DP{F}^{{\lambda }_2}$. The post-stimulus undershoot decreased as $DP{F}^{{\lambda }_2}$ increased, and the increment of decrease became small as the $DP{F}^{{\lambda }_2}$ values increased (Figure 10). Similarly, a slight change in post-stimulus undershoot was observed as $DP{F}^{{\lambda }_1}$ increased (Figure 6). Case 4: The optical density related to wavelength λ₁ was constrained to be negative, and the one related to wavelength λ₂ as positive. In this case, no effect was observed on post-stimulus undershoot for either $DP{F}^{{\lambda }_1}$ or $DP{F}^{{\lambda }_2}$ variation. One of the main reasons for this is that negative output was not possible in this case. It was observed that the peak of the response decreased as $DP{F}^{{\lambda }_2}$ increased and that there was a slight change in the peak of HRF as $DP{F}^{{\lambda }_1}$ increased (Figures 7, 11).

FIGURE 5

Figure 5. Results for both negative cases (Case: 2) with variation in value of $DP{F}^{{\lambda }_1}$ (3–8) under stimulation St₁-St₅ (top left, top middle, top right, bottom left, and bottom middle).

FIGURE 6

Figure 6. Results for first positive and second negative cases (Case: 3) with variation in value of $DP{F}^{{\lambda }_1}$ (3–8) under stimulation St₁-St₅ (top left, top middle, top right, bottom left, and bottom middle).

FIGURE 7

Figure 7. Results for first negative and second positive cases (Case: 4) with variation in value of $DP{F}^{{\lambda }_1}$ (3–8) under stimulation St₁-St₅ (top left, top middle, top right, bottom left, and bottom middle).

FIGURE 8

Figure 8. Results for both positive cases (Case: 1) with variation in value of $DP{F}^{{\lambda }_2}$ (3–8) under stimulation St₁-St₅ (top left, top middle, top right, bottom left, and bottom middle).

FIGURE 9

Figure 9. Results for both negative cases (Case: 2) with variation in value of $DP{F}^{{\lambda }_2}$ (3–8) under stimulation St₁-St₅ (top left, top middle, top right, bottom left, and bottom middle).

FIGURE 10

Figure 10. Results for first positive and second negative cases (Case: 3) with variation in value of $DP{F}^{{\lambda }_2}$ (3–8) under stimulation St₁-St₅ (top left, top middle, top right, bottom left, and bottom middle).

FIGURE 11

Figure 11. Results for first negative and second positive cases (Case: 4) with variation in value of $DP{F}^{{\lambda }_2}$ (3–8) under stimulation St₁-St₅ (top left, top middle, top right, bottom left, and bottom middle).

The real data on the healthy human subjects were collected using fNIRS in a visuo-motor experiment. The data was fed to Equation (16), and the best possible HRF was estimated using cost function J₂. The results (Figure 12) indicated that the peak of HRF increased as $DP{F}^{{\lambda }_1}$ was increased from lower to higher values. Also, there was an increment in the full width at half maximum (FWHM) of the HRF. Similarly, in contrast to the $DP{F}^{{\lambda }_1}$ variation, a decrement in the peak of HRF was observed as $DP{F}^{{\lambda }_2}$ was increased, resulting also in a decrement of FWHM (Figure 13). The changes in the peak values of HRF with variation of $DP{F}^{{\lambda }_1}$ and $DP{F}^{{\lambda }_2}$ are plotted in Figure 14. The results clearly show that a change in DPF corresponding to any wavelength significantly effects HRF shape. It is well known fact that in most BCI-fNIRS applications, features are observed from HRF attributes. Thus, a slight change in HRF on the basis of the DPF can lead to incurred error, misleading results, and reduced accuracy.

FIGURE 12

Figure 12. Results for fNIRS data set with variation in $DP{F}^{{\lambda }_1}$.

FIGURE 13

Figure 13. Results for fNIRS data set with variation in $DP{F}^{{\lambda }_2}$.

FIGURE 14

Figure 14. Effects of variations in $DP{F}^{{\lambda }_1}$ and $DP{F}^{{\lambda }_2}$ on peak of HRF for all subjects.

It was evident from the diffusion equation relating to the transportation of light through the homogeneous semi-infinite medium that the DPF depends on the reduced scattering coefficient, the absorption coefficient, and the source-detector separation (Scholkmann and Wolf, 2013). The human brain tissue causes high scattering of NIR light. The difference in structure and behavior of tissue is also possible among different subjects. Therefore, each human being respond differently to NIR light. That is due to the variation in the optical properties of the path-length of light photons. It also causes the change in the traveled distance of light photons. Talukdar et al. (2013) concluded that multiple layer of human scalp, skull, and brain tissues can effect scattering and thus DPF, which leads to systematic errors. Thus, the variation of DPF among the subjects is obvious. Similarly, the MBLL can interpret the dependency of HRF on DPF (Maikala, 2010). Therefore, it is worthwhile to analyze the effects of DPF change on NIR light based on human brain signals, particularly given that different brain areas have different scattering and absorption properties according to a subject's specific physiology and condition. A possible means of analyzing the effect of DPF on HRF is to split the HRF into three sections, namely pre-stimuli, main response peak, and post-stimuli undershoot. Therefore, a generic overview of the effects of DPF on HRF was presented in this paper.

Conclusion

In this study, the effects of the DPF were analyzed for fNIRS-observed data. The observed optical densities were utilized to formulate an iterative optimization problem for estimation of the best possible HRFs for different values of DPF related to wavelengths λ₁ and λ₂. Different simulated data sets were generated for a combination of different stimuli and activation times. The simulated HRF was fed into an optimization problem in order to split it into two different optical densities for further analyses. Later, different values of DPF were used to regenerate actual data and its effects. In addition to the simulation of data, the motor cortices of six healthy human subjects were scanned for a DPF analysis of measured fNIRS signals. It was concluded that the DPF affects the attributes of the HRF and that correct values of DPF are required for accurate results and their proper interpretation.

Author Contributions

MAK has proposed the methodology for analysis, did simulations and wrote the first draft of the manuscript. MM has done literature survey, implementation of methodology, did simulations and helped in editing of manuscript before submission. MJ has supervised the research.

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 National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIP) (No. 2015R1A5A1037668).

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

Cavazza, A., Aranyi, G., and Charles, F. (2017). BCI control of heuristic search algorithms. Front. Neuroinform. 11:6. doi: 10.3389/fninf.2017.00006

CrossRef Full Text | Google Scholar

Chatterjee, S., Phillips, J. P., and Kyriacou, P. A. (2015). “Differential pathlength factor estimation for brain-like tissue from a single-layer Monte Carlo model,” in Engineering in Medicine and Biology Society (EMBC), 2015 37th Annual International Conference of the IEEE (Milan), 3279–3282.

Google Scholar

Delpy, D. T., Cope, M., van der Zee, P., Arridge, S. R., Wray, S., and Wyatt, J. S. (1988). Estimation of optical pathlength through tissue from direct time of flight measurements. Phys. Med. Biol. 33, 1433–1442. doi: 10.1088/0031-9155/33/12/008

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. doi: 10.1203/00006450-199605000-00025

PubMed Abstract | CrossRef Full Text | Google Scholar

Duncan, A., Meek, J. H., Clemence, M., Elwell, C. E., Tyszczuk, L., Cope, M., et al. (1995). Optical pathlength measurements on adult head, calf and forearm and the head of newborn infants using phase resolved optical spectroscopy. Phys. Med. Biol. 40, 295–304. doi: 10.1088/0031-9155/40/2/007

CrossRef Full Text | Google Scholar

Duncan, A., Whitlock, T. L., Cope, M., and Delpy, D. T. (1993). A multiwavelength wideband intensity modulated optical spectrometer for near-infrared spectroscopy and imaging. Proc. SPIE 1888, 248–257. doi: 10.1117/12.154641

CrossRef Full Text | Google Scholar

Essenpreis, M., Elwell, C. E., Cope, M., van der Zee, P., Arridge, S. R., and Delpy, D. T. (1993). Spectral dependence of temporal point spread functions in human tissue. Appl. Opt. 32, 418–424. doi: 10.1364/AO.32.000418

CrossRef Full Text | Google Scholar

Ferrari, M., Wei, Q., Carraresi, L., De Blasi, R., and Zaceanti, G. (1992). Time resolved spectroscopy of human forearm. J. Photochem. Photobid. Biol. 16, 141–153. doi: 10.1016/1011-1344(92)80005-G

PubMed Abstract | CrossRef Full Text | 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. doi: 10.1006/nimg.1997.0306

CrossRef Full Text

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 Map. 2, 189–210. doi: 10.1002/hbm.460020402

CrossRef Full Text | Google Scholar

García-Prieto, J., Bajo, R., and Pereda, E. (2017). Efficient computation of functional brain networks: toward real-time functional connectivity. Front. Hum. Neurosci.11:8. doi: 10.3389/fninf.2017.00008

CrossRef Full Text | Google Scholar

Haftka, R. T., Gurdal, Z., and Kamat, M. P. (1990). Element of Structural Optimization. Dordrecht: Kluwer Academic Publisher.

Hiraoka, M., Firbank, M., Essenpreis, M., Cope, M., Arridge, S., Van Der Zee, P., et al. (1993). A Monte Carlo investigation of optical pathlength in inhomogeneous tissue and its application to near-infrared spectroscopy. Phys. Med. Biol. 38, 1859–1876. doi: 10.1088/0031-9155/38/12/011

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

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

Kamran, M. A., Jeong, M.-J., and Mannan, M. M. N. (2015). Optimal hemodynamic response model for functional near-infrared spectroscopy. Front. Behav. Neurosci. 9:151. doi: 10.3389/fnbeh.2015.00151

PubMed Abstract | CrossRef Full Text | Google Scholar

Kamran, M. A., Mannan, M. M., and Jeong, M. Y. (2016). Cortical signal analysis and advances in functional near-infrared spectroscopy signal: a review. Front. Hum. Neurosci. 10:261. doi: 10.3389/fnhum.2016.00261

CrossRef Full Text | Google Scholar

Kohl, M., Nolte, C., Heekeren, H. R., Horst, S., Scholz, U., Obrig, H., et al. (1998). Determination of the wavelength dependence of the differential pathlength factor from near-infrared pulse signals. Phys. Med. Biol. 43, 1771–1782. doi: 10.1088/0031-9155/43/6/028.

PubMed Abstract | CrossRef Full Text | Google Scholar

Koray, C., Bulent, S., Yasemin, P. K., and Ata, 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

CrossRef Full Text | Google Scholar

Lackowicz, I. R., and Berndt, K. (1990). Frequency-domain measurements of photon migration in tissue. Chem. Phys. Lett. 166, 246–225. doi: 10.1016/0009-2614(90)80024-8

CrossRef Full Text

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Maikala, R. (2010). Modified Beer's law – historical perspectives and relevance in near-infrared monitoring of optical properties of human tissue. Int. J. Ind. Ergon. 40, 125–134. doi: 10.1016/j.ergon.2009.02.011

CrossRef Full Text | Google Scholar

Mannan, M. M. N., Jeong, M. Y., and Kamran, M. A. (2016a). Hybrid ICA- regression: automatic identification and removal of ocular artifacts from electroencephalographic signals. Front. Hum. Neurosci. 10:193. doi: 10.3389/fnhum.2016.00193

CrossRef Full Text | Google Scholar

Mannan, M. M. N., Kim, S., Jeong, M. Y., and Kamran, M. A. (2016b). Hybrid EEG-eye tracker: automatic identification and removal of eye movement and blink artifacts from electroencephalographic signal. Sensors 16:241. doi: 10.3390/s16020241

CrossRef Full Text | Google Scholar

Meszlényi, R. J., Buza, K., and Vidnyánszky, Z. (2017). Resting State fMRI functional connectivity-based classification using a convolutional neural network architecture. Front. Hum. Neurosci. 11:61. doi: 10.3389/fninf.2017.00061

CrossRef Full Text | Google Scholar

Naseer, N., and Hong, K.-S. (2015). fNIRS-based brain-computer interfaces: a review. Front. Hum. Neurosci. 9:3. doi: 10.3389/fnhum.2015.00003

CrossRef Full Text | Google Scholar

Piao, D. D., Barbour, R., Graber, H., and Lee, D. (2015). “How does the differential pathlength factor for steady-state near-infrared spectroscopy of homogeneous medium vary with geometry?,” in Optical Trapping Applications: Optical Society of America (Vancouver, BC).

Google Scholar

Prince, S. J. D., 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/00319155/48/11/301.

PubMed Abstract | CrossRef Full Text | Google Scholar

Scholkmann, F., and Wolf, M. (2013). General equation for the differential pathlength factor of the frontal head depending on wavelength and age. J. Biomed. Opt. 18:105004. doi: 10.1117/1.JBO.18.10.105004

CrossRef Full Text | Google Scholar

Schroeter, M. L., Bucheler, M. M., Muller, K., Uludag, K., Villringer, A., and Cramon, D. Y. (2004). Towards standard analysis of functional near-infrared imaging. Neuroimage 21, 283–290. doi: 10.1016/j.neuroimage.2003.09.054

PubMed Abstract | CrossRef Full Text | Google Scholar

Schroeter, M. L., Zysset, S., Kruggel, F., and Cramon, D. Y. (2003). Age dependency of the hemodynamic response as measured by functional near-infrared spectroscopy. Neuroimage 19, 555–564. doi: 10.1016/S1053-8119(03)00155-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Talukdar, T., Moore, J. H., and Diamond, S. G. (2013). Contineous correction of differential path length factor in near-infrared spectroscopy. J. Biomed. Opt. 18:56001. doi: 10.1117/1.JBO.18.5.056001

CrossRef Full Text

van der Zee, P., Arridge, S. R., Cope, M., and Delpy, D. T. (1990). The effect of optode positioning on optical pathlength in near infrared spectroscopy of brain. Adv. Exp. Med. Bid. 227, 79–84. doi: 10.1007/978-1-4684-8181-5_11

CrossRef Full Text | Google Scholar

van der Zee, P., Cope, M., Arridge, S. R., Essenpreis, M., Potter, L. A., Edwards, A. D., et al. (1992). Experimentally measured optical pathlengths for the adult head, calf and forearm and the head of the newborn infant as a function of inter optode spacing. Adv. Exp. Med. Biol. 316, 143–153. doi: 10.1007/978-1-4615-3404-4_17

PubMed Abstract | 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

Zhao, H., Tanikawa, Y., Gao, F., Onodera, Y., Sassaroli, A., Tanaka, K., et al. (2002). Maps of optical differential pathlength factor of human adult forehead, somatosensory motor and occipital regions at multi-wavelength in NIR. Phys. Med. Biol. 47, 2075–2093. doi: 10.1088/0031-9155/47/12/306

CrossRef Full Text | Google Scholar