Algorithms for bioluminescence tomography incorporating anatomical information and reconstruction of tissue optical properties

Mohamed A. Naser; Michael S. Patterson

doi:10.1364/BOE.1.000512

1. Introduction

Optical molecular imaging in small animals has been attracting much interest in recent years as a noninvasive method of visualization of disease progression and response to treatment [1–4]. Small animals can be labeled with optical molecular probes which emit light that can be detected at the surface giving information about the distributions of specific genes or proteins associated with these sources and so realizing a noninvasive investigation of the small animals [5]. The direct image captured by the CCD camera does not reflect the actual internal source distribution as in the visible and near infrared range the light transported inside the tissue suffers from multiple scattering and absorption processes before escaping. The goal of bioluminescence tomography (BLT) is to resolve the molecular process inside the tissue by measuring the light over the animal surface and reconstructing the correct distribution of the bioluminescence sources.

For wavelengths where scattering dominates absorption, the diffusion model has been extensively used to describe light propagation inside tissue [6–16]. Diffuse optical tomography (DOT) involves reconstructing the spatial distribution of the optical properties of tissue from the measured light fluence rate at the tissue surface. The measured data are used to estimate the scattering and absorption coefficients of different tissue types inside the animal at a specific wavelength. In this study, the reconstructed image for the optical properties for wavelengths in the range 580-620 nm is then used with the measured light at the boundary to solve the BLT problem and reconstruct the correct distribution of the bioluminescence sources.

The reconstruction of the tissue optical properties and BLT are ill-posed problems [17]. The points at the surface boundary are the only accessible data points and their number is much less than the number of variables representing the absorption and scattering coefficients at all points inside the tissue. This precludes a unique solution to the problem [11]. Therefore, the solution should include some priori information about the distribution of the optical properties inside the tissue and the characteristics of the bioluminescence sources. Using the anatomical information of the x-ray CT image, different regions inside the animal can be distinguished due to their differences in x-ray attenuation. By making use of that, and by assuming that each identified region has the same optical properties, the number of variables will be significantly reduced and a unique solution can be obtained [18]. For the BLT problem, different researchers have used sparsity regularization [19–21], to describe the bioluminescence sources as localized and sparse. Making use of these source characteristics, an iterative minimization algorithm can be used to get the source distribution using an updated permissible region. The permissible region is initially the entire domain and is then reduced iteratively to get the best localization of the reconstructed sources.

In this paper, a novel approach to solving the BLT problem is described. The algorithm uses the anatomical information of the CT image to guide the cw DOT reconstruction which is then used for solving the BLT problem. The model is applied to 2D cylindrical objects to estimate the optical properties of different regions inside the object assuming that every element of each region has the same optical properties. The light fluence rate at the boundary is written as a Taylor series expansion around an initial guess corresponding to a homogenous object with the same optical properties in all tissues. The expansion is approximated by the first three terms in the series which include the first and the second order derivatives and the higher orders are neglected. The first and second order derivatives represent the Jacobian and Hessian and they give the change of light fluence rate at the boundary due to small change in the tissue optical properties. The derivatives are obtained by the direct. The Jacobian and Hessian are used in an iterative algorithm to reconstruct the tissue optical properties. The reconstructed optical properties are close to the actual values even in the presence of noise. The advantage of the proposed iterative DOT algorithm in this paper over the previous Levenberg-Marquardt (LM) approach is that there is no need of any regularization parameter as the second order derivative is already given and the symmetric matrix formed from the Jacobian is not ill-posed due to the small number of variables representing the tissue optical properties. In addition, the algorithm does not require the initial guess of the tissue optical properties to be very close to the actual values as in the case of linear approximation without considering the second order derivatives. In the presence of noise, the inclusion of the second order derivatives reduces the error in the reconstructed optical properties compared to those recovered using only the first order approximation.

The second part of the model developed in the paper is used to solve the BLT for the 2D object based on the reconstructed optical properties and measured light data. The Green’s function of the system at specific wavelengths is calculated numerically using the reconstructed optical properties. A minimization problem is solved where the objective function is the first norm of the difference between the actual measured light fluence rate and the values calculated using the Green’s function. The solution of the minimization problem gives an estimate of the source strength and distribution that results in the best fit to the measured data. The objective function is minimized iteratively such that in every iteration, the permissible region in which the source is contained is reduced, removing the points corresponding to small values and making the search region smaller, hence reducing the number of variables in the ill-posed problem. The main advantage of the proposed algorithm over previous algorithms that used the permissible region technique is that we do not need to know from the beginning where the source should be in order to fix the region where we search. Instead we start with the whole object as the initial permissible region which is then shrunk iteratively until the BL sources are properly reconstructed. The objective function uses the first norm in addition to a zero norm penalty term to localize the source distribution. The choice of the first norm for the objective function instead of the second norm has some advantages, as it shows better stability in the presence of noise and it is faster to calculate. The reason for the preference of the L1 norm over the L2 norm for noisy data is that the L1 norm is robust to outliers and large amplitude anomalies, see for example [22]. The time required to solve the minimization problem is greatly reduced if the first and second derivatives of the objective function are provided. Using the first norm will make the second derivatives vanish and avoid the huge memory required to store the calculated second derivative matrix. The two models of optical properties reconstruction and BLT rely on solving the diffusion equation numerically using the method of finite elements (FE). The numerical technique described in the paper to solve the diffusion equation is described for the 2D problem, but the upgrade to solve 3D problems will not change the algorithms. It will, however, change the basis function of the finite elements and the formulation that converts the differential equation to a matrix equation.

