Abstract
Background
The bioluminescent enzyme firefly luciferase (Luc) or variants of green fluorescent protein (GFP) in transformed cells can be effectively used to reveal molecular and cellular features of neoplasia in vivo. Tumor cell growth and regression in response to various therapies can be evaluated by using bioluminescent imaging. In bioluminescent imaging, light propagates in highly scattering tissue, and the diffusion approximation is sufficiently accurate to predict the imaging signal around the biological tissue. The numerical solutions to the diffusion equation take large amounts of computational time, and the studies for its analytic solutions have attracted more attention in biomedical engineering applications.
Methods
Biological tissue is a turbid medium that both scatters and absorbs photons. An accurate model for the propagation of photons through tissue can be adopted from transport theory, and its diffusion approximation is applied to predict the imaging signal around the biological tissue. The solution to the diffusion equation is formulated by the convolution between its Green's function and source term. The formulation of photon diffusion from spherical bioluminescent sources in an infinite homogeneous medium can be obtained to accelerate the forward simulation of bioluminescent phenomena.
Results
The closed form solutions have been derived for the timedependent diffusion equation and the steadystate diffusion equation with solid and hollow spherical sources in a homogeneous medium, respectively. Meanwhile, the relationship between solutions with a solid sphere source and ones with a surface sphere source is obtained.
Conclusion
We have formulated solutions for the diffusion equation with solid and hollow spherical sources in an infinite homogeneous medium. These solutions have been verified by Monte Carlo simulation for use in biomedical optical imaging studies. The closed form solution is highly accurate and more computationally efficient in biomedical engineering applications. By using our analytic solutions for spherical sources, we can better predict bioluminescent signals and better understand both the potential for, and the limitations of, bioluminescent tomography in an idealized case. The formulas are particularly valuable for furthering the development of bioluminescent tomography.
Keywords:
Diffusion equation; Green's function; analytical solution; Monte Carlo simulation; bioluminescent imagingBackground
Light propagation in scattering media has attracted attention in many areas of physics, biology, medicine, and engineering. This process can be described by the radiative transfer equation [1], which is a complex integrodifferential equation. Under strong scattering conditions, the radiative transfer equation can be reduced to the diffusion equation, but this still does not permit exact solutions in most cases. Recently, various approximate methods for the radiative transfer equation have been presented [25]. Arridge et al. [2] presented a method for handling nonscattering regions within diffusing domains. The method develops from an iterative radiositydiffusion approach that uses a finite element method. The fundamental idea is to introduce extra equations into the standard diffusion finite element method to represent nondiffusive light propagation across a nonscattering region. Aydin et al. [3] introduced a twodimensional, finite element spherical harmonics radiation transport method for the simulation of light propagation in tissue. They derived a photontransport forward model in turbid media that treats weak inhomogeneities through a Born approximation of the Boltzmann radiative transfer equation. This model more conveniently replaces the commonly used diffusion approximation in optical tomography [4]. Analytic solutions to the diffusion equation in some simplified cases have been explored previously [57]. Boas et al. [5] presented a spherical harmonic solution in the case of spherical inhomogeneity in an infinite medium. In this instance, a closedform infinite series solution was obtained. Additionally, they compared their calculations to experimental results involving a perfect absorber in a large diffuse region. They observed good agreement between the experimental and the analytical data. Walker et al. [6] extended the above spherical results to cylindrical inhomogeneity. Fishkin et al. [7] derived analytical expressions based on the diffusion approximation that describe the photon density in a uniform, infinite, strongly scattering medium that contains a sinusoidally intensitymodulated point source of light. They also reported good agreement between their analytical and experimental results. However, these solutions to the diffusion equation were derived for a point source only.
Recently, light propagation in scattering media has attracted increased attention in biomedical applications [8], especially in bioluminescent imaging as applied in bioluminescent tomography. The bioluminescent enzyme firefly luciferase (Luc) or variants of green fluorescent protein (GFP) in transformed cells can be effectively used to reveal molecular and cellular features of neoplasia in vivo. Tumor cell growth and regression in response to various therapies can be evaluated by using bioluminescent imaging. With a number of key collaborators, we are developing a bioluminescent tomography system for the 3D reconstruction of a bioluminescent source distribution in a living animal marked by bioluminescent reporter luciferases from diffuse signals detected on a spherical surface containing the animal. In bioluminescent experiments, the density of the bioluminescent signals decay exponentially with time. Actually, bioluminescent imaging is a timedependent problem. In this paper, we formulate solutions to the timedependent diffusion equation and the steadystate diffusion equation with solid and hollow spherical sources, respectively, in an infinite homogeneous medium. These formulas are then verified in a Monte Carlo simulation performed using a commercial software package, TracePro. The solution to the diffusion model in a homogeneous medium can be used to achieve reasonably good qualitative estimates of the signal level and can be applied to estimate rapidly the positions of light sources for bioluminescent imaging. The closed form solution is highly accurate and more computationally efficient in biomedical engineering applications.
Methods
The radiative transport equation in a media that scatters and absorbs photons is
L(r, s, t) is the radiance (the power per unit area and unit solid angle) at position r, traveling in direction s, at time t, in the unit of W/(m^{2}•sr). C is the speed of light in the medium. μ_{t }= μ_{s }+ μ_{a }is the optical transport coefficient, where μ_{s }is the scattering coefficient, and μ_{a }is the absorption coefficient. S(r, s, t) is the spatial and angular distribution of the source in the unit of W/(m^{3}•sr). The normalized phase function p(s', s) represents the probability of scattering in a direction s' from direction s. The photon fluence is given by
The photon flux or current density is given by
Both the fluence and the flux are in units of W m^{2}. The transport equation can be regarded as a conservation equation for the radiance; it neglects coherence and polarization effects.
In vivo experimental light propagates in a highly scattering tissue in which the scattering coefficient is much greater than the absorption coefficient, and the diffusion approximation is sufficiently accurate. In the P_{1 }approximation, the radiance is expanded in spherical harmonics truncated to the first order, and the radiance equation can be reduced to the diffusion equation [5,9]
where Φ(r, t) is the photon fluence rate; D = {3[μ_{a }+ (1  g)μ_{s}]}^{1 }is the diffusion coefficient; g is the anisotropy coefficient. Analytic solutions to the diffusion equation are difficult to obtain in most cases and its numerical solutions take large amounts of computational time. By setting S(r, t) = δ(r, t), the timedomain Green's function for an infinite homogeneous medium is given by
The Green's function (5) helps to further simplify the computation of the diffusion equation (4). The solution to the diffusion equation for a general source term S(r, t) is given by the convolution between S(r, t) and the Green's function Γ(r, t):
In the following section, we analytically solve for the photon diffusion with spherical sources embedded in an infinite homogeneous medium.
Solution for a spherical surface source
A spherical source with photon power uniformly distributed on the surface is defined by
where p is the source power; R is the radius of the spherical source; and ω is the decay coefficient of the source power related to time.
In the spherical coordinates, the convolution integral (6) becomes
Integrating with respect to the angular variables θ and φ, Eq. (8) becomes
Letting r ≥ R, , we have
where , and .
Solution for a spherical solid source
A spherical solid source with a radially and exponentially decreasing power density is defined by
where p is the source power; R is the radius of the spherical source; ω is the decay coefficient of the source power related to time; and β is the exponential decay coefficient related to the radius of the spherical source.
The convolution integral (6) becomes
Integrating with respect to the angular variables θ and φ, and changing the order of integration, Eq. (12) becomes
Setting r ≥ R, , Φ(r, t) can be expressed as
Performing an integration around variable t', we obtain
Integrating with respect to variables r', we obtain
where
and E(R, k, a, b) is defined as
Setting β = 0, we obtain the solution for a spherical solid source with a constant power density as follows:
where
Solutions to the steadystate diffusion equation
Similarly, expressing the solution of the steadystate diffusion equation as the convolution between S(r, t) and its Green's function, and simplifying the convolution, the solutions to the steadystate diffusion equation can be easily derived for spherical surface and solid sources, respectively. In this case, the solution for a spherical surface source can be shown to be
and the solution for a spherical solid source can be expressed as
where p is the source power; R is the radius of the spherical source; and μ_{eff }denotes the effective attenuation coefficient defined as μ_{eff }= (μ_{a }/D)^{1/2 }= [3 μ_{a }(μ_{a }+ μ_{s }(1  g))]^{1/2}.
Results
To verify the correctness of the above analytic results, two sets of computer simulation experiments were performed using the wellknown MatLab and a commercial software package, TracePro (Lambda Research Corp., Littleton, MA, USA; http://www.lambdares.com/products/tracepro/index.phtml webcite). The first set of tests examined the relationship between the solutions for spherical surface and solid sources. The other set of tests demonstrated good agreement between the analytically computed and the Monte Carlo simulated data.
Comparison between solutions for spherical solid and surface sources
The numerical tests were conducted with various parameters using solutions for spherical surface and solid sources. In the first experimental setting, a spherical solid source of radius 1 mm with power 0.0002513 (W) was embedded in an infinite homogeneous medium. The medium was defined by μ_{s }= 40.3 mm^{1}, μ_{a }= 0.15 mm^{1}, n = 1, and g = 0.994. In the second setting, a spherical solid source of radius 0.3 mm with power 0.0002513 (W) was assumed, with a medium of μ_{s }= 10.27 mm^{1}, μ_{a }= 0.082 mm^{1}, n = 1.37, and g = 0.96. Figs. 1 and 2 present some typical analytically calculated outcomes. These results indicate that the directly computed data using the solution formula for a spherical solid source are in an excellent agreement with their counterparts indirectly computed using the solution for a spherical surface source. In other words, the two solutions are consistent.
Figure 1. Comparison of flux computed on a spherical surface in the first experimental setting according to the solutions for spherical surface and solid sources, respectively.
Figure 2. Comparison of flux computed on a spherical surface computed in the second experimental setting according to the solutions for spherical surface and solid sources, respectively.
Comparison between analytic solution and Monte Carlo simulation
In the following two numerical experimental settings, the optical parameters of the medium were set at μ_{s }= 23 mm^{1}, μ_{a }= 0.35 mm^{1}, n = 1, and g = 0.94, the same as the optical parameters of a mouse lung. The source, a sphere with a radius of 1.3 mm and power 0.00025 (W), was embedded in an infinite homogeneous medium. The power incidence on the spherical surfaces with different radii was computed by using the solution formulas and the Monte Carlo method as implemented in TracePro, respectively. Comparison of the two is shown in Fig. 3. In the other numerical experiment, the optical parameters of the medium were μ_{s }= 16 mm^{1}, μ_{a }= 0.20 mm^{1}, n = 1.37, and g = 0.85, the same as the optical parameters of a mouse heart. The source is a sphere of radius 1 mm with power 0.00025 (W). Fig. 4 presents a comparison of the analytic results and the simulated data from TracePro. It can be observed that there is good agreement between the computed and the simulated data.
Figure 3. Comparison of flux computed on a spherical surface in the first experimental setting based on the analytic solution and the Monte Carlo simulation, respectively.
Figure 4. Comparison of flux computed on a spherical surface in the second experimental setting based on the analytic solution and the Monte Carlo simulation, respectively.
Discussions and conclusion
The primary utility of our analytic solutions is to accelerate the forward simulation of bioluminescent phenomena which will facilitate studies on bioluminescent tomography. The spherical sources we assumed are the most fundamental building blocks of the physical phantoms used in the evaluation of image quality. By using our analytic solutions for spherical sources, we can better predict bioluminescent signals and better understand both the potential for, and the limitations of, bioluminescent tomography in an idealized case.
In conclusion, we have formulated solutions for the diffusion equation with solid and hollow spherical sources in an infinite homogeneous medium. These solutions have been verified by Monte Carlo simulation for use in biomedical optical imaging studies. The formulas are particularly valuable for furthering the development of bioluminescent tomography.
APPENDIX: Integral formulas
Acknowledgment
This work was supported in part by the NIH grant R21/R33 EB001685.
References

Case KM, Zweifel PF: Linear Transport Theory. Mass: AddisonWesley; 1967.

Arridge SR, Deghani H, Schweiger M, Okada E: The finite element model for the propagation of light in scattering media: a direct method for domains with nonscattering regions.
Med Phys 2000, 27:252264. PubMed Abstract  Publisher Full Text

Aydin ED, Oliveira CR, Goddard AJ: A comparison between transport and diffusion calculations using a finite elementspherical harmonicsradiation transport method.
Med Phys 2002, 29:20132023. PubMed Abstract  Publisher Full Text

Xu M, Cai W, Lax M, Alfano RR: Photon transport forward model for imaging in turbid media.

Boas DA: Diffuse Photon Probes of Structural and Dynamical Properties of Turbid Media, PHD. Thesis, University of Pennsylvania, Philadelphia; 1996.

Walker SA, Boas DA, Gratton E: Photon density waves scattered from cylindrical inhomogeneities: theory and experiments.

Fishkin JB, Gratton E: Propagation of photondensity waves in strongly scattering media containing an absorbing semiinfinite plane bounded by a straight edge.
J Opt Soc Am A 1993, 10:127140. PubMed Abstract

Rice BW, Cable MD, Nelson MB: In vivo imaging of lightemitting probes.
J of Bio Opt 2001, 6:432440. Publisher Full Text

Arridge SR, Schweiger M, Hiraoka M, Delpy DT: A finite element approach for modeling photon transport in tissue.
Med Phys 1993, 20:299309. PubMed Abstract  Publisher Full Text

Wang LH, Jacques SL, Zheng LQ: MCMLMonte Carlo modeling of photon transport in multilayered tissues.
Comput Methods Programs Biomed 1995, 47:131146. PubMed Abstract  Publisher Full Text

Patterson MS, Chance B, Wilson BC: Time resolved reflectance and transmittance for the noninvasive measurement of optical properties.

Markel VA, Schotland JC: Inverse problem in optical diffusion tomography. 1 FourierLaplace inversion formulas.
J Opt Soc Am A Opt Image Sci Vis 2001, 18:13361347. PubMed Abstract

Haskell RC, Svaasand LO, Tsay T, Feng T, McAdams MS, Tromberg BJ: Boundary conditions for the diffusion equation in radiative transfer.
J Opt Soc Am A 1994, 11:27272741. PubMed Abstract

Flock ST, Patterson MS, Wilson BC, Wyman DR: Monte Carlo modeling of light propagation in highly scattering tissues 1: Model predictions and comparison with diffusion theory.
IEEE Trans Biomed Eng 1989, 36:11621168. PubMed Abstract  Publisher Full Text

Sassaroli A, Blumetti C, Martelli F, Alianelli L, Contini D, Ismaelli A, Zaccanti G: Monte Carlo procedure for investigating light propagation and imaging of highly scattering media.