1.

Introduction

Light propagation in tissues is governed by the theory of radiative transport.¹ The radiative transport equation takes into account absorption and scattering due to inhomogeneities in tissues. Analytical solutions to this integrodifferential equation are known only for relatively simple problems.^{2 3}

For light that has propagated deeply into an optically thick medium, the radiance becomes nearly isotropic due to multiple scattering. For that case the transport equation can be replaced by the diffusion equation.¹ Solving the diffusion equation is much easier than solving the transport equation. However, the diffusion approximation is limited to situations in which the direction dependence of the radiance is nearly negligible. Hence, it does not approximate well the solution to the transport equation near collimated sources and interfaces with significant refractive index mismatch.⁴

The limitations of the diffusion approximation are problematic. Diagnostic measurements at small source-detector separations yield valuable site-specific information that is not available at larger separations.⁵ Moreover, noninvasive measurements are made at the boundary of a medium where the refractive index may be significantly different from that inside the medium. Nonetheless, it is used extensively to model diagnostic data.

Recently, there has been much progress in acquiring diverse sets of diagnostic measurements of spatial, spectral, angular, and polarization variations of light scattered by tissues.⁶ Those data show that these diverse diagnostic measurements are rich with physiological information. However, this information content is not understood well enough that this diversity can be leveraged as a gain in biomedical applications. Modeling this data accurately necessitates a more sophisticated theory than the diffusion approximation. To overcome the limitations inherent to the diffusion approximation, we seek a method to solve the transport equation.

Typically, one solves the transport equation using Monte Carlo simulations. Monte Carlo simulations trace individual photons as they propagate in and interact with the medium. Statistical averages of data collected from following a large number of these photons yield desired results. Although this method gives exact results up to statistical errors, the convergence rate is O(N^1/2) with N denoting the number of collected photons. This slow convergence rate motivates us to seek a direct numerical solution method as an alternative.

A two-layered medium is a useful model for tissues because it allows for different optical properties to be prescribed for the superficial and deep regions of tissues. For example, this difference between optical properties is necessary to model accurately light propagation through tissues consisting of epithelial and stromal tissue layers. Because most cancers originate in the epithelial layer,^{6 7} it is important to determine how light has interacted with this layer. However, light also penetrates into the stromal layer, multiply scatters within it, and interacts back and forth with the epithelial layer before emerging at the surface. The coupling between these two layers is nontrivial since photons can scatter between them any number of times.

There are several studies that have used the diffusion approximation to calculate the reflectance due to a two-layer medium.^{8 9 10 11 12} Because the diffusion approximation is not accurate for describing light near sources, its use in the top layer can be problematic. A recent study by Chang etal.¹³ has proposed using the Beer–Lambert law for the epithelial tissue layer and the diffusion approximation in the stromal layer. Alternatively, hybrid methods that solve the transport equation in the top layer using Monte Carlo simulations and solve the diffusion approximation in the bottom one have been used to study this problem.^{14 15} These methods overcome the limitations of the diffusion approximation. However, both of them use heuristic arguments, and are not derived systematically from the transport equation. Rather, they have been shown empirically to agree with the solution of the transport equation.

From the experimental point of view, different gating techniques have been proposed to discriminate between the light coming from the superficial and deep tissue layers. Time gating¹⁶ uses picosecond time-resolved techniques to select the signal arriving early to the detector thereby excluding the multiply scattered photons that have penetrated deeply. Polarization gating^{17 18} is based on the fact that single scattering in the top layer retains the polarization state of the incident light, while the degree of polarization of multiply scattered light in the deep tissue is nearly negligible. Therefore, the information contained in backscattered light from epithelial tissues is almost completely contained in the photons that retain their original polarization state.

Even though these sophisticated experimental techniques can be used to study epithelial tissue layers, the need to analyze, interpret, and predict diverse data from diagnostic measurements accurately necessitates the development of sophisticated models. We address this need by solving the transport equation. Since we shall use a direct numerical solution method, it is more complicated mathematically than solving the problem using the diffusion equation or Monte Carlo simulations. However, a direct numerical solution method provides a tool to investigate the information content available in diverse diagnostic data.

