Analytical reconstruction of the bioluminescent source with priors

Hakan Erkol; Aytac Demirkiran; Nasire Uluc; Mehmet B. Unlu

doi:10.1364/OE.22.019758

1. Introduction

Excitement over the development of stand-alone bioluminescence tomography systems has been decreasing for the last few years. This has been partly due to the frustration caused by the lack of high resolution images. Unfortunately, the multiple scattering of near infrared photons inside tissue results poor images. Another contributing factor is the non-uniqueness of the inverse source problem for measurements done in one wavelength. Increasing the number of wavelengths compromises the simplicity of the tomographic system and increases the measurement time as well as the duration of the image reconstruction. As a result, time decay of the bioluminescent source becomes an issue to be taken into account in small animal studies. The current direction is the development of multimodal systems. Combined X-ray computed tomography and bioluminescence imaging systems and sometimes with the addition of diffuse optical tomography component have been built in the literature [1–5].

The structural information obtained from the computed tomography is routinely used to increase the spatial resolution and contrast of the bioluminescence tomography images. Whereas the diffuse optical tomography can provide the functional priors needed for quantification.

Photon transport in a diffusive medium can be described exactly by the radiative transfer equation [6, 7]. However, since the radiative transfer equation is very difficult to solve analytically, it is generally reduced to the diffusion equation [8, 9]. It is well known that tissue is a diffusive medium so that light propagation in tissue can be described well by the diffusion equation [10].

There are several optical imaging techniques dealing with the problem of light diffusion in tissue such as diffuse optical tomography (DOT) and bioluminescence tomography (BLT). In DOT, known external light sources are employed for the purpose of reconstruction of optical parameters of a diffusive medium [11, 12] whereas in BLT, unknown internal bioluminescent light distributions in a diffusive medium are tried to be quantified from boundary measurements [13–17].

Mathematical framework for DOT and BLT involves two conventional steps called as the forward and inverse problems [18]. The forward problem of both DOT and BLT aims to calculate the photon fluence rate within the tissue and the outward flux through the tissue boundary for a given source distribution and tissue optical properties. On the other hand, the inverse problem utilizes the forward problem to reconstruct the distributions of the optical properties of the tissue in DOT and the unknown bioluminescent source within the tissue in BLT using the measured boundary flux. Therefore, both DOT and BLT demand to solve the forward problem both quickly and accurately.

It is very difficult to obtain exact analytical solutions to the diffusion equation describing light propagation in tissue because of the structural inhomogeneities and irregular geometries. However, analytical methods to the forward problem, which are more readily applicable in practice, can be used for some definite geometries such as infinite or semi-infinite slabs, spherical, parallel-epiped, circular and cylindrical regions [10, 18–34]. Analytic solutions to the diffusion equation for some layered geometries are also found in which all layers are homogeneous [21, 29, 31, 35–39]. It is important to note that studies of light diffusion in layered geometries are very important since many imaging applications consist of layered structures [21, 40–47]. In addition, in BLT, bioluminescent source is inside the tissue and the detectors are placed on the boundary surface. If the tissue is very thin or long enough, then 3 dimensional polar coordinates can be reduced (or approximated) to 2 dimensional polar coordinates. In the Diffuse optical tomography, the detectors are also placed on the boundaries of a circular tissue domain and the approximation is also two dimensional. Therefore, these regular models can generate insight into the complex geometric models of tissues.

Solutions to the diffusion equation using Green’s function method for various homogenous geometries including slabs, cylinders, spheres for both the time and frequency domains can be found in the literature [20, 22]. Arridge et al. solved the diffusion equation for an infinite line source in an infinite cylinder and a point source in a sphere in the frequency domain using the Green’s function method [22]. Sikora et al. solved the diffusion equation for concentric spheres by using a series expansion method [38]. Kienle et al. obtained an analytical solution to the diffusion equation in the continuous wave, frequency and time domains for a multiple layered finite cylinder. Kienle et al. followed the Green’s function methodology by applying the extrapolated boundary conditions [39]. In addition, Zhang et al. solved the steady state equation in cylindrical coordinates for a point source by using modified bessel functions and applying the extrapolated boundary conditions [33].

We solve the diffusion equation analytically in a novel way for a point source represented by the Dirac delta function with the strength of α. We first employ detector readings and then combine them with the Robin boundary condition. Next, we utilize the continuities of the photon fluence rate and its derivative at the surface of the inner region to determine the photon fluence rate between the inner and the outer regions. Later, by making use of the continuity of the photon fluence rate and the discontinuity of its derivative at the source position, we obtain the photon density for the inner region in terms of the boundary measurements, the diffusion coefficients, and the associated Bessel functions.

In seeking analytical solutions to the forward problem, the source term is generally assumed to be the point source of unit strength in the literature mentioned above. However, determining the unknown strength of internal bioluminescent point sources using boundary measurements is essential in BLT applications. As another novelty of this work, the source term is not of unit strength but some value described by α. We concentrate on the determination of the strength as well as the location of the source inside a homogeneous and inhomogeneous (a two-layered) concentric circular turbid media with different optical properties by using the detector readings performed on the outer surface. As a result, we obtain analytical expressions for the strength in terms of the source position, detector readings and the diffusion coefficients.

In BLT, if both of the strength and position of the internal source are unknowns, they cannot be determined with good accuracy from boundary measurements due to the ill-posed nature of the problem. Numerical methods applied in literature are also in need of a priori knowledge. This a priori knowledge is achieved from other imaging modalities providing spatial priors (e.g. the exact position of the internal source) such as magnetic resonance imaging, X-ray micro computed tomography or ultrasound computed tomography [4, 16, 48]. BLT together with some of the imaging system providing structural information can be utilized as a multimodality system.

