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 () 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 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 (time-segmented ) was estimated from the slope of each segment. The difference between the estimated and the actual in the lower layer then depends on the ratio of the actual in the upper layer to that in the lower layer, which was linearly related to the ratio of the estimated in an early time segment to that in a later segment. Using this relationship, it was possible to correct the obtained from the later time segment to derive the 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 under several conditions in which and the reduced scattering coefficient () of each layer were changed, (2) the difference between the time-segmented at a later time and the actual of the bottom layer, and its relation to the ratio of the time-segmented 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 values. Then, the experimental results are reported: (4) the determination of of the bottom layer of a four-layered epoxy resin-based phantom with the “time segment analysis.”
Methods for Determining the of the Bottom Layer
The details of the method have been described in a previous paper,26 and the following is a brief summary.
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 for the reference and the object is expressed by and ,27,28 where is the light velocity in a medium and is a scattering function which is dependent on the of the medium in the photon diffusion regime. The light attenuation, , 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 and of a reference , , can be expressed by
In the case of an object that is a layered medium in which varies with layer (i.e., depth) and the reference is homogeneous, is not a linear function of even though can be ignored and , that is, is different at each detection time. The method introduces a time-dependent apparent absorption coefficient, , given by , where the subscript refers to the ’th layer in the object medium, is the stochastic time-dependent partial path length in the ’th layer for the photons detected at time , and is the total path length of photons detected at time , which is defined as . Then is expressed as follows:
Time Segment Analysis
To simplify the process of the analysis, the curve within a short time range (e.g., a few hundred picoseconds, ) is considered a straight line, and a term for the mean absorption differences, , is introduced, the value of this is obtained by , rather than . In this study, the 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 () by linear regression. The mean of each time segment is referred to as the time-segmented , . Then, assuming that a reference and an object have the same and refractive index, this gives the following equation,
The time-segmented was calculated by substituting , , and in Eq. (3).
The segments of the curve were numbered as 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, () denotes the in segment .
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 , anisotropy parameter ), were compared with those in an isotropic setting () 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 (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 .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 () 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 value mainly represented statistical bias and statistical variations that were attributed to the Monte Carlo simulation. The values in all the mathematical models described in Sec. 3.2 (Table 1) were (). Thus, it was concluded that the statistical error of the Monte Carlo simulation was negligible.
Optical properties and layer thickness of four layered models used in simulation.
|Layer (#)||μs′ (mm−1)||μa (mm−1)||d (mm)|
|1||1.3||0.0013–0.02||(d1,d2)=(3,5), (3, 7),|
|2||1.0||0.011–0.014||(4, 7), (5, 7), (5, 5), (5, 9)|
|1||1.1||0.0013–0.02||(d1,d2)=(3,7), (5, 7)|
Four-Layered Models for Simulation
We used a four-layered semi-infinite model as a human head model. The layers are numbered , and 4 in order of depth from the outside, and denoted by #: the layers #1 to #4 correspond to the scalp, skull, cerebrospinal fluid (CSF), and brain, respectively. The , , and the thickness of each layer are denoted by , , and , 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 and values: 1.3 and for Group A, 1.1 and for Group B, and 1.8 and for Group C. The and values were set to 0.1 and in all three groups. In Group A, was varied in the range of 0.013 to in layer #1, 0.011 to in layer #2, and 0.01 to in layer #4 separately. In Groups B and C, only and were varied and was set to . In all the models, was set to be larger than as the of the scalp is generally larger than that of the skull in adult heads,4344.–45 and was . The thicknesses of layers #1 and #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 and were set as 1 and 90 mm, respectively. We mainly evaluated temporal variations in and differences between and the actual for the Group A models. The models in Groups B and C were used to investigate the effects of the variation in the upper layers on . The of the homogeneous reference was changed within the range of 1.1 to as described below, while the and the thickness were set to and 100 mm, respectively.
Determination of for the Reference Mediums
In the previous study,26 it was found that the results of our method depend on the value of the reference, which was ideally the same as that of the upper layer of a two-layered object medium. In general, the value of a multilayered medium determined by the curve fitting procedure as described in Sec. 3.1 is closer to the 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 value but independent of the value, when the optical properties are in the range of the published data for the biological tissue.46 In this article, the 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), depends on the temporal change in in the time range of each segment in each layer. Therefore, to interpret the calculated , 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 [, , , ; and , , , ], were perpendicularly partitioned into 1 mm thick laminae and the attenuation changes () caused by altering in a lamina by of the initial value were calculated. The time-independent mean partial path length of photons in the lamina where is changed can be estimated by the 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 to () obtained by dividing the calculated in each time segment by . 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, (50 for the reference) mm, were prepared. The base of the phantoms was epoxy resin (Beuhler, Lake Bluff, Illinois). We adjusted the and to 0.1 to 1.77 and 0.0033 to 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 values of the objects determined in the two time ranges used in the fitting (, 250 to 1000 ps; and , 0 to 4500 ps). The and of the reference were 0.001 and , respectively (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 ]. Reflected light is collected by a fiber bundle (3-mm diameter, ). 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 , corresponding to a refractive index of 1.58. Before the calculation of , the profiles of the time-resolved reflectance measured by TRS were deconvoluted by the instrument function by using Bayesian deconvolution.47
Temporal Variations in
Figure 1 shows the temporal variations in obtained from the nine models in Group A. Figure 1(a) shows results in the mediums where both and were larger than (squares) and where was larger but was equal to (triangles). Figure 1(b) shows results in the mediums where was smaller than and , while were larger than (diamonds), equal to (crosses), or smaller than (circles). The values are expressed as the ratio to . When both and were larger than [squares in Fig. 1(a)], the ratios of the to the 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 . When was equal to [triangles in Fig. 1(a)], is nearest 1 in the time segments later than IV. For all the models in Fig. 1(a), the upper layer thickness () influenced the values of in segments-I and -II but not those in the later segments.
In the models where was smaller than [Fig. 1(b)], the in all the time segments later than segment-II were smaller than 1. When was equal to or larger than , the was almost constant in segments-II to -VI, and a little smaller than the values. The values of in the later time segments were independent of the value of . When both and were smaller than , the values were much smaller than one in all the time segments. The decreased between segments I and II, then increased slightly in the later time segments, and remained about 20% smaller than . For all the models in Fig. 1(b), the differences between and in the later time segments were larger in the medium with than in that with (3, 7). Further, the 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 was observed in other models with different sets of and thicknesses.
Differences Between the and the Actual in the Bottom Layer
In the previous study of two-layered models, it was found that the ratio of (V) to the actual of the bottom layer was expressed by a function of the ratio of in an earlier time segment to that in a later one. In this article we examined the relationship between the ratio of to and the ratio of in segment-I to that in segment-V () obtained from Group A (Fig. 2). All data obtained from the various mediums with different sets of () were plotted together in the graph. The plot shows that the could be expressed as one polynomial function of the with an adjusted coefficient of determination . The deviation of from the regression curve was less than 6% (mean error was 1.6%).
Figure 3 is a plot of the versus the for Groups B and C, where the values of layers #1 and #2 are different from those in Group A. The solid line represents the regression curve of the relation between and 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 ). The deviation of 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 and the can be expressed by irrespective of the values of the upper layers.
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 was larger than [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 was smaller than [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.
Figures 5(a) and 5(b) show the profiles of the time-resolved reflectances of the object  and reference  at 760 nm, which were deconvoluted by the incident light, and the curve, respectively. The profiles of were reliable in the time range of 0.5 to 3 ns (segments-I to -V). From Fig. 5(b), and were estimated at each measurement position. Substituting these values into the function , was estimated to be (). This value was very close to the actual of the object ().
Characterization of in View of the Change in s-mPPL in Each Layer
The temporal variations in in four-layered mediums (Fig. 1) largely agreed with the previous results for a two-layered model.26 The changed with the time segment, reflecting the differences in in the depth direction. The main finding is that the values in later time segments exclusively depend on the differences between and . That is, the accurately represents the in the mediums where was larger than or equal to [Fig. 1(a)], whereas in mediums where was smaller than , the was smaller than even in the later time segments irrespective of values [Fig. 1(b)]. In addition, the dependence of on the upper-layer thickness varied with the difference between and (black versus grey symbols in Fig. 1).
As defined in Sec. 2.2, represents the mean in a segment. The value depends on the mean temporal change in for each time segment in the individual layers [Eq. (2)]. Therefore, a consideration of the contribution of the mean change in to the change in (; 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 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 was smaller than , 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 to the change in was high enough to strongly influence the even in later time segments. This dominant change in can explain the finding that the at the later time represented values smaller than in the mediums where is smaller than . The dominance of the change in also accounts for the difference in at later time segments between the mediums with and with (5, 9). The reason why there was little effect of on the temporal pattern of changes in may be ascribed to the smaller temporal change in s-mPPL in layer #1.
When was larger than , at segments-V and -VI mainly represented the [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 to the change in in these segments was smaller than that of . Therefore, it was concluded that the size difference between the and the is critical for the change in in a later time segment, which determines values.
Determination of in the Bottom Layer by Correction of
Like in two-layered models, when the value of the bottom layer was larger than that of the upper layer, the was smaller than the . With two-layered models, the ratio of (IV) or (V) to the real of the lower layer is expressed by the affine function of the , which enables estimation of the of the lower layer.26 The difference between and in the four-layered models also depends on the , though the was not expressed by an affine function but by a polynomial function of the (Figs. 2 and 3). This relationship between the and was independent of the values and thicknesses of the layers #1 and #2. In addition to the high 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 and 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 between the upper layer of the object and the reference significantly influenced the in the time segments earlier than 2500 ps.26 However, this difference did not influence the after 2500 ps, and the difference in between the lower layer of the object and the reference rarely affected the values in all the time segments. The present study used a curve fitting procedure to determine the value of the reference. The effect of the time range of the curve fitting on the within the range of 250 to 3000 ps was investigated. When the of the reference was determined by fitting over the wider time ranges, the plots of versus 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 of the upper layers. However, there were no differences in 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 , 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 () and actual values () was maximally about (about 5%). Considering the deviation of from the regression curve (Figs. 2 and 3) that was less than 8%, the estimation error is acceptable. Further, since the relationship between the and was independent of the values and thicknesses of the layers #1 and #2 in the range used in this study ( 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 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 was expressed by a unique function of the on the condition that the thickness of the upper layer was less than 15 mm. Since the 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).