1.

Introduction

The highly scattering nature of tissue prevents embedded inclusions from being adequately resolved in optical images.¹ Hence, differing theoretical constructs have been proposed to separate the effects of scattering and absorption, enabling one to quantify optical coefficients as a spectroscopic signature of abnormal tissue embedded in thick, otherwise normal tissue. The most widely used of these theoretical constructs are the diffusion-like models based on the diffusion approximation of the transport equation and models based on random walk theory (RWT).² These theories can be used to fit experimental data and retrieve optical properties or can be implemented in computationally intensive inverse algorithms to map the optical properties of the tissue. However, it is well established that those inverse algorithms are ill posed. Thus, results obtained from such algorithms often differ from nominal values of the medium under interrogation. For this reason, we believe it is necessary to test the forward problem with independent data sets obtained from actual experiments or Monte Carlo (MC) simulations in which the parameters involved (e.g., size of the target or the amplitude of the perturbations) are varied. This allows one to test the robustness of the theory and therefore apprehend its limits. If inverse algorithms are applied to actual experiments in the absence of such tests, the origin of observed discrepancies between the nominal values and the constructed values will be unclear: are the discrepancies due to the ill poseness of the inverse method or the inaccuracy of the theory?

For example, we have used RWT to derive the absorptive perturbation kernel of a small abnormal inclusion embedded in a tissue-like turbid medium.³ However, when we have used the standard method of integrating our perturbation kernel over the volume of a nonlocalized target, we found that simple extrapolation (i.e., integration over the volume) of the perturbation method yields inaccurate results for large inclusions. This was the motivation for deriving a corrective factor which takes into account effects of each absorptive unit (e.g., voxel) on other units and vice versa. In other words, assuming that each voxel is independent can result in significant errors. It should be noted that the perturbation kernel derived from RWT is similar to that obtained from the diffusion approximation.³ Hence, the present work could possibly solve discrepancies encountered in reconstructed images using diffusion approximation-based approaches.

In time-resolved transillumination experiments, the temporal distribution of photon paths inside the tissue is critically dependent on the optical coefficients of the tissue background, as well as those of the abnormal target(s) under investigation. The mean time that photons have spent inside an abnormal target can be seen as a quantitative parameter for imaging purposes. It is worth noting that this parameter was first introduced theoretically as a “photon hitting density” by Schotland, Haselgrove, and Leigh.⁴ The authors reduced the problem of calculating the photon hitting density to that of determining the diffusion Green’s functions. However, no comparisons of this model with experimental data or Monte Carlo simulations have been published. Though it was not analyzed, the factor would be most significant in the case of relatively large, nonlocalized absorbing targets.

Recently, Zaccanti et al.;⁵ have proposed a method to measure the mean flight time 〈t_i〉 spent by photons inside a given volume of a highly diffusing medium. The authors have substantiated their approach with Monte Carlo simulations and experiments. In the most recent paper by the same group,⁶ Monte Carlo simulations and experimental data were used to assess the accuracy of a perturbation approximation to the diffusion model for scattering and absorbing inhomogeneities.

Using the data, obtained from those researchers at the University of Florence, we show that the time dependent contrast functions proposed by Gandjbakhche et al.;⁷ could provide a direct and quantitative analysis of the mean time spent inside an abnormal inclusion. Similarly, Monte Carlo simulations and experiments aimed at an estimate of 〈t_i〉 can substantiate assumptions of the random walk (RW) model of the time-resolved contrast functions C_a(x,Δt) for absorptive inclusions,³ in particular, for the case of relatively large nonlocalized inclusions.⁸ To test the robustness of our analysis, we also present in Sec. 3 an analysis of another independent set of time-resolved experimental data on absorptive inclusions of three different sizes obtained by Painchaud et al.;⁹

2.

Monte Carlo Simulations and Experimental Setup

Both experiments and Monte Carlo (MC) simulations employ the standard transillumination geometry for a highly scattering slab of thickness T=40 mm (Figure 1). Sizes and amplitudes of absorptive perturbations were varied over a broad range. Setup details have been described recently elsewhere.⁵ This setup has been employed to experimentally determine the mean time spent by photons inside a given inclusion before being detected on the side of the slab opposite the source. The measurements were made using the original methodology of the paper,⁵ i.e., by measuring the relative variation of received time-resolved power P(Δt,x,μ_a+δμ_a)/P(Δt,x,μ_a) due to a small variation of the absorption coefficient δμ_a inside the inclusion

