Perturbative forward solver software for small localized fluorophores in tissue.

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 r, and that the fluorescent light at wavelength λ, 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 μ, μ′, and n, and at λ as μ, μ′, and n. We assume that a fluorophore is distributed inside a volume V′ with absorption coefficient and reduced scattering coefficient μ and μ′ at λ, and μ and μ′ at λ.

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(μ′ + μ′)] and D = 1/[3(μ′ + μ′)] are the diffusion coefficients at excitation and emission wavelengths. Equations (1) and (2) are coupled by the excitation fluence, Φ, which determines the source term q(r,t) of Eq. (2). In particular, the probability dP that fluorescent photons are generated in the time interval dt′ at time t′, in the wavelength interval dλ around λ, and in the volume element d3r′ at r′, is where η(λ,λ) is the quantum yield or quantum efficiency per unit wavelength for emission at λ 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 η→ = η(λ,λ)dλ quantum efficiency of the fluorophore in the wavelength range dλ around λ, i.e. the ratio of the number of photons emitted in the wavelength range dλ to the number of absorbed photons. According to Green’s function theorem, the fluence rate at emission wavelength, Φ (r,t), due to the fluorescent photons generated in the interval dλ, is with G Green’s function of the medium for the fluence at the emission wavelength. We stress that the Green’s functions G and G 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 Φ0 with the exponential decay of the fluorophore

It 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 Φ0 (r,t) and Φ (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 G 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 G 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 G 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. q0 and q0 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) = Φ(r)+Φ(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 Φ0 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 Φ0 can be obtained by using the CW solution, provided that the coefficient μ is replaced with the term [μ – (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/μ′). The absorption coefficient of the fluorophore at the excitation wavelength is μ = 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 μ = 0.002 mm−1. The refractive index of the medium and of the external is 1.4.

Fig. 3

Identical to Fig. 2 except for the Fluorescence CW reflectance.

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 μ of the fluorophore at the excitation wavelength were considered (for each case we have V′μ = 0.01 mm2): V′ = 10 mm3, μ = 0.001 mm−1; V′ = 100 mm3, μ = 0.0001 mm−1; V′ = 500 mm3, μ = 0.00002 mm−1, and V′ = 1000 mm3, μ = 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 μ, the RTE hybrid model provides the same results since Eq. (11) depends only on the product V′μ 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′μ = 0.01 mm2): V′ = 10 mm3, μ = 0.001 mm−1; V′ = 100 mm3, μ = 0.0001 mm−1; V′ = 500 mm3, μ = 0.00002 mm−1, and V′ = 1000 mm3, μ = 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 G 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.