1. Introduction

Gene therapy is a breakthrough in the modern medicine, which promises to cure diseases by modifying gene expressions. A key for the development of gene therapy is to monitor the in vivo gene transfer and its efficacy in the mouse model. Traditional biopsy methods are invasive, insensitive, inaccurate, inefficient, and limited in the extent. To map the distribution of the administered gene, reporter genes such as those producing luciferase are used to generate light signals within a living mouse. Bioluminescent imaging (BLI) is an optical technique for sensing gene expression, protein functions and other biological processes in mouse models by reporters such as luciferases that generate internal biological light sources [1, 2]. The light emitted within the mouse can be captured externally using a highly sensitive CCD camera [3]. “The collection times for BLI are relatively short compared to non-optical modalities, an advantage of optical imaging methods in general” [1].

BLI has great potentials in various biomedical applications, including regenerative medicine, developmental therapeutics, treatment of residual minimal disease, and studies on cancer

stem cells. “Additionally, time-lapse whole body imaging of animals bearing xenografts of appropriately labeled cancer cells can provide new information on tumor growth dynamics and metastasis patterns that it is not possible to obtain by invasive experimental approaches” [4]. “Bioluminescence imaging (BLI) is a highly sensitive tool for visualizing tumors, neoplastic development, metastatic spread, and response to therapy” [5]. Although the spatial resolution is limited when compared with other imaging modalities, the advantages of BLI in vivo imaging are sensitivity, speed, non-invasiveness, cost, ease, and low background noise (in contrast to fluorescence, PET and other non-optical modalities) [1]. BLI could be applied to study almost all diseases in small animal models [6, 1, 7, 8, 9, 10, 2].

However, the current BLI technology has not explored the full potential of this approach. It only works in 2D imaging mode and is incapable of 3D imaging of the light source location associated with specific organs and tissues [2]. Since its first introduction in 2003 [11, 12], the bioluminescence tomography (BLT) has been undergone a rapid development [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 4, 22, 2, 23, 24]. BLT is to address the needs for 3D localization and quantification of an internal bioluminescent source distribution in a small animal. With BLT, optimal analyzes on a bioluminescent source distribution become feasible inside a living mouse. Although there are currently several approaches for this technique [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 4, 22, 2, 23, 24], the originally proposed approach in [11, 12] is still one of the promising techniques in this field. With this technique, the complete knowledge on the optical properties of anatomical components is assumed to be available from an independent tomography scan, such as a CT/micro-CT, MRI scan and/or diffuse optical tomography (DOT), by image segmentation and optical property mapping. That is, we can segment the image volume into a number of anatomical structures, and assign optical properties to each component using a database of the optical properties, or use the DOT technique for this purpose [11, 12].

Traditionally, optical tomography sends visible or near infra-red light to probe a scattering object, and reconstructs the distribution of internal optical properties, such as absorption and scattering coefficients [25, 26, 27]. In contrast to this active imaging modality, BLT reconstructs an internal bioluminescent source distribution generated by luciferase activated with reporter genes, from BLI measurements at the object surface. Mathematically, BLT is a source inversion problem based on the boundary measurement, and hence is a highly ill-posed inverse problem per se. To obtain satisfactory results for the BLT, prior knowledge must be utilized to regularize the problem. The tomographic feasibility and the solution uniqueness have been theoretically studied in [13]. It was proved that the uniqueness of the BLT does not hold in general. In our previous studies, we utilized constraints such as the non-negativity and source support [14, 15, 18, 19]. Other constraints such as the range of the source intensity may be effective as well. Another approach is to utilize spectral resolved measurement such as the hyper-spectral or multi-spectral measurement [4, 22, 24, 28] and multi-spectral source information [29] to improve the BLT reconstruction.

In previous studies, complete data measured on the full object boundary is required to reconstruct an internal bioluminescent source distribution. In practice, the measured data for BLT is often incomplete due to physical limitations as in X-ray CT [30, 31, 26]. The BLT in this situation is similar to CT image reconstruction from angle-limited data. In this work, we propose an approach for BLT from incomplete or partial boundary measurement and establish two iterative reconstruction algorithms based on the diffusion approximation. First we generalize the results on the solution uniqueness in [13] to the partial measurement case. It is proved that similar results still hold given practical optical signals on the mouse body surface, although the solution characterization is more complicated. A new important formulation is the definition of the Dirichlet-to-Neumann (or Steklov-Poincaré) map in (21) for the partial measurement

case. Then we extend two of our reconstruction methods for BLT in [14, 15, 18] to the partial measurement case. The first algorithm is a variant of the well-known expectation-maximization (EM) algorithm [32, 33, 34]. The second one is based on the projected Landweber scheme [35, 36, 37, 38]. Both methods allow the incorporation of knowledge-based constraints. Two practical constraints, the source non-negativity and support constraints, are introduced to regularize the BLT problem and produce stability. The initial choice of both methods and its influence on the regularization and stability are also discussed. Both algorithms are evaluated and validated with intensive numerical simulation and a physical phantom. Finally algorithmic issues are investigated, especially the method to avoid the inverse crime in numerical simulations.

