Open Access

Formulation of photon diffusion from spherical bioluminescent sources in an infinite homogeneous medium

BioMedical Engineering OnLine20043:12

DOI: 10.1186/1475-925X-3-12

Received: 15 October 2003

Accepted: 04 May 2004

Published: 04 May 2004

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 time-dependent diffusion equation and the steady-state 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 imaging

Background

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 radiosity-diffusion 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 two-dimensional, finite element spherical harmonics radiation transport method for the simulation of light propagation in tissue. They derived a photon-transport 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 closed-form 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 intensity-modulated 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 time-dependent problem. In this paper, we formulate solutions to the time-dependent diffusion equation and the steady-state 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

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equa_HTML.gif

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/(m2•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/(m3•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

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equb_HTML.gif

The photon flux or current density is given by

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equc_HTML.gif

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 P1 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]

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equd_HTML.gif

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 time-domain Green's function for an infinite homogeneous medium is given by

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Eque_HTML.gif

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

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equf_HTML.gif

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

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equg_HTML.gif

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

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equh_HTML.gif

Integrating with respect to the angular variables θ and φ, Eq. (8) becomes

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equi_HTML.gif
Letting r ≥ R,
https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equj_HTML.gif
, we have
https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equk_HTML.gif
where
https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equl_HTML.gif
, and https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_IEq1_HTML.gif .

Solution for a spherical solid source

A spherical solid source with a radially and exponentially decreasing power density is defined by

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equm_HTML.gif

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

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equn_HTML.gif

Integrating with respect to the angular variables θ and φ, and changing the order of integration, Eq. (12) becomes

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equo_HTML.gif
Setting r ≥ R,
https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equp_HTML.gif
, Φ(r, t) can be expressed as
https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equq_HTML.gif

Performing an integration around variable t', we obtain

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equr_HTML.gif

Integrating with respect to variables r', we obtain

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equs_HTML.gif

where

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equt_HTML.gif

and E(R, k, a, b) is defined as

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equu_HTML.gif

Setting β = 0, we obtain the solution for a spherical solid source with a constant power density as follows:

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equv_HTML.gif

where

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equw_HTML.gif

Solutions to the steady-state diffusion equation

Similarly, expressing the solution of the steady-state diffusion equation as the convolution between S(r, t) and its Green's function, and simplifying the convolution, the solutions to the steady-state 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

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equx_HTML.gif

and the solution for a spherical solid source can be expressed as

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equy_HTML.gif

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 μaa + μ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 well-known MatLab and a commercial software package, TracePro (Lambda Research Corp., Littleton, MA, USA; http://www.lambdares.com/products/tracepro/index.phtml). 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.
https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Fig1_HTML.jpg
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.

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Fig2_HTML.jpg
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.
https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Fig3_HTML.jpg
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.

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Fig4_HTML.jpg
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

https://static-content.springer.com/image/art%3A10.1186%2F1475-925X-3-12/MediaObjects/12938_2003_Article_39_Equz_HTML.gif

Declarations

Acknowledgment

This work was supported in part by the NIH grant R21/R33 EB001685.

Authors’ Affiliations

(1)
CT/Micro-CT Laboratory University of Iowa Iowa City
(2)
Optical Imaging Laboratory Department of Biomedical Engineering Texas A&M University College Station

References

  1. Case KM, Zweifel PF: Linear Transport Theory. Mass: Addison-Wesley 1967.Google Scholar
  2. 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: 252–264. 10.1118/1.598868View ArticleGoogle Scholar
  3. Aydin ED, Oliveira CR, Goddard AJ: A comparison between transport and diffusion calculations using a finite element-spherical harmonicsradiation transport method. Med Phys 2002, 29: 2013–2023. 10.1118/1.1500404View ArticleGoogle Scholar
  4. Xu M, Cai W, Lax M, Alfano RR: Photon transport forward model for imaging in turbid media. Opt Lett 2001, 26: 1066–1068.View ArticleGoogle Scholar
  5. Boas DA: Diffuse Photon Probes of Structural and Dynamical Properties of Turbid Media, PHD. Thesis, University of Pennsylvania, Philadelphia 1996.Google Scholar
  6. Walker SA, Boas DA, Gratton E: Photon density waves scattered from cylindrical inhomogeneities: theory and experiments. Appl Opt 1998, 37: 1935–1944.View ArticleGoogle Scholar
  7. Fishkin JB, Gratton E: Propagation of photon-density waves in strongly scattering media containing an absorbing semi-infinite plane bounded by a straight edge. J Opt Soc Am A 1993, 10: 127–140.View ArticleGoogle Scholar
  8. Rice BW, Cable MD, Nelson MB: In vivo imaging of light-emitting probes. J of Bio Opt 2001, 6: 432–440. 10.1117/1.1413210View ArticleGoogle Scholar
  9. Arridge SR, Schweiger M, Hiraoka M, Delpy DT: A finite element approach for modeling photon transport in tissue. Med Phys 1993, 20: 299–309. 10.1118/1.597069View ArticleGoogle Scholar
  10. Wang LH, Jacques SL, Zheng LQ: MCML-Monte Carlo modeling of photon transport in multilayered tissues. Comput Methods Programs Biomed 1995, 47: 131–146. 10.1016/0169-2607(95)01640-FView ArticleGoogle Scholar
  11. Patterson MS, Chance B, Wilson BC: Time resolved reflectance and transmittance for the non-invasive measurement of optical properties. Appl Opt 1989, 28: 2331–2336.View ArticleGoogle Scholar
  12. Markel VA, Schotland JC: Inverse problem in optical diffusion tomography. 1 Fourier-Laplace inversion formulas. J Opt Soc Am A Opt Image Sci Vis 2001, 18: 1336–1347.MathSciNetView ArticleGoogle Scholar
  13. 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: 2727–2741.View ArticleGoogle Scholar
  14. 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: 1162–1168. 10.1109/TBME.1989.1173624View ArticleGoogle Scholar
  15. 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. Appl Opt 1998, 37: 7392–7400.View ArticleGoogle Scholar

Copyright

© Cong et al; licensee BioMed Central Ltd. 2004

This article is published under license to BioMed Central Ltd. This is an Open Access article: verbatim copying and redistribution of this article are permitted in all media for any purpose, provided this notice is preserved along with the article's original URL.