Figure 1

Geometric scheme assumed for theoretical model, Monte Carlo simulations, and corresponding experimental setup.

From MC simulations, both the temporal distribution of photons that have traveled from a source to a detector aligned on the opposite side of a slab, and the mean time spent by photons inside an absorptive inclusion have been derived. As is normally done, the trajectories of individual photons are simulated according to statistical rules that govern photon interactions within the medium. A sophisticated MC model based on scaling relationships that has been developed at the University of Florence¹⁰ has also been used in this research. The simulation results reported here have been obtained by assuming an asymmetry factor g equal to zero, corresponding to Raleigh scatterers. Compared to the appropriately scaled results for nonzero g values, this approach produces only a very slight discrepancy near the source and for photons with very short flight times as expected.

3.

Comparison of the RW Model With Experimental Data and MC Simulations

To compare the results of experiments and MC simulations with predictions of the RW model, we take advantage of an evident relationship between the mean time that photons spend inside the inclusion 〈t_i〉(Δt) and the corresponding absorptive contrast C_a(x,Δt) which is valid, if the extra optical depth inside inclusion Δτ=Δμ_ac〈t_i〉(Δt)⩽1 is

RW formulas for C_a(x,Δt) were recently derived and used for quantification of inclusions in tissue-like turbid slabs from time of flight experimental data.^{3 8} Within the RW model framework, for a relatively small (size d_i≪L, where L is slab thickness) abnormally absorbing inclusion, we can represent the inclusion by a set of N_a ³ independent absorbers located on a cubic lattice with spacing d_s=√2/μ_s ^′, where μ_s ^′ is the transport-corrected scattering coefficient of the background. The contrast function is⁸

Here, p(Δt), is the probability that, in the absence of the inclusion, an injected photon would reach the detector with time delay Δt, W(x,Δt) is the probability that a photon, having visited the inclusion, will reach the detector with time delay Δt, and η _eff is the effective absorptivity of an elementary absorber⁸

Taking into account that N_a ³=V(μ_s ^′/√2)³, and d_i ²=V^2/3 for relatively small inclusions with extra optical depth τ due to nonzero absorption perturbation, Δμ_a, not exceeding unity, (τ=Δμ_aμ_s ^′V^2/3/2⩽1), we obtain

Comparison of Eq. (5) with the results of MC simulations shows very good agreement, as presented in Figures 2(a) and 2(b), for the case of small spherical inclusions of volume V=524 mm ³ and positive or negative absorption perturbations embedded halfway between the source and detector in a scattering medium with μ_s ^′=0.5 mm ⁻¹. Figures 2(a) and 2(b) correspond to the positive absorption perturbation of Δμ_a=0.0025 mm ⁻¹ (μ_a ⁽⁰⁾=0) and the negative absorption perturbation of Δμ_a=−0.001 mm ⁻¹ (μ_a ⁽⁰⁾=0.001 mm ⁻¹), respectively. The exponential factor in Eq. (5) accounts for a small difference (∼6) between the results of the MC simulations for two cases. It is worth noting that the RW model predicts an increase in 〈t_i(x,Δt)〉 with decreasing Δt for the shortest time delays, resulting from the index (number of RW steps) n+1 of the probability W instead of n in the numerator of Eq. (5).⁷ This kind of behavior at small Δt is not explained by the standard diffusion model of time-resolved contrast functions as discussed, for example, by Morin, et al.;¹¹ Researchers using diffusion models generally introduce an ad hoc correction that is equivalent to one extra RW step.^{11 12}

Figure 2

Time-resolved mean time 〈t_i〉(Δt) spent by photons inside an inclusion of the volume V_i=524 mm ³ with the absorption perturbation positioned midway between the source and detector (slab thickness L=40 mm, μ_s ^′=0.5 mm ⁻¹, refraction index n=1.33 ); (a) μ_a ⁽⁰⁾=0, Δμ_a=0.0025 mm ⁻¹, (b) μ_a ⁽⁰⁾=0.001 mm ⁻¹, Δμ_a=−0.001 mm ⁻¹.

To further test the accuracy of the RW approach for estimates of absorptive contrast, we look at a correction factor for nonlocalized absorbing inclusions. We use the results of MC simulations and experiments that provide the time-integrated (CW) values of the mean time spent by photons inside an absorptive inclusion. To quantitatively analyze these results, we take into account the reduction in the average number of RW steps inside an inclusion due to the background absorption μ_a ⁽⁰⁾ (if it is not negligible). This reduction is derived, assuming, in accordance with the Central Limit Theorem, that the distribution of the random walk number of steps n _RW is Gaussian

