## 1.

## Introduction

Near-infrared spectroscopy (NIRS) has several advantages compared to other functional measurement methods, such as functional magnetic resonance imaging (fMRI), positron emission tomography, and electroencephalography (EEG). The advantages include good temporal resolution and measurement of both oxygenated $\left(\mathrm{Hb}{\mathrm{O}}_{2}\right)$ and deoxygenated (HbR) hemoglobin, portability, low restraints of subjects, and low cost of the measurement equipment. The disadvantages are low spatial resolution and potentially unstable measurement because of the instability of an optical probe contact that might be induced by body motion and hair absorption. The systemic physiological interference is an another difficult problem in NIRS measurement.

The well-known systemic physiological interference is caused by heart beat and respiration. The suppression of this kind of systemic interference is relatively easy because the correlation between the interference and the functional response of brain functional activity is usually low. For example, systemic interference caused by heart beat can be effectively removed using a low pass filtering technique. However, low pass filtering may not be appropriate for respiration because the frequency spectra of respiration and functional activity may overlap each other. Even in this case, it can be suppressed by using the low correlation property between the interference (respiration) and brain functional activity. Zhang solved this problem by using an adaptive filtering technique.^{1, 2} Saager and Berger applied the least-squares fit technique to this problem, and gave a detailed analysis of its performance.^{3} Both methods are based on the multidistance measurement and have an advantage that they can remove not only the interference in the superficial layer (scalp and skull layers), but also the interference included in the cerebral layer, as far as the interference is low correlated with the brain functional signal.

There is an another form of interference, which is observed synchronously with the heart rate increase during execution of a task (for example, finger opposition task).^{4, 5} This kind of interference is highly correlated with the functional response itself. Thus, it is difficult to remove it by the methods using the low correlation property between them. We consider the removal of this kind of interference in this paper.

Two approaches have been proposed to isolate and eliminate this kind of interference. The first approach is to use the spatial locality difference between the brain functional activity and the systemic interference. It is well known through fMRI studies that the brain functional activity is highly localized. Conversely, systemic interference is generally global. Thus, this difference can be used to differentiate them. Pioneering work of this approach was reported by Franceschini
^{4} They assumed that the NIRS signal measured at an area where no functional activity is expected should be similar to the systemic interference, and isolated the functional response of the activated region by subtracting this signal from the measured signal. This approach was refined as a method that uses principal component analysis.^{5, 6} A similar approach based on independent component analysis instead of PCA has also been proposed.^{7} This approach has an advantage in that a conventional continuous-wave NIRS system is available, and the method is easy to implement. However, it has a tendency to decrease the amplitude of the functional response in the activated regions and to propagate noise from noisy channels to all other channels.^{5} Another disadvantage is that it requires a large-scale NIRS system equipping many measurement channels because it is based on the spatial uniformity of the systemic interference.

The second approach relies on the multidistance measurement. If we assume a multi layered slab model as a human adult head, the temporal absorbance change can be approximated as a linear sum of the products of a partial pathlength and an absorption coefficient change of each layer.
^{8, 9} Because the partial pathlength depends on the distance between a light source and a detector, absorption coefficient change at each layer can be computed easily if several detectors with different source-detector separations are used and partial pathlength parameters of each source-detector pair are known.
^{8, 10} Thus, if such an interference occurs mainly in the superficial layer, multidistance measurement method can separate it from the functional response.

The problem with the second approach is that it is very difficult to determine, experimentally, the partial pathlength of each layer. Fabbri
^{10} pointed out that the ratio of the partial pathlength of the superficial layer at different detectors can be calculated as the ratio of the absorbance changes, if the hemoglobin changes in the cerebral layer are negligible. This approach is effective. However, requirement of static hemodynamics in the cerebral layer may be a strong constraint in some situations.

A partial pathlength of each layer can be predicted computationally by Monte Carlo simulation. Monte Carlo simulation requires a lot of physiological parameters of a head model, such as the thickness and transport scattering coefficient of each tissue. Tissue thickness can be measured by MRI, however, doing so may spoil the advantages of low restraint of subjects and low cost provided by NIRS. For example, although NIRS offers advantages for brain functional measurement of infants and children because of its low restraint property, its applicability would be limited if MRI measurement is required. The transport scattering coefficient of each layer is much more difficult to measure.

Leung
^{11} have provided a theoretical analysis of the second approach by the finite element method of the diffusion equation based on assumed standard physiological parameters. They showed how exactly the hemoglobin change at gray matter can be estimated using this approach and gave the optimum detector arrangement. In this paper, deviations of blood volume, oxygen saturation, and thickness of a superficial layer are considered, and a range of data were simulated corresponding to different thicknesses and physiological conditions of each layer. The optimum regression parameters for brain hemodynamic change estimation were calculated directly from the simulated data without using the modified Beer-Lambert law. They used the partial least-squares method to obtain these parameters. In their paper, however, only the thickness parameter of a superficial layer was considered as a deviation, and the transport scattering parameter of all layers and thickness parameter of the other layers were not considered.

Here, we analyzed the multidistance measurement method to separate the brain functional response from the systemic interference in a superficial layer caused by task execution. Other interferences, evoked by heart beat, respiration, etc., were assumed to be removed beforehand by preprocessing. We used the modified Beer-Lambert law, and partial optical path lengths were estimated by Monte Carlo simulation. We supposed a population of subjects and assumed that physiological parameters, the transport scattering parameters, and thicknesses of layers of each subject were distributed around the assumed standard values. An estimation method giving a good average separation performance for the population was derived by multilinear regression.

## 2.

## Theory

Our study used two light detectors ( ${d}_{1}$ and ${d}_{2}$ ) and two wavelengths ( ${\lambda}_{1}$ and ${\lambda}_{2}$ ) of light. A five-layer slab model consisting of scalp, skull, cerebrospinal fluid (CSF), gray matter, and white matter was used as a human adult head model.

## 2.1.

### NIRS Observation Model Based on Partial Pathlength

The temporal absorbance change at the detector
${d}_{i}$
$(i=1,2)$
of wavelength
${\lambda}_{j}$
$(j=1,2)$
is represented as follows:^{8, 9}

## Eq. 1

$\Delta {A}_{{d}_{i},{\lambda}_{j}}={l}_{{d}_{i},{\lambda}_{j}}^{\mathrm{sc}}\Delta {\mu}_{\mathrm{a},{\lambda}_{j}}^{\mathrm{sc}}+{l}_{{d}_{i},{\lambda}_{j}}^{\mathrm{sk}}\Delta {\mu}_{\mathrm{a},{\lambda}_{j}}^{\mathrm{sk}}+{l}_{{d}_{i},{\lambda}_{j}}^{\mathrm{csf}}\Delta {\mu}_{\mathrm{a},{\lambda}_{j}}^{\mathrm{csf}}+{l}_{{d}_{i},{\lambda}_{j}}^{\mathrm{gm}}\Delta {\mu}_{\mathrm{a},{\lambda}_{j}}^{\mathrm{gm}}+{l}_{{d}_{i},{\lambda}_{j}}^{\mathrm{wm}}\Delta {\mu}_{\mathrm{a},{\lambda}_{j}}^{\mathrm{wm}},$## Eq. 2

$\Delta {A}_{{d}_{i},{\lambda}_{j}}={l}_{{d}_{i},{\lambda}_{j}}^{\mathrm{sp}}\Delta {\mu}_{\mathrm{a},{\lambda}_{j}}^{\mathrm{sp}}+{l}_{{d}_{i},{\lambda}_{j}}^{\mathrm{gm}}\Delta {\mu}_{\mathrm{a},{\lambda}_{j}}^{\mathrm{gm}},$## Eq. 3

${l}_{{d}_{i},{\lambda}_{j}}^{\mathrm{sp}}={l}_{{d}_{i},{\lambda}_{j}}^{\mathrm{sc}}+{l}_{{d}_{i},{\lambda}_{j}}^{\mathrm{sk}},$## Eq. 4

$\Delta {\mu}_{\mathrm{a},{\lambda}_{j}}^{\mathrm{sp}}=\Delta {\mu}_{\mathrm{a},{\lambda}_{j}}^{\mathrm{sc}}=\Delta {\mu}_{\mathrm{a},{\lambda}_{j}}^{\mathrm{sk}}.$## Eq. 5

$\mathit{a}={(\Delta {A}_{{d}_{1},{\lambda}_{1}},\Delta {A}_{{d}_{2},{\lambda}_{1}},\Delta {A}_{{d}_{1},{\lambda}_{2}},\Delta {A}_{{d}_{2},{\lambda}_{2}})}^{T},$## Eq. 6

$\mathbf{\mu}={(\Delta {\mu}_{\mathrm{a},{\lambda}_{1}}^{\mathrm{sp}},\Delta {\mu}_{\mathrm{a},{\lambda}_{1}}^{\mathrm{gm}},\Delta {\mu}_{\mathrm{a},{\lambda}_{2}}^{\mathrm{sp}},\Delta {\mu}_{\mathrm{a},{\lambda}_{2}}^{\mathrm{gm}})}^{T}.$## Eq. 8

$L=\left(\begin{array}{cccc}{l}_{{d}_{1},{\lambda}_{1}}^{\mathrm{sp}}& {l}_{{d}_{1},{\lambda}_{1}}^{\mathrm{gm}}& 0& 0\\ {l}_{{d}_{2},{\lambda}_{1}}^{\mathrm{sp}}& {l}_{{d}_{2},{\lambda}_{1}}^{\mathrm{gm}}& 0& 0\\ 0& 0& {l}_{{d}_{1},{\lambda}_{2}}^{\mathrm{sp}}& {l}_{{d}_{1},{\lambda}_{2}}^{\mathrm{gm}}\\ 0& 0& {l}_{{d}_{2},{\lambda}_{2}}^{\mathrm{sp}}& {l}_{{d}_{2},{\lambda}_{2}}^{\mathrm{gm}}\end{array}\right).$Next we give a vector representation of $\mathrm{Hb}{\mathrm{O}}_{2}$ and HbR concentration change in superficial and gray matter layers

## Eq. 9

$\mathit{x}={(\Delta \mathrm{Hb}{\mathrm{O}}_{2}^{\mathrm{sp}},\Delta \mathrm{Hb}{\mathrm{O}}_{2}^{\mathrm{gm}},\Delta \mathrm{Hb}{\mathrm{R}}^{\mathrm{sp}},\Delta \mathrm{Hb}{\mathrm{R}}^{\mathrm{gm}})}^{T}.$## Eq. 11

$E=\left(\begin{array}{cccc}{\epsilon}_{\mathrm{Hb}{\mathrm{O}}_{2},{\lambda}_{1}}& 0& {\epsilon}_{\mathrm{Hb}\mathrm{R},{\lambda}_{1}}& 0\\ 0& {\epsilon}_{\mathrm{Hb}{\mathrm{O}}_{2},{\lambda}_{1}}& 0& {\epsilon}_{\mathrm{Hb}\mathrm{R},{\lambda}_{1}}\\ {\epsilon}_{\mathrm{Hb}{\mathrm{O}}_{2},{\lambda}_{2}}& 0& {\epsilon}_{\mathrm{Hb}\mathrm{R},{\lambda}_{2}}& 0\\ 0& {\epsilon}_{\mathrm{Hb}{\mathrm{O}}_{2},{\lambda}_{2}}& 0& {\epsilon}_{\mathrm{Hb}\mathrm{R},{\lambda}_{2}}\end{array}\right).$## 2.2.

### Optimizing the Estimation of $\Delta \mathrm{Hb}{\mathrm{O}}_{2}$ and $\Delta \mathrm{Hb}\mathrm{R}$

Hemoglobin concentration change can be estimated from the measured temporal absorbance change
$\mathit{a}$
and Eq. 12 if we know
$E$
and
$L$
.
$E$
can be obtained from prior studies.^{12} However, partial pathlength matrix
$L$
is determined by physiological parameter values, which depend on a subject and a measurement position, and it is difficult to obtain experimentally. We describe in this paper a separation method that provides a good average separation performance for a variety of subjects and measurement positions.

The various partial pathlength matrices are written as $\mathcal{L}=\{{\mathrm{L}}_{1},{\mathrm{L}}_{2},\dots ,{\mathrm{L}}_{K}\}$ corresponding to a variety of subjects and measurement positions. $\mathcal{L}$ indicates that various values of the temporal absorbance change ${\mathit{a}}_{k}={L}_{k}E\mathit{x}$ are observed even if the hemoglobin concentration change $\mathit{x}$ is fixed. We introduce an inverse transformation matrix ${E}^{-1}W$ , which estimates hemoglobin concentration change from temporal absorbance change values. $W$ estimates absorption coefficient change from the temporal absorbance change, and ${E}^{-1}$ transforms absorption coefficient change to hemoglobin concentration change. Here we assume that the absorption coefficient changes of two wavelengths are estimated independently. Thus, $W=\mathrm{diag}({W}_{1},{W}_{2})$ , where ${W}_{1}$ and ${W}_{2}$ are estimation matrices for ${\lambda}_{1}$ and ${\lambda}_{2}$ , respectively. Because $W$ includes coefficients to estimate absorption coefficient changes, we call $W$ a measurement weighting.

When we use a measurement weighting $W$ , the hemoglobin concentration change when $L={L}_{k}$ is estimated as follows:

The estimation error is## Eq. 15

${\widehat{\mathit{x}}}_{k}\left(t\right)-\mathit{x}\left(t\right)=({E}^{-1}W{L}_{k}E-I)\mathit{x}\left(t\right).$## Eq. 16

${J}_{k}=\frac{1}{T}{\int}_{0}^{T}{[{\widehat{\mathit{x}}}_{k}\left(t\right)-\mathit{x}\left(t\right)]}^{T}[{\widehat{\mathit{x}}}_{k}\left(t\right)-\mathit{x}\left(t\right)]\mathrm{d}t,$## Eq. 17

