Mathematical Study and Numerical Simulation of Multispectral Bioluminescence Tomography

Multispectral bioluminescence tomography (BLT) attracts increasingly more attention in the area of optical molecular imaging. In this paper, we analyze the properties of the solutions to the regularized and discretized multispectral BLT problems. First, we show the solution existence, uniqueness, and its continuous dependence on the data. Then, we introduce stable numerical schemes and derive error estimates for numerical solutions. We report some numerical results to illustrate the performance of the numerical methods on the quality of multispectral BLT reconstruction.


INTRODUCTION
Genetically engineered mice are popular laboratory models for numerous studies directly relevant to the human healthcare. Mouse imaging allows in vivo evaluation of physiological and pathological processes under controlled conditions simulating therapeutic interventions. One of the most promising mouse imaging modalities is based on luciferase report gene [1,2]. When the corresponding substrate is administered into a mouse, the resultant biochemical reaction generates bioluminescent light at visible and infrared wavelengths. These bioluminescent photons carry important information about tumor burden, micro-metastases, infectious loci, therapeutic gene delivery, and so on. Bioluminescent imaging (BLI) utilizes sensitive photon detection techniques to observe bioluminescent signals on the body surface of the mouse.
Funded by NIH/NIBIB, we built the first bioluminescence tomography (BLT) prototype [3], and developed the BLT theory and methods [4][5][6][7][8]. Quickly, BLT has grown into a hot area, in which several groups are actively producing valuable results [9][10][11][12][13][14]. Currently, we are improving our BLT prototype with more features and better performance in hope that BLT will become a universal and powerful tool for molecular and cellular imaging in the near future. While multispectral BLT results were already reported in the literature [6,12,14], there is a critical need to establish a mathematical theory of multispectral BLT.
It is well known that all the common luciferase enzymes, from firefly (Fluc), click beetle (CBGr68, CBRed) and Renilla reniformins (HRLuc), share a wide spectral range, roughly from 400 to 750 nm [15]. Furthermore, the newly developed tricolor reporter generates bioluminescent light that is rich in green and red spectral bands. Because the optical properties of the tissues depend on the wavelength, the intensity and spectrum of bioluminescent light measured on the body surface of a mouse are a function of the spectrum of the bioluminescent source, its 3D distribution, and the individual anatomy of the mouse. Clearly, multispectral data are more informative than mixed data, and must be analyzed for the optimal BLT performance.
This paper provides a theoretical study of the multispectral BLT model. The spectrum is divided into certain numbers of bands, say i 0 bands Λ 1 , . . . , Λ i0 , with Λ i = [λ i 1 , λ i ), is a partition of the spectrum range. Denote by p the bioluminescent source distribution function. Then the bioluminescent source distribution within the band Λ i is ω i p. In this paper, we always understand the range of the index i to be 2 International Journal of Biomedical Imaging 1, . . . , i 0 . The weight ω i 0, and i0 i=1 ω i = 1. Let Ω R d (d 3) be the biological medium with the boundary Γ. Let Ω 0 be a measurable subset of Ω (Ω 0 = Ω is allowed) where the light source exists. The set Ω 0 is known as the permissible region. For each spectral band Λ i , 1 i i 0 , we use the following diffusion equation to describe the photon density (1) , μ a,i (x) and μ ¼ s,i (x) are the absorption coefficient and the reduced scattering coefficient within the band Λ i , χ Ω0 is the characteristic function of Ω 0 , that is, its value is 1 in Ω 0 , and is 0 in ΩÒΩ 0 . The bioluminescent imaging experiments are usually performed in a dark environment so that the natural boundary condition takes the form [16] Here, ∂/∂ν stands for the outward normal derivative, with γ being the refractive index of the medium. Let Γ 0 Γ with meas(Γ 0 ) > 0 be a part of the boundary where measurement takes place; Γ 0 = Γ is allowed. With the emission filter of bandpass Λ i , the measured quantity is the outgoing flux density [16]: The pointwise formulation of the BLT problem is to determine the source function from (1)-(4) for 1 i i 0 . As noted in [8], the pointwise formulation (1)-(4) for 1 i i 0 is ill-posed: (1) in general, there are infinite many solutions; (2) when the form of the source function is specified, generally there are no solutions; (3) the source function does not depend continuously on the data. We study the multispectral BLT problem through Tikhonov regularization as follows.
Denote by Q ad the admissible set of the source functions. We assume Q ad is a closed convex subset of L 2 (Ω 0 ). Examples include Q ad = L 2 (Ω 0 ), or the subset of L 2 (Ω 0 ) of nonnegatively valued functions, or certain finite dimensional subspace or subset of L 2 (Ω 0 ). For a weak formulation of the boundary value problem (1)-(2), we need the space V = H 1 (Ω). For any q ¾ L 2 (Ω 0 ), and any i = 1, . . . , i 0 , define u i = u i (q) ¾ V to be the solution of the problem From the assumptions on the data, we can apply the wellknown Lax-Milgram lemma [17,18] to conclude that the solution u i (q) exists and is unique. We note that the mapping q u i (q) is linear. Choose weights w i > 0, 1 i i 0 , with i0 i=1 w i = 1. Define a Tikhonov regularization functional [19,20] In the definition of J ε (¡), other norms than ¡ L 2 (Γ) and ¡ L 2 (Ω0) may be used if it is required to do so based on practical considerations. In this paper, we use (6) for the definition of J ε (¡) that leads to easier implementation. The following formulas for Gâteaux derivatives of J ε (¡) hold: for p, q ¾ L 2 (Ω 0 ). Thus, J ε (¡) is strictly convex when ε > 0.
We then introduce the following multispectral BLT problem.
This paper is on a study of the multispectral BLT Problem 1. In the next section, we focus on the solution existence, uniqueness, and continuous dependence on the data. In Section 3, we discuss numerical methods for the multispectral BLT reconstruction and derive error estimates for the numerical solutions. In Section 4, we include some numerical examples to show the performance of the numerical methods. We end the paper by a concluding remark summarizing the main contributions of the paper.

