2.1.

Material and Structure of Phantom

We used reports by Fukui et al.⁴⁶ and Wang et al.⁴⁷ to set the absorption coefficient (${\mu }_a$) and reduced scattering coefficient (${\mu }_s^{\prime }$) values of the scalp and skull and the gray matter of a human adult head to the range of 0.012 to 0.036 and 0.73 to $2.3{\mathrm{mm}}^{-1}$, respectively. In this study, the ${\mu }_a$ and ${\mu }_s^{\prime }$ values of the two layers of absorber were chosen from around this range. A white polyoxymethylene (POM) was used as the base material of the phantom because it has a comparable ${\mu }_s^{\prime }$ value ($0.9{\mathrm{mm}}^{-1}$) to biological tissue and is commercially available and easily worked, although it does not absorb light very much (${\mu }_a=0.002{\mathrm{mm}}^{-1}$). Three absorber materials with ${\mu }_a$ values of 0.01, 0.02, and $0.03{\mathrm{mm}}^{-1}$ were thus prepared. The ${\mu }_s^{\prime }$ value of the absorbers was $0.9{\mathrm{mm}}^{-1}$ to match that of the POM. The absorbers were made by mixing an epoxy resin (MY753, Aeropia Chemical Supplies) and a hardener (XD716, Aeropia Chemical Supplies) at a $3:1$ weight ratio.⁴⁸ An infrared dye, Projet 830 (Avecia Ltd.), and titanium dioxide (SuperWhite, Alec Tiranti Ltd.) were mixed into the absorber materials to control the absorption coefficient and the reduced scattering coefficient of the absorbers, respectively.

Using these materials, a phantom with two absorber layers (upper and lower) in low-absorption and high-scattering base material was developed. The phantom size was $100\times 120(\text{horizontal})\times 60(\text{vertical})\mathrm{mm}$. The absorbers ($75\times 60\times 5\mathrm{mm}$) were held by absorber-holing frames (made of white POM) and were inserted into two cavities ($100\times 70\times 5\mathrm{mm}$) at depths of 2 to 7 mm and 12 to 17 mm. The hexagon socket set screw was used to attach the absorber to the absorber-holding frames. These frames, along with an overall view of the developed dynamic phantom, are shown in Fig. 1. One absorber-holding frame was attached to each automatic stage with an aluminum spacer between the aluminum frame and the stage.

Fig. 1

Photos of the developed phantom. (a) Upper layer and (b) lower layer of absorber-holding frame. (c) Overall view of dynamic phantom.

In our phantom model, to uphold the assumption that an absorber’s position change in a layer region is equivalent to a homogeneous ${\mu }_a$ change in the corresponding layer region, the thickness of each absorber (and absorber-holding frame) should be sufficiently larger than photon-transport length scales (approx. $1/{\mu }_s^{\prime }$). The depth (12 mm) and thickness (5 mm) of the lower layer was chosen from around the values of reported scalp-cortex distance (12 to 18 mm) and thickness of gray matter (4 mm).^49,50 The thickness of the upper layer (5 mm) was also chosen from the reported value of scalp thickness (3.0 to 7.1 mm).^49,50 The upper-layer simulating scalp should be as shallow as possible, but we put a 2-mm-thick superficial layer (scatterer) over the upper layer to strengthen the structure of the phantom. In this way, we determined the depth and thickness of each layer not simply based on anatomical considerations but also on the design of a phantom for industrial purposes.

2.2.

Driving Mechanism of Absorbers

To dynamically change the position of absorbers in the phantom, two layers of absorbers were held in place by absorber-holding frames using two one-axis automatic stages (SGSP20-85, Sigma Koki Co., Ltd.). The stage and absorber-holding frame was connected to an aluminum frame. The two stages were controlled by a PC that transmitted commands to the stage controller via RS232C. Each absorber-holding frame was thus independently controlled by a corresponding one-axis automatic stage.