The organization of the paper is as follows. We introduce the problem of BLT from partial measurement in §2, reformulate it in an operator form by generalizing the Dirichlet-to- Neumann map in §3, and extend our results on the solution uniqueness of BLT in §4. Then, we present our iterative BLT algorithms in §6. We report our numerical and experimental results in §7. Finally we discuss technical problems and research topics with the current BLT techniques in §8 and conclude the paper in §9.

2. Formulation of BLT

Let Ω be a bounded smooth domain in the three–dimensional Euclidean space R ³ that contains an object to be imaged. BLT is to reconstruct the source q from measurement on the boundary of Ω enclosing the source. Let u(x,θ) be the radiance in direction ϕ∈S ² at x∈Ω, where S ² is the unit sphere. A general model for light migration in a random medium is the radiative transfer equation (RTE) [39, 25, 26]

for t>0, and x∈Ω, where c denotes the particle speed, µ=µ_a+µ_s with µ_a and µ_s being the absorption and scattering coefficients respectively, the scattering kernel η is normalized such that ∫S ²η(θ·θ′)dθ′=1, and q is the internal light source. In (1), the radiance u(x,θ,t) is in Wcm^-2 sr^-1, the source term q(x,θ,t) in Wcm^-3 sr^-1, the scattering coefficient µ_s and the absorption coefficient µ_a both are given in cm^-1, and the scattering phase function η is in sr^-1 [40]. To find a solution for (1), we need the initial condition

and the boundary condition for u

where ν is the exterior normal on the boundary Γ of Ω. The forward problem (1), (2) and (3) admits a unique solution under natural assumptions on µ, µ_a and η [41]. The homogeneous condition g-(x,θ,t)=0 specifies that no photons travel in an inward direction at the boundary, except for source terms [25, p. R50], which is the case for BLT. In terms of the radiance u(x,θ,t), the measurement at a boundary points is given by

With the RTE, it is quite complex to reconstruct the light source q from the measurement g and with the above initial-boundary problem being as the forward process. This is largely due to the difficulty in computing the solution u for the forward problem (1), (2) and (3).

Typical values of the optical parameters for biological tissues are µ_a=0.1–1.0mm^-1, µ_s= 100–200mm^-1, respectively. This means that the mean free path of the particles is between 0.005 and 0.01 mm, which is very small compared to a typical object. Thus the predominant phenomenon in optical tomography is scattering rather than transport. Therefore, the diffusion approximation has been widely applied to simplify the RTE (1) in optical tomography [25, 26, 27]. The diffusion approximation assumes that the radiance u(x,θ,t) is isotropic and takes the average radiance u₀(x,t) to approximate the radiance u(x,θ,t)

$\mathrm{Let}\overline{\eta }=\frac{1}{4\pi }{\int }_{S^2}\theta ·\theta \prime \eta \left(\theta ·\theta \prime \right)d\theta \prime \mathrm{and}$

The forward process in the diffusion regime is then given by the following initial-boundary problem [39, 25, 26]

For diffusion approximation based BLT, the measured quantity at x∈Γ for t≥0 is given by

Because the internal bioluminescence distribution is relatively stable, we can use the stationary version of (9) and (11) as the forward model for BLT. By discarding all the time dependent terms in (9) and (11), the stationary forward model is given by the following boundary value problem (BVP)

In practice, it is difficult to obtain all the measurement on the boundary Γ. We consider the case in which the measurement is conducted on some disjoint patches Γ_j⊂Γ for j=1, ⋯,J, each of which is smooth, closed and connected. Let Γ_P=∪^J_j=1 Γ_j. The BLT problem from partial measurement can be stated as follows: Given the incoming light g- on Γ and outgoing radiance g on Γ_P, find a source q₀ with the corresponding diffusion approximation u₀ such that

In a typical BLT setting, g-=0, because there is no incoming light. The optical parameters D and µ_a in the above formulations can be established point-wise as mentioned above [11, 12].

3. Reformulation of BLT(P)

The uniqueness property of the BLT problem was already studied [13]. For the BLT(P) problem (15), we analyze its solution uniqueness in this section. The reader is reminded that the presentation in its mathematically accurate form will require rather technical and tedious assumptions on the domain and the coefficients. Here we will make the mathematical presentations as precise as possible while keeping a reasonable readability. For details, please refer to [13, 14] and the references therein.

In additions to the assumptions for the the BLT(P) problem (15), we assume that the parameters D>D₀>0 for some positive constant D₀ and that µ_a=0 are bounded functions. We further assume that D is sufficiently regular near Γ. For example, D is equal to a constant near Γ. Let γ₀ and γ₁ be the boundary value maps

and L be the differential operator

Given $f\in H^{\frac{1}{2}}\left({\Gamma }_P\right)$, let w ₁∈H ¹(Ω) be the solution of the following mixed boundary value problem (MBVP) [42, 43]

We define a linear operator $N_{{\Gamma }_P}$ from $H^{\frac{1}{2}}\left({\Gamma }_P\right)$ to $H^{-\frac{1}{2}}\left({\Gamma }_P\right)$(ΓP) by

On the other hand, for q₀∈L ²(Ω), we consider the following MBVP

and define another linear operator ${\Lambda }_{{\Gamma }_P}$ from L ²(Ω) to $H^{\frac{1}{2}}\left({\Gamma }_P\right)$ by

Both operators NΓ_P and ΛΓ_P are extensions of the well-known Dirichlet-to-Neumann (or Steklov-Poincareé) map ([44]) and an relevant operator Λ in the complete measurement case ([14, 15, 18]) to the partial measurement case, respectively.

