Abstract
This paper presents a review of the phenomena regarding light-tissue interactions, especially absorption and scattering. The most important mathematical approaches for modeling the light transport in tissues and their domain of application: “first-order scattering,” “Kubelka-Munk theory,” “diffusion approximation,” “Monte Carlo simulation,” “inverse adding-doubling” and “finite element method” are briefly described.
1. Introduction
When tissues are exposed to light reflection, refraction, absorption, or scattering can occur, which lead to energy losses in the incident beam.
Refraction is not significant in biomedical applications, except for laser irradiation of transparent media, such as cornea tissue; in opaque media the most important phenomena are scattering and absorption, depending on the material type of the tissue and the incident wavelength. Knowledge of absorbing and scattering properties of the tissues is needed for predicting success of laser surgery treatment.
Direct measurement methods simply use the Beer attenuation law, but they need corrections when surface reflections occur due to the mismatched refractive indexes.
Indirect techniques use theoretical models for the scattering phenomena; the indirect noniterative methods need simple equations to connect optical properties to the measured quantities, while the indirect iterative methods can develop sophisticated models in which the optical properties are iterated until the computed reflection and transmission match the measured values.
2. Basic Phenomena Regarding Light and Tissues
Reflection means the electromagnetic waves return from surfaces upon they are incident, generally being boundary surfaces between two materials of different refractive indexes, such as air and tissue. The simple law of reflection states that the reflection angle equals the incidence angle, while the surface is supposed to be smooth, having small irregularities compared to the radiation wavelength. The real tissues do not act like optical mirrors; the roughness of the reflecting surface leads to multiple beam reflections (diffuse reflection).
Refraction means a displacement of the transmitted beam through the surface that separates two media with different refractive indexes, and it originates from the change of the speed of light passing through the surface.
Refraction usually occurs together with reflection; the reflectivity of a surface is a measure of the amount of reflected radiation, and it is the ratio of the reflected and incident electric field amplitudes.
The reflectance is the ratio of the reflected and incident intensities; thus, it is equal to the square of the reflectivity. At normal incidence the reflectance of a surface separating air and water is about 2% [1, 2], but in several cases the reflectance can be important, with values that cannot be neglected; in these cases proper eye protection when using lasers is needed.
Even one might expect that the intensities of the reflected and of the refracted beams would be complementary such as their sum would be equal to the incident intensity, which is not exact; the intensity is the ratio between the power and the area unit, while the cross-sections of the incident, reflected and refracted beams do not have the same value, except for the normal incidence. But the total energy in these beams is conserved.
Absorption means that the incident beam intensity decreases when passing through the tissue, because the beam energy is partly converted into heat motion or into vibration of the absorbent medium molecules.
The laws of absorption describe the effect of the tissue thickness (Lambert law): or the effect of the concentration (Beer law): where is the beam intensity at the distance , is the intensity of the incident beam, is the absorption coefficient of the medium, is the concentration of the absorbing medium, and depends on other parameters of the substance than concentration. The inverse of is named the length of absorption, and it measures the distance at which the incident intensity drops by .
Absorption in biological tissues is mainly determined by the water molecules, especially in the IR region of the spectrum, and by protein and pigments in the UV and visible range; melanin is the basic pigment of the skin and hemoglobin is component of vascularized tissues. Most biomolecules have complex absorption band structures in the 400–600 nm range, but the spectral range of 600–1200 nm (the therapeutic window) is free of absorption phenomena; neither water nor macromolecules absorb near IR, so that the light penetrates biological tissues with little loss and enables treatments of profound structures. The skin is the highest absorber in the visible domain, while the cornea is totally transparent.
There is an almost perfect match between the absorption peaks of skin melanin and hemoglobin (Figure 1) and the green and yellow wavelengths at 531 nm and 568 nm of krypton ion lasers, meaning these lasers can be used for coagulating blood vessels. Sometimes the original absorption of the tissue is increased using special dyes and inks before laser treatments, so that the treatment efficiency would be higher and the neighborhood tissues are less affected.
Scattering can be elastic (when the incident photon energy has the same value as the scattered photon energy) or inelastic (when a fraction of the incident photon energy is converted into forced vibrations of the medium particles). Rayleigh scattering is of elastic type, where the scattering coefficient decreases with the fourth power of the wavelength.
Brillouin scattering is inelastic and is caused by the acoustic waves that induce inhomogeneities of the refractive index; it can be seen like a kind of optical Doppler Effect, as the scattered photon frequency is shifted up or down when the scattering particles move towards or away from the light source.
The decrease of the intensity during scattering is given by a similar law as absorption: In most tissues (turbid media) both absorption and scattering occur simultaneously; their entire attenuation coefficient is
2.1. Unscattered (Coherent) Transmission
Considering an unscattered beam, with no surface reflections, incident on a slab of tissue having the thickness d, the transmission is exponentially attenuated (Beer’s law); the unscattered (collimated) transmission is given by Thus, the total attenuation coefficient can be determined: In the presence of the mismatched surface reflections, corrections are required; for instance, when a tissue sample is placed between glass slides, the collimated beam is reflected at the air-slide, slide-tissue, tissue-slide, and slide-air surfaces. If the sample is only a few optical depths (8) thick, one must consider multiple internal reflections. Thus, a net reflection coefficient is given [3] by where and are the Fresnel reflections at the air-glass, respectively, glass-tissue interfaces; then the measured transmission is solved to obtain , in order to calculate .
The optical depth is defined by where is the element of the optical path and is the total length of the optical path. If the attenuation coefficient is constant, we meet the case of homogeneous attenuation; thus, .
In the literature [1, 3], we also find reduced scattering and attenuation coefficients: which include the particular situation when (just forward scattering) when the intensity is not attenuated. The linear transport coefficient describes the inverse of the mean free path between two interaction events in a strong scattering medium.
No theory completely explains why the photons are preferably scattered in the forward direction, nor the wavelength dependencies of the phenomena. A more convenient approach is to work with the photon probability function (or phase function) to be scattered by an angle , which can be fitted to experimental measurements. According to this dependence, the scattering can be isotropic, which means does not depend on , or else anisotropic. One can evaluate the anisotropy of the scattering using the coefficient of anisotropy defined as the average cosine of the scattering angle : where is the solid angle element. Isotropic scattering will use , while means totally forward scattering and means totally backward scattering. For most biological tissues we can assume that values are in the range of 0.7–0.99 [1] because the most frequent scattering angles range between 8° and 45°.
The probability (phase) function is normalized: The best fit to the experiments (Figure 2) was given by the Henyey-Greenstein phase function: This is equivalent with the series: where are Legendre polynomials.
The absorption coefficients are generally obtained subtracting transmitted, reflected, and scattered intensities from the incident intensity; experimental methods determine either the total attenuation coefficient or both absorption and scattering coefficients. Rotating the detector, one can obtain the angular dependence of the scattered intensity, thus having the anisotropy coefficient .
3. General Physical Model for Light Propagation in Tissue
Mathematical description of absorption and scattering can be made in two ways. The most fundamental approach is the analytical theory based on Maxwell’s equations, but their complexity and the inhomogeneities of biological tissues limit the possibility to obtain the exact analytical solutions, thus limiting the applicability of the theory.
The second approach is the photon transport theory, which deals with photon beams passing through absorbing and scattering media, without considering Maxwell’s equations; it was extensively used for laser-tissue interactions where its predictions were satisfactory in many cases, though it is a less strict theory compared to analytical theories.
Radiance J is the power flux density in a given direction within the solid angle unit (W·cm−2·sr−1); it decreases due to scattering and absorption but increases by the light scattered from directions into direction : where is the phase function of a scattered photon from to direction, is the attenuation coefficient which includes scattering and absorption coefficients, and denotes the infinitesimal path length. The differential equation (15) for radiance is called the radiative transport equation. For symmetric scattering about the optical axis, , where is the scattering angle.
Intensity is the measurable optical property; it is obtained by integrating radiance over the solid angle: Its value can be measured within the experiments, so one can express radiance by where means a solid angle delta function oriented in the direction.
When laser propagates inside a turbid medium, the radiance can be separated in two terms depicting a coherent (unscattered) component and a diffuse (scattering) component: The coherent radiance can be calculated from which has the solution where is the incident intensity and the dimensionless parameter is given by (5). The unscattered term contains all the light that did not interact with the tissue and is characterized by exponential decaying.
Evaluation of the second term, the diffuse radiance, is an important problem as it contains all light that has been scattered at least once; one cannot exactly determine the path of all the scattered photons so adequate approximations or statistical approach have to be made depending on whether absorption or scattering is dominant. Complexity of the approach is related to its accuracy and to the computing time it needs.
4. Mathematical Methods
Mathematical methods are based upon assumptions regarding the incident light sources and boundary conditions. They are referred to as “first-order scattering,” “Kubelka-Munk theory,” “diffusion approximation,” “Monte Carlo simulation,” and “inverse adding-doubling.”
4.1. First-Order Scattering
Scatering means that multiple scattering is not considered; in some cases the diffuse radiance is much smaller than the coherent radiance, and one can assume it can be neglected: and the intensity is given by Lambert law: Although it is a simple solution, it is often inapplicable for the therapeutic window (600–1200 nm).
4.2. Kubelka-Munk Theory
Unlike the previous assumption, Kubelka and Munk supposed the radiance to be diffuse () and defined two coefficients for the absorption and for scattering of diffuse radiation, different from and that referred only to coherent radiation.
Considering diffuse radiance through a one-dimensional isotropic slab with no reflection at the boundaries, this approach is equivalent to the diffusion model with two-phase functions, peaked forward and backward [3].
A flux in the incident direction and a backscattered flux are shown in Figure 3; two differential equations can be written, each states that in both directions the radiance encounters two losses (scattering and absorption) and a gain from opposite direction scattered photons: where is the mean direction of the incident radiation. Their general solutions are with Considering average values of scattered and coherent path lengths [1] and since diffuse scattering implies that does not depend on the scattering angle, one can find the connection between the absorption and scattering coefficients:
The Kubelka-Munk expressions for reflection and transmission of diffuse radiance on a slab of thickness are [3] so that the scattering and absorption coefficients can be expressed in terms of the measured values of and : where the parameters and can be found using This theory is a simple model to measure the optical properties of the tissues, but it is restricted by the assumption that the incident light is already diffuse, the isotropic scattering and matched boundary refractive indexes, which are atypical for laser-tissue interaction. It can be extended considering more fluxes, but an important disadvantage is the extended computing time.
4.3. Diffusion Approximation
When the scattering phenomena dominate absorption, the diffuse term in (18) can be expanded in a series: where is the diffuse intensity and the vector flux is given by The diffuse intensity satisfies where denotes the diffusion parameter () that is an approximation of the measured effective attenuation coefficient of diffuse light: where denotes the effective diffusion length and is the term for the source of the scattered photons. It was shown [1] that where is the incident flux and is the optical depth given by (9).
Finally, the diffusion approximation states that with . Different sets of values for , , and provide similar radiances in diffusion approximation calculus.
4.4. Monte Carlo Simulations
The Monte Carlo method essentially runs a computer simulation based upon a numerical approach to the transport equation (15). The statistical approach implies the simulation of a number of photons random walk; the statistical accuracy of the results is proportional to , so that a valuable approximation has to take into account a large number of photons. This method has become a powerful tool for many disciplines and it requires large computers or networks.
The main idea of Monte Carlo method for the absorption and scattering phenomena is to follow the optical path of a photon through a turbid medium. The distance between two collisions is chosen from a logarithmic distribution using a random computer generated number.
Absorption is depicted as a decrease of a weight attributed to the photon during propagation; scattering is provided by choosing a new direction of propagation, according to a given phase function and another random number. The whole procedure continues until the photon weight reaches a minimum cut-off value, or the photon escapes from the considered region.
Monte Carlo simulation (Figure 4) needs five steps [1].(1)Photons are generated at a surface of the considered region, so that their distribution can be fitted to a given light source (i.e., Gaussian beam).(2)Pathway generation: the distance to the first collision is computed supposing that absorbing and scattering particles are randomly distributed: a random number is generated, so the distance where is the particle density and is their scattering cross-section [1]; thus, represents the mean free path and a scattering point is obtained. Then a second random number is generated to determine the scattering angle according to the phase function; then the third random number is generated to get the azimuth angle: (3)Absorption makes the photon weight decrease by , being the absorption coefficient.(4)Elimination of the photon when its attributed weight reaches a certain value, then a new photon is launched and proceeds with step 1.(5)Detection: after repeating steps 1–4 for a sufficient number of photons, the computer has stored a map of pathways, so one can make statistical predictions upon the fraction of the incident photons being absorbed by the medium and the spatial and angular distribution of the transmitted photons.
As the simulation accuracy increases with the larger numbers of photons, the necessity of extending calculations makes them time consuming. If earlier results are stored in the computer memory and used again if needed with the same phase function, computing time will be saved.
4.5. Inverse Adding-Doubling Method
The name of the method comes from reversing the usual process when calculating reflectance and transmittance from the optical properties; this method assumes that reflection and transmission of incident light at a certain angle are known. If we need the same properties for a layer twice as thick, we divide it into two equal slabs, and then add the reflection and transmission contributions of either slab. These properties can be obtained for an arbitrary slab of tissue starting with a thin, known slab, and doubling it until the necessary thickness is achieved.
The adding method can be extended to simulate layered tissues with different optical properties.
4.6. The Finite Element Method: Boundary and Source Conditions
The model to calculate light transport in strongly scattering materials was needed as a component in image reconstruction schemes that localize the optical properties of the tissues from boundary measurements by solving the inverse problem. Diffusion approximation to the radiative transfer equation (15) is an adequate model for scatter-dominant materials for it assumes that scattering dominates absorption (generally true for biological tissues) and the anisotropic propagation is weak, that is not the case near sources and boundaries.
When the only sources of light are monochromatic lasers, a frequency-independent model of light transport is suitable; multi-wavelength systems will be applied sequentially, so that the change of frequency is expressed by the change in the optical properties.
The finite element method for the propagation of light in scattering media aims [5, 6] to find the photon density and the radiance in a strong scattering domain , together with the outward current (existence) through the boundary , using the diffusion approximation model to the radiative transport equation (15). The existence through boundary at the point is defined as where is the normal to at .
The diffusion equation (32) can be written as where is the photon density at and denotes the diffusion coefficient (34). The solution we seek is a continuous, linear approximation of . The considered domain Ω is partitioned into nonoverlapping elements , in each element is assumed linear. Nodes are attached to the element vertices and at each point within the element is the linear interpolation of nodal values where are linear nodal shape functions [5] with support over all elements which have the node in the position as a vertex, and . Equation (39) becomes The solution of (32) requires appropriate boundary conditions. The Dirichlet condition (DBC) means that the medium around is a perfect absorber; photons are absorbed when crossing , so that outside the domain the photon density equals zero. A more realistic boundary condition would be which is named [5] a Robin boundary condition (RBC) that constrains a linear combination of the photon densities and the current at . It can be modified to incorporate a mismatch of the refractive indexes within and in the surrounding medium; we assume so that (41) becomes where is the parameter governing internal reflection at the boundary that can be fitted from experimental curves. A modified Robin condition [5] can be with , or if we use a different approach [6] to derive from Fresnel’s laws: with the critical angle and (i.e., for we get ).
While RBC is a more accurate physical situation than the DBC, the mathematical approach is simple with DBC. To get a compromise, we can introduce an extrapolated boundary at a certain distance from the physical boundary at which DBC apply. For an arbitrary two-dimensional domain , the extrapolated boundary condition is obtained simply by adding a border of thickness around . Different expressions [5] of extrapolating can be applied to arbitrary geometries if the radius of the boundary curvature is small compared to the mean free path, both for matched and mismatched refractive index.
There are two possibilities to model the light sources incident at a point on the boundary: collimated or diffuse sources.(1)The diffusion equation cannot describe correctly collimated sources, so we can represent a collimated pencil beam by an isotropic source at the depth that is accurate at distances larger than the mean free path from the source but breaks down close to the source. The implementation of the collimated source needs a delta-shaped term in (39) or other models have been using the analogy to nuclear engineering with cylinder sources and exponentially decaying.(2)Diffuse sources on the surface can be regarded as an inward directed diffuse photon current, distributed over the illuminated boundary segment . The Dirichlet condition (40) is written for diffuse sources where , is the source current strength, is a weighting function, and is the outward normal to at .
The inclusion of the source as a photon current through boundaries modifies the Robin condition (43) along to Experimental data [5, 6] lead to the conclusion that the Robin model and the extrapolated boundary models are equivalent, but the Dirichlet model shows significant lack of accuracy that should prevent its use in most applications.
Figure 5 shows intensity distributions calculated [1] with either method inside a turbid medium, assuming isotropic scattering; denotes the ratio of absorption and scattering coefficients (optical albedo): The values and mean that scattering is predominant. The ordinate shows diffuse intensity in units of incident intensity.
5. Conclusions
Different methods for solving the transport equation were discussed; the most important are the Kubelka-Munk theory, the diffusion approximation, and Monte Carlo simulations.
The Kubelka-Munk theory deals only with diffuse radiation and it is limited to the cases where scattering dominates absorption; an important disadvantage is the one-dimensional geometry.
The diffusion approximation is not restricted to diffuse radiation, but it also applies to predominant scattering phenomena, which is a powerful tool.
Monte Carlo simulations provide most accurate solutions, since many input parameters may be taken into account in specially developed computer programs; the method allows two-dimensional and three-dimensional evaluations although they require a long computing time.
The finite element method was compared [5] to Monte Carlo corresponding methods and it was found to be in good agreement with Robin boundary conditions; the diffusion approximation, although not strictly valid near boundaries, is still able to produce correct data when using appropriate boundary conditions, which matches best a given experimental situation. It can be used when the Monte Carlo computing speed is significantly reduced.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.