The automatic stages and the phantom were attached to an aluminum bottom plate. The structure of the phantom and the driving mechanism is shown in Fig. 2. Figure 2(a) shows a top view of the phantom showing one source point, the phosphor detector point, an outline of the “banana” width of the photon spatial distribution, and one of the absorber slabs (superimposed). Figure 2(b) shows a front view of the phantom. Each absorber-holding frame is independently controlled by a corresponding one-axis automatic stage.

Fig. 2

Structure of phantom and driving mechanism of absorbers. (a) Top view. (b) Front view. Each absorber-holding frame is independently controlled by a corresponding one-axis automatic stage.

Two absorbers were placed on the opposite side from the center line (scanning line) so that sensitivity to the absorption change of each layer would not be much influenced by the absorption change of another layer. The opposite-side placement minimizes the number of photons that propagates in both absorbers, whereas same-side placement maximizes it. Therefore, the opposite-side placement of absorbers minimizes any interaction effects between two absorbers that could lead to nonlinearity between the position and the $\mathrm{\Delta }\mathrm{OD}$.

The target position was converted to the stroke in stage coordinate value and sent as a command from the PC to the stage controller. The origin and positive or negative direction of the position coordinate axes were set so that the partial optical path length in the absorber was increased (i.e., the optical density was increased) when the stage position was positive in the defined coordinate axis (hereafter, “$x$” used for the place coordinates for both layers). Therefore, the directions of the position coordinates for the upper and lower layers were set to be opposite. Optical density change (i.e., effective $\mathrm{\Delta }{\mu }_a$) was controlled by the partial effective path length in absorbers that was determined by the stage-position change ($\mathrm{\Delta }x$), not by the $\mathrm{\Delta }{\mu }_a$ of the absorber.

4.1.

Calibration Method

A calibration curve of the relationship between the stage position and changes in optical density ($\mathrm{\Delta }\mathrm{OD}$) is necessary in order to use the phantom for NIRS measurement. To calibrate the phantom, its $\mathrm{\Delta }\mathrm{OD}$ was measured as the inserted absorber at either of the upper or lower layers driven by a one-axis automatic stage. The layer that is not driven is replaced by base material (white POM) by setting the stage positions to $x=-37$ and $-43\mathrm{mm}$ for the upper and lower layers, respectively. $\mathrm{\Delta }\mathrm{OD}$ was calculated by

4.2.

Relationship Between Stage Position and $\mathrm{\Delta }\mathrm{OD}$

To ensure a linear relationship between the stage-position change ($\mathrm{\Delta }x$) and $\mathrm{\Delta }\mathrm{OD}$ for both layers, and to minimize the interaction effect between two absorber layers, the absorption coefficient of the upper layer should be smaller. Furthermore, in terms of phantom design, to make the $\mathrm{\Delta }\mathrm{OD}$ values at a typical S-D distance (for example, 29.5 mm) caused by the movements of two absorber layers as equal as possible, the absorption coefficient of the lower layer should be higher than that of the upper layer. Among all combinations of ${\mu }_a$ values, we selected 0.01 and $0.03{\mathrm{mm}}^{-1}$ for the upper and lower layers, respectively.

The relationship between stage position (the coordinates of which are defined in Sec. 2.2) and $\mathrm{\Delta }\mathrm{OD}$ for 10 different S-D distances is shown in Fig. 4(a) and 4(b). Figure 4(a) describes $\mathrm{\Delta }\mathrm{OD}$ changes as only the upper layer being driven (${\mu }_a$ of absorber: $0.01{\mathrm{mm}}^{-1}$), and Fig. 4(b) describes $\mathrm{\Delta }\mathrm{OD}$ changes as only the lower layer being driven (${\mu }_a$ of absorber: $0.03{\mathrm{mm}}^{-1}$). In Fig. 4(a), the inflection point of $\mathrm{\Delta }\mathrm{OD}$ is deviated to the minus side because the photon distribution is changed due to the movement of the upper-layer absorber. In contrast, the deviation from the center line of the inflection point of $\mathrm{\Delta }\mathrm{OD}$ for the lower layer is relatively smaller than that for the upper layer because the ratio of partial effective path length in the lower absorber to total path length is smaller than that in the upper absorber, and the influence of stage position on the photon distribution in the phantom is relatively small.

