Reflectance measurements of layered media with diffuse photon-density waves: a potential tool for evaluating deep burns and subcutaneous lesions

The basic principles of a non-contact, near-infrared technique for the mapping of layered tissues are discussed theoretically and verified experimentally. The propagation properties of diffuse photon-density waves in tissues depend on the optical properties of the tissue. When a layered medium is irradiated by amplitude modulated light, the difference in optical properties between the layers is evident in the phase and amplitude of the diffuse reflection coefficient, which is a result of the interference of the partial waves propagating in the different layers. Thus, diffuse photon-density waves are applicable to the analysis of the structure of layered tissue. The probing depth is determined by the modulation frequency of the incident light. For modulation frequencies between several hundred megahertz and a few gigahertz, this allows us to analyse the properties of muscle tissue of up to 4-8 mm below the surface. Experimental results based on chicken breast muscle are given. As an example, the technique might be of use for evaluating the depth of necrosis and the blood volume fraction in deep burns.


Introduction
Light impinging on a turbid medium is only partly reflected at the surface. The major body of photons penetrates into the medium, thus probing it to a certain depth, before being absorbed or scattered out. Measuring the reflectance from a turbid medium, therefore, not only yields information about the very surface, but also down to a certain depth in the material. This effect forms the basis for optical tomography techniques used for the investigation of the composition, structure and properties of biological tissue. As these techniques can be applied non-invasively, they are ideally suited for diagnostic purposes in medicine. In the wavelength range where haemoglobin absorption is low, i.e. for red/infrared light, photons can travel up to several millimetres through tissue, gathering information about its constituency. Optical tomography techniques can be divided into two main groups, according to how light scattering in tissue is handled: techniques based on spatial coherence, such as confocal microscopy (Schmitt et al 1994) and low-coherence tomography (Schmitt and Knüttel 1997), and techniques based on diffusive propagation of photons, such as time-resolved (Patterson et al 1989, Arridge et al 1991, Schotland et al 1993, Jacques et al 1995 and diffuse photon-density wave methods (Fishkin et al 1991, Svaasand et al 1993, Patterson 1995. 0031-9155/99/030801+13$19.50 © 1999 IOP Publishing Ltd Techniques belonging to the former group show a good spatial resolution, typically in the range of 5-20 µm for low-coherence microscopy (Fercher 1996), and in the submicron range for confocal microscopy. The resolution in the transverse direction is limited by the spatial coherence of the light. The depth resolution limit of confocal microscopy is given by diffraction, whereas in the case of low-coherence microscopy it depends on the temporal coherence of the source. On the other hand, the depth in tissue that can be analysed by both lowcoherence and confocal microscopy is limited by the mean free path length between scattering events, as scattering leads to a loss of spatial coherence. The transport mean free path, i.e. the depth where an average photon has completely lost its information about the original direction of motion, typically lies in the range of 0.2-1 mm for near-infrared light (Cheong et al 1990), and gives the order of an upper limit for the probing depth by coherence techniques.
The other group of optical tomography techniques uses the diffusion of photons in a highly scattering environment explicitly to investigate the tissue. One subgroup, the time-resolved techniques, measures the time of flight of the photons on their way through the tissue. Of particular interest are those photons that have crossed the tissue on an almost straight path, the so-called ballistic and quasi-ballistic photons. Their measurement is extremely difficult; firstly because scattering reduces their number rapidly, and secondly because an extremely high time resolution is necessary (Jacques et al 1995). The transient of the time distribution for later arriving photons, that means photons which have undergone multiple scattering, varies more slowly, and its measurement can be used for the characterization of tissue (Chance 1995). The Fourier transform relates the time-resolved techniques to a second subgroup, which measures the propagation properties of so-called diffuse photon-density waves. These are excited by harmonically modulated light, and information about the tissue can be gathered by analysing the amplitude and phaseshift relative to the source signal after the light has travelled over a certain distance. The penetration depth for such waves has typical values of 3-9 mm for light in the near-infrared range (Cheong et al 1990). The spatial resolution corresponds to the scattering coefficient of the tissue, and is comparable to the transport mean free path. Techniques based on photon diffusion thus probe far deeper into the tissue than methods based on spatial coherence, but at the cost of considerably reduced spatial resolution. Imaging of spherical objects by diffuse photon-density waves has been investigated, first for objects embedded in a homogeneous medium (Boas et al 1993) and later for objects in a layered medium of spherical symmetry (Yao et al 1996). There is hope that this technique can be applied to the diagnosis of breast carcinoma (Fishkin et al 1997). Recently, interest in the analysis of plane layered media has evolved (Alexandrakis et al 1998, Kienle et al 1998, Spott et al 1998.
The investigations of the propagation properties in layered tissue published so far have concentrated on infinitely thin pencil and point sources. In this study, in contrast, a plane wave source is assumed, which is a suitable model when the dimension of the homogeneously irradiated area on the surface of the tissue is large in comparison with the depth to which the light penetrates into the skin. While the necessary size of the irradiated area naturally puts a strong limit to the spatial resolution in the lateral direction, such measurements are less susceptible to local variations in the tissue, and thus yield more precise global values, as was shown in previous work with steady-state measurements (Svaasand et al 1995, Spott et al 1997. The mathematical formalism used here is a straightforward extension of the steady-state model into the frequency domain. The underlying equations are assembled in the appendix. The rest of the main part of this paper is structured as follows. In section 2 the observable phenomena, namely variations of the amplitude and phase of the diffusely reflected light, depending on the properties of the layered structure, are investigated theoretically, based on the formulae given in the appendix. It shows that the different propagation properties of the photon-density waves in the layers with different optical properties lead to interesting interference phenomena. These become manifest in that neither the amplitude nor the phaseshift of the diffusely reflected light necessarily lie in the range of the corresponding values for uniform, semi-infinite media with the same properties as the different layers. Section 3 presents the results of a simple ex vivo experiment. Reflectance measurements are taken from partially coagulated chicken breast muscle and prove the theoretical predictions to be qualitatively correct. In the last section of this paper, the results are discussed and a conclusion is given.

Theory
A harmonically modulated light wave, which propagates through a turbid medium, sets up an energy density = DC + AC . DC is the steady-state value of the energy density, and AC a harmonically varying contribution. When dealing with diffuse photon-density waves, it is this latter fraction which is of interest, because it describes the wave-like character of the energy density. In the appendix the governing equations based on the diffusion approximation to the Boltzmann transport equation are given. The diffusion approximation leads to a second-order differential equation of Helmholtz type (A.7) in terms of the fluence rate ϕ. The fluence rate is defined as the integral of the radiance L over all solid angles (International Organization of Standardization 1992): In a similar manner to the energy density, the fluence rate can be split up into a steady-state and a time-varying fraction. As the energy density and the fluence rate are proportional to each other, so are their time-varying fractions: where the proportionality factor c is the speed of light in the medium. A solution of the diffusion equation (A.7) for a plane wave travelling in the positive x-direction in a homogeneous, sourcefree medium is given by where ϕ 0 is the fluence rate at x = 0 and δ c the effective complex penetration depth. For a harmonic wave ϕ AC ∼ exp(−iωt), δ c takes a complex value, and the solution (3) can be rewritten as The exponential term in equation (4) is characterized by the penetration depth δ AC (ω): and determines the attenuation of the amplitude of the oscillation. The sine term in equation and describes the propagation of the wave. The propagation properties of diffuse photon-density waves are thus characterized by the penetration depth δ AC and the group velocity v. Both quantities depend on the modulation frequency, and it is instructive to study their behaviour by an investigation of their Taylor series expansions. This results in a decisive quantity for their behaviour, given by the product ωτ a , where τ a = (cµ a ) −1 is the loss relaxation time, i.e. the mean time of flight between two absorption interactions; µ a is the absorption coefficient. For small values of ωτ a , the penetration depth can be approximated by (ωτ a 1) where the penetration depth at zero modulation frequency, δ DC , is defined by Thus, the higher the modulation frequency, the lower the penetration depth and thereby the higher the attenuation of the amplitude of the oscillation. For the phase velocity the following expansion is found: where the value in the zero-frequency limit, v DC , is given by The phase velocity therefore increases with increasing modulation frequency. The series expansions lead to the conclusion that changes in the attenuation and the phase velocity, and thereby changes in the propagation properties of the diffuse wave, become of importance once ωτ a ≈ 1. This condition motivates the definition of the loss relaxation frequency f a : For low-loss tissues with a red/near-infrared absorption coefficient µ a = 10-20 m −1 , which applies to galliform muscle tissue; for example, a loss relaxation frequency of f a = 0.3-0.7 GHz is found. The absorption coefficient values for green/yellow light are µ a = 130-350 m −1 with a corresponding relaxation frequency f a = 4-12 GHz (Cheong et al 1990). Figure 1 shows the tissue structure investigated in this study. A plane layer of thickness d sits on top of a semi-infinite medium of different optical properties, and is bounded by air on the upper side. When this structure is irradiated by amplitude modulated light, a diffuse photon-density wave is excited in the tissue. The different propagation properties in the tissue layers become visible in the light which is scattered out of the tissue, i.e. diffusely reflected, due to the interference of the partial waves propagating in the two layers. The amount of diffusely reflected light, the diffusely reflected flux, is given by −j | x=0 + , i.e. the negative of the flux in the tissue next to the boundary. The negative sign indicates that the reflected light is directed out of the tissue, in the negative x-direction. The diffuse reflection coefficient γ is defined as the ratio between the diffusely reflected flux and the incident flux P 0 and can, as a complex quantity, be written as where the amplitude |γ | of the reflection coefficient is the ratio of the amplitude of the backscattered to the incident flux, and the phase γ of the reflection coefficient is the phaseshift of the reflected relative to the incident signal.
As an example of the frequency dependence, the amplitude and the phase of the reflection coefficient of a two-layered structure versus the thickness of the upper layer and the modulation frequency are given in figure 2. For the upper layer, an absorption coefficient of µ a,1 = 10 m −1 and a reduced scattering coefficient of µ s,1 = 700 m −1 were assumed. These values correspond to the optical properties of chicken breast muscle for near-infrared light (Cheong et al 1990,  Marquez et al 1998). The same absorption coefficient was chosen for the semi-infinite bottom layer, i.e. µ a,2 = µ a,1 , but a higher reduced scattering coefficient of µ s,2 = 2400 m −1 . This structure should thus correspond to a layering of relatively low-scattering non-coagulated tissue on top of relatively highly scattering coagulated tissue. It should therefore be suitable to characterize the reflectance of partly coagulated tissue. This will be investigated in detail in the next section. As can be observed from the graph of the amplitude response in figure 2, and in more detail for the three modulation frequencies of f = 100 MHz, 1 GHz and 2 GHz in figure 3, the amplitude, i.e. the reflectivity, decreases as the highly scattering layer is shifted away from the surface, and flattens out at a depth approximately 30% higher than the penetration depth of the upper layer. For a low modulation frequency of 100 MHz, the penetration depth of the upper layer is δ AC,1 = 6.8 mm; the reflectivity thus stabilizes once the bottom layer is more than 9 mm below the surface. At 1 GHz, the penetration depth is reduced to δ AC,1 = 4.9 mm, and the amplitude of the reflection coefficient changes only slightly once the upper layer is thicker than 7 mm. The same behaviour is noticed at a modulation frequency of 2 GHz.   The lower part of figure 2 shows the phaseshift γ of the reflected light relative to the incident light. The phase increases strongly as the highly scattering bottom layer starts moving away from the surface; the maximum phaseshift is reached at an upper layer thickness of 3 mm. This maximum is a clear sign for the interference of the partial waves in the two layers. Once the upper layer has become thicker than 3 mm, the phase drops off and tends asymptotically towards a constant value. The frequency response of the phase is characterized by two different behaviours: below the absorption relaxation frequency of f a = (2πτ a ) −1 = 477 MHz, the phase shift increases almost linearly with the frequency; at higher frequencies, the shift grows more slowly.
The corresponding results for the reversed sequence of the two layers, i.e. the highly scattering layer with µ a,1 = 10 m −1 and µ s,1 = 2400 m −1 on top of a low-scattering layer with µ a,2 = 10 m −1 and µ s,2 = 700 m −1 , is given in figure 4. As the upper layer grows in thickness, it again shows that the reflectivity gets close to a constant value more quickly for higher than for lower modulation frequencies. This corresponds to the behaviour of the penetration depths of the upper layer: 3.7 mm at 100 MHz, 2.6 mm at 1 GHz and 2 mm at 2 GHz. Contrary to the results presented in figure 2, the phase is found to drop sharply as the upper layer starts to grow, reaching a minimum at a thickness of approximately 1 mm. As the thickness continues to increase, the phase grows and quickly tends to a constant value. Again it is observed that for modulation frequencies below the absorption relaxation frequency of f a = 477 MHz, the phaseshift grows more strongly with increasing modulation frequency than at frequencies above f a .

Experiment
In order to verify the simulation results experimentally, a simple reflectance experiment was carried out on tissue slabs consisting of two layers of distinctively different optical properties. The layering could be achieved by placing a slab of chicken breast muscle on a hot plate. Thus, a coagulated layer of gradually growing thickness was covered by a top layer of diminishing thickness, consisting of tissue not yet denatured by the temperature. This order, with a lowscattering layer on top of a highly scattering layer, corresponds to an inverse layering of burnt tissue and has the experimental advantage that the amplitude and the phase of the reflection coefficient change less rapidly compared with the inverse order, as can be seen by comparing figure 2 with figure 4. The power of the hot plate was adjusted such that it gave a slight boiling at the interface to the chicken muscle. The heat propagated into the muscle by thermal conduction. The temperature rise T at distance x from the interface hot plate/meat slab and time t after switching on the heat approximately follows the law (Carslaw and Jaeger 1959) The temperature of the heater/tissue interface was T heater = 100 • C and the room temperature T room = 25 • C. The thermal diffusivity of the muscle tissue is given by χ = 1.5 × 10 −7 m 2 s −1 (Duck 1990). The temperature distribution versus distance from the interface and elapsed time is shown in figure 5. The horizontal plane at 65 • C corresponds to the protein denaturation temperature, i.e. describes when coagulation occurs. The intersection line between the two graphs, also shown as a projection into the x-t plane, marks the position of the coagulation zone depending on time.
The correctness of this model was validated by heating slabs of tissue for certain amounts of time and determining the position of the coagulation front by cutting through the slab. This inspection revealed that the coagulated zone was very well defined as long as it was thinner than 4 mm; at greater thicknesses, the transition to the uncoagulated part of the tissue became more diffuse. The time required for the coagulation zone to propagate 10 mm into the slab was measured to be about 15 min, in good agreement with the theoretically predicted values.
To carry out the reflectance measurement, the slab of chicken breast muscle, positioned on top of the hot plate, was irradiated by harmonically modulated laser light of 674 nm wavelength; see figure 6 for a sketch of the experimental set-up. A diode laser was modulated at frequencies from 100 MHz to 2 GHz by the frequency generator of a network analyser (impedance meter). The light was coupled into a 100 µm core graded-index fibre and guided to a combination of a collimating lens and a microscope lens, which produced an irradiated area of 5 cm on the sample when the lens combination was kept at a distance of approximately 20 cm. The reflected light was collected with a 600 µm core step-index fibre. This fibre was kept slightly askew relative to the source fibre to avoid picking up specularly reflected light from the surface of the sample, and its tip was positioned such that the field picked up was well within the irradiated spot. The signal was detected with an avalanche photodiode with a roll-off frequency of 1 GHz. The amplitude and phaseshift relative to the phase of the incident light could then be recorded by the network analyser. the time stamp taken at every measurement by inverse application of equation (13), assuming a coagulation temperature of 65 • C. The amplitude was normalized to coincide with the simulated values for a zero-thickness top layer. For the relative phaseshift, the measurement of the fully coagulated slab had to be used as a reference at every frequency, as the frequency response of the system was undetermined. When this normalization is disregarded, the measured phaseshift follows the prediction of the simulation (cf figure 2). With increasing thickness of the non-coagulated layer, the phaseshift increases rapidly, to find its maximum for an upper layer of 3 mm thickness. Thereafter, it sags slightly and settles at a constant value, which is higher than for the completely coagulated slab. The measured amplitude shows, in accordance with the simulated values in figure 2, a slight decrease as the modulation frequency increases. In contrast to the simulation, the measured values drop faster as the coagulation zone moves away from the surface.

Discussion and conclusions
The diffuse reflection coefficient of a turbid tissue is given by the flux at the upper surface of the slab, which is, according to the boundary condition (A.13), proportional to the fluence rate in this region. the fluence rate is influenced by different partial waves: the diffuse wave in the upper layer of a layered tissue interferes with a wave propagating in the bottom layer at a different speed, due to the difference in optical properties. Moreover, the diffuse light waves are subject to internal reflection at the tissue/air interface. The interference of these partial waves leads to marked maxima and minima in the phase of the diffusely reflected wave, depending on the optical properties of the layers. Moreover, the dependence of the penetration depth on the modulation frequency allows us to manipulate the probing depth of the diffuse wave. Thus, the technique investigated in this study permits the gathering of information about the structure of layered tissue, and should therefore be applicable to the investigation of subcutaneous lesions and the necrotization depth of burns.
The mathematical model proposed in this study has shown to be in comparatively good, in case of the phaseshift almost quantitative, agreement with the experimental results. As soon as more experimental data become available, further improvements can be expected with respect to the model. Moreover, in vivo measurements in a clinical environment will be important for an assessment of the applicability of the method.
where q describes the source distribution. For the case of a harmonically varying field, ϕ ∼ exp(−iωt), equation (A.6) can be rewritten into a second-order equation of Helmholtz type where the effective complex penetration depth is given by The source distribution in a two-layered one-dimensional model as sketched in figure 1, i.e. a layer of thickness d on top of a semi-infinite medium, is expressed by an approximate source function of the following form: where P 0 is the part of the incident light which is transmitted through the surface. The subscripts 1 and 2 refer to the upper and lower layers respectively. The phaseshift of the collimated source beam at a depth x, relative to the surface, is given by the factor exp (−iωx/c). This phaseshift is taken into account by introducing complex transport coefficients µ ω,tr,1 = µ tr,1 (1 + iωτ tr,1 ) 0 < x d µ ω,tr,2 = µ tr,2 (1 + iωτ tr,2 ) d <x<∞. where the constants A 1 . . . A 3 are determined by application of the boundary conditions at the interface between the two layers and the interface tissue/air. The condition at the boundary separating the two layers can be expressed by the continuity in fluence rate and diffuse photon flux ϕ 1 | x=d − = ϕ 2 | x=d + j 1 | x=d − = j 2 | x=d + . (A.12) This boundary condition is, however, invalid at the tissue/air interface, as this is an interface between a scattering and a non-scattering medium. Here, the partial-current boundary condition is applicable, which relates the diffuse photon flux out of the medium to the fluence rate right below the interface (Haskell et al 1994). It can be expressed in the form where R eff is the effective reflection coefficient. The value of R eff can be calculated by integrating the Fresnel reflection coefficient for unpolarized light over all angles of incidence, and is found to be 0.431 for an index of refraction of n = 1.33, and 0.493 for n = 1.40.
The diffuse reflection coefficient is given by the ratio of the diffusely reflected light, i.e. the light flux emerging out of the tissue, to the incident flux P 0 γ = − j| x=0 P 0 . (A.14) When the incident light is time varying, the reflection coefficient takes complex values. For steady-state light it is real. For the two-layer model, the reflection coefficient can be expressed by an explicit formula, as follows from equations (A.3), (A.11) and (A.14):