Estimating optical properties in layered tissues by use of the Born approximation of the radiative transport equation.

We use the Born approximation of the radiative transport equation to recover simultaneously the absorption and scattering coefﬁcients in a single layer of a two-layer tissue sample from reﬂectance data. This method reduces the estimation of both optical properties to a single linear, least-squares problem. It is valid over length scales smaller than a transport mean free path and hence is useful for epithelial tissue layers. We demonstrate the accuracy of this method by using spatially resolved reﬂectance data computed with Monte Carlo simulations. © 2006 Optical Society of America OCIS codes: 170.3660, 030.5620, 100.3190 .

Determining optical properties of layered tissue structures is important for biomedical diagnostics. Several methods to determine layered tissue properties have been developed, but the vast majority of them use the diffusion approximation of the radiative transport equation (RTE). [1][2][3] These methods are limited to systems whose reduced scattering coefficients are much larger than their absorption coefficients ͑ s Јӷ a ͒ and that have a minimum layer thickness of the order of a few transport mean free paths ͓l * ϵ 1/͑ a + s Ј͔͒. However, several tissue systems exist that do not satisfy these requirements. Of prime importance is epithelial tissue, which consists of a cellular layer supported by an underlying stroma. The cellular layer has a characteristic thickness of Շ500 m, which is significantly smaller than a transport mean free path. Hence methods based on the diffusion approximation cannot be applied credibly for evaluation of epithelial-stromal transformations. Specialized techniques such as elastic scattering spectroscopy and differential path-length spectroscopy have been developed to measure properties of the epithelium alone. 4,5 While they are valuable, such techniques are unable to measure and isolate concomitant changes in the underlying stroma.
Using the Born approximation (BA) of the RTE, we formulate a novel method with which to determine optical properties of layered media. This method applies over a broad range of optical properties and to layer thicknesses of the order of a single scattering mean free path. Using the BA, we recover simultaneously the absorption and scattering coefficients in a single layer of a two-layer medium from spatially resolved reflectance data. In the calculations that follow, we assume that the optical properties of the other layer are known a priori. However, this method extends readily for more general problems. The estimates for the optical properties are the solution of an N ϫ 2 linear, least-squares problem with N ജ 2 denot-ing the number of measurements. Using Monte Carlo simulations for measured data, we show that this BA method is effective over a large range of parameter values.
Continuous-wave light transport in tissues is governed by the time-independent RTE The radiance ⌿ depends on direction , a vector on unit sphere S 2 , and position r. Here, a and s are the absorption and scattering coefficients, respectively. Scattering operator L is defined as The scattering phase function f gives the fraction of light traveling initially in direction Ј that scatters into direction . In Eq. (1), Q denotes an interior source. Green's function G͑ , r ; Ј , rЈ͒ for half-space D = ͕z Ͼ 0͖ composed of a uniform turbid medium is the solution of Eq. (1) with a and s constants and Q͑ , r͒ = ␦͑ − Ј͒␦͑r − rЈ͒, with r, rЈ D. In addition, we prescribe that ⌿ satisfy the boundary condition Here R͑͒ denotes the directional variation of the Fresnel reflection coefficient that is due to the refractive-index mismatch at the boundary. For this study we use the analytical representation of this half-space Green's function that is computed in terms of plane-wave solutions. 6,7 When the absorption and scattering coefficients are spatially heterogeneous, we represent them as a,s ͑r͒ = a0,s0 + ␦ a,s ͑r͒,