Fig. 4

Calibration results. Relationship between stage position and $\mathrm{\Delta }\mathrm{OD}$ for 10 conditions of S-D distance when (a) only upper layer (${\mu }_a$ of absorber: $0.01{\mathrm{mm}}^{-1}$) was driven and (b) only lower layer (${\mu }_a$ of absorber: $0.03{\mathrm{mm}}^{-1}$) was driven. Gradient of $\mathrm{\Delta }\mathrm{OD}$ at each stage position for 10 conditions of S-D distance when (c) only upper layer (${\mu }_a$ of absorber: $0.01{\mathrm{mm}}^{-1}$) was driven and (d) only lower layer (${\mu }_a$ of absorber: $0.03{\mathrm{mm}}^{-1}$) was driven.

Moreover, in both conditions (a) and (b), the initial rise of $\mathrm{\Delta }\mathrm{OD}$ while increasing the stage position from $x=-30\mathrm{mm}$ is earlier for the longer S-D distance because the spatial distribution of photons in the phantom is larger for the longer S-D distance probe, and the probe is sensitive to the absorber that is further from the center line ($x=0$ in the position coordinate).

To quantify the deviation of the inflection point of $\mathrm{\Delta }\mathrm{OD}$ from the center point, we calculated the gradient of $\mathrm{\Delta }\mathrm{OD}$. The relationship between the stage position and the stage-axial spatial gradient of the $\mathrm{\Delta }\mathrm{OD}$ [$\partial (\mathrm{\Delta }\mathrm{OD})/\partial x$] at each S-D distance is shown in Fig. 4(c) and 4(d). Figure 4(c) describes $\partial (\mathrm{\Delta }\mathrm{OD})/\partial x$ changes as only the upper layer being driven (${\mu }_a$ of absorber: $0.01{\mathrm{mm}}^{-1}$) and Fig. 4(d) describes $\partial (\mathrm{\Delta }\mathrm{OD})/\partial x$ changes as only the lower layer being driven (${\mu }_a$ of absorber: $0.03{\mathrm{mm}}^{-1}$). The Gaussian curves fitted to the $\mathrm{\Delta }\mathrm{OD}$ gradient of the measurement data are shown in the figures. The mean stage positions at the maximal $\mathrm{\Delta }\mathrm{OD}$ gradient were $x=-1.9\mathrm{mm}$ (upper layer) and $x=-0.7\mathrm{mm}$ (lower layer). The mean values and standard deviations of full width at half maximum (FWHM) of the $\mathrm{\Delta }\mathrm{OD}$ gradient were $9.4\pm 2.3\mathrm{mm}$ (upper layer) and $19\pm 1.5\mathrm{mm}$ (lower layer). The minimum FWHM values were 6.0 mm at 9.5-mm S-D distance and 16.8 mm at 13.5-mm S-D distance for the upper and lower layers, respectively. The $\mathrm{\Delta }\mathrm{OD}$ gradient at the 5.5-mm S-D distance was not used to calculate the mean position or FWHM for the lower layer because the $\mathrm{\Delta }\mathrm{OD}$ gradient was almost equal to zero.

