2.1.

Optical Properties and Transport Parameters

This tutorial follows the notation of Welch and van Gemert.¹ 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}∕{\mathrm{cm}}^2$ and $T$ is dimensionless.

Table 1

Tissue optical properties.

Table 2

Transport properties used in diffusion theory.

Table 3

Parameters for light transport.

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.² 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∕[{\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.

Fig. 1

Monte Carlo simulation of photons launched as collimated beams from optical fibers. The light appears to diffuse from a central point located one ${\mathit{MFP}}^{\prime }=1∕[{\mu }_a+{\mu }_s(1-g)]$ in front of the delivery fiber (red dot at $r=0$ , $z=0.3\mathrm{cm}$ ). Iso-concentration contours based on Monte Carlo are drawn as black dots and black line. Dashed red circular lines centered around the center of diffusion indicate the prediction of diffusion theory. Near the source, diffusion theory does not agree with the Monte Carlo simulation, but distant from the source agreement is good. The mismatched air/medium boundary at $z=0$ slightly distorts the iso-concentration curves of the Monte Carlo simulation, so Monte Carlo and dashed lines of diffusion theory do not exactly match. (Optical properties: ${\mu }_a=0.1{\mathrm{cm}}^{-1}$ , ${\mu }_s=100{\mathrm{cm}}^{-1}$ , $g=0.90$ , $n=1.4$ .) (Color online only.)

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.³ 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}∕{\mathrm{cm}}^2$ per W], which is the transport $T_{\mathit{MC}}(x,y,z)$ $[1∕{\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\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)∕{\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)∕{\mu }_a=10$ is 0.252 (based on Monte Carlo simulations), where the mismatched tissue/air surface boundary has $n_{\text{tissue}}∕n_{\text{air}}=1.4$ . A study of reflectance versus ${\mu }_s(1-g)∕{\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)∕{\mu }_a\text{equal}\sim 20]$ .⁴ So the reflectance of a tissue should exceed at least 25% and preferably 40%, or ${\mu }_s(1-g)∕{\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)∕{\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)∕{\mu }_a⪡10$ . 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}∕{\mathrm{cm}}^2]$ is proportional to the concentration of optical energy $C$ $[\mathrm{J}∕{\mathrm{cm}}^3]$ :

Eq. 1

$$\varphi =cC,$$