Assume that q0 is a solution of the BLT(P) problem (15) with one corresponding radiance diffusion approximation solution u, given the observed g on G P and assumed g- on Γ. Then u is the unique solution of the following MBVP

On the other hand, let w₁ be defined as in (18)–(20) with f=g-+2g, w ₂ be defined as in (22)–(24), and ν=w ₁+w ₂. Then we have

By the solution uniqueness of this MBVP [42], it follows that v satisfies (26)–(28). Hence, u=ν is the required radiance u that generates the measurement on Γ_P. The measurement equation (12) implies that

Hence, q₀ satisfies the following equation

Conversely, if there exists a q₀ satisfying (33), we can construct u=v=w₁ +w₂ as above. It follows that u satisfies the forward model and the measurement equation. In summary we have

Proposition 3.1. q₀ is a solution to the (BLT(P)) problem (15) if and only if it is a solution to the equation (33).

This is an extension of the result for the BLT problem in our previous studies.[13, 14, 15, 18]

4. Solution structure of BLT(P)

In practice the BLT(P) problem is expected to have at least one solution. Therefore, we will not discuss the existence of the BLT(P) problem but focus on the uniqueness of the BLT(P) solution. The uniqueness of the BLT solution was discussed in our previous study [13]. In what follows we will demonstrate how to extend the results to the BLT(P) problem.

We need the following notations from functional analysis [45]. Let A be a linear operator from a Banach space X to a Banach space Y. The kernel or null space of A is defined as 𝒩 [A]={x∈X : A[x]=0}, and the range of A is 𝓡[A]= {y∈Y: y=A[x] for some x∈X}. For a subspace M of a Hilbert space H, M⊥ is the set of all y∈H, such that 〈y,x〉=0 for all x∈M.

By Theorem 3.1, the uniqueness of BLT(P) solution can be characterized by determining the kernel $\mathit{𝒩}\left[{\Lambda }_{{\Gamma }_P}\right]$ of the operator ${\Lambda }_{{\Gamma }_P}:L^2\left(\Omega \right)=H^{\frac{1}{2}}\left({\Gamma }_P\right)\subset L^2\left({\Gamma }_P\right)$. We need the Green’s formula [42, 43]

For $\psi \in H^{\frac{1}{2}}\left({\Gamma }_P\right)$, let

$T_{{\Gamma }_P}$ be defined by

as the unique solution in H ¹(Ω)⊂L ²(Ω) of the boundary problem

Then, we can prove with the Green’s formula that for the operators ${\Lambda }_{{\Gamma }_P}:L^2\left(\Omega \right)\to H^{\frac{1}{2}}\left({\Gamma }_P\right)\subset L^2\left({\Gamma }_P\right)$ and $T_{{\Gamma }_P}:H^{\frac{1}{2}}\left({\Gamma }_P\right)\subset L^2\left({\Gamma }_P\right)\to L^2\left(\Omega \right)$,

i.e., they are adjoint to each other,

Hence, the kernel of ${\Lambda }_{{\Gamma }_P}$ is [45]

To obtain an explicit characterization of $\mathit{𝒩}\left[{\Lambda }_{{\Gamma }_P}\right]$, let

Then we have

Proposition 4.1.

Proof. If $q\in L\left[H_{0,{\Gamma }_P}^2\left(\Omega \right)\right]$ with q=L[p] for some $p\in H_{0,{\Gamma }_P}^2\left(\Omega \right)$, then for $v=T_{{\Gamma }_P}\left[\psi \right]\in 𝓡\left[T_{{\Gamma }_P}\right]$, by the Green’s formula (34),

because ${{\gamma }_0\left[p\right]\mid }_{{\Gamma }_P}=0$, ${{\gamma }_1\left[p\right]\mid }_{{\Gamma }_P}=0$, and the boundary integral on Γ\Γ_P is equal to zero. Hence, $q\perp 𝓡\left[T_{{\Gamma }_P}\right]$. Therefore, $L\left[H_{0,{\Gamma }_P}^2\left(\Omega \right)\right]\subset 𝓡{\left[T_{{\Gamma }_P}\right]}^{\perp }$.

Conversely, assume that $q\in 𝓡{\left[T_{{\Gamma }_P}\right]}^{\perp }=\mathit{𝒩}\left[{\Lambda }_{{\Gamma }_P}\right]$. We have, by (22)–(25), there exists w ₂ such that

We have w₂∈H ²(Ω) by the regularity theory for second order elliptic partial differential equations [42, 43]. The above boundary conditions imply that $w_2\in H_{0,{\Gamma }_P}^2\left(\Omega \right)$. Hence, $q=L\left[w_2\right]\in L\left[H_{0,{\Gamma }_P}^2\left(\Omega \right)\right]$. The conclusion follows immediately.

By Proposition 3.1, the BLT(P) problem is equivalent to the linear equation (33) with q ₀ as the unknown to be found. All the solutions q ₀ to (33) form a convex set in L ²(Ω). Hence, there exists one unique solution of the minimal L ²-norm among those solutions [45]. Let this minimal norm solution be denoted as qH. Then, all the solutions can be expressed as qH+𝒩[Λ].

We summarize the above results into the following theorem.

