## 1.

## Introduction

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.

## 2.

## Analytic Modeling of Diffuse Light Transport

## 2.1.

### 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
$P$
and light transport
$T$
differ for a point source, a line source, and a planar source, but the fluence rate
$\varphi =PT$
always has the same units. For a point source,
$P$
has units of
$W$
and
$T$
has units of
${\mathrm{cm}}^{-2}$
. For a line source,
$P$
has units of W/cm and
$T$
has units of
${\mathrm{cm}}^{-1}$
. For a planar source,
$P$
has units of
$\mathrm{W}\u2215{\mathrm{cm}}^{2}$
and
$T$
is dimensionless.

## Table 1

Tissue optical properties.

Quantity | Symbol | Units |
---|---|---|

Absorption coefficient | ${\mu}_{a}$ | ${\mathrm{cm}}^{-1}$ |

Scattering coefficient | ${\mu}_{s}$ | ${\mathrm{cm}}^{-1}$ |

Anisotropy of scattering | $g$ | dimensionless |

Refractive index | $n$ | dimensionless |

## Table 2

Transport properties used in diffusion theory.

Quantity | Symbol | Units |
---|---|---|

Reduced scattering coefficient | ${\mu}_{s}^{\prime}={\mu}_{s}(1-g)$ | ${\mathrm{cm}}^{-1}$ |

Transport mean free path | ${\mathit{MFP}}^{\prime}=1\u2215({\mu}_{a}+{\mu}_{s}^{\prime})$ | cm |

Diffusion length | $D={\mathit{MFP}}^{\prime}\u22153$ | cm |

Optical penetration depth | $\delta =\sqrt{D\u2215{\mu}_{a}}$ | cm |

## Table 3

Parameters for light transport.

Quantity | Symbol | Units |
---|---|---|

Radiant power | $P$ | W, W/cm, $\mathrm{W}\u2215{\mathrm{cm}}^{2}$ |

Radiant energy | $Q$ | J |

Transport | $T$ | ${\mathrm{cm}}^{-2}$ , ${\mathrm{cm}}^{-1}$ , dimensionless |

Fluence rate | $\varphi =PT$ | $\mathrm{W}\u2215{\mathrm{cm}}^{2}$ |

Fluence | $\psi =QT$ | $\mathrm{J}\u2215{\mathrm{cm}}^{2}$ |

Speed of light in tissue | $c={c}_{0}\u2215n$ | cm/s |

Concentration of radiant energy | $C=\varphi \u2215c$ | $\mathrm{J}\u2215{\mathrm{cm}}^{3}$ |

## 2.2.

### 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
$1\u2215[{\mu}_{a}+{\mu}_{s}(1-g)]$
in front of the fiber. This distance is called the transport mean free path, or
${\mathit{MFP}}^{\prime}$
(sometimes called
${l}^{*}$
). The figure illustrates diffusion theory and the Monte Carlo approach agreement at distances exceeding one
${\mathit{MFP}}^{\prime}$
from the tip of the fiber and exceeding one
${\mathit{MFP}}^{\prime}$
from the apparent source of diffusion located one
${\mathit{MFP}}^{\prime}$
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,
${\varphi}_{\mathit{MC}}(x,y,z)$
[
$\mathrm{W}\u2215{\mathrm{cm}}^{2}$
per W], which is the transport
${T}_{\mathit{MC}}(x,y,z)$
$[1\u2215{\mathrm{cm}}^{2}]$
When a photon reaches a depth
$z$
that exceeds
${\mathit{MFP}}^{\prime}$
, 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
$P(x,y,z)$
[W per voxel]. After launching many photons, the Monte Carlo simulation yields both a
${\varphi}_{\mathit{MC}}(x,y,z)$
and a distributed source
$P(x,y,z)$
. A point spread function or Green’s function,
$G({x}^{\prime},{y}^{\prime},{z}^{\prime},x,y,z)$
, based on diffusion theory is convolved against
$P({x}^{\prime},{y}^{\prime},{z}^{\prime})$
to yield the fluence rate,
${\varphi}_{\mathit{DT}}(x,y,z)=G\otimes P$
, throughout the medium due to diffusion. The final answer is
$\varphi ={\varphi}_{\mathit{MC}}+{\varphi}_{\mathit{DT}}$
.

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 ${\mathit{MPF}}^{\prime}$ corresponds to less than 1 or $2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ . 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
${\mu}_{s}(1-g)\u2215{\mu}_{a}$
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
${\mu}_{s}(1-g)$
rather than on the specific values of
${\mu}_{s}$
and
$g$
. The total diffuse reflectance from a tissue that has
${\mu}_{s}(1-g)\u2215{\mu}_{a}=10$
is 0.252 (based on Monte Carlo simulations), where the mismatched tissue/air surface boundary has
${n}_{\text{tissue}}\u2215{n}_{\text{air}}=1.4$
. A study of reflectance versus
${\mu}_{s}(1-g)\u2215{\mu}_{a}$
where
${\mu}_{s}$
and
$g$
were varied but
${\mu}_{s}(1-g)$
was conserved, indicated that reflectance becomes independent of
$g$
when reflectance exceeds
$\sim 0.40$
$[{\mu}_{s}(1-g)\u2215{\mu}_{a}\phantom{\rule{0.3em}{0ex}}\text{equal}\sim 20]$
.^{4} So the reflectance of a tissue should exceed at least 25% and preferably 40%, or
${\mu}_{s}(1-g)\u2215{\mu}_{a}$
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
${\mu}_{s}(1-g)\u2215{\mu}_{a}$
, 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 ${\mu}_{s}(1-g)\u2215{\mu}_{a}\u2aa110$ . 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.

## 2.3.

### Diffusion Theory

## 2.3.1.

#### Fluence rate and concentration of optical energy

The fluence rate $\varphi $ $[\mathrm{W}\u2215{\mathrm{cm}}^{2}]$ is proportional to the concentration of optical energy $C$ $[\mathrm{J}\u2215{\mathrm{cm}}^{3}]$ :

where $c$ is the speed of light in the medium, $c={c}_{o}\u2215n$ [cm/s], ${c}_{o}=2.998\times {10}^{10}\phantom{\rule{0.3em}{0ex}}\mathrm{cm}\u2215\mathrm{s}$ , and $n$ is the refractive index of the medium. For example, consider a uniform beam of $1\phantom{\rule{0.3em}{0ex}}\mathrm{W}\u2215{\mathrm{cm}}^{2}$ irradiance passing through a $1\text{-}{\mathrm{cm}}^{3}$ cube of water (see Fig. 2 ). If the fluence rate $\varphi $ in the nonscattering water $(n=1.33)$ equals $1\phantom{\rule{0.3em}{0ex}}\mathrm{W}\u2215{\mathrm{cm}}^{2}$ at $500\text{-}\mathrm{nm}$ wavelength, then the corresponding concentration of photons is $0.18\phantom{\rule{0.3em}{0ex}}\mathrm{pM}$ $(0.18\times {10}^{-12}\phantom{\rule{0.3em}{0ex}}\text{moles}\u2215\text{liter})$ :## Eq. 2

