Estimating the absorption coefficient of the bottom layer in four-layered turbid mediums based on the time-domain depth sensitivity of near-infrared light reflectance

Abstract. Expanding our previously proposed “time segment analysis” for a two-layered turbid medium, this study attempted to selectively determine the absorption coefficient (μa) of the bottom layer in a four-layered human head model with time-domain near-infrared measurements. The difference curve in the temporal profiles of the light attenuation between an object and a reference medium, which are obtained from Monte Carlo simulations, is divided into segments along the time axis, and a slope for each segment is calculated to obtain the depth-dependent μa(μaseg). The reduced scattering coefficient (μs′) of the reference is determined by curve fitting with the temporal point spread function derived from the analytical solution of the diffusion equation to the time-resolved reflectance of the object. The deviation of μaseg from the actual μa is expressed by a function of the ratio of μaseg in an earlier time segment to that in a later segment for mediums with different optical properties and thicknesses of the upper layers. Using this function, it is possible to determine the μa of the bottom layer in a four-layered epoxy resin-based phantom. These results suggest that the method reported here has potential for determining the μa of the cerebral tissue in humans.


Introduction
In 1977, Jöbsis 1 was the first to report near-infrared measurements of tissue oxygenation in living animals, and near-infrared spectroscopy (NIRS) has been developing into a useful tool for monitoring of cerebral hemodynamic changes.Today, NIRS is used in a wide variety of research fields and in clinical medicine, while quantitative and selective extraction of signals originating in the brain remains difficult.Recently, there have been reports that the skin blood flow influences NIRS signals, an issue that cannot be ignored in NIRS studies. 2,3Over the past 35 years, the NIRS community has focused on solving these issues and has developed a wide range of types of NIRS measurements.][6][7] Propagation of photons in living tissue is described by the radiative transport equation (RTE); however, the RTE, an integro-differential equation, is not easy to solve even by numerical methods.Therefore, a diffusion approximation to the RTE is widely used instead, since light propagation can be con-sidered isotropic at the macroscopic level, when the scattering interaction is much larger than the absorption and the point of interest is far from light sources and boundaries. 83][14] Other methods have also been proposed to estimate the optical properties of layered mediums, such as the multilayered (time-dependent) DE, [15][16][17] also a method with spatially and time-resolved reflectance, 18 and a method considering time-dependent mean partial path lengths. 19However, none of these methods have been applied to practical human head measurements.
Diffuse optical tomography (DOT) is a technique for reconstructing images of optical properties using multiple light sources and detectors 20 and is thought to be the most promising technique for the quantitative detection of focal changes in cerebral hemodynamics.DOT can be performed with timedomain, 21 frequency-domain, 22 and continuous-wave (CW) measurement instruments. 23Recently, CW high-density DOT with high spatial resolution has been developed and applied to functional hemodynamic maps of the adult human visual cortex. 24However, this CW-DOT only provides qualitative images of the measured changes as it is based on linear image reconstruction.Nonlinear iterative reconstruction schemes are required to obtain quantitative images, and such procedures are still under development. 25e have previously proposed a novel approach, "time segment analysis," for determination of the μ a of each layer in a two-layered turbid medium with a single-channel TRS instrument. 26In that approach, first a time-attenuation difference curve, derived by the subtraction of the time-resolved reflectance of a reference medium from that of an object, was divided into segments along the time axis (e.g., 500 ps time width).Then, a depth-dependent μ a (time-segmented μ a ) was estimated from the slope of each segment.The difference between the estimated μ a and the actual μ a in the lower layer then depends on the ratio of the actual μ a in the upper layer to that in the lower layer, which was linearly related to the ratio of the estimated μ a in an early time segment to that in a later segment.Using this relationship, it was possible to correct the μ a obtained from the later time segment to derive the μ a of the lower layer in twolayered mediums.An advantage of this approach is that it involves a single-distance measurement, which is more suitable than a multidistance measurement for measurements of the head with its heterogeneous superficial layer and curvature.A single distance measurement is also straightforward to expand into multichannel measurements.To apply our method to human head measurements, the validity and reliability of the method must be confirmed with three or more layered models and with other models with more complex structures.In the present study, the investigation focuses on the applicability of our method 26 to a four-layered slab model with a Monte Carlo simulation and phantom experiments.
The article here first reports the results of simulations: (1) time-segmented μ a under several conditions in which μ a and the reduced scattering coefficient (μ s 0 ) of each layer were changed, (2) the difference between the time-segmented μ a at a later time and the actual μ a of the bottom layer, and its relation to the ratio of the time-segmented μ a in an earlier time segment to that of a later one, and (3) the mean partial path length of photons detected in each time segment to be able to interpret the time-segmented μ a values.Then, the experimental results are reported: (4) the determination of μ a of the bottom layer of a four-layered epoxy resin-based phantom with the "time segment analysis." 2 Methods for Determining the μ a of the Bottom Layer The details of the method have been described in a previous paper, 26 and the following is a brief summary.