2. DOT and BLT forward model

In the forward model, the light flux from the tissue surface must be calculated assuming that we know the tissue optical properties and source distribution. Light propagation in a dense elastic-scattering medium like tissue is similar to neutral-particle transport which can be described by the radiative transport equation (RTE) - a specific form of the generalized Boltzmann equation [23,24]. In the visible and near infrared, some approximations can be applied to the RTE to avoid its high computation cost and detailed knowledge of tissue optical properties, and the simple diffusion equation can be used to describe the light propagation. The diffusion equation has been extensively used to model light propagation in tissues where the scattering of light is dominant over absorption. The diffusion equation at a specific wavelength in the frequency domain and in 2D space is given by [23].

[- (\frac{\partial}{\partial x} κ (x, y; λ) \frac{\partial}{\partial x} + \frac{\partial}{\partial y} κ (x, y; λ) \frac{\partial}{\partial y}) + (μ_{a} (x, y; λ) + \frac{i ω}{c})] φ (x, y; λ; ω) = q (x, y; λ; ω),

where

φ (x, y; λ; ω)

is the light fluence rate at position

(x, y)

and wavelength λ and the light modulation frequency is given by

ω = 2 π f

, and the speed of light in tissue is given by c. The isotropic source distribution is given by

q (x, y; λ; ω)

. The spatial distribution of the tissue optical properties at wavelength λ are given by the absorption coefficient

μ_{a} (x, y; λ)

and the diffusion coefficient

κ (x, y; λ)

, where the diffusion coefficient is given by