Finally, the corresponding formula for 〈t_i〉_cw is obtained by a summation over the number of steps in the numerator and denominator of Eq. (5), and a modification of the index of the exponent

In Figure 3 we present a comparison of the calculated dependency of 〈t_i〉_cw on the lateral shift x of the inclusion center relative to the source-detector axis with the experimentally observed values and results from MC simulations presented by Zaccanti et al.;⁵ Agreement between the model and data proves to be good. It is worth noting that the shape of the lateral distribution is close to Gaussian as has been previously shown.¹³

Figure 3

Lateral distribution of the CW mean time 〈t_i〉(x), spent by photons inside an inclusion of the volume V_i=592 mm ³ with the absorption perturbation Δμ_a=0.0069 mm ⁻¹ positioned in the slab midplane; x is the distance between the inclusion center the source-detector axis.

Figure 4 shows 〈t_i〉_cw as a function of the background absorption coefficient μ_a ⁽⁰⁾. MC simulations for this case confirm the RW predictions, illustrating the accuracy of the Green functions estimates based on the RW model.

Figure 4

Dependence of the CW mean time 〈t_i〉_cw on the background absorption coefficient μ_a ⁽⁰⁾ ( μ_s ^′=0.5 mm ⁻¹, inclusion volume V_i=524 mm ³, absorption perturbation Δμ_a=0.000 05 mm ⁻¹ ).

Two more dependencies of the mean time 〈t_i〉_cw on the parameters of the inclusion itself, i.e., the target volume, V, and the amplitude of absorption perturbation, Δμ_a, are presented in Figures 5 and 6, respectively. Figure 5 shows the match between the new theory, MC simulations, and experimental data.⁵ In Figures 6(a) and 6(b), only MC results were available. Figure 6(a) corresponds to the case without the background absorption (μ_a ⁽⁰⁾=0) and positive absorption perturbations Δμ_a, while Figure 6(b) corresponds to the case of μ_a ⁽⁰⁾=0.02 mm ⁻¹ and the broader range of perturbation amplitudes, including negative values Δμ_a. For relatively small values of inclusion volume and/or absorption perturbation, both relationships are well described by the RW formulas for quasipoint absorptive inclusions,³ while for larger abnormalities good agreement is only achieved if one uses Eq. (9) with the exponential correction factor [for small μ_a ⁽⁰⁾ this factor is presented in Eq. (4)].

Figure 5

Dependence of the CW mean time 〈t_i〉_cw on the inclusion volume V (absorption perturbation Δμ_a=0.0069 mm ⁻¹ ); both the experimental data and results of MC simulations for the same conditions are shown.

Figure 6

Dependence of the CW mean time 〈t_i〉_cw on the value of the absorption perturbation Δμ_a (inclusion volume V=524 mm ³, background μ_s ^′=0.5 mm ⁻¹ ); (a) μ_a ⁽⁰⁾=0, (b) μ_a ⁽⁰⁾=0.02 mm ⁻¹ .

Additional comparisons between the results of MC simulations for the mean photon pathlengths 〈l_i〉_CW=c〈t_i〉_CW inside the absorptive inclusion of the volume V=524 mm ³ with predictions of the RW model [Eq. (9)] are presented in Table 1 for a wide range of μ_a ⁽⁰⁾ and Δμ_a. Errors from the RW expression were less than 6 (and generally much less) for optical depths τ⩽1. The error increases up to ∼10–12 for τ≈1.3 (the last row of Table 1). The dependence of absorptive contrast amplitude (which is proportional to the mean time spent inside the inclusion 〈t_i(t)〉 ) on the size of the inclusion is further tested by using independent experimental data from Painchaud et al.;⁹ Time-resolved contrast functions for three similar, purely absorptive inclusions (cylinders with equal diameter and height d) that differ only in size (d=3.6, 5.0, 7.0 mm) were presented. All measurements were reported to have been made in transillumination geometry with the inclusion halfway between the source and detector in a slab of thickness L=20 mm (scattering and absorption coefficients of the background were μ_s ^′=1.13 mm ⁻¹ and μ_a ⁽⁰⁾=0.002 mm ⁻¹ ). In Figure 7(a) we present experimental contrast amplitudes as a function of the cubed size of the inclusion for the time delay Δt=1000 ps. As expected, for the relatively large inclusion size d=7.0 mm when a simple perturbation analysis [corresponding to no correction as in Eq. (3) with η _eff =η₀=Δμ_a/μ_s ^′ ] fails, absorptive contrast deviates from proportionality to the inclusion volume V∝d³. However, if one plots the same contrast amplitude as a function of parameter Q=d³ exp(−Δμ_aμ_s ^′V^2/3/2) [Figure 7(b)] that takes into account the exponential correction factor in Eq. (4), the results are very close to linear and pass through the origin, as we expect from our model [see Eqs. (2) and (5)]. From these equations, taking into account that the RW ratio (W_n/p_n)≈0.16 for the number of steps, n, corresponding to Δt=1000 ps, we were able to estimate μ_a inside the inclusion to be μ_a=0.024±0.001 mm ⁻¹ for all considered inclusion sizes. This is close to the nominal value of μ_a=0.029 mm ⁻¹.