Theoretical Basis
Two kinds of mediums, a reference medium and an object medium for the measurements, are considered.When a light impulse is irradiated on the surface of the mediums, the reflected light intensity at time t for the reference and the object is expressed by I R ðtÞ ¼ S R ðtÞ expð−μ a;R • c R • tÞ and IðtÞ ¼ SðtÞ expð−μ a • c • tÞ, 27,28 where c is the light velocity in a medium and SðtÞ is a scattering function which is dependent on the μ s 0 of the medium in the photon diffusion regime.The light attenuation, A, is defined as the logarithm of the inverse of the reflectance.According to the time-resolved Beer-Lambert law, the difference between the attenuation of an object ½AðtÞ and of a reference ½A R ðtÞ, A diff ðtÞ, can be expressed by where S diff ðtÞ ¼ ln½S R ðtÞ∕SðtÞ.
In the case of an object that is a layered medium in which μ a varies with layer (i.e., depth) and the reference is homogeneous, A diff ðtÞ is not a linear function of t even though S diff ðtÞ can be ignored and c ¼ c R , that is, μ a is different at each detection time.The method introduces a time-dependent apparent absorption coefficient, μ a ðtÞ, given by P i μ ai l i ðtÞ∕LðtÞ, where the subscript i refers to the i'th layer in the object medium, l i ðtÞ is the stochastic time-dependent partial path length in the i'th layer for the photons detected at time t, and LðtÞ is the total path length of photons detected at time t, which is defined as LðtÞ ¼ ct.Then dA diff ðtÞ∕dt is expressed as follows: where dA diff ðtÞ∕dt expresses the variation of μ a in the depth direction in the object.

Time Segment Analysis
To simplify the process of the analysis, the A diff − t curve within a short time range (e.g., a few hundred picoseconds, Δt) is considered a straight line, and a term for the mean absorption differences, hdA diff ðtÞ∕dti, is introduced, the value of this is obtained by ΔA diff ðtÞ∕Δt, rather than dA diff ðtÞ∕dt.In this study, the A diff − t curve in the range 500 to 3500 ps is divided into six time segments of 500 ps, and the slope of each time segment is estimated (hdA diff ðtÞ∕dti) by linear regression.The mean μ a ðtÞ of each time segment is referred to as the time-segmented μ a , μ seg a .Then, assuming that a reference and an object have the same μ s 0 and refractive index, this gives the following equation, The time-segmented μ a ðμ seg a Þ was calculated by substituting hdA diff ðtÞ∕dti, μ a ; R, and c in Eq. (3).
The segments of the A diff − t curve were numbered as K ¼ I; II; : : : ; VI in time order, so the segment of a time range from 500 to 1000 ps is referred to as segment-I and that from 2500 to 3000 ps was as segment-V, μ seg a (K) denotes the μ seg a in segment K. To calculate the reflectances from four-layered (object) and from homogeneous (reference) mediums, we used the Monte Carlo code developed by Wang and Jacques. 29The code was modified to fit our measurement system.Source photons were perpendicularly irradiated onto the surface of the semi-infinite mediums, and photons that were emitted from a detector position on the surface were all detected with a time step of 10 ps, and the source-detector distance was 30 mm.Reflectance profiles simulated in an anisotropic setting (scattering coefficient μ s ¼ 10 mm −1 , anisotropy parameter g ¼ 0.9), were compared with those in an isotropic setting (μ s ¼ 1.0 mm −1 ) and it was confirmed that the differences between the two were small enough to be disregarded in these simulating parameter conditions.Thus, an isotropic scattering was assumed, reducing the calculation time, and the calculations were repeated until the number of detected photons reached 1,000,000.A light velocity of 0.219 mm∕ps (refractive index 1.37) was used in all the calculations.