Theorem 4.2. Assume that the BLT(P) problem is solvable. For any couple (g-,g) such that

there is one special solution qH for the BLT(P) problem (15), which is of the minimal L²-norm among all the solutions. Then, any solution can be expressed as q₀=q_H+L[p], for some $p=H_{0,{\Gamma }_P}^2\left(\Omega \right)$.

5. Uniqueness of the BLT(P) solution in RBF

It is remarkable that the uniqueness results in [13] for sources consisting of radial base functions (RBF) still hold for the BLT(P) problem after examining the assumptions and proofs therein. We assume the following conditions on Ω, D, µ_a and q₀.

C1: Ω is a bounded C ² domain of R ³ and partitioned into non-overlapping sub-domains Ω_i, i=1,2,…, I;

C2: Each Ω_i is connected with piecewise C ² boundary Γ_i;

C3: D and µ_a are C ² near the boundary of each sub-domain.

C4: D>D ₀>0 for some positive constant D ₀ is Lipschitz on each sub-domain; µ_a≥0 and µ_a∈L^p(Ω) for some p>N/2;

The first uniqueness result is for sources consisting of bioluminescent point impulses

where each a_s is a constant, and ys the location of a point source inside some Ω_i, for s=1,⋯,S.

Theorem 5.1. ([13]) Assume the conditions C1–C4 hold. If $q_0(y)={\displaystyle {\sum }_{s=1}^S{a}_s\delta (y-y_s)}$ and $q_0^{\prime }(y)={\displaystyle {\sum }_{s=1}^{S\prime }{a}_s^{\prime }\delta (y-y_s^{\prime })}$ are two solutions to the BLT(P) problem (15) with the same measurement on ΓP, then S=S′ and there is a permutation τ of {1,⋯,S} such that a′s=aτ(s) and y′s=yτ_(s).

The second uniqueness result is for sources as a linear combination of RBFs, which includes the previous type of sources in (45) as the limiting case,

where each $g_s\in L^2\left(B_{r_0^s,r_1^s}\left(x_s\right)\right)$ is continuous, the source centers {xs} are distinct, and each source support $B_{r_0^s,r_1^s}\left(x_s\right)\subset \subset {\Omega }_i$ for some 1≤i≤I. For each 0≤r₀<r₁<∞, x ₀∈R³, $B_{r_0,r_1}\left(x_0\right)$ denotes a hollow ball specified by r ₀<‖x-x ₀‖ < r₁ for r ₀>0 or a solid ball specified by ‖x-x₀‖ < r₁ for r₀=0. Let φ(r) be defined as follows. For µ_a=0,

and for µ_a>0,

For the second result we need to replace the condition C4 with the following condition C4* and a new condition C5:

C4*: D and µ_a are piecewise constant in the sense that there exist constants D ₁,⋯,D _I>0 and µ ₁,⋯,µ ₁=0 such that D(x)≡D_i and µ_a(x)=µ_i,∀x∈ Ω_i.

C5: For each sub-domain Ω_m, 1≤m≤ I, there exists a sequence of indices 1≤i ₁,i ₂,…, i _k≤I with the following connectivity property: the intersection ${\Gamma }_P\cap {\Gamma }_{i_1}$ contains a smooth C ² open patch and ${\Gamma }_{i_j}\cap {\Gamma }_{i_{j+1\mathrm{}}}$ contains a smooth C ² open patch, for j=1,…,k-1, and ${\Omega }_{i_k}={\Omega }_m$;

Note that the condition C4* is a special case of the condition C4. The second uniqueness result can be stated as follows.

Theorem 5.2. ([13]) Assume the conditions C1–C3, C4* and C5 hold. If $q_1\left(y\right)=\sum _{s=1}^S{g}_s\left(\parallel y-y_s\parallel \right)\chi B_{r_0^s,r_1^s}\left(y_s\right)$ and $q_2\left(y\right)=\sum _{s=1}^{S\prime }{G}_s\left(\parallel y-Y_s\parallel \right)\chi B_{R_0^s,R_1^s}\left(X_s\right)$ are two solutions to the BLT(P) problem (15) with the same measurement on ΓP, then S=S′ and there exist a permutation τ of {1, ⋯ ,S} and a map C: {1,⋯,S} → [1,I] such that Y_s = y_τ(s) ∈ Ω_C(s) and

where φ_j is given by (48) or (47) with D=D_j and µ_a=µ_j.

6. Reconstruction methods

In our previous studies, the BLT problem with complete measurement is reformulated as a linear equation by the Dirichlet-to-Neumann map and solved with the proposed EM-like and constrained Landweber (CL) iterative algorithms therein. In this section we follow the same approach to establish reconstruction methods for the BLT(P) problem by solving the linear equation (33).

Let $b=N_{{\Gamma }_P}\left[g^{-}+2g\right]+g\in H^{\frac{1}{2}}\left({\Gamma }_P\right)$. By (33), the BLT solution q₀ satisfies the following equation

The extended EM and CL iterative algorithms for solving the BLT(P) problem (50) are established in the following.

6.1. Method I: EM method

Based on the formulation in [26, 14], let