WELL-POSEDNESS
We first address the existence and uniqueness issue. For this purpose, we make some assumptions on the given data. We assume Ω R d is a nonempty, open, bounded set with a Lip- Weimin Han et al.

3
If Q ad L 2 (Ω 0 ) is a subspace, then the solution p ε ¾ Q ad is characterized by a variational equation Proof. For ε > 0, Problem 1 is a constrained optimization for a strictly convex objective functional over a closed convex set. Thus, it has a unique solution (see, e.g., [17,Theorem 3.3.12]), and the solution p ε is characterized by the relation J ¼ ε (p ε )(q p ε ) 0, for all q ¾ Q ad (see [17,Theorem 5.3.19]), which is the variational inequality (8). When Q ad L 2 (Ω 0 ) is a subspace, the inequality (8) reduces to (9) by a standard argument.
We then consider the continuous dependence of the solution on the data. Proof. The solution p ε depends continuously on all the data. In order to maintain the length of the proof, in the following we show the continuous dependence of p ε on ε, D i , μ a,i , Q i and ω i , 1 i i 0 . To simplify the notation, let p ε ¾ Q ad be the solution of Problem 1 with the data ε + δ ε with δ ε ε/2, in Ω, and ω i + δ ωi with δ ωi ω 0 /2. Also, in this proof, we use c for constants that may depend on ε, D 0 , ω 0 , Q i , and μ a,i for 1 i i 0 , but are independent of δ ε , δ ωi , δ μa,i , δ Di , and δ Qi for 1 i i 0 .
Similar to (8), p ε is characterized by the inequality Take q = p ε in this inequality, q = p ε in (8), and use these two inequalities to obtain After some manipulations, By the definitions of u i (q) and u i (q), we obtain the equality Thus, We can bound p ε L 2 (Ω0) by p ε L 2 (Ω0) + p ε p ε L 2 (Ω0) . Then from (12), Hence, the solution depends continuously on the data.
Next, we consider the solution behavior when ε 0+.
Note that a solution p ¾ Q ad of Problem 1 with ε = 0 is 4 International Journal of Biomedical Imaging characterized by the inequality Let S 0 Q ad be the solution set of Problem 1 with ε = 0. As in [21], the following result holds.