Simulation and Experiment
In the previous study, 26 we confirmed that temporal profiles of detected light intensity in homogeneous and two-layered phantoms, which consisted of two layers with optical properties the same as those of the homogeneous one, were almost in agreement.Thus, the effect of the interface between two layers on the pattern of photon propagation was considered to be negligible in the Monte Carlo simulation and phantom measurement under present study conditions.
The statistical error of the Monte Carlo simulation was estimated by fitting a theoretical TPSF for a semi-infinite homogenous medium (based on the solution of the DE derived by employing the extrapolated boundary condition of the reflectance mode [30][31][32] ) to the TPSF obtained from the Monte Carlo simulation (time range of 0 to 4500 ps).The photon diffusion coefficient was expressed by D ¼ 1∕3μ s 0 . 33In the fitting procedure, a weighted least squares fitting method based on the Levenberg-Marquardt Method, 34,35 which is an iterative improved method, was used.In case of practical TRS measurements of inhomogeneous objects, such as the human head, when the chi square nu (χ 2 ν ) value is 0.7 to 1.3, it is considered that the theoretical TPSF is well fitted to the observed one. 36,37In the present study, it was assumed that the best fitting TPSF was the reference true solution, and that the χ 2 ν value mainly represented statistical bias and statistical variations that were attributed to the Monte Carlo simulation.The χ 2 ν values in all the mathematical models described in Sec.3.2 (Table 1) were 1.09 AE 0.05 (mean AE SD).Thus, it was concluded that the statistical error of the Monte Carlo simulation was negligible.

Four-Layered Models for Simulation
We used a four-layered semi-infinite model as a human head model.The layers are numbered i ¼ 1; 2; 3, and 4 in order of depth from the outside, and denoted by #i: the layers #1 to #4 correspond to the scalp, skull, cerebrospinal fluid (CSF), and brain, respectively.The μ s 0 , μ a , and the thickness d of each layer are denoted by μ si 0 , μ ai , and d i , respectively.Table 1 shows the values used in this paper.The optical properties of each layer were chosen from the published data. 7,38,39][42] The models were classified into three groups according to the μ s1 0 and μ s2 0 values: 1.3 and 1.0 mm −1 for Group A, 1.1 and 0.8 mm −1 for Group B, and 1.8 and 1.2 mm −1 for Group C. The μ s3 0 and μ s4 0 values were set to 0.1 and 1.7 mm −1 in all three groups.In Group A, μ a was varied in the range of 0.013 to 0.02 mm −1 in layer #1, 0.011 to 0.014 mm −1 in layer #2, and 0.01 to 0.02 mm −1 in layer #4 separately.In Groups B and C, only μ a1 and μ a4 were varied and μ a2 was set to 0.012 mm −1 .In all the models, μ a1 was set to be larger than μ a2 as the μ a of the scalp is generally larger than that of the skull in adult heads, [43][44][45] and μ a3 was 0.0033 mm −1 .The thicknesses of layers #1 and #2, (d 1 ; d 2 ) (mm), were set to (3, 5), (3, 7), (4, 7), (5, 7), (5, 5), or (5, 9) in Group A, (3, 7) or (5, 7) in Group B, and (5, 7) in Group C, whereas d 3 and d 4 were set as 1 and 90 mm, respectively.We mainly evaluated temporal variations in μ seg a and differences between μ seg a and the actual μ a for the Group A models.The models in Groups B and C were used to investigate the effects of the μ s 0 variation in the upper layers on μ seg a .The μ s 0 of the homogeneous reference was changed within the range of 1.1 to 1.8 mm −1 as described below, while the μ a and the thickness were set to 0 mm −1 and 100 mm, respectively.