To proportionally change the $\mathrm{\Delta }\mathrm{OD}$ by the change of the stage position, the stage should be driven within the range where the spatial slope of $\mathrm{\Delta }\mathrm{OD}$ is constant. When an approximate-linearity of the relationship between changes in stage position and $\mathrm{\Delta }\mathrm{OD}$ was defined as the stage positions where the change in $\mathrm{\Delta }\mathrm{OD}$ gradient is within 10% from the maximum (at the inflection point), the calculated ranges were $x=-1.9\pm 1.8$ and $x=-0.7\pm 3.6\mathrm{mm}$ for the upper and lower layers, respectively. The ranges were calculated by using the mean FWHM values of the $\mathrm{\Delta }\mathrm{OD}$ gradient as representative values. The waveform of $\mathrm{\Delta }\mathrm{OD}$ was approximately the same as that of the stage position when the stages were driven within the range; therefore, the $\mathrm{\Delta }\mathrm{OD}$ change profile can be easily designed. The waveforms of the driving stages used for the following experiments were designed so that the stage positions were all within the above range.

4.3.

Comparison with Results by Monte Carlo Simulation

To investigate the consistency between the developed phantom and a simulation model where the same ${\mu }_a$ and ${\mu }_s^{\prime }$ were assumed, $\mathrm{\Delta }\mathrm{OD}$ measured while each layer region was switched from uniform POM to uniform absorber ($\mathrm{\Delta }{\mathrm{OD}}_{\mathrm{sw}}$) was compared with the linear-fitted partial effective path length ($L_{\mathrm{eff}}$) of photons calculated by a three-dimensional (3-D) Monte Carlo simulation^53–55 in the following way.

From experiments, $\mathrm{\Delta }{\mathrm{OD}}_{\mathrm{sw}}$ was obtained by the $\mathrm{\Delta }\mathrm{OD}$ measured at the 30-mm position where each layer region underneath the measurement optodes was switched from uniform POM to uniform absorber, while another layer was uniform POM. Monte Carlo simulation was used to calculate the $L_{\mathrm{eff}}^{\text{upper}}$ and $L_{\mathrm{eff}}^{\text{lower}}$ values by using absorption coefficients $({\mu }_a^{\text{upper}},{\mu }_a^{\text{lower}})=(0.01,0.002)$ and $(0.002,0.03)({\mathrm{mm}}^{-1})$, respectively. When the change in the product of absorption coefficient and partial effective path length is assumed to be equal to $\mathrm{\Delta }\mathrm{OD}$,⁵⁴ $\mathrm{\Delta }\mathrm{OD}$ is expressed as:

The data obtained at the longest S-D distance (41.5 mm) was not used for the fitting because of a low signal-to-noise ratio. The obtained scaling factors were 0.010 and $0.023{\mathrm{mm}}^{-1}$ for upper- and lower-layer $L_{\mathrm{eff}}$, respectively. The experimental $\mathrm{\Delta }{\mathrm{OD}}_{\mathrm{sw}}$, the computed $\mathrm{\Delta }\mathrm{OD}$ (linear fit of the simulated $L_{\mathrm{eff}}$) of each layer, and the ratio of $\mathrm{\Delta }\mathrm{OD}$ at each S-D distance are plotted in Fig. 5(a). The goodness of fit between the experimental and computed $\mathrm{\Delta }\mathrm{OD}$ was tested using the chi-square test. Results showed that the validity of the model (spline function from computed $L_{\mathrm{eff}}$), including S-D distance dependency of $\mathrm{\Delta }\mathrm{OD}$ for each layer, was warranted [reduced chi-square: ${\chi }^2=0.0003$ (upper), 0.017 (lower) $<1$]. According to Eq. (2), when $L_{\mathrm{eff}}$ values are constant and the $\mathrm{\Delta }{\mu }_a$ of either the upper or lower layer is equal to zero, the obtained scaling factor should be comparable to $\mathrm{\Delta }{\mu }_a$. The scaling factors for the upper and lower layers are therefore expected to be 0.008 (0.010 to 0.002) and 0.028 (0.030 to 0.002) ${\mathrm{mm}}^{-1}$, respectively, which are the values of the assumed absorption-coefficient changes. The differences between the scaling factor and acual experimental $\mathrm{\Delta }{\mu }_a$ were 20% and 22%.

Fig. 5

