Open Access
Add to BibSonomyAbstract
In this paper a forward solver software for the time domain and the CW domain based on the Born approximation for simulating the effect of small localized fluorophores embedded in a non-fluorescent biological tissue is proposed. The fluorescence emission is treated with a mathematical model that describes the migration of photons from the source to the fluorophore and of emitted fluorescent photons from the fluorophore to the detector for all those geometries for which Green’s functions are available. Subroutines written in FORTRAN that can be used for calculating the fluorescent signal for the infinite medium and for the slab are provided with a linked file. With these subroutines, quantities such as reflectance, transmittance, and fluence rate can be calculated.
© 2011 Optical Society of America
1. Introduction
Fluorescence spectroscopy and imaging using near infrared light are methods of analysis that have become emerging and promising techniques when applied in medicine [1–7]. Applications rely on measurements of the steady-state emission spectrum when continuous wave (CW) sources are used, on measurements of the time-resolved (TR) fluorescence signal when pulsed laser sources are employed (time domain), or on measurements of frequency domain signals (amplitudes and phase-shift measurements) when modulated laser sources are employed [8]. The development of fluorescence applications has been accompanied by the elaboration of several forward solvers in order to obtain the fluorescence signal. In the literature we have examples of forward solvers based on solutions of the diffusion equation (DE) [9–12], on the random walk theory [13], on the spherical harmonics approximation [14], and on numerical solutions of the radiative transfer equation (RTE) obtained using the Monte Carlo method [15,16] or with the finite-difference discrete-ordinates method [17]. The Born approximation applied to the DE has also been used for measurements in the frequency domain [18, 19]. For the inverse problem O’Leary et al. [18] used a first-order perturbation model for the diffuse photon density waves emitted from a localized fluorophore inside an infinite medium. Ntziachristos and Weissleder [19] used a normalized Born approximation of the DE to develop an algorithm for the tomographic reconstruction of small fluorescent targets inside an infinite medium. Li et al. [20] have presented solutions for the fluorescent diffuse photon density waves of a fluorescent target embedded inside an infinite or a semi-infinite fluorescent diffusive medium. These solutions for the frequency domain are not analytical in a proper sense since numerical evaluations of integrals are required.
We propose a perturbative time domain and the CW domain forward solver software based on the Born approximation in order to simulate the effect of localized fluorophores embedded in a non fluorescent biological tissue. The fluorescence emission is treated with a mathematical model that is capable of describing the migration of photons from the source to the fluorophore and of emitted fluorescent photons from the fluorophore to the detector, for all those geometries for which Green’s functions can be calculated. The method is implemented here for the infinite medium and for the slab for which some subroutines written in FORTRAN are provided with the linked file (Media 1) . Quantities such as reflectance, transmittance, and fluence rate can be calculated by using these subroutines and the software package distributed with a recently published book [21]. The forward solver software has been developed for the time domain and the CW domain, however guidelines that can extend its use also to the frequency domain have been provided as well.
The software presented in this paper is a tool for photon migration studies that can be utilised when a fluorescent inclusion is present in tissue. For instance, this is the case of optical molecular imaging that uses near-infrared fluorescent probes [14, 17]. In optical molecular imaging, a fluorescent biochemical marker is injected into tissue and will emit near-infrared light upon excitation from an external light source. By using measurements of the light intensity on the tissue surface, it is possible to determine the spatial concentration distribution of the marker inside tissue by recording tomographical data sets and employing appropriate tomographic image reconstruction schemes [22]. The inverse fluorescent source problem in optical molecular imaging has been solved by using the DE [19, 23, 24] or the RTE [17, 25] as a forward model. Thus, the availability of fast and precise forward solvers is fundamental for the inverse problem in optical molecular imaging.
The theory is presented in section 2; the software is described in section 3 and in the body text of the subroutines provided with the linked file (Media 1). In section 3, some examples of use of the software are also provided together with a comparison with the results of Monte Carlo simulations. Lastly, in section 4 a summary of the paper and a discussion of possible further developments of the software are provided.
2. Theory
The forward solver which we propose considers the effects of fluorescence on photon migration by the probability of having fluorescent photons emitted inside the medium. The fluorescent signal is affected by multiple scattering at both the excitation and the emission wavelengths. In Fig. 1, a diagram shows that the excitation light at wavelength λ is an isotropic point source at rs, and that the fluorescent light at wavelength λe, which is emitted from a fluorophore at r′, is collected at point r. The absorption coefficient, the reduced scattering coefficient, and the refractive index of the background medium at λ are denoted, respectively, as μa, μ′s, and n, and at λe as μae, μ′se, and ne. We assume that a fluorophore is distributed inside a volume V′ with absorption coefficient and reduced scattering coefficient μaf and μ′sf at λ, and μafe and μ′sfe at λe.
Fig. 1 Diagram of a photon path from the source to the fluorophore and from the fluorophore to the detector.
A similarity can be observed between the behavior of an absorbing inclusion and that of a fluorescent inclusion when these are inside a scattering medium. Both of them subtract energy from propagation in accordance to the phenomena of photons absorption. Therefore, the expression of the perturbed fluence due to an absorbing defect obtained with the Born approximation (see, for instance, chapter 7 in [21]) can be used to describe the photons absorbed by a fluorescent inclusion. The signal at emission wavelength can then be calculated by using the probability function that the fluorophore re-emits the absorbed photons.
2.1. Perturbative model based on the DE
Light propagation through biological tissue can be described with the time-dependent DE [21]. Therefore, the phenomena of fluorescence can be described by coupled diffusion equations at excitation and emission wavelengths that are respectively, [10]:
where D = 1/[3(μ′s + μ′sf)] and De = 1/[3(μ′se + μ′sef)] are the diffusion coefficients at excitation and emission wavelengths. Equations (1) and (2) are coupled by the excitation fluence, Φ, which determines the source term qe(r,t) of Eq. (2). In particular, the probability dPf that fluorescent photons are generated in the time interval dt′ at time t′, in the wavelength interval dλe around λe, and in the volume element d3r′ at r′, is where η(λ,λe) is the quantum yield or quantum efficiency per unit wavelength for emission at λe given excitation at λ, and G is the Green’s function of the medium for the fluence at the excitation wavelength. Therefore, the coupling of Eqs. (1) and (2) is given by the term where τ is the fluorescent lifetime of the fluorophore (i.e. the average time that the fluorophore spends in the excited state prior to returning to the ground state). By substituting Eq. (4) in Eq. (2), the following equation at emission wavelength is obtained: with η→e = η(λ,λe)dλe quantum efficiency of the fluorophore in the wavelength range dλe around λe, i.e. the ratio of the number of photons emitted in the wavelength range dλe to the number of absorbed photons. According to Green’s function theorem, the fluence rate at emission wavelength, Φe (r,t), due to the fluorescent photons generated in the interval dλe, is with Ge Green’s function of the medium for the fluence at the emission wavelength. We stress that the Green’s functions G and Ge are the Green’s function of the same medium, so that the different notation is maintained only to distinguish between the excitation and the emission wavelengths. If all the fluorescence photons were promptly emitted, the fluence at emission wavelength would simply be given by Equation (6) is thus the convolution of Φe0 with the exponential decay of the fluorophoreIt can be noted that the integral of Eq. (7) is similar to the expression of the perturbed fluence due to an absorbing defect placed inside a turbid medium obtained using the Born approximation (see, for instance, chapter 7 in [21], Eq. (7.11) for the DE or Eq. (7.63) for the RTE). Therefore, the Born approximation can also be used to calculate the fluorescence signal. That is, if the fluorophore occupies a small volume, and if at excitation wavelength it perturbs the fluence rate of the background medium of a small quantity then the following approximations can be made:
and the calculation of the TR expressions for the fluence Φe0 (r,t) and Φe (r,t) can simply be obtained as These equations are perturbative solutions for a fluorescence emission in which the integral on V′ is omitted, and are valid only for small fluorescent inclusions.Equations (10) and (11) were obtained using the DE, but can also be obtained by using the RTE [21]. We stress that they can be calculated for any geometry for which analytical or numerical Green’s functions for the DE or the RTE are available. For G and Ge the solutions and the software presented in [21] can be used and Eqs. (10) and (11) can be therefore calculated for all the geometries analyzed in [21].
For the sake of clarity we stress that, for the theory presented in this section as well as the models proposed in [18,19], the Born approximation is applied twice: i.e., a first time at excitation wavelength, and a second time at emission wavelength. In fact, both the Green’s functions G and Ge are approximated to that of the background medium. Therefore, these models are similar but not exactly equal to the Born approximation.
The following model is expected to fail when the fluorophore significantly perturbs the fluence at excitation wavelength and thus when the approximations of Eqs. (9) are no longer valid and so the actual expressions of G and Ge should be used. In this case, a Green’s function of the medium obtained with the high order perturbation theory could be used to improve the forward solver [26–28].
The theory has been developed for time-dependent sources. It is easy to extend this approach to CW sources. By eliminating the time dependence in Eq. (1), a CW perturbative solution similar to Eq. (10) can be written as follows:
Last, we considered the frequency-domain: When modulated sources with a DC and a AC component were considered [29], the source term of Eq. (1) could be written as
where ω = 2πf is the pulsation and f is the frequency of modulation. q0DC and q0AC are the DC and AC parts of the source strength. The solution for the fluence can be written as the sum of the DC and AC solutions, i.e., Φ(r,t) = ΦDC(r)+ΦAC(r,t). The AC solution of the DE with the source term given by Eq. (13) must be oscillating at the same pulsation ω as the source. Thus, the AC solution of Eq. (1) will have the form: If we substitute the fluence in Eq. (1) with Eq. (14), we obtain that ΦAC0 can be described by the equation Equation (15) does not have a time-dependent term, and it is a CW diffusion equation. Thus, the general solution for ΦAC0 can be obtained by using the CW solution, provided that the coefficient μa is replaced with the term [μa – (iω)/v].2.2. Hybrid approach for the slab
The theory described above is obtained by using the DE and is not restricted to any particular geometry. We should mention that Eqs. (11) and (12) can also be obtained by using the RTE [21]. In general, solutions based on the DE are affected by the intrinsic approximations of this theory, which show different consequences for steady-state and time-dependent sources. For steady-state sources, solutions based on the DE provide an insufficient description of photons detected at short distances [21]. For time-dependent sources, DE solutions provide an insufficient description of early received photons [21]. The limitations of the DE are a serious problem in modeling photon migration, especially when small volumes of diffusive media are considered. To overcome these limitations, solutions that are an improvement of the DE solutions [30–34], hybrid solutions [35–37], or RTE solutions [38–40] have been provided.
To evaluate the fluorescence signal with Eqs. (11) and (12), analytical expressions for the TR and CW Green’s functions at excitation and emission wavelengths are required. The Green’s function for the slab is calculated by using the hybrid approach described in [21], which makes use of RTE solutions for the infinite medium, together with the DE framework approach, which is based on the method of images. For the time domain (Eq. (11)), the solutions for the infinite medium that can be used are the solution of the DE, the solution of the telegrapher equation, and the solution of the RTE for isotropic scattering proposed by Paasschens [38]. For the CW domain (Eq. (12)), these solutions for the infinite medium are considered: the solution of the DE and the improved solution of the DE (IDE) for isotropic scattering proposed by Graff and Rinzema [32]. The solution by Graff and Rinzema will be denoted inside the software section as IDE solution in the CW domain.
The limitations of the solutions implemented for the slab, i.e. of Eqs. (11) and (12) calculated using the hybrid approach adopted in [21], can be ascribed to: I) approximations due to the assumptions contained in Eqs. (9) that are similar to the Born approximation; II) approximations of the hybrid Green’s functions for the slab described in [21]. Improvements on this theory could be obtained by implementing: A) a better approximation for the Green’s function describing the fluence at the fluorophore site, which could be obtained by using high-order solutions of the perturbation theory [26,27,42] or exact RTE solutions [39]; B) the whole integral on the volume V′.
3. Software and examples of use
With this paper (see Media 1), we provide some subroutines written in FOR-TRAN that together with the subroutines described in [21] can be used for calculating Eqs. (11) and (12) for the infinite medium and for the homogeneous slab. The Homogeneous_Fluence_Fluorescence_Infinite subroutine is the basis for calculating the fluorescence in the infinite medium. The calculation can be implemented in the time-domain by using Eq. (11) and in the CW domain by using Eq. (12). The subroutines necessary for the calculation are RTE_Fluence_Infinite_TR (solution by Paasschens [38]), DE_Fluence_Infinite_TR, DE_Fluence_Infinite_CW, and IDE_Fluence_Infinite_CW (solution by Graaff and Rinzema [32]). The code of these subroutines can be found in [21] or in the linked file (Media 1). The Homogeneous_Fluence_Fluorescence_Slab subroutine is the basis for calculating the fluorescence in the slab. The subroutines necessary for the calculation are the same mentioned above.
Some instructions for the use of the Homogeneous_Fluence_Fluorescence_Infinite and Homogeneous_Fluence_Fluorescence_Slab subroutines are provided within the body of the text on the subroutines. The reflectance, R(ρ,t), that is due to the fluorophore can be calculated by using Fick’s law or the partial current boundary condition [21]:
with D being the diffusion coefficient, A the coefficient for the Fresnel Reflections [21], and Φ the fluorescent fluence calculated with the proposed subroutines.We stress that when Eqs. (11) and (12) are calculated using RTE Green’s functions for the infinite medium, the limit of the calculated solution tends to become exact as the volume of the fluorescent inclusion approaches to zero. Therefore, the proposed forward solver can be exact for small fluorescent inclusions. By making use of the recent RTE Green’s function obtained by Liemert and Kienle [40] or by Machida et al. [39], it is possible to build an exact forward solver for small localized fluorescent inclusions even for anisotropic scattering. Last, we stress that subroutines for calculating the fluorescence signal emitted by a localized fluorescent inclusion placed inside a two-layer medium, a sphere, a parallelepiped, and a cylinder, can be obtained by using the software package provided in [21] and by implementing for these geometries a similar approach to that used for the infinite medium and the slab.
In order to cite an example of use of the proposed software, we consider a fluorophore embedded in a slab 2000 mm thick at ρ =15 mm. The fluorophore is assumed by a sphere centered at (7.5, 0, 7.5) mm with a volume V′ = 1.25 mm3. The absorption and the reduced scattering coefficient of the background medium are equal to 0.01 and 1. mm−1, respectively, at both the emission and the excitation wavelengths. The refractive index of the medium and of the external is 1.4 at both the emission and the excitation wavelengths. An excitation spatial and temporal isotropic Dirac delta source of unitary strength is at (0, 0, 1/μ′s). The absorption coefficient of the fluorophore at the excitation wavelength is μaf = 0.002 mm−1 and 0 at the emission wavelength and its quantum efficiency is unitary. The same scattering properties as those for the background medium were present in the fluorophore. In Fig. 2 the TR reflectance fluorescence signal is obtained with the Homogeneous_Fluence_Fluorescence_Slab subroutine (that calculate Eq. (11) with the hybrid model for the slab) and with the partial current boundary condition. The hybrid model is calculated by using Paasschens RTE solution for the infinite medium [38]. The TR reflectance shown in Fig. 2 is the probability per unit time and per unit surface area that a photon emitted by the fluorophore exits from the medium at distance ρ and time t. In Fig. 2 we have assumed that fluorescence photons are emitted with a τ =1 ns lifetime. The computation time for calculating Eq. (11) depends on the number of temporal intervals used for calculating the time convolution integral. This parameter can be set to the right value inside the Homogeneous_Fluence_Fluorescence_Slab subroutine (Integer parameter denoted nd). By using 200 time windows (it is sufficient to have a negligible error in the calculated fluence due to the numerical evaluation of the integral): the computation time for calculating Eq. (11) for a single time is about 0.1 s with an Intel Centrino2 processor. In Fig. 3 a map of the CW reflectance fluorescence signal (Eq. (12)) is plotted that was obtained using the hybrid model (Homogeneous_Fluence_Fluorescence_Slab subroutine) based on the IDE solution by Graaff and Rinzema and by using the same case as Fig. 2. The unit of the plotted CW reflectance is mm−2, and represents the probability per unit surface area that a photon emitted by the fluorophore exits the from slab at a certain position on the z = 0 plane. The computation time of the whole map of Fig. 2 (mapping density of 1600 points) is about 0.01 s.
Fig. 2 Fluorescence TR reflectance signal received from a slab 2000 mm thick at ρ =15 mm. The absorption and the reduced scattering coefficient of the background medium are equal to 0.01 and 1. mm−1, respectively, at both the emission and the excitation wavelengths. The fluorophore is a sphere centered at (7.5, 0, 7.5) mm and with a volume V′ = 1.25 mm3. The absorption coefficient of the fluorophore at excitation wavelength is μaf = 0.002 mm−1. The refractive index of the medium and of the external is 1.4.
In Fig. 4 the TR reflectance fluorescence signal obtained with the Homogeneous_Fluence_Fluorescence_Slab subroutine and with the partial current boundary condition is compared with the reflectance simulated with a Monte Carlo code for a slab 2000 mm thick with the receiver at ρ =30 mm. For obtaining the Monte Carlo results we have modified a previously developed Monte Carlo code for describing the effects of absorbing and scattering defects [43]. The reduced scattering coefficient of the background medium is assumed equal to 0.5 mm−1 at both the emission and the excitation wavelengths. The fluorophore is assumed by a sphere centered at (15, 0, 15) mm. The same scattering properties as those for the background medium were present in the fluorophore. In the MC simulations four different couples of values for the volume V′ of the sphere and for the absorption μaf of the fluorophore at the excitation wavelength were considered (for each case we have V′μaf = 0.01 mm2): V′ = 10 mm3, μaf = 0.001 mm−1; V′ = 100 mm3, μaf = 0.0001 mm−1; V′ = 500 mm3, μaf = 0.00002 mm−1, and V′ = 1000 mm3, μaf = 0.00001 mm−1. The absorption coefficient of the fluorophore at emission wavelength is 0. The absorption coefficient of the background medium is 0.01 mm−1 in Fig. 4a, and 0 in Fig. 4b at both the emission and the excitation wavelengths. The refractive index of the medium and of the external is 1.4 at both the emission and the excitation wavelengths. In the comparison of Fig. 4 we have assumed that fluorescence photons are promptly emitted. Thus in Eq. (11), and consequently inside the Homogeneous_Fluence_Fluorescence_Slab subroutine, a τ =10−10 ps is set. The Monte Carlo results are obtained with an isotropic scattering function. The hybrid model is calculated by using in Eq. (11) and in the Homogeneous_Fluence_Fluorescence_Slab subroutine the RTE solution for the infinite medium. For the four different combinations of the volume V′ and of the absorption μaf, the RTE hybrid model provides the same results since Eq. (11) depends only on the product V′μaf that is constant. These results show that for the smallest volume V′ the model is in pretty good agreement with the Monte Carlo results, while some deviations can be appreciated for the higher values of the volume V′. This fact confirms the expectations that the model provide a non sufficient description of photon migration for large volumes of the fluorophore. Similar results were also obtained by calculating the hybrid model for the slab with the DE solution for the infinite medium.
Fig. 4 Comparison between hybrid model and Monte Carlo results for the fluorescence TR reflectance signal received from a slab 2000 mm thick with the receiver at ρ = 30 mm. The reduced scattering coefficient of the background medium is assumed at the emission and the excitation wavelength equal to 0.5 mm−1. The absorption coefficient of the background medium is 0.01 mm−1 in Fig. a), and 0 in Fig. b). In the MC simulations four different values of the volume and of the absorption of the fluorophore have been considered (for each case we have V′μaf = 0.01 mm2): V′ = 10 mm3, μaf = 0.001 mm−1; V′ = 100 mm3, μaf = 0.0001 mm−1; V′ = 500 mm3, μaf = 0.00002 mm−1, and V′ = 1000 mm3, μaf = 0.00001 mm−1. The fluorophore is assumed by a sphere centered at (15, 0, 15) mm. The refractive index of the medium and of the external is 1.4. The hybrid model is calculated with Eq. (11) and with the partial current boundary condition, and it has been implemented using the RTE solution for the infinite medium.
4. Summary and discussion
We proposed a perturbative forward solver software for the time domain and the CW domain based on the Born approximation in order to simulate the effect of small localized fluorophores embedded inside a non-fluorescent biological tissue. The fluorescence emission is treated with a model that describes the migration of photons from the source to the fluorophore and of emitted fluorescent photons from the fluorophore to the detector, for all those geometries for which Green’s functions can be calculated. The method has been implemented for the infinite medium and for the slab for which two subroutines (Homogeneous_Fluence_Fluorescence_Infinite and Homogeneous_Fluence_Fluorescence_Slab) written in FORTRAN language are provided with the linked file (Media 1). Quantities such as reflectance, transmittance, and fluence rate can be calculated using these subroutines and the software package distributed with a recently-published book [21]. Similar subroutines can be also obtained for other geometries such as the two-layer slab, the sphere, the parallelepiped, and the cylinder, by using the software distributed in [21] and by implementing for these geometries similar subroutines to those provided for the infinite medium and the slab. The software has been developed for the time and the CW domains, however it can also be extended to the frequency domain by following the guidelines provided in Sect. 2. The computation time of the said forward solver is significantly less than any elementary Monte Carlo program dedicated to the same kind of calculation. Typically, the computation time for calculating the TR reflectance (at a single value of time) emitted by a fluorescent inclusion is about 0.1 s and the computation time for calculating the CW reflectance is less than 10−5 s. The accuracy of the proposed forward solver is excellent for small localized fluorescent sources, while it decreases when larger fluorescent inclusions are considered as it is shown in Fig. 4. The main use of forward solver software that we have proposed is for photon migration through biological tissues with embedded small localized fluorescent inclusions, such as those encountered in optical molecular imaging.
The proposed forward solver is not meant to replace sophisticated numerical algorithms that are developed for solving partial differential equations in very general conditions. Nonetheless, we can apply this method to a certain number of geometries for which the Green’s function can be expressed as a closed-form solution or it can be calculate numerically, and it has the advantage of being fast and precise for small localized fluorescent perturbations. In the case of larger fluorescent inclusions, it is possible to extend the validity of the theory to a wider range of fluorescent inclusions by implementing the perturbation theory beyond the limit of the Born approximation [26–28,41,42]. In particular, the solutions obtained by Sassaroli et al. [28,41,42] for the higher order perturbation theory for CW, time and frequency domains have shown to be sufficiently accurate to describe correctly absorbing relative perturbations of intensities of up to 50–60%. For this reason, the forward solver that we propose could be further improved by using the perturbed solutions proposed in [28, 41, 42] for the Green’s functions G and Ge of Eqs. (11) and (12).
In recent decades, several analytical and numerical solutions of the RTE have been reported. For instance, Machida et al. numerically obtained the Green’s function for the RTE in the slab geometry for anisotropic scattering [39], and Liemert and Kienle obtained a CW analytical Green’s function of the RTE for the infinite medium which is valid for anisotropic scattering [40]. We wish to emphasize that the availability of an analytical Green’s function valid for anisotropic scattering makes it possible to further extend the perturbative approach presented in this paper.
Acknowledgments
This work has been partially supported with funding from the European Communitys Seventh Framework Program (FP7/2007-2013) under grant agreement No. HEALTH-F5-2008-201076.
References and links
1. S. Andersson-Engels and B. C. Wilson, “In vivo fluorescence in clinical oncology: fundamental and practical issues,” J. Cell. Pharmacol. 3, 66–79 (1992).
2. D. R. Braichotte, J. F. Savary, P. Monnier, and H. E. van den Bergh, “Optimizing light dosimetry in photodynamic theraphy of early stage carcinomas of the esophagus using fluorescence spectroscopy,” Laser Surg. Med. 19, 340–346 (1996). [CrossRef]
3. I. J. Bigio and J. R. Mourant, “Ultraviolet and visible spectroscopies for tissue diagnostics: fluorescence spectroscopy and elastic-scattering spectroscopy,” Phys. Med. Biol. 42, 803–814 (1997). [CrossRef] [PubMed]
4. G. A. Wagnieres, W. M. Star, and B. C. Wilson, “In vivo fluorescence spectroscopy and imaging for oncological applications,” Photochem. Photobiol. 68, 603–632 (1998). [PubMed]
5. A. Liebert, H. Wabnitz, H. Obrig, R. Erdmann, M. Möller, R. Macdonald, H. Rinneberg, A. Villinger, and J. Steinbrink, “Non-invasive detection of fluorescence from exogeneous chromophores in the adult human brain,” NeuroImage 31, 600–608 (2006). [CrossRef] [PubMed]
6. R. Cubeddu, D. Comelli, C. D’Andrea, P. Taroni, and G. Valentini, “Time-resolved fluorescence imaging in biology and medicine,” J. Phys. D: Appl. Phys. 35, R61–R76 (2002). [CrossRef]
7. V. Ntziachristos, C. Bremer, and R. Weissleder, “Fluorescence imaging with near-infrared light: new technological advances that enable in vivo molecular imaging,” Eur. Radiol. 13, 195–208 (2002).
8. J. R. Lakowicz, Principles of Fluorescence Spectroscopy (Springer, New York, 2006). [CrossRef]
9. M. Patterson and B. Pogue, “Mathematical model for time-resolved and frequency-domain fluorescence spectroscopy in biological tissues,” Appl. Opt. 33, 1963–1974 (1994). [CrossRef] [PubMed]
10. C. L. Hutchinson, J. R. Lakowicz, and E. Sevick-Muraca, “Fluorescence lifetime-based sensing in tissues: A computational study,” NeuroImage 68, 1574–1582 (1995).
11. D. E. Hyde, T. J. Farrel, M. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved fluorescence from depth-dependent fluorophore concentrations,” Phys. Med. Biol. 46, 369–383 (2001). [CrossRef] [PubMed]
12. M. Sadoqi, P. Riseborough, and S. Kumar, “Analytical models of time resolved fluorescence spectroscopy in tissues,” Phys. Med. Biol. 46, 2725–2743 (2001). [CrossRef] [PubMed]
13. D. Hattery, V. Chernomordik, M. Loew, I. Gannot, and A. Gandjbakhche, “Analytical solutions for time-resolved imaging in a turbid medium such as tissue,” Opt. Express 16, 13188–13202 (2001).
14. Y. Lu, B. Zhu, H. Shen, J. C. Rasmussen, G. Wang, and E. M. Sevick-Muraca, “A parallel adaptive finite element simplified spherical harmonics approximation solver for frequency domain fluorescence molecular imaging,” Phys. Med. Biol. 55, 4625–4645 (2010). [CrossRef] [PubMed]
15. J. Swartling, A. Pifferi, A. M. K. Enejder, and S. Andersson-Engels, “Accelerated Monte Carlo models to simulate fluorescence spectra from layered tissues,” J. Opt. Soc. Am. A 20, 714–727 (2003). [CrossRef]
16. A. Liebert, H. Wabnitz, N. Żołek, and R. Macdonald, “Monte Carlo algorithm for efficient simulation of time-resolved fluorescence in layered turbid media,” Opt. Express 16, 13188–13202 (2008). [CrossRef] [PubMed]
17. A. D. Klose, V. Ntziachristos, and A. H. Hielscher, “The inverse source problem based on the radiative transfer equation in optical molecular imaging,” J. Comput. Phys. 202, 323–345 (2005). [CrossRef]
18. M. A. O’Leary, D. A. Boas, X. D. Li, B. Chance, and A. G. Yodh, “Fluorescence lifetime imaging in turbid media,” Opt. Lett. 21, 158–160 (1996). [CrossRef]
19. V. Ntziachristos and R. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation,” Opt. Lett. 26, 893–895 (2001). [CrossRef]
20. X. D. Li, M. A. O’Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Fluorescent diffuse photon density waves in homogeneous and heterogenous turbid media: analytic solutions and applications,” Appl. Opt. 35, 3746–3758 (1996). [CrossRef] [PubMed]
21. F. Martelli, S. Del Bianco, A. Ismaelli, and G. ZaccantiLight Propagation through Biological Tissue and Other Diffusive Media: Theory, Solutions, and Software (SPIE, Bellingham, 2010). [CrossRef] [PubMed]
22. A. B. Milstein, S. Oh, K. J. Webb, C. A. Bouman, Q. Zhang, D. A. Boas, and R. P. Millane, “Fluorescence optical diffusion tomography,” Appl. Opt. 42, 3081–3094 (2003). [CrossRef] [PubMed]
23. H. Jiang, “Frequency-domain fluorescent diffusion tomography: a finite-element-based algorithm and simulations,” Appl. Opt. 37, 5337–5343 (1998). [CrossRef]
24. L. Zhang, G. Gao, H. He, and H. Zhao, “Three-dimentional scheme for time-domain fluorescence molecular tomography based on Laplace transforms with noise-robust factors,” Opt. Express 16, 7214–7222 (2008). [CrossRef] [PubMed]
25. A. D. Klose and A. H. Hielscher, “Fluorescence tomography with simulated data based on the equation of radiative tranfer,” Opt. Lett. 28, 1019–1021 (2003). [CrossRef] [PubMed]
26. B. Wassermann, “Limits of high-order perturbation theory in time-domain optical mammography,” Phys. Rev. E 74, 031908 (2006). [CrossRef]
27. D. Grosenick, A. Kummrow, R. Macdonald, P. M. Schlag, and H. Rinneberg, “Evaluation of higher-order time-domain perturbation theory of photon diffusion on breast equivalent phantoms and optical mammograms,” Phys. Rev. E 76, 061908 (2007). [CrossRef]
28. A. Sassaroli, F. Martelli, and S. Fantini, “Perturbation theory for the diffusion equation by use of the moments of the generalized temporal point-spread function. I. Theory,” J. Opt. Soc. Am. A 48, 2105–2118 (2006). [CrossRef]
29. T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. 73, 076701 (2010). [CrossRef]
30. A. Kienle and M. S. Patterson, “Improved solution of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium,” J. Opt. Soc. Am. A 14, 246–524 (1997). [CrossRef]
31. V. Venugopalan, J. S. You, and B. J. Tromberg, “Radiative transport in the diffusion approximation: An extention for highly absorbing media and small source-detector separations,” Phys. Rev. E 58, 2395–2407 (1998). [CrossRef]
32. R. Graaff and K. Rinzema, “Practical improvements on photon diffusion theory: application to isotropic scattering,” Phys. Med. Biol. 46, 3043–3050 (2001). [CrossRef] [PubMed]
33. M. Chu, K. Vishwanath, A. D. Klose, and H. Dehghani, “Light transport in biological tissue using three-dimensional frequency-domain simplified spherical harmonics equations,” Phys. Med. Biol. 54, 2493–2509 (2009). [CrossRef] [PubMed]
34. A. Liemert and A. Kienle, “Analytical solutions of the simplified spherical harmonics equations,” Opt. Lett. 35, 3507–3509 (2010). [CrossRef] [PubMed]
35. P. N. Reinersman and K. L. Carder, “Hybrid numerical method for solution of the radiative transfer equation in one, two, or three dimensions,” Appl. Opt. 43, 2734–2743 (2004). [CrossRef] [PubMed]
36. T. Tarvainen, M. Vauhkonen, V. Kolehmainen, S. R. Arridge, and J. Kaipio, “Coupled radiative transfer equation and diffusion approximation model for photon migration in turbid medium with low-scattering and non-scattering regions,” Phys. Med. Biol. 50, 4913–4930 (2005). [CrossRef] [PubMed]
37. F. Martelli, A. Sassaroli, A. Pifferi, A. Torricelli, L. Spinelli, and G. Zaccanti, “Heuristic Green’s function of the time dependent radiative transfer equation for a semi-infinite medium,” Opt. Express 15, 18168–18175 (2007). [CrossRef] [PubMed]
38. J. C. J. Paasschens, “Solution of the time-dependent Boltzmann equation,” Phys. Rev. E 56, 1135–1141 (1997). [CrossRef]
39. M. Machida, G. Y. Panasyuk, J. C. Schotland, and V. A. Markel, “The Greens function for the radiative transport equation in the slab geometry,” J. Phys. A: Math. Theor. 43, 065402 (2010). [CrossRef]
40. A. Liemert and A. Kienle, “Analytical solution of the radiative transfer equation for the infinite-space fluence,” Phys. Rev. A 83, 015804 (2011). [CrossRef]
41. A. Sassaroli, F. Martelli, and S. Fantini, “Higher-order perturbation theory for the diffusion equation in heterogeneous media: application to layered and slab geometries,” Appl. Opt. 48, D62–D73 (2009). [CrossRef] [PubMed]
42. A. Sassaroli, F. Martelli, and S. Fantini, “Perturbation theory for the diffusion equation by use of the moments of the generalized temporal point-spread function. III. Frequency-domain and time-domain results,” J. Opt. Soc. Am. A 27, 1723–1742 (2010). [CrossRef]
43. A. Sassaroli, C. Blumetti, F. Martelli, L. Alianelli, D. Contini, A. Ismaelli, and G. Zaccanti, “Monte Carlo procedure for investigating light propagation and imaging of highly scattering media,” Appl. Opt. 37, 7392–7400 (1998). [CrossRef]