$$\varphi \frac{n}{{c}_{0}}\frac{\lambda}{h{c}_{0}}\frac{1000}{{N}_{\mathrm{Av}}}=(1\phantom{\rule{0.3em}{0ex}}\mathrm{W}\u2215{\mathrm{cm}}^{2})\frac{1.33}{{(3\times {10}^{10}\phantom{\rule{0.3em}{0ex}}\mathrm{cm}\u2215\mathrm{s})}^{2}}\frac{(500\times {10}^{-7}\phantom{\rule{0.3em}{0ex}}\mathrm{cm})}{(6.63\times {10}^{-34}Js)}\frac{1000\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{3}\u2215\text{liter}}{(6.02\times {10}^{23}\phantom{\rule{0.3em}{0ex}}{\text{mole}}^{-1})}=0.18\times {10}^{-12}\mathrm{M}.$$*in vacuo*speed of light, indicates the number of photons per J energy. The factor 1000 indicates the conversion $1000\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{3}\u2215\text{liter}$ , and ${N}_{\mathrm{Av}}$ is Avagadro’s number for the molecules per mole.

## 2.3.2.

#### 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 $J$ $[\mathrm{W}\u2215{\mathrm{cm}}^{2}]$ down a concentration gradient $\partial C\u2215\partial x$ $[(\mathrm{J}\u2215{\mathrm{cm}}^{3})\u2215\mathrm{cm}]$ is

where $\chi $ is the diffusivity $[{\mathrm{cm}}^{2}\u2215\mathrm{s}]$ equal to $cD$ for light, with the diffusion length $D={\mathit{MFP}}^{\prime}\u22153\phantom{\rule{0.3em}{0ex}}\left[\mathrm{cm}\right]$ . In general, $\chi $ can refer to the diffusion of heat, molecules, etc., or radiant energy. For light, $C=\varphi \u2215c$ and $\chi =cD$ , and the product $\chi C=\left(cD\right)(\varphi \u2215c)=D\varphi $ . 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 $x$ 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 $dA$ at $x=0$ , $z=0$ , oriented such that the normal to the window is aligned with the $z$ 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 ${\mu}_{s}^{\prime}$ $\left[{\mathrm{cm}}^{-1}\right]$ . The fluence rate upstream from the window along the gradient encounters the scatterers within an incremental volume $dV$ and a power is scattered ${\mu}_{s}^{\prime}\varphi dV$ [W]. A fraction $\mathrm{cos}\left(\theta \right)dA\u2215\left(4\pi {r}^{2}\right)$ of this scattered power is directed toward the window, which is a distance $r$ from the volume $dV$ , where $\theta $ 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 $\mathrm{exp}(-r\u2215{\mathit{MFP}}^{\prime})$ , 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 $dA$ yields the net flux density $J$ $[\mathrm{W}\u2215{\mathrm{cm}}^{2}]$ passing through the incremental window:

## Eq. 4

$$J={\int}_{z=-\infty}^{0}{\int}_{x=0}^{\infty}{\mu}_{s}^{\prime}\varphi (x,z)\frac{\mathrm{cos}\left(\theta \right)}{4\pi {r}^{2}}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-r\u2215MF{P}^{\prime})\mathrm{sign}(-z)2\pi xdxdz.$$Fick’s first law holds if the gradient of fluence rate is linear within the two hemispherical regions of integration, i.e., $\varphi \left(z\right)=\varphi \left(0\right)+az$ , where $z$ is the position along the axis perpendicular to the window and $a$ is a constant. This situation is called diffusion theory, and is sometimes called the ${P}^{1}$ approximation. If, however, there is too much curvature in the gradient of $\varphi $ , then the $\varphi $ near the window must be described by higher order terms, $\varphi \left(z\right)=\varphi \left(0\right)+az+b{z}^{2}+\dots $ . An example of this is the so-called ${P}^{3}$ 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 ${P}^{3}$ approximation is a more appropriate method to describe the transport. However, if the source and point of observation are separated by more than a ${\mathit{MFP}}^{\prime}$ , the ${P}^{1}$ 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 $\varphi $ across the boundary will be continuous, although the rate of energy deposition ${\mu}_{a}\varphi $ can be discontinuous.

## 2.3.3.

#### Time-resolved diffusion

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:

## Eq. 5

$$C(r,t)=Q\frac{\mathrm{exp}[-{r}^{2}\u2215\left(4\chi t\right)]}{{\left(4\pi \chi t\right)}^{3\u22152}}.$$To describe the diffusion of fluence rate $\varphi $ after deposition of a point source of radiant energy, $Q$ $\left[J\right]$ , the speed of light is introduced, since $\varphi =cC$ . Also, the diffusivity $\chi $ is replaced by $cD$ . Finally, the role of absorption is added by including the term $\mathrm{exp}(-{\mu}_{a}ct)$ , where $ct$ is the path length of individual photons at time $t$ after the impulse, and ${\mu}_{a}$ $\left[{\mathrm{cm}}^{-1}\right]$ is the absorption coefficient. The result is:

## Eq. 6

$$\varphi (r,t)=cQ\frac{\mathrm{exp}[-{r}^{2}\u2215\left(4cDt\right)]}{{\left(4\pi cDt\right)}^{3\u22152}}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-{\mu}_{a}ct).$$## 2.3.4.

#### Steady-state diffusion

Integration of the prior time-resolved $\varphi (r,t)$ in Eq. 6 over all time yields the radiant exposure $\psi \left(r\right)$ $[\mathrm{J}\u2215{\mathrm{cm}}^{2}]$ , which does not depend on time. The direct integration of Eq. 6 using rules of integration is not immediately obvious. Intuition indicates that $\psi $ should fall exponentially with distance from a source. A numerical integration of $\varphi (r,t)$ over all time at each position $r$ (not shown here) confirms that at large $r$ , the radiant exposure $\psi \left(r\right)$ falls as $\mathrm{exp}(-r\u2215\text{constant})$ . The constant is given the symbol $\delta $ [cm] and is called the optical penetration depth. The factor $1\u2215\delta $ is called the effective attenuation coefficient ${\mu}_{\mathrm{eff}}$ . The integration becomes easier by considering conservation of energy. The proper form of $\psi $ is such that the integral of energy deposition over an infinite volume equals the total energy deposited $Q$ . The energy deposition ${\mu}_{a}\psi $ $[\mathrm{W}\u2215{\mathrm{cm}}^{3}]$ is integrated over the infinite spatial domain:

## Eq. 7

$${\int}_{\infty}{\mu}_{a}\psi \left(r\right)dV={\int}_{r=0}^{\infty}{\mu}_{a}\frac{\mathrm{exp}(-r\u2215\delta )}{\text{factor}}4\pi {r}^{2}dr=Q,$$## Eq. 8

$$\psi \left(r\right)={\int}_{t=0}^{\infty}\varphi (r,t)dt=Q\frac{\mathrm{exp}(-r\u2215\delta )}{4\pi {\mu}_{a}{\delta}^{2}r}=Q\frac{\mathrm{exp}(-r\u2215\delta )}{4\pi Dr}.$$This $\psi \left(r\right)$ is the impulse response of radiant exposure to an impulse of energy deposition $Q$ $\left[J\right]$ . If one has a repetitive sequence of impulses at a frequency $f$ , $[1\u2215\mathrm{s}]$ the time averaged power is $P=fQ$ [W]. In the limit of $f\to \text{infinity}$ and $Q\to 0$ such that $P$ is maintained at a constant average power, the superposition of the time-resolved responses $\varphi (r,t)$ to each impulse $Q$ blurs into a steady-state fluence rate proportional to $\psi \left(r\right)$ , with $Q$ replaced by $P$ :

This equation describes steady-state diffusion of light from a point source of power $P$ . Figure 4 gives an example of steady-state diffusion of light.The transport factor $T$ equals $\varphi \u2215P$ and characterizes the light transport in a tissue independent of the strength of the source. The following equations list the transport $T$ for 1-D planar diffusion from a planar source ${P}_{\text{planar}}$ $[\mathrm{W}\u2215{\mathrm{cm}}^{2}]$ , for 2-D cylindrical diffusion from a line source ${P}_{\mathrm{cyl}}$ [W/cm], and for 3-D spherical diffusion from a point source, ${P}_{\mathrm{sph}}$ [W], derived as outlined before from conservation of energy:

## 2.3.5.

#### Modulated diffusion

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:

## Eq. 11

$$P\left(t\right)={P}_{o}[1+{M}_{o}\phantom{\rule{0.2em}{0ex}}\mathrm{sin}\left(\omega t\right)],$$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) is

## 12.

## 14.

## 15.

^{5}The consequence of these expressions is that the fluence rate seen at the detector has a steady-state component ${P}_{o}{T}_{ss}\left(r\right)$ , plus a modulated component characterized by a modulation $M$ [dimensionless] and a phase $\Phi $ [radians]:

## 16.

For very low frequencies of modulation as $\omega \to 0$ , the factor $a$ approaches 1 and the factor $b$ approaches 0. Therefore, ${k}^{\u2033}$ approaches $1\u2215\delta $ and ${k}^{\prime}$ approaches 0. $M$ approaches ${M}_{o}$ , and $\Phi $ approaches 0. So the transport becomes

## Eq. 17

$$\varphi (r,t)={P}_{o}{T}_{ss}\left(r\right)[1+{M}_{o}\phantom{\rule{0.2em}{0ex}}\mathrm{sin}\left(\omega t\right)]=P\left(t\right){T}_{ss}\left(r\right),$$For very high frequencies of modulation as $\omega \u2aa2{\mu}_{a}c$ , the factor $a$ approaches ${[\omega \u2215\left({\mu}_{a}c\right)]}^{1\u22152}$ , and the factor $b$ approaches $\pi \u22154$ . Since $\mathrm{cos}(\pi \u22154)=\mathrm{sin}(\pi \u22154)=1\u2215\sqrt{2}$ , the values ${k}^{\u2033}$ and ${k}^{\prime}$ approach the same limiting value:

## Eq. 18

$${k}^{\u2033}={k}^{\prime}=\frac{1}{\delta \sqrt{2}}{\left(\frac{\omega}{{\mu}_{a}c}\right)}^{1\u22152}.$$## 2.4.

### 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 $P$ of light is launched isotropically at a position ${x}_{s},{y}_{s},{z}_{s}$ , the transport to a detector at position ${x}_{d},{y}_{d},{z}_{d}$ over a distance ${r}_{s\to d}$ is

and the fluence rate at the detector is ${\varphi}_{d}=P{T}_{s\to d}$ . The background optical properties are ${\mu}_{ao}$ and ${\mu}_{so}^{\prime}$ , yielding a background diffusion constant $D=(1\u22153)\u2215({\mu}_{ao}+{\mu}_{so}^{\prime})$ , and a background optical penetration depth $\delta ={(D\u2215{\mu}_{ao})}^{1\u22152}$ .Now consider the presence of a small absorbing object of volume $V$ at position ${x}_{p},{y}_{p},{z}_{p}$ with an extra absorption $\Delta {\mu}_{a}$ , such that its absorption is ${\mu}_{a}={\mu}_{\mathrm{ao}}+\Delta {\mu}_{a}$ . The fluence rate at the object is ${\varphi}_{\mathrm{obj}}=P{T}_{s\to p}$ , where ${T}_{s\to p}$ is computed using a source-object distance ${r}_{s\to p}={[{({x}_{p}-{x}_{s})}^{2}+{({y}_{p}-{y}_{s})}^{2}+{({z}_{p}-{z}_{s})}^{2}]}^{1\u22152}$ . The extra energy deposition in the absorbing object, above the background absorption, is $\Delta {\mu}_{a}{\varphi}_{\mathrm{obj}}$ $[\mathrm{W}\u2215{\mathrm{cm}}^{3}]$ . The total extra power deposited in the object is $\Delta {\mu}_{a}{\varphi}_{\mathrm{obj}}V$ [W], which has the units of power. Because the object is absorbing, the object behaves as a source of negative power ${P}_{p}$ :