Proposition 1.
Assume S 0 is nonempty. Then S 0 is closed and convex. Moreover, where p 0 ¾ S 0 is the solution of Problem 1 with minimal Proof. The closedness and convexity of S 0 are easy to deduce.
Here we only prove (18). Take q = p 0 in (8), q = p ε in (17) for p = p 0 , and add the two inequalities to obtain , and p ε is uniformly bounded. Let p ε ¼ be a subsequence of p ε , converging weakly to p. Since S 0 is weakly closed, p 0 L 2 (Ω0) . Now p 0 is the unique element in S 0 with minimal L 2 (Ω 0 ) norm, p = p 0 . Thus the limit p = p 0 does not depend on the subsequence selected; consequently, p ε converges weakly to p 0 in L 2 (Ω 0 ) as ε 0. Using the relation we conclude the strong convergence p ε p 0 L 2 (Ω0) 0 as ε 0.
As a simple consequence of Proposition 1, if the solution set S 0 = p is a singleton, then p ε p in L 2 (Ω), as ε 0. Note that if Q ad is a bounded set, then S 0 is nonempty. This follows from applying a standard result on convex minimization, for example, [17,Theorem 3.3.12]. As in the first part of the proof of Theorem 1, L 2 (Ω 0 ) is a Hilbert space, Q ad L 2 (Ω 0 ) is convex and closed, J ε ε=0 : Q ad R is convex and continuous. Since Q ad is assumed to be bounded, Problem 1 with ε = 0 has a solution. Without further information on Q ad , though, we cannot ascertain uniqueness of a solution when ε = 0.

NUMERICAL APPROXIMATIONS
In this section, we discretize Problem 1 and derive error estimates for the numerical solutions. First we need to discretize the boundary value problem (5). Let T h (h: meshsize) be a regular family of finite element partitions of Ω into triangular/tetrahedral elements such that each element at the boundary Γ has at most one nonstraight face (for a three-dimensional domain) or side (for a two-dimensional domain). For each triangulation T h = K , let V h V be the linear element space. For any q ¾ L 2 (Ω 0 ), denote by The solution u h i (q) exists and is unique. Let The admissible source function space Q ad may or may not need to be discretized. In general, let Q ad,1 Q ad be nonempty, closed and convex. Later in the section, we will consider two possible choices of Q ad,1 . We then introduce the following discretization of Problem 1.
Similar to Theorem 1 and Proposition 1, we have the following result.

Proposition 2.
For ε > 0, Problem 2 has a unique solution p h ε ¾ Q ad,1 , which is characterized by the discrete variational inequality: If Q ad,1 is a subspace of L 2 (Ω 0 ), then p h ε is characterized by a variational equation: If Q ad,1 is a bounded set, then S h 0 is nonempty. In concrete situations, it is possible to show the nonemptyness of the solution set S h 0 directly. We then turn to error estimation. For this purpose, we further assume The assumptions (25) are made to ensure the validity of the solution regularity (26) used below in error estimation. Without the solution regularity property, error estimates with lower convergence orders can still be derived. Let We have the finite element interpolation error estimate: This error estimate is usually proved when Ω is a polyhedral/polygonal domain so that each element K in a finite element partition T h has straight faces/sides on its boundary (e.g., [23,24]). For applications in bioluminescence tomography, Ω is a smooth domain, and is not polyhedral. In such an application, the error estimate (27) still holds [8]. For the finite element solution of (21), there is a constant c > 0 independent of h and ε such that This error bound is shown in [8] in a somewhat different setting. Now we distinguish two cases in the discretization of the admissible set Q ad,1 . We first consider the case Q ad,1 = Q ad . This is the natural choice when Q ad is a finite dimensional subspace or subset of linear combinations of specified functions such as the characteristic functions of certain subsets of Ω. We have the following error bound. Theorem 3. With Q ad,1 = Q ad , there is a constant c > 0 independent of ε and h such that Proof. We choose q = p ε in (23), q = p h ε in (8), and use the two inequalities to obtain The error bound (29) follows from this inequality together with (28).
Further error bounds require more information on the data. We present two sample results as consequences of Theorem 3.
Assume Q ad is a bounded set in L 2 (Ω). Then p ε L 2 (Ω0) and p h ε L 2 (Ω0) are uniformly bounded with respect to ε and h. By (26), u i (p ε ) H 2 (Ω) , and hence u i (p ε ) L 2 (Γ) as well, is uniformly bounded. So from (29), we see that there is a constant c > 0 independent of ε and h such that Next, we assume the data are compatible in the sense that there exists p 1 ¾ Q ad such that u i (p 1 ) = 2AQ i on Γ 0 for 1 and so for all ε > 0, p ε L 2 (Ω0) + ε 1/2 u i (p ε ) 2AQ i L 2 (Γ0) 2 p 1 L 2 (Ω0) . Thus from (29), The first term on the right-hand side is bounded as follows: Therefore, we conclude that for some constant c > 0 independent of ε and h, International Journal of Biomedical Imaging We now consider the situation where Q ad is a general admissible set and needs to be discretized. In addition to the regular family of finite element partitions T h of Ω, let T 0,H be a regular family of finite element partitions of Ω 0 such that each element at the boundary ∂Ω 0 has at most one nonstraight face (for a threedimensional domain) or side (for a two-dimensional domain). The partitions T h and T 0,H do not need to be related; however, T h can be constructed based on T 0,H . Let Q H L 2 (Ω 0 ) be the piecewise constant space. Define Q ad,1 = Q H ad Q H Q ad . We denote the solution of Problem 2 by p h,H ε . Denote by Π H : L 2 (Ω 0 ) Q H the orthogonal projection opeartor: for q ¾ L 2 (Ω 0 ), We will use the following properties: By the elementwise formula since Q ad L 2 (Ω 0 ) is convex, we see that Π H : Q ad Q H ad , that is, for q ¾ Q ad , its piecewise constant orthogonal projection Π H q ¾ Q H ad . We first present a preparatory result.