Figure 7

Time-resolved contrast amplitude, as a function of size of inclusion; comparison with experimental data from Ref. 9 for three absorptive inclusions, positioned at the midplane of the highly scattering slab (μ_s ^′=1.13 mm ⁻¹) of the thickness L=20 mm that differ only in size (d=3.6,5,7 mm). (a) Parameter Q related to size that corresponds to x axis of the graph, is simply d³. (b) Parameter Q includes in addition to d³ an exponential correction factor, described in Eq. (4).

Table 1

4.

Discussion

Comparisons of theoretical predictions, based on the RW model, with MC simulations and two independent experimental data sets show that the mean time spent by photons inside an absorptive nonlocalized inclusion 〈t_i〉 cannot be accurately calculated without introducing a size-dependent exponential correction factor. This factor takes into account the probability of additional absorption of the photon at the RW lattice points inside the inclusion.

Starting with the simple relationship between 〈t_i〉 and the corresponding absorptive contrast function [Eq. (2)], and using the RW derived formulas for mean time, 〈t_i〉, yields a practical framework in which the results of MC simulations and experiments may be compared. The excellent match between the simulations and experiments in this framework substantiates the use of this recently developed RW methodology^{3 8} to quantify optical characteristics and sizes of abnormalities inside a turbid slab for relatively large (nonlocalized) absorptive inclusions.

This methodology is a multi-step analysis of time-resolved data. From images observed at differing flight times, we construct the time-dependent contrast functions, fit our theoretical expressions, and compute the optical properties of the background, and those of the inclusion along with its size. This method was successfully applied to reconstruct the optical characteristics and sizes of inclusions with increased scattering and/or absorption inside several independent tissue-like phantoms.^{3 8} We used early time of flight data to estimate amplitudes of the scattering perturbations and sizes of the abnormality, while the long delayed photons are used to estimate absorption perturbations. Our exponential correction factor has been shown, for relatively large absorptive inclusions (both positive and negative), to improve the accuracy of estimated inclusion sizes and may be easily incorporated into our existing technique for analysis of absorptive contrast.³ It is worth noting that our RW-based technique for analysis of the scattering perturbations^{3 8} permits one to consider the case of the relatively large (nonlocalized) scattering inclusions.

Recently, in a pilot study, the same RW methodology was used to estimate optical parameters of the breast tumors from in vivo time-resolved data, obtained at two wavelengths.¹⁴ From the absorption coefficients at both wavelengths blood oxygen saturation was found for the tumors and surrounding tissue. We plan to continue our work on quantification of tumor characteristics from in vivo time-resolved measurements in order to better characterize different types of abnormalities, and in particular, relatively large tumors.

It should be noted that in our approach, we do not try to retrieve the whole distributions of scattering and absorption coefficients over the tissue volume, as is usually the goal of diffuse optical tomography. Instead, we consider a more limited problem, i.e., to estimate under reasonable assumptions optical parameters and sizes of a single tissue abnormality with an absorptive and/or scattering perturbation inside a quasihomogeneous region of interest inside turbid tissue.

Our analysis suggests that for accurate quantification of the absorption coefficient of an embedded target, a corrective factor which is size dependent must be included in the forward model based on diffusion-like theories. Unfortunately, this inclusion makes the expression of absorptive contrast nonlinear. As a result, development of algorithms for diffuse optical tomography, already cumbersome, becomes more challenging.