1.

Introduction

Functional near-infrared spectroscopy (fNIRS) is a functional neuroimaging technique similar to functional magnetic resonance imaging (fMRI) and positron emission tomography (PET). It detects cerebral functional hemodynamic changes, but not the cerebral neural activation itself. fNIRS potentially has some advantages over other techniques: portability, inexpensive equipment, and availability for simultaneous use with other measurement devices including fMRI, PET, electroencephalography (EEG), and magnetoencephalography (MEG). However, undesirable artifacts are often observed in fNIRS measurements. Systemic physiological changes, physical movements, and even posture changes in a subject can generate such artifacts. These artifacts are observed because the detected light is influenced not only by the hemodynamic change in gray matter, but also by scattering and absorption changes in all tissues through which the light propagates. It is not known exactly how scattering and absorption changes occur in each tissue with a subject’s physiological change. Thus, various statistical methods, such as conventional filtering and multivariate data analysis, have been used to remove such artifacts from fNIRS measurements.

Artifacts induced by physiological activity, such as cardiac pulsation, respiration, and a blood pressure change, can be effectively removed by low-pass filtering or a moving average method. These techniques work well in a real-time environment because they are very simple to implement. However, they tend to overflatten the signal. Thus, for example, a fast-rising signal of brain activation may be deleted if they are applied. To solve this problem, Zhang, Brown, and Strangman^{1, 2} proposed an adaptive filtering technique based on low correlation between artifacts and cerebral functional activity. Saager and Berger³ applied the least square fit technique to this problem and conducted a detailed analysis of its performance. Kohno ⁴ used a powerful technique, independent component analysis, to separate a statistically independent signal from data contaminated with various artifacts.

These techniques and a conventional block averaging technique are useful for suppressing task-uncorrelated noises. However, if the data are contaminated with an artifact that is strongly correlated with a task sequence, these methods cannot work effectively. For example, Izzetoglu ⁵ pointed out that a subject’s head movement generates a blood flow change; thus, if the subject moves his head synchronously with a task, an artifact induced by this blood flow change will be observed. Since this kind of artifact is tightly correlated with the task sequence, the block averaging technique cannot remove it.

In general, to remove task-correlated artifacts, we adopt an experimental design including a reference task that is considered to generate the same signal as a target task except for a functional component. The artifacts can be canceled, but the functional component remains when we calculate the difference between the target and reference measurements. Such a design, however, usually requires a longer execution time, and it is very difficult to verify the equivalence of artifact components in target and reference tasks.

Another approach to this problem is multivariate data analysis. This approach is based on the empirical knowledge that the artifact signal is observed globally, and the functional signal is observed locally. Principal component analysis is used to extract the functional signal from multichannel measurement data. This method was reported to be effective.^{6, 7} However, the use of many optodes, most of which have to be positioned at areas unrelated to the target, is a big problem for practical measurement. In addition, if several activations at different areas occur simultaneously, separation of these signals will be difficult.

We do not know precisely how task-correlated artifacts are generated. Presumably they originate mainly from something other than activity in the cortical layer. Thus, it may be useful for functional detection if we can separate the signal of the cortical layer from that of other layers. Hence, several methods using multidistance probe arrangements for separating NIRS signals in lower layers from other layers have been proposed.

A NIRS algorithm using a multidistance probe arrangement, NIR spatially resolved spectroscopy (NIR-SRS), is known.^{8, 9} NIR-SRS was originally developed to measure the tissue oxygen index (TOI), and Al-Rawi, Smielewski, and Kirkpatrick reported that an extracarotid artery clamp had little influence on the TOI measured by NIR-SRS.¹⁰ This result may suggest that NIR-SRS can remove the interference of extracranial blood flow from brain functional activity signals. However, we cannot theoretically validate its ability to reduce artifacts of a superficial layer, because NIR-SRS is based on the theory of light propagation in a monolayered medium.¹¹

Improved approaches of multidistance probe arrangement based on two-layer models were reported. Choi ¹² used an analytical solution of a two-layer model that contains an additional superficial layer on a basal semi-infinite medium. They obtained the hemoglobin oxygenation etc. of extra- and intracranial layers by a theoretical curve fitting to experimental data. Fabbri ¹³ conducted Monte Carlo simulations of a two-layered tissue model and showed that a multidistance probe arrangement is effective in detecting the absorption change in the lower layer (cortical layer), even when optical changes arise in both lower and upper layers. This is a viable technique for separating the signal at the cortical layer from artifacts, regardless of whether they are task correlated or uncorrelated.

Validity of these algorithms greatly depends on the optical models assumed. It should be examined by theoretical or simulation approaches whether such a simplified model of a few layers provides a sufficient condition to separate cortical signals from others. In addition, only absorption change in tissues has been considered in most algorithms based on the modified Lambert-Beer law. In a recent study, Tomita, Ohtomo, and Suzuki reported that hemodynamic changes could be accompanied by scattering changes caused by blood cell aggregation.¹⁴ Since scattering also causes optical attenuation changes, we think that analysis considering both absorption and scattering is important.

In this work, we present a theoretical analysis and Monte Carlo simulations of a five-layered tissue model in which both scattering and absorption changes occur, and show that a multidistance probe arrangement is effective in extracting absorption changes in the gray matter layer. The simulation result is used to determine the appropriate probe arrangement. The proposed method is compared with a conventional method through actual experiments on nonfunctional and functional tasks. Artifacts induced by body tilting, head nodding, and breath holding were clearly observed when a conventional method was used, and were appreciably reduced by the proposed method. Signals evoked by single-sided finger movements were observed globally at both hemispheres when we used a conventional method. On the other hand, localized signals at the primary motor area were observed by the proposed method. A statistically significant increase in oxygenated hemoglobin (HbO) and decrease in deoxygenated hemoglobin (HbR) were simultaneously observed only at each contralateral primary motor area.