where $c$ is the speed of light in the medium, $c=c_o∕n$ [cm/s], $c_o=2.998\times {10}^{10}\mathrm{cm}∕\mathrm{s}$ , and $n$ is the refractive index of the medium. For example, consider a uniform beam of $1\mathrm{W}∕{\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\mathrm{W}∕{\mathrm{cm}}^2$ at $500\text{-}\mathrm{nm}$ wavelength, then the corresponding concentration of photons is $0.18\mathrm{pM}$ $(0.18\times {10}^{-12}\text{moles}∕\text{liter})$ :

Eq. 2

$$\varphi \frac{n}{c_0}\frac{\lambda }{h{c}_0}\frac{1000}{N_{\mathrm{Av}}}=(1\mathrm{W}∕{\mathrm{cm}}^2)\frac{1.33}{{(3\times {10}^{10}\mathrm{cm}∕\mathrm{s})}^2}\frac{(500\times {10}^{-7}\mathrm{cm})}{(6.63\times {10}^{-34}Js)}\frac{1000{\mathrm{cm}}^3∕\text{liter}}{(6.02\times {10}^{23}{\text{mole}}^{-1})}=0.18\times {10}^{-12}\mathrm{M}.$$

The factor $\lambda ∕hc$ , where $\lambda$ is wavelength, $h$ is Planck’s constant, and $c$ is the in vacuo speed of light, indicates the number of photons per J energy. The factor 1000 indicates the conversion $1000{\mathrm{cm}}^3∕\text{liter}$ , and $N_{\mathrm{Av}}$ is Avagadro’s number for the molecules per mole.

Fig. 2

The concentration of photons in a cube of water irradiated with $1\mathrm{W}∕{\mathrm{cm}}^2$ is $0.18\mathrm{pM}$ .

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{\exp [-r^2∕\left(4\chi t\right)]}{{\left(4\pi \chi t\right)}^{3∕2}}.$$

The $Q$ 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 $(r=0)$ at time zero $(t=0)$ . The second factor has units of ${\mathrm{cm}}^{-3}$ and is referred to as Green’s function. Hence, the concentration $C$ at a distance $r$ from the source at a time $t$ after the initial deposition of $Q$ is expressed as $[\text{units}∕{\mathrm{cm}}^3]$ .

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 $\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{\exp [-r^2∕\left(4cDt\right)]}{{\left(4\pi cDt\right)}^{3∕2}}\exp (-{\mu }_{a}ct).$$

This is the time-resolved 3-D diffusion equation for light in a medium with scattering and absorption. Figure 3 gives an example of time-resolved diffusion of light from an impulse.

Fig. 3

Time-resolved diffusion of light from an impulse of light, $Q=1\mathrm{J}$ , delivered as a point source at $r=0$ , $t=0$ . The time-resolved fluence rate $\varphi (r,t)$ is shown at $t=0.1$ , 0.2, and $0.5\mathrm{ns}$ .

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}∕{\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 $\exp (-r∕\text{constant})$ . The constant is given the symbol $\delta$ [cm] and is called the optical penetration depth. The factor $1∕\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}∕{\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{\exp (-r∕\delta )}{\text{factor}}4\pi r^{2}dr=Q,$$

where $dV$ is expressed as an incremental shell volume $4\pi r^{2}dr$ . The denominator called factor needs to be specified. The factor needs the term $1∕Q$ to yield the final value $Q$ . The factor needs ${\mu }_a$ to cancel the ${\mu }_a$ in the numerator. The factor needs $4\pi$ to cancel the $4\pi$ in the numerator. The factor needs $r$ to reduce the $r^2$ in the numerator to simply $r$ , so that $\exp (-r∕\delta )rdr$ can be integrated by parts. A factor ${\delta }^2$ is generated by this integration, so the factor needs a ${\delta }^2$ to cancel this term. In summary, the factor must equal $4\pi {\mu }_a{\delta }^{2}r∕Q$ to allow conservation of energy in Eq. 7. Therefore, the expression for the radiant exposure $\psi \left(r\right)$ in response to an impulse of energy $Q$ at $r=0$ and $t=0$ is

Eq. 8

$$\psi \left(r\right)={\int }_{t=0}^{\infty }\varphi (r,t)dt=Q\frac{\exp (-r∕\delta )}{4\pi {\mu }_a{\delta }^{2}r}=Q\frac{\exp (-r∕\delta )}{4\pi Dr}.$$

The factor ${\mu }_a{\delta }^2$ is identical to $D$ (derivation not shown here), and the expression using $D$ is also shown in Eq. 8 and is commonly used.

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∕\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$ :

Eq. 9

$$\varphi \left(r\right)=P\frac{\exp (-r∕\delta )}{4\pi Dr}.$$

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.

Fig. 4

Steady-state diffusion of light from a point source $P=1\mathrm{W}$ at $r=0$ . The dashed line indicates diffusion with no absorption. The data near $r=0$ are not correct because diffusion theory is inaccurate near the source where the gradient of fluence rate has significant curvature $(\mid {\partial }^2\varphi ∕\partial r^2\mid ⪢0)$ . (Top) Fluence rate, $\varphi$ . (Bottom) Residual error expressed as $\left(\mathit{MC}\text{-}\mathit{DT}\right)∕DT$ , where $\mathit{MC}$ is the Monte Carlo result and $\mathit{DT}$ is the diffusion theory result.

The transport factor $T$ equals $\varphi ∕P$ 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}∕{\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:

10.

$$T_{\text{planar}}\left(z\right)=\frac{\exp (-z∕\delta )}{2{\mu }_a\delta },$$

$$T_{\mathrm{cyl}}\left(r\right)\approx \sqrt{\frac{2\delta }{\pi r}}\frac{\exp (-r∕\delta )}{2\pi {\mu }_a{\delta }^2},$$

$$T_{\mathrm{sph}}\left(r\right)=\frac{\exp (-r∕\delta )}{4\pi {\mu }_a{\delta }^{2}r}.$$

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\sin \left(\omega t\right)],$$