which is the log likelihood function when the measured data b is subject to Poissonian distribution. By the maximum principle of elliptic partial differential equations [42], it follows that ${{\rm A}}_{{\Gamma }_P}\left[q_0\right]\ge 0$ and b≥0 if q₀≥0, q₀≠0, g≥0, g≠0 and g-≥0. We try to find a solution for the BLT(P) problem by performing the following optimization

We first assume that q₀>0 is a minimizer of F. The case of q₀≥0 can be handled similarly as the limiting case. We need to find the Frechet derivative of F. Let

where v is an arbitrary bounded function of L ²(Ω), and compute

${\frac{d}{\mathrm{dt}}f\left(t\right)\mid }_{t=0}={\int }_{{\Gamma }_P}\left\{b\frac{1}{{\Lambda }_{{\Gamma }_P}\left[q_0\right]}-1\right\}{\Lambda }_{{\Gamma }_P}\left[v\right]d\Gamma ={\int }_{\Omega }{\Lambda }_{{\Gamma }_P}^{*}\left[\frac{b}{{\Lambda }_{{\Gamma }_P}\left[q_0\right]}-1\right]v\mathrm{dx}.$

Hence, the Frechet derivative of F is

If q₀>0 is a solution of (52), it follows that F′[q₀]=0. The general case of q₀=0 is given by the following Kuhn-Tucker condition [26]

Let ${\varphi }_1={\Lambda }_{{\Gamma }_P}^{*}\left[1\right]=T_{{\Gamma }_P}\left[1\right]$, i.e., the solution of the following MBVP, by (40) and (36)–(38),

It follows from the maximum principle of elliptic partial differential equations that 0 <ϕ ₁=1 [46]. Because ${\Lambda }_{{\Gamma }_P}^{*}=T_{{\Gamma }_P}$ by (40), the Kuhn-Tucker condition (55) can be rewritten as

Then we obtain the following EM formula for BLT from partial measurement

6.2. Method II: CL method

Assume that we have some prior knowledge about the source represented as a convex set

which is a closed convex subset of L ²(Ω). Let P_𝓒 be the orthogonal projection operator from L ²(Ω) to 𝓒. Then the CL scheme, or projected Landweber scheme, for solving (50), is given as follows

where λ_n is a relaxation parameter and ${\Lambda }_{{\Gamma }_P}^{*}=T_{{\Gamma }_P}$ by (40). The convergence property of the CL scheme was studied in [47, 48] and improved in [37, 38]. Conditions for the relaxation parameter depend on the operator norm $\parallel {\Lambda }_{{\Gamma }_P}^{*}{\Lambda }_{{\Gamma }_P}\parallel ={\parallel {\Lambda }_{{\Gamma }_P}\parallel }^2$ For example, $0<{\lambda }_n<\frac{2}{{\mid {\Lambda }_{{\Gamma }_P}\mid }^2}$ [37, 38]. The limit of q⁽ⁿ⁾0 is a solution to the constrained least-squares problem

In this work, the non-negativity and source support constraints are applied. The non-negativity constraint implies that the source is non-negative and is utilized by setting the negative parts of iterated sources to be zero. The source support constraint assumes that the source is non-zero only within some sub-region Ω0 of the region Ω, and is utilized by setting the iterated sources to be zero outside Ω₀. The choice of Ω₀ is discussed in §6.3.2. Both constraints are applied after each iteration.

6.3. Relevant issues

6.3.1. Computational environment

A common feature of the proposed EM and CL methods is that they are of an iterative nature. At each step, they require the evaluation of the operators ${\Lambda }_{{\Gamma }_P}$ by (25) based on the MBVP (22)–(24) and $T_{{\Gamma }_P}$ by (35) based on the MBVP (36)–(38). Both MBVPs are of the same type can be solved with the finite element method (FEM) [49]. In this work, both MBVPs are solved with the FEM software Comsol_™ and Matlab^™. Note b in (50) for both methods can be computed by (21) after solving the MBVP (18)™(20). The computer is a Dell^™ workstation, Precision 670, with dual Intel^™ Xeon CPUs of main frequency 2.8GHz and 6GB memory. The operating system is Microsoft^™ Windows XP Professional X64 Edition. Other details are reported in §7.

6.3.2. Choice of q⁽⁰⁾₀

In all iterative image restoration methods, we have to initialize the process. We propose the following method to choose an initial guess q⁽⁰⁾₀. By Green’s formula, let v=u and p=1 in (34), we have,

If we replace q₀ with q⁽⁰⁾₀, we obtain

Because u=w₁ by the the maximum principle of elliptic partial differential equations [42], we have

In practice, the source q₀ is compactly supported on a subset Ω₀ inside Ω. We choose q₀ such that it is equal to a positive constant in its support Ω₀ and zero otherwise. Hence

where Q₀ is a constant, and χΩ₀ (x)=1 on Ω0 and is zero otherwise. The support Ω₀ of the source q₀ is part of the prior knowledge, which was termed the permissible region in and could be inferred from the measured data [19]. Please refer to §8 for further discussions. Because only the data g on ΓP is available, the final estimate is

where |Ω₀| is the volume of Ω₀.

6.3.3. Convergence criteria

The convergence criteria for both algorithms may include (1) when the iteration number n reaches an assumed maximum number; (2) when the successive incremental |q₀⁽ⁿ⁺¹⁾-q⁽ⁿ⁾ ₀| is smaller than an assumed error level. In this work, the convergence criterion is by manually setting the iteration numbers to fixed numbers for both methods, respectively.