Determination of μ s 0 for the Reference Mediums
In the previous study, 26 it was found that the results of our method depend on the μ s 0 value of the reference, which was ideally the same as that of the upper layer of a two-layered object medium.In general, the μ s 0 value of a multilayered medium determined by the curve fitting procedure as described in Sec.3.1 is closer to the μ s 0 of the upper layer than that of the lower layer.This can be explained by the fact that the slope of the rising phase of the TPSF is dependent on the μ s 0 value but independent of the μ a value, when the optical properties are in the range of the published data for the biological tissue. 46In this article, the μ s 0 value of the reference was determined by curve fitting with the TPSF derived from an analytical solution of the DE for a semi-infinite homogenous medium to the profile of the time-resolved reflectance of the object medium in the time range of 250 to 1000 ps.The reason why we chose this time range is described in Discussion (Sec.5.2).

Calculation of Mean Partial Path Lengths
As shown in Eqs. ( 2) and ( 3), μ seg a depends on the temporal change in l i ðtÞ in the time range of each segment in each layer.Therefore, to interpret the calculated μ seg a , the mean partial path lengths in each layer for a four-layered medium were also estimated by the Monte Carlo simulation.Two models in Group A [μ a1 ¼ 0.015, μ a2 ¼ 0.012, μ a4 ¼ 0.02, ðd 1 ; d 2 Þ ¼ ð5; 5Þ; and μ a1 ¼ 0.015, μ a2 ¼ 0.012, μ a4 ¼ 0.01, ðd 1 ; d 2 Þ ¼ ð5; 5Þ], were perpendicularly partitioned into 1 mm thick laminae and the attenuation changes (ΔA) caused by altering μ a in a lamina by þ10% of the initial value were calculated.The time-independent mean partial path length of photons in the lamina where μ a is changed can be estimated by the ΔA∕Δμ a of the whole time range of the time-resolved reflectance.The present study estimated a time-dependent mean partial path length given by the time segment unit [denoted as the time-segment dependent mean partial path length (s-mPPL)] at the depth of z to z þ 1 mm (z ¼ 0; 1; 2; : : : ; 19) obtained by dividing the ΔA calculated in each time segment by Δμ a .We also calculated a time-segment dependent mean total path length (s-mTPL) in the medium, the product of the light velocity and the mean time of flight of photons in each time segment.

Epoxy Resin-Based Phantoms
Four-layered (object) and homogeneous (reference) phantoms, 140 × 140 × 63 (50 for the reference) mm, were prepared.The base of the phantoms was epoxy resin (Beuhler, Lake Bluff, Illinois).We adjusted the F s 0 and μ a to 0.1 to 1.77 and 0.0033 to 0.015 mm −1 at 760 nm, respectively, by adding titanium oxide (Tian Kogyo, Ube, Japan) as scatterers and ink (Greenish brown: Chugai Pharmaceutical, Tokyo, Japan) as an absorber.The thickness of each layer and the optical properties of the object are shown in Table 2.These optical properties are theoretical values based on the concentrations of titanium oxide and ink.Unlike the mathematical phantoms, there were no differences in the μ s 0 values of the objects determined in the two time ranges used in the fitting (1.520 mm −1 , 250 to 1000 ps; and 1.502 mm −1 , 0 to 4500 ps).The μ a and μ s 0 of the reference were 0.001 and 1.52 mm −1 , respectively (Table 2).