This perturbing negative power propagates in a spherically symmetric manner throughout the tissue and reaches the detector, and yields a perturbing fluence rate at the detector, ${P}_{p}{T}_{p\to d}$ . The net fluence rate at the detector isThe first term on the right is positive. The second term on the right is negative for an absorbing object with $\Delta {\mu}_{a}>0$ . If $\Delta {\mu}_{a}<0$ , then ${P}_{p}$ is positive and the perturbation ${P}_{p}{T}_{p\to d}$ 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 ${\mathbf{P}}_{j}$ lists all sources. A vector ${\mathbf{P}}_{i}$ lists only the virtual sources, and is a subset of ${\mathbf{P}}_{j}$ . All the ${N}_{0}$ real sources are called ${\mathbf{P}}_{j}(j=1:{N}_{0})$ , and the ${N}_{i}$ virtual sources due to the perturbing objects are called ${\mathbf{P}}_{i}(i=1:{N}_{i})={\mathbf{P}}_{j}(j={N}_{0}+1:{N}_{j})$ . A perturbation matrix $\mathbf{M}$ is defined with ${N}_{j}$ columns to match ${\mathbf{P}}_{j}$ and ${N}_{i}$ rows to match ${\mathbf{P}}_{i}$ . The $\mathbf{M}$ matrix consists of elements that scale each source in ${\mathbf{P}}_{j}$ to yield a contribution to an updated value of ${\mathbf{P}}_{i}$ :

or## Eq. 24

$${\mathbf{P}}_{j}(j={N}_{0}+1:{N}_{j})={\mathbf{P}}_{i}(i=1:{N}_{i})+[\mathbf{M}{\mathbf{P}}_{j}-i(i=1:{N}_{i})]\u22153.$$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 $\mathbf{S}$ and the elements in the $\mathbf{M}$ matrix. One is due to absorption and one is due to scattering. Consider a spherical object of radius $R$ that is both absorbing and scattering. The volume of the object can be subdivided into a set of ${N}_{a}$ subvolumes that fill the object’s volume, and each is called a “virtual absorbing volume source.” The vector ${P}_{j}(j={N}_{0}+1:{N}_{0}+{N}_{a})$ is filled with these virtual absorbing volume sources, in this example expressed as small spheres, each of radius $a$ and subvolume ${V}_{i}=(4\u22153)\pi {a}^{3}$ . The radius $a$ is usually chosen to approximately equal ${\mathit{MFP}}^{\prime}\u22152$ , which sets the number of small spheres required to fill the object’s volume ${N}_{a}$ , but the total volume of all the small spheres must equal the volume of the object $V={N}_{a}(4\u22153)\pi {a}^{3}$ , so $a={[3V\u2215\left(4\pi {N}_{a}\right)]}^{1\u22153}$ . The elements of the $\mathbf{M}$ matrix for $i=1:{N}_{a}$ are specified for each $j$ ’th source by:

## Eq. 25

$$\mathbf{M}(j,i)=-[\Delta {\mu}_{a}\left(i\right)+\Delta {\mu}_{s}^{\prime}\left(i\right)\frac{{\mu}_{ao}}{{\mu}_{so}^{\prime}}]V\left(i\right)T(j,i).$$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 $\Delta {\mu}_{si}^{\prime}$ . The surface is divided into a set of ${N}_{s}$ incremental surfaces, each with a surface area $\Delta {A}_{i}$ . The sum of the incremental surfaces equals the total surface area of the object. The vector ${\mathbf{P}}_{j}(j={N}_{0}+{N}_{a}+1:{N}_{j})$ , where ${N}_{j}={N}_{o}+{N}_{a}+{N}_{s}$ , is filled with these virtual scattering surface sources. The elements of the $\mathbf{M}$ matrix for $i={N}_{a}+1:{N}_{i}$ , where ${N}_{i}={N}_{a}+{N}_{s}$ , is specified for each $j$ source by:

## Eq. 26

$$\mathbf{M}(j,i)=\frac{\Delta {\mu}_{a}+\Delta {\mu}_{s}^{\prime}}{{\mu}_{so}^{\prime}}\Delta A\left(i\right)\frac{D}{d}[{T}_{1}(j,i)-{T}_{2}(j,i)],$$*outside*the object at the $i$ ’th virtual source position, and ${T}_{2}(j,i)$ is the transport from the $j$ ’th source to a position $d\u22152$ just

*inside*the object at the $i$ ’th virtual source position. The distance $d$ is usually chosen to equal ${\mathit{MFP}}^{\prime}$ , 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 $(D\u2215d)({T}_{1}-{T}_{2})$ 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 ${\mathbf{P}}_{i}$ 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 ${\mathbf{P}}_{i}$ 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 $\mathbf{M}$ , a special case occurs when $i$ equals $j$ , and the self-perturbation by a virtual source must be evaluated. This term is here called ${\mathbf{M}}_{j,j}$ . For the scattering surface sources, the value of ${\mathbf{M}}_{j,j}$ is zero, because a point source cannot exert a gradient across itself. For the absorbing volume sources of radius $a$ , the value of ${\mathbf{M}}_{j,j}$ is calculated:

## Eq. 27

$${\mathbf{M}}_{j,j}={\int}_{r=0}^{a}\Delta {\mu}_{a}\frac{\mathrm{exp}(-r\u2215\delta )}{4\pi Dr}2\pi {r}^{2}dr=\frac{\Delta {\mu}_{a}}{{\mu}_{a0}}[1-(1+\frac{a}{\delta})\mathrm{exp}(-a\u2215\delta )].$$Figure 9 shows an example of a $0.4\text{-}\mathrm{cm}$ -diam. spherical object with adding absorber ( $\Delta {\mu}_{a}=1\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , $\Delta {\mu}_{s}^{\prime}=0$ , ${\mu}_{a0}=0.1\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , ${\mu}_{s}^{\prime}=10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ ). The background light field ${\varphi}_{0}$ is shown in Fig. 9a. The light source is located at $(x,y,z)=(-1,0,0)$ . The perturbation itself, ${\varphi}_{p}$ , is shown in Fig. 9b. The object casts a spherical zone of negative perturbation around the object, about $-10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-2}$ effect ( $\mathrm{W}\u2215{\mathrm{cm}}^{2}$ per W of the source) within the object and $-5\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-2}$ effect in the region surrounding the object. Figure 9c shows the fractional perturbation, ${\varphi}_{p}\u2215{\varphi}_{0}$ [dimensionless], illustrating that the perturbation is about $-0.30$ within the object and cases a shadow of perturbation about $-0.10$ to $-0.20$ 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 ( $\Delta {\mu}_{s}^{\prime}=10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , $\Delta {\mu}_{a}=0$ , ${\mu}_{a0}=0.1\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , ${\mu}_{s}^{\prime}=10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ ). 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 ${\varphi}_{p}\u2215{\varphi}_{0}$ 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}

## 2.5.

