Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues

We introduce a novel and efficient method to provide solutions to inverse photon migration problems in hetero- geneous turbid media. The method extracts derivative information from a single Monte Carlo simulation to permit the rapid determination of rates of change in the detected photon signal with respect to perturbations in background tissue optical properties. We then feed this derivative information to a nonlinear optimization algorithm to determine the optical properties of the tissue heterogeneity under examination. We demonstrate the use of this approach to solve rapidly a two-region inverse problem of photon migration in the transport regime, for which diffusion-approximation-based approaches are not applicable.

Recently there has been heightened interest 1 -3 in quantifying both the structure and the function of biological tissue on millimeter to submillimeter spatial scales. For tissues, such dimensions are intermediate between the coherent regime, covered by confocal and optical coherence microscopies, and the diffuse regime, covered by standard photon migration techniques. In this intermediate regime the coherent character of the light is partially lost, but the photon directions are not suff iciently random to permit the application of the diffusion approximation to radiative transport theory. Such a situation arises when one is probing tissues (a) in which optical absorption is comparable with or greater than scattering and (or) (b) over length scales that are comparable with the transport mean free path. However, if appropriate optical models can be employed in this regime, which we term the transport regime, the potential for determining optical scattering and absorption properties on a spatial scale of hundreds of micrometers to a few millimeters can be realized. Such an advance would provide the basis for developing quantitative noninvasive diagnostic tools that are complementary to methods that operate in either the diffuse or the coherent regime.
In this Letter we present a new approach to solving inverse photon migration problems that is applicable to both transport and diffuse regimes. The novelty of this approach lies in the rapid extraction of derivative information by use of the perturbation Monte Carlo (pMC) method. 4,5 We use this derivative information within a gradient-based nonlinear optimization algorithm to determine the optical properties of the tissue heterogeneity under consideration. This technique thus allows for the resolution of inverse photon migration problems in computational times that are only marginally larger than for a single Monte Carlo (MC) simulation alone. This gain in computation efficiency is achieved because pMC analyses derived from a single MC simulation replace the repetitive MC simulations that would otherwise be required for providing the derivative information that is necessary for resolution of the inverse problem. In addition, the positive correlation in the pMC estimates captures differential effects with much greater accuracy than if independent MC simulations were used. We stress that the pMC approach, as with any other MC simulation of light transport, permits the accommodation of arbitrary geometric conf igurations and eliminates the limitations that are inherent in the diffusion approximation. As a result, this inversion scheme is applicable down to length scales that approach the single-scattering mean free path.
Our approach is depicted schematically in Fig. 1. First we adopt a f ixed source -detector geometry and background tissue optical properties with which to generate and save a basic set of biographies for detected photons by use of a conventional MC simulation. Second, we employ a pMC algorithm that rapidly predicts the signal that would be detected as a result of any number of specified changes in the spatial distribution of optical properties [e.g.,m a ͑r͒,m s ͑r͒,ĝ͑r͒] relative to the background [e.g., m a ͑r͒, m s ͑r͒, g͑r͒]. We obtain such predictions by reweighting the photon biographies, using ratios of the perturbed to the unperturbed probability-density functions. If we consider the background problem to be a homogeneous sample with optical absorption and scattering coeff icients m a and m s , respectively ͑m t m a 1 m s ͒, we generate photon biographies by sampling intercollision distances s from the probability-density function m t exp͑2m t s͒. The weight of each photon is reduced by the factor ͑m s ͞m t ͒ at each collision point, and thus j A͑m s ͞m t ͒ k is the weight of each photon that reaches a detector after k collisions, where A is a factor that accounts for the refractive-index mismatch between the tissue and the detector. If a perturbing region replaces a when the detected photon experiences j collisions and traverses a total length S in the perturbed region. We note that Eq. (1) can be modif ied to accommodate other characteristics of the perturbed region including the location of its boundaries and moments of the scattering phase function. Thus these other parameters can, in principle, also be determined by the pMC method. The pMC algorithm predicts the measured signalĵ that results from a perturbation in the background optical properties, and our goal is to use this information to guide an optimization algorithm to resolve the inverse problem. We accomplish this by demonstrating that MC estimates of ͑≠ĵ͞≠a a ͒, ͑≠ĵ͞≠a s ͒, and higher-order derivatives are correctly estimated by a closed-form differentiation of Eq. (1). This derivative information is used by a gradient-based nonlinear optimization algorithm for the determination of the optical properties in the perturbed region that best match the detected photon signal (in a x 2 sense). The idea for estimating derivatives based on a single MC simulation for neutron transport applications was described previously by both Dejonghe 6 and Rief. 7 The technique has been used principally for performing sensitivity analyses in forward MC simulations, although Hall has applied "differential Monte Carlo" to the adjustment of material properties to fit reaction rates in a specific neutron transport problem. 8 To demonstrate the performance of the technique, we consider optical property changes within a layered tissue model. When layered tissues are probed by use of noninvasive optical methods, the goal is to identify the optical properties of each layer from a measured signal. The tissue that we consider in this demonstration is composed of a thin (0.5-mm) upper layer above a much thicker (19.5-mm) base layer. We treat spatially homogeneous perturbations in either the top or the bottom layer relative to a background tissue that possesses spatially homogeneous optical absorption, scattering, and anisotropy coefficients of m a 0.034 ͑1͞mm͒, m s 6.11 ͑1͞mm͒, and g 0.9, respectively. We use the Henyey -Greenstein phase function to describe the angular scattering distribution. These properties are representative of normal cervical stromal tissue at l 849 nm, as reported by Hornung and co-workers. 9 A MC database that describes the transport of 3 3 10 6 photons through this homogeneous background tissue is generated by use of these properties. The perturbed region is assumed to be either one of the two layers, and data required for determiningĵ are tallied and saved. The measured data are obtained with a two-region MC simulation that provides real and imaginary components of the simulated measurement signal in the frequency domain. 10 Measurements are produced at 20 equispaced frequencies in 0-1000 MHz for three detectors that span 0 -3 mm, and a Box-Müller transformation is utilized to add 2% noise to these simulated measurements.
Using an initial guess equal to the background optical properties, we make MC estimates of the signal derivative with respect to m a and m s and use them to guide a two-parameter Levenberg -Marquardt algorithm to minimize the difference between the measured data and the perturbation approximation. Although simultaneous perturbations of both m a and m s can be determined by our solution method, we chose test cases that consisted of different relative changes of either m a or m s in one of the two layers. These test cases highlight the effectiveness of the two-parameter optimization algorithm in identifying changes in one of the two parameters. For all test cases the pMC approach converges in fewer than 11 iterations, with completion of each iteration requiring ϳ3 min on a computer with a 500-MHz Pentium III processor.  respectively. Both cases produce excellent decoupling of the changes in m a from m s . Specifically, although m a varies by a factor of 9 over the whole range tested, the retrieved value for m s never deviates from the true value by more than 63%. Moreover, the results are surprisingly good for the top-layer perturbation, even though its thickness is only 0.5 mm (fewer than three scattering mean free paths). The top layer exhibits larger error bars because the percentage of collisions that occur in the upper layer is only ഠ17% of the total as determined from the photon biographies. Figures 3(a) and 3(b) display results of similar quality for m s perturbations ranging from 70% to 130% of the background optical properties in the top and bottom layers, respectively. These results also exhibit excellent decoupling of the changes in m s from m a .
In summary, we have developed a general method for solving important inverse problems in biomedical optics. Its novelty lies in the use of a pMC method to perform a sensitivity analysis of the photon migration process with respect to spatially heterogeneous tissue optical properties. The use of the pMC method eliminates the need for running many separate, independent MC simulations to extract the dependence of the solution on the optical parameters under study. The latter approach would not only be orders of magnitude more costly but would be much less accurate as well. It is important to appreciate that the model geometry considered here dictates a solution in the transport regime, as it uses data that are retrieved with source -detector separations smaller than 3 mm, for which the signal captured cannot be adequately described by standard diffusion-approximation-based approaches. Moreover, whereas the frequency-domain test case was chosen to mimic best the experimental study conducted by Hornung and co-workers, 9 results of equivalent quality have been obtained in the case of continuous illumination and detection. We believe that this technique will prove to be valuable in resolving a wide range of inverse problems in photon migration that are cumbersome or inaccessible when conventional approaches that make use of full transport solutions are used.