(a) Experimental $\mathrm{\Delta }{\mathrm{OD}}_{\mathrm{sw}}$, the computed $\mathrm{\Delta }\mathrm{OD}$ (linear fit of the simulated $L_{\mathrm{eff}}$) of each layer, and the ratio of $\mathrm{\Delta }\mathrm{OD}(R_{\mathrm{\Delta }\mathrm{OD}})$ at each S-D distance. (b) Experimental (exp.) and simulated (sim.) $L_{\mathrm{eff}}$ of upper and lower layers. Experimental $L_{\mathrm{eff}}$ was calculated using $\mathrm{\Delta }{\mathrm{OD}}_{\mathrm{sw}}$ and $\mathrm{\Delta }{\mu }_a$. When $\mathrm{\Delta }{\mathrm{OD}}_{\mathrm{sw}}$ was measured, absorption coefficients of upper and lower layers were changed from 0.002 to $0.01{\mathrm{mm}}^{-1}$ and from 0.002 to $0.03{\mathrm{mm}}^{-1}$, respectively. In the Monte Carlo simulation, absorption coefficients of the upper and lower layers were set to 0.01 and $0.03{\mathrm{mm}}^{-1}$, respectively.

The $L_{\mathrm{eff}}$ values computed by Monte Carlo simulation and the experimental $L_{\mathrm{eff}}$ values ($=\mathrm{\Delta }{\mathrm{OD}}_{\mathrm{sw}}/\mathrm{\Delta }{\mu }_a$) are shown in Fig. 5(b). Experimental $L_{\mathrm{eff}}$ values were calculated using expected $\mathrm{\Delta }{\mu }_a$ values. The differences between the $L_{\mathrm{eff}}$ values from the measurement and the simulation were, on average, 19.0% and 9.4% ($L_{\mathrm{eff}}$ at a S-D distance of 5.5 mm for the lower layer was not used because the $L_{\mathrm{eff}}$ from the simulation was almost equal to zero, so the error was 240%). The difference of the values between experiment and simulation could have been caused by 1. an inaccurate measurement of the absorber’s optical properties, 2. the difference between the photon distributions (and subsequent $L_{\mathrm{eff}}$ change) under the condition when the lower layer was filled with absorber and the condition when it was not, or 3. measurement noise during the experiments.

5.1.

Synthesis of Waveforms

We used the phantom to evaluate the multidistance analytical methods for discriminating the effects of superficial tissue in the following way. The waveform of the absorption change at each layer can be easily reproduced by controlling the 1-axis stage with digital-base commands. While the upper and lower layers were driven under the following two conditions, $\mathrm{\Delta }\mathrm{OD}$ was measured by the multidistance NIRS system described in Sec. 3. Measurements were also conducted while only one layer was driven under the two conditions.

Under condition 1, the stage position of the upper layer was driven so as to simulate white noise while the stage position of the lower layer was driven so as to simulate a change in cerebral blood volume (CBV)⁵⁶ evoked by neural activation. Under condition 2, the stage position of the upper layer was driven so as to simulate a biological fluctuation,^57,58 while the stage position of the lower layer was driven so as to simulate a change in CBV response. All three simulated waveforms could be obtained with the in vivo NIRS measurement.

Three types of time-course data, along with the power spectral densities of the stage-position changes, are shown in Fig. 6(a) white noise, 6(b) biological fluctuation, and 6(c) CBV response. The time-course data were normalized to zero-mean and unit-variance data when the power spectrum densities were calculated. Power spectral densities were calculated on the assumption that one point of time-course data is equivalent to 1 sec. At frequencies lower than 0.1 Hz, the waveforms in 6(b) and 6(c) commonly have a lot of power compared to those at higher frequencies. These two waveforms cannot thus be discriminated with frequency filtering. The correlation coefficients between the waveforms in 6(a) and 6(c) and those in 6(b) and 6(c) were 0.015 and $-0.190$, respectively.

Fig. 6