Instrumentation
A two-channel TRS instrument (TRS-20, Hamamatsu Photonics, Hamamatsu, Japan) was employed.The TRS-20 consists of three pulsed laser diodes with different wavelengths (760, 800, and 830 nm) with durations of around 150 ps [full width at half maximum (FWHM)] at the repetition rate of 5 MHz, a high-speed photomultiplier tube for single photon detection, and a circuit for time-resolved measurements based on a time correlated single photon counting method.Incident light is delivered to the sample through an optical light guide [GI type, 200 μm core diameter, numerical aperture ðNAÞ ¼ 0.25].Reflected light is collected by a fiber bundle (3-mm diameter, NA ¼ 0.26).Minimum data acquisition time is 100 ms.

Time-Resolved Measurements of a Phantom
The incident and detecting light guides were placed on the upper surface of a phantom, separated by 30 mm.Measurements with an accumulated time of 10 s were conducted at five different positions within an area 40 mm from the edge of the phantom to avoid the distortion of photon propagation due to the edge.Specular reflection and light leakage were prevented with a black light-guide holder.Instrument responses of TRS-20 were measured by placing the incident fiber opposite the detecting fiber with a neutral density filter between the two.The instrument response of was around 250 ps FWHM at all wavelengths.
The optical properties of the phantoms were determined by curve fitting with the TPSF derived from an analytical solution of the DE with the extrapolated boundary condition, convoluted by the instrument function to the measured profile of the timeresolved reflectance as described in Sec.3.1.The light velocity in the phantoms was assumed to be 0.19 mm∕ps, corresponding to a refractive index of 1.58.Before the calculation of A diff ðtÞ, the profiles of the time-resolved reflectance measured by TRS were deconvoluted by the instrument function by using Bayesian deconvolution. 47Results

Temporal Variations in μ seg a
Figure 1 shows the temporal variations in μ seg a obtained from the nine models in Group A. Figure 1(a) shows results in the mediums where both μ a1 and μ a2 were larger than μ a4 (squares) and where μ a1 was larger but μ a2 was equal to μ a4 (triangles).Figure 1(b) shows results in the mediums where μ a2 was smaller than FFFF a4 and μ a1 , while μ a1 were larger than (diamonds), equal to (crosses), or smaller than μ a4 (circles).The μ seg a values are expressed as the ratio to μ a4 .When both μ a1 and μ a2 were larger than μ a4 [squares in Fig. 1(a)], the ratios of the μ seg a to the μ a4 ½μ seg a ðkÞ∕μ a4 term were above 1 in all the time segments.The values gradually decreased to become close to 1 as the time segments become later, and reached a value slightly larger than μ a4 .When μ a2 was equal to μ a4 [triangles in Fig. 1(a)], μ seg a ðkÞ∕μ a4 is nearest 1 in the time segments later than IV.For all the models in Fig. 1(a), the upper layer thickness (d 1 ; d 2 ) influenced the values of μ seg a in segments-I and -II but not those in the later segments.
In the models where μ a2 was smaller than μ a4 [Fig.1(b)], the μ seg a ðkÞ∕μ a4 in all the time segments later than segment-II were smaller than 1.When μ a1 was equal to or larger than μ a4 , the μ seg a was almost constant in segments-II to -VI, and a little smaller than the μ a4 values.The values of μ seg a in the later time segments were independent of the value of μ a1 .When both μ a1 and μ a2 were smaller than μ a4 , the μ seg a ðkÞ∕μ a4 values were much smaller than one in all the time segments.The μ seg a decreased between segments I and II, then μ seg a increased slightly in the later time segments, and remained about 20% smaller than μ a4 .For all the models in Fig. 1(b), the differences between μ a4 and μ seg a in the later time segments were larger in the medium with ðd 1 ; d 2 Þ ¼ ð5; 9Þ than in that with (3, 7).Further, the μ seg a values of the medium with (5, 9) were smaller than those of the medium with (3, 7) in all segments.The same tendency of the temporal variations of μ seg a was observed in other models with different sets of μ ai and thicknesses.In the previous study of two-layered models, it was found that the ratio of μ seg a (V) to the actual μ a of the bottom layer was expressed by a function of the ratio of μ seg a in an earlier time segment to that in a later one.In this article we examined the relationship between the ratio of μ seg a ðVÞ to μ a4 ½μ seg a ðVÞ∕μ a4 and the ratio of μ seg a in segment-I to that in segment-V (μ seg a ratio I∕V ) obtained from Group A (Fig. 2).All data obtained from the various mediums with different sets of (d 1 ; d 2 ) were plotted together in the graph.The plot shows that the μ seg a ðVÞ∕μ a4 could be expressed as one polynomial function of the μ seg a ratio I∕V ½μ seg 35 with an adjusted coefficient of determination R 2 ¼ 0.976.The deviation of μ seg a ðVÞ∕μ a4 from the regression curve was less than 6% (mean error was 1.6%).
Figure 3 is a plot of the μ seg a ðVÞ∕μ a4 versus the μ seg a ratio I∕V for Groups B and C, where the μ s 0 values of layers #1 and #2 are different from those in Group A. The solid line represents the regression curve of the relation between μ seg a ðVÞ∕μ a4 and μ seg a ratio I∕V for the Group A results, as also shown in Fig. 2. The data points from Groups B and C are in good agreement with the regression curve obtained from the data sets of Group A (an adjusted R 2 ¼ 0.975).The deviation of μ seg a ðVÞ∕μ a4 in Groups B and C from the regression curve was less than 8% (mean error was 2.7%).This result suggests the conclusion that the relation between μ seg a ðVÞ∕μ a4 and the μ seg a ratio I∕V can be expressed by 35 irrespective of the μ s 0 values of the upper layers.(a) μ seg a in the models where both μ a1 and μ a2 are larger than μ a4 (black filled square, gray filled square; μ a1 ¼ 0.015 mm −1 , μ a2 ¼ 0.012 mm −1 , μ a4 ¼ 0.01 mm −1 ), and μ a1 is larger than and μ a2 is equal to μ a4 (black filled triangle, gray filled triangle;