7. Experimental results

7.1. Numerical experiments

For inverse problems, numerical tests of reconstruction methods usually make use of simulated data from the numerical solution of the forward problem. One typical issue is coined as the inverse crimes in [50]. This happens especially when insufficient rough discretization or the same discretization are used for the forward and inverse process, because “it is possible that the essential ill-posedness of the inverse problem may not be evident” [50, p. 304]. Hence, the results could be overly optimistic and unreliable. “Unfortunately, not all of the numerical reconstructions which have appeared in the literature meet with this obvious requirement” [50, p. 133]. As suggested, “it is crucial that the synthetic data be obtained by a forward solver which has no connection to the inverse solver under consideration” [50, p. 133].

One approach to avoid the inverse crime is to use different discretizations in the forward and the inverse process [50, p. 304]. Because our formulation of the reconstruction methods is analytical and independent of any discretization, we can use different finite elements (shape functions and meshes) in Comsol^™ for the forward process and reconstruction algorithms, respectively. Moreover, we can change the mesh sizes with the adaptive mesh technique at each iteration step for the reconstruction algorithms to solve the MBVPs (22)–(24) and (36)–(38). During the iteration intermediate results at different meshes are interpolated to the required nodes with the built-in bilinear interpolation method in Comsol^™, when values at the nodes are required.

We have carries out intensive numerical simulation for various numerical phantoms. Various algorithmic settings such as the initial choices, iteration numbers, shape functions and meshes of the FEM have also been evaluated. Because the algorithms depend on multiple factors, the optimal setting is still under investigation. Nevertheless both methods can reliably reconstruct the sources in most settings. Due to the limited space, we report one representative result with the CL method in the following.

In this experiment, a numerical phantom with the same geometrical structure as the physical phantom in §7.2 is used. This is a cylindrical heterogeneous mouse chest phantom (Fig. 1(a)) of 30mm height and 30mm diameter. Its structure is shown in Fig. 1(b). Three sources of the form

for i=1,2,3, are set up in this simulation, where Ω_i is a ball centered at x_i: Ω_i={x : ‖x-_xi‖<r_i}. The radius r_i are all set to 1mm. The sources are centered at x ₁=(-0.9cm, 0.25cm, 0cm), x ₂=(-0.9cm, -0.25cm, 0cm) and x ₃=(0.9cm, 0.25cm, 0cm), with the intensity values A ₁=25.1µW/cm³, A ₂=23.3µW/cm³ and A ₃=25.1µW/cm³, respectively. These intensity values are set according to the total source power of the physical phantom in §7.2 so that the total source power of each source is equal to one of the physical sources in §7.2. The optical coefficients of the phantom are set as in Table 1.

Table 1. Optical parameters for the phantom

View Table | View all tables in this article

As discussed above, we use different finite elements (shape functions and meshes) for the forward process and inverse process, respectively. The information of the finite elements is

presented in Table 2, which is part of the output of the command meshinfo of Comsol^™.

 

Fig. 1. (a) A heterogeneous mouse phantom consisting of bone (B), heart (H), lungs (L), and muscle (M). (b) A cross-section through two luminescent sources in the left lung and another source in the right lung. The four arrows show the directions of the CCD camera for data measurement.

Download Full Size | PDF

Table 2. Finite element information for the simulation

View Table | View all tables in this article

Figure 2 shows the results obtained with the CL method. The initial support or the permissible region is set to Q₀={(x,y,z) : 0.8<(x ²+y ²)^1/2 <1.2, -0.15<z<0.15}. The relaxation coefficient λ is manually set to λ=20. The computational overhead for the cases of complete measurement and partial measurement is about the same. It takes about 6 hours for 70 iterations to get the results in Fig. 2. The figures are generated with the commands postplot, geomplot and meshplot of Comsol^™ with manually adjusted parameters at different views, respectively.

In the case of the partial measurement from the front view, the source in the right lung close to the back view is not reconstructed, see Fig. 2 (c) and (d). This is reasonable because that source is far from the measurement surface in terms of the mean-free path. When complete measured data is used for reconstruction, this source could be reliably reconstructed and is shown in Fig. 2 (a) and (b). Quantitative results about the location accuracy are compiled into Table 3. The absolute error is defined by $\sqrt{{\left(x_{i,r}-x_i\right)}^2+{\left(y_{i,r}-y_i\right)}^2+{\left(z_{i,r}-z_i\right)}^2}$. (x_i,r,y_i,r,z_i,r) is the reconstructed center of each source and is estimated interactively from the reconstructed

source distribution from different views. The relative error is defined by the absolute error divided by $\sqrt{x_i^2+y_i^2+z_i^2}$. The results demonstrate that the same location accuracy for the left two sources can be obtained with only the partial measurement in the front view.

 

Fig. 2. Reconstructed results by the CL method and a cross-section at z=0cm. (a) and (b) are results from data measured at the four views. (c) and (d) are from data measured at the front view only.

Download Full Size | PDF

Table 3. Quantitative results for the reconstructed locations of the three sources at S1=(-0.90, 0.25, 0), S2=(-0.90, -0.25, 0) and S3=(0.90, 0.25, 0), respectively. The unit is cm.

View Table | View all tables in this article

