1.

Introduction

Photon migration through a turbid medium, e.g. tissue, depends on the optical properties of that medium. These properties are represented by the absorption (μ_a) and scattering (μ_s) coefficients, the probability distribution of the scattering angles, p(θ), and the tissue refractive index. Recently, large efforts have been spent on in vivo determination of the tissue optical properties, using affordable techniques, such as measurement of the spatially resolved diffuse reflectance, that can be implemented in instruments for clinical use.^{1 2 3 4 5 6 7 8 9} The optical properties can be used directly for diagnostics, since, e.g., brain and breast tumors have higher absorption than surrounding, non-neoplastic tissue.^{1 10} Furthermore, knowledge of tissue bulk absorption at distinct wavelengths can be used for analysis of the tissue content of major chromophores, such as hemoglobin, melanin, fat, and water.¹⁰

By using two wavelengths, the relative amounts of oxygenated and reduced hemoglobin can be estimated, as in pulse oximetry. However, the latter method is based on the assumption that the migrated photon pathlength is equivalent at the two different wavelengths⁶—an assumption that can be questioned, considering the fact that optical properties display a wavelength dependence.^{9 11} Laser Doppler flowmetry (LDF) is an established method for monitoring microvascular perfusion in vivo. However, one problem with the method is the difficulty in making inter- and even intraindividual comparisons, due to varying tissue optical properties. It has been shown, that for a constant perfusion, the LDF output signal is affected both by changes in scattering and absorption of the turbid medium.¹² The generated perfusion estimate is based on the interaction between photons and moving scatterers, an estimate sensitive to both the velocity and concentration¹³ of moving scatterers (mainly red blood cells) as well as the photon pathlength.^{14 15} The origin of the pathlength variations is not only found in the tissue optical properties, but also in the source detector separation, ρ, of the LDF probe (normally 0.25–1 mm).^{12 15}

The tight relationship between optical properties and photon pathlength suggests that similar approaches could be used in estimating the two; hence, possibly enabling us to eliminate the pathlength influence on the LDF readings. To our knowledge, no author has previously presented a method based on local reflectance at small ρ for determination of the photon pathlength in a turbid medium. Such a method could also be used for further development of pathlength dependent methods, e.g. pulse oximetry and reflectance spectroscopy. Therefore, the aim of the present study is to develop a method for local estimation of photon pathlength, using the spatially resolved diffuse reflectance profile in the 0–2 mm range.

We will show in this paper, that for a wide range of optical properties, applicable to human skin, the average pathlength can vary by almost a factor 6 for a source detector separation of 2 mm. We have devised methods based on a one-layer homogenous tissue model, that can predict the average pathlength of the photons at various source detector separations, with a root-mean-square (rms) error of about 5.

2.

Materials and Methods

2.1.

Simulation Model

For small source detector separations (0–2 mm) and optical properties where ρμ_s ^′<10 [μ_s ^′=μ_s(1−〈cos θ〉), the reduced scattering coefficient, 〈cos θ〉=average cosine of the scattering angles], the diffusion approximation of the transport equation is not generally applicable.^{2 16} Further, the diffusion approximation is only valid when μ_s ^′≫μ_a. ^{16 17} Therefore, with the setup used in the present study, photon pathlengths must be statistically determined by means of Monte Carlo simulations. The Monte Carlo simulation software used in this study was developed by de Mul et al.; (MontCarl 2001, version 20.01 a).¹⁸

For all simulations, a homogenous semi-infinite slab with a thickness of 100 mm with different optical properties according to Tables 1 and 2, was used. A low concentration of homogeneously distributed moving scatterers (corresponding to μ_s=0.1 mm ⁻¹ ) with a constant velocity, v=1.0 mm/s, parallel to the slab surface, was introduced in the model. All scattering events, due to both moving and static scatterers, were modeled with the same phase function. In order to estimate the LDF perfusion, all Doppler shifts due to photon interactions with moving scatterers were stored. The refractive index of the ambient air was set to n=1.0, whereas the refractive index of the slab was set to n=1.44, which is considered a relevant value for human skin.^{19 20} A divergent circular light beam (NA=0.37) with a diameter of 0.2 mm and a rectangular intensity distribution impinged on the slab surface. All photons exiting the upper slab surface at a radial distance (ρ) from the center of the source in the range 0.13⩽ρ⩽2.17 mm were detected. This geometry was subsequently transformed mathematically into mimicking a linear array of ten fibers, with one transmitting, and nine receiving fibers (Figure 1). The fibers were located adjacent to each other, with a center-to-center separation of 230 μm, each fiber having a core diameter of 200 μm, and a surrounding cladding with a thickness of 15 μm. The specific geometrical and optical properties of the simulated fiber optic probe is similar to a probe we have previously used for LDF measurements.¹² In all simulations, 500 000 photons were detected. The number of emitted photons in each simulation ranged from less than 800 000 to more than 16 million for the different setups of optical properties. In order to simulate a light source with a constant intensity, the number of detected photons was normalized by the number of emitted photons.