2.

Theory

2.1.

Optical Multilayer Model

A typical probe arrangement of fNIRS is illustrated in Fig. 1 (source and detector 1). In this arrangement, light power is detectable only when scattering exists in a medium. Light from a source of power $I_{0,\lambda }$ propagates in a medium according to the spatial distribution of the scattering and absorption coefficients ( ${\mu }_{s,\lambda }^{\prime }$ and ${\mu }_{a,\lambda }$ ), and reaches a detector at a distance $d_1$ with power $I_{\lambda }$ . Here, $\lambda$ indicates the wavelength of the incident light. The optical attenuation is defined as $A_{d,\lambda }=-\ln (I_{\lambda }∕I_{0,\lambda })$ . Since the input power $I_{0,\lambda }$ is difficult to measure, the optical attenuation change from a certain standard condition is used in an ordinary fNIRS measurement. When the change in scattering and absorption coefficients ( $\Delta {\mu }_{s,\lambda }^{\prime }$ and $\Delta {\mu }_{a,\lambda }$ ) from the standard condition is small, the optical attenuation change $\Delta A_{d,\lambda }$ can be represented as

Eq. 1

$\Delta A_{d,\lambda }={\phantom{|}\frac{\partial A_{d,\lambda }}{\partial {\mu }_a}|}_{\mathrm{sv}}\Delta {\mu }_{a,\lambda }+{\phantom{|}\frac{\partial A_{d,\lambda }}{\partial {\mu }_s^{\prime }}|}_{\mathrm{sv}}\Delta {\mu }_{s,\lambda }^{\prime }.$

Here, ${\phantom{|}X|}_{\mathrm{sv}}$ denotes the value of $X$ at the standard condition. The coefficients ${\phantom{|}\partial A_{d,\lambda }∕\partial {\mu }_a|}_{\mathrm{sv}}$ and ${\phantom{|}\partial A_{d,\lambda }∕\partial {\mu }_s^{\prime }|}_{\mathrm{sv}}$ have the dimension of length.

Fig. 1

Schematic illustration of multidistance probe arrangement on the human head.

If the second term of Eq. 1 is negligible, this equation corresponds to the modified Lambert-Beer law, which is the theoretical basis of a conventional fNIRS algorithm. In this case, the coefficient ${\phantom{|}\partial A_{d,\lambda }∕\partial {\mu }_a|}_{\mathrm{sv}}$ is called the mean optical path length. In general, however, the second term cannot be neglected. For instance, Tomita, Ohtomo, and Suzuki experimentally demonstrated in vitro that the scattering property of a blood volume changes when its flow velocity changes.¹⁴ Thus, in this work we set the starting point of our theoretical discussion at Eq. 1.

We use a five-layered optical slab model of an adult head, which consists of layers of scalp (sc), skull (sk), cerebrospinal fluid (csf), gray matter (gm), and white matter (wm). Supposing the temporal change in scattering and absorption to be uniform in each layer, Eq. 1 is transformed to

Eq. 2

$\Delta A_{d,\lambda }=\sum _L{l}_{d,\lambda }^L\Delta {\mu }_{a,\lambda }^L+\sum _L{m}_{d,\lambda }^L\Delta {\mu }_{s,\lambda }^{\prime L},$

where $v^L$ denotes the parameter $v$ at the layer $L$ . The summations are accomplished with $L∊\{\mathrm{sc},\mathrm{sk},\mathrm{csf},\mathrm{gm},\mathrm{wm}\}$ , and $l_{d,\lambda }^L$ and $m_{d,\lambda }^L$ denote

Eq. 3

$\{\begin{array}{c}{l}_{d,\lambda }^L={\phantom{|}\frac{\partial A_{d,\lambda }}{\partial {\mu }_a^L}|}_{\mathrm{sv}}\\ m_{d,\lambda }^L={\phantom{|}\frac{\partial A_{d,\lambda }}{\partial {\mu }_s^{\prime L}}|}_{\mathrm{sv}}\end{array}\phantom{\}}.$

Instead of the analytical derivation of the coefficients $l_{d,\lambda }^L$ and $m_{d,\lambda }^L$ , we numerically obtained these values by a Monte Carlo simulation. The coefficients are calculated as follows:

Eq. 4

$\{\begin{array}{l}{l}_{d,\lambda }^L=\{A_{d,\lambda }({\overline{\mu }}_a^L+{\mu }_a^L,{\overline{\mathbf{\mu }}}_o)-A_{d,\lambda }({\overline{\mu }}_a^L,{\overline{\mathbf{\mu }}}_o)\}∕{\mu }_a^L\\ m_{d,\lambda }^L=\{A_{d,\lambda }({\overline{\mu }}_s^{\prime L}+{\mu }_s^{\prime L},{\overline{\mathbf{\mu }}}_o^{\prime })-A_{d,\lambda }({\overline{\mu }}_s^{\prime L},{\overline{\mathbf{\mu }}}_o^{\prime })\}∕{\mu }_s^{\prime L}\end{array}\phantom{\}},$

