This report is a tutorial introduction to diffuse light transport in biological tissues. Section 1 presents the basics of diffuse light transport, showing the simple equations for time-resolved, steady-state, and modulated light transport. The perturbation method for handling slight heterogeneities in optical properties is introduced. The treatment of an air/tissue surface boundary condition is considered. Section 2 describes numerical methods for simulating light transport in complex tissues. The goal of this work is to provide the novice in biomedical optics with an introduction to diffuse light transport, and the underpinnings of how basic approaches to solving light transport problems can be solved.
Analytic Modeling of Diffuse Light Transport
Optical Properties and Transport Parameters
This tutorial follows the notation of Welch and van Gemert.1 The tissue optical properties are shown in Table 1 . These properties yield the transport properties used in diffusion theory, as given in Table 2 . In discussing light transport, the parameters given in Table 3 are used. The units of power and light transport differ for a point source, a line source, and a planar source, but the fluence rate always has the same units. For a point source, has units of and has units of . For a line source, has units of W/cm and has units of . For a planar source, has units of and is dimensionless.
Tissue optical properties.
|Anisotropy of scattering||dimensionless|
Transport properties used in diffusion theory.
|Reduced scattering coefficient|
|Transport mean free path||cm|
|Optical penetration depth||cm|
Parameters for light transport.
|Radiant power||W, W/cm,|
|Transport||, , dimensionless|
|Speed of light in tissue||cm/s|
|Concentration of radiant energy|
Diffuse Versus Nondiffuse Regimes
The movement of light can be discussed in several ways. For example, 1. light moves as an electromagnetic wave in a medium (used in interferometry), 2. light moves as ballistic photons, each with a direction of travel that can be redirected by scattering (used in Monte Carlo simulations), or 3. light can be treated as a concentration of optical energy that diffuses down a concentration gradient (used in diffusion theory simulations). This work is devoted to discussing the third example, the movement of optical energy by diffusion down concentration gradients.
Before proceeding, it is important to understand when diffusion theory is not appropriate. Diffusion assumes that the quantity that is diffusing, in our case optical or radiant energy, does not have a preferential direction of travel. Hence, the net movement of that quantity is by diffusion down a concentration gradient, following Fick’s first law of diffusion (discussed later). When photons are delivered as a collimated beam into a medium, they definitely have a direction of movement. As the photons are scattered by interaction with a tissue, they lose their directionality and hence become eligible for diffusion.
Monte Carlo simulations provide a method for specifying the movement of ballistic photons.2 Figure 1 compares the Monte Carlo and diffusion theory descriptions of the distribution of light in an infinite medium where light is delivered as a collimated beam by an optical fiber submerged in the light-scattering medium. The figure illustrates that in a local region near the source, diffusion theory fails to accurately describe the light distribution. However, distant from the source, diffusion theory is quite good, as if the collimated photons had been thrown forward by the optical fiber and the diffusion process had emanated from this location, a distance in front of the fiber. This distance is called the transport mean free path, or (sometimes called ). The figure illustrates diffusion theory and the Monte Carlo approach agreement at distances exceeding one from the tip of the fiber and exceeding one from the apparent source of diffusion located one beyond the fiber tip.
The idea of a hybrid Monte Carlo/diffusion theory has been introduced, in which a Monte Carlo simulation launches the photons and diffusion theory propagates them after they have multiply scattered.3 A Monte Carlo simulation launches photons into a scattering medium and propagates the photons according to probability density functions for the step size between photon/tissue interaction sites and the angles of deflection. As the photon propagates initially, the fluence rate in response to a 1-W incident power rate is recorded in an array, [ per W], which is the transport When a photon reaches a depth that exceeds , the Monte Carlo simulation of that photon movement is halted and the photon power is deposited into the local voxel of an array that accumulates the photons as a distributed source [W per voxel]. After launching many photons, the Monte Carlo simulation yields both a and a distributed source . A point spread function or Green’s function, , based on diffusion theory is convolved against to yield the fluence rate, , throughout the medium due to diffusion. The final answer is .
If one ignores the region near a source, then diffusion theory is a very useful and reasonably accurate method for describing the distribution of light in a light-scattering light-absorbing medium. If the observation point is more than one or two transport mean free paths from a source, diffusion theory is usually accurate to within a few percent. Typically in the visible or near-infrared wavelength range, a distance of one corresponds to less than 1 or . So diffusion theory is quite appropriate and useful when discussing light transport over several millimeters or more.
Diffusion theory also assumes that photons are able to participate in a random walk. Hence, photons should be able to undergo several scattering events before being absorbed. The commonly cited criterion, which is only a rough rule of thumb, is that the ratio should exceed 10. If there are sufficient scattering events before an absorption event occurs, then the average photon propagation after many scattering events becomes dependent on rather than on the specific values of and . The total diffuse reflectance from a tissue that has is 0.252 (based on Monte Carlo simulations), where the mismatched tissue/air surface boundary has . A study of reflectance versus where and were varied but was conserved, indicated that reflectance becomes independent of when reflectance exceeds .4 So the reflectance of a tissue should exceed at least 25% and preferably 40%, or should exceed 10 and preferably 20, if one wishes to use diffusion theory reliably. Fortunately, this is usually the case for tissues when using visible to near-infrared wavelengths of light. In tissues with lower reflectance or ratio , diffusion theory can still be useful but caution as to the accuracy of diffusion theory is advised.
In summary, diffusion theory is very quick and useful. It is not good near sources. This caution extends to locations near boundaries between regions of differing optical properties, especially the air/tissue surface, and to locations near strongly absorbing objects. Diffusion theory is not so good in tissues with strong absorption, where . For much work in the visible to near-infrared wavelength range, diffusion theory is quite good, and should be understood by every student of biomedical optics.
Fluence rate and concentration of optical energy
The fluence rate is proportional to the concentration of optical energy :is the speed of light in the medium, [cm/s], , and is the refractive index of the medium. For example, consider a uniform beam of irradiance passing through a cube of water (see Fig. 2 ). If the fluence rate in the nonscattering water equals at wavelength, then the corresponding concentration of photons is : , where is wavelength, is Planck’s constant, and is the in vacuo speed of light, indicates the number of photons per J energy. The factor 1000 indicates the conversion , and is Avagadro’s number for the molecules per mole.
Flux and concentration gradients
Fluence rate can be regarded as a measure of concentration, and the diffusion of light as the movement of photons down concentration gradients. This gradient-driven movement of light is described by Fick’s first law of diffusion, in which the flux down a concentration gradient isis the diffusivity equal to for light, with the diffusion length . In general, can refer to the diffusion of heat, molecules, etc., or radiant energy. For light, and , and the product . Therefore, the speed of light is canceled in Eq. 3.
The description of the flux of optical energy down a gradient of fluence rate along the dimension follows Fick’s law, which relates the local gradient to the local flux [Eq. 3]. The description was originally developed to describe neutron scattering in nuclear reactor cores, where the role of absorption was neglected. Consider a small incremental window of cross sectional area at , , oriented such that the normal to the window is aligned with the axis. The flux of light through the window is due to the scatter of light from two small hemispherical regions on either side of the window, which scatter light through the window from both directions. For this example, isotropic scatterers are assumed, described by the reduced scattering coefficient . The fluence rate upstream from the window along the gradient encounters the scatterers within an incremental volume and a power is scattered [W]. A fraction of this scattered power is directed toward the window, which is a distance from the volume , where is the angle between the normal of the window and the direction of scatter toward the window. The fraction of this power aimed at the window that reaches the window is , since scattering and absorption attenuate the light. This attenuation limits the radius of the hemispheres from which the window accepts scattered light. The balance between these two fluxes through the window, the downstream flux minus the upstream flux, yields the net flux through the window. Normalizing the net flux by yields the net flux density passing through the incremental window:is the incremental volume expressed as an annular ring. The factor causes the integrand to be positive when (upstream of the gradient) and to be negative when (downstream of the gradient). This expression yields a value for the flux density given a particular gradient . The diffusion length is thereby specified as , using Fick’s first law.
Fick’s first law holds if the gradient of fluence rate is linear within the two hemispherical regions of integration, i.e., , where is the position along the axis perpendicular to the window and is a constant. This situation is called diffusion theory, and is sometimes called the approximation. If, however, there is too much curvature in the gradient of , then the near the window must be described by higher order terms, . An example of this is the so-called approximation. Such curvature occurs when the point of interest is close to a source of light or a sink for light (an absorber), or near a boundary of mismatched optical properties. Hence, diffusion theory does not work well near discontinuities, and the approximation is a more appropriate method to describe the transport. However, if the source and point of observation are separated by more than a , the approximation is usually sufficient. The presence of a boundary is discussed later (see Boundaries in see Sec. 2.5).
The use of Fick’s law has two important consequences that should always be kept in mind. First, there is continuity of the fluence rate across any boundary, even if the tissues on either side of the boundary have different optical properties. Second, the flux incident on a boundary equals the flux leaving the boundary, in other words photons do not accumulate on a boundary surface. If two tissues with differing absorption and scattering properties are in contact and light diffuses through their common boundary, the fluence rate across the boundary will be continuous, although the rate of energy deposition can be discontinuous.
Now that photon movement is understood as a diffusion process, consider the standard time-resolved diffusion equation, which is applicable to anything that can diffuse but is now applied to the diffusion of radiant energy. The generic solution to the diffusion equation, for diffusion from a point source, or delta function source, in both space and time is:is the point source of whatever substance or energy [units] that will diffuse, for example heat, molecules, photons, etc., which is deposited at the origin at time zero . The second factor has units of and is referred to as Green’s function. Hence, the concentration at a distance from the source at a time after the initial deposition of is expressed as .
To describe the diffusion of fluence rate after deposition of a point source of radiant energy, , the speed of light is introduced, since . Also, the diffusivity is replaced by . Finally, the role of absorption is added by including the term , where is the path length of individual photons at time after the impulse, and is the absorption coefficient. The result is:3 gives an example of time-resolved diffusion of light from an impulse.
Integration of the prior time-resolved in Eq. 6 over all time yields the radiant exposure , which does not depend on time. The direct integration of Eq. 6 using rules of integration is not immediately obvious. Intuition indicates that should fall exponentially with distance from a source. A numerical integration of over all time at each position (not shown here) confirms that at large , the radiant exposure falls as . The constant is given the symbol [cm] and is called the optical penetration depth. The factor is called the effective attenuation coefficient . The integration becomes easier by considering conservation of energy. The proper form of is such that the integral of energy deposition over an infinite volume equals the total energy deposited . The energy deposition is integrated over the infinite spatial domain:is expressed as an incremental shell volume . The denominator called factor needs to be specified. The factor needs the term to yield the final value . The factor needs to cancel the in the numerator. The factor needs to cancel the in the numerator. The factor needs to reduce the in the numerator to simply , so that can be integrated by parts. A factor is generated by this integration, so the factor needs a to cancel this term. In summary, the factor must equal to allow conservation of energy in Eq. 7. Therefore, the expression for the radiant exposure in response to an impulse of energy at and is is identical to (derivation not shown here), and the expression using is also shown in Eq. 8 and is commonly used.
This is the impulse response of radiant exposure to an impulse of energy deposition . If one has a repetitive sequence of impulses at a frequency , the time averaged power is [W]. In the limit of and such that is maintained at a constant average power, the superposition of the time-resolved responses to each impulse blurs into a steady-state fluence rate proportional to , with replaced by :. Figure 4 gives an example of steady-state diffusion of light.
The transport factor equals and characterizes the light transport in a tissue independent of the strength of the source. The following equations list the transport for 1-D planar diffusion from a planar source , for 2-D cylindrical diffusion from a line source [W/cm], and for 3-D spherical diffusion from a point source, [W], derived as outlined before from conservation of energy:
If the power of a point source is modulated sinusoidally, then a population of photons will be periodically injected into the medium. The time-resolved source is:is the modulation of the source, [dimensionless], and is the angular frequency of modulation [radians/s], , where f is the frequency in [Hz].
These injected photons will diffuse down the gradient from the source toward some point of observation at a detector. The diffusion theory expression for transport from such a time-resolved modulated source (a point source within a medium with no nearby boundaries) isis the steady-state transport: and are defined:5 The consequence of these expressions is that the fluence rate seen at the detector has a steady-state component , plus a modulated component characterized by a modulation [dimensionless] and a phase [radians]: is swept over a broad range, then the values of and versus frequency constitute a frequency domain description of light transport.
For very low frequencies of modulation as , the factor approaches 1 and the factor approaches 0. Therefore, approaches and approaches 0. approaches , and approaches 0. So the transport becomes.
For very high frequencies of modulation as , the factor approaches , and the factor approaches . Since , the values and approach the same limiting value:is a little below or a little above , the behavior of and varies versus , each in a unique manner, and therefore and observed at a particular distance from the source encode the optical properties and . At very low , the phase is too low to reliably measure. At high , the becomes insensitive to . So the use of modulated light is typically in the intermediate region where is a little less than . For a tissue with and , the value is , and the frequency is . Typically, modulated diffusion measurements are made with frequencies in the range , which correspond to , respectively. Figure 5 shows and versus for a series of modulation frequencies. Figure 6 shows and versus modulation frequency for a series of positions.
Modeling Heterogeneities by the Perturbation Method
Diffusion theory can often provide a quick assessment of the influence of a perturbing object within a tissue. Such a perturbing object can be a region of tissue with different absorption and/or scattering properties, perhaps a tumor with an elevated density of blood vessels due to angiogenesis, a portwine stain lesion with abnormal vasculature, a fibrous growth with increased scattering, or a fluid-filled cyst with decreased scattering.
The perturbation method can use the Green’s function of diffusion theory to assess the influence of a perturbing object on a distant point of observation. The method is introduced using the steady-state diffusion equation, but the method can be developed using frequency-domain transport, as shown in Ref. 6, or time-resolved transport. The method need not be based on diffusion theory, any transport function will suffice.
Consider the example in Fig. 7 . If a power of light is launched isotropically at a position , the transport to a detector at position over a distance is. The background optical properties are and , yielding a background diffusion constant , and a background optical penetration depth .
Now consider the presence of a small absorbing object of volume at position with an extra absorption , such that its absorption is . The fluence rate at the object is , where is computed using a source-object distance . The extra energy deposition in the absorbing object, above the background absorption, is . The total extra power deposited in the object is [W], which has the units of power. Because the object is absorbing, the object behaves as a source of negative power :. The net fluence rate at the detector is . If , then is positive and the perturbation is positive.
This is the basic idea of perturbation theory. More details can be found in Ref. 6.
If one has many real sources and many virtual sources due to many perturbing objects, the process is treated as a matrix problem. A vector lists all sources. A vector lists only the virtual sources, and is a subset of . All the real sources are called , and the virtual sources due to the perturbing objects are called . A perturbation matrix is defined with columns to match and rows to match . The matrix consists of elements that scale each source in to yield a contribution to an updated value of :is calculated by . This first calculation of is called the Born approximation. The values of are used to update the virtual source values in :22, 23 is iterated, until the values of no longer change. This occurs rapidly for an absorbing object, usually less than nine iterations. For a scattering object (see later), convergence usually takes more iterations, perhaps 20 to 50, but is still rapid. For scattering objects, it is usually very helpful to let each update be a partial update, for example only a step toward the calculated update value, in other words:
If an object is rather large, then it is better to subdivide the object volume into smaller subvolumes, each acting as an independent virtual source. For a scattering object, virtual sources, some positive and some negative, are placed on the surface of the object to create an effective dipole that mimics the behavior of the scattering object. The following discussion outlines these procedures. Figure 8 illustrates the virtual sources described, Fig. 8a showing the absorbing subvolumes and Fig. 8b showing the surface scattering sources that yield a dipole with positive sources on the hemisphere facing upstream into the gradient of fluence rate toward the real source, and negative sources on the hemisphere facing downstream down the gradient of fluence rate away from the real source.
There are two kinds of perturbation that contribute to the source vector and the elements in the matrix. One is due to absorption and one is due to scattering. Consider a spherical object of radius that is both absorbing and scattering. The volume of the object can be subdivided into a set of subvolumes that fill the object’s volume, and each is called a “virtual absorbing volume source.” The vector is filled with these virtual absorbing volume sources, in this example expressed as small spheres, each of radius and subvolume . The radius is usually chosen to approximately equal , which sets the number of small spheres required to fill the object’s volume , but the total volume of all the small spheres must equal the volume of the object , so . The elements of the matrix for are specified for each ’th source by:denotes the effect of the ’th source on the ’th virtual source, and is the same as the introductory example presented earlier [see Eq. 20]. The second term, is a scattering-enhanced absorption due to the extra scattering that causes light to be trapped within the object, thereby extending its path length within the object.
The surface of the object is studded with point sources called “virtual scattering surface sources,” which mimic the scattering of the object due to its . The surface is divided into a set of incremental surfaces, each with a surface area . The sum of the incremental surfaces equals the total surface area of the object. The vector , where , is filled with these virtual scattering surface sources. The elements of the matrix for , where , is specified for each source by:is the transport from the ’th source to a position just outside the object at the ’th virtual source position, and is the transport from the ’th source to a position just inside the object at the ’th virtual source position. The distance is usually chosen to equal , such that these two sources inside and outside the object’s surface are sufficiently distant from each other to interact using diffusion theory. Positions 1 and 2 are aligned along the normal to the surface. Consequently, the factor is the flux crossing into the object normal to the surface at the incremental area. The strength of the virtual surface scattering source is proportional to this flux. If flux gradient is driving light into the object, as occurs on the upstream side of the object closest to the real source, then the value of is positive. If the flux gradient is driving light out of the object, as occurs on the downstream side of the object away from the real source, then the value of is negative. The positive and negative surface sources on the object create a dipole that radiates a dipolar perturbation oriented along the gradient of flux. Light scatters upstream toward the true light source, and a shadow is cast downstream behind the object.
In the matrix , a special case occurs when equals , and the self-perturbation by a virtual source must be evaluated. This term is here called . For the scattering surface sources, the value of is zero, because a point source cannot exert a gradient across itself. For the absorbing volume sources of radius , the value of is calculated:22, 23], a stable vector is achieved. The fluence rate at any observation point, for example at a detector, can now be calculated: is the transport from the ’th source to the point of observation. All the sources , both real and virtual, contribute to the observed fluence rate.
Figure 9 shows an example of a -diam. spherical object with adding absorber ( , , , ). The background light field is shown in Fig. 9a. The light source is located at . The perturbation itself, , is shown in Fig. 9b. The object casts a spherical zone of negative perturbation around the object, about effect ( per W of the source) within the object and effect in the region surrounding the object. Figure 9c shows the fractional perturbation, [dimensionless], illustrating that the perturbation is about within the object and cases a shadow of perturbation about to behind the object, downstream from the source. The object’s perturbation can significantly affect the light field downstream from the object, but it cannot overcome the strong light field upstream toward the source.
Figure 10 shows the same spherical object, this time with added scattering ( , , , ). This time, the object backscatters a positive perturbation upstream toward the source and casts a negative shadow downstream behind the object. Figure 10c shows that the fractional perturbation is a little different than the absorption only example [Fig. 9c], but essentially presents the same behavior of casting a shadow downstream behind the object away from the source. These examples illustrate that both increased absorption and increased scattering cast shadows downstream, but have marginal effect upstream. A detector is more sensitive to the presence of a perturbing object if the detector is placed in the shadow cast by the object.
In summary, the perturbation method is rather simple and quick, and provides a good estimate of the effect of a perturbing object on the field of light at some distant point of observation, such as a detector. The perturbation method has been used to interpret the gastric endoscopic images of reflectance to determine the size and depth of a subsurface blood vessel, to aid the decision about how to handle a gastric bleeding site,7 and to assess the ability to measure the oxygen saturation of a neonatal brain while still in the birth canal by delivery of light to and collection of escaping optical flux from the maternal abdomen, to aid the decision to deliver the neonate by Caesarean section.8
A word on boundaries is worthwhile, since most biophotonic applications involve delivery and/or collection of light through an air/tissue surface. The boundary condition is modeled using the method of images adapted to solve the relationship between fluence rate at the surface and flux escaping the surface. The air/tissue boundary is replaced by an infinite tissue with a source of negative light located above (or outside) the original tissue boundary, which serves to draw light out of the original tissue region, thereby mimicking the loss of light due to escape from the tissue.
Figure 11a schematically depicts light diffusing toward an air/tissue surface, labeled as a flux, and diffusing away from the surface, labeled as an flux. The light is shown as a zig-zag trajectory, at relative to the axis, because on average each photon travels a distance for every movement along the axis. Consider a test point just below the surface within the tissue, . The fluence rate at this test point will be twice the value of flux passing along the axis contributed by both the and fluxes, in other words . As the light hits the surface, there is internal reflectance of about half the light . So the value of is . Hence, . The escaping flux is the observable reflectance . Hence, . Therefore, the fluence rate at the surface isevaluated at . However, the boundary condition is mimicked by assuming an infinite medium and placing a surface of symmetry outside the true tissue surface at a position , where is a negative value. For each true source of light within the tissue, one places a negative image source at a position symmetrically opposite the true source, outside the original tissue above this plane of symmetry. For example, if a point source of light is located at a depth within the tissue, then a negative point source of equal power is placed at a position .
The example situation illustrated in Fig. 11b uses 1-D planar diffusion theory to better illustrate the gradient of between a real source at and a negative image source at . When 3-D diffusion theory is used to consider the response to a point source that is close to the surface, the gradient in the region near the point of light delivery may not be perfectly linear, and the solution may not be accurate in this region. However, the solution yields an accurate solution distant from the point source, and the boundary condition discussed here is applicable to 3-D diffusion from a point source at .
The choice of is such that the gradient from the tissue surface to this virtual surface of symmetry is . Therefore, the escaping flux is equal to29, 30 and eliminating yields30 describes the boundary condition. The value of internal reflectance specifies the value . The value can be calculated: is the Fresnel reflectance as a function of the angle of incidence relative to the surface normal and the refractive index mismatch at the boundary. The calculation of is 0.4722 for an air/water interface and for an air/tissue interface . A very good approximation is is , calculated by Egan and Hilgeman9 as cited by Groenhuis Ferwerda, and Ten Bosch.10 It is very good for .
Therefore, for a generic air/tissue interface, the value of is , and the value of is . Therefore, . Hence, the position of the negative image source is placed at .
Once the positions and have been specified, the superposition of the two sources yields the net fluence rate in the tissue:, this gradient will drive a flux out of the tissue that predicts the escaping reflectance as a function of radial position : distant from the source at .11
A common situation is a narrow collimated beam of light illuminating a tissue. The collimated beam launches light into the tissue, and the light acts as if there is a point source at . Equations 34, 35 yield good approximations to the fluence rate and escaping flux at locations greater than distant from the origin at and greater than from the source at ( , ). Monte Carlo simulations are recommended for regions close to the source.
Diffusion Implementation in Numerical Methods
Partial Differential Equation Representation of Diffusion
Solution to the diffusion problem can be accomplished through formalized approaches to solving the partial differential equation itself. This can be done numerically or analytically, and the solutions can be derived from the time-domain equation:11, 12, 13 However, the validity and ease of analytic methods tend to break down where the exterior or interior boundaries become complex, thereby requiring more elaborate region integral approaches. Additionally, if the goal of the problem is the combination of clinical or biological imaging tissue volumes where both the boundaries are complex and the interior property values vary, then the optimal choice may transition away from hybrid analytic/numerical solutions toward purely numerical approaches.
Numerical methods to solve the diffusion equation are numerous, including finite difference, finite element and boundary element approaches.14, 15, 16, 17 The major difference between these approaches is in the formalized approach to solving the partial differential equation; however, the choice also leads to major differences in how the solution must be calculated. From a practical standpoint one of the key differences separating the numerical methods is the shape and origin of the mesh of points on which the region to be simulated is discretized. These differences can be seen in the representative meshes displayed in Fig. 12 .
Often the trickiest part of solving the solution in these domains is how to deal with the boundary between the diffusing regime, and what is outside this domain. Boundary conditions can make solving for analytic solutions more complex, and also the conditions must be validated against either measurements or a Monte Carlo type of simulation.
The choice of numerical approach is often driven by the experience with the method, as each requires an expert to create the code, yet often applications of the codes are better suited to different applications. In Table 4 , a summary of how the mesh and memory used scales with the size of the problem is shown. Additionally, the key advantages and weaknesses of each method are listed briefly.
The three major numerical methods for solving the diffusion equation are listed with associated advantage and disadvantage interms of how they scale as the problem gets larger. Here, N is the number of nodes (unknowns) in the mesh, and B is the bandwidth of the mesh. The memory and size of FEM and FDM assume that the mesh is banded (in 2-D), a banded solver such as LU decomposition is used,14, 15 and a sparse matrix solver is used in 3-D.
|Integratedanalyticmethod||Arbitrary||Small memory||Fast,conceptually theeasiest approach||Integration over highlyshaped regions andheterogeneties are difficult|
|Finitedifferencemethod||Cartesian||2-D: 3-D:||2-D: 3-D:||Ease of constructionand discretization,sparse/bandedmatrices||Difficult to handleirregular boundaries|
|Finite elementmethod||Arbitrary(triangular elementsin 2-D; tetrahedronsin 3-D)||2-D: 3-D:||2-D: 3-D:||Arbitrary shapes,flexible,sparse/bandedmatrices||Volume meshing in 3-Dcan be difficult forarbitrary shapes|
|Boundaryelementmethod||Arbitrary(line segments in2-D; triangularelements in 3-D)||Arbitrary shapes,nodes only onboundary||Full matrices;assumption of piece-wise constant regions|
Numerical Methods to Solve the Diffusion Equation
Finite difference method
Finite difference methods (FDM) approaches are routinely solved on a discretized Cartesian mesh, where the volume is arbitrary. However, since the mesh is Cartesian, then the volume must be discretized finely enough to allow accurate representation of the tissue boundaries. Usually this regular Cartesian grid used in discretization is the major limitation of finite difference methods, in that most complex boundaries cannot be discretized finely enough to represent smooth curved boundaries. The finite difference solution is obtained by solving the numerical molecule for a diffusive term, on the nodes given in the mesh.14, 15 The time-domain diffusion equation has a major limitation in this approach, in that the finite difference approach has large time requirements due to the iterative approach to converging on the solution, and so the added need to iterate each time step makes this method intractable for very large discretizations. Nonetheless, for highly complex partial differential equations, and very large problems, the finite difference approach has always been the standard to initiate test solutions, mainly because of its ease of construction and intuitive approach to discretization. The diffusion equation is simply discretized into points for the coordinates, and the fluence at each location calculated by a straightforward approximation of the derivative by a finite difference, such that it becomes:7 and that the optical properties and vary spatially and so must be discretized as well as the fluence). This finite difference molecule can then be solved sequentially at each node, and by iterating this over the entire mesh until there is little change between iterations, the solution is reached. This is the most intuitive manner to solve this, called the Jacobi method. Since the Jacobi finite difference method is plagued by slow convergence issues, many alternative methods have been produced to speed up the convergence. Most notable is the implicit solution, where the problem is solved as a matrix inversion, rather than iterative relaxation approach. This requires developing a full matrix representation of how each node in the mesh affects all other nodes and summing up these contributions numerically into a matrix that can be solved. Then the solution is reached simply by inverting this large matrix. This approach obviously has large memory demands for reasonably complex problems, and in contrast the explicit solver requires little memory. Several hybrid approaches have been used, with one notable one being called the Crank-Nicholson alternating difference implicit method. This approach cycles between explicit steps and implicit solution steps to speed up the convergence of the solution, but requires more work to create the initial solution. Most finite difference solvers for regular equations are freely available through academic shareware libraries on the Internet.
Multigrid solver methods have been developed as more elaborate approaches to solving the finite difference method faster.18 The central concept in these approaches is to solve the relaxation on a very coarse grid, and then allow rapid extrapolation of the solution from a coarse grid to a finer grid. The optimally designed approach will actually use several difference mesh grid resolutions, and start out solving the most course grid, then systematically solve the solution on further and further refinements of the mesh. Then, accurate solutions will typically require cycling back and forth on the course grid and fine grid a few times to allow rapid convergence with maximal efficiency. A uniquely high speed solver was developed by the researchers at National Oceanic and Atmospheric Administration (NOAA) and used for diffuse spectroscopy for more than a decade.19, 20, 21
Finite element method
The FEM approach to solve partial differential equations is based on solving the equation as integrated on arbitrary discrete elements, called basis functions , which map pieces of the entire domain between discrete points, called nodes,15, 16 as shown in Fig. 12b. This is often called the Galerkin method, and is developed in a formalized approach to discretizing any partial differential equation, resulting in a linear system of equations, which can then be solved by matrix methods. The diffusion equation is usually solved in this approach, in what is called the weak formulation, which is a theoretical framework simplifying the solution found, such that it can be solved in the linear formulation.
In the finite element formalism, the frequency domain solution can be discretized onto the basis functions multiplied by weighting coefficients , where is the discrete fluence at each node with weights , which are linear mapping functions. The discrete values of the fluence are determined as part of the solution process, given simple fixed weight functions. In the Galerkin method of weighted residuals, Eq. 37 is multiplied by an identical set of weighting functions and integrated over the entire problem domain to give:simply denotes integration, such that , and is the entire region of the imaging field over which the integrand is calculated. To avoid the need for second order differentiation of the basis functions, this equation is manipulated further through Green’s theorem to arrive at a form without second derivatives: is the surface normal from the surface being integrated by . This representation has two advantages: (1) the second derivative terms have been eliminated, and (2) the surface integral provides a vehicle for implementing the boundary conditions. In fact, its integrand is exactly the normal component of the flux through the boundary surface. At all node points within the finite element mesh that are not on its outer boundary, this term will be zero. In solutions of the prior equation, the integrand in this flux term has been typically approximated by: therefore defined as: is the normal direction from the surface, along the direction of , at all boundary locations . In this derivation, has a theoretical value near , for a tissue-air interface.22 While analytical forms of can be determined, it must be recognized that this coefficient is, at best, a rough estimate of the relative fluence rate gradient proportionality. In reality, diffusion theory may not match the measurements near the source or at the boundaries, due to the irradiance being highly anisotropic. So with the limitations in the accuracy of diffusion theory at the boundary, this parameter may not have any meaningful physical value that would allow a good match to measured data of photon transport in tissue. In practice, this coupling parameter is often chosen to allow a good match between the forward simulations and measured data from tissue phantoms.
Using the weak form of the earlier equation, and discretizing the diffusion coefficient and the absorption coefficient in the same manner as was done with the fluence, then the equation becomes:, , and denote the different weighting functions for fluence and tissue properties. The approach is used to create a numerical matrix map of the diffusion connectivity between the basis functions, where each element of is given by: are linear or nearly linear interpolation functions between neighboring nodes (i.e., linearly varying from 1 at the node to 0 at all neighboring nodes), so that all significant interactions only take place between neighboring nodes in this integration. Thus, the resulting matrix is then a continuous map of the diffusion equation, with values represented at the discrete node points, and with only local connectivity between nodes. The basis functions are used to interpolate the values to all points between the nodes. More complex basis function representations exist, but are often not used in simpler solvers.15
The approach to simplifying this equation is too lengthy to explain here, but is well described in Refs. 16, 17, 23. The Galerkin method results in a forward solution that reduces to a matrix equation representation of the diffusion equation:is the vector of fluence values at all nodes and is the vector of sources at each node, driving the predicted fluence. This fluence is then solved simply by inverting the connectivity matrix : is always a banded matrix with many zero elements, and so this can be solved quite efficiently using a specialized banded matrix solver. They key to efficiency here is then using a mesh that has the lowest bandwidth possible. The bandwidth is simply a measure of the node numbering, and indicates the maximum difference between the node numbers that are attached to one another, via being in the same triangular or tetrahedral element. Bandwidth reduction of a generated mesh is often required for efficient solution with a banded matrix solver. However, in 3-D, sparse matrix storage schemes are used for efficiency, and the computational burden increase is usually formidable. Yet several working solutions are freely available and commonly used. Generally though, mesh resolutions above 20,000 to 40,000 nodes can become too large to solve on a single computational engine, and must be solved with parallel computer processor implementation or via commercially optimized solvers.
The boundary conditions are applied by adding to the discretized version of this equation into the matrix equation.24, 25, 26, 27 This only has additive values at the boundary nodes, and is equal to zero at all internal nodes. Zero fluence boundary conditions are often used only in the simplest of simulations, and are easily implemented by setting the off-diagonal matrix elements to zero along the rows corresponding to boundary nodes, and then allowing the diagonal value to be unity. Extrapolated boundary conditions are implemented by adding additional nodes to the mesh exterior and having measurement surface points within the mesh at the location of the true boundary.
Once the source and boundary conditions are set up in the matrix and vector, the result from inversion of is a discretized map of the optical fluence pattern , which is accurate if the domain has been discretized finely enough. Again, this approach has strict limits on the discretization and can run into memory problems for arbitrarily large sized domains.
Boundary element method
The boundary element method (BEM) is a less explored area of numerical methods, but is quite useful for problems where the interior properties are known to be homogeneous. This method provides a unique approach where the exterior shape can be complex, and the fluence at the boundaries can then be estimated given knowledge of the homogeneous interior regions. Because of the homogeneity constraints of this approach, it is best implemented when the exterior shapes of tissues and the boundaries of interior organs are well known a priori from some structural imaging modality, and then the light fluences at the boundaries of these tissues can be calculated. This has been shown to be fundamentally useful in electrical impedance problems,28 and has been extended to optical diffusion problems in shape recovery,29 fluorescence prediction30 and spectral region estimation.31 The type of mesh used is simply the exterior boundary mesh of points to discretize the region, as shown in Fig. 12c.
The boundary element solution is actually quite similar to the method of analytic integration, in that its derivation begins by assuming that there is a Green’s function solution for the partial differential equation to be solved, and that this solution could then be used to solve for an arbitrary solution by integration over the boundary shape. Under homogeneity conditions, the diffusion equation is simplified into the modified Helmholtz equation for which the Green’s function can be found to be:, and is constant in subdomain .
The final version of this expansion for the diffusion equation leads to:31and in this equation have elements: is the source distribution is the discrete integral over the boundary angular space, and is the solid angle enclosed by the boundary at node . The complete context of the BEM equation derivation is beyond the scope of this review, but can be found in the papers referenced earlier. The BEM provides a major benefit in obtaining the discretization of the domain, especially in 3-D by requiring only surface meshes. Its main weakness is the use of full matrices, which are required, as compared to the banded or sparse matrices used for FEM calculations.
Applications of Diffusion Modeling of Light Transport
Estimation of Tissue Properties
There are several classes of image or property recovery problems that are solved with these forward diffusion models. The commonality in them is that usually there is some measured dataset , which is a function of the fluence at the boundary . This can be simulated by a diffusion theory calculation , such as the phase, logarithm of the intensity, relative intensity, intensity ratio, mean transit time, etc. These measurements are then used to predict the properties within some unknown region of the tissue. The terminology commonly used is that the forward problem involves obtaining the fluence given a known set of optical properties. The problem of predicting the optical properties given some measure of the fluence is then called the inverse problem. In the inverse problem, regions to be recovered could be:
1. the entire region, assumed to be homogeneous
2. layers of tissue
3. locally embedded regions within the bulk
4. recovery of properties at every point in the tissue.
The mathematical formulation would be as follows. Assuming that the initial guess is close to the true vector of property values, the Newton method would assume that the goal is to minimize the square difference between the measured and calculated datasets:is the matrix with elements defined by the derivative of at each measurement point with respect to the property values or at each node, giving matrix elements: , such that: ; however, since is rarely a square matrix (i.e., the number of nodes is almost never exactly equal to the number of measurements), then the update equation becomes: is to allow it to be regularized, using a scalar regularization parameter , which is either empirically or algebraically derived: largely depend on the size and how well posed the matrix is. In the smallest case, this could be a single region problem with multiple measurements , leading to a matrix that is , and readily solved without regularization. In situations where the regions are in layers or where there are “lumps” to be recovered, then the matrix is still small, and may or may not require regularization. An example of this is shown in Fig. 13 , where the interior of a rat cranium is being imaged with NIR light. Figure 13a shows the MRI image of the cranium, where the brain is visible at the bottom as a lighter gray region, the skull is seen as dark black regions, and the muscle is seen as darker gray through the cranium. Figure 13b is a FEM mesh created to allow diffusion simulation of this tissue, separating the muscle and brain into two distinct regions. The diffuse reconstructed values of the tissue absorption coefficient are shown in Fig. 13d where the differences in absorption between brain and muscle tissue are apparent, even in this ill-posed diffuse tomography image. In the full image reconstruction problem, where there are nodes and measurements, the matrix is and can require considerable effort to regularize and invert, and still will lead to diffuse blurry images, of the type shown in Fig. 13d. Diffuse tomography always leads to these low resolution types of images, because the forward diffusion transport process makes the image reconstruction problem highly ill-posed and nonunique. Virtually all diffuse tomography imaging is limited in this manner, in that the abilty to recover sharp features is significantly diminished. Resent work in NIR imaging has been in the area of using diffuse spectroscopy as an adjuvant to standard imaging systems, such as MRI or CT. An example of how this would be done is shown in Fig. 13c, where the same rodent brain as before is recovered as a predefined region in the absorption coefficient. In this approach, it is possible to obtain quantitatively accurate bulk tissue values, because the inversion problem is no longer one dominated by image reconstruction, since the recovery only involves two regions instead of regions. It has been shown that image-guided recovery of NIR measurements allows much more accurate spectroscopy of tissue, if there is a meaningful way to combine NIR with an imaging system to give the prior information.
Recovery of tissue values at a single wavelength are readily achieved, and computer codes such as TOAST (University College London34) and NIRFAST (Dartmouth College and University of Exeter35) are freely available for this problem.
Diffuse spectroscopy of tissue is often the primary goal in diffuse optics, and the nature of the problem can be transmission spectroscopy for absorption or scatter, or luminescence spectroscopy for bioluminescence, fluorescence, or phosphorescence. Any problem that requires multiple wavelength fitting to a predetermined extinction or absorption spectrum can fall into the secondary inversion problem of spectroscopy.
In absorption spectroscopy, the goal is to then fit the measured or computed absorption coefficients (e.g., from solving the inverse problem in previous section) to the known molar absorption spectra of the features that absorb in tissue, such as hemoglobins, breakdown products of heme, water, lipids, cytochromes, and potentially other absorbers that are endogenous or exogenous. The problem is formulated as:, with elements of molar absorption for each wavelength. The inversion is then:14a , to show the trough in absorption between hemoglobin and water, bounding the diffuse regime on the red and NIR regions, respectively. Outside the range of , it is difficult to believe that diffusion theory would be an acceptable approximation, because in many tissues the absorption coefficient is of similar value to the transport scattering coefficient, as seen in this graph.
Spectroscopy of tissue with this approach to fitting the spectra has been shown to be enormously successful in breast tumor spectroscopy,36 brain tissue monitoring,37 and muscle physiology,38, 39 to name a few applications. Image-guided diffuse optical spectroscopy has recently been established as a useful tool to image cancer tumors, and uses the principles of region-based estimation to reduce the size of the Jacobian and allow more accurate recovery of the interior regions. Additionally, incorporating the spectral inversion into the image recovery problem has been demonstrated to allow more accurate estimation of the tissue values.40, 41, 42 It is generally assumed that addition of any constraints into the inversion problem, in the way of spectral or spatial prior knowledge, will allow a well-posed inverse problem and potentially more accurate estimation of tissue properties.
Elastic Scatter Spectroscopy
In scatter spectroscopy, the nature of the scatter is elastic scattering, and known to fit to Mie theory under certain circumstances. Several research groups have shown the ability to fit the scattered spectrum to the effective particle size and number density, under the assumption of a Mie-like scattering model. Alternatively, the spectra is thought to fit to a Mie-like spectra with some Rayleigh scatter influence, which has been shown to be near:or , and ultimately these parameters are not intuitive, as the values for and are arbitrary from an analytic approximation to the Mie theory. The important spectra to fit comes from the Mie theory: is the particle size histogram, with diameter and index ratio , with scattering efficiency from Mie theory43 and calculated anisotropy parameter .44 An illustration of the change in the scattering spectrum with size of the scatterer is shown in Fig. 14b.
Fitting remitted or transmitted scatter spectra to models of the scatterer particle size has become an extremely active area of research, with extensive preclinical and clinical data indicating that scattering spectra can be used to estimate scatterer composition and size in certain cases.45, 46, 47, 48, 49, 50
Luminescence Spectroscopy: Fluorescence, Phosphorescence, and Bioluminescence
Fitting of luminescence spectroscopy is less common, but becoming increasingly important for recovery of a unique concentration value,51 and for applications where the spectral shape is altered on remission from the tissue, providing some measure of depth from which the signal is coming.52 This fitting is identical to the absorption spectra fitting, although often there is need for addition of a background spectrum or contaminant spectra that must be included to allow accurate fitting. Examples of fluorescence spectra from relevant fluorophores, which are used routinely in tissue imaging, are shown in Fig. 14c. The amplitudes of each are normalized to unity at their maximum, but the yield of fluorescence is directly proportional to (1) the molar extinction coefficient at the excitation light wavelength, (2) the intensity of the excitation light, (3) the fluorescence quantum yield of the fluorophore within tissue, and (4) the concentration of the fluorophore present. A reduction in any of these four things will lead to reduced fluorescence from the tissue sampled.
Fluorescence tomography for small animal imaging has been largely established by Zacharakis and by Ntziachristos 53, 54, 55, 56 in recent years. Yet the majority of small animal luminescence imaging with commercial instrumentation is done without much consideration of diffusion issues or compensation for diffusion, which can suffice for relative comparisons of fluorophore production or uptake. However, when the location of the fluorophore changes, or when the tissue shape change is relevant, then more extensive modeling is required. Additionally, it has become well established in recent years that if surface measurements are sufficient, then the reflectance geometry of having the source and detector on the same side of the tissue is good enough for imaging. However, if the fluorophore is deeper in the tissue, then the transmission geometry provides a superior signal-to-background value.57 This transmission-based measurement for uptake of fluorophores in tumor tissue will become more utilized as this mode of measurement becomes standard in commercial small animal imaging systems.
This work focuses on the fundamental concepts leading to the mathematical and numerical estimations of light diffusion in tissue. While the theory can be developed from a number of different perspectives, the focus here is on phenomenological description, allowing physical insights into the problem. Analytic derivations are kept to a minimum while outlining the overall concepts. In the numerical methods presented, the focus is on providing a practical outlines of the strengths and weaknesses of the different methods. This review is a brief introduction to what has become a well developed mature field of diffuse light transport modeling. The references listed are a good introduction to the fundamentals and applications that have emerged in the past few decades of work.