Figure 1

Probe design. Ten fibers are arranged in a linear array, and each fiber has a core diameter of 0.2 mm (0.23 mm including cladding).

Table 1

Optical properties of the reference space. 〈cos θ〉: mean cosine of the photon scattering angle, θ, calculated from the combination of 96 Henyey–Greenstein phase function, and 4 isotropic scattering.
μa (mm−1)	0.01	0.05	0.10	0.15	0.20	0.25
μs (mm−1)	5	10	20	30	40	50
〈cos θ〉	0.80	0.85	0.90	0.95

Table 2

Optical properties of the validation space. 〈cos θ〉: mean cosine of the photon scattering angle, θ, calculated from the combination of 96 Henyey–Greenstein phase function, and 4 isotropic scattering.
μa (mm−1)	0.03	0.075	0.125	0.175	0.225
μs (mm−1)	7.5	15	25	35	45
〈cos θ〉	0.825	0.875	0.925

The Henyey–Greenstein phase function²¹

Eq. (1)

$${\mathit{p}}_{\text{HG}}(\theta )=\frac{1}{4\pi }\frac{1-{\mathit{g}}_{\text{HG}}^2}{{(1+{\mathit{g}}_{\text{HG}}^2-2_{\mathit{g}\text{HG}}\text{cos}\theta )}^{\mathrm{3/2}}}$$

was used in conjunction with an isotropic phase function to describe the photon scattering, where θ is the deflection angle and g _HG is the anisotropy factor. The proportions used were 96 Henyey–Greenstein and 4 isotropic phase function. Hence, the resulting phase function can be expressed as

Eq. (2)

$$p(\theta )=\frac{\beta }{4\pi }+(1-\beta ){\mathit{p}}_{\text{HG}}(\theta )\text{,}$$

where β=0.04. This combination of a highly forward scattering component, such as the Henyey–Greenstein phase function (g _HG ⩾0.7), and an isotropic component, has been found to adequately describe light scattering in biological tissues.^{1 2}

A large reference space, used to develop pathlength estimation methods, was defined. It consisted of equidistant optical properties, with the intention to encompass the range of values from in vitro and in vivo estimations of human epidermis and dermis at λ=632 nm. ^{11 19 20 22} Therefore, all combinations of parameters in Table 1 were chosen as input parameters in the Monte Carlo model, thus, requiring 6×6×4=144 simulations.

Second, a validation space for evaluation of the accuracy of the pathlength estimation methods was defined. It was setup to maximize the distances to the nearest combination of optical properties in the reference space (Table 2). The validation space consisted of 5×5×3=75 simulations.

3.

Results

The variations in the average photon pathlength as a function of the source detector separation at discrete distances, ρ_i, for the 144 different combinations of optical properties in the reference space is depicted in Figure 2. The mean and SD of the pl_i basically increase linearly with the source detector separation. For the ninth detector fiber (ρ_i=2.07 mm), the average pl_i varies almost sixfold (5.0–28.8 mm). Even for the first detector fiber (ρ_i=0.23 mm), the ratio between the longest and the shortest pl_i is close to 2.6 (range 0.94–2.41 mm).

Figure 2

Variations of pathlength (pl_i) as a function of source detector separation (ρ_i). Mean (x)±one SD (error bars) are given for each fiber. The triangles denote the maximum and minimum values.

The result of predicting pl_i of the validation space using the average pl_i from the reference space is presented in Figures 3 and 4. This method overestimates pl_i by 9–19 (mean) with a SD of 18–32, and a rms of the relative error $\left[\left(\widehat{p{l}_i}-p{l}_i\right)/p{l}_i\right]$of 20–37 (hereafter referred to as rms error). The range of pl ratios $\left(\widehat{p{l}_i}/p{l}_i\right)$ for the ninth fiber was 0.52–1.81.

Figure 3