where ${\mathbf{\mu }}_o$ represents all parameters except ${\mu }_a^L$ , and ${\mathbf{\mu }}_o^{\prime }$ represents all parameters except ${\mu }_s^{\prime L}$ . The bar symbol over a parameter $\mu$ $\left(\overline{\mu }\right)$ denotes the standard condition value. Optical attenuations such as $A_{d,\lambda }({\overline{\mu }}_a^L+{\mu }_a^L,{\overline{\mathbf{\mu }}}_o)$ and $A_{d,\lambda }({\overline{\mu }}_s^{\prime L}+{\mu }_s^{\prime L},{\overline{\mathbf{\mu }}}_o^{\prime })$ were obtained by a Monte Carlo simulation, where ${\mu }_a^L$ and ${\mu }_s^{\prime L}$ were set to satisfy ${\mu }_a^L∕{\overline{\mu }}_a^L=0.1$ and ${\mu }_s^{\prime L}∕{\overline{\mu }}_s^{\prime L}=0.1$ , respectively. Parameter values at the standard condition and layer thickness were based on the literature (Okui and Okada¹⁵) and are shown in Table 1 . The program tMCimg, provided to the public by the Photon Migration Imaging Laboratory at Massachusetts General Hospital (see http://www.nmr.mgh.harvard.edu/PMI/index.htm), was used in the simulation.

Table 1

Hemodynamic parameters, thickness, and optical properties for each layer of an adult head model.

Tissue type	Thickness[mm]	μa [mm−1] at 670∕690∕720∕750∕780∕830 [nm]
Scalp	3	$0.031∕0.026∕0.022∕0.022∕0.020∕0.019$
Skull	7	$0.027∕0.017∕0.013∕0.016∕0.016∕0.017$
CSF	2	$0.0044∕0.0029∕0.003∕0.0045∕0.0044∕0.0056$
Gray matter	4	$0.048∕0.039∕0.036∕0.038∕0.036∕0.041$
White matter	34	$0.024∕0.018∕0.014∕0.016∕0.016∕0.018$
Tissue type		${\mu }_s^{\prime }$ $\left[{\mathrm{mm}}^{-1}\right]$ at $670∕690∕720∕750∕780∕830$ [nm]
Scalp		$2.53∕2.38∕2.24∕2.11∕2.00∕1.84$
Skull		$2.33∕2.13∕1.92∕1.79∕1.66∕1.47$
CSF		$0.35∕0.32∕0.29∕0.27∕0.25∕0.22$
Gray matter		$2.69∕2.57∕2.46∕2.39∕2.31∕2.10$
White matter		$10.49∕10.03∕9.78∕9.62∕9.25∕8.82$

Figure 2 shows a 2-D plot of optical attenuation changes when parameters ${\mu }_a^L$ and ${\mu }_s^{\prime L}$ increase by 10% from standard values in each layer $L$ . Five source-detector distances (10, 15, 20, 25, and $30\mathrm{mm}$ ) were assumed, and an incident light of wavelength $780\mathrm{nm}$ was used in the simulation. This figure shows the following.

• • Scattering and absorption changes at the scalp and skull layers strongly affect the optical attenuation change. Interestingly, the attenuation change increases (skull layer) and decreases (scalp layer) proportionally if we increase the source-detector distance.

• • The optical attenuation change induced by an absorption change in the csf layer is small, and the attenuation change due to a scattering change is comparable to that due to an absorption change in the gray matter layer. The observed attenuation change in the csf layer is thought to be negligible because that layer contains no blood cell, and thus scattering and absorption changes due to blood flow changes should be very small.

• • The absorption change in the gray matter layer influences the attenuation change considerably. This change increases proportionally with the source-detector distance. Remarkably, the optical attenuation change caused by a scattering change in the gray matter layer is very small.

• • The optical attenuation change induced by absorption and scattering changes in the white matter layer is very small.

Consequently, we assume that both the scattering and absorption terms of the csf and white matter layers and the scattering term of the gray matter layer barely contribute to the optical attenuation change. We omit these terms from Eq. 2.

Eq. 5

$\Delta A_{d,\lambda }=(l_{d,\lambda }^{\mathrm{sc}}\Delta {\mu }_{a,\lambda }^{\mathrm{sc}}+m_{d,\lambda }^{\mathrm{sc}}\Delta {\mu }_{s,\lambda }^{\prime \mathrm{sc}})+(l_{d,\lambda }^{\mathrm{sk}}\Delta {\mu }_{a,\lambda }^{\mathrm{sk}}+m_{d,\lambda }^{\mathrm{sk}}\Delta {\mu }_{s,\lambda }^{\prime \mathrm{sk}})+l_{d,\lambda }^{\mathrm{gm}}\Delta {\mu }_{a,\lambda }^{\mathrm{gm}}.$

Concerning the skull layer, diploic veins are found in an area of cancellous bone between the inner and outer compact bone layers in the skull, and blood is drained in this area. However, we do not exactly know how the blood in the skull layer changes its volume and flow with systemic activity changes. Here, we assume two contrasting conditions of minimally and maximally changing hemodynamics in the skull layer. In a minimal condition, hemodynamic change in the skull layer can be negligible by supposing that blood flow in the skull layer is barely affected by systemic activity changes. In this case, we assume $\Delta {\mu }_{a,\lambda }^{\mathrm{sk}}=0$ and $\Delta {\mu }_{s,\lambda }^{\prime \mathrm{sk}}=0$ , thus, the skull term in Eq. 5 is omitted and the scalp term remains as a single superficial layer term. In contrast, we suppose a maximal condition that hemodynamic change in the skull is very similar to that in the scalp layer. In this case, we assume $\Delta {\mu }_{a,\lambda }^{\mathrm{sk}}=\Delta {\mu }_{a,\lambda }^{\mathrm{sc}}$ and $\Delta {\mu }_{s,\lambda }^{\prime \mathrm{sk}}=\Delta {\mu }_{s,\lambda }^{\prime \mathrm{sc}}$ , and scalp and skull terms in Eq. 5 are merged into a single superficial layer term. In either condition, Eq. 5 is simplified,

Eq. 6

$\Delta A_{d,\lambda }=l_{d,\lambda }^{\mathrm{sp}}\Delta {\mu }_{a,\lambda }^{\mathrm{sp}}+m_{d,\lambda }^{\mathrm{sp}}\Delta {\mu }_{s^{\prime },\lambda }^{\mathrm{sp}}+l_{d,\lambda }^{\mathrm{gm}}\Delta {\mu }_{a,\lambda }^{\mathrm{gm}},$

where

Eq. 7

$\{\begin{array}{l}{l}_{d,\lambda }^{\mathrm{sp}}=l_{d,\lambda }^{\mathrm{sc}}\\ m_{d,\lambda }^{\mathrm{sp}}=m_{d,\lambda }^{\mathrm{sc}}\end{array}\phantom{\}},\text{if}\{\begin{array}{c}\Delta {\mu }_{a,\lambda }^{\mathrm{sk}}=0\\ \Delta {\mu }_{s,\lambda }^{\prime \mathrm{sk}}=0\end{array}\phantom{\}},$

and

Eq. 8

$\{\begin{array}{l}{l}_{d,\lambda }^{\mathrm{sp}}=l_{d,\lambda }^{\mathrm{sc}}+l_{d,\lambda }^{\mathrm{sk}}\\ m_{d,\lambda }^{\mathrm{sp}}=m_{d,\lambda }^{\mathrm{sc}}+m_{d,\lambda }^{\mathrm{sk}}\end{array}\phantom{\}},\text{if}\{\begin{array}{c}\Delta {\mu }_{a,\lambda }^{\mathrm{sk}}=\Delta {\mu }_{a,\lambda }^{\mathrm{sc}}\\ \Delta {\mu }_{s,\lambda }^{\prime \mathrm{sk}}=\Delta {\mu }_{s,\lambda }^{\prime \mathrm{sc}}\end{array}\phantom{\}}.$

Here, sp denotes the superficial layer. As a result, the five-layer model was simplified to a two-layer model.

Fig. 2

2-D plot of optical attenuation changes where the parameters ${\mu }_a^L$ and ${\mu }_s^{\prime L}$ increase by 10% from the standard values in each layer $L∊\{\mathrm{sc},\mathrm{sk},\mathrm{csf},\mathrm{gm},\mathrm{wm}\}$ . Five source-detector distances (10, 15, 20, 25, and $30\mathrm{mm}$ ) were assumed, and an incident light of wavelength $780\mathrm{nm}$ was used in the simulation.

2.2.

Artifact Cancellation by Multidistance Probe Arrangement

If we use two detectors at different source-detector distances $d_1$ and $d_2$ $(d_1>d_2)$ , the changes in optical attenuation are written as follows.

Eq. 9

$\{\begin{array}{c}\Delta A_{d_1,\lambda }=l_{d_1,\lambda }^{\mathrm{sp}}\Delta {\mu }_{a,\lambda }^{\mathrm{sp}}+m_{d_1,\lambda }^{\mathrm{sp}}\Delta {\mu }_{s,\lambda }^{\prime \mathrm{sp}}+l_{d_1,\lambda }^{\mathrm{gm}}\Delta {\mu }_{a,\lambda }^{\mathrm{gm}}\\ \Delta A_{d_2,\lambda }=l_{d_2,\lambda }^{\mathrm{sp}}\Delta {\mu }_{a,\lambda }^{\mathrm{sp}}+m_{d_2,\lambda }^{\mathrm{sp}}\Delta {\mu }_{s,\lambda }^{\prime \mathrm{sp}}+l_{d_2,\lambda }^{\mathrm{gm}}\Delta {\mu }_{a,\lambda }^{\mathrm{gm}}\end{array}\phantom{\}}.$

Here, we introduce the ratios $k_{a,\lambda }$ and $k_{s,\lambda }$ as follows.

Eq. 10

$\{\begin{array}{l}{k}_{a,\lambda }=l_{d_1,\lambda }^{\mathrm{sp}}∕l_{d_2,\lambda }^{\mathrm{sp}}\\ k_{s,\lambda }=m_{d_1,\lambda }^{\mathrm{sp}}∕m_{d_2,\lambda }^{\mathrm{sp}}\end{array}\phantom{\}}.$

If detector distances can be arranged to satisfy $k_{a,\lambda }=k_{s,\lambda }$ , we can calculate the absorption change at the gray matter layer $\Delta {\mu }_{a,\lambda }^{\mathrm{gm}}$ as follows.

Eq. 11

$\Delta {\mu }_{a,\lambda }^{\mathrm{gm}}=\frac{\Delta A_{d_1,\lambda }-k_{\lambda }\Delta A_{d_2,\lambda }}{l_{d_1,\lambda }^{\mathrm{gm}}-k_{\lambda }{l}_{d_2,\lambda }^{\mathrm{gm}}},$

where $k_{\lambda }=k_{a,\lambda }=k_{s,\lambda }$ .

Figure 3 shows simulated values of $k_{a,\lambda }$ and $k_{s,\lambda }$ with various combinations of source-detector distances. In many cases of both the minimal and maximal conditions of the superficial layer, $k_{a,\lambda }$ and $k_{s,\lambda }$ had very similar values. However, when $d_2$ was $10\mathrm{mm}$ , the values did not coincide sufficiently.

Fig. 3

2-D plots of $k_{a,\lambda }$ and $k_{s,\lambda }$ at various combinations of source-detector distances in two contrasting conditions on the superficial layer. (a) Minimal condition ( $l_{d,\lambda }^{\mathrm{sp}}=l_{d,\lambda }^{\mathrm{sc}}$ , $m_{d,\lambda }^{\mathrm{sp}}=m_{d,\lambda }^{\mathrm{sc}}$ ). (b) Maximal condition ( $l_{d,\lambda }^{\mathrm{sp}}=l_{d,\lambda }^{\mathrm{sc}}+l_{d,\lambda }^{\mathrm{sk}}$ , $m_{d,\lambda }^{\mathrm{sp}}=m_{d,\lambda }^{\mathrm{sc}}+m_{d,\lambda }^{\mathrm{sk}}$ ). A wavelength of $780\mathrm{nm}$ was assumed.

We simulated $\Delta A_{d_1,\lambda }-k_{a,\lambda }\Delta A_{d_2,\lambda }$ for various combinations of detector distances. The absorption coefficient at each activation volume of the gray matter layer along the line connecting a source and a detector is assumed to increase by 10% from its standard value. In fMRI studies on human adult brains, functional activations are often observed at a volume of several tens of ${\mathrm{mm}}^3$ . On the other hand, if we set a smaller activation volume (larger photon numbers); a longer computational time will be required to obtain sufficient quality of the simulation result. From this practical reason, we fixed a $5\times 5\times 4{\mathrm{mm}}^3$ volume as an activation volume in this simulation. The result is shown in Fig. 4, which shows the sensitivity distribution for detecting a local absorption change in the gray matter layer. The sensitivities obtained under the minimal and maximal skull layer conditions are very similar to each other. These tend to become high when $d_1$ is large and $d_2$ is small. Thus, $d_1$ should be large and $d_2$ should be small as much as possible to realize good sensitivity. However, if the source-detector distance is increased, the measurement noise would become large, since fewer photons can reach the detector. In addition, if vasculature and hemodynamic activity in the tissues are not sufficiently homogeneous, a larger distance between two detectors may violate our fundamental assumption that each tissue layer is a uniform slab.

Fig. 4

The detection sensitivity of various combinations of source-detector distances, where the absorption coefficient at each $5\times 5\times 4{\mathrm{mm}}^3$ volume of the gray matter layer along the line connecting a source and a detector is assumed to increase by 10% from its standard value. In each frame, the simulated values of $\Delta A_{d_1,\lambda }-k_{a,\lambda }$ $\Delta A_{d_2,\lambda }$ at various positions of the activated volume were drawn by blue and red lines, which correspond with the minimal condition ( $l_{d,\lambda }^{\mathrm{sp}}=l_{d,\lambda }^{\mathrm{sc}}$ , $m_{d,\lambda }^{\mathrm{sp}}=m_{d,\lambda }^{\mathrm{sc}}$ ) and the maximal condition ( $l_{d,\lambda }^{\mathrm{sp}}=l_{d,\lambda }^{\mathrm{sc}}+l_{d,\lambda }^{\mathrm{sk}}$ , $m_{d,\lambda }^{\mathrm{sp}}=m_{d,\lambda }^{\mathrm{sc}}+m_{d,\lambda }^{\mathrm{sk}}$ ) of superficial layer, respectively. Green and two blue vertical lines indicate the positions of the source and two detectors, respectively. A wavelength of $780\mathrm{nm}$ was assumed.

Consequently, an effective probe arrangement should satisfy the following constraints.

• • The source-detector distance combination should realize the condition $k_{a,\lambda }=k_{s,\lambda }$ as much as possible.

• • $d_1$ should be large and $d_2$ should be small to realize good sensitivity to local absorption changes in the gray matter layer.

• • To keep the measurement noise small, $d_1$ should not be too large.

• • The difference between $d_1$ and $d_2$ should be small enough to satisfy the assumption of uniform tissue layers.

Considering Figs. 3 and 4, we chose the source-detector distance combination $d_1=30$ and $d_2=20\mathrm{mm}$ as a compromise meeting these conditions.

2.3.

Estimation of Partial Path Length Ratio $k$

To calculate the absorption change in the gray matter layer according to Eq. 11, we need to know $k_{\lambda }$ . We may be able to use a fixed value for $k_{\lambda }$ ,¹⁶ however, it is better to estimate it before each experiment. We may be able to obtain it as $k_{\lambda }=\Delta A_{d_1,\lambda }^{\text{rest}}∕\Delta A_{d_2,\lambda }^{\text{rest}}$ ,¹³ where $\Delta A_{d_1,\lambda }^{\text{rest}}$ and $\Delta A_{d_2,\lambda }^{\text{rest}}$ are optical attenuation changes at the resting state of a subject. Here, it is assumed that the absorption change in gray matter is small at the resting state. However, because the absolute value of $\Delta A_{d_2,\lambda }^{\text{rest}}$ could be very small, possibly zero, the ratio $\Delta A_{d_1,\lambda }^{\text{rest}}∕\Delta A_{d_2,\lambda }^{\text{rest}}$ can be very unstable, even when the measurement noise is small. We proposed another estimation method that is stable and gives more accurate values in noisy conditions. The details of the method are described in Appendix { label needed for app[@id='xA'] } in Sec. 6.

We use a triple wavelength measurement, and $k_{\lambda }$ is estimated for each wavelength. Monte Carlo simulations under both the minimal and maximal conditions of the superficial layer indicate that the wavelength dependence of $k_{a,\lambda }$ and $k_{s,\lambda }$ is weak when $d_1=30\mathrm{mm}$ and $d_2=20\mathrm{mm}$ (Fig. 5 ). To statistically stabilize the estimation, we used the average $k$ of three estimated values when using Eq. 11.

Fig. 5

2-D plot of wavelength dependence of $k_{a,\lambda }$ and $k_{s,\lambda }$ in two contrasting conditions on the superficial layer. Blue and red colored markers correspond with the minimal condition ( $l_{d,\lambda }^{\mathrm{sp}}=l_{d,\lambda }^{\mathrm{sc}}$ , $m_{d,\lambda }^{\mathrm{sp}}=m_{d,\lambda }^{\mathrm{sc}}$ ) and the maximal condition ( $l_{d,\lambda }^{\mathrm{sp}}=l_{d,\lambda }^{\mathrm{sc}}+l_{d,\lambda }^{\mathrm{sk}}$ , $m_{d,\lambda }^{\mathrm{sp}}=m_{d,\lambda }^{\mathrm{sc}}+m_{d,\lambda }^{\mathrm{sk}}$ ) of superficial layers, respectively. Source-detector distances of $d_1=30\mathrm{mm}$ and $d_2=20\mathrm{mm}$ were assumed.

3.

Methods

A healthy adult male volunteer participated in the experiment. This study was approved by the Institutional Review Board of National Institute of Advanced Industrial Science and Technology (Japan), and the participant gave his written informed consent.

For fNIRS measurement, we used a near-infrared oximeter with a multiprobe adapter system (NIRO-200 with C9866, Hamamatsu Photonics KK), which has two sources, each with wavelengths 776, 809, and $850\mathrm{nm}$ , and eight detectors. The sampling frequency was $0.5\mathrm{Hz}$ . Each recorded signal was high-pass filtered off-line at a cut-off frequency of $0.01\mathrm{Hz}$ . Based on the simulation results, we set distances from the source to the detectors at $d_1=30\mathrm{mm}$ and $d_2=20\mathrm{mm}$ . Using an originally designed holder, the two sources and eight detectors were fixed on the scalp as shown in Figs. 6e and 6f . The detectors were numbered as shown in Fig. 6c. We implemented the conventional method using detectors 1, 3, 5, and 7 and the proposed method using paired detectors 1–2, 3–4, 5–6, and 7–8.

Fig. 6

(a) Axial and (b) coronal sections of the functional volume image of left-hand finger opposition. (c) and (d) Intersections of solid lines give the projections of activated areas on the scalp surface. Paired detectors are located at the position where their middle point corresponds to the projected point. Markers numbered 1 to 8 indicate the detector positions. Markers S1 and S2 indicate source positions. (e) and (f) Photographs of the probe set on a subject model.

Paired detectors 1–2 and 7–8 were positioned on the activation area of the finger opposition task as follows. To locate the activation area during the task, T1-weighted images and echo planar images were measured using MR equipment (MR Signa 3.0T, GE Yokogawa Medical Systems KK). The subject was instructed to tap his thumb with his index finger at a frequency of $4\mathrm{Hz}$ . The block conditions were alternated using visual cues in this order: left-hand finger stimulation periods $\left(20\mathrm{s}\right)$ , rest periods $\left(20\mathrm{s}\right)$ , right-hand finger stimulation periods $\left(20\mathrm{s}\right)$ , and rest periods $\left(20\mathrm{s}\right)$ . A complete session consisted of an initial reference rest period $\left(20\mathrm{s}\right)$ and five trials of the block condition (rest, left finger task, rest, right finger task) without any interruption. T-contrast images of left- and right-hand finger opposition against rest periods were obtained using SPM5 (see http://www.fil.iion.ucl.ac.uk/spm/). The images of left-hand activity versus rest appear as colored pixels in Figs. 6a and 6b. The detector was located directly above the activation area, and the position of the activation area was identified by the overlying T1-weighted functional images. The Cz position of an EEG 10–20 system, which is defined as the midpoint of the medial sagittal surface curve connecting nasion and inion, was also identified using the T1-weighted image. The paired detectors 1–2 and 7–8 were located relative to the Cz position.

Optical attenuation changes at the resting state were measured for $10\min$ to estimate the ratio $k$ of each pair of detectors. Table 2 shows the correlation coefficients of the optical attenuation changes in each pair of detectors during the resting state. These values show that the attenuation changes at the paired detectors are highly correlated with each other. Based on the measured optical attenuation change, the $k$ value for each pair was calculated and is shown in the last column of Table 2.

Table 2

Correlation coefficients of ΔAd,λ between adjacent detectors and k values.

To verify the effect of the proposed method on the cancellation of artifacts caused by posture change, body motion, and respiration, the subject, sitting in a chair, was instructed to perform three tasks: tilting his upper body forward by about $30\mathrm{deg}$ , nodding his head by about $30\mathrm{deg}$ with a frequency of $0.25\mathrm{Hz}$ , and holding his breath. In each task, the block conditions were alternated through task periods $\left(20\mathrm{s}\right)$ and rest periods $\left(20\mathrm{s}\right)$ using auditory cues. A complete session consisted of an initial reference rest period $\left(20\mathrm{s}\right)$ followed by a block of five task/rest sequences.

The finger opposition task was also performed to verify that the proposed method can detect a brain functional signal. To evaluate the performance of the proposed and the conventional (using a single detector and the modified Lambert-Beer law) methods for detecting brain functional activity in the finger opposition task, we analyzed statistical significance by the following procedure. The total task sequence comprises five trials of rest-task pairs. We averaged $\Delta \mathrm{Hb}\mathrm{O}$ and $\Delta \mathrm{HbR}$ at each task and rest periods (each has a duration of $20\mathrm{s}$ ). We adopted a temporal offset of $6\mathrm{s}$ for each period to take into account transient phases of hemodynamic change. Then we calculated the difference between the rest and task averages of five trials. We finally applied the paired $t$ -test to the resultant five differences.

4.

Results

The conventional and proposed methods are compared for cases of body tilting, head nodding, and breath holding in Fig. 7 . $\Delta \mathrm{Hb}\mathrm{O}$ and $\Delta \mathrm{HbR}$ of five task/rest sequences are block-averaged. In all cases using the conventional method, $\Delta \mathrm{Hb}\mathrm{O}$ and $\Delta \mathrm{HbR}$ significantly changed during the task. The magnitude of these changes was comparable to or larger than that of the finger opposition task (shown in Fig. 8 ). The proposed method effectively reduced the magnitude of the changes.

Fig. 7

Block averages and standard deviations of hemoglobin changes during three tasks: body tilting, head nodding, and breath holding. Red, blue, and green lines show $\Delta \mathrm{Hb}\mathrm{O}$ , $\Delta \mathrm{HbR}$ , and the task period, respectively. For each task, frames on left and right sides correspond to the conventional and proposed methods, respectively. Four frames in each column correspond to detector positions 1, 3, 5, and 7 for the conventional method and 1–2, 3–4, 5–6, and 7–8 for the proposed method.

Fig. 8

Hemoglobin changes during the finger opposition task. Red, blue, and green lines show $\Delta \mathrm{Hb}\mathrm{O}$ , and $\Delta \mathrm{HbR}$ , and the task period, respectively. The two left and two right columns show block-averaged hemoglobin changes and their standard deviations for the conventional and proposed methods, respectively. In each method, the left and right frames correspond to left and right finger oppositions, respectively. Four frames in each column correspond to detector positions 1, 3, 5, and 7 for the conventional method, and 1–2, 3–4, 5–6, and 7–8 for the proposed method.

The observed $\Delta \mathrm{Hb}\mathrm{O}$ and $\Delta \mathrm{HbR}$ for the finger opposition task using the two methods are shown in Fig. 8. In the conventional method, increases in $\Delta \mathrm{Hb}\mathrm{O}$ were observed at almost all detectors regardless of the side on which finger opposition was performed. The changes in $\Delta \mathrm{HbR}$ were smaller than those in $\Delta \mathrm{Hb}\mathrm{O}$ . On the other hand, significant simultaneous and opposing changes of $\Delta \mathrm{Hb}\mathrm{O}$ and $\Delta \mathrm{HbR}$ were observed only at the detector positions of 1-2 and 7-8 if we used the proposed method. These changes were observed at the contralateral sides in the finger opposition task. The changes in $\Delta \mathrm{Hb}\mathrm{O}$ were almost equal in magnitude to those in $\Delta \mathrm{HbR}$ .

The differences between the rest and task period averages of five trials are shown in the first (conventional method) and third (proposed method) rows of Fig. 9 . In this figure, averaged differences and their standard errors are given. The results of paired $t$ -tests are shown in the second (conventional method) and fourth (proposed method) rows of Fig. 9. Their numerical values are also given in Table 3 . The figure and table show the following. In the conventional method, the detector positions where $\Delta \mathrm{Hb}\mathrm{O}$ changed over a level of statistical significance of $p=0.05$ were not lateralized during finger opposition. $\Delta \mathrm{HbR}$ showed statistical significance for right finger opposition; however, it showed no significance for left finger opposition. These results suggest that the conventional method failed to detect a cortical functional activation. On the other hand, in the proposed method, the statistical significance of $\Delta \mathrm{Hb}\mathrm{O}$ and $\Delta \mathrm{HbR}$ were clearly lateralized. Statistical significances of the $\Delta \mathrm{Hb}\mathrm{O}$ increase and the $\Delta \mathrm{HbR}$ decrease were simultaneously observed for only two cases: 1 at the 1-2 detector position when left finger opposition was performed, and 2 at the 7-8 detector position when right finger opposition was performed. These positions correspond to locations directly above the activation area identified by fMRI. These results suggest that the proposed method works better than the conventional method in detecting cortical activation.

Fig. 9

Statistical analysis of cerebral activation during the finger opposition task. Upper and lower subfigures show the analysis results of the conventional and proposed methods, respectively. The upper row of each subfigure indicates the difference between the rest and task averages and its standard error. The lower row of each subfigure indicates $t$ -values of paired $t$ -tests. In each frame, red and blue bars indicate values of HbO and HbR, respectively. The left and right parts of each frame correspond to task condition, left- and right-finger opposition, respectively. Asterisks in some frames of $t$ -value indicate statistical significance over the level of $p=0.05$ .

Table 3

T -values of the paired t -test for the finger opposition task. Asterisks indicate statistical significance level of p=0.05 .

5.

Discussion

Artifacts caused by physical and systemic physiological activities have been a serious problem in fNIRS measurement. We confirmed the existence of artifacts in conventional fNIRS measurement using three tasks: body tilting, head nodding, and breath holding. Their magnitude was comparable to or larger than that observed as a functional signal during finger opposition.

The instability of the optical contact between probes and scalp may have caused some of these artifacts. However, these kinds of signal changes are also observed even when we use another type of probe (A9782, Hamamatsu Photonics KK), which can be attached very stably on the forehead of a subject (data not shown here). Thus, we think that physical movements and systemic physiological activity changes strongly influence conventional fNIRS measurement and may generate serious artifacts, even if we can realize stable optical contacts. Experimental results showed that such artifacts can be eliminated effectively by the proposed method.

Some researchers have reported that NIRS signals were observed at both the contralateral and ipsilateral hemisphere when a subject performs a one-sided finger opposition.^{6, 21} We confirmed the existence of such a global change in the conventional fNIRS signal. Franceschini suggested that some systemic response accompanying the finger movement causes such global changes in $\Delta \mathrm{Hb}\mathrm{O}$ and $\Delta \mathrm{HbR}$ .⁶ This component cannot be easily removed by conventional data analysis, and it certainly interferes with detection of brain function. However, if these hemodynamic changes mainly occur in the scalp or skull layers, the interference can be reduced by the proposed method, and in fact, they were reduced well in our experiment (see Fig. 8).

There are some arguments about the correlation between the fNIRS signal and brain functional activity. Some studies claim that a strong correlation exists between $\Delta \mathrm{HbR}$ decrease and the blood oxygen level dependent signal of fMRI.²² In the proposed method, the localization of $\Delta \mathrm{HbR}$ decrease coincides well with the localization of fMRI signals. Therefore, this observation supports a consistent relationship between fNIRS and fMRI measurements.

Many researchers consider that an increase in $\Delta \mathrm{Hb}\mathrm{O}$ corresponds to brain activity. In our result using the proposed method, the increase in $\Delta \mathrm{Hb}\mathrm{O}$ was observed simultaneously with the decrease in $\Delta \mathrm{HbR}$ ; consequently, it is localized in the area where the fMRI signal was observed. On the other hand, in the conventional method, a certain global component of the $\Delta \mathrm{Hb}\mathrm{O}$ increase at a superficial layer seems to be superposed on the $\Delta \mathrm{Hb}\mathrm{O}$ increase in the gray matter. Apparently this component is not directly related to brain activity, because it is well known that the cerebral region associated with motor function is highly localized.²³ If we consider such a contaminated signal of $\Delta \mathrm{Hb}\mathrm{O}$ increase as a “functional” signal, an artificial correlation of $\Delta \mathrm{Hb}\mathrm{O}$ to the task sequence may lead fNIRS into false-positive detection of brain activity.

In our multidistance probe arrangement, the distances between the source and each detector were set to 20 and $30\mathrm{mm}$ . However, some studies use another arrangement in which one detector was positioned close to the source to detect the signal from the scalp exclusively, while the other detector is positioned further away to detect the signal from gray matter. Leung, Elwell, and Delpy calculated the measurement accuracy and optimized distance conditions using various numbers of detectors.²⁴ In the case using two detectors, they showed that a distance combination of 5 and $50\mathrm{mm}$ is optimal to detect attenuation change in brain. Saager and Berger also used this type of arrangement and showed by Monte Carlo simulations that it is possible to remove artifacts uncorrelated with a functional signal.³ Our method is different from their studies in the following three points. First, we intend to remove task-correlated or function-correlated artifacts. Thus, we cannot use the low correlation property between artifacts and functional signals. Second, we intend to remove artifacts caused by both the scattering and absorption changes in superficial layers. If a source-detector distance less than $10\mathrm{mm}$ is selected, $k_a=k_s$ does not hold sufficiently (see Fig. 3), and it can counteract this type of artifact cancellation. Third, if we use their arrangements, the distance between two detectors becomes greater than $20\mathrm{mm}$ . This can violate the assumption of uniform slab in a two-layer model, because vasculature and hemodynamics at the actual superficial layer may not be sufficiently homogeneous.

We adopted the simplification that only functional hemodynamics occur in the gray matter layer. Under this simplification, detected signals from gray matter directly indicate a functional signal. However, systemic physiological activities can affect hemodynamics in the gray matter layer as well as in the superficial layer. Therefore, we discuss in Appendix { label needed for app[@id='xB'] } in Sec. 7 how systemic hemodynamics in gray matter, if it exists, affects the detection of brain function. In brief, the influence of systemic hemodynamics in gray matter can also be removed by the proposed method; see Appendix { label needed for app[@id='xB'] } in Sec. 7 for details.

In this work, data obtained from one single subject were shown. We implemented similar experiments on several subjects. Also in these cases, we confirmed effective artifact reduction and functional signal localization; however, it may require further validation with more subjects. Through further trials of the proposed method, $\Delta \mathrm{Hb}\mathrm{O}$ and $\Delta \mathrm{HbR}$ changes by cerebral activation of human adults, infants, and pathological cases may be revealed in the future. We expect that these can be adopted as cerebral hemodynamic response functions (HRF) in the statistical approaches based on the general linear model.^{25, 26} Such a statistical parameter mapping (SPM)-like analysis based on a precise HRF will increase the reliability of fNIRS measurement.