͑4͒
respectively. According to the BA, with ⌿ 0 denoting the solution of the unperturbed problem. For a function u͑ , r͒ we define the operation G D ͓u͔ used in relation (5) as We now consider the half-space D to be composed of two layers. The top and bottom layers are denoted D t = ͕0 Ͻ z Ͻ z 0 ͖ and D b = ͕z Ͼ z 0 ͖, respectively. We assume that layer thickness z 0 is known. Suppose that either the top layer or the bottom layer has absorption and scattering coefficients a,s = a0,s0 + ␦ a,s , respectively, with ␦ a and ␦ s representing constant perturbations. We assume that the other layer is unperturbed, so its optical properties are a0 and s0 . Using relation (5), we model the reflectance for the two cases, respectively. The quantities F 0 and F a,s t,b are defined as In Eqs. (8), S NA 2 denotes the set of directions within the numerical aperture of the detector. Notice that in Eqs. (8b) and (8c) the spatial integration is restricted to either the top or the bottom layer.
Suppose that data are collected at N ജ 2 detectors located at 1 ,¯, N . Rearranging terms in Eq. (7), we obtain the following linear system: . N ϫ 1 vector y has entries y n = F data ͑ n ͒ − F 0 ͑ n ͒. In fact, only N = 2 measurements are needed to compute the perturbations. However, we incorporate more measurements to help stabilize this system under the influence of instrument noise. Hence we seek the least-squares solution of Eq. (9) to obtain estimates for both ␦ a and ␦ s .
To demonstrate the utility of our inversion scheme, we generate measured data F data ͑͒ with Monte Carlo simulations of radiative transport in a twolayer model system representative of cervical epithelium. The epithelium is represented by a 500 m thick cellular layer supported by a semi-infinite stromal layer (Fig. 1). The baseline optical properties in both layers are taken to be those of normal cervical stromal tissue at = 849 nm: a0 = 0.034/ mm and s0 Ј = s0 ͑1−g͒ = 0.611/ mm. 8 The Henyey-Greenstein scattering phase function is used with asymmetry parameter g = 0.7. With these optical properties, the cellular layer has a thickness of less than l * /3. We employ six detectors, each with diameter 200 m, located at source-detector separations of 1 = 0.7 mm, 2 = 1.1 mm, 3 = 2.0 mm, 4 = 2.5 mm, 5 = 3.0 mm, and 6 = 3.7 mm. The Monte Carlo simulations were run with 30ϫ 10 6 photons and 2% Gaussian noise added to the data to typify instrument noise. For the simulations of measured data, we perturbed either a or s in either the top or the bottom layer. We assume that the optical properties in the unperturbed layer are known and are given by a0 and s0 . Changing only one parameter at a time in the simulated data tested the ability of our twoparameter inversion [Eq. (9)] to decouple and recover both ␦ a and ␦ s . The ranges of a and s perturbations used here are representative of absorption and scattering properties measured from ten patients who exhibited high-grade squamous intraepithelial lesions. 8 For a and s perturbations in the top layer [ Fig.  1(a)], the recovered optical properties are shown in Fig. 2 by filled circles and by open circles, respec-  tively. In Fig. 2(a) we show results for seven simulated measurements in which s is held fixed and a varies from 33% to 300% of a0 . The diagonal line indicates the true values of a , and the horizontal line indicates the true values of s . The results show that the recovered values of a and s for all seven test cases have relative errors within 21% and 10%, respectively. In Fig. 2(b) we show results for five simulated measurements in which a is held fixed and s varies from 60% to 140% of s0 . The diagonal line indicates the true s values, and the horizontal line indicates the true values of a . For the data in which s is perturbed, 80% to 120% of s0 , the recovered values of a and s have relative errors within 15% and 18%, respectively. Outside this range, the recovered values for s are of similar quality. However, the recovered values of a are poorer, with the maximum relative error exceeding 40%.
For a and s perturbations in the bottom layer [ Fig. 1(b)], the recovered optical properties are shown in Fig. 3. In Fig. 3(a) we show results from seven simulated measurements in which s is held fixed and a varies from 33% to 300% of a0 . The recovered values of a and s have relative errors within 25% and 5%, respectively. The recovered values of a degrade systematically as the perturbation grows. In Fig. 3(b) we show results from five simulated measurements in which a is held fixed and s varies from 60% to 140% of the baseline value. The recovered values of a and s have relative errors within 25% and 3%, respectively. The bottom-layer results are poorer than those of the top layer. This is so because the BA is less able to account for perturbations distributed over a large volume. Nonetheless, application of the BA to recover the optical properties of the bottom layer provides accurate trends of the per-turbed variables. Moreover, it does so without significant coupling between the optical parameters.
We have demonstrated the use of the Born approximation of the radiative transport equation to estimate optical properties in a single layer of a twolayer medium. This method applies to systems with length scales smaller than l * because it is based on the RTE rather than on the diffusion approximation. The BA of the RTE has been used in imaging applications of photon migration. 9,10 Our implementation of the BA of the RTE for the two-layer problem uses an analytical representation of Green's function and yields an N ϫ 2 linear, least-squares problem. Hence the associated computational costs are small.
This method requires knowledge of the optical properties of some reference medium that is close to the true medium. This a priori information is contained implicitly in the computation of the half-space Green's function. Nonetheless, this method yields results that are accurate for a broad range of optical parameters ͑6 Շ s Ј/ a Շ 25͒ that are not close to the reference medium. These data suggest that the method does not depend strongly on this a priori knowledge. Moreover, we have assumed here knowledge of asymmetry parameter g and layer thickness z 0 and of which layer is perturbed. However, the BA is not restricted by these assumptions and can be generalized to take them into account. These issues require more attention and will be subjects of future research that includes validation by use of experimental measurements of layered tissue phantoms.
The work of A. D. Kim (adkim@ucmerced.edu) is supported by the National Science Foundation through grant DMS-0504858. C. Hayakwa and V. Venugopalan are supported by the National Institutes of Health through grants R01-EB00345 and P41-RR01192.