Temporal Variations in the s-mPPL in Each Layer
in the models where μ a2 is smaller than μ a4 , while μ a1 is larger than (black filled diamond, gray filled diamond; , and smaller than (black filled circle, gray filled circle; μ a1 ¼ 0.015 mm −1 , μ a2 ¼ 0.012 mm −1 , μ a4 ¼ 0.02 mm −1 ) μ a4 .The vertical axis represents the ratio of μ Fig. 3 The ratio of the μ seg a (V) to μ a4 in reference to the ratio of the μ seg a (I) to (V) (μ seg a ratio I∕V ) predicted by Monte Carlo simulation in all the mediums in Groups B and C, where μ s 0 values of layers #1 and #2 were different from those in Group A. Solid line, which is the same as that in Fig. 2, represents the regression curve obtained from the data sets of Group A.
[Fig.4(a)], the s-mPPLs in layers #2 and #4 markedly increased with time, and the s-mPPL in layer #2 was the longest in the latter time segments.The temporal changes in s-mPPL in layer #1, which was longer than that in layer # 4, was small compared to those in layers #2 and #4.The s-mPPL in layer #3, which was the shortest, also showed small increments with time.
In the medium where μ a4 was smaller than μ a2 [Fig.4(b)], the s-mPPL in layer #4 markedly increased with time, with a steeper slope than that in the medium in Fig. 4(a).The s-mPPLs in layer #2 increased with time in time segments I-IV, with amplitudes very similar to those in Fig. 4(a).In time segment-V, however, the s-mPPL in layer #2 did not increase.The temporal variations in the s-mPPLs in layers #1 and #3 were relatively small, especially after time segment-III.

Phantom Experiment
Figures 5(a) and 5(b) show the profiles of the time-resolved reflectances of the object [IðtÞ] and reference [I R ðtÞ] at 760 nm, which were deconvoluted by the incident light, and the A diff − t curve, respectively.The profiles of A diff ðtÞ were reliable in the time range of 0.5 to 3 ns (segments-I to -V).From Fig. 5(b), μ seg a ðIÞ and μ seg a ðVÞ were estimated at each measurement position.Substituting these values into the function 35, μ a4 was estimated to be 0.0156 AE 0.002∕mm (mean AE SD).This value was very close to the actual μ a4 of the object (0.015∕mm).The temporal variations in μ seg a in four-layered mediums (Fig. 1) largely agreed with the previous results for a two-layered model. 26The μ seg a changed with the time segment, reflecting the differences in μ a in the depth direction.The main finding is that the μ seg a values in later time segments exclusively depend on the differences between μ a2 and μ a4 .That is, the μ seg a ðVÞ accurately represents the μ a4 in the mediums where μ a2 was larger than or equal to μ a4 [Fig.1(a)], whereas in mediums where μ a2 was smaller than μ a4 , the μ seg a ðVÞ was smaller than μ a4 even in the later time segments irrespective of μ a1 values [Fig.1(b)].In addition, the dependence of μ seg a on the upperlayer thickness varied with the difference between μ a2 and μ a4 (black versus grey symbols in Fig. 1).As defined in Sec.2.2, μ seg a represents the mean μ a ðtÞ in a segment.The μ seg a value depends on the mean temporal change in l i ðtÞ for each time segment in the individual layers [Eq.( 2)].Therefore, a consideration of the contribution of the mean change in l i ðtÞ to the change in LðtÞ (¼ cΔt; it is constant in all the segments here) in each time segment is helpful to interpret the findings mentioned above.From Figs. 4(a) and 4(b), we can predict the mean temporal change in l i ðtÞ and its amplitude within the segment (500 ps) by interpolation of s-mPPL between the time segments.Here, it should be noted that the temporal patterns of changes in the s-mPPLs provide more important information than the amplitude of the changes as described below.

Discussion
When μ a2 was smaller than μ a4 , s-mPPL in layer #2 increased considerably in all the time ranges [Fig.4(a)].This implies that the contribution of the mean change in l 2 ðtÞ to the change in LðtÞ was high enough to strongly influence the μ seg a even in later time segments.This dominant change in l 2 ðtÞ can explain the finding that the μ seg a at the later time represented values smaller than μ a4 in the mediums where μ a2 is smaller than μ a4 .The dominance of the change in l 2 ðtÞ also accounts for the difference in μ seg a at later time segments between the mediums with ðd 1 ; d 2 Þ ¼ ð3; 7Þ and with (5, 9).The reason why there was little effect of μ a1 on the temporal pattern of changes in μ seg a may be ascribed to the smaller temporal change in s-mPPL in layer #1.
When μ a2 was larger than μ a4 , μ seg a at segments-V and -VI mainly represented the μ a4 [Fig.1(a)].This is explained by the finding that there was no difference in s-mPPL in layer #2 between segment-IV and segment-V, while s-mPPL in layer #4 in segment-V was longer than that in segment-IV [Fig.4(a)], which indicates that the contribution of the mean change in l 2 ðtÞ to the change in LðtÞ in these segments was smaller than that of l 4 ðtÞ.Therefore, it was concluded that the size difference between the μ a2 and the μ a4 is critical for the change in l i ðtÞ in a later time segment, which determines μ seg a values.Like in two-layered models, when the μ a value of the bottom layer was larger than that of the upper layer, the μ seg a ðVÞ was smaller than the μ a4 .With two-layered models, the ratio of μ seg a (IV) or (V) to the real μ a of the lower layer is expressed by the affine function of the μ seg a ratio I∕V , which enables estimation of the μ a of the lower layer. 26The difference between μ seg a ðVÞ and μ a4 in the four-layered models also depends on the μ seg a ratio I∕V , though the μ seg a ðVÞ∕μ a4 was not expressed by an affine function but by a polynomial function of the μ seg a ratio I∕V (Figs. 2 and 3).This relationship between the μ seg a ðVÞ∕μ a4 and μ seg a ratio I∕V was independent of the μ s 0 values and thicknesses of the layers #1 and #2.In addition to the high R 2 values, experimental measurements with the epoxy resin-based phantom support the validity of using this function, although there is no functional explanation why the relationship between μ seg a ðVÞ∕μ a4 and μ seg a ratio I∕V in various object mediums with different conditions of the optical properties and thicknesses of the upper layers can be expressed by a single polynomial function.

Determination of an Adequate Reference in Time Segment Analysis
We have reported that the difference in μ s 0 between the upper layer of the object and the reference significantly influenced the μ seg a in the time segments earlier than 2500 ps. 26However, this difference did not influence the μ seg a after 2500 ps, and the difference in μ s 0 between the lower layer of the object and the reference rarely affected the μ seg a values in all the time segments.The present study used a curve fitting procedure to determine the μ s 0 value of the reference.The effect of the time range of the curve fitting on the μ seg a within the range of 250 to 3000 ps was investigated.When the μ s 0 of the reference was determined by fitting over the wider time ranges, the plots of μ seg a ðVÞ∕μ a4 versus μ seg a ratio I∕V did not converge on a single curve.Therefore, we used the time range of 250 to 1000 ps as the most appropriate range for fitting to determine the averaged μ s 0 of the upper layers.However, there were no differences in μ s 0 values of the epoxy resin-based phantom of the two fitting ranges (250 to 1000 ps and 0 to 4500 ps).In practical measurements, an inappropriate convolution process by the incident light pulse and the time lag of the pulse launching time may induce errors when deriving the μ s 0 , whereas our simulation did not include noise.This may account for the discrepancy between the simulations and the experimental results.

Possibility of Application of the Present Method to Human Head Measurements
As described above, practical measurements are, in general, accompanied by noises arising from various causes, such as a change in coupling between a light guide and the object's surface.In this article, we examined the applicability of the present method to practical measurements by using the four layered slab phantom.The difference between the estimated (0.0156 AE 0.002∕mm) and actual μ a values (0.015∕mm) was maximally about 0.008∕mm (about 5%).Considering the deviation of μ seg a ðVÞ∕μ a4 from the regression curve (Figs. 2  and 3) that was less than 8%, the estimation error is acceptable.Further, since the relationship between the μ seg a ðVÞ∕μ a4 and μ seg a ratio I∕V was independent of the μ s 0 values and thicknesses of the layers #1 and #2 in the range used in this study (μ s 0 of 0.8 to 1.8/mm, thickness of less than 15 mm), it is expected that the applicability of the present approach could be extended to in vivo measurements.
In the adult human head, the total thickness of scalp, skull, and CSF is 11 to 15 mm in the forehead, and about 12 mm in the temporal and occipital areas. 32,33,42Therefore, it is likely that the μ a of the brain can be estimated by the present method on the assumption that the human head is a four-layered slab model.It is however difficult to simply apply this method to measurements in the parietal region, where the total thickness is 15 to 25 mm.Further, the four-layered slab would be too simple for a human head model.As the next step, thus, we will investigate the applicability of our "time segment analysis" to more realistic human head models with the total thickness of the extracerebral tissues in the range of 8 to 25 mm.

Conclusion
This study demonstrated that the "time segment analysis" of time resolved reflectance could be applied to four-layered slab models.The μ seg a ðVÞ∕μ a4 was expressed by a unique function of the μ seg a ratio I∕V on the condition that the thickness of the upper layer was less than 15 mm.Since the μ a of the bottom layer can be determined by this function, the approach reported here has the potential to selectively and quantitatively measure hemoglobin concentrations in the human brain.

3. 1
Monte Carlo Simulations -Calculation of the Time-Resolved Reflectances

4. 2
Differences Between the μ seg a and the Actual μ a in the Bottom Layer

Figures 4 (Fig. 1
Figures 4(a) and 4(b) show temporal variations in the s-mPPL in the four layers.In the medium where μ a4 was larger than μ a2

5. 1
Characterization of μ seg a in View of the Change in s-mPPL in Each Layer

5. 2
Determination of μ a in the Bottom Layer by Correction of μ seg a

Table 1
Optical properties and layer thickness of four layered models used in simulation.

Table 2
Optical properties at 760 nm and layer thickness of epoxy resin-based phantoms, where μ s 0 and F a are theoretical values based on the concentration of titanium oxide and ink.