κ (x, y; λ) = \frac{1}{3 (μ_{a} (x, y; λ) + {μ^{'}}_{s} (x, y; λ))},

where

{μ^{'}}_{s} (x, y; λ)

is the reduced scattering coefficient which is defined as

{μ^{'}}_{s} = μ_{s} (1 - g)

, where

μ_{s}

is the scattering coefficient and g is the anisotropy factor. To solve the differential equation in (1), the boundary condition at the air-tissue boundary should be defined. The boundary condition considered is a Robin-type boundary condition [7], in which the light fluence rate at the boundary is equal to the gradient of fluence rate at the boundary multiplied by a factor to account for internal reflection due to refractive index mismatch. At the boundary, the light fluence rate satisfies the equation

φ (x, y; λ; ω) + 2 κ (x, y; λ) A \hat{n} (x, y) . \nabla φ (x, y; λ; ω) = 0,

where

\hat{n} (x, y)

is a unit vector pointed outwardly normal to the surface, and A is derived from Fresnel’s law as [25]

A = \frac{(2 / (1 - R_{0})) - 1 + {| \cos θ_{c} |}^{3}}{1 - {| \cos θ_{c} |}^{2}},

where the critical angle

θ_{c} = \sin^{- 1} (1 / n)

and

R_{0} = {(n - 1)}^{2} / {(n + 1)}^{2}

and n is the tissue refractive index at the boundary while the air refractive index is assumed 1. Equation (1) can be solved numerically using the method of finite elements (FE), where the space is discretized using small triangular elements. Assuming that the space is divided using a mesh which contains M triangles and N nodes, the differential equation can be transformed to a matrix equation [26]

{[K]}_{N \times N} φ_{N \times 1} = {[q]}_{N \times 1},

where the light fluence rate can then be obtained from Eq. (5) by a simple matrix inversion.

3. Optical properties Reconstruction

The forward model described in the previous section can be used to calculate the light fluence rate at the boundary due to any source distribution inside the tissue assuming known optical properties. There are actually two inverse problems that have been solved in this paper; the first is the DOT problem which is to reconstruct unknown optical properties assuming known sources and the second one is the BLT problem which is to reconstruct unknown sources assuming known optical properties. Therefore, estimating the unknown source distribution is dependent on first reconstructing the unknown optical properties of in-vivo tissues by solving the first inverse problem. Since the first inverse problem is to estimate the optical properties of many points of the mesh inside the tissue from few data points at the boundary, it is an ill-posed problem and the solution is not unique [11]. One way to overcome this is to reduce the number of variables by constraining the nodes corresponding to the same region to have the same optical properties making use of the information provided by CT imaging. In this case the number of variables will be equal to twice the number of regions, as each region will have two parameters: the scattering and absorption coefficients.

Assuming that the light fluence rate data captured by the detectors at a specific wavelength around the object can be written as a Taylor series expansion such that

φ_{d} = φ_{d}^{0} + \sum_{i} \frac{\partial φ_{d}}{\partial u_{i}} δ μ_{i} + \frac{1}{2} \sum_{i, j} \frac{\partial^{2} φ_{d}}{\partial μ_{i} \partial μ_{j}} δ μ_{i} δ μ_{j} + \dots,

where

φ_{d}

is the light fluence rate value at a detector d and d varies from 1 to the total number of detectors and

φ_{d}^{0}

is the value corresponding to an initial estimate of the optical properties and known sources.

δ μ_{i}

gives the change in the optical properties for region i and i = 1, 2, …., 2N and N is the total number of regions identified from the x-ray CT image. The vector

δ μ

is given by

δ μ = [\begin{matrix} δ μ_{s 1} \\ ⋮ \\ δ μ_{s N} \\ δ μ_{a 1} \\ ⋮ \\ δ μ_{s N} \end{matrix}],

and

μ_{s i}

and

μ_{a i}

are the scattering and absorption coefficient of region i. The Jacobian (

J_{d}

) or the sensitivity row vector

\frac{\partial φ_{d}}{\partial μ_{i}}, i = 1, 2, \dots, 2 N

gives the rate of change of the light fluence rate due to small changes in the tissue optical properties. The Hessian (

H_{d}

)

\frac{\partial^{2} φ_{d}}{\partial μ_{i} \partial μ_{j}}, i, j = 1, 2, \dots, 2 N

gives the second derivative of the light fluence rate with respect to small changes of the optical properties. By considering only terms up to the second order, Eq. (6) can be written in the compact form

φ_{d} = φ_{d}^{0} + J_{d} δ μ + \frac{1}{2} δ μ^{T} H_{d} δ μ

Multiplying both sides of (8) by $J_{d}^{T}$ , and with some rearrangement we get

δ μ = {(J^{T} J)}^{- 1} (J^{T} (δ φ - \frac{1}{2} δ μ^{T} H_{d} δ μ)),

where

δ φ_{d} = φ_{d} - φ_{d}^{0}

. Equation (9) is a nonlinear equation which gives the variation of the optical properties in terms of the variation in the light fluence rate at the boundary. For

N_{d}

detectors, Eq. (9) is given by

δ μ = {(J^{T} J)}^{- 1} (J^{T} (δ φ - \frac{1}{2} \sum_{d = 1}^{N_{d}} δ μ^{T} H_{d} δ μ)),

where

J = [\begin{matrix} \frac{\partial φ_{1}}{\partial μ_{s 1}} & \dots & \frac{\partial φ_{1}}{\partial μ_{s N}} & \frac{\partial φ_{1}}{\partial μ_{a 1}} & \dots & \frac{\partial φ_{1}}{\partial μ_{a N}} \\ ⋮ & \dots & ⋮ & ⋮ & \dots & ⋮ \\ \frac{\partial φ_{N_{d}}}{\partial μ_{s 1}} & \dots & \frac{\partial φ_{N_{d}}}{\partial μ_{s N}} & \frac{\partial φ_{N_{d}}}{\partial μ_{a 1}} & \dots & \frac{\partial φ_{N_{d}}}{\partial μ_{a N}} \end{matrix}],

and for every detector d, where d varies from 1 to N_d

H_{d} = [\begin{matrix} \frac{\partial^{2} φ_{d}}{\partial μ_{s 1}^{2}} & \dots & \frac{\partial^{2} φ_{d}}{\partial μ_{s 1} \partial μ_{s N}} & \frac{\partial^{2} φ_{d}}{\partial μ_{s 1} \partial μ_{a 1}} & \dots & \frac{\partial^{2} φ_{d}}{\partial μ_{s 1} \partial μ_{a N}} \\ ⋮ & \dots & ⋮ & ⋮ & \dots & ⋮ \\ \frac{\partial^{2} φ_{d}}{\partial μ_{s N} \partial μ_{s 1}} & \dots & \frac{\partial^{2} φ_{d}}{\partial μ_{s N}^{2}} & \frac{\partial^{2} φ_{d}}{\partial μ_{s N} \partial μ_{a 1}} & \dots & \frac{\partial^{2} φ_{d}}{\partial μ_{s N} \partial μ_{a N}} \\ \frac{\partial^{2} φ_{d}}{\partial μ_{a 1} \partial μ_{s 1}} & \dots & \frac{\partial^{2} φ_{d}}{\partial μ_{a 1} \partial μ_{s N}} & \frac{\partial^{2} φ_{d}}{\partial μ_{a 1}^{2}} & \dots & \frac{\partial^{2} φ_{d}}{\partial μ_{a 1} \partial μ_{a N}} \\ ⋮ & \dots & ⋮ & ⋮ & \dots & ⋮ \\ \frac{\partial^{2} φ_{d}}{\partial μ_{a N} \partial μ_{s 1}} & \dots & \frac{\partial^{2} φ_{d}}{\partial μ_{a N} \partial μ_{s N}} & \frac{\partial^{2} φ_{d}}{\partial μ_{a N} \partial μ_{a 1}} & \dots & \frac{\partial^{2} φ_{d}}{\partial μ_{a N}^{2}} \end{matrix}]

The reconstruction algorithm is summarized in the flowchart in Fig. 1 .

Fig. 1 Flowchart for iterative solution of the optical properties reconstruction.

Download Full Size | PDF

The derivatives of the light fluence rate at the detectors corresponding to the scattering and absorptions coefficients in the Jacobian and Hessian matrices can be calculated using the direct method. For N regions within the mesh, a change or perturbation for either $μ_{s}$ or $μ_{a}$ or both of them depending on whether we calculate the first derivative or the second derivative is made in each region by a small fraction such as 1% and these data are used to calculate the corresponding light fluence rate at detectors to get the change. This process is repeated for every source. For example, the first and second derivatives of detector number k due to changes in the optical properties of regions i and j are given by

{\frac{\partial φ_{k}}{\partial μ_{i}} |}_{μ_{i}} = \frac{φ_{k} (μ_{i} + Δ μ_{i}) - φ_{k} (μ_{i} - Δ μ_{i})}{2 Δ μ_{i}}, μ_{i} = μ_{s i} or μ_{a i},

and

\begin{array}{l} {\frac{\partial^{2} φ_{k}}{\partial μ_{i}^{2}} |}_{μ_{i}} = \frac{φ_{k} (μ_{i} + Δ μ_{i}) + φ_{k} (μ_{i} - Δ μ_{i}) - 2 φ_{k} (μ_{i})}{Δ μ_{i}^{2}} \\ {\frac{\partial^{2} φ_{k}}{\partial μ_{i} \partial μ_{j}} |}_{μ_{i}, μ_{j}} = \frac{1}{4 Δ μ_{i} Δ μ_{j}} (\begin{array}{l} φ_{k} (μ_{i} + Δ μ_{i}, μ_{j} + Δ μ_{j}) + φ_{k} (μ_{i} - Δ μ_{i}, μ_{j} - Δ μ_{j}) \\ - φ_{k} (μ_{i} + Δ μ_{i}, μ_{j} - Δ μ_{j}) - φ_{k} (μ_{i} - Δ μ_{i}, μ_{j} + Δ μ_{j}) \end{array}) \end{array}

Due to the large range in the values of the light fluence rate at different detectors, it is better and more accurate to use the logarithm of the absolute value instead of the magnitude of the light fluence rate [14]. So in equations from 14 to 20, $φ_{k} \to \log | φ_{k} |$ , and since the modulation frequency is set to zero, the phase of the light fluence rate at all detectors is equal and all information is only in the magnitude. The reason for using zero modulation frequency is that the same camera used for the bioluminescence measurement is employed for optical properties reconstruction. The camera gives information about magnitude but not the phase of the light fluence rate.

The 2D object (Fig. 2 ) used in the simulations is 25 mm in diameter and represents an approximate transverse view of the abdomen of a mouse. It is a simplified version of Fig. 6(a) in Ref [27]. It consists of 6 different regions representing different organs and tissues in the mouse. The regions from 1 to 6 represent the skin (E1), the bowel and gut (E2), the kidneys (E3, E4), the bone (E5) and adipose tissue (E6). The different regions are assumed to be distinguished through CT x-ray imaging due to the differences in their attenuation coefficients and all elements within each region are assumed to have the same optical properties. In Fig. 2, the source-detectors setup used in the simulation is shown. There are 16 sources uniformly distributed around the object and for every source, the light fluence rate is calculated at 48 detectors on the opposite side facing the source and constituting a 270° field of view as shown in Fig. 2. A full 360° is not practical to achieve with a camera and external sources. The total number of detector readings $N_{d} = 16 * 48 = 768$ , so the size of the Jacobian is 768 × 12 and the size of the Hessian matrix $H_{d}$ for every one of the 768 detector readings is 12 × 12.

Fig. 2 Schematic of the object used in the simulation illustrating 6 different regions: E1 is the skin, E2 is the bowel, E3 and E4 are the kidneys, E5 is the bone and E6 is adipose tissue. The source-detectors setup is also shown in the figure; 16 sources (one is indicated by the star) uniformly distributed are used and associated with every source 48 detectors forming an approximately 270° of view; the total number of data readings are 16x48 which is 768.

Download Full Size | PDF

Fig. 6 BLT reconstruction of one bioluminescence Gaussian source localized in the left kidney; 3% Gaussian noise is added. (a) The actual source with a total power of 0.4845 nW. (b) The reconstructed source with a total power 0.2397 nW. (c) The reconstructed source with a total power 0.4975 nW when the actual source power used is 0.969 nW.

Download Full Size | PDF

4. Optical properties calculations

The absorption and scattering coefficient for each organ and tissue in the heterogeneous object are estimated using empirical formulas describing the spectral variation in each tissue. The spectral variation of the scattering coefficient can be described by the empirical formula [28]

{μ^{'}}_{s} (λ) = a \times λ^{- b},

where a and b are fitting parameters describing the spectral variation in each tissue. The absorption coefficient is calculated assuming that light absorption in tissues depend only on the percentage concentration of oxy-hemoglobin (HbO₂), deoxy-hemoglobin (Hb), and water (W) in the tissue and is given by [29]

μ_{a} (λ) = S_{B} (x μ_{aHb} (λ) + (1 - x) μ_{{aHbO}_{2}} (λ)) + S_{W} μ_{aW} (λ),

where

x = {HbO}_{2} / ({HbO}_{2} +Hb)

and

S_{B}

and

S_{W}

are fitting parameters to match the absorption data available in the literature and as shown in [30] they were consistent with the blood and water fractions in different mouse organs. The scattering and absorption coefficient for different tissues in the heterogeneous object are calculated at three different wavelengths using Eqs. (15) and (16) and they are shown in Table 1 . The values for the parameters a, b, x,

S_{B}

and

S_{W}

in Eqs. (15) and (16) are reported in [28] and the absorption coefficients for the hemoglobin and water used are given in [29].

Table 1. Optical properties for each region in the heterogeneous object calculated at three different wavelengths. The values of the scattering and absorption coefficients are given in (mm⁻¹)

View Table

5. Bioluminescence tomography

The Green’s function of the system gives the response at any point due to a unit excitation at any other point. For the diffusion equation, the Green’s function in 2D space is defined as

[- (\frac{\partial}{\partial x} κ (x, y; λ) \frac{\partial}{\partial x} + \frac{\partial}{\partial y} κ (x, y; λ) \frac{\partial}{\partial y}) + μ_{a} (x, y; λ)] G (x, x^{'}; y, y^{'}; λ) = δ (x - x^{'}) δ (y - y^{'})

The Green’s function defined in (17) gives the light fluence rate at the point $(x, y)$ due to a point source excitation at the point $(x^{'}, y^{'})$ at zero modulation frequency. For any other excitation of the form $s (x, y; λ)$ , the light fluence rate can be obtained from the Green’s function and is given by

φ (x, y; λ) = \iint_{x^{'}, y^{'}} G (x, x^{'}; y, y^{'}; λ) s (x^{'}, y^{'}; λ) d x^{'} d y^{'}

Equation (17) can be obtained in a matrix form by discretizing the space using, for example, the method of finite elements and the differential equation in Eq. (17) can be transformed to

D G = I \Rightarrow G = D^{- 1},

where D is the discretized version of the diffusion operator shown in Eq. (17) and I is the identity matrix. The light fluence rate due to any source distribution is given by the simple matrix multiplication of the Green’s function and the source vector.

The bioluminescence tomography problem is to obtain the estimate of the source strength and distribution that gives the best agreement with the light fluence rate captured by the detectors around the object. So the bioluminescence sources should satisfy the equation

G_{p} s = φ,

where

G_{p} = G (p, :)

is a matrix that gives the Green’s function values at the boundary points corresponding to delta function excitations at any other point. The index p corresponds to the index of the detectors points at the boundary and if the detector points do not coincide with the edge points of the FE mesh, interpolations are used to get the values of the Green’s function at the detector positions. The direct solution to the problem by inverting the Green’s function matrix such that

s = {(G_{p}^{T} G_{p})}^{- 1} G_{p}^{T} φ

does not give a correct solution. This is because the square matrix

G_{p}^{T} G_{p}

is an ill-posed matrix because the number of data points obtained from the detectors around the object surface is generally much smaller than the number of variables represented by the values of the sources at each point inside the object. In order to solve this problem, we can make use of information regarding the source characteristics. It is not likely to find the source very close to the tissue surface (the diffusion approximation itself may not be suitable to describe the photon transport in this case and other models should be considered), so we can exclude the region within about one transport length of the surface from the permissible region. Thus the initial guess for the permissible region is the whole object area except a small shell of one transport length close to edge. We can also make use of the fact that, for many applications, the sources are sparse in space, so we are looking for a few localized sources with a total number of source points less than or equal to the number of measured light fluence rate data points. In this case, an iterative solution can update the permissible region to shrink the domain by considering the most likely points where the sources can be found and excluding points of small, zero or negative values from the permissible domain. Every iteration includes solving a minimization problem which has the form

\begin{array}{l} min \sum_{λ} {‖ G_{p} (λ, R) s (R) - φ (λ) ‖}^{1} \\ s . t . 0 \leq s \leq s_{\max} \\ s . t . R \in Permissible Region \end{array}

The minimization is done to the superposition at different wavelengths of the first norm of the difference between the multiplication of the Green’s function and the estimated sources (which gives the calculated light fluence rate at the boundary points) and the measured light fluence rate. It is assumed that the source strengths at different wavelengths are the same, but this is not necessary if the source emission spectrum is known. In this case, the source in Eq. (22) will be $s (λ, R) \to s (R) e (λ)$ , where $e (λ)$ is the source emission spectrum. After every iteration, the permissible region is reduced by dividing the number of points in the permissible region by a number greater than one to get the number of points in the new permissible region by removing points corresponding to the lowest source values. The choice of this number is optional depending on the choice of the number of iterations and the initial and final number of points required in the permissible region for the final solution. The simplest way to do this is to choose the points corresponding to the largest source values to be the new permissible region. A better method that ensures no important points are lost is to choose the points such that the point is selected when the source values at the point and the surrounding neighbor points are high. This is achieved by enclosing every point in the permissible region with a circle of radius equal to one transport length to account for the source values of the neighbor points as well and calculating the total power of the source inside the circle. This will give priority to points in the center of the circular domain over points at the edge of the domain which is desirable when we have point or Gaussian sources to reconstruct. The minimization of the objective function is performed with the requirement that the source values lie between zero and a maximum estimated value. The algorithm terminates when the number of points in the permissible region becomes smaller than the number of detectors.

One of the practical issues when solving the minimization problem in Eq. (22) is to add a penalty term to the objective function to direct the minimizer to choose localized source distributions rather than distributed ones. For example, a point source at a specific location inside the tissue and sources forming a ring centered at the same point can give the same light fluence rate at the boundary points and are hence degenerate. In order to remove this degeneracy and form localized sources, the objective function can be modified to have the form

min \sum_{λ} {‖ G_{p} (λ) s - φ (λ) ‖}^{1} + α {‖ s ‖}^{0},

where α is a negative number. This penalty is chosen to be the zero norm of the source term and it gives the largest positive source values that minimize the difference between the calculated and measured values of the boundary light fluence rate. This will lead to the choice of strong sources deeper inside the tissue instead of weak sources near the boundary. When the number of points in the permissible region become equal to or smaller than the number of detectors, α is set to zero to give the correct solution in the last few iterations.

6. Results and discussions

6.1 Optical properties reconstruction

The optical properties reconstruction algorithm described in the flowchart in Fig. 1 is applied to the heterogeneous object shown in Fig. 2. All elements within each region in the object are assumed to have the same optical properties. We are assuming that we start with an accurate segmented image provided by x-ray CT. It could also be provided by MR or by a deformable atlas of mouse anatomy. The question of whether CT will work in this regard remains to be determined in a prototype instrument. The contrast obtainable with x-ray CT is highly dependent on the radiation dose that is acceptable in a given application. In the computer simulation, 16 sources and 48 detectors associated with each source have been used. The external light source is simulated as a point isotropic source located one transport length, $1 / {μ^{'}}_{s}$ , from the boundary. The forward model algorithm is run to generate the light fluence rate at 48 detectors covering 270° field of view on the opposite side of the source as shown in Fig. 2. The forward model subroutine is called 16 times for the 16 sources and the total number of detector readings is 16 × 48 = 768. Since the modulation frequency is zero, the light fluence rate will be a real number.

The Jacobian and Hessian can be calculated by performing a small perturbation to the optical properties around an initial guess. The initial values for the scattering and absorption coefficients can be taken to be the same in all regions (homogeneous structure) and a small perturbation in one of the optical properties is made to calculate the change in the light fluence rate due to this change. The initial value for the scattering and absorption coefficients for all regions are 2 and 0.02 mm⁻¹, respectively. For every region, both the scattering and absorption coefficients are perturbed by ±1% to calculate the partial derivatives of the light fluence rate with respect to the scattering and absorption coefficients of that region as shown in Eq. (13), and therefore the forward model subroutine is called $16 (sources) \times 4 (\pm Δ μ_{s, a}) \times 6 (regions) = 384$ times to calculate the Jacobian. As an example,

J_{11} = {\frac{\partial \log (φ_{1})}{\partial {μ^{'}}_{s 1}} |}_{\begin{array}{l} {μ^{'}}_{s} = 2 \\ μ_{a} = 0.02 \end{array}} = [\begin{array}{l} \log (φ_{1} (2.01, 2, \dots, 2; 0.02, \dots, 0.02)) \\ - \log (φ_{1} (1.99, 2, \dots, 2; 0.02, \dots, 0.02)) \end{array}] / (2 \times 0.01)

To calculate the Hessian, the change in the light fluence rate due to perturbation of two optical properties whether in the same region or different regions should be obtained as described in Eq. (14). For every detector reading, the assigned Hessian matrix is a 12 × 12 matrix. The diagonal terms call the forward model 16(sources) × [4( ± Δµ_s,a) × 6(regions) + 1(unperturbed)] = 400 times, and the off diagonal terms require 16(sources) × [4( ± Δµ_s,a) × (11 + 10 + … + 1)] = 4224 times. Only half of the matrix needs to be calculated due to the symmetry. Making use of calculating the Jacobian, the diagonal terms require only 16 calls to the forward model for the unperturbed term.

Figure 4 shows a comparison between the actual and reconstructed scattering and absorption coefficients and the corresponding relative errors in all regions. The actual optical properties were calculated using Eqs. (15) and (16) and the results are shown in Table 1. These optical properties are used to calculate the light fluence rate which will represent the measured data. Gaussian noise with standard deviation of 1% of the signal is added to these data to simulate an actual experiment and provide the reference data $φ_{m}$ for the reconstruction algorithm. The result in Fig. 4 shows that the reconstructed optical properties match very well with relative error <1.5% for a noise power of 1% using the first order approximation.

Fig. 4 Comparisons between the actual and reconstructed optical properties of different tissue regions, see Fig. 2, at 600 nm in the heterogeneous object and percentage relative error. Gaussian noise with standard deviation of 1% is added and the reconstruction model considers the 1st order approximation (the Jacobian). (a) and (b) show the actual (blue line with circles) and reconstructed (red line with stars) scattering and absorption coefficients in different regions, and (c) and (d) show the relative errors in the reconstructed scattering and absorption coefficients.

Download Full Size | PDF

Figure 5 shows a comparison between the results of the reconstruction algorithm when using the first order and second order approximation at higher level of noise power. Using the second order approximation rather than the first order one can reduce the maximum reconstruction relative error from 11% to < 7% at noise power of 3% as in Fig. 5.

Fig. 5 (a) and (b) show comparisons between the actual (blue lines with circles) and reconstructed (red lines with stars) scattering and absorption coefficients of different tissue regions at 600 nm in the heterogeneous object using the first order approximation (the Jacobian); the Gaussian noise added is 3%. (c) and (d) Show the same results using the second order approximation (the Hessian).

Download Full Size | PDF

6.2 BLT reconstruction

The reconstructed optical properties in the presence of 3% noise at the wavelength 580 nm, 600 nm, and 620 nm were used for the BLT reconstruction as this is a realistic test of the complete process. The actual optical properties shown in Table 1 were used in the forward model to simulate the measured data due to the BLT sources inside the object. Gaussian noise of 3% standard deviation was added to these data to simulate experimental measurements. 150 detectors were uniformly distributed around the cylindrical object. This is the same as the number of edge points of the mesh used in the simulations. The light fluence rate at the detector is calculated by interpolation of the values at the neighboring edge points. The Green’s functions at different wavelengths are calculated from Eq. (19) using the reconstructed optical properties. The source powers at the three wavelengths are assumed to be the same. An iterative solution to the minimization problem as described in the flowchart in Fig. 3 is used to get the estimate of the source values that gives the best fit to the measured light fluence rate at the detectors. The objective function is minimized using the fmincon function of Matlab. The source is constrained between 0 and 1, where it is assumed that its maximum value does not exceed 1 nw/mm². The first and second derivatives of the objective function are provided to enhance the minimization speed. The fist derivative used in the model for the objective function is given by $grad= \sum_{λ} C^{T} (λ) * G_{p} (λ) + α$ , where $C (λ) = sign (G_{p} (λ) s - φ (λ))$ is a column vector of size equal to the number of detectors and its elements values equal 1 for $G_{p} (λ) s > φ (λ)$ , −1 for $G_{p} (λ) s < φ (λ)$ , and 0 for $G_{p} (λ) s = φ (λ)$ where the first norm is not differentiable. The second order derivative used is zero. At every iteration, both the source vector and the Green’s function are changed such that $s \leftarrow s (R)$ and $G_{p} (λ) \leftarrow G (index_p, R; λ)$ , where index_p is the index corresponding to the detector points and R is the index of points corresponding to the permissible region. The permissible region is reduced in every iteration by removing points with lower power and keeping the points that have high power.

Fig. 3 Flowchart for iterative solution of the bioluminescence sources reconstruction algorithm.

Download Full Size | PDF

Figure 6 shows the result of the BLT reconstruction algorithm for a single Gaussian source with total power of 0.4845 nW and FWHM of 1 mm located in the center of the left kidney. Figure 6(b) shows the reconstruction after 30 iterations where there are 15 points in the permissible region. The total simulation time is 89 seconds corresponding to average time per iteration of 2.967 seconds. The simulation is slower at the beginning because of the large number of points in the permissible region and accelerates as the region shrinks. The permissible region here starts at 1389 points and the algorithm is terminated when the number of points becomes 15 after 30 iterations. The number of points in the permissible region is reduced by a factor of (1389/15) ^ (1/30) = 1.1629 at each iteration. This factor can be altered by the user and affects the number of iterations required for solution. The choice here of the final number of points in the permissible region, which should be less than the number of detectors, is optional. However, the reconstruction results are improved when the number of points of the reconstructed sources becomes closer to the actual sources. The total power of the reconstructed source is 0.2397 nW, and it becomes 0.4975 nW if the actual source power is doubled to 0.969 nW. The total power of the reconstructed source is calculated by integrating the source density over a circle that encloses the non-zero values of the source points. Although the total power of the reconstructed source is only about half its actual value, the maximum values of the actual and reconstructed sources are comparable. The maximum values of the actual and reconstructed sources for a total power of 0.4845 of the actual source are 0.0739 and 0.0637 nW/mm², respectively and for a total power of 0.969 of the actual source are 0.1478 and 0.1152, respectively. The reason for the differences in the reconstructed total power is that the minimizer algorithm source constraint leads to having few localized sources in the permissible region; however the localized source does not have a specific shape in space. Therefore, if the localized source values are distributed uniformly, the edge points that are close to the detectors will have higher contributions to the light fluence rate and so reduce the total source power required to fit the boundary data. Although the reconstructed source power is only half the actual source power, the algorithm conserves the linearity of the source power reconstruction. The results in Fig. 6 show the linearity of the reconstruction as the power of the reconstructed source is doubled when the actual source power is doubled.

Figure 7 shows the results of the BLT reconstruction algorithm if we have a source in each kidney. The reconstruction can detect the position of the sources with good accuracy and also the relative strength of the two sources.

Fig. 7 BLT reconstruction of two bioluminescence Gaussian sources localized in the left and right kidneys with added 3% Gaussian noise. The actual sources are shown in (a), (c), and (e), while the corresponding reconstructed sources are shown in (b), (d), and (f) respectively. The two sources have the same strength in (a) and (b) and the total powers of the actual left and right sources are 0.96907 and 0.97743 nW and the left and right reconstructed sources are 0.260 nW and 0.22 nW, respectively. (c) and (d) Show the results when the left source power is the double the right source power; and the total powers of the actual left and right sources are 0.96907 and 0.4887 nW and the left and right reconstructed sources are 0.236 nW and 0.117 nW, respectively. (e) and (f) Show the case when the right source power is double the left one; The total powers of the actual left and right sources are 0.4845 and 0.9774 and the left and right reconstructed sources are 0.1431 nW and 0.2329 nW, respectively.

Download Full Size | PDF

7. Conclusion

Two numerical algorithms based on finite element methods for optical properties reconstructions and BLT in an object with known structure have been developed. A second order equation for the light fluence rate dependence on the scattering and absorption coefficients for all regions inside the object has been formulated and used to calculate the Jacobian and Hessian employed in the iterative solution of the reconstruction. For low noise, the first order approximation is sufficient to give good results with small relative error. For higher noise power, the second order term shows better performance and faster reconstruction. The BLT algorithm developed in the paper is based on calculating the Green’s function of the system at different wavelengths. An adaptive modification of the permissible region is performed to favor the recovery of a few strong sources. The results show good agreement between the localization of the actual and reconstructed sources in addition to the ability to estimate relative source magnitudes.

Acknowledgments

This work was supported in part by grants from the Ontario Ministry of Research and Innovation postdoctoral fellowship program.

References and links

1. V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol. 23(3), 313–320 (2005). [CrossRef] [PubMed]

2. R. Weissleder and U. Mahmood, “Molecular imaging,” Radiology 219(2), 316–333 (2001). [PubMed]

3. J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discov. 7(7), 591–607 (2008). [CrossRef] [PubMed]

4. J. Tian, J. Bai, X. P. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, “Multimodality molecular imaging,” IEEE Eng. Med. Biol. Mag. 27(5), 48–57 (2008). [CrossRef] [PubMed]

5. C. H. Contag and M. H. Bachmann, “Advances in in vivo bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. 4(1), 235–260 (2002). [CrossRef] [PubMed]

6. S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. 20(2), 299–309 (1993). [CrossRef] [PubMed]

7. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995). [CrossRef] [PubMed]

8. H. Jiang, “Frequency-domain fluorescent diffusion tomography: a finite-element-based algorithm and simulations,” Appl. Opt. 37(22), 5337–5343 (1998). [CrossRef] [PubMed]

9. A. H. Hielscher, A. D. Klose, and K. M. Hanson, “Gradient-based iterative image reconstruction scheme for time-resolved optical tomography,” IEEE Trans. Med. Imaging 18(3), 262–271 (1999). [CrossRef] [PubMed]

10. X. Gu, Q. Zhang, L. Larcom, and H. Jiang, “Three-dimensional bioluminescence tomography with model-based reconstruction,” Opt. Express 12(17), 3996–4000 (2004). [CrossRef] [PubMed]

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

12. W. Cong, G. Wang, D. Kumar, Y. Liu, M. Jiang, L. V. Wang, E. A. Hoffman, G. McLennan, P. B. McCray, J. Zabner, and A. Cong, “Practical reconstruction method for bioluminescence tomography,” Opt. Express 13(18), 6756–6771 (2005). [CrossRef] [PubMed]

13. N. V. Slavine, M. A. Lewis, E. Richer, and P. P. Antich, “Iterative reconstruction method for light emitting sources based on the diffusion equation,” Med. Phys. 33(1), 61–68 (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(3), 365–367 (2006). [CrossRef] [PubMed]

15. D. C. Comsa, T. J. Farrell, and M. S. Patterson, “Quantification of bioluminescence images of point source objects using diffusion theory models,” Phys. Med. Biol. 51(15), 3733–3746 (2006). [CrossRef] [PubMed]

16. C. Kuo, O. Coquoz, T. L. Troy, H. Xu, and B. W. Rice, “Three-dimensional reconstruction of in vivo bioluminescent sources based on multispectral imaging,” J. Biomed. Opt. 12(2), 024007 (2007). [CrossRef] [PubMed]

17. X. Song, D. Wang, N. Chen, J. Bai, and H. Wang, “Reconstruction for free-space fluorescence tomography using a novel hybrid adaptive finite element algorithm,” Opt. Express 15(26), 18300–18317 (2007). [CrossRef] [PubMed]

18. P. K. Yalavarthy, B. W. Pogue, H. Dehghani, and K. D. Paulsen, “Weight-matrix structured regularization provides optimal generalized least-squares estimate in diffuse optical tomography,” Med. Phys. 34(6), 2085–2098 (2007). [CrossRef] [PubMed]

19. Y. Lu, X. Zhang, A. Douraghy, D. Stout, J. Tian, T. F. Chan, and A. F. Chatziioannou, “Source reconstruction for spectrally-resolved bioluminescence tomography with sparse a priori information,” Opt. Express 17(10), 8062–8080 (2009). [CrossRef] [PubMed]

20. P. Mohajerani, A. A. Eftekhar, J. Huang, and A. Adibi, “Optimal sparse solution for fluorescent diffuse optical tomography: theory and phantom experimental results,” Appl. Opt. 46(10), 1679–1685 (2007). [CrossRef] [PubMed]

21. N. Cao, A. Nehorai, and M. Jacobs, “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Opt. Express 15(21), 13695–13708 (2007). [CrossRef] [PubMed]

22. J. F. Claerbout and F. Muir, “Robust modeling with erratic data,” Geophysics 38(5), 826–844 (1973). [CrossRef]

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

24. A. D. Klose, “Transport-theory-based stochastic image reconstruction of bioluminescence sources,” J. Opt. Soc. Am. A 24(6), 1601–1608 (2007). [CrossRef]

25. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng. 25(6), 711–732 (2009). [CrossRef] [PubMed]

26. L. J. Volakis, A. Chatterjee, and L. C. Kempel, Finite Element Method for Electromagnetics, (IEEE Press, New York, 1998).

27. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Effect of optical property estimation accuracy on tomographic bioluminescence imaging: simulation of a combined optical-PET (OPET) system,” Phys. Med. Biol. 51(8), 2045–2053 (2006). [CrossRef] [PubMed]

28. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. 50(17), 4225–4241 (2005). [CrossRef] [PubMed]

29. S. A. Prahl, http://omlc.ogi.edu/spectra/index.html (Oregon Medical Laser Clinic) (2001).

30. R. P. Brown, M. D. Delp, S. L. Lindstedt, L. R. Rhomberg, and R. P. Beliles, “Physiological parameter values for physiologically based pharmacokinetic models,” Toxicol. Ind. Health 13(4), 407–484 (1997). [PubMed]

	580 nm		600 nm		620 nm
	${μ^{'}}_{s}$	$μ_{a}$	${μ^{'}}_{s}$	$μ_{a}$	${μ^{'}}_{s}$	$μ_{a}$
Skin	2.6006	0.1299	2.5054	0.0391	2.4167	0.0178
Bowel	1.3741	0.0202	1.3175	0.0073	1.2650	0.0041
Kidney	2.8013	0.1216	2.6615	0.0372	2.5329	0.0176
Bone	3.0847	0.1042	2.9347	0.0328	2.7966	0.0146
Adipose	1.3037	0.0077	1.2805	0.0031	1.2584	0.0022

Algorithms for bioluminescence tomography incorporating anatomical information and reconstruction of tissue optical properties

Abstract

1. Introduction

2. DOT and BLT forward model

3. Optical properties Reconstruction

4. Optical properties calculations

5. Bioluminescence tomography

6. Results and discussions

6.1 Optical properties reconstruction

6.2 BLT reconstruction

7. Conclusion

Acknowledgments

References and links

Cited By

Figures (7)

Tables (1)

Equations (24)

Biomedical Optics Express