Time course and power spectrum density (PSD) of stage-position change simulating NIRS signal when one point of time course data was equivalent to 1 sec: (a) white noise, (b) biological fluctuation, and (c) cerebral blood volume (CBV) response evoked by neuronal activation.

5.2.

Extraction of Deep-Layer Signal with Subtraction Method

Using our dynamic phantom, the performance of a subtraction method with short-separation regression,^29–31 which extracts the deep absorption change from a signal that is contaminated by the absorption change of the surface layer, was investigated while both absorbers of the upper and lower layers were independently driven. In the subtraction method, the signal from a short S-D distance channel was linearly fitted to that from a long-distance channel.^29–31 In the present study, we compare the result from the method with the correct waveform that was obtained while only the lower layer was driven. Performance was defined as the amplitude ratio of the subtraction result to actual lower layer change (target).

When the subtraction method is used, lower-layer signal amplitude at a longer S-D distance ($\mathrm{\Delta }{\mathrm{OD}}_{\mathrm{SD}1}^{\text{lower}}$) is also subtracted by a lower-layer signal that is included in the shorter S-D distance signal (regressor channel) ($\mathrm{\Delta }{\mathrm{OD}}_{\mathrm{SD}2}^{\text{lower}}$), where the scaling factor is the ratio of the upper-layer $\mathrm{\Delta }\mathrm{OD}$ of longer S-D distance to shorter S-D distance ($\mathrm{\Delta }{\mathrm{OD}}_{\mathrm{SD}1}^{\text{upper}}/\mathrm{\Delta }{\mathrm{OD}}_{\mathrm{SD}2}^{\text{upper}}$) because subtraction is performed so as to eliminate the upper-layer signal. If the performance of the subtraction method is defined as the amplitude ratio (AR) of the subtraction result to only the lower-layer change (target), AR can be expressed as:

5.3.

Results of Performance Evaluation of Subtraction Method

The $\mathrm{\Delta }\mathrm{OD}$ measured at each S-D distance when the upper and lower layers were driven in white-noise and a CBV waveform, respectively (condition 1), is shown in Fig. 7. Mixed signals of white-noise and CBV waveforms were obtained, but almost the same shape of CBV signal can be observed in $\mathrm{\Delta }\mathrm{OD}$ signals at longer S-D distances. The subtraction results with the short-distance regression of $\mathrm{\Delta }\mathrm{OD}$ measured while both layers were driven under condition 1, along with $\mathrm{\Delta }\mathrm{OD}$ that was measured at a long distance while only the lower layer was singly driven (target), are shown in Fig. 8. The noise that is commonly present in $\mathrm{\Delta }\mathrm{OD}$ signals at all S-D distances was eliminated, and a CBV signal of the lower layer was clearly obtained when the shortest S-D distance (5.5 mm) channel was used.

Fig. 7

$\mathrm{\Delta }\mathrm{OD}$ measured at each S-D distance as upper and lower layers driven in white-noise and CBV waveform, respectively (condition 1). S-D distances: (a) 5.5 mm, (b) 13.5 mm, (c) 25.5 mm, and (d) 37.5 mm.

Fig. 8

Subtraction result with short-distance regression of $\mathrm{\Delta }\mathrm{OD}$ measured while both layers are driven in condition 1, and $\mathrm{\Delta }\mathrm{OD}$ measured at a long distance while only lower layer (target) is driven. Combination of S-D distances: (a) 37.5 and 5.5 mm, (b) 37.5 and 13.5 mm, (c) 37.5 and 25.5 mm, (d) 25.5 and 5.5 mm, (e) 25.5 and 13.5 mm, and (f) 13.5 and 5.5 mm.

The results under condition 2 are shown in Figs. 9 and 10. As seen in Fig. 10, even if the tendency of the CBV response was not clearly visible from the raw $\mathrm{\Delta }\mathrm{OD}$ signal, the CBV signal at the lower layer was obtained when the shortest S-D distance channel was used. The subtraction method was effective even under the condition that the two waveforms had power in the same frequency range (Fig. 6).