### Boundaries

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 $J$ flux, and diffusing away from the surface, labeled as an $I$ flux. The light is shown as a zig-zag trajectory, at $60\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$ relative to the $z$ axis, because on average each photon travels a distance $L=2\Delta z$ for every $\Delta z$ movement along the $z$ axis. Consider a test point just below the surface within the tissue, $z=0$ . The fluence rate $\varphi (z=0)$ at this test point will be twice the value of flux passing along the $z$ axis contributed by both the $I$ and $J$ fluxes, in other words $\varphi \left(0\right)=2(I+J)$ . As the light hits the surface, there is internal reflectance ${r}_{i}$ of about half the light $({r}_{i}\approx 0.50)$ . So the value of $I$ is ${r}_{i}J$ . Hence, $\varphi \left(0\right)=2(I+J)=2({r}_{i}J+J)=2(1+{r}_{i})J$ . The escaping flux is the observable reflectance $R=(1-{r}_{i})J$ . Hence, $J=R\u2215(1-{r}_{i})$ . Therefore, the fluence rate at the surface is

This is the first key equation. Next, consider Fick’s first law of diffusion. The escaping flux is $R=-D(\partial \varphi \u2215\partial z)$ evaluated at $z\approx 0$ . 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 $z={z}_{b}$ , where ${z}_{b}$ 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 $1\text{-}\mathrm{W}$ point source of light is located at a depth ${z}_{s}$ within the tissue, then a negative point source of equal power is placed at a position ${z}_{i}=2{z}_{b}-{z}_{s}$ .The example situation illustrated in Fig. 11b uses 1-D planar diffusion theory ${T}_{\text{planar}}\left(z\right)=\mathrm{exp}(-z\u2215\delta )\u2215\left(2{\mu}_{\alpha}\delta \right)$ to better illustrate the gradient of $\varphi $ between a real source at ${z}_{s}$ and a negative image source at ${z}_{i}$ . When 3-D diffusion theory is used to consider the response to a point source that is close to the surface, the $\varphi $ 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 $\left(1\phantom{\rule{0.3em}{0ex}}{\mathit{MFP}}^{\prime}\right)$ from the point source, and the boundary condition discussed here is applicable to 3-D diffusion from a point source at ${z}_{s}$ .

The choice of ${z}_{b}$ is such that the gradient from the tissue surface to this virtual surface of symmetry is $-\partial \varphi \u2215\partial z=\varphi \left(0\right)\u2215{z}_{b}$ . Therefore, the escaping flux is equal to

This is the second key equation. Combining Eqs. 29, 30 and eliminating $R$ yieldsEquation 30 describes the boundary condition. The value of internal reflectance ${r}_{i}$ specifies the value ${z}_{b}$ . The value ${r}_{i}$ can be calculated:## Eq. 32

$${r}_{i}={\int}_{\theta =0}^{\pi \u22152}{R}_{\text{Fresnel}}\left(\theta \right)2\pi \phantom{\rule{0.2em}{0ex}}\mathrm{sin}\left(\theta \right)d\theta ,$$^{9}as cited by Groenhuis Ferwerda, and Ten Bosch.

^{10}It is very good for $1.0<n<2.2$ .

Therefore, for a generic air/tissue interface, the value of ${r}_{i}$ is $\sim 0.5$ , and the value of $(1+{r}_{i})\u2215(1-{r}_{i})$ is $\sim 3$ . Therefore, ${z}_{b}\approx -\left(2\right)\left(3\right)({z}_{s}\u22153)\approx -2{z}_{s}$ . Hence, the position of the negative image source is placed at ${z}_{i}\approx -5{z}_{s}$ .

Once the positions ${z}_{s}$ and ${z}_{i}$ have been specified, the superposition of the two sources yields the net fluence rate in the tissue:

## Eq. 34

$$\varphi (r,z)=\frac{\mathrm{exp}(-{r}_{1}\u2215\delta )}{4\pi D{r}_{1}}-\frac{\mathrm{exp}(-{r}_{2}\u2215\delta )}{4\pi D{r}_{2}},$$## Eq. 35