Another quantitative index for BLT is the reconstructed source power compared to its orig-inal value. As reported in the literature, the source power was estimated as the source integral ∫q(x)dx of the source intensity over its support [19, 51]. Another approach for the source power estimation for a RBF source is based on the result of Theorem 5.2 by computing the source moment ∫φ_s(‖x-X_s‖)G_s(‖x-X_s‖)dx over its support, which is equal to $4\pi {\int }_{R_0^s}^{R_1^s}{r}^2{\phi }_s\left(r\right)G_s\left(r\right)\mathrm{dr}$. The reconstructed source intensities and sizes are estimated interactively from orthogonal views of the reconstructed source distributions. Table 4 presents the results of the reconstructed source integrals and its absolute and relative errors with the true value, respectively, for each source. Table 5 shows the corresponding results for the source moments. The differences of the results in both tables are discussed in §8.

Table 4. Quantitative results for the reconstructed source integrals of the sources. The sources are listed in the order as in Table 3. Their true values are 105.1, 97.4 and 105.1, respectively. The unit is nW.

View Table | View all tables in this article

Table 5. Quantitative results for the reconstructed source moments of the sources. The sources are listed in the order as in Table 3. Their true values are 125.5, 116.5 and 125.5, respectively. The unit is nW.

View Table | View all tables in this article

7.2. Physical experiment

We use the same physical phantom in our previous study [19]. As described in §7.1, a cylindrical heterogeneous mouse chest phantom (Fig. 1(a)) of 30mm height and 30mm diameter was designed and fabricated. It consisted of four different materials high-density polyethylene (8624K16), nylon 6/6 (8538K23), delrin (8579K21) and polypropylene (8658K11) (McMaster- Carr supply company, Chicago, IL, US) to represent muscle (M), lungs (L), heart (H) and bone (B), respectively, as shown in Fig. 1(b). A luminescent light stick (Glowproducts, Canada) was selected as the testing source. The stick consisted of a glass vial containing one chemical solution and a larger plastic vial containing another solution with the former being embedded in

the latter. By bending the plastic vial, the glass vial can be broken to mix the two solutions after which luminescent light was emitted. The particular dye in the chemical solution was for red light, and it could last for approximately 4 hours at an emission wavelength range between 650nm and 700nm, being close to that of the red spectral region of the luciferase. Two small holes of diameter 0.6mm and height 3mm were drilled in the phantom with their centers at (-0.9cm, 0.15cm, 0.0cm) and (-0.9cm, -0.15cm, 0.0cm) in the left lung region of the phantom, respectively. Two red luminescent liquid filled catheter tubes of 1.9mm height and 0.56mm diameter were placed inside the two small holes, respectively. We measured the total power of the red luminescent liquid filled polythene tubes with the CCD camera. They were 105.1 nW and 97.4 nW, respectively.

 

Fig. 3. (a) A cross-section through two hollow cylinders for hosting luminescent sources in the left lung. The four arrows show the direction of the CCD camera during data acquisition. (b) The measurements at the four views combined along the phantom side surface with unit µW/cm²

Download Full Size | PDF

In our bioluminescent imaging, a CCD camera (Princeton Instruments VA 1300B, Roper Scientific, Trenton, NJ) was used for data acquisition on the phantom surface. The collected bioluminescent views were transformed from grey-scale pixel values into equivalent numbers in physical units. We used the same calibration approach described in the previous study [19]. The measured data on the CCD were transformed to the form of the measurement equation

(12) through a geometric optics mapping of light beams. The required optical parameters of the phantom in Table 1 were the same as in the previous study [19].

Due to the limited space, representative results from the physical phantom are given in Fig. 4(a) and Fig. 4(b) using the EM algorithm from only the data measured in the four views of the side surface of the phantom. We conducted experiments for reconstruction from the data measured in one or several of the four views. Figure 4(c) and 4(d) are the results reconstructed using the EM algorithm from only the data measured in the front view. The results from the measured data in other single views or combinations of these three views were not encouraging, because the signal-to-noise ratios were too low. For the results in Fig. 4 (a)–(d), the iteration number for the EM algorithm was set to 50 along with an initial source support region, Q₀={(x,y,z) : 0.8<(x²+y²)1/2<1.2,-0.15<z<0.15}. It takes about 5 hours for the 50 iterations of the EM algorithm. For the CL algorithm, similar results were obtained but the separation of the two sources was not as good as that obtained using the EM algorithm.

 

Fig. 4. Representative results reconstructed by the EM algorithm. (a) and (c) are the sources reconstructed by the EMalgorithm from the data measured in the four views and in the front view only, respectively. (b) and (d) are cross-sections at z=0 cm of the sources in (a) and (c), respectively.

Download Full Size | PDF

Quantitative results about the location accuracy and source power computed by the source integral are summarized in Table 6 and Table 7. The source powers estimated by source moments are not available because the original source is not of RBF sources, please refer to discussions in §8. The reconstructed source intensities and sizes are estimated interactively from orthogonal views of the reconstructed source distributions. The information of the finite element in this

experiment is presented in Table 8.

Table 6. Quantitative results for the reconstructed locations of the two sources at S1=(-0.90,0.15,0) and S2=(-0.90,-0.15,0), respectively. The unit is cm.

View Table | View all tables in this article