Fig. 9

$\mathrm{\Delta }\mathrm{OD}$ measured at each S-D distance as upper and lower layers were driven in waveform of biological fluctuation and CBV, respectively (condition 2). S-D distances: (a) 5.5 mm, (b) 13.5 mm, (c) 25.5 mm, and (d) 37.5 mm.

Fig. 10

Subtraction result with short-distance regression of $\mathrm{\Delta }\mathrm{OD}$ measured while both layers were driven in condition 2, and $\mathrm{\Delta }\mathrm{OD}$ measured at a long distance while only lower layer (target) was driven. Combination of S-D distances: (a) 37.5 and 5.5 mm, (b) 37.5 and 13.5 mm, (c) 37.5 and 25.5 mm, (d) 25.5 and 5.5 mm, (e) 25.5 and 13.5 mm, and (f) 13.5 and 5.5 mm.

Table 1(a)–(c) shows the performance of the subtraction method defined as the amplitude ratios (AR) of the subtraction result to only the lower-layer change (target) under three conditions: (a) AR from computed $\mathrm{\Delta }\mathrm{OD}$, (b) experimental AR under condition 1, and (c) experimental AR under condition 2. The mean error and its standard deviation of experimental AR from computed AR was $23\pm 22\%$, but under four conditions of S-D distance combinations, the errors were less than 10%. For example, when 25.5- and 5.5-mm S-D distances were used under condition 1, the error was only 2.5%. These errors might have been caused by several factors, such as waveform combination or method of calculating AR. We demonstrated that the approximate performance of the subtraction method can be estimated by the computed AR and that a shorter S-D distance is better for a regressor.

Table 1

Performance of subtraction method in extraction of lower-layer signal. Amplitude ratio (AR) of subtraction result to only lower-layer change (target) under three conditions: (a) AR from computed ΔOD; (b) experimental AR under condition 1; and (c) experimental AR under condition 2.

A CBV response with good reproducibility was obtained three successive times. When a 25.5-mm channel was used as a regressor, the lower signal was not accurately obtained (AR values were 0.13 and 0.21 under conditions 1 and 2, respectively) because this channel contained a considerable amount of the lower signal ($R_{\mathrm{\Delta }\mathrm{OD}}=2.5$). In the phantom measurement, up to now, it has been difficult to change the absorption coefficient (or subsequent optical density) with complete reproducibility because conventional dynamic phantoms are based on a two-state control. Our dynamic phantom, in contrast, can be applied to a variety of signal discrimination methods.

5.4.

Signal Discrimination with ICA

ICA is a signal discrimination method that extracts independent components from multiple signals without knowledge of the obtained signal by utilizing the statistical independence of the source components.⁵⁹ This method is effective for analyzing signals that have multiple signal sources and need to be measured at multiple points.^27,28

ICA can be used to extract multiple independent components and weights at each measurement position. Original data can be reconstructed by totaling the products of independent components and weights. In this study, we used the time-delayed decorrelation (TDD)-ICA algorithm, which assumes that the time-delayed cross-correlation between independent components can vanish at any time.^28,60–62 Multiple signals used in this analysis are, in this case, the time series recorded at multiple S-D distances.

We tested whether or not we could discriminate between absorption changes of the upper and lower layers from the obtained $\mathrm{\Delta }\mathrm{OD}$ time course data with ICA when the absorption of both layers of the phantom were simultaneously changed (condition 2, stated above). The obtained components were compared with the $\mathrm{\Delta }\mathrm{OD}$ time-course data obtained when the absorption of only one layer was driven (actual signals). Performance was quantitatively evaluated on the basis of the correlation coefficient between the obtained component and the actual signal.

5.5.

Results of Performance Evaluation of ICA Method