$$R\left(r\right)={\phantom{\mid}-D\frac{\partial \varphi}{\partial z}\mid}_{z=0}=\frac{1}{4\pi}{z}_{s}(\frac{1}{\delta}+\frac{1}{{r}_{1}})\frac{\mathrm{exp}(-\mid {r}_{1}\mid \u2215\delta )}{4\pi D{r}_{1}^{2}}+({z}_{s}+2{z}_{b})(\frac{1}{\delta}+\frac{1}{{r}_{2}})\frac{\mathrm{exp}(-\mid {r}_{2}\mid \u2215\delta )}{4\pi D{r}_{2}^{2}}.$$^{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 ${z}_{s}=1\phantom{\rule{0.3em}{0ex}}{\mathit{MFP}}^{\prime}$ . Equations 34, 35 yield good approximations to the fluence rate $\varphi (z,r)$ and escaping flux $R\left(r\right)$ at locations greater than $1\phantom{\rule{0.3em}{0ex}}{\mathrm{MPF}}^{\prime}$ distant from the origin at $(r=0)$ and greater than $1\phantom{\rule{0.3em}{0ex}}{\mathrm{MPF}}^{\prime}$ from the source at ( $r=0$ , $z={z}_{s}$ ). Monte Carlo simulations are recommended for regions close to the source.

## 3.

## Diffusion Implementation in Numerical Methods

## 3.1.

### 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:

## Eq. 36

$$\frac{1}{c}\frac{\partial \varphi (r,t)}{\partial t}-\nabla \cdot D\nabla \varphi (\mathbf{r},t)+{\mu}_{a}\varphi (\mathbf{r},t)=S(\mathbf{r},t),$$## Eq. 37

$$-\nabla \cdot D\nabla \varphi (\mathbf{r},\omega )+({\mu}_{a}+\frac{i\omega}{c})\varphi (\mathbf{r},\omega )=S(\mathbf{r},\omega ).$$## Eq. 38

$$-\nabla \cdot D\nabla \varphi \left(\mathbf{r}\right)+{\mu}_{a}\varphi \left(\mathbf{r}\right)=S\left(\mathbf{r}\right).$$^{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.

## Table 4

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.

Method | Mesh/grid | Memory | Time | Advantages | Disadvantages |
---|---|---|---|---|---|

Integratedanalyticmethod | Arbitrary | Small memory | $N$ | Fast,conceptually theeasiest approach | Integration over highlyshaped regions andheterogeneties are difficult |

Finitedifferencemethod | Cartesian | 2-D: ${B}^{*}N$ 3-D: $N$ | 2-D: ${B}^{2*}N$ 3-D: $N$ | Ease of constructionand discretization,sparse/bandedmatrices | Difficult to handleirregular boundaries |

Finite elementmethod | Arbitrary(triangular elementsin 2-D; tetrahedronsin 3-D) | 2-D: ${B}^{*}N$ 3-D: $N$ | 2-D: ${B}^{2*}N$ 3-D: $N$ | Arbitrary shapes,flexible,sparse/bandedmatrices | Volume meshing in 3-Dcan be difficult forarbitrary shapes |

Boundaryelementmethod | Arbitrary(line segments in2-D; triangularelements in 3-D) | ${N}^{2}$ | ${N}^{3}$ | Arbitrary shapes,nodes only onboundary | Full matrices;assumption of piece-wise constant regions |

## 3.2.

### Numerical Methods to Solve the Diffusion Equation

## 3.2.1.

#### 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
$i,j,k$
for the
$x,y,z$
coordinates, and the fluence at each location calculated by a straightforward approximation of the derivative by a finite difference, such that it becomes:

## Eq. 39

$${D}_{i,j,k}(\frac{{\varphi}_{i+1,j,k}-{\varphi}_{i-1,j,k}}{\Delta {L}_{i}^{2}}+\frac{{\varphi}_{i,j+1,k}-{\varphi}_{i,j-1,k}}{\Delta {L}_{j}^{2}}+\frac{{\varphi}_{i,j,k+1}-{\varphi}_{i,j,k-1}}{\Delta {L}_{k}^{2}})+[\frac{({D}_{i+1,j,k}-{D}_{i,j,k})({\varphi}_{i+1,j,k}-{\varphi}_{i,j,k})}{\Delta {L}_{i}^{2}}+\frac{({D}_{i,j+1,k}-{D}_{i,j,k})({\varphi}_{i,j+1,k}-{\varphi}_{i,j,k})}{\Delta {L}_{j}^{2}}+\frac{({D}_{i,j,k+1}-{D}_{i,j,k})({\varphi}_{i,j,k+1}-{\varphi}_{i,j,k})}{\Delta {L}_{k}^{2}}]+{\mu}_{a}{\varphi}_{i,j,k}={S}_{ijk}.$$^{7}and that the optical properties $D$ and ${\mu}_{a}$ 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}

## 3.2.2.

#### 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
${\phi}_{i}$
, 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 $\varphi (r,\omega )$ can be discretized onto the basis functions multiplied by weighting coefficients $\varphi ={\sum}_{j=1}^{N}{\Phi}_{j}{\varphi}_{j}$ , where $\Phi $ is the discrete fluence at each node $j$ with weights $\phi $ , 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 ${\varphi}_{i}$ and integrated over the entire problem domain to give:

## Eq. 40

$$\u27e8-\nabla \cdot D\nabla \Phi {\phi}_{i}\u27e9+\u27e8({\mu}_{a}+\frac{i\omega}{c})\Phi {\phi}_{i}\u27e9=\u27e8S{\phi}_{i}\u27e9,$$## Eq. 41

$$\u27e8D\nabla \Phi \cdot \nabla {\phi}_{i}\u27e9-\oint \widehat{s}\cdot D\nabla \Phi {\phi}_{i}ds+\u27e8({\mu}_{a}+\frac{i\omega}{c})\Phi {\phi}_{i}\u27e9=\u27e8S{\phi}_{i}\u27e9.$$## Eq. 43

$$\alpha =\frac{\widehat{s}\cdot D\nabla \varphi}{\varphi}=\frac{D}{\varphi}{\phantom{\mid}\frac{\partial \varphi}{\partial r}\mid}_{r=R},$$^{22}While analytical forms of $\alpha $ 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 $(D=\sum {D}_{i}{\phi}_{i})$ and the absorption coefficient $({\mu}_{a}=\sum {\mu}_{ai}{\phi}_{i})$ in the same manner as was done with the fluence, then the equation becomes:

## Eq. 44

$$\u27e8D\Phi {\phi}_{k}\nabla {\phi}_{j}\cdot \nabla {\phi}_{i}\u27e9+\u27e8{\mu}_{a}\Phi {\phi}_{l}{\phi}_{i}\u27e9=\u27e8S{\phi}_{i}\u27e9+\alpha \Phi {\phi}_{i}.$$## Eq. 45

$${a}_{ij}=\u27e8-\sum {D}_{i}{\phi}_{k}\nabla {\phi}_{j}\nabla {\phi}_{i}-\sum {\mu}_{ai}{\phi}_{i}{\phi}_{l}\u27e9,$$^{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:

where $\mathbf{\Phi}$ is the vector of fluence values at all nodes $N$ and $\mathbf{S}$ is the vector of sources at each node, driving the predicted fluence. This fluence is then solved simply by inverting the connectivity matrix $\mathbf{A}$ :The diffusion matrix $\mathbf{A}$ 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
$\left[\mathbf{A}\right]$
$\mathbf{\Phi}=\mathbf{S}$
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 $\mathbf{A}$ matrix and $\mathbf{S}$ vector, the result from inversion of $\mathbf{A}$ is a discretized map of the optical fluence pattern $\mathbf{\Phi}$ , 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.

## 3.2.3.

#### 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 prediction^{30} 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:

## Eq. 48

$${G}_{i}(r,{r}_{i})=\frac{\mathrm{exp}\left(\frac{-{k}_{l}\mid r-{r}_{i}\mid}{\sqrt{{D}_{l}}}\right)}{4\pi {D}_{l}\mid r-{r}_{i}\mid},\phantom{\rule{1em}{0ex}}\text{in}\phantom{\rule{0.3em}{0ex}}3\text{-}\mathrm{D},$$The final version of this expansion for the diffusion equation leads to:^{31}

## Eq. 49

$$\left[\mathbf{A}\right]\mathbf{\Phi}-\left[\mathbf{B}\right]\left\{\mathbf{D}(\partial \Phi \u2215\partial n)\right\}=\mathbf{S},$$## 50.

## 4.

## Applications of Diffusion Modeling of Light Transport

## 4.1.

### 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 ${y}_{d}$ , which is a function of the fluence at the boundary $\varphi $ . This can be simulated by a diffusion theory calculation ${y}_{d}=f\left[\varphi ({\mu}_{a},{\mu}_{s}^{\u2215})\right]$ , 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.

^{32, 33}However, inversion with models that are not for infinite media do not lead to an explicit inversion formula.

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:

## Eq. 51

$$\text{minimize}\phantom{\rule{0.3em}{0ex}}f({\mu}_{a},{\mu}_{s}^{\u2215})={\{y-{y}_{d}({\mu}_{a},{\mu}_{s}^{\u2215})\}}^{2}.$$## Eq. 53

$${j}_{ik}=\{\partial {y}_{di}\u2215d{\mu}_{ak};\partial {y}_{di}\u2215d{\mu}_{sk}^{\u2215}\}.$$## Eq. 54

$$0=\{y-{y}_{d}({\mu}_{a0},{\mu}_{s0}^{\u2215})\}+\mathbf{J}\{\Delta {\mu}_{a},\Delta {\mu}_{s}^{\u2215}\}+\cdots \phantom{\rule{0.2em}{0ex}},$$## Eq. 55

$$\mathbf{J}\{\Delta {\mu}_{a},\Delta {\mu}_{s}^{\u2215}\}=\{y-{y}_{d}({\mu}_{ai},{\mu}_{si}^{\u2215})\}.$$## Eq. 56

$$\{\Delta {\mu}_{a,I},\Delta {\mu}_{si}^{\u2215}\}={\left[{\mathbf{J}}^{T}\mathbf{J}\right]}^{-1}{\mathbf{J}}^{T}\{y-{y}_{d}({\mu}_{ai},{\mu}_{si}^{\u2215})\}.$$## Eq. 57

$$\{\Delta {\mu}_{a,I},\Delta {\mu}_{si}^{\u2215}\}={[{\mathbf{J}}^{T}\mathbf{J}+\lambda \mathbf{I}]}^{-1}\left[{\mathbf{J}}^{T}\right]\{y-{y}_{d}({\mu}_{ai},{\mu}_{si}^{\u2215})\}.$$Recovery of tissue values at a single wavelength are readily achieved, and computer codes such as TOAST (University College London^{34}) and NIRFAST (Dartmouth College and University of Exeter^{35}) are freely available for this problem.

## 4.2.

### Absorption Spectroscopy

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:

## Eq. 59

$$\mathbf{C}={\left[{\mathbf{\epsilon}}^{T}\mathbf{\epsilon}\right]}^{-1}{\epsilon}^{T}{\mu}_{a}.$$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.

## 4.3.

### 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:

Since this is a nonlinear equation, it is again not easily inverted for $a$ or $c$ , and ultimately these parameters are not intuitive, as the values for $a$ and $b$ are arbitrary from an analytic approximation to the Mie theory. The important spectra to fit comes from the Mie theory:## Eq. 61

$${\mu}_{s}^{\prime}\left(\lambda \right)={N}_{o}\sum _{i=1}^{p}f\left({a}_{i}\right)\left(\pi {a}_{i}^{2}\right){Q}_{\mathrm{scat}}(m,{a}_{i},\lambda )[1-g(m,{a}_{i},\lambda )],$$^{43}and calculated anisotropy parameter $g$ .

^{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}

## 4.4.

### 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.

## 5.

## Conclusion

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.

## References

**,” Comput. Methods Programs Biomed., 47 (2), 131 –146 (1995). https://doi.org/10.1016/0169-2607(95)01640-F 0169-2607 Google Scholar**

*MCML-Monte-Carlo modeling of light transport in multilayered tissues***,” J. Opt. Soc. Am. A Opt. Image Sci. Vis, 10 (8), 1746 –1752 (1993). 1084-7529 Google Scholar**

*Hybrid model of Monte-Carlo simulation and diffusion-theory for light reflectance by turbid media***,” Opt. Eng., 32 (2), 258 –266 (1993). https://doi.org/10.1117/1.2170990 0091-3286 Google Scholar**

*Tissue characterization and imaging using photon density waves***,” J. Opt. Soc. Am. A Opt. Image Sci. Vis, 14 (1), 255 –261 (1997). 1084-7529 Google Scholar**

*Perturbation theory for diffuse light transport in complex biological tissues***,” Gastrointest. Endosc., 41 (3), 218 –224 (1995). 0016-5107 Google Scholar**

*Infrared video imaging of subsurface vessels—a feasibility study for the endoscopic management of gastrointestinal-bleeding***,” J. Biomed. Opt., 5 (3), 277 –282 (2000). https://doi.org/10.1117/1.429996 1083-3668 Google Scholar**

*Modeling photon transport in transabdominal fetal oximetry***,” Appl. Opt., 22 (16), 2456 –2462 (1983). 0003-6935 Google Scholar**

*Scattering and absorption of turbid materials determined from reflectaion measurements. 1: Theory***,” Appl. Opt., 28 (12), 2331 –2336 (1989). 0003-6935 Google Scholar**

*Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical-properties***,” Phys. Med. Biol., 37 (7), 1531 –1560 (1992). https://doi.org/10.1088/0031-9155/37/7/005 0031-9155 Google Scholar**

*The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis***,” Phys. Med. Biol., 39 1157 –1180 (1994). https://doi.org/10.1088/0031-9155/39/7/008 0031-9155 Google Scholar**

*Frequency domain optical absorption spectroscopy of finite tissue volumes using diffusion theory***,” Med. Phys., 22 (6), 691 –701 (1995). https://doi.org/10.1118/1.597488 0094-2405 Google Scholar**

*Spatially varying optical property reconstruction using a finite element diffusion equation approximation***,” Med. Phys., 20 (2), 299 –309 (1993). https://doi.org/10.1118/1.597069 0094-2405 Google Scholar**

*A finite-element approach for modeling photon transport in tissue***,” Phys. Med. Biol., 40 (10), 1709 –1729 (1995). https://doi.org/10.1088/0031-9155/40/10/011 0031-9155 Google Scholar**

*Initial assessment of a simple system for frequency-domain diffuse optical tomography***,” J. Appl. Math. Comput., 34 I13 –46 (1989). Google Scholar**

*MUDPACK: multigrid portable fortran software for the efficient solution of linear elliptic differential equations***,” Phys. Med. Biol., 40 1709 –1729 (1995). https://doi.org/10.1088/0031-9155/40/10/011 0031-9155 Google Scholar**

*Initial assessment of a simple system for frequency domain diffuse optical tomography***,” Appl. Opt., 36 2260 –2272 (1997). 0003-6935 Google Scholar**

*Imaging of fluorescent yield and lifetime from multiply scattered light re-emitted from tissues and other random media***,” Appl. Opt., 40 (4), 588 –600 (2001). https://doi.org/10.1364/AO.40.000588 0003-6935 Google Scholar**

*Three-dimensional simulation of near-infrared diffusion in tissue: boundary condition and geometry analysis for finite element image reconstruction***,” Inverse Probl., 15 (2), R41 –R93 (1999). https://doi.org/10.1088/0266-5611/15/2/022 0266-5611 Google Scholar**

*Optical tomography in medical imaging***,” J. Opt. Soc. Am. A, 11 (10), 2727 –2741 (1994). 0740-3232 Google Scholar**

*Boundary conditions for the diffusion equation in radiative transfer***,” J. Opt. Soc. Am. A, 12 (11), 2532 –2539 (1995). 0740-3232 Google Scholar**

*Boundary-conditions for diffusion of light***,” Phys. Med. Biol., 40 1957 –1975 (1995). https://doi.org/10.1088/0031-9155/40/11/013 0031-9155 Google Scholar**

*Influence of boundary conditions on the accuracy of diffusion theory in time-resolved reflectance spectroscopy of biological tissues***,” Appl. Opt., 28 2331 –2336 (1989). 0003-6935 Google Scholar**

*Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties***,” IEEE Trans. Med. Imaging, 25 (9), 1180 –1188 (2006). https://doi.org/10.1109/TMI.2006.879957 0278-0062 Google Scholar**

*Electrode boundary conditions and experimental validation for BEM-based EIT forward and inverse solutions***,” Phys. Med. Biol., 51 497 –516 (2006). https://doi.org/10.1088/0031-9155/51/3/003 0031-9155 Google Scholar**

*Diffuse photon propagation in multilayered geometries***,” J. Comp. Physiol., 210 109 –132 (2005). 0373-0859 Google Scholar**

*Fluorescence photon migration by the boundary element method***,” Med. Phys., 0094-2405 Google Scholar**

*A boundary element approach for image-guided near-infrared absorption and scatter estimation***,” J. Opt. Soc. Am. B, 11 (10), 2128 –2138 (1994). 0740-3224 Google Scholar**

*Semi-infinite-geometry boundary-problem for light migration in highly scattering media—a frequency-domain study in the diffusion-approximation***,” J. Biomed. Opt., 5 (2), 182 –93 (2000). 1083-3668 Google Scholar**

*Calibration of near infrared frequency-domain tissue spectroscopy for absolute absorption coefficient quantitation in neonatal head-simulating phantoms***,” Breast Cancer Res. Treat., 7 (6), 279 –285 (2005). 0167-6806 Google Scholar**

*Imaging in breast cancer: diffuse optics in breast cancer: detecting tumors in pre-menopausal women and monitoring neoadjuvant chemotherapy***,” Neuroimage, 23 S275 –S288 (2004). 1053-8119 Google Scholar**

*Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy***,” Philos. Trans. R. Soc. London, Ser. B, 352 (1354), 727 –735 (1997). 0962-8436 Google Scholar**

*Measurements of scattering and absorption changes in muscle and brain***,” J. Biomed. Opt., 5 (1), 102 –105 (2000). https://doi.org/10.1117/1.429975 1083-3668 Google Scholar**

*Quantification of ischemic muscle deoxygenation by near infrared time-resolved spectroscopy***,” Appl. Opt., 44 (10), 1948 –1956 (2005). https://doi.org/10.1364/AO.44.001948 0003-6935 Google Scholar**

*Optimal linear inverse solution with multiple priors in diffuse optical tomography***,” Acad. Radiol., 1076-6332 Google Scholar**

*In vivo hemoglobin and water concentration, oxygen saturation, and scattering estimates from near-infrared breast tomography using spectral reconstruction***,” Proc. Natl. Acad. Sci. U.S.A., 103 (23), 8828 –8833 (2006). https://doi.org/10.1073/pnas.0509636103 0027-8424 Google Scholar**

*Imaging breast adipose and fibroglandular tissue molecular signatures using hybrid MRI-guided near-infrared spectral tomography***,” Appl. Opt., 30 (31), 4507 –4514 (1991). 0003-6935 Google Scholar**

*Light scattering in intralipid-10% in the wavelength range of $400\u20131100\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$***,” Nature (London), 406 (6791), 35 –36 (2000). https://doi.org/10.1038/35017638 0028-0836 Google Scholar**

*Detection of preinvasive cancer cells***,” Nat. Med., 7 (11), 1245 –1248 (2001). 1078-8956 Google Scholar**

*Imaging human epithelial properties with polarized light-scattering spectroscopy***,” Cancer Res., 63 (13), 3556 –3559 (2003). 0008-5472 Google Scholar**

*In situ detection of neoplastic transformation and chemopreventive effects in rat esophagus epithelium using angle-resolved low-coherence interferometry***,” Gastrointest. Endosc., 63 (2), 257 –261 (2006). 0016-5107 Google Scholar**

*Elastic scattering spectroscopy for the diagnosis of colonic lesions: initial results of a novel optical biopsy technique***,” J. Biomed. Opt., 11 (4), 041106 (2006). https://doi.org/10.1117/1.2342747 1083-3668 Google Scholar**

*Image reconstruction of effective mie scattering parameters of breast tissue in vivo with near-infrared tomography***,” Opt. Lett., 32 (8), 933 –935 (2007). https://doi.org/10.1364/OL.32.000933 0146-9592 Google Scholar**

*Image-guided optical spectroscopy provides molecular-specific information in vivo: MRI-guided spectroscopy of breast cancer hemoglobin, water, and scatterer size***,” Cytom. A., 69 (8), 748 –758 (2006). Google Scholar**

*Spectral imaging in biology and medicine: slices of life***,” Opt. Lett., 31 (3), 365 –367 (2006). https://doi.org/10.1364/OL.31.000365 0146-9592 Google Scholar**

*Spectrally-resolved bioluminescence optical tomography***,” IEEE Trans. Med. Imaging, 24 (7), 878 –885 (2005). https://doi.org/10.1109/TMI.2004.843254 0278-0062 Google Scholar**

*Fluorescent protein tomography scanner for small animal imaging***,” Proc. Natl. Acad. Sci. U.S.A., 101 (33), 12294 –12299 (2004). https://doi.org/10.1073/pnas.0401137101 0027-8424 Google Scholar**

*Visualization of antitumor treatment by means of fluorescence molecular tomography with an annexin V-Cy5.5 conjugate***,” Nat. Med., 8 (7), 757 –760 (2002). 1078-8956 Google Scholar**

*Fluorescence molecular tomography resolves protease activity in vivo***,” Mol. Imaging, 1 (2), 82 –88 (2002). https://doi.org/10.1162/153535002320162732 1535-3508 Google Scholar**

*In vivo tomographic imaging of near-infrared fluorescent probes***,” J. Biomed. Opt., 10 064007 (2005). https://doi.org/10.1117/1.2136148 1083-3668 Google Scholar**

*Planar fluorescence imaging using normalized data***,” Thayer School of Engineering, 195 (2005) Google Scholar**

*MRI-coupled broadband near-infrared tomography for small animal brain studies***,” J. Biomed. Opt., 10 (5), 051704-1--10 (2005). https://doi.org/10.1117/1.2098607 1083-3668 Google Scholar**

*Approximation of Mie scattering parameters from near-infrared tomography of health breast tissue in vivo*