In this paper we shall study the scalar, steady-state transport equation. In particular, we shall calculate the solution of the transport equation for a two-layered medium in which a finite slab resides on top of a half space. Using the theory of Green’s functions, we replace the two-layer problem by an equivalent slab problem. The key to this equivalent slab problem is the prescription of an alternate boundary condition at the bottom of the slab that takes into account exactly the multiple scattering of light in the lower layer. This boundary condition is given in terms of the surface Green’s function for the half space. The surface Green’s function is related directly to the volume Green’s function. We use the method due to Kim¹⁹ to compute the volume Green’s function for the half space. By solving this equivalent slab problem rather than the two-layer problem, we reduce the computational domain to the top layer. Hence, this result is useful for developing methods to probe epithelial tissue layers.

For the case in which the finite slab is optically thin, we compute an asymptotic solution to the equivalent slab problem. This analytical solution gives a systematic representation of interaction between the top and bottom layer. We shall show that this thin layer approximation works well for all angles except those that are near-grazing.

The paper is organized as follows. In Sec. 2 we discuss the two-layer problem and its direct numerical solution in one spatial dimension. In Sec. 3 we discuss the equivalent problem in the top slab using the alternate boundary condition at the bottom of the slab. Included in this discussion is the computation of the surface Green’s function for the lower medium and the direct numerical solution of the problem. In Sec. 4 we compute an asymptotic approximation to the equivalent slab problem for the case in which the slab is optically thin. In Sec. 5 we discuss the extension of this work to three-dimensional problems. Section 6 contains numerical results demonstrating the validity of this theory. Section 7 is the conclusions. The Appendix gives the method used to calculate plane wave solutions that are used throughout this discussion.

2.

The Two-Layer Problem

The radiance Ψ is the radiant power per unit solid angle per unit area perpendicular to the direction of propagation. It depends on the unit direction vector ω and the position vector r. The radiance in tissues is governed by the transport equation

We seek the solution of (1) in a two-layered medium in which a finite slab 0<z<z₀ is situated on top of the half space z>z₀. In particular, we wish to calculate the light backscattered by this medium. The optical properties: the absorption coefficient μ_a, the scattering coefficient μ_s, and the scattering phase function f, are homogeneous in each layer. However, the properties in one layer may be different from the other layer. To acknowledge this difference, we denote the optical properties in the top layer by (μ_a1,μ_s1,f₁), and those in the lower layer by (μ_a2,μ_s2,f₂). For simplicity, we assume index-matched boundaries. This discussion extends readily for the case in which boundaries are index mismatched.²⁰

2.1.

The One-Dimensional Problem

To begin this discussion, we study the one-dimensional problem. For this problem, a plane wave that illuminates normal to the top boundary of the slab is the only source. A sketch of this problem is shown in Fig. 1.

Figure 1

A sketch of the two-layer problem.

For this problem the solution Ψ is independent of the transverse spatial coordinates x and y. In addition, because we consider a plane wave incident normally on the medium, the solution is axisymmetric and hence, independent of φ. We represent the solution Ψ of this one dimensional, two-layer problem as

Eq. (3)