$=\frac{1}{T}{\int}_{0}^{T}\mathit{x}{\left(t\right)}^{T}{({E}^{-1}W{L}_{k}E-I)}^{T}({E}^{-1}W{L}_{k}E-I)\mathit{x}\left(t\right)\mathrm{d}t,$## Eq. 19

${\Sigma}_{x}=\frac{1}{T}{\int}_{0}^{T}\mathit{x}\left(t\right){\mathit{x}}^{T}\left(t\right)\mathrm{d}t.$The optimum weighting $W=\mathrm{diag}({W}_{1},{W}_{2})$ to minimize $J$ and its derivation detail is given in the Appendix. The optimum weighting [Eq. 43] gives a good average separation performance for a partial path-length matrix distribution $\mathcal{L}$ on average. Thus, if we apply this weighting to the multidistance measurement data, it is expected to give a better separation performance than the weighting determined without considering a parameter distribution. This will be illustrated later in Fig. 5.

## 2.3.

### Random Generation of Partial Path-length Matrices

Because a partial path-length matrix $L$ is a block diagonal matrix, we write $L=\mathrm{diag}({\Lambda}^{1},{\Lambda}^{2})$ . ${\Lambda}^{j}$ corresponds to a measurement of ${\lambda}_{j}$ light. The value of a partial path-length depends on the thickness and transport scattering coefficient of each layer of a head. Thus, we use ${\Lambda}^{j}\left({\mathbf{u}}^{j}\right)$ to emphasize this. The set of parameters ${\mathbf{u}}^{j}$ are as follows:

## Eq. 21

${\mathbf{u}}^{j}={({u}_{1}^{j},{u}_{2}^{j},{u}_{3}^{j},{u}_{4}^{j},{u}_{5}^{j},{u}_{6}^{j},{u}_{7}^{j},{u}_{8}^{j},{u}_{9}^{j})}^{T}$## Eq. 22

$={({\mu}_{\mathrm{s},{\lambda}_{j}}^{\prime \mathrm{sc}},{\mu}_{\mathrm{s},{\lambda}_{j}}^{\prime \mathrm{sk}},{\mu}_{\mathrm{s},{\lambda}_{j}}^{\prime \mathrm{csf}},{\mu}_{\mathrm{s},{\lambda}_{j}}^{\prime \mathrm{gm}},{\mu}_{\mathrm{s},{\lambda}_{j}}^{\prime \mathrm{wm}},{t}^{\mathrm{sc}},{t}^{\mathrm{sk}},{t}^{\mathrm{csf}},{t}^{\mathrm{gm}})}^{T}.$A set of partial path-length matrices can be obtained if we execute a Monte Carlo simulation based on various multilayered model parameters. This simulation, however, requires significant computation time; thus, execution of many sets of parameters is difficult. Instead, we assumed that ${\Lambda}^{j}({\mathbf{u}}_{0}^{j}+\Delta {\mathbf{u}}^{j})$ can be approximated in a linear form if $\Delta {\mathbf{u}}^{j}$ is small

## Eq. 23

${\Lambda}^{j}({\mathbf{u}}_{0}^{j}+\Delta {\mathbf{u}}^{j})={\Lambda}^{j}\left({\mathbf{u}}_{0}^{j}\right)+\sum _{m=1}^{9}{\phantom{|}\frac{\partial {\Lambda}^{j}}{\partial {u}_{m}}|}_{{\mathbf{u}}_{0}^{j}}\Delta {u}_{m}^{j}.$Many partial path-length matrices can be generated based on Eq. 23. The first-order coefficients of Eq. 23 indicate how much the partial path length changes when thicknesses and transport scattering coefficients of layers change; we term this the sensitivity coefficient of the partial pathlength. The approximate value of the sensitivity coefficient can be obtained if we investigate the partial pathlength change produced by a small increase in the parameter value ${u}_{m}$ through simulations.

We computed the partial-pathlength and sensitivity coefficient of standard parameters ${u}_{m,0}^{j}$ by Monte Carlo simulation. We assumed that a parameter ${u}_{m}^{j}$ is a random variable and obeys a normal distribution $N({u}_{m,0}^{j},{\sigma}_{m}^{j})$ . If we determine a parameter value ${u}_{m}^{j}$ by random sampling of $N({u}_{m,0}^{j},{\sigma}_{m}^{j})$ , then many partial pathlength matrices can be generated by Eq. 23.

## 3.

## Method

## 3.1.

### Monte Carlo Simulation

The partial path-length and sensitivity coefficient based on standard parameter values were computed by a Monte Carlo simulation. The simulation program, tMCimg, provided in public by the Photon Migration Imaging Laboratory at Massachusetts General Hospital was used. A five-layer model of an adult head consisting of scalp, skull, CSF, gray matter, and white matter was used (Fig. 1 ). The size of simulated tissues was $100\times 100\times 50\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ . The wavelengths of the light source were 800 and $840\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ . The source-detector separation was assumed from $10\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}40\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ every $5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ .

The standard parameters used in the simulation are given in Table 1
. Thickness and optical properties (
${\mu}_{\mathrm{a}}$
and
${\mu}_{\mathrm{s}}^{\prime}$
) for near-infrared light of
$800\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
for each layer were chosen from reported data.^{13}
$\mathrm{Hb}{\mathrm{O}}_{2}$
and HbR concentrations were calculated from the absorption coefficients at
$800\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
with the assumption that oxygen saturation is 70%. The absorption coefficient calculated from water absorption coefficients given in the paper^{14} was also used in this case.^{15} These values were used to calculate the absorption coefficients at
$840\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
. Because it is known that the transport scattering coefficient decreases if the wavelength of infrared light increases,^{11, 16} transport scattering coefficients at
$840\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
were set to be 95% of those at
$800\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
.

## Table 1

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

Tissue type | Baselineconcentration | Thickness(mm) | Absorptioncoefficient μa 800nm∕840nm (mm−1) | Transport scatteringcoefficient μs′ 800nm∕840nm (mm−1) | |
---|---|---|---|---|---|

HbO2 (mM) | HbR(mM) | ||||

Scalp | 0.064 | 0.027 | 3 | $0.018\u22150.021$ | $1.90\u22151.81$ |

Skull | 0.057 | 0.024 | 7 | $0.016\u22150.019$ | $1.60\u22151.52$ |

CSF | 0.014 | 0.006 | 2 | $0.004\u22150.005$ | $0.24\u22150.23$ |

Gray matter | 0.128 | 0.055 | 4 | $0.036\u22150.042$ | $2.20\u22152.09$ |

White matter | 0.050 | 0.021 | 34 | $0.014\u22150.017$ | $9.10\u22158.65$ |

## 3.2.

### Model Time Course of Hemodynamic Change

It is very difficult to know real values $\mathit{x}={(\Delta \mathrm{Hb}{\mathrm{O}}^{\mathrm{sp}},\Delta \mathrm{Hb}{\mathrm{O}}_{2}^{\mathrm{gm}},\Delta \mathrm{Hb}{\mathrm{R}}^{\mathrm{sp}},\Delta \mathrm{Hb}{\mathrm{R}}^{\mathrm{gm}})}^{T}$ in superficial and gray-matter layers. Thus, we use a simple model of hemodynamic change time course. In this model, $\Delta \mathrm{Hb}{\mathrm{O}}_{2}^{\mathrm{gm}}$ and $\Delta \mathrm{Hb}{\mathrm{R}}^{\mathrm{gm}}$ are functional responses of the executed task, and $\Delta \mathrm{Hb}{\mathrm{O}}^{\mathrm{sp}}$ and $\Delta \mathrm{Hb}{\mathrm{R}}^{\mathrm{sp}}$ are systemic interference correlated with the task execution.

The hemodynamic responses in both layers were defined as the convolution of the stimulation
$s\left(t\right)$
, [
$s\left(t\right)=0$
for rest period, and 1 for stimulation] and a prototypical hemodynamic impulse response
$h\left(t\right)$
,^{1, 17} in our model

## Eq. 24

$s\left(t\right)=\{\begin{array}{ll}1& \text{if}\phantom{\rule{0.3em}{0ex}}t\u220a\text{rest}\\ 0& \text{if}\phantom{\rule{0.3em}{0ex}}t\u220a\text{stimulation}\end{array}\phantom{\}},$^{18}show that $\Delta \mathrm{Hb}{\mathrm{O}}^{\mathrm{gm}}$ and $\Delta \mathrm{Hb}{\mathrm{R}}^{\mathrm{gm}}$ have a similar time course with a very little delay and that the functional response is slower than the systemic response. The parameter concerning the impulse response width was set at $b=8$ for all responses. Figure 2a shows the stimulation $s\left(t\right)$ $(0\u2a7dt\u2a7d40)$ , and Fig. 2b shows the hemodynamic impulse response of each layer. The peak values of all responses are normalized to be unity.

The evoked hemodynamic response $x\left(t\right)$ was the convolution of $s\left(t\right)$ and $h\left(t\right)$

where $C$ is a scaling parameter. Figure 2c shows the simulated evoked hemodynamic response of each layer. Peak values of $\Delta \mathrm{Hb}{\mathrm{O}}^{\mathrm{gm}}$ and $\Delta \mathrm{Hb}{\mathrm{R}}^{\mathrm{gm}}$ were set to be 0.017 and $-0.005\phantom{\rule{0.3em}{0ex}}\mathrm{mM}$ , respectively. These values were about 13% and 9% of the assumed baseline conditions.^{19}On the other hand, peak values of the superficial layer were set to be $1\u221510$ of the gray-matter layer based on the following findings: The temporal absorbance changes of the systemic interference of a finger-tapping task and functional response of brain functional activity are considered to be comparable