Using the data obtained under condition 2 (both layers driven; upper layer: biological fluctuation, lower layer: CBV), ICA was used to extract independent components and their weights at each S-D distance. Figure 11(a) plots the time course of independent components extracted with ICA and Fig. 11(b) plots weights of each component at each S-D distance. The top-two largest components in the contributing ratio were extracted. Components 1 and 2 correspond to the absorption change of the lower and upper layers, respectively. The correlation coefficients between the independent components (extracted when both layers were driven) and the $\mathrm{\Delta }\mathrm{OD}$ time course (obtained at a 25.5-mm S-D distance in advance while only each layer was driven) were 0.988 (upper layer) and 0.996 (lower layer). (The channel of the 41.5-mm S-D distance was not used because of a low signal-to-noise ratio.)

Fig. 11

ICA results. (a) Independent components extracted with ICA and (b) weights for each component at each S-D distance.

The dependence of a component’s weight on S-D distance was similar to that of $\mathrm{\Delta }\mathrm{OD}$ and $L_{\mathrm{eff}}$ (see Fig. 5). The reasons for the slight difference between these dependences are 1. the interference of one layer absorption change to the sensitivity to that of other layer’s absorption change, 2. the correlation between two-layer waveforms, and 3. the measurement noise.

According to Eq. (2), $\mathrm{\Delta }\mathrm{OD}$ is the product of the change in absorption coefficient and partial effective path length, and therefore the weights at each S-D distance in Fig. 11(b) reflect the changes in partial effective path length depending on the S-D distance shown in Fig. 5. This result shows that the phantom simulates the absorption change in the entire layer by sliding the position of the absorber (the coefficient of which is fixed).

The effectiveness of ICA for discriminating between surface- and deep-layer effects was experimentally demonstrated under a condition in which the time-course changes in the absorption of the upper and lower layers were different (preferably, statistically independent). Such validation is difficult to ensure with human measurement because the exact amount of absorption change cannot be precisely known.

6.1.

Summary

We developed a dynamic phantom for simulating tissue with two absorption layers that are independently driven with two one-axis automatic stages that can cause the absorption change of any waveform. When several simulated NIRS signals—such as white noise, physiological fluctuation, and CBV response—were created by the driving stages, the time-course data of $\mathrm{\Delta }\mathrm{OD}$ were obtained with a good reproducibility and linear relationship to stage-position changes. These results show that $\mathrm{\Delta }\mathrm{OD}$ (i.e., $\mathrm{\Delta }{\mu }_a$) is caused by the stage-position change ($\mathrm{\Delta }x$).

We used this phantom to evaluate the performance of a subtraction method with a short-distance regressor. When a short S-D distance channel that is not sensitive to lower-layer absorption change is used as a regressor, the signal that resulted from subtraction was in better agreement with the actual lower signal. Furthermore, to demonstrate the advantage of simulating absorption coefficient waveforms with our phantom, we applied an ICA method to the signals obtained while two layers were concurrently driven with different time courses. When the extracted components were compared to the signal obtained when the absorber of each layer was singly driven, the correlation coefficients were over 0.98 for both layers.

These results demonstrate that our dynamic phantom can be used to evaluate methods for discriminating between the effects of scalp and cerebral blood flow in NIRS signals.

6.2.

Future Work

It has been reported that systemic fluctuations related to blood pressure or heart rate—in particular, because of dynamic cerebral autoregulation^63,64 or veins on the brain surface³⁴—can influence the NIRS signal.^65,66 In these cases, to extract cerebral blood including systemic signal, we need to use a signal discrimination method that considers the dependence of signal amplitude on S-D distance instead of just the waveform characteristics of the signal. The dynamic phantom we developed in the present study should contribute to a direct validation of other methods to extract desired signals (e.g., a CBV signal) from mixed signals. It should be noted that, however, the results from this phantom cannot be easily generalized because the structure and optical properties of it are not the same as human. If this phantom model is extended to more realistic geometries with less spatial homogeneity, the result from the phantom would be more reliable and would contribute to the signal discrimination. In particular, the effects of correlation or covariance, frequency characteristics, and waveforms on signal-processing methods can be directly and quantitatively investigated.