Prediction of pathlength based on average pathlength. Predicted divided by simulated pathlength $\left(\widehat{p{l}_i}/p{l}_i\right)$ is plotted as a function of source detector separation (ρ_i). Mean (x)±one SD (error bars) are given for each fiber. The triangles denote the maximum and minimum values.

Figure 4

Prediction of pathlength based on all of the various methods presented in this paper, given as the rms of the relative error $[(\widehat{p{l}_i}-p{l}_i)/p{l}_i]$ for each method vs ρ_i. Predictions based on: ○, average pathlength, +, KNN method using m₂ and m₃ values based on all nine fibers, ▵, KNN method using R_i values from fibers 3 and 7, ×, KNN method using m₁, m₂, and m₃ values based on all nine fibers, and □, MPR method using m₁, m₂, and m₃ values based on all nine fibers.

The more sophisticated prediction methods KNN, LIP, and MPR yielded results with comparable accuracy. However, the LIP method was not able to predict all pl_i, due to the fact that some predictor points in the validation space fell outside the convex hull of the triangulated reference space. Therefore, this method was not further used. As for the KNN method, the data set was preprocessed using autoscaling, since this in general improved accuracy. The KNN method was evaluated for K in the range 1–16 and the optimal K value was chosen individually for each fiber (range 3–10), where the optimal K is defined as the K value resulting in the prediction with the smallest rms error.

One previously suggested method,²⁴ is to predict pl_i based on two R_i values. The R_i from fibers 3 and 7 ( ρ_i=0.69 mm and ρ_i=1.61 mm ) gave the most accurate $\widehat{p{l}_i}$ (Figure 4), when applying the KNN method. On average, this method yielded a 1.6–5.6 overestimation of $\widehat{p{l}_i},$with a SD of 3.6–6.9 and a rms error of 3.9–8.7. The range of pl ratios for the worst case, fiber nine, was 0.91–1.36.

Extending the KNN approach to three different R_i yielded slightly more accurate and precise results, compared to the case with two different R_i. The R_i values from fibers 2, 4, and 7 appeared to give the most accurate $\widehat{p{l}_i}.$ The $\widehat{p{l}_i}$ is overestimated by 1.5–4.5, with a SD of 4.0–6.4 and a rms error of 4.2–7.6. The range of pl ratios for fiber 6 (worst case) is 0.83–1.24. For the sake of clarity, these data were excluded from Figure 4.

Estimation of pl_i based on m₂ and m₃ (derived from all nine fibers) using the KNN method is also presented in Figure 4. The $\widehat{p{l}_i}$ is overestimated by 2.4–7.3, with a SD of 10.3–28.1 and a rms error of 10.5–28.8. The worst case range of pl ratios was for fiber 9 (0.41–2.01).

The KNN method using m₁, m₂, and m₃, predicted pl_i to within −0.09–1.1, with a SD of 3.8–6.9 and a rms error of 3.8–7.0 (Figure 4). The worst case range of pl ratios (fiber 9) was between 0.82 and 1.16. Decreasing the number of detectors in the fitting to Eq. (3) results in a decreased accuracy and precision of the $\widehat{p{l}_i}$ using the KNN method. In all cases, the best predictions are obtained using all nine detectors. However, generally, accuracy and precision was maintained down to using only the five detector fibers closest to the source. After removal of more fibers, accuracy unequivocally deteriorated.

Applying the MPR method to the same m₁, m₂, and m₃ data set, slightly improved the predictions (Figures 4 and 5). The $\widehat{p{l}_i}$ mean error ranged between −0.01 and 0.78 of the simulated pl_i using all nine detectors, with a SD of 3.0–5.5 and a rms error of 3.0–5.4. The range of pl ratios for fiber 6 (worst case) was 0.86–1.15.

Figure 5

Prediction of pathlength based on calculated m₁, m₂, and m₃ values using the MPR method and all nine fibers in the data fitting to Eq. (3). Predicted divided by simulated pathlength is plotted as a function of ρ_i. Mean (x)±one SD (error bars) are given for each fiber. The triangles denote the maximum and minimum values.

For a fixed 〈cos θ〉=0.875, the calculated LDF perfusion [Eq. (7)] increases almost linearly with pl_i (Figure 6). The perfusion estimate ranges from 69–694 a.u., considering all source detector separations. Even within one fiber, the perfusion estimate can vary substantially. At ρ_i=0.46 mm, the perfusion estimate varies between 109 and 237 a.u. (Figure 7), and increases linearly with pl_i (ranging from 1.75 to 3.41 mm).

