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}

Next we give a vector representation of $\mathrm{Hb}{\mathrm{O}}_2$ and HbR concentration change in superficial and gray matter layers

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.¹² 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 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.

Fig. 5

Mean-squared errors of $\Delta \mathrm{Hb}{\mathrm{O}}^{\mathrm{gm}}$ and $\Delta \mathrm{Hb}{\mathrm{R}}^{\mathrm{gm}}$ estimated by $W_{\mathrm{opt}}$ (filled bars) and $L_{\mathrm{std}}^{-1}$ (unfilled bars). Standard deviations of normal distributions of model parameters were set to be 10, 20, and 30% of standard parameter values shown in Table 1. (a) and (c) MSEs of $\Delta \mathrm{Hb}{\mathrm{O}}^{\mathrm{gm}}$ , (b) and (d) MSEs of $\Delta \mathrm{Hb}{\mathrm{R}}^{\mathrm{gm}}$ , and (a) and (b) standard artifact case. Changes in the superficial layer is $1∕10\text{th}$ of changes in the gray-matter layer. (c) and (d) Large artifact case. Changes in the superficial layer is $1∕5$ of changes in the gray-matter layer.

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:

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

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.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\mathrm{mm}$ . The wavelengths of the light source were 800 and $840\mathrm{nm}$ . The source-detector separation was assumed from $10\text{to}40\mathrm{mm}$ every $5\mathrm{mm}$ .

Fig. 1

Geometry of a five-layer adult head model for Monte Carlo simulation.

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\mathrm{nm}$ for each layer were chosen from reported data.¹³ $\mathrm{Hb}{\mathrm{O}}_2$ and HbR concentrations were calculated from the absorption coefficients at $800\mathrm{nm}$ with the assumption that oxygen saturation is 70%. The absorption coefficient calculated from water absorption coefficients given in the paper¹⁴ was also used in this case.¹⁵ These values were used to calculate the absorption coefficients at $840\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\mathrm{nm}$ were set to be 95% of those at $800\mathrm{nm}$ .

Table 1

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

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

Fig. 2

Model time course of hemodynamic change: (a) Stimulation paradigm, (b) normalized hemodynamic impulse response, and (c) evoked hemodynamic change.

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

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

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.

Fig. 3

Partial path length of superficial and gray-matter layer based on the standard model parameters. Solid lines indicate the partial pathlength of ${\lambda }_1$ , and dashed lines indicate the partial pathlength of ${\lambda }_2$ .

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.

Fig. 4

Relative (percent) change of partial path length in superficial and gray-matter layers, where model parameters were increased by 10% from the standard values: (a) Cases in which the transport scattering coefficient of each layer was increased and (b) cases in which the thickness of each layer was increased. $X$ -axis corresponds to the source-detector separation.

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\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\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∕5\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}}$ .

Fig. 6

(a, b) Hemoglobin changes estimated by two weightings ( $W_{\mathrm{opt}}$ and $L_{\mathrm{std}}^{-1}$ ) and dual-detector NIRS algorithm. (c) Hemoglobin changes estimated by an ordinary single detector NIRS algorithm. (d–f) Normalized hemoglobin changes of (a–c). Ten different partial path-length matrices were used in the simulation, and they were randomly selected from the test data. In the test data generation, standard deviation of model parameter distribution was assumed to be 20% of the standard values. The estimated hemoglobin changes are shown as light red (HbO) and light blue (HbR) solid lines. The assumed original hemoglobin changes are also shown in the figures.

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\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.²⁰

Next, we normalized the scale of each response in Figs. 6a, 6b, 6c by making its value at $t=15\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.