Lemma 1. There is a constant c > 0 independent of h and H such that
Proof. From the definitions of u i (q) and u h i (Π H q), we have The sum of the last two integrals on the right-hand side can be replaced by the following, with the help of (43), )) dx, which is bounded by (40). Then after some algebraic manipulations we obtain We use the error bound inf v h ¾V h u i (q) v h H 1 (Ω) ch u i (q) H 2 (Ω) and the regularity bound (26) in (45) to obtain (42).
We now prove the following error estimate.
Weimin Han et al.

7
Proof. Similar to the first part of the proof of Theorem 3, there holds where we used the property (p h,H ε , Π H p ε p ε ) L 2 (Ω0) = 0. Then, Applying the error bound (28) and Lemma 1, we can then deduce (46). Similar to Theorem 3, we present two sample results as consequences of Theorem 4. If Q ad is bounded in L 2 (Ω 0 ), then there is a constant c > 0 independent of ε, h, and H such that If the data are compatible, then there is a constant c > 0 independent of ε, h, and H such that These error bounds involve an approximation error term p ε Π H p ε L 2 (Ω0) . This term can always be bounded by p ε L 2 (Ω0) + Π H p ε L 2 (Ω0) 2 p ε L 2 (Ω0) c. Moreover, as is shown in [8], if S 0 = , then p ε Π H p ε L 2 (Ω0) 0 as H, ε 0, and if p ε ¾ H 1 (Ω 0 ), then p ε Π H p ε L 2 (Ω0) cH p ε H 1 (Ω0) .
We underline that the above theoretical results on the numerical solutions with the second choice of Q ad,1 are still valid if Q H L 2 (Ω 0 ) is a general finite element space containing piecewise constants. The proofs of the results are the same as long as we define Π H to be the orthogonal projection operator in L 2 (Ω 0 ) onto the space of piecewise constants.
When the regularization parameter ε is chosen related to the discretization parameters h and H, we may express the error bounds in terms of the discretization parameters only. For example, from (36), max 1 i i0 u i (p ε ) u h i (p h ε ) L 2 (Γ0) ch 3/2 , and if ε = ch β , 0 < β < 3 in (36), then p ε p h ε L 2 (Ω0) Finally, we comment that when the solution set S 0 is nonempty, convergence of the numerical solution p h ε to the minimal energy solution p 0 ¾ S 0 follows from the triangle together with (18) and the convergence of p h ε to p ε in L 2 (Ω 0 ). A similar statement holds for the convergence of p h,H ε to p 0 ¾ S 0 .

NUMERICAL SIMULATION
For a numerical simulation of source reconstruction, we consider a cylindrical phantom with a diameter 20 mm and a height 20 mm. We choose the coordinate system so that the phantom is represented as The measurement boundary is the entire lateral side of the cylinder: Based on the results of the bioluminescent spectral analysis experiment in [15], we split the spectrum into three regions: Λ 1 = [400 nm, 530 nm), Λ 2 = [530 nm, 630 nm), Λ 3 = [630 nm, 750 nm], and quantify the energy distribution weights to be ω 1 = 0.29 in Λ 1 , ω 2 = 0.48 in Λ 2 , and ω 3 = 0.23 in Λ 3 . The optical parameters are assigned as follows (unit: mm 1 ): The refractive index A is 1.37 in the biological tissue. The phantom is discretized into 58161 tetrahedral elements and 10778 nodes. A total of 2170 datum nodes are distributed along Γ 0 , and simulated measurement data of photon density at datum nodes are generated from the diffusion approximation model. Two uniform light sources are embedded into the phantom. The first light source has a power 2.2 nano-Watts and is distributed in those elements whose vertices have a distance less than or equal to 2 mm from ( 4, 3, 10) T . The second light source has a power 2.4 nano-Watts and is distributed in those elements whose vertices have a distance less than or equal to 2 mm from ( 4, 3, 10) T . Denote the region where the true light sources are distributed to be D tr . Figure 1 shows the finite element mesh, the 3D and 2D views of the true source distribution in the phantom. Here and in the following figures, the 2D view is for the midsection (x 3 = 10) of the cylinder.
We choose the permissible region and use piecewise constants to approximate the source density function. In simulating the difference between discretized solution u h and solution u of the boundary value problem (1)-(2), we introduce random noises of the sizes 5%, 10%, and 20% in the datasets on the lateral surface covering phantom. Then the regularization method with ε = 1.0 ¢ 10 8 is applied to reconstruct source distribution. A modified Newton method and an active set strategy with simple constraint are used in solving the discrete Problem 2.
The reconstructed results show that with random noises of sizes 5%, 10%, and 20%, we have 80%, 77%, and 76% of the total power of the reconstructed light sources located in D tr , the region of the true light source. Figures 2-4 illustrate 3D and 2D views of the reconstructed light source distribution with 5%, 10%, and 20% random noises in surface flux density, respectively. The numerical results show that the numerical method is computationally efficient, stable, and robust with respect to noise in the measured data.

CONCLUDING REMARK
Multispectral bioluminescence tomography is a new development in optical imaging and has a great potential Weimin Han et al.  to advance the field of biomedical imaging. In this paper, we have established a mathematical framework for studies of multispectral bioluminescence tomography. We have analyzed the theoretical properties of the multispectral imaging model including the solution existence, uniqueness, and continuous dependence on the data. We have rigorously demonstrated the convergence and error bounds for the discrete source functions that are obtained by minimizing the discretized objective functions subject to the PDE constraints. Numerical examples have illustrated the performance of the numerical scheme for multispectral bioluminescence tomography. Weimin Han received his B.S. degree from Fudan University in 1983, M.S. degree from Chinese Academy of Sciences in 1986, and Ph.D. degree from University of Maryland in 1991, all in mathematics. He has been on the faculty of the Department of Mathematics, University of Iowa since 1991, as an Assistant Professor in 1991-1996, Associate Professor in 1996-1999, and Professor starting from 1999. His research interests are mathematical analysis and numerical approximations of problems arising in applications. He has authored/coauthored 7 books and more than 90 journal papers in numerical analysis, computational mechanics, and inverse problems. Then, he was Professor with Departments of Radiology, Biomedical Engineering, Civil Engineering, and Mathematics, and Director of the Center for X-Ray and Optical Tomography, University of Iowa, until 2006. Currently, he is Samuel Reynolds Pritchard Professor of Engineering, and Director of Biomedical Imaging Division, VT-WFU School of Biomedical Engineering and Sciences, Virginia Tech. His interests include computed tomography, bioluminescence tomography, and systems biomedicine. He has authored or coauthored over 300 journal articles and conference papers, including the first paper on spiral/helical cone-beam CT and the first paper on bioluminescence tomography. He is the Editorin-Chief for International Journal of Biomedical Imaging, and Associate Editor for IEEE Transactions on Medical Imaging and Medical Physics. He is an IEEE Fellow and an AIMBE Fellow. He is also recognized by a number of awards for academic achievements.