Figure 6

Perfusion estimates (Perf_i) as a function of average pathlength, depicted for all ρ_i and 〈cos θ〉=0.875. The span of optical properties is defined in Table 2.

Figure 7

Perfusion estimates vs pl_i for ρ_i=0.46 mm and 〈cos θ〉=0.875. The span of optical properties is defined in Table 2.

4.

Discussion

The aim of this study was to devise a method that could estimate the average pathlength, pl_i, migrated by photons in turbid media, at discrete source detector separations, ρ_i. The method has to be robust and applicable to a wide range of optical properties, in order to be useful in a clinical setting. It should be easy to measure and calculate the predictors, allowing the investigator to monitor pathlength variations in vivo, preferably real time. The tissue volumetric resolution of the optical properties and, hence, pathlength determination, depends on the source detector separation, ρ. Using a small ρ, local tissue inhomogeneities in the mm³ range can be revealed.² By using a larger ρ, the influence of deeper tissue structures in the cm³ domain is detected.¹⁰ We present a method, based on the spatially resolved diffuse reflectance in the 0–2 mm range, that yields $\widehat{p{l}_i}$ with a rms error of approximately 5.

The significance of being able to predict pl is illustrated in Figure 2. The longest pathlength migrated by a photon can be almost six times greater than the shortest one with a source detector separation of 2.07 mm. Source detector separations up to about one mm can yield more than three-fold variations in pl (2.8–9.2 mm).

A number of different methods for predicting the pl_i are presented in this paper. All but one are based on measuring the diffuse reflectance, at two to nine discrete detector locations. The simplest method finds $\widehat{p{l}_i}$ based on the average of the simulated photon pathlengths as a function of ρ_i; thus, rendering any measurements unnecessary. However, these predictions have a systematic overestimation of 9–19 and large SDs, probably due to the choice of validation data set, making this method less useful.

Of the three remaining estimation methods evaluated in this study, only one is independent of the absolute magnitude of the diffuse reflectance, R_i, thereby avoiding problems with absolute intensity calibration. This method uses a combined preprocessing technique of linearization and data fitting to an analytic expression, originating from diffusion theory [Eq. (3)], and results in two parameters ( m₂ and m₃ ) that merely describe the shape of R_i vs ρ_i, but not the absolute magnitude. Obviously, information is lost this way, but estimations are still more accurate than the previous method in predicting pl_i. However, the method is far too imprecise, with rms errors up to 29. This result is in agreement with those of Dam et al.; who found that predicting optical properties based on relative reflectance profiles was 5–10 times less accurate than using absolute reflectance.³

One intuitive approach previously suggested by us, is using the diffuse reflectance detected in two (neighboring) fibers as predictors.²⁴ The choice of fibers will strongly influence the accuracy of the prediction algorithm. We found slight overestimations of pl_i of up to 5.6, and rms errors up to 8.7. Obviously, the method can be expanded to incorporate more than two fibers in the analysis. The case of three different fibers improved results slightly, with overestimations ranging up to 4.5. A major drawback with this method is the inherent sensitivity to disturbances in any one of the channels/fibers during measurement.⁴

In order to condense the information from measurements, and to reduce sensitivity to data collection disturbance in one or several fibers, it is practical to fit diffuse reflectance data to an analytic expression [Eq. (3)], as mentioned previously. Using all three parameters, m₁, m₂, and m₃, in describing the R_i, both the shape ( m₂ and m₃ ) and absolute magnitude (m₁) are considered. The best accuracy in predicting pl_i based on m₁, m₂, and m₃ was seen if all nine simulated fibers were used in the data fitting. The KNN method was considered an intuitive approach in finding $\widehat{p{l}_i},$but required relatively high K values (3–10) to yield accurate estimates. Better estimates were found using a multiple polynomial regression model of the third degree. This way, $\widehat{p{l}_i}$ based on m₁, m₂, and m₃ became very accurate (overestimations up to 0.8) and precise (SD and rms error up to 5.5). The worst case range of $\widehat{p{l}_i}/p{l}_i$ was 0.86–1.15 (fiber 6).

The accuracy and precision should be viewed in light of the choice of optical parameters in the reference and validation space (Tables 1 and 2), where the parameters in the validation space were chosen to maximize the distances to the nearest combinations of optical properties in the reference space. Thus, the accuracy listed earlier for the various prediction methods are expected to be the worst case scenarios.