Our calculations are based on the fact that if source strength is known a priori, then source location could be calculated or vice versa: if source location is known a priori, source strength could be calculated. For this reason, firstly; the location (or strength) of the source has to be determined. Furthermore, it is also crucial that in our analytical calculations, a single wavelength measurement is sufficient to obtain an accurate source reconstruction regardless of the uniqueness problem in BLT inverse source problem [1, 14]. Hence, this method may also be used to reduce the imaging time.

2. Method

2.1. Theoretical method

The diffusion equation in time domain is given by the following partial differential equation [9]

\frac{\partial Φ (r, t)}{c \partial t} + μ_{a} Φ (r, t) - \nabla \cdot [D \nabla Φ (r, t)] = S (r, t)

where Φ, c, D, μ_a, and S represent the photon fluence rate at position r and time t, the speed of light, the photon diffusion coefficient, the absorption coefficient and source distributions of photons, respectively. Here, the diffusion coefficient is given as

D (r) = \frac{1}{3 [μ_{a} (r) + μ_{s}^{'} (r)]}

where μ′_s is the reduced scattering coefficient.

Using the Fourier transform, Φ(r, t) = ∫Φ(r, ω)exp(iωt)dω, the diffusion equation in frequency domain takes the following form

(\frac{i ω}{c} + μ_{a} - D \nabla^{2}) Φ (r, ω) = S (r, ω)

where D is assumed to be constant.

If the source term is treated as a point source represented by the Dirac delta function, the Fourier transform of the diffusion equation in frequency domain can be written as

\nabla^{2} Φ (r, ω) - (\frac{i ω}{c} + μ_{a}) \frac{1}{D} Φ (r, ω) = - \frac{1}{D} α_{i} δ (r - r_{i})

where r_i and α_i are the position and the strength of the point source, respectively. For bioluminescence applications, ω is taken to be zero so that Eq. (3) reduces to

\nabla^{2} Φ (r) - \frac{μ_{a}}{D} Φ (r) = - \frac{1}{D} α_{i} δ (r - r_{i}) .

2.1.1. Light diffusion in a circular homogenous tissue

From now on, we focus on two dimensional cylindrical coordinates and make a change of variable r → ρ. In Fig. 1, we consider a polar tissue domain in which a point source is placed at (ρ = ρ_i, θ = 0). In other words, it is located in the x axis a distance of ρ_i away from the origin. The center of the circular domain is chosen as the origin and there is no any other light source except from the point source. Here, the radius and the diffusion coefficient of the domain are represented by R and D, respectively.

Fig. 1 The circular region with an inclusion represented by a point source located at ρ_i.

Download Full Size | PDF

For polar coordinates, Eq. (4) takes the following form

\nabla^{2} Φ (ρ, ϕ) - \frac{μ_{a}}{D} Φ (ρ, ϕ) = - \frac{1}{D} α_{i} δ (ρ - ρ_{i}, ϕ - ϕ_{i})

where ρ_i and α_i are the position and the strength of the point source, respectively. If (ρ, ϕ) ≠ (ρ_i, ϕ_i), assuming that Φ(ρ, ϕ) = Φ(ρ, −ϕ) (as explained in the discussion and results part, the assumption is quite realistic), the solution of Eq. (5) is

Φ (ρ, ϕ) = \sum_{m = 0}^{\infty} (A_{m} J_{m} (k ρ) + B_{m} Y_{m} (k ρ)) cos (m ϕ)

where

k = \sqrt{- \frac{μ_{a}}{D}}

, J_m and Y_m are the Bessel functions of the first kind and the second kind, respectively.

Now we consider a detector located on the boundary of the tissue domain to obtain a solution near the detector. At the location of the detector ρ = R, defining Φ(ρ, ϕ) ≡ Φ₁(ρ, ϕ), we take the general solution given by Eq. (6) between the boundary of the domain and the point source in Fig. 1, i.e. r ∈ [ρ_i, R]. Hence, the transformation of the measured data Φ₁(R, t) in the continuous wave domain (ω = 0) [49], Φ₁(kR, ϕ), can be written as

Φ_{1} (k R, ϕ) = \sum_{m = 0}^{\infty} (A_{m} J_{m} (k R) + B_{m} Y_{m} (k R)) cos (m ϕ) .

Multiplying Eq. (7) with cos(m′ϕ) and then integrating over ϕ between 0 and 2π yields

A_{m} J_{m} (k R) + B_{m} Y_{m} (k R) = \frac{1}{π} \int_{0}^{2 π} Φ_{1} (R, ϕ) cos (m ϕ) d ϕ

where the orthonormality relation of cosine function is used

\int_{0}^{2 π} cos (m ϕ) cos (m^{'} ϕ) d ϕ = π δ_{m, m^{'}}

for m, m′ ≠ 0. If m and m′ are zero, the integral in Eq. (9) is equal to 2π. To obtain A_m and B_m, the Robin boundary condition is introduced [8, 50, 51]:

Φ_{1} (R) + {2 ζ D \frac{\partial Φ_{1} (ρ)}{\partial ρ} |}_{ρ = R} = 0

where ζ and D are constants representing the internal reflection of light due to the index of refraction mismatch for tissue-air and the diffusion coefficient, respectively. Substituting Φ₁(kR, ϕ) into Eq. (10) yields

A_{m} = \frac{R}{2} I_{m} (\frac{Y_{m} (k R)}{2 ζ D} + Y_{m}^{'} (k R))

and

B_{m} = - \frac{R}{2} I_{m} (\frac{J_{m} (k R)}{2 ζ D} + J_{m}^{'} (k R))

where

I_{m} = \int_{0}^{2 π} Φ_{1} (R, ϕ) cos (m ϕ) d ϕ,

prime stands for differentation with respect to ρ and the Wronskian of the Bessel functions,

W (J_{m} (k R), Y_{m} (k R)) = \frac{2}{π R}

, is used.

The photon fluence rate between the origin and the source, ρ ∈ [0, ρ_i], is

Φ_{2} (ρ, ϕ, k) = \sum_{m = 0}^{\infty} (C_{m} J_{m} (k ρ) + D_{m} Y_{m} (k ρ)) cos (m ϕ)

but since Y_m(kρ) blows up as ρ goes to 0, D_m must be zero so that

Φ_{2} (ρ, ϕ, k) = \sum_{m = 0}^{\infty} C_{m} J_{m} (k ρ) cos (m ϕ) .

The continuity of the photon fluence rate at ρ = ρ_i (the source position) is

{Φ_{2} (ρ, ϕ, k) |}_{ρ = ρ_{i}} = {Φ_{1} (ρ, ϕ, k) |}_{ρ = ρ_{i}} .

For a point source, the diffusion equation for the inner region in Fig. 1 is

\begin{array}{r} \sum_{m = 0}^{\infty} (\frac{1}{ρ} \frac{d}{d ρ} (ρ \frac{d}{d ρ} g_{m} (k ρ)) cos (m ϕ) + \frac{1}{ρ^{2}} \frac{d^{2}}{d ϕ^{2}} (cos (m ϕ)) g_{m} (k ρ) \\ - \frac{μ_{a}}{D} Φ (k ρ, ϕ)) = - \frac{α_{i}}{D} δ (ρ, ρ_{i}) \end{array}

where g_m(kρ) = [J_m(kρ), Y_m(kρ)] [52–54]. In order to obtain a discontinuity relation for the derivative of the photon fluence rate, we multiply Eq. (17) with ρ cos(m′ϕ) and integrate between (ρ = ρ_i − ε, ϕ = 0) and (ρ = ρ_i + ε, ϕ = 2π) with respect to ρ and ϕ. Considering ε → 0⁺ and applying the continuity of the photon fluence rate gives

\begin{array}{r} \sum_{m = 0}^{\infty} (\int_{ρ_{i} - ε}^{ρ_{i} + ε} (\frac{d}{d ρ} (ρ \frac{d g_{m} (k ρ)}{d ρ})) d ρ \int_{0}^{2 π} cos (m ϕ) cos (m^{'} ϕ) d ϕ \\ = - \frac{1}{D} α_{i} \int_{ρ_{i} - ε}^{ρ_{i} + ε} \int_{0}^{2 π} cos (m^{'} ϕ) δ (ρ, ρ_{i}) ρ d ρ d ϕ \end{array}

where D is assumed to be constant through the domain. Substituting

δ (ρ, ρ_{i}) = \frac{δ (ρ - ρ_{i}) δ (ϕ - ϕ_{i})}{ρ}

and

\int_{0}^{2 π} cos (m ϕ) cos (m^{'} ϕ) d ϕ = π δ_{m m^{'}}

(notice that m, m′ ≠ 0, if both of them are zero, the orthonormality relation is 2π) into Eq. (18) as ε → 0⁺ and utilizing the continuity of the photon fluence rate yields the following discontinuity relation

{\frac{d g_{m} (k ρ)}{d ρ} |}_{ρ = ρ_{i} + ε} {- \frac{d g_{m} (k ρ)}{d ρ} |}_{ρ = ρ_{i} - ε} = - \frac{α_{i}}{π D ρ_{i}} cos (m ϕ_{i}) .

Writing Φ₂(ρ, ϕ, k) as the photon fluence rate between the origin and the point source and Φ₁(ρ, ϕ, k) as the photon fluence rate between the point source and the boundary of the domain (see Fig. 1) into Eqs. (16) and (21) yields

A_{m} J_{m} (k ρ_{i}) + B_{m} Y_{m} (k ρ_{i}) = C_{m} J_{m} (k ρ_{i})

and

A_{m} J_{m}^{'} (k ρ_{i}) + B_{m} Y_{m}^{'} (k ρ_{i}) = C_{m} J_{m}^{'} (k ρ_{i}) - \frac{α_{i}}{π D ρ_{i}} cos (m ϕ_{i}) .

Solving Eqs. (22) and (23) for C_m and α_i leads to

C_{m} = A_{m} + B_{m} \frac{Y_{m} (k ρ_{i})}{J_{m} (k ρ_{i})}

and

α_{i} = \frac{{DRI}_{m}}{J_{m} (k ρ_{i}) cos (m ϕ_{i})} (\frac{J_{m} (k R)}{2 ζ D} + J_{m}^{'} (k R)),

respectively. Note that I_m corresponds to a specific m value which is dependent on the angular profile of photon fluence rate due to the source. For this reason, the strength of the source takes just one physical value corresponding to that m value. In fact, it is impossible in practice to obtain an infinite number of detector measurements on the boundary. Hence, as it will be explained in numerical part, the surface integral, I_m, can be approximated by the summation.

2.1.2. Light diffusion in a circular heterogonous tissue

In Fig. 2, we consider a two-layered concentric tissue domain in which a point source is placed at ρ = ρ_i, between the center of the domain and the inner boundary.

Fig. 2 Concentric two-layered domain.

Download Full Size | PDF

Here, the radii and the diffusion coefficients of the inner and the outer regions are represented by R_in, R_out, D_in and D_out, respectively.

Following the similar steps performed in the previous section, the solutions for the outer region is

Φ_{out} (ρ, ϕ) = \sum_{m = 0}^{\infty} (A_{m} J_{m} (k_{out} ρ) + B_{m} Y_{m} (k_{out} ρ)) cos (m ϕ)

where

k_{out} = \sqrt{- \frac{μ_{a, out}}{D_{out}}}

and R_in ≤ ρ ≤ R_out. Using detector readings and the Robin boundary condition, constants A_m and B_m are expressed in terms of detector readings as it follows

A_{m} = \frac{R_{out}}{2} I_{m} (\frac{Y_{m} (k_{out} R_{out})}{2 ζ_{t a} D_{out}} + Y_{m}^{'} (k_{out} R_{out}))

and

B_{m} = - \frac{R_{out}}{2} I_{m} (\frac{J_{m} (k_{out} R_{out})}{2 ζ_{t a} D_{out}} + J_{m}^{'} (k_{out} R_{out}))

where

I_{m} = \int_{0}^{2 π} Φ (R_{out}, ϕ) cos (m ϕ) d ϕ,

prime stands for differentiation with respect to ρ and ζ_ta is the internal reflection of light due to the index of refraction mismatch for tissue-air.

The solution between the point source and the inner boundary has the following form

Φ_{in} (ρ, ϕ) = \sum_{m = 0}^{\infty} (C_{m} J_{m} (k_{in} ρ) + D_{m} Y_{m} (k_{in} ρ)) cos (m ϕ)

where

k_{in} = \sqrt{- \frac{μ_{a, in}}{D_{in}}}

and 0 ≤ ρ ≤ R_in. Now, in order to determine the constants C_m and D_m, we use the following boundary conditions of Φ and the photon flux −D. ∇Φ at the boundary of the inner layer (ρ = R_in) [55–57]

{{(\frac{n_{2}}{n_{1}})}^{2} Φ_{in} (k_{in} ρ) |}_{ρ = R_{in}} {- Φ_{out} (k_{out} ρ) |}_{ρ = R_{in}} {= c (\frac{n_{2}}{n_{1}}) (- D . \nabla Φ) |}_{ρ = R_{in}}

and

{- D_{in} \frac{\partial Φ_{in}}{\partial ρ} |}_{ρ = R_{in}} = {- D_{out} \frac{\partial Φ_{out}}{\partial ρ} |}_{ρ = R_{in}},

respectively. For tissue-tissue boundary, n₁ ≃ n₂ and

c (\frac{n_{2}}{n_{1}}) ≃ 0

[58,59] so that Eq. (31) becomes

{Φ_{in} (k_{in} ρ) |}_{ρ = R_{in}} = {Φ_{out} (k_{out} ρ) |}_{ρ = R_{in}} .

Solving Eqs. (32) and (33) for C_m and D_m leads to

\begin{array}{r} C_{m} = \frac{π R_{in}}{2} {A_{m} [J_{m} (k_{out} R_{in}) Y_{m}^{'} (k_{in} R_{in}) - \frac{D_{out}}{D_{in}} J_{m}^{'} (k_{out} R_{in}) Y_{m} (k_{in} R_{in})] \\ + B_{m} [Y_{m} (k_{out} R_{in}) Y_{m}^{'} (k_{in} R_{in}) - \frac{D_{out}}{D_{in}} Y_{m}^{'} (k_{out} R_{in}) Y_{m} (k_{in} R_{in})]} \end{array}

and

\begin{array}{r} D_{m} = - \frac{π R_{in}}{2} {A_{m} [J_{m} (k_{out} R_{in}) J_{m}^{'} (k_{in} R_{in}) - \frac{D_{out}}{D_{in}} J_{m}^{'} (k_{out} R_{in}) J_{m} (k_{in} R_{in})] \\ + B_{m} [Y_{m} (k_{out} R_{in}) J_{m}^{'} (k_{in} R_{in}) - \frac{D_{out}}{D_{in}} Y_{m}^{'} (k_{out} R_{in}) J_{m} (k_{in} R_{in})]}, \end{array}

respectively.

From now on, following the same steps performed for the homogeneous case, using the continuity of the Φ and the discontinuity of its derivative at the source position, the strength of the source is obtained

α_{i} = - \frac{2 D_{in}}{J_{m} (k_{in} ρ_{i}) cos (m ϕ_{i})} D_{m} .

Therefore, the strength is turned out to be dependent on D_m and hence detector readings.

2.2. Numerical method

In this section, in order to validate the robustness of the analytical method presented in the previous section, numerical calculations are performed for realistic tissue parameters. The detector readings obtained from the circular boundary are used in the calculations for both homogeneous and heterogeneous domains, respectively. Firstly, the strength of the bioluminescent source is calculated with a known source position, analytically. The change of the source strength with respect to the source position is also investigated. Next, the source position is obtained with a known source strength. The behavior of the position with increasing strength values is examined. The calculated values of the strength and position are compared with the true values. The analytical solutions are obtained based on the detector readings obtained from the boundary. For this reason, the tissue-diffusion equation is solved by using the finite element method (FEM) to introduce synthetic detector measurements. The tissue-diffusion equation is solved on a mesh of 1761 nodes and 3264 first order triangular elements to obtain the photon fluence rate. Figure 3 shows the finite element mesh used in the calculations.

Fig. 3 Finite element mesh used in the simulations.

Download Full Size | PDF

The detector readings are obtained from 128 different points from the boundary of the circular domain. The distance between the successive detector positions is the same through the boundary. In other words, each detector reading is measured for every 360/128= 2.8125 degree.

The surface integral on the boundary given in the previous section, I_m,

I_{m} = \int_{0}^{2 π} Φ_{1} (R, ϕ) cos (m ϕ) d ϕ

states that the strength has to get only one value corresponding to a specific m since it depends on I_m. In other words, the summation over all m values does not give a realistic result. As it will be explained in the next section, the calculations show that the calculated strength and position of the source are very close to their true values for just one m value. Obtaining the strength or position of the source for a specific m value makes it possible to get rid of the non-uniqueness problem. While calculating the strength or position of the source, if both of them are unknowns then it leads to the ill-posedness. For this reason, in our calculations, one of them is calculated on the condition that the other one is known.

Note that, in BLT, it is not possible to obtain an infinite number of detector readings from the boundary. Thus, only a finite number of measurements can be collected so that in our calculations, we write I_m as the following summation over the boundary readings instead of the integral

I_{m} = \int_{0}^{2 π} Φ (R, ϕ) cos (m ϕ) d ϕ = \sum_{i = 1}^{128} Φ (R, ϕ_{i}) cos (m ϕ_{i}) .

Note that the distance between the successive detector measurements is s = Rθ = 2πR/128 and the first measurement, i = 1 is chosen such that the corresponding detector position is on the x axis at the boundary of the domain, (ρ = R, θ = 0).

3. Discussion and results

Our numerical results show that when m = 1, the calculated strength and position values are very close to the true values. When m ≠ 1, the calculated results are not comparable with the true ones. Hence, based on the detector readings, the photon transport in tissue is described well for m = 1.

As well as the analytical calculations, 2D simulations are also performed using the finite element method. For the point source approximation, the radius of the circular boundary is taken 30 times higher than the radius of the point source. Here, the diameter of the circular domain is 40 mm. The tissue-diffusion equation parameters; absorption and reduced scattering coefficients are chosen as μ_a = 0.00281 mm⁻¹ and μ′_s = 1.6667 mm⁻¹ [12, 60], respectively.

In the theoretical method, the tissue-diffusion equation was solved based on the assumption Φ(ρ, ϕ) = Φ(ρ, −ϕ) leading to cos(mϕ) solutions. This approximation is quite realistic because in our calculations, the radius of domain is greater than the radius of the source so that the source can be treated as a point source. If the point source is located on the x-axis, the photon densities for the upper and lower regions of the circular domain are expected to be the same. The results of simulations show that they are really the same. These results obtained by FEM can be seen in Fig. 4 showing the distribution of the photon fluence rate along the measurement points on the boundary for the homogeneous case. Here, the photon fluence rate due to the point source located on the x-axis has its maximum value so that the measurement point corresponding to this maximum reading is set as the first measurement.

Fig. 4 Photon fluence rate (number of photons per mm⁻²) vs. number of measurement for a point source located x = 10 mm and y = 0 where the radius of the circular domain R=20 mm.

Download Full Size | PDF

To obtain the strength and position, the Robin boundary condition is used and the constant ζ representing the internal reflection of light due to the index of refraction mismatch is taken as 2.74 [12, 60].

For the homogeneous domain, the strength given in the theoretical part is

α_{m} = \frac{{DRI}_{m}}{J_{m} (k ρ_{i}) cos (m ϕ_{i})} [\frac{J_{m} (k R)}{2 ζ D} + J_{m}^{'} (k R)]

where

I_{m} = \sum_{i = 1}^{128} Φ (R, ϕ_{i}) cos (m ϕ_{i})

, R is the radius of the circular domain, (ρ_i, ϕ_i) is the position of the point source and prime stands for differentiation with respect to ρ. We perform our calculations for m = 1 among the possible m = 0, 1, 2... values since m = 1 value gives a very good realistic result compared with the others. In our calculations, we take the radius of the circular domain as R = 20 mm. After setting the true value of the source and obtaining the detector readings on the boundary at 128 different angles using the FEM, we insert these measurements in I_m and hence in Eq. (39). First of all, we calculate the strength α at various source positions inside the domain and show our results in Fig. 5.

Fig. 5 Source strength α vs. source position ρ where ⋄ and + represent the calculated value and the true value, respectively; the radius of the circular domain, R = 20 mm and the true value of the source strength, α_true = 1mm⁻¹.

Download Full Size | PDF

Figure 5 shows that the calculated value is fluctuating very slightly around the true value. On the other hand, while getting closer to the boundary or origin, the tissue-diffusion approximation begins to lose its feasibility to describe the photon transport in the tissue, and hence; I_m obtained from the numerical values on the boundary leads to a larger difference between the measurement based analytical result and the true value. Therefore, the calculated results are in a very good accordance with the true ones unless the source is very close to the origin or the boundary as expected. We also obtain the position of the source by solving Eq. (39) for various strength values. Before solving Eq. (39), we obtain detector readings and then I_m value by setting the true value of strength. We plot the position with respect to the strength in Fig. 6. The calculated position is around 9.67 mm and very close to the true value. The calculated position changes very slightly with the strength. The reason for this accuracy is that the tissue-diffusion equation works well at the position of the source. When the source is located near the origin or boundary, likewise the calculated strength, the calculated value of the position deviates considerably from its true value as expected.

Fig. 6 Position of the source ρ vs. strength of the source α for the radius of the circular domain, R = 20 mm where + and ⋄ represent the calculated and the true value, respectively.

Download Full Size | PDF

For the heterogenous domain, similar numerical calculations are performed to obtain the strength given by Eq. (36) in the theoretical part, $α_{i} = - \frac{2 D_{in}}{J_{m} (k_{in} ρ_{i}) cos (m ϕ_{i})} D_{m}$ . Here, D_m is given by Eq. (35) in terms of A_m and B_m or the detector measurements.

Here, likewise the homogeneous case, our numerical calculations show that the calculated results are very close to the true values for m = 1 case. Therefore, Φ must have one specific m value in order that the solution to be physical for the heterogeneous case.

The heterogeneous domain consists of two different concentric circular regions. The absorption and the reduced scattering coefficients for the inner region are μ_a,in = 0.00281 mm⁻¹ and μ_s′,in = 1.6667 mm⁻¹ [12,60], respectively. For the outer tissue region, the parameters are taken as μ_a,out = 3μ_a/4 mm⁻¹ and μ_s′,out = 3μ_s′/4 mm⁻¹. When the Robin boundary condition is used for the outer boundary, the mismatch index of tissue-air is taken as ζ = 2.74 [12, 60, 61]. For biological applications, the differences between the diffusion parameters for the inner and the outer regions lead a negligible change in the refractive index so that we take the same refractive index for both regions, n_in = n_out.

Likewise the homogeneous case, the point source is located on the x-axis. Figure 7 shows the distribution of the photon fluence rate along the measurement points on the outer boundary obtained by FEM.

Fig. 7 Photon fluence rate (number of photons per mm⁻²) vs. number of measurement for a point source located x = 10 mm and y = 0 where the radii of the outer and inner circular domain R=20, 15mm, respectively.

Download Full Size | PDF

The measurement point corresponding to the maximum reading is set as the first measurement. In Fig. 8, for various source positions, we calculated the strength of the source for R_in = 15 mm and R_out = 20 mm.

Fig. 8 Source strength α vs. source position ρ for the layered domain. Here, ⋄ and + represent the calculated values and the true value, respectively; the radii of the outer and inner domain R = 20, 15mm and the true value of the source strength α_true = 1mm⁻¹.

Download Full Size | PDF

When the source is getting closer to the origin, the difference between the calculated and the true value is increasing because the detector readings obtained from the numerical solutions of the diffusion equation are getting worse. Thus, this deterioration affects the analytical results. On the other hand, when the source is very close to the inner layer, this time the boundary conditions coincide with the discontinuity relation of the Dirac delta function. For this reason, the calculated vale of the strength deviates from its true value considerably since the method does not work when source position is in the vicinity of the inner boundary. The method works very well unless the source position is very close to the inner boundary or the origin. In Fig. 9, we also calculate the source position for various strength values and see that it is around 8.85 mm and changing very slightly with respect to the strength.

Fig. 9 Source position ρ vs. source strength α for the layered circular domain. Here, the outer and inner radii R = 20, 15 mm where ⋄ and + represent the calculated and the true value, respectively.

Download Full Size | PDF

The true value of the position is 10 mm. As it can be seen from Figs. 8 and 9, the method works well in the heterogeneous medium. When the calculated position values for heterogeneous case are compared with the homogeneous case’s results (Figs. 6 and 9), it can be easily seen that the latter approaches the true value better as expected.

Therefore, we have presented an analytical reconstruction approach to obtain the strength or position of a bioluminescence point source described by the Dirac delta function for homogeneous and heterogeneous circular turbid media. Main goals of BLT are to determine the strength and location of an internal source. These goals are generally achieved by means of numerical techniques. In spite of the fact that numerical techniques can be used for complex geometries, they are computationally expensive in practice. On the other hand, analytical approaches can positively affect the computational cost and give more accurate results compared with the numerical approaches. Although these approaches are useful for just some definite geometries, they can generate insight into more complex geometries. For example, a sufficiently long cylindrical medium can be approximated by a 2D polar coordinate system. This kind of system can be investigated by utilizing the methodology presented in this work. Therefore, our method offers a different analytical point of view to BLT. Because of the fact that our method is based on an analytical approach, it gives more accurate results compared with some other numerical methods (FEM or FDM). The method can save computing time since it does not require any simulations.

Acknowledgments

This research is supported in part by Marie Curie Reintegration Grant 268287, Bogazici University Research funding BAP 7126 and BAP 6033, and TUBITAK Grant 112T253.

References and links

1. G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,” Med. Phys. 31(8), 2289–2299 (2004). [CrossRef] [PubMed]

2. C. Darne, Y. Lu, and E. M. Sevick-Muraca, “Small animal fluorescence and bioluminescence tomography: a review of approaches, algorithms and technology update,” Phys. Med. Biol. 59, R1–R64 (2014). [CrossRef]

3. J. A. Guggenheim, H. R. A. Basevi, J. Frampton, I. B. Styles, and H. Dehghani, “Multi-modal molecular diffuse optical tomography system for small animal imaging,” Meas. Sci. Technol. 24, 105405 (2013). [CrossRef] [PubMed]

4. H. Yan, Y. Lin, W. C. Barber, M. B. Unlu, and G. Gulsen, “A gantry-based tri-modality system for bioluminescence tomography,” Rev. Sci. Instrum. 83, 043708 (2012). [CrossRef] [PubMed]

5. A. Corlu, R. Choe, T. Durduran, K. Lee, M. Schweiger, S. R. Arridge, E. M. C. Hillman, and A. G. Yodh, “Diffuse optical tomography with spectral constraints and wavelength optimization,” Appl. Opt. 44(11), 2082–2093 (2005). [CrossRef] [PubMed]

6. A. Ishimaru, Wave Propagation and Scattering in Random Media (Academic, 1978), Vol 1.

7. W. Han, J. A. Eichholz, and G. Wang, “On a family of differential approximations of the radiative transfer equation,” J. Math. Chem. 50, 689–702 (2011). [CrossRef]

8. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, 41–93 (1999). [CrossRef]

9. L. H. V. Wang and H. Wu, Biomedical Optics: Principles and Imaging Hoboken (Wiley, 2007).

10. M. S. Patterson and S. J. Madsen, “Diffusion equation representation of photon migration in tissue,” IEEE MTT-S Digest 2, 905–908 (1991).

11. M. Schweiger and S. R. Arridge, “Optical tomographic reconstructon in a complex head model using a priori region boundary condition,” Phys. Med. Biol. 44, 2703–2721 (1999). [CrossRef] [PubMed]

12. M. B. Unlu, O. Birgul, R. Shafiha, G. Gulsen, and O. Nalcıoğlu, “Diffuse optical tomographic reconstruction using multifrequency data,” J. Biomed. Opt. 11, 054008 (2006). [CrossRef]

13. G. Wang, W. Cong, K. Durairaj, X. Qian, H. Shen, P. Sinn, E. Hoffman, G. Mclennan, and M. Henry, “In vivo mouse studies with bioluminescence tomography,” Opt. Express 14, 7801–7809 (2006). [CrossRef] [PubMed]

14. H. Dehghani, S. C. Davis, S. Jiang, B. W. Pogue, K. D. Paulsen, and M. S. Patterson, “Spectrally resolved bioluminescence optical tomography,” Opt. Lett. 31, 365–367 (2006). [CrossRef] [PubMed]

15. H. Dehghani and B. W. Pogue, “Spectrally resolved bioluminescence optical tomography using the reciprocity approach,” Med. Phys. 35, 4863–4871 (2008). [CrossRef] [PubMed]

16. M. A. Naser and M. S. Patterson, “Algorithms for bioluminescence tomography incorporating anatomical information and reconstruction of tissue optical properties,” Biomed. Opt. Express 1, 512–526 (2010). [CrossRef]

17. Q. Zhang, L. Yin, Y. Tan, Z. Yuan, and H. Jiang, “Quantitative bioluminescence tomography guided by diffuse optical tomography,” Opt. Express 16, 1481–1486 (2008). [CrossRef] [PubMed]

18. S. R. Arridge, “A finite element approach for modelling photon transport in tissue,” Med. Phys. 20, 299–309 (1993). [CrossRef] [PubMed]

19. M. Schweiger, S. R. Arridge, and D. T. Delpy, “Application of the finite element method for the forward and inverse model in optical tomography,” J. Mat. Imaging Vis. 3, 263–283 (1993). [CrossRef]

20. M. S. Patterson, “Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties,” Appl. Opt. 28, 2331–2336 (1989). [CrossRef] [PubMed]

21. J. M. Schmitt, G. X. Zhou, E. C. Walker, and R. T. Wall, “Multilayer model of photon diffusion in skin,” J. Opt. Soc. Am. A 7, 2141–2153 (1990). [CrossRef] [PubMed]

22. S. R. Arridge, M. Cope, and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol. 37, 1531–1560 (1992). [CrossRef] [PubMed]

23. T. J. Farrel, M. S. Patterson, and B. C. Wilson, “A diffusion theory model of specially resolved steady state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys. 19, 879–888 (1992). [CrossRef]

24. D. A. Boas, M. A. O’leary, B. Chances, and A. G. Yodh, “Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: Analytic solution and applications,” Proc. Natl. Acad. Sci. U.S.A. 91, 4887–4891 (1994). [CrossRef] [PubMed]

25. S. A. Walker, D. A. Boas, and E. Gratton, “Photon density waves scattered from cylindrical inhomogeneities: theory and experiments,” Appl. Opt. 37, 1935–1944 (1998). [CrossRef]

26. B. W. Pogue and M. S. Patterson, “Frequency-domain optical absorption spectroscopy of finite tissue volumes using diffusion theory,” Phys. Med. Biol. 39, 1157–1180 (1994). [CrossRef] [PubMed]

27. D. Contini, F. Martelli, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation I. Theory,” Appl. Opt. 36, 4587–4599 (1997). [CrossRef] [PubMed]

28. A. Kienle and M. S. Patterson, “Improved solutions 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–254 (1997). [CrossRef]

29. F. Martelli, A. Sassaroli, Y. Yamada, and G. Zaccanti, “Analytical approximate solutions of the timedomain diffusion equation in layered slabs,” J. Opt. Soc. Am. A. 19, 71–80 (2002). [CrossRef]

30. A. Kienle, “Light diffusion through a turbid parallelepiped,” J. Opt. Soc. Am. A 22, 1883–1888 (2005). [CrossRef]

31. F. Martelli, A. Sassaroli, S. D. Bianco, and G. Zaccanti, “Solution of the time-dependent diffusion equation for a three-layer medium: application to study photon migration through a simplified adult head model,” Phys. Med. Biol. 52, 2827–2843 (2007). [CrossRef] [PubMed]

32. A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. I. Homogeneous case,” Opt. Express 18, 9456–9473 (2010). [CrossRef] [PubMed]

33. A. Zhang, D. Piao, C. F. Bunting, and B. W. Pogue, “Photon diffusion in a homogeneous medium bounded externally or internally by an infinetely long circular cylindrical applicator. I. Steady-state theory,” J. Opt. Soc. Am. A 27, 648–662 (2010). [CrossRef]

34. W. Cong, L. V. Wang, and G. Wang, “Formulation of photon diffusion from spherical bioluminescent sources in an infinite homogeneous medium,” Biomed. Eng. Online 3, 1–6 (2004). [CrossRef]

35. S. Takatani and M. D. Graham, “Theoretical analysis of diffuse reflectance from a two-layer tissue model,” IEEE Trans. Biomed. Eng. 26, 656–664 (1979). [CrossRef] [PubMed]

36. A. Liemert and A. Kienle, ‘Light diffusion in a turbid cylinder. II. Layered case,” Opt. Express 18, 9266–9279 (2010). [CrossRef] [PubMed]

37. I. Dayan, S. Havlin, and G. H. Weiss, “Photon migration in a two-layer turbid medium A diffusion analysis,” J. Mod. Opt. 39, 1567–1582 (1992). [CrossRef]

38. J. Sikora, A. Zacharopoulos, A. Douiri, M. Schweiger, L. Horesh, S. R. Arridge, and J. Ripoll, “Diffuse photon propagation in multilayered geometries,” Phys. Med. Biol. 51, 497–516 (2006). [CrossRef] [PubMed]

39. A. Liemert and A. Kienle, “Light diffusion in N-layered turbid media: frequency and time domains,” J. Biomed. Opt. 15, 025002 (2010). [CrossRef] [PubMed]

40. E. Okada, M. Firbank, M. Schwiger, S. R. Arridge, M. Cope, and D. T. Delpy, “Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head,” Appl. Opt. 36, 21–31 (1997). [CrossRef] [PubMed]

41. L. O. Svaasand, T. Spott, J. B. Fishkin, T. Pham, B. J. Tromberg, and M. W. Berns, “Reflectance measurements of layered media with diffuse photon-density waves: a potential tool for evaluating deep burns and subcutaneous lesions,” Phys. Med. Biol. 44, 801–813 (1999). [CrossRef] [PubMed]

42. A. Kienle and T. Glanzmann, “In vivo determination of the optical properties of muscle with time-resolved reflectance using a layered model,” Phys. Med. Biol. 44, 2689–2702 (1999). [CrossRef] [PubMed]

43. B. W. Pogue, T. O. McBride, U. L. Osterberg, and K. D. Paulsen, “Comparison of imaging geometries for diffuse optical tomography of tissue,” Opt. Express 4, 270–286 (1999). [CrossRef] [PubMed]

44. H. Dehghani and D. T. Delpy, “Near-infrared spectroscopy of the adult head: effect of scattering and absorbing obstructions in the cerebrospinal fluid layer on light distribution in the tissue,” Appl. Opt. 39, 4721–4729 (2000). [CrossRef]

45. H. Dehghani, B. A. Brooksby, P. W. Pogue, and K. D. Paulsen, “Effects of refractive index on near-infrared tomography of the breast,” Appl. Opt. 44, 1870–1878 (2005). [CrossRef] [PubMed]

46. P. K. Yalavarthy, H. Dehghani, B. W. Pogue, and K. D. Paulsen, “Critical computational aspects of near infrared circular tomographic imaging: Analysis of measurement number, mesh resolution and reconstruction basis,” Opt. Express 14, 6113–6127 (2006). [CrossRef] [PubMed]

47. M. Schweiger and S. R. Arridge, “The finite-element method for the propagation of light in scattering media: frequency domain case,” Med. Phys. 24, 895–902 (1997). [CrossRef] [PubMed]

48. J. Zhang, D. Chen, J. Liang, H. Xue, J. Lei, Q. Wang, D. Chen, M. Meng, Z. Jin, and J. Tian, “Incorporating MRI structural information into bioluminescence tomography: system, heterogeneous reconstruction and in vivo quantification,” Biomed. Opt. Express 5, 1861–1876 (2014). [CrossRef] [PubMed]

49. M. A. Anastasio, J. Zhang, D. Modgil, and P. J. L. Rivière, “Application of inverse source concepts to photoacoustic tomography,” Inverse Problems 23, 21–35 (2007). [CrossRef]

50. H. Gao, H. Zhao, W. Cong, and G. Wang, “Bioluminescence tomography in Gaussian prior,” Biomed. Opt. Express 1, 1259–1277 (2010). [CrossRef]

51. B. W. Pogue, S. Geimer, T. O. McBride, S. Jiang, U. L. Österberg, and K. D. Paulsen, “Three-dimensional simulation of near-infrared diffusion in tissue: boundary condition and geometry analysis for finite-element image reconstruction,” Appl. Opt. 40, 588–600 (2001). [CrossRef]

52. S. G. Mykhlin, Mathematical Physics: An Advanced Course (North-Holland, 1970).

53. Y. V. Egorov and M. A. Shubin, Partial Differential Equations (Springer, 1992).

54. E. Demiralp and H. Beker, “Properties of bound states of the Schrödinger equation with attractive Dirac delta potentials,” J. Phys. A 36, 7449–7459 (2003). [CrossRef]

55. M. Shendeleva, “Radiative transfer in a turbid medium with a varying refractive index: comment,” J. Opt. Soc. Am. A 21, 2464–2467 (2004). [CrossRef]

56. R. Aronson, “Boundary conditions for diffusion of light,” J. Opt. Soc. Am. A 12, 2532–2539 (1995). [CrossRef]

57. G. W. Faris, “Diffusion equation boundary conditions for the interface between turbid media: a comment,” J. Opt. Soc. Am. A 19, 519–520 (2002). [CrossRef]

58. T. J. Farrel and M. S. Patterson, “Experimental verification of the effect of refractive index mismatch on the light fluence in a turbid medium,” J. Biomed. Opt. 6, 468–473 (2001). [CrossRef]

59. J. M. Tualle, J. Prat, E. Tinet, and S. Avrillier, “Real-space Greens function calculation for the solution of the diffusion equation in stratified turbid media,” J. Opt. Soc. Am. A 17, 2046–2055 (2000). [CrossRef]

60. H. Erkol and M. B. Unlu, “Virtual source method for diffuse optical imaging,” Appl. Opt. 52, 4933–4970 (2013). [CrossRef] [PubMed]

61. T. Farrel and M. S. Patterson, “Experimental verification of the effect of refractive index mismatch on the light fluence in a turbid medium,” J. Biomed. Opt. 6, 468–473 (2001). [CrossRef]

Analytical reconstruction of the bioluminescent source with priors

Abstract

1. Introduction

2. Method

2.1. Theoretical method

2.1.1. Light diffusion in a circular homogenous tissue

2.1.2. Light diffusion in a circular heterogonous tissue

2.2. Numerical method

3. Discussion and results

Acknowledgments

References and links

Cited By

Figures (9)

Equations (39)

Optics Express