^{4}and Monte Carlo simulation shows that the partial path length of a superficial layer is more than 10 times longer than that of a gray-matter layer. Thus, we assumed that hemoglobin concentration changes of a superficial layer are $1\u221510$ of those of gray-matter layer.

The covariance matrix ${\Sigma}_{x}$ [Eq. 19] was calculated based on these simulated hemodynamic responses and used for the optimization process.

## 4.

## Results and Discussion

## 4.1.

### Partial Path Length and Its Sensitivity Coefficient

We computed partial path lengths of superficial and gray-matter layers at various source-detector separations by Monte Carlo simulation (Fig. 3 ). Solid lines and dashed lines in Fig. 3 correspond to the light sources of wavelengths ${\lambda}_{1}$ and ${\lambda}_{2}$ , respectively. Figure 3 shows that partial path lengths of different wavelength have similar values, especially for gray matter. It also shows that the partial path length of the superficial layer is more than 10 times greater than that of gray matter. This difference is pronounced when the source-detector separation is small, which means that blood flow change in a superficial layer may have a strong influence on the measurement.

Figure 4 shows relative (percent) change of partial pathlength of gray-matter and superficial layers, that is, the ratio of partial path-length change to its baseline value, where model parameters were increased by 10% from the standard values. We showed these values instead of sensitivity coefficients because they may be more comprehensible. Figure 4 indicates that the thickness parameter of scalp and skull layers has a significant influence on the partial path length of the gray-matter and superficial layers. The transport scattering coefficient of the skull layer has also a comparable influence on the partial pathlength of the gray-matter layer. On the other hand, partial path lengths of both layers are insensitive to all other parameters.

The validity of the first-order approximation of the partial path-length matrix [Eq. 23] was checked by another Monte Carlo simulation. The simulation result (data not shown) shows Eq. 23 holds approximately in the range of $\pm 30\%$ of the transport scattering coefficient change and $\pm 1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ of the layer thickness change.

## 4.2.

### Optimum Weighting and Its Performance

To obtain a measurement weighting that provides a good average performance for the subject population, we generated a set $\mathcal{L}$ of partial path-length matrices based on its linear approximation model [Eq. 23] and the assumption that model parameters are distributed normally around the standard parameter values. The source-detector separations of two detectors were set to be 20 and $30\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ . The optimum weighting ${W}_{\mathrm{opt}}$ for $\mathcal{L}$ was derived by Eq. 43. Standard deviations of normal distributions of model parameters were set to be 10, 20, and 30% of standard parameter values shown in Table 1. Two thousand partial path-length matrices were randomly generated for each case, and half of them were used as training data to compute ${W}_{\mathrm{opt}}$ and the other half as test data to evaluate the performance of the obtained ${W}_{\mathrm{opt}}$ .

When we use a simulated evoked hemodynamic response of Fig. 2, the mean-squared error (MSE) of $\Delta \mathrm{Hb}{\mathrm{O}}^{\mathrm{gm}}$ and $\Delta \mathrm{Hb}{\mathrm{R}}^{\mathrm{gm}}$ estimated by two measurement weightings, ${W}_{\mathrm{opt}}$ and ${L}_{\mathrm{std}}^{-1}$ are shown in Figs. 5a and 5b. ${W}_{\mathrm{opt}}$ is the optimum weighting obtained by the training data, and its MSEs are shown by filled bars. ${L}_{\mathrm{std}}^{-1}$ is the inverse of the standard partial path-length matrix ${L}_{\mathrm{std}}$ , where ${L}_{\mathrm{std}}$ is a partial path-length matrix determined by the standard model parameters. This means that ${L}_{\mathrm{std}}^{-1}$ is a weighting not to consider the model parameter distribution. Its MSEs are shown by unfilled bars. Figures 5a and 5b show that MSE of both $\Delta \mathrm{Hb}{\mathrm{O}}^{\mathrm{gm}}$ and $\Delta \mathrm{Hb}{\mathrm{R}}^{\mathrm{gm}}$ increases almost linearly when deviations of model parameters increases. It also shows that MSEs of the optimum weighting are much less than those obtained by ${L}_{\mathrm{std}}^{-1}$ , especially when the deviation of model parameter distribution is large. Thus, the proposed optimization method is effective to minimize the estimation errors on average.

Figures 5c and 5d show MSEs of the case where the hemoglobin changes in the superficial layer is larger than the case of Fig. 2. In this case, the peak values of $\Delta \mathrm{Hb}{\mathrm{O}}^{\mathrm{sp}}$ and $\Delta \mathrm{Hb}{\mathrm{R}}^{\mathrm{sp}}$ were set $1\u22155\text{th}$ of the gray matter layer values. Figures 5c and 5d show that MSEs of the optimum weighting slightly increase in this case; however, MSEs of ${L}_{\mathrm{std}}^{-1}$ are significantly larger than the small artifact case (Fig. 2 case). This result shows that the optimum weighting can give a better and robust estimation of brain hemodynamic response on average.