All the aforementioned results are based on a fixed relation between the isotropic and the anisotropic component of the phase function, but with varying 〈cos θ〉. The importance of the phase function in determining the reflectance profile, depends on the source detector separation, expressed in terms of ρμ_s ^′. ² For a high albedo and ρμ_s ^′>10, the diffusion approximation holds and the reflectance only depends on μ_a and μ_s ^′. The diffusion approximation can be extended to smaller distances, but only if the phase function is known.² Generally, for separations in the range 0.5<ρμ_s ^′<10, the reflectance depends on the first and second moment of the phase function and for smaller distances on even higher moments of the phase function.² In our study, the range of optical properties is wide in the sense that source-detector separations of 2 mm yield ρμ_s ^′ in the range 0.5–20. Since our model almost encompasses all three regimes mentioned earlier, Monte Carlo simulations were chosen to determine the reflectance. One limitation of this study is, that we did not consider other combinations of phase functions. To do this, however, the range of optical properties should be divided into smaller sets corresponding to the ρμ_s ^′ regimes mentioned earlier.

It has been proposed that μ_a and μ_s ^′ should be measured in the diffusive region and the phase function at small distances.^{1 2} This and other studies have shown that limiting the separation reduces the estimation accuracy.²⁵ Therefore, there is a trade off between spatial resolution and the estimation accuracy. Our results were obtained with a single layer model. To predict the optical properties and, hence, the photon pathlength for a layered model, requires some a priori knowledge of the layer structure and the optical properties of the layers.^{6 26}

From Figure 6, the need for correction of pathlength-related variations in the LDF perfusion estimate is obvious. While the true perfusion through the slab is kept constant, and the optical properties are varied according to Table 2 (fixed 〈cos θ〉=0.875 ), the estimated perfusion value displays a ten-fold difference between its lowest and highest value, looking at all source detector separations. Hence, at least part of an increase in the LDF perfusion signal with increasing ρ_i, traditionally attributed to sampling deeper more highly perfused areas, could simply be related to a longer pl_i, and, thus, greater probability of interaction with moving scatterers. Further, it is logical to assume that the probability of multiple Doppler shifts and nonlinear homodyne effects increase with longer pathlengths. Thus, the slope of the perfusion is expected to be somewhat underestimated for increasing pl_i. ²⁷ In a fixed fiber, at a distance relevant to LDF (ρ_i=0.46 mm), the perfusion estimate varied between 109 and 237 a.u. (Figure 7), and increased linearly with pl_i (ranging from 1.75 to 3.41 mm). It is evident that a compensation for this pathlength-related perfusion variation could improve LDF accuracy.

5.

Conclusion

In this study, we have demonstrated the substantial variations in pathlength traversed by individual photons through a turbid medium, with optical properties relevant to human skin. We present a multiple polynomial regression method that, based on spatially resolved diffuse reflectance, can predict the average pathlength as a function of source detector separation (up to 2 mm) with a rms error of about 5. The K-nearest neighbor method was the other key approach investigated. It yielded slightly less precise predictions with a rms error of approximately 7. If no preprocessing was carried out on the reflectance data, the accuracy deteriorated, but the precision was essentially retained. Finally, the average pathlength as a function of source detector separation was also predicted based on the average of all simulated photon paths, yielding gross overestimations and a rms error of up to almost 40. The results implicate that LDF perfusion estimates can be improved by assessing the pathlength. Other possible applications are reflectance spectroscopy and pulse oximetry.

Acknowledgments

This study was supported by the European Commission through the SMT4-CT97-2148 contract. Under this contract, a cooperation runs between the Universities of Twente and Groningen (the Netherlands), Linko¨ping and Malmo¨ (Sweden), Toulouse (France), the companies Perimed AB (Sweden), Moor Instruments and Oxford Optronix (UK), and the Institute of Biocybernetics and Biomedical Engineering in Warsaw (Poland). Financial support was also received from the Swedish Heart Lung Foundation, Project No. 200141208 and MedComp AB, Sweden.

Appendix A

Figure 8

(Appendix A) Top: Detail of probe design, with four adjacent fibers, including cladding. ρ_i=distance from center of source fiber to center of ith detector fiber, ρ_i,n=distance from center of source fiber to center of nth subpartition in ith detector fiber, ρ=radial distance from center of source fiber to arbitrary point in the ith detector fiber, r=core radius of fiber. Bottom: Ring-to-fiber weighting factor, W(ρ_i,n), as a function of the distance to the center of the source fiber.

Appendix B