where $M_o$ is the modulation of the source, $0<M_o<1$ [dimensionless], and $\omega$ is the angular frequency of modulation [radians/s], $\omega =2\pi f$ , 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) is

12.

$$\varphi (r,t)=P_0[T_{ss}\left(r\right)+M_o\frac{\exp (-r{k}^{″})}{4\pi Dr}\sin (\omega t-k^{\prime }r)],$$

$$=P_0{T}_{ss}\left(r\right)[1+M_o\exp [-r(k^{″}-1∕\delta )]\sin (\omega t-k^{\prime }r)],$$

where $T_{ss}\left(r\right)$ $[1∕{\mathrm{cm}}^2]$ is the steady-state transport:

Eq. 13

$$T_{ss}\left(r\right)=\frac{\exp (-r∕\delta )}{4\pi Dr},$$

the factors $k\prime$ and $k″$ are defined:

14.

$$k^{″}=\frac{1}{\delta \sqrt{2}}{\{{[1+{\left(\frac{\omega }{{\mu }_{a}c}\right)}^2]}^{1∕2}+1\}}^{1∕2}=\frac{a}{\delta }\cos \left(b\right),$$

$$k^{\prime }=\frac{1}{\delta \sqrt{2}}{\{{[1+{\left(\frac{\omega }{{\mu }_{a}c}\right)}^2]}^{1∕2}-1\}}^{1∕2}=\frac{a}{\delta }\sin \left(b\right),$$

where

15.

$$a={[1+{\left(\frac{\omega }{{\mu }_{a}c}\right)}^2]}^{1∕4},$$

$$b=\frac{1}{2}\arctan \left(\frac{\omega }{{\mu }_{a}c}\right).$$

These above expressions are consistent with the appendix of Svaasand ⁵ 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.

$$\text{Modulation}M=M_o\exp [-r(k^{″}-1∕\delta )],$$

$$\text{Phase}\Phi =k^{\prime }r.$$

If the angular frequency of modulation $\omega$ is swept over a broad range, then the values of $M$ and $\Phi$ versus frequency constitute a frequency domain description of light transport.

For very low frequencies of modulation as $\omega \to 0$ , the factor $a$ approaches 1 and the factor $b$ approaches 0. Therefore, $k^{″}$ approaches $1∕\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\sin \left(\omega t\right)]=P\left(t\right)T_{ss}\left(r\right),$$

which is simply the steady-state response to the slowly varying light source $P\left(t\right)$ .

For very high frequencies of modulation as $\omega ⪢{\mu }_{a}c$ , the factor $a$ approaches ${[\omega ∕\left({\mu }_{a}c\right)]}^{1∕2}$ , and the factor $b$ approaches $\pi ∕4$ . Since $\cos (\pi ∕4)=\sin (\pi ∕4)=1∕\sqrt{2}$ , the values $k^{″}$ and $k^{\prime }$ approach the same limiting value:

Eq. 18

$$k^{″}=k^{\prime }=\frac{1}{\delta \sqrt{2}}{\left(\frac{\omega }{{\mu }_{a}c}\right)}^{1∕2}.$$

In the intermediate regime where $\omega$ is a little below or a little above ${\mu }_{a}c$ , the behavior of $M$ and $\Phi$ varies versus $\omega$ , each in a unique manner, and therefore $M$ and $\Phi$ observed at a particular distance $r$ from the source encode the optical properties ${\mu }_a$ and ${\mu }_s(1-g)$ . At very low $\omega ⪡{\mu }_{a}c$ , the phase $\Phi$ is too low to reliably measure. At high $\omega >{\mu }_{a}c$ , the $\Phi$ becomes insensitive to ${\mu }_a$ . So the use of modulated light is typically in the intermediate region where $\omega$ is a little less than ${\mu }_{a}c$ . For a tissue with $n=1.37$ and ${\mu }_a=0.1{\mathrm{cm}}^{-1}$ , the value $\omega ={\mu }_a{c}_0∕n$ is $1.65\times {10}^{-9}[\mathrm{rad}∕\mathrm{s}]$ , and the frequency $f=\omega ∕\left(2\pi \right)$ is $262\mathrm{MHz}$ . Typically, modulated diffusion measurements are made with frequencies in the range $f=25\text{to}1000\mathrm{MHz}$ , which correspond to ${\mu }_a=\omega n∕c_0=0.007\text{to}0.29{\mathrm{cm}}^{-1}$ , respectively. Figure 5 shows $M$ and $\Phi$ versus $r$ for a series of modulation frequencies. Figure 6 shows $M$ and $\Phi$ versus modulation frequency for a series of $r$ positions.