Figure 6 shows how much the estimated hemoglobin changes deviate from the assumed original changes (Fig. 2) when model parameter changes. The standard deviation of model parameter distribution was assumed to be 20% of the standard values. Ten partial path-length matrices were randomly selected from the generated test data, and two measurement weightings ( ${W}_{\mathrm{opt}}$ and ${L}_{\mathrm{std}}^{-1}$ ) were used to estimate the hemoglobin changes based on Eq. 14. The estimated hemoglobin changes are light red (HbO) and blue (HbR) solid lines. The assumed original hemoglobin changes are also shown in the figures by thick solid lines. The estimation result by ${W}_{\mathrm{opt}}$ and ${L}_{\mathrm{std}}^{-1}$ are shown in Figs. 6a and 6b, respectively, which show that the estimated changes by ${L}_{\mathrm{std}}^{-1}$ vary widely compared to those by ${W}_{\mathrm{opt}}$ .

Figure 6c is the estimation result where we applied an ordinary single-detector NIRS algorithm to the test data under the same conditions. The source-detector separation was set to be
$30\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
in this case. The sum of the partial path lengths of all layers was considered as the total path length in this case. The obtained hemoglobin changes were significantly smaller than the original changes because the total path length used in the algorithm is longer than the partial path length of the gray-matter layer. This is called partial volume effect.^{20}

Next, we normalized the scale of each response in Figs. 6a, 6b, 6c by making its value at
$t=15\phantom{\rule{0.3em}{0ex}}\mathrm{s}$
, 1 for HbO and
$-1$
for HbR, to analyze the response shapes without being affected by their scale variation. The original changes were also normalized in this way. The results are shown in Figs. 6d, 6e, 6f. Figures 6d and 6e show that the shapes of the responses estimated by multidistance measurement are close to the original hemoglobin changes in the gray-matter layer, although the rising and decaying phases of the response are distorted. These distortions are considered to be caused by cross talk.^{15, 20} On the other hand, the shapes estimated by an ordinary single-detector NIRS algorithm [Fig. 6f] are different from the original changes in gray-matter layer but close to those in the superficial layer. These results show that the multidistance measurement is effective to extract the correct hemoglobin changes in the gray-matter layer.

## 5.

## Conclusion

The multidistance measurement method can separate the functional response of brain activity from the systemic physiological interference caused by task execution. However, because measurement weighting is difficult to determine experimentally, we estimated it by a Monte Carlo simulation. A five-layered slab model was used to represent a human adult head. Because model parameters such as the thickness and transport scattering coefficient depend on a subject and a measurement position and are difficult to measure, we assumed that model parameters obey normal distributions around standard parameter values. The measurement weighting that produces a good average separation performance for the model parameter distribution was determined.

Partial path lengths of superficial and gray-matter layers were computed based on standard parameters by a Monte Carlo simulation, and their sensitivity coefficients, which indicate how much the partial path length changes when model parameter changes, were also given. The results indicate that the partial path length of the superficial layer is much longer than that of the gray-matter layer. The thickness parameter of scalp and skull layers has a significant influence on the partial path length of the gray-matter and superficial layers. The transport scattering coefficient of the skull layer has also a comparable influence on the partial path length of the gray-matter layer. Partial path lengths are insensitive to all other parameters.

Using the obtained partial path lengths and their sensitivity coefficients, a measurement weighting designed to yield the average optimum separation performance for model parameter distribution was derived. The obtained optimum weighting is robust to model parameter deviation and realizes smaller MSEs compared to a weighting that is determined without considering model parameter distribution.

When optimizing a measurement weighting, the covariance matrix
${\Sigma}_{x}$
of hemoglobin changes is required. In this paper, we assumed a simple model of a hemodynamic change time course according to the recent study.^{18}
${\Sigma}_{x}$
was calculated by this model. Although we considered a statistical distribution of layer parameters, the variations of hemodynamic change was not considered in this paper to simplify the analysis.

A significant systemic physiological interference, which is highly correlated with the functional response, is observed for a motor task (finger opposition), and only a small interference is observed for a cognitive and a visual tasks.^{5} Thus, the proposed method is effective mainly for a motor task experiment. However, this method can also suppress the systemic interference in the extracerebral layer caused by heart beat and respiration, although the theoretical analysis is given without considering these interferences in this paper to make the analysis simple. Thus, the proposed method is expected to be effective for many cases.