Table 7. Quantitative results for the reconstructed source integrals of the sources. The sources are listed in the order as in Table 6. Their true values are 105.1 and 97.4, respectively. The unit is nW.

View Table | View all tables in this article

Table 8. Finite element information for the physical phantom experiment

View Table | View all tables in this article

8. Discussions

Because there is no unique solution to the BLT(P) problem in the general case by Theorem 4.2, one may consider to utilize the minimal norm solution qH as the solution of the BLT(P) problem. The minimal norm solution qH is unique and also called the minimal energy solution, and advocated in other applications [38, 52]. However, it can be similarly demonstrated as in our previous study [14] that the minimal energy source solution is not favorable for the BLT(P) problem in general. Because sources of compact supports are commonly encountered in practice, it can be proved in the same way that such sources cannot be found as the minimal norm solution for the BLT(P) problem [14]. It has been reported that adequate prior knowledge must be utilized to obtain a physically favorable BLT solution [11, 12, 13, 14, 15, 16, 18, 19, 20]. From the physical perspective, the multi-spectral technique is also a promising approach to resolve this issue, although it needs extra hardware and exposure time at a compromised signal-to-noise ratio [4, 22, 24, 28, 29]. The methods for BLT from partial measured data in this work can be combined with the multi-spectral technique as well.

The iterative approach provides a mechanism for incorporating prior knowledge based constraints and has been widely used in practice. In this paper we have established two iterative algorithms for BLT from partial measurement. In this work, two constraints, the non-negativity and source support constraints, are applied. The EM algorithm induces the non-negativity of the source because of the maximum principle of elliptic partial differential equations, if the initial choice is non-negative. The source support with the EM algorithm is automatically implied because the support will not increase during the iteration due to its multiplication operation in (60). Hence, it seems that the EM algorithm is more preferable than the CL algorithm with those two constraints. However, the CL algorithm is more flexible to add other constraints and has the established convergence property. Other constraints highlighted by Theorem 5.1 and 5.2 are to restrict the solution space to a sub-space of bioluminescent source distributions so that the solution uniqueness holds to a practical extent. Nevertheless, the source support constraint helps to resolve the non-uniqueness because of the property of both the EM and CL algorithms.

Another advantage of the proposed iterative methods for BLT is its regularization property. For inverse problems, popular approaches include the Tikhonov [53] and the iterative regularization methods [54]. For the latter the regularization is achieved by setting the iteration number. Both of our proposed methods are iterative and impose regularization to the BLT problem by setting iteration numbers. It was proved in [55] that the CL and Tikhonov regularization methods are equivalent. As discussed in the previous paragraph, the source support of the initial choice also helps to regularize the problem and produce stability.

There are remaining mathematical issues with the proposed algorithms. For the EM algorithm, its convergence has not been established, though it converges in our experiments. For the CL method, there is no guide in choosing the relaxation coefficient λ_n. This parameter depends on the operator norm $\parallel {\Lambda }_{{\Gamma }_P}^{*}{\Lambda }_{{\Gamma }_P}\parallel ={\parallel {\Lambda }_{{\Gamma }_P}\parallel }^2$, which is equivalent to find the minimal eigenvalue of ${\Lambda }_{{\Gamma }_P}^{*}{\Lambda }_{{\Gamma }_P}$ or the minimal singular value of ${\Lambda }_{{\Gamma }_P}$. This can be reduced to a boundary eigenvalue problem of partial differential equations. Both problems are left for future investigation.

There are also remaining physical issues with the proposed approach. We have used a geometric optics mapping of light beams to transform the measured data by a CCD detector to the surface measurement equation (12). More advanced techniques should be used to improve this process, such as the recently proposed non-contact measurement technique for fluorescence tomography [56]. Furthermore, the diffusion approximation can be improved by the radiative transfer equation [51].

There are two methods in quantifying the source power, namely, by the source integral and source moment. By (49), the source integral is equal to the source moment when the absorption coefficient µ_a=0. As presented in Table 4 and 5, the error in source power computed by the source moment is smaller than that by the source integral in numerical simulations. This is in accordance with the theoretical result in Theorem 5.2. The results of source moments presented in this work are based on estimated source sizes interactively from orthogonal views of the reconstructed source distributions. More sophisticated techniques of RBF approximations should be applied to get more accurate estimations [57]. Although estimating source powers by the source moment is more accurate than by the source integral, it is unclear how the source moment is related to the total biological activities and how to extend it to a general source without the RBF structure. Both methods of computing source powers will be investigated in real biological experiments in our future work.

The computational overhead reported in this work is about 5–6 hours. Most of the computation is spent on solving MBVPs, for which there are sophisticated parallel algorithms. We are working to improve the computing performance with the Linux cluster technology, the most popular parallel computing technology at present.

9. Conclusions

To overcome the limitations in the previous BLT studies, we have formulated a mathematical framework for BLT from partial boundary measurement, extended our results on the solution uniqueness, and proposed the two iterative reconstruction algorithms based on the diffusion approximation. Either of the methods is suitable for incorporating practical constraints. Two practical constraints, the source non-negativity and support constraints, have been introduced to regularize the BLT problem and produce stability. The initial choice of both methods and its influence on the regularization and stability have also been discussed. Numerical and physical phantoms have been used to evaluate and validate the proposed algorithms with quantitative results. Technical problems and research topics with the BLT technique have also been discussed.