$$\mathrm{\Psi }(\nu ,z)=\{\begin{array}{cc}{\mathrm{\Psi }}_1(\nu ,z),\hfill & 0<z<z_0,\hfill \\ {\mathrm{\Psi }}_2(\nu ,z),\hfill & z>z_0,\hfill \end{array}$$

with Ψ₁ satisfying

Eq. (4)

$$\nu {\partial }_z{\mathrm{\Psi }}_1+{\mu }_{a1}\mathrm{\Psi }=-{\mu }_{s1}{\mathrm{\Psi }}_1+{\mu }_{s1}{{\displaystyle \int }}_{-1}^1{p}_1(\nu ,{\nu }^{\prime })\times {\mathrm{\Psi }}_1({\nu }^{\prime },z)d{\nu }^{\prime },0<z<z_0,$$

and Ψ₂ satisfying

Eq. (5)

$$\nu {\partial }_z{\mathrm{\Psi }}_2+{\mu }_{a2}{\mathrm{\Psi }}_2=-{\mu }_{s2}{\mathrm{\Psi }}_2+{\mu }_{s2}{{\displaystyle \int }}_{-1}^1{p}_2(\nu ,{\nu }^{\prime })\times {\mathrm{\Psi }}_2({\nu }^{\prime },z)d{\nu }^{\prime },z>z_0\text{.}$$

The functions p_1,2 are related to f_1,2 by

Eq. (6)

$$p_{1,2}(\nu ,{\nu }^{\prime })={{\displaystyle \int }}_{-\pi }^{\pi }{f}_{1,2}[\nu {\nu }^{\prime }+\sqrt{(1-{\nu }^2)(1-{\nu }^{\prime 2})}\times \text{cos}(\phi -{\phi }^{\prime })]d(\phi -{\phi }^{\prime })\text{.}$$

At z=0 we prescribe the boundary condition

Eq. (7)

$${\mathrm{\Psi }}_1(\nu ,0)=\frac{F}{2\pi }\delta (\nu -1),0<\nu \le 1\text{.}$$

This boundary condition corresponds to a plane wave with flux F incident normally on the boundary at z=0. At the interface located at z=z₀, we prescribe that the radiance must be continuous over all directions

Eq. (8)

$${\mathrm{\Psi }}_1(\nu ,z_0)={\mathrm{\Psi }}_2(\nu ,z_0),-1\le \nu \le 1\text{.}$$

In addition, we prescribe that Ψ₂→0 as z→∞.

3.

The Equivalent Problem in the Slab

To compute the solution to the two-layer problem, one must solve two transport equations for the slab and for the half space. Those two solutions must satisfy the matching condition at z=z₀ given by (8). Because most early precancerous tissues develop in the epithelial layer, we would like to focus our analysis on the top layer. However, diagnostic measurements taken at the boundary z=0 include light that has penetrated through the top layer into the bottom one, and scattered from there up into the slab. From there, light scatters back and forth between the two layers any number of times before emerging as backscattered light. This interaction is taken into account by the coupling of Ψ₁ and Ψ₂ in (8).

Suppose that the optical properties of the bottom layer are known. Under this assumption, we introduce a boundary condition at z=z₀ that takes into account light backscattered by the bottom layer. In this way we are able to replace the two-layer problem by an equivalent problem for a slab corresponding to the top layer.

Consider the half space problem in which the top layer has been removed. A light source Ψ(ν,z₀) for 0<ν⩽1 is incident on the half space at the boundary z=z₀. There are no other sources. The radiance backscattered by this half space is given in terms of the surface Green’s function G_s as³

Equation (14) gives the radiance backscattered by the lower medium due to light incident from above it. We can interpret this formula as a diffuse reflection law (see Fig. 2). Rather than solve the problem in the lower layer, we can use (14) as a boundary condition for Ψ₁ at z=z₀ and for −1⩽ν<0. Hence, we replace the original two-layer problem by one in the finite slab 0<z<z₀:

Figure 2

A sketch of the equivalent slab problem in which we replace the transport equation in the lower half space by an alternate boundary condition given in terms of the surface Green’s function G_s.

The boundary condition given by (18) can be written formally also in terms of the Chandrasekhar’s scattering S function.² For the case that scattering is isotropic, this boundary condition can be written in closed-form in terms of Chandrasekhar’s H-function² or using the Green’s function given by Case and Zweifel.³

4.

Thin Layer Approximation

Suppose the top layer is very thin compared with the scattering mean free path z₀≪l_s=1/μ_s. Then, we can compute an analytical approximation to the solution of the equivalent slab problem. We introduce the nondimensional variables

The thin layer approximation corresponds to the limit in which ε→0. In addition, we assume here that α is much smaller than ε. This assumption is not necessary, but simplifies the following analysis. Furthermore, it is consistent with typical optical properties for tissues (e.g., the absorption coefficient is approximately one order of magnitude smaller than the scattering coefficient).

We represent Ψ as the asymptotic expansion

To O(ε) we obtain

One can continue this computation in a similar manner to obtain higher order corrections. For example, one can consider the limit in which α=O(ε²). Hence, by computing the asymptotic expansion of Ψ up to O(ε²), one takes into account absorption in the top layer. This computation involves only a sequence of elementary calculations. However, the analytical expressions become complicated, so we do not discuss them here.

The thin layer approximation yields a power series in z¯. This power series is related to the successive scattering or Neumann series. For that case the leading order term in that series is given by the Beer–Lambert law.¹ Hence, it is proportional to a exponential function decaying with respect to depth. Consequently, higher order corrections involve exponential functions also. By expanding these exponential functions as Taylor series and truncating them in a manner that is consistent with the balance of scales given in (22), we can show that the thin layer approximation and the successive order scattering series are equivalent asymptotically. Hence, we have shown that the model due to Chang etal.¹³ in which they use the Beer–Lambert law for the top layer follows readily from this systematic treatment of the transport equation. Moreover, this perturbation analysis provides an estimate for the error incurred by applying this approximation. Hence, this analysis provides additional insight into the range of validity of that model.

5.

Extension to Three-Dimensional Problems

Even though we have discussed and demonstrated this theory for the one-dimensional problem, it extends easily to the three-dimensional problem. For that case the equivalent slab problem is

Rather than solve the slab problem in the physical domain, we solve for its Fourier transform with respect to the transverse spatial coordinates, x and y. The Fourier transform of Ψ with respect to x and y is

Kim¹⁶ gives a method to compute the volume Green’s function for the half space in this Fourier domain. Using that result we compute G^_s. Hence, we are able to solve the equivalent slab problem for each transverse spatial frequency pair (k_x,k_y) over a discrete spectrum. In fact, the equivalent slab problem for each (k_x,k_y) pair is a one-dimensional problem just like the one we have discussed earlier. Once that spectrum is computed, we recover the solution in the physical domain using discrete inverse Fourier transforms. Because the slab problem for each spatial frequency pair is decoupled from the others, this computation is easily implemented on a concurrent computing system thereby increasing the efficiency of this calculation.

6.

Numerical Results

We have used the Henyey–Greenstein scattering phase function to compute our numerical solutions. It is given by

Its only parameter is the anisotropy factor g. By computing (6) using (43), we determine the corresponding scattering function in the one-dimensional problem to be

For the top layer, we have set the absorption coefficient to be μ_a1=0.02 mm ⁻¹ and the scattering coefficient to be μ_s1=6.50 mm ⁻¹. For the bottom layer, we have set the absorption coefficient to be μ_a1=0.01 mm ⁻¹ and the scattering coefficient to be μ_s1=6.00 mm ⁻¹. For both layers, we have set the anisotropy factor to be g=0.80. We have taken these values from the study by Kienle etal.⁹

6.1.

Plane Wave Source

In Fig. 3 we show numerical results for the case in which a plane wave is incident normally on the two-layer medium. The plane wave solutions that we calculated numerically used M=64 Gauss–Legendre quadrature points. We show the solutions calculated for the two-layer problem (“x” symbols), the equivalent slab problem (solid curves) and the thin layer approximation (dashed curves). These solutions have been calculated for different top layer thicknesses: (a) z₀=0.1 mm, (b) z₀=0.5 mm, (c) z₀=2.0 mm, and (d) z₀=6.0 mm.

Figure 3

The radiance backscattered by a two-layer medium due to a plane wave normally incident on the medium. The optical properties of the top slab are (μ_a1,μ_s1,g₁)=(0.02 mm ⁻¹,6.5 mm ⁻¹,0.80) and the optical properties of the bottom half space are (μ_a2,μ_s2,g₂)=(0.01 mm ⁻¹,6.0 mm ⁻¹,0.80). The thickness of the top slab is (a) z₀=0.1 mm, (b) z₀=0.5 mm, (c) z₀=2.0 mm, and (d) z₀=6.0 mm. The x symbols are for the solution to the two-layer problem, the solid curve is for the solution of the equivalent slab problem and the dashed curved is the thin layer approximation.

In all of these results we observe that the solutions to the two-layer problem and the equivalent slab problem are indistinguishable. Because the theory underlying the equivalent slab problem is exact, the only quantitative differences between these results are below the double precision arithmetic used to calculate them. The thin layer approximation gives good agreement for z₀=0.1 mm corresponding to ε=0.65. For the cases with larger slab thicknesses, we observe that the thin layer approximation is not accurate. The accuracy of this approximation breaks down especially for near-grazing angles for which ν≈0. This break down occurs because higher order scattering events, which are not taken into account by the thin layer approximation, are significant for near-grazing angles.

6.2.

Gaussian Beam Source

Using the extension of the theory discussed in Sec. 5, we have computed the solution of the two-layer problem and the equivalent slab problem due to a Gaussian beam incident normally on the medium. For that case, we can reduce the number of independent variables by working in cylindrical coordinates. Because the beam is incident normally on the medium, the problem is axisymmetric. Hence, the reflectance is defined by the Hankel transform

Eq. (45)

$$R(\rho )=\frac{1}{2\pi }{{\displaystyle \int }}_0^{\infty }\hat{F}(k)J_0(k\rho )kdk,$$

with

Eq. (46)

$$\hat{F}(k)={{\displaystyle \int }}_{-\pi }^{\pi }{{\displaystyle \int }}_{-1}^0\hat{\mathrm{\Psi }}(\nu ,\phi ,z=0;k)\nu d\nu d\phi \text{.}$$

Here ρ=(x²+y²)^1/2 and k=(k_x ²+k_y ²)^1/2. Additional details on the implementation of the direct numerical solution including the numerical calculation of (45) are contained in Ref. 23.

In Fig. 4 we show the radial dependence of the reflectance due to a Gaussian beam incident normally on the medium. The beam width has been set to 2 mm. The plane wave solutions that we calculated numerically used M=16 for the Gauss product quadrature rule. We show the solutions calculated for the two-layer problem (x symbols) and the equivalent slab problem (solid curves). These solutions have been calculated for different top layer thicknesses: (a) z₀=0.5 mm, (b) z₀=2.0 mm, and (c) z₀=6.0 mm. For all cases we observe that the solutions to the two-layer problem are indistinguishable from those to the equivalent slab problem.

Figure 4

The reflectance by a two-layer medium due to a Gaussian beam normally incident on the medium. All parameters are the same as for Fig. 3. The beam width is 2.0 mm. The thickness of the top slab is (a) z₀=0.5 mm, (b) z₀=2.0 mm, and (c) z₀=6.0 mm. The x symbols are for the solution to the two-layer problem, and the solid curve is for the solution of the equivalent slab problem.

7.

Conclusions

We have discussed the direct numerical solution of the transport equation for a two-layer medium. Using the theory of Green’s functions, we have also reduced the two-layer problem to an equivalent slab problem. This equivalent slab problem uses an alternate boundary condition given in terms of the surface Green’s function for the lower layer. Hence, it depends only on the optical properties of the lower layer. The theory behind the equivalent slab problem is exact.

We have discussed in detail the one-dimensional problem in which a plane wave is incident normally on the top layer. Numerical results confirm that the solutions to the two-layer problem and the equivalent slab problem are equal as expected. This theory extends easily to three-dimensional problems through the use of Fourier transform methods. We have shown numerical results for the problem in which Gaussian beam is incident normally on the medium. Again, the solutions given by the two-layer problem and the equivalent slab problem are indistinguishable.

For the case in which the top layer is optically thin, we have derived an asymptotic approximation. The backscattered radiance given by the thin layer approximation agrees well with that given by the numerical solution to the two-layer problem. However, this approximation breaks down when the thickness of the top layer increases. Our numerical results confirm this. The accuracy breaks down most dramatically for near-grazing angles. This discrepancy is because backscattered light at near-grazing angles is made up of multiply scattered light that is not taken into account by this asymptotic theory.

These results suggest that this equivalent slab problem should be extremely useful for studying the epithelial tissue layer when the optical properties of the stromal layer are known.

Appendix: Plane Wave Solutions

Plane wave solutions are a general class of solutions that take the form Ψ(ν,z)=e^λzV(ν). Substituting this solution form into the transport equation yields

We normalize the plane wave solutions according to

The set of eigenfunctions {V_j(ν)} are complete over the full range. Hence, any arbitrary function can be expanded in terms of plane wave solutions.

We cannot calculate the plane wave solutions analytically, in general. Hence, we calculate them numerically using the discrete ordinate method. We approximate the integral operation in (A1) using the Gauss–Legendre quadrature rule

Because of the Gauss–Legendre quadrature rule is symmetric about the origin, the numerically calculated plane wave solutions obey the symmetry property. Hence, λ_−j=−λ_j and V_−j(ν_m)=V_j(−ν_m)=V_j(ν_M−m+1) for j=1,…,M. We order and index the eigenvalues as

For the three-dimensional problem given in (40), one must include the φ dependence of the radiance. The method to calculate plane wave solutions follows in a similar way as shown above for the one-dimensional problem. We use a product Gaussian quadrature rule to approximate the integral operations: