20 September 2013 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
Author Affiliations +
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 (μ seg a) . 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 μ seg a from the actual μa is expressed by a function of the ratio of μ seg a 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.



In 1977, Jöbsis1 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,3 Over the past 35 years, the NIRS community has focused on solving these issues and has developed a wide range of types of NIRS measurements. Among these, time-domain measurement [time-resolved spectroscopy (TRS)] is a very promising approach to these issues, as TRS provides a temporal profile of the detected light intensity [temporal point spread function (TPSF)], which carries information about depth-dependent attenuation based on the correlation of detection time and penetration depth of photons (time-domain depth sensitivity).45.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 considered 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.8 In conventional time-domain measurements of the human head, the absorption coefficient (μa) of the cerebral tissue has commonly been determined by fitting the TPSF derived from an analytical solution of the photon diffusion equation (DE) to the profile of the time-resolved reflectance on the assumption that the head is semi-infinite and homogeneous.910.11 This approach selectively measures cerebral optical properties, but the selective and quantitative accuracies are inadequate because of limitations of the diffusion approximation to the RTE for inhomogeneous biological tissue.1213.14 Other methods have also been proposed to estimate the optical properties of layered mediums, such as the multilayered (time-dependent) DE,1516.17 also a method with spatially and time-resolved reflectance,18 and a method considering time-dependent mean partial path lengths.19 However, 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 detectors20 and is thought to be the most promising technique for the quantitative detection of focal changes in cerebral hemodynamics. DOT can be performed with time-domain,21 frequency-domain,22 and continuous-wave (CW) measurement instruments.23 Recently, CW high-density DOT with high spatial resolution has been developed and applied to functional hemodynamic maps of the adult human visual cortex.24 However, 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.25

We 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.26 In 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 two-layered 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 method26 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) 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.”


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 IR(t)=SR(t)exp(μa,R·cR·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 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 [AR(t)], Adiff(t), can be expressed by


where Sdiff(t)=ln[SR(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, Adiff(t) is not a linear function of t even though Sdiff(t) can be ignored and c=cR, that is, μa is different at each detection time. The method introduces a time-dependent apparent absorption coefficient, μa(t), given by iμaili(t)/L(t), where the subscript i refers to the i’th layer in the object medium, li(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 dAdiff(t)/dt is expressed as follows:


where dAdiff(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 Adifft 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, dAdiff(t)/dt, is introduced, the value of this is obtained by ΔAdiff(t)/Δt, rather than dAdiff(t)/dt. In this study, the Adifft 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 (dAdiff(t)/dt) by linear regression. The mean μa(t) of each time segment is referred to as the time-segmented μa, μaseg. Then, assuming that a reference and an object have the same μs and refractive index, this gives the following equation,



The time-segmented μa(μaseg) was calculated by substituting dAdiff(t)/dt, μa,R, and c in Eq. (3).

The segments of the Adifft 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, μaseg (K) denotes the μaseg in segment K.


Simulation and Experiment


Monte Carlo Simulations – Calculation of the Time-Resolved Reflectances

To calculate the reflectances from four-layered (object) and from homogeneous (reference) mediums, we used the Monte Carlo code developed by Wang and Jacques.29 The 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=10mm1, anisotropy parameter g=0.9), were compared with those in an isotropic setting (μs=1.0mm1) 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.219mm/ps (refractive index 1.37) was used in all the calculations.

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 mode3031.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.33 In 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,37 In 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±0.05 (mean±SD). Thus, it was concluded that the statistical error of the Monte Carlo simulation was negligible.

Table 1

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

Layer (#)μs′ (mm−1)μa (mm−1)d (mm)
Group A
11.30.0013–0.02(d1,d2)=(3,5), (3, 7),
21.00.011–0.014(4, 7), (5, 7), (5, 5), (5, 9)
Group B
11.10.0013–0.02(d1,d2)=(3,7), (5, 7)
Group C


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, μa, and the thickness d of each layer are denoted by μsi, μai, and di, 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 The depth from the scalp to the cortical surface was set to 9.0 to 15.0 mm in the forehead region.4041.42

The models were classified into three groups according to the μs1 and μs2 values: 1.3 and 1.0mm1 for Group A, 1.1 and 0.8mm1 for Group B, and 1.8 and 1.2mm1 for Group C. The μs3 and μs4 values were set to 0.1 and 1.7mm1 in all three groups. In Group A, μa was varied in the range of 0.013 to 0.02mm1 in layer #1, 0.011 to 0.014mm1 in layer #2, and 0.01 to 0.02mm1 in layer #4 separately. In Groups B and C, only μa1 and μa4 were varied and μa2 was set to 0.012mm1. 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,4344.45 and μa3 was 0.0033mm1. The thicknesses of layers #1 and #2, (d1,d2) (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 d3 and d4 were set as 1 and 90 mm, respectively. We mainly evaluated temporal variations in μaseg and differences between μaseg 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 variation in the upper layers on μaseg. The μs of the homogeneous reference was changed within the range of 1.1 to 1.8mm1 as described below, while the μa and the thickness were set to 0mm1 and 100 mm, respectively.


Determination of μs for the Reference Mediums

In the previous study,26 it was found that the results of our method depend on the μs 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 value of a multilayered medium determined by the curve fitting procedure as described in Sec. 3.1 is closer to the μs 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 value but independent of the μa value, when the optical properties are in the range of the published data for the biological tissue.46 In this article, the μs 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), μaseg depends on the temporal change in li(t) in the time range of each segment in each layer. Therefore, to interpret the calculated μaseg, 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, (d1,d2)=(5,5); and μa1=0.015, μa2=0.012, μa4=0.01, (d1,d2)=(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+1mm (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 Fs and μa to 0.1 to 1.77 and 0.0033 to 0.015mm1 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 values of the objects determined in the two time ranges used in the fitting (1.520mm1, 250 to 1000 ps; and 1.502mm1, 0 to 4500 ps). The μa and μs of the reference were 0.001 and 1.52mm1, respectively (Table 2).

Table 2

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

Layer (#)μs′ (mm−1)μa (mm−1)d (mm)



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 time-resolved reflectance as described in Sec. 3.1. The light velocity in the phantoms was assumed to be 0.19mm/ps, corresponding to a refractive index of 1.58. Before the calculation of Adiff(t), the profiles of the time-resolved reflectance measured by TRS were deconvoluted by the instrument function by using Bayesian deconvolution.47




Temporal Variations in μaseg

Figure 1 shows the temporal variations in μaseg 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 FFFFa4 and μa1, while μa1 were larger than (diamonds), equal to (crosses), or smaller than μa4 (circles). The μaseg 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 μaseg to the μa4 [μaseg(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)], μaseg(k)/μa4 is nearest 1 in the time segments later than IV. For all the models in Fig. 1(a), the upper layer thickness (d1,d2) influenced the values of μaseg in segments-I and -II but not those in the later segments.

Fig. 1

Temporal variations of μaseg of four-layered mediums (Group A) with two different sets of (d1,d2) predicted by Monte Carlo simulation. (d1,d2)=(3,7) (black symbols) and (d1,d2)=(5,9) (gray symbols). (a) μaseg in the models where both μa1 and μa2 are larger than μa4 (black filled square, gray filled square; μa1=0.015mm1, μa2=0.012mm1, μa4=0.01mm1), and μa1 is larger than and μa2 is equal to μa4 (black filled triangle, gray filled triangle; μa1=0.015mm1, μa2=0.012mm1, μa4=0.012mm1). (b) μaseg in the models where μa2 is smaller than μa4, while μa1 is larger than (black filled diamond, gray filled diamond; μa1=0.018mm1, μa2=0.012mm1, μa4=0.015mm1), equal to (μa1=0.015mm1, μa2=0.012mm1, μa4=0.015mm1), and smaller than (black filled circle, gray filled circle; μa1=0.015mm1, μa2=0.012mm1, μa4=0.02mm1) μa4. The vertical axis represents the ratio of μaseg to μa4 (μaseg(k)/μa4).


In the models where μa2 was smaller than μa4 [Fig. 1(b)], the μaseg(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 μaseg was almost constant in segments-II to -VI, and a little smaller than the μa4 values. The values of μaseg in the later time segments were independent of the value of μa1. When both μa1 and μa2 were smaller than μa4, the μaseg(k)/μa4 values were much smaller than one in all the time segments. The μaseg decreased between segments I and II, then μaseg 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 μaseg in the later time segments were larger in the medium with (d1,d2)=(5,9) than in that with (3, 7). Further, the μaseg 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 μaseg was observed in other models with different sets of μai and thicknesses.


Differences Between the μaseg and the Actual μa in the Bottom Layer

In the previous study of two-layered models, it was found that the ratio of μaseg (V) to the actual μa of the bottom layer was expressed by a function of the ratio of μaseg in an earlier time segment to that in a later one. In this article we examined the relationship between the ratio of μaseg(V) to μa4[μaseg(V)/μa4] and the ratio of μaseg in segment-I to that in segment-V (μasegratioI/V) obtained from Group A (Fig. 2). All data obtained from the various mediums with different sets of (d1,d2) were plotted together in the graph. The plot shows that the μaseg(V)/μa4 could be expressed as one polynomial function of the μasegratioI/V [μaseg(V)/μa4=1.70(μasegratioI/V)2+4.82(μasegratioI/V)2.35] with an adjusted coefficient of determination R2=0.976. The deviation of μaseg(V)/μa4 from the regression curve was less than 6% (mean error was 1.6%).

Fig. 2

The ratio of the μaseg (V) to μa4 in reference to the ratio of the μaseg (I) to (V) (μasegratioI/V) predicted by Monte Carlo simulation in all the mediums of Group A. Different symbols represent six different sets of (d1,d2): (3, 5), (+); (3, 7), (∇); (4, 7), (×); (5, 7), (triangle); (5, 5), (filled circle); (5, 9), (diamond). Solid line represents the polynomial regression curve obtained from all the data.


Figure 3 is a plot of the μaseg(V)/μa4 versus the μasegratioI/V for Groups B and C, where the μs values of layers #1 and #2 are different from those in Group A. The solid line represents the regression curve of the relation between μaseg(V)/μa4 and μasegratioI/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 R2=0.975). The deviation of μaseg(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 μaseg(V)/μa4 and the μasegratioI/V can be expressed by [1.70(μasegratioI/V)2+4.82(μasegratioI/V)2.35] irrespective of the μs values of the upper layers.

Fig. 3

The ratio of the μaseg (V) to μa4 in reference to the ratio of the μaseg (I) to (V) (μasegratioI/V) predicted by Monte Carlo simulation in all the mediums in Groups B and C, where μs 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.



Temporal Variations in the s-mPPL in Each Layer

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 [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.

Fig. 4

Temporal variations of s-mPPL predicted by the Monte Carlo simulation at each layer of #1−#4 in two four-layered mediums in Group A. (a) μa1=0.015mm1, μa2=0.012mm1, μa4=0.02mm1, (d1,d2)=(5,5); (b) μa1=0.015mm1, μa2=0.012mm1, μa4=0.01mm1, (d1,d2)=(5,5)

The s-mPPL in each layer was calculated by (the sum of the mPPL/mTPL in each depth within the corresponding layer) × (the midpoint time in each time segment) × (light velocity in the mediums).


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 [IR(t)] at 760 nm, which were deconvoluted by the incident light, and the Adifft curve, respectively. The profiles of Adiff(t) were reliable in the time range of 0.5 to 3 ns (segments-I to -V). From Fig. 5(b), μaseg(I) and μaseg(V) were estimated at each measurement position. Substituting these values into the function [μaseg(V)/μa4=1.70(μasegratioI/V)2+4.82(μasegratioI/V)2.35], μa4 was estimated to be 0.0156±0.002/mm (mean±SD). This value was very close to the actual μa4 of the object (0.015/mm).

Fig. 5

Time-resolved reflectance and Adifft curve obtained by phantom measurements. (a) Temporal profiles of detected light through a reference [IR(t)] and an object medium [I(t)] on a natural logarithmic vertical scale. The time of the peak of the instrument function was taken as 0 s. The data of detected light intensity are deconvoluted by the temporal profile of the incident light. (b) Adifft curve obtained from the data in (a) in the time period of 0 to 3.5 ns. The Roman numbers show segment numbers of the curve parted every 0.5 ns to obtain the μaseg.





Characterization of μaseg in View of the Change in s-mPPL in Each Layer

The temporal variations in μaseg in four-layered mediums (Fig. 1) largely agreed with the previous results for a two-layered model.26 The μaseg changed with the time segment, reflecting the differences in μa in the depth direction. The main finding is that the μaseg values in later time segments exclusively depend on the differences between μa2 and μa4. That is, the μaseg(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 μaseg(V) was smaller than μa4 even in the later time segments irrespective of μa1 values [Fig. 1(b)]. In addition, the dependence of μaseg on the upper-layer thickness varied with the difference between μa2 and μa4 (black versus grey symbols in Fig. 1).

As defined in Sec. 2.2, μaseg represents the mean μa(t) in a segment. The μaseg value depends on the mean temporal change in li(t) for each time segment in the individual layers [Eq. (2)]. Therefore, a consideration of the contribution of the mean change in li(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 li(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.

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 l2(t) to the change in L(t) was high enough to strongly influence the μaseg even in later time segments. This dominant change in l2(t) can explain the finding that the μaseg at the later time represented values smaller than μa4 in the mediums where μa2 is smaller than μa4. The dominance of the change in l2(t) also accounts for the difference in μaseg at later time segments between the mediums with (d1,d2)=(3,7) and with (5, 9). The reason why there was little effect of μa1 on the temporal pattern of changes in μaseg may be ascribed to the smaller temporal change in s-mPPL in layer #1.

When μa2 was larger than μa4, μaseg 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 l2(t) to the change in L(t) in these segments was smaller than that of l4(t). Therefore, it was concluded that the size difference between the μa2 and the μa4 is critical for the change in li(t) in a later time segment, which determines μaseg values.


Determination of μa in the Bottom Layer by Correction of μaseg

Like in two-layered models, when the μa value of the bottom layer was larger than that of the upper layer, the μaseg(V) was smaller than the μa4. With two-layered models, the ratio of μaseg(IV) or (V) to the real μa of the lower layer is expressed by the affine function of the μasegratioI/V, which enables estimation of the μa of the lower layer.26 The difference between μaseg(V) and μa4 in the four-layered models also depends on the μasegratioI/V, though the μaseg(V)/μa4 was not expressed by an affine function but by a polynomial function of the μasegratioI/V (Figs. 2 and 3). This relationship between the μaseg(V)/μa4 and μasegratioI/V was independent of the μs values and thicknesses of the layers #1 and #2. In addition to the high R2 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 μaseg(V)/μa4 and μasegratioI/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 between the upper layer of the object and the reference significantly influenced the μaseg in the time segments earlier than 2500 ps.26 However, this difference did not influence the μaseg after 2500 ps, and the difference in μs between the lower layer of the object and the reference rarely affected the μaseg values in all the time segments. The present study used a curve fitting procedure to determine the μs value of the reference. The effect of the time range of the curve fitting on the μaseg within the range of 250 to 3000 ps was investigated. When the μs of the reference was determined by fitting over the wider time ranges, the plots of μaseg(V)/μa4 versus μasegratioI/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 of the upper layers. However, there were no differences in μs 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, 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±0.002/mm) and actual μa values (0.015/mm) was maximally about 0.008/mm (about 5%). Considering the deviation of μaseg(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 μaseg(V)/μa4 and μasegratioI/V was independent of the μs values and thicknesses of the layers #1 and #2 in the range used in this study (μs 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,42 Therefore, 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.



This study demonstrated that the “time segment analysis” of time resolved reflectance could be applied to four-layered slab models. The μaseg(V)/μa4 was expressed by a unique function of the μasegratioI/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.


The authors very much appreciate the helpful discussion with Y. Yamashita and M. Oda in Hamamatsu Photonics KK. This research was supported by Grants-in Aid for Scientific Research (C) from the Japan Society for the Promotion of Science (No. 16500468).


1. F. F. Jöbsis, “Noninvasive infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science 198(4323), 1264–1267 (1977).SCIEAS0036-8075 http://dx.doi.org/10.1126/science.929199 Google Scholar

2. T. Takahashiet al., “Influence of skin blood flow on near-infrared spectroscopy signals measured on the forehead during a verbal fluency task,” NeuroImage 57(3), 991–1002 (2011).NEIMEF1053-8119 http://dx.doi.org/10.1016/j.neuroimage.2011.05.012 Google Scholar

3. E. Kirilinaet al., “The physiological origin of task-evoked systemic artifacts in functional near infrared spectroscopy,” NeuroImage 61(1), 70–81 (2012).NEIMEF1053-8119 http://dx.doi.org/10.1016/j.neuroimage.2012.02.074 Google Scholar

4. S. Del BiancoF. MartelliG. Zaccanti, “Penetration depth of light re-emitted by a diffusive medium: theoretical and experimental investigation,” Phys. Med. Biol. 47(23), 4131–4144 (2002).PHMBA70031-9155 http://dx.doi.org/10.1088/0031-9155/47/23/301 Google Scholar

5. J. SelbD. K. JosephD. A. Boas, “Time-gated optical system for depth-resolved functional brain imaging,” J. Biomed. Opt. 11(4), 044008 (2006).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.2337320 Google Scholar

6. A. Torricelliet al., “Time-resolved reflectance at null source-detector separation: improving contrast and resolution in diffuse optical imaging,” Phys. Rev. Lett. 95(7), 078101 (2005).PRLTAO0031-9007 http://dx.doi.org/10.1103/PhysRevLett.95.078101 Google Scholar

7. G. Zaccantiet al., “Optical properties of biological tissues,” Proc. SPIE 2389, 513–521 (1995).PSISDG0277-786X http://dx.doi.org/10.1117/12.210000 Google Scholar

8. A. Ishimaru, “Diffusion of a pulse in densely distributed scatters,” J. Opt. Soc. Am. 68(8), 1045–1050 (1978).JOSAAH0030-3941 http://dx.doi.org/10.1364/JOSA.68.001045 Google Scholar

9. S. Ijichiet al., “Quantification of cerebral hemoglobin as a function of oxygenation using near-infrared time-resolved spectroscopy in a piglet model of hypoxia,” J. Biomed. Opt. 10(2), 024026 (2005).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.1899184 Google Scholar

10. V. Quaresimaet al., “Bilateral prefrontal cortex oxygenation responses to a verbal fluency task: a multichannel time-resolved near-infrared topography study,” J. Biomed. Opt. 10(1), 011012 (2005).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.1851512 Google Scholar

11. Y. Hoshiet al., “Resting hypofrontality in schizophrenia: a study using near-infrared time-resolved spectroscopy,” Schizophr. Res. 84(2–3), 411–420 (2006).SCRSEH0920-9964 http://dx.doi.org/10.1016/j.schres.2006.03.010 Google Scholar

12. A. LiebertA. Kienle, “Light diffusion in N-layered turbid media: frequency and time domains,” J. Biomed. Opt. 15(2), 025002 (2010).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.3368682 Google Scholar

13. A. Kienleet al., “Investigation of two-layered turbid media with time-resolved reflectance,” Appl. Opt. 37(28), 6852–6862 (1998).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.37.006852 Google Scholar

14. A. H. Hielscheret al., “Time-resolved photon emission from layered turbid media,” Appl. Opt. 35(4), 719–728 (1996).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.35.000719 Google Scholar

15. A. Kienleet al., “Noninvasive determination of the optical properties of two-layered turbid media,” Appl. Opt. 37(5), 779–791 (1998).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.37.000779 Google Scholar

16. F. Martelliet al., “Analytical approximate solutions of the time-domain diffusion equation in layered slabs,” J. Opt. Soc. Am. A 19(1), 71–80 (2002).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.19.000071 Google Scholar

17. F. Martelliet al., “Phantom validation and in vivo application of an inversion procedure for retrieving the optical properties of diffusive layered media from time-resolved reflectance measurements,” Opt. Lett. 29(17), 2037–2039 (2004).OPLEDP0146-9592 http://dx.doi.org/10.1364/OL.29.002037 Google Scholar

18. M. ShimadaY. HoshiY. Yamada, “Simple algorithm for the measurement of absorption coefficients of a two-layered medium by spatially resolved and time-resolved reflectance,” Appl. Opt. 44(35), 7554–7563 (2005).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.44.007554 Google Scholar

19. J. Steinbrinket al., “Determining changes in NIR absorption using a layered model of the human head,” Phys. Med. Biol. 46(3), 879–896 (2001).PHMBA70031-9155 http://dx.doi.org/10.1088/0031-9155/46/3/320 Google Scholar

20. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15(6), R41–R93 (1999).INPEEY0266-5611 http://dx.doi.org/10.1088/0266-5611/15/2/022 Google Scholar

21. F. GaoP. PouletY. Yamada, “Simultaneous mapping of absorption and scattering coefficients from a three-dimensional model of time-resolved optical tomography,” Appl. Opt. 39(31), 5898–5910 (2000).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.39.005898 Google Scholar

22. B. W. Pogueet al., “Initial assessment of a simple system for frequency domain diffuse optical tomography,” Phys. Med. Biol. 40(10), 1709–1729 (1995).PHMBA70031-9155 http://dx.doi.org/10.1088/0031-9155/40/10/011 Google Scholar

23. D. A. Boaset al., “The accuracy of near infrared spectroscopy and imaging during focal changes in cerebral hemodynamics,” NeuroImage 13(1), 76–90 (2001).NEIMEF1053-8119 http://dx.doi.org/10.1006/nimg.2000.0674 Google Scholar

24. B. W. Zeffet al., “Retinotopic mapping of adult human visual cortex with high-density diffuse optical tomography,” Proc. Natl. Acad. Aci. U. S. A. 104(29), 12169–12174 (2007). Google Scholar

25. A. CharetteJ. BoulangerH. K. Kim, “An overview on recent radiation transport algorithm development for optical tomography imaging,” JQSRT 109(17–18), 2743–2766 (2008).JQSRAE0022-4073 http://dx.doi.org/10.1016/j.jqsrt.2008.06.007 Google Scholar

26. C. Satoet al., “Extraction of depth-dependent signals from time-resolved reflectance in layered turbid media,” J. Biomed. Opt. 10(6), 064008 (2005).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.2136312 Google Scholar

27. Y. NomuraO. HazekiM. Tamura, “Exponential attenuation of light along nonlinear path through the biological model,” Adv. Exp. Med. Biol. 248, 77–80 (1989).AEMBAP0065-2598 http://dx.doi.org/10.1007/978-1-4684-5643-1 Google Scholar

28. Y. TsuchiyaT. Urakami, “Photon migration model for turbid biological medium having various shapes,” Jpn. J. Appl. Phys. 34, L79–L81 (1995).JJAPA50021-4922 http://dx.doi.org/10.1143/JJAP.34.L79 Google Scholar

29. L. WangS. L. Jacques, Monte Carlo modeling of light transport in multi-layered tissues in standard C, University of Texas/M.D. Anderson Cancer Center, Houston (1992). Google Scholar

30. T. J. FarrelM. S. Patterson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo.” Med. Phys. 19(4), 879–888 (1992).MPHYA60094-2405 http://dx.doi.org/10.1118/1.596777 Google Scholar

31. R. C. Haskellet al., “ Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A Opt. Image Sci. Vis. 11(10), 2727–2741 (1994).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.11.002727 Google Scholar

32. A. KienleM. S. Patterson, “Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium,” J. Opt. Soc. Am. A Opt. Image Sci. Vis. 14(1), 246–254 (1997).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.14.000246 Google Scholar

33. K. FurutsuY. Yamada, “Diffusion approximation for a dissipative random medium and the applications,” Phys. Rev. E 50(1), 3634–3640 (1994).PLEEE81063-651X http://dx.doi.org/10.1103/PhysRevE.50.3634 Google Scholar

34. K. Levenberg, “A method for the solution of certain problems in least squares, Quart. Appl. Math. 2, 164–168 (1944).QAMAAY0033-569X Google Scholar

35. D. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” SIAM J. Appl. Math. 11(2), 431–441 (1963).SMJMAP0036-1399 http://dx.doi.org/10.1137/0111030 Google Scholar

36. P. B. Bevington, Data Reduction and Error Analysis for the Physical Science, McGraw-Hill, New York (1969). Google Scholar

37. V. D. O ConnorD. Phillips, Time-Correlated Single-Photon Counting, Academic Press, London (1985). Google Scholar

38. F. Bevilacquaet al., “In vivo local determination of tissue optical properties: applications to human brain,” Appl. Opt. 38(22), 4939–4950 (1999).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.38.004939 Google Scholar

39. J. H. Choiet al., “Noninvasive determination of the optical properties of adult brain: near-infrared spectroscopy approach,” J. Biomed. Opt. 9(1), 221–229 (2004).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.1628242 Google Scholar

40. M. Okamotoet al., “Three-dimensional probabilistic anatomical cranio-cerebral correlation via the international 10–20 system oriented for transcranial functional brain mapping,” Neuroimage 21(1), 99–111 (2004).NEIMEF1053-8119 http://dx.doi.org/10.1016/j.neuroimage.2003.08.026 Google Scholar

41. H. Zhaoet al., “Maps of optical differential pathlength factor of human adult forehead, somatosensory motor, occipital regions at multi-wavelengths in NIR,” Phys. Med. Biol. 47(1), 2075–2093 (2002).PHMBA70031-9155 http://dx.doi.org/10.1088/0031-9155/47/12/306 Google Scholar

42. A. Liebertet al., “Time-resolved multidistance near-infrared spectroscopy of the adult head: intracerebral and extracerebral absorption changes from moments of distribution of times of flight of photons,” Appl. Opt. 43(15), 3037–3047 (2004).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.43.003037 Google Scholar

43. D. A. BoasA. M. Dale, “Simulation study of magnetic resonance imaging-guided cortically constrained diffuse optical tomography of human brain function,” Appl. Opt. 44(9), 1957–1968 (2005).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.44.001957 Google Scholar

44. Y. Hoshiet al., “Reevaluation of near-infrared light propagation in the adult human head: implications for functional near-infrared spectroscopy,” J. Biomed. Opt. 10(6), 064032 (2005).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.2142325 Google Scholar

45. E. OkadaD. T. Delpy, “Near-infrared light propagation in an adult head model. I. Modeling of low-level scattering in the cerebrospinal fluid layer,” Appl. Opt. 42(16), 2906–2914 (2003).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.42.002906 Google Scholar

46. J. Johanssonet al., “Time-resolved NIR/Vis spectroscopy for analysis of solids: Pharmaceutical tablets.” Appl. Spectro. 56(6), 725–731 (2002).APSPA40003-7028 http://dx.doi.org/10.1366/000370202760077676 Google Scholar

47. T. J. KennettW. V. PrestwichA. Robertson, “Baysian deconvolution II: noise properties.” Nucl. Instr. Meth. 151, 293–301 (1978).NUIMAL0029-554X http://dx.doi.org/10.1016/0029-554X(78)90503-7 Google Scholar

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Chie Sato, Chie Sato, Miho Shimada, Miho Shimada, Yukari Tanikawa, Yukari Tanikawa, Yoko Hoshi, Yoko Hoshi, } "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," Journal of Biomedical Optics 18(9), 097005 (20 September 2013). https://doi.org/10.1117/1.JBO.18.9.097005 . Submission:

Back to Top