In our method, the hemoglobin concentration change is estimated by applying an inverse of molar absorption coefficient matrix and an optimized measurement weighting to the absorbance change obtained by the multidistance measurement [Eq. 13]. The applicability of the method strongly depends on the accuracy of the assumed model and its parameters. Thus, its naive application to a real measurement may be inappropriate. As shown in Fig. 6, however, because an ordinary single-detector NIRS measurement may lead to large estimation errors if there exists a hemoglobin concentration change in a superficial layer, we think that the proposed method is valuable even when we have only a rough knowledge of the model and its parameter values.

## Appendices

### Appendix: Derivation of Optimum Measurement Weighting, $\mathit{W}$

The criterion $J$ of Eq. 20 is represented as follows:

## Eq. 27

$J=\mathrm{tr}\left\{F{W}^{T}GW\right\}-2\phantom{\rule{0.3em}{0ex}}\mathrm{tr}\left(H{W}^{T}\right)+\mathrm{tr}{\Sigma}_{x},$## Eq. 30

$H={\left({E}^{-1}\right)}^{T}{\Sigma}_{x}{E}^{T}\left(\frac{1}{K}\sum _{k=1}^{K}{L}_{k}^{T}\right).$## Eq. 34

$J=\mathrm{tr}\left\{{F}_{11}{W}_{1}^{T}{G}_{11}{W}_{1}\right\}+\mathrm{tr}\left\{{F}_{12}{W}_{2}^{T}{G}_{21}{W}_{1}\right\}+\mathrm{tr}\left\{{F}_{21}{W}_{1}^{T}{G}_{12}{W}_{2}\right\}+\mathrm{tr}\left\{{F}_{22}{W}_{2}^{T}{G}_{22}{W}_{2}\right\}-2\phantom{\rule{0.3em}{0ex}}\mathrm{tr}\left({H}_{11}{W}_{1}^{T}\right)-2\phantom{\rule{0.3em}{0ex}}\mathrm{tr}\left({H}_{22}{W}_{2}^{T}\right)+\mathrm{tr}{\Sigma}_{n}.$## Eq. 35

$\frac{\partial J}{\partial {W}_{1}}=2{G}_{11}{W}_{1}{F}_{11}+{G}_{21}^{T}{W}_{2}{F}_{12}^{T}+{G}_{12}{W}_{2}{F}_{21}-2{H}_{11}$## Eq. 37

$\frac{\partial J}{\partial {W}_{2}}=2{G}_{21}{W}_{1}{F}_{12}+2{G}_{22}{W}_{2}{F}_{22}-2{H}_{11}.$## Eq. 40

$\{{F}_{11}\otimes {G}_{11}\}\mathrm{cs}{W}_{1}+\{{F}_{21}^{T}\otimes {G}_{12}\}\mathrm{cs}{W}_{2}=\mathrm{cs}{H}_{11},$## Eq. 41

$\{{F}_{12}^{T}\otimes {G}_{21}\}\mathrm{cs}{W}_{1}+\{{F}_{22}\otimes {G}_{22}\}\mathrm{cs}{W}_{2}=\mathrm{cs}{H}_{11}.$## Eq. 42

$\left(\begin{array}{cc}\{{F}_{11}\otimes {G}_{11}\}& \{{F}_{21}^{T}\otimes {G}_{12}\}\\ \{{F}_{12}^{T}\otimes {G}_{21}\}& \{{F}_{22}\otimes {G}_{22}\}\end{array}\right)\left(\begin{array}{c}\mathrm{cs}{W}_{1}\\ \mathrm{cs}{W}_{2}\end{array}\right)=\left(\begin{array}{c}\mathrm{cs}{H}_{11}\\ \mathrm{cs}{H}_{22}\end{array}\right).$## Eq. 43

$\left(\begin{array}{c}\mathrm{cs}{W}_{1}\\ \mathrm{cs}{W}_{2}\end{array}\right)={\left(\begin{array}{cc}\{{F}_{11}\otimes {G}_{11}\}& \{{F}_{21}^{T}\otimes {G}_{12}\}\\ \{{F}_{12}^{T}\otimes {G}_{21}\}& \{{F}_{22}\otimes {G}_{22}\}\end{array}\right)}^{-1}\left(\begin{array}{c}\mathrm{cs}{H}_{11}\\ \mathrm{cs}{H}_{22}\end{array}\right).$## Acknowledgment

The authors thank the Photon Migration Imaging Laboratory at Massachusetts General Hospital for the Monte Carlo simulation code (available at http://www.nmr.mgh.harvard.edu/PMI/index.htm).

## References

*In vivo*measurements of the wavelength dependence of tissue-scattering coefficients between 760 and $900\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ measured with time-resolved spectroscopy,” Appl. Opt., 36 386 –396 (1997). https://doi.org/10.1364/AO.36.000386 0003-6935 Google Scholar