Fig. 5

Modulation ( $M$ [dimensionless]) and phase ( $\Phi$ [radians]) versus distance ( $r$ [cm]) of observation point from source, for various frequencies ( $f$ [Hz]) of full source modulation, $M_o=1$ . Optical properties of medium are absorption ${\mu }_a=0.1{\mathrm{cm}}^{-1}$ , and reduced scattering ${\mu }_s(1-g)=10{\mathrm{cm}}^{-1}$ .

Fig. 6

Modulation $\left(M\right)$ and phase $\left(\Phi \right)$ versus frequency of modulation, observed at $r=1$ , 2, and $3\mathrm{cm}$ from source. The source modulation and tissue properties are the same as in Fig. 5.

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

Fig. 7

A $0.2\text{-}\mathrm{cm}$ -diam. absorbing object perturbs the field of light from a source, as seen by a detector. The transport of light from the source to the perturbing object $T_{s\to p}$ is absorbed by the volume $\left(V\right)$ of the sphere, which has an incremental absorption $\Delta {\mu }_a$ above the background tissue optical properties, but the same scattering properties as the background tissue. The power absorbed by the sphere, $\Delta {\mu }_a{T}_{s\to p}V$ , has units of power [W], and acts as a negative optical power, $P_p=-\Delta {\mu }_a{T}_{s\to p}V$ , that propagates throughout the tissue and subtracts from the background light supplied by the source. The detector sees a local fluence rate ${\varphi }_d=P{T}_{s\to d}+P_p{T}_{p\to d}$ . Hence, the signal at the detector slightly drops in the presence of the perturbing object.

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∕2}$ . The extra energy deposition in the absorbing object, above the background absorption, is $\Delta {\mu }_a{\varphi }_{\mathrm{obj}}$ $[\mathrm{W}∕{\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 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$ :

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.

Fig. 8

The virtual sources for perturbation by an object. (a) Subvolumes within an object, accounting for object absorption and scattering-enhanced absorption. (b) Surface sources on the object’s surface, creating a dipole oriented along the gradient of fluence rate.

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∕3)\pi a^3$ . The radius $a$ is usually chosen to approximately equal ${\mathit{MFP}}^{\prime }∕2$ , 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∕3)\pi a^3$ , so $a={[3V∕\left(4\pi N_a\right)]}^{1∕3}$ . The elements of the $\mathbf{M}$ matrix for $i=1:N_a$ are specified for each $j$ ’th source by:

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:

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:

Figure 9 shows an example of a $0.4\text{-}\mathrm{cm}$ -diam. spherical object with adding absorber ( $\Delta {\mu }_a=1{\mathrm{cm}}^{-1}$ , $\Delta {\mu }_s^{\prime }=0$ , ${\mu }_{a0}=0.1{\mathrm{cm}}^{-1}$ , ${\mu }_s^{\prime }=10{\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{\mathrm{cm}}^{-2}$ effect ( $\mathrm{W}∕{\mathrm{cm}}^2$ per W of the source) within the object and $-5{\mathrm{cm}}^{-2}$ effect in the region surrounding the object. Figure 9c shows the fractional perturbation, ${\varphi }_p∕{\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.

Fig. 9

An object with increased absorption ( ${\mu }_{a0}=0.1{\mathrm{cm}}^{-1}$ , ${\mu }_{s{0}^{\prime }}=10{\mathrm{cm}}^{-1}$ , $\Delta {\mu }_a=1{\mathrm{cm}}^{-1}$ , $\Delta {\mu }_s^{\prime }=0$ ). (a) The background fluence distribution ${\varphi }_0(x,z)$ at $y=0[\mathrm{W}∕{\mathrm{cm}}^2]$ , diffusing from the source. (b) The perturbing fluence ${\varphi }_p$ due to the object. (c) The fraction perturbation ${\varphi }_p∕{\varphi }_0$ . The object asserts a spherical zone of negative perturbation, which casts a perturbing shadow downstream but has little effect upstream toward the source of light.

Figure 10 shows the same spherical object, this time with added scattering ( $\Delta {\mu }_s^{\prime }=10{\mathrm{cm}}^{-1}$ , $\Delta {\mu }_a=0$ , ${\mu }_{a0}=0.1{\mathrm{cm}}^{-1}$ , ${\mu }_s^{\prime }=10{\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∕{\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.

Fig. 10

An object with increase scattering ( ${\mu }_{a0}=0.1{\mathrm{cm}}^{-1}$ , ${\mu }_{s{0}^{\prime }}=10{\mathrm{cm}}^{-1}$ , $\Delta {\mu }_a=0$ , $\Delta {\mu }_s^{\prime }=10{\mathrm{cm}}^{-1}$ ). (a), (b), and (c) are the same as in Fig. 9. The object asserts a dipole pattern of perturbation, in which light is scattered upstream toward the source, and a shadow is cast downstream. However, the perturbation has little effect upstream toward the source of light.

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,⁷ 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.⁸

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\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∕(1-r_i)$ . Therefore, the fluence rate at the surface is

Fig. 11

Schematic explanation of mismatched boundary condition. (a) The photons diffuse in a J flux toward the surface and an $I$ flux away from the surface, and on average travel a distance $L=2\Delta z$ per $\Delta z$ movement along the $z$ axis. The internal reflectance $r_i$ at the surface yields $I=r_{i}J$ . The escaping reflectance $R=(1-r_i)J$ . Hence, the fluence rate at the surface is $\varphi \left(0\right)=2R(1+r_i)∕(1-r_i)$ , Eq. 29. (b) The 1-D fluence rate distribution due to a planar source at $z=z_s$ and an air/tissue boundary at $z_0$ . A surface of symmetry is placed as $z_b$ . If a true positive source of light is present at $z_s$ within the tissue, a negative image source is placed at $z_i=2z_b-z_s$ outside the tissue. The gradient of light between $z_0$ and $z_b$ drives a flux out of the surface, $R=D\varphi \left(z_0\right)∕\mid z_b\mid$ , Eq. 30. Combining Eqs. 29, 30 yields the position $z_b=-2(1+r_i)∕(1-r_i)D$ , which is Eq. 31.

The example situation illustrated in Fig. 11b uses 1-D planar diffusion theory $T_{\text{planar}}\left(z\right)=\exp (-z∕\delta )∕\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{\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 ∕\partial z=\varphi \left(0\right)∕z_b$ . Therefore, the escaping flux is equal to

Therefore, for a generic air/tissue interface, the value of $r_i$ is $\sim 0.5$ , and the value of $(1+r_i)∕(1-r_i)$ is $\sim 3$ . Therefore, $z_b\approx -\left(2\right)\left(3\right)(z_s∕3)\approx -2z_s$ . Hence, the position of the negative image source is placed at $z_i\approx -5z_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:

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{\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{\mathrm{MPF}}^{\prime }$ distant from the origin at $(r=0)$ and greater than $1{\mathrm{MPF}}^{\prime }$ from the source at ( $r=0$ , $z=z_s$ ). Monte Carlo simulations are recommended for regions close to the source.

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:

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 .

Fig. 12

Example meshes are shown for (a) finite difference solver in Cartesian coordiates with adaptive mesh refinement in each of the $x,y,z$ directions; (b) finite element solver using triangular elements in two dimensions or tetrahedral elements in three dimensions; and (c) a boundary element solver where the mesh is just composed of nodes at the boundaries outlining homogeneous regions.

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.

3.2.

Numerical Methods to Solve the Diffusion Equation