## Abstract

Parameterizing the bioluminescent source globally in Gaussians provides several advantages over voxel representation in bioluminescence tomography. It is mathematically unique to recover Gaussians [Med. Phys. **31**(8), 2289 (2004)] and practically sufficient to approximate various shapes by Gaussians in diffusive medium. The computational burden is significantly reduced since much fewer unknowns are required. Besides, there are physiological evidences that the source can be modeled by Gaussians. The simulations show that the proposed model and algorithm significantly improves accuracy and stability in the presence of Gaussian or non- Gaussian sources, noisy data or the optical background mismatch. It is also validated through *in vivo* experimental data.

© 2010 OSA

## 1. Introduction

Bioluminescence tomography (BLT) [1–3] is an emerging molecular imaging modality, which can be used to monitor physiological and pathological activities at molecular levels, specially visualize primary tumor growth and tumor cell metastasis. In BLT [4–13], one attempts to recover both location and intensity of bioluminescent sources from boundary measurements.

Conventionally the bioluminescent source is represented in voxels. For example, in a triangulation {*τ _{i}, i≤N*} of the domain

*Ω*, the source in piecewise constants is

*1*is the basis function supported on the subdomain

_{i}*τ*that is 1 on

_{i}*τ*, and 0 otherwise, and

_{i}*q*is the average intensity value on

_{i}*τ*to be recovered.

_{i}The accurate source reconstruction in voxels is highly challenging, which intrinsically comes from its non-uniqueness and severe illposedness [3]. For example, even with a large number of boundary data, which is practically much less than the number of voxels, the recovery of deeply embedded objects still has a poor resolution. An attempt to achieve the high-resolution reconstruction in the voxel-based BLT was made through the modeling by radiative transfer equation (RTE) and the use of angularly-resolved data in the medium of a few mean free paths [14,15], although it may take quite an effort to acquire the angularly-resolved data in practice. Please note that the popular diffusion approximation (DA) cannot model the angularly-resolved measures.

In this work, instead of voxel representation, we model the source in Gaussians in reconstructions. A motivation is that because the light of interest is highly diffusive, a localized bioluminescent tumor does look like a Gaussian at a sufficiently large distance. Therefore it is practically sufficient to model various shapes as Gaussians in the diffusive medium, especially those deeply embedded objects. With this new representation, the reconstruction itself is unique [3] and the BLT is reduced to the problem with a few parameters. Moreover, the new model can be justified from physiological reasons. Please see the next section for the detailed explanation of the benefits for such a parameterization.

The outcome of the Gaussian model is formulated as a minimization problem with respect to the number of Gaussians, intensity, the center, the radius and the aligned direction of each Gaussian.

Just to clarify, although the number of Gaussians may be fixed in a specific reconstruction process, the exact number of actual Gaussian sources is not assumed. If the number of Gaussians in reconstruction is larger than the actual number, the source can still be recovered; otherwise, the actual number can be approached from a smaller initial guess in a combinatorial fashion through minimizing the data discrepancy. This can be illustrated in the following example. Suppose there are three Gaussians. In the reconstruction, we can set the number of Gaussians to four. With the proposed reconstruction algorithm, three recovered Gaussians would correspond to actual ones while the other one has the zero or nearly zero intensity values. On the other hand, if the assumed number is smaller than the actual number, e.g., only one Gaussian, the data discrepancy cannot be effectively reduced, and in this case the problem can be formulated as combinatorial optimization problem in which the number of Gaussians is increased iteratively until three so that data discrepancy can be reduced to the desired tolerance. On the other hand, a rough estimate of the actual number of Gaussians does improve both accuracy and efficiency of the reconstruction. In summary, the actual number of Gaussians is not assumed and it is practically sufficient to fix the number of Gaussians in reconstructions since the sources are localized and limited in quantity.

Please notice that, in [16] on diffuse optical imaging, modeling the absorption coefficient as a single Gaussian, the authors also demonstrated the reconstruction improvement when using *a priori* information from ultrasound images for the center of the Gaussian. In this work, we do not assume *a priori* information from other imaging modalities although it can be used for further improvement. For example, we do not assume the center of Gaussians that is to be recovered. However, we do impose some geometric constraints to avoid the apparent non-uniqueness of geometric parameterization. For example, we constrain on the distance between any two Gaussians so that they do not overlap, or use the size of the phantom as the bound of the center. Please refer to the next section for the details on the applied constraints.

Last, although we model the source in Gaussians in the reconstruction, the actual source distribution in simulations is not necessarily Gaussian. However the reconstruction is better if the actual source is in Gaussians as well. Please refer to the result section for details.

Here is the outline of the paper. In section 2, we introduce the new model of bioluminescent source by Gaussians, various intrinsic geometric constraints, and then the minimization algorithm; in section 3, the proposed model and algorithm is validated through simulations and experimental data.

## 2. Methods

#### 2.1. Forward Modeling

RTE is the most accurate among realistic models for light propagation in tissue [17,18]. However, due to its large dimensionality, an efficient solver of RTE is non-trivial [19]. Many efforts have been devoted in the development of RTE solver. Please see [20] and references therein for numerical methods for solving the deterministic RTE, and [21] and references therein for stochastic approaches by Monte Carlo method. The alternative modeling strategy is to simplify RTE into appropriate approximation according to particular optical regimes of interest [17,22,23]. Among all the approximate models, the most popular one is DA, which is valid in the scattering dominant regime. For simplicity, we will use DA in this proof-of-concept study. Please note that the modeling error in practice comes from not only the chosen model, but also the approximation of optical properties for the underlying medium.

In the following DA with Robin boundary condition, we use *f* for the light intensity, *q* for the bioluminescent source, *μ _{a}* for the absorption coefficient,

*μ*for the reduced scattering coefficient,

_{s}^{'}*D*for the diffusion coefficient, i.e.,

*D*= 1/ [3(

*μ*+

_{a}*μ*)], $\widehat{n}$ for the outgoing normal at the medium boundary,

_{s}^{'}*A*for the boundary constant [24].

Here we choose piecewise linear finite element method (FEM) to solve Eq. (2). After discretization, Eq. (2) is reduced to a linear system with the discretized light intensity *Φ*, the discretized source *Q*, and the system matrix *F* dependent on *μ _{a}* and

*D*, i.e.,

We refer the readers to [24] for the details of FEM solution of Eq. (2). After solving Eq. (3), we have the following boundary fluxes

#### 2.2. Source Representation by Gaussians

As a quantitative imaging method, BLT is expected to provide the source intensity or power accurately. However, there are some intrinsic aforementioned difficulties in achieving high-resolution reconstruction with the voxel-based BLT. Therefore, instead of local voxel representation of the bioluminescent source, we consider the global representation by characteristic shapes of the tumor, e.g., the Gaussians in this study. A motivation is that because the light of interest is highly diffusive, a localized bioluminescent tumor does look like a Gaussian if it is sufficiently distant from the boundary. It turns out that this new representation with significantly reduced degrees of freedom is more robust than the voxel-based representation.

Mathematically, it was proven in [3]: the solution of BLT is non-unique without the incorporation of effective *a priori* knowledge on the source distribution; in particular the uniqueness can be established when source is assumed to be composed of solid balls with the known intensities. However, using solid balls is not effective in representing general source distribution, especially features like anisotropy and fast decay of source intensity which are quite common in practice. In this study, we use summation of Gaussians to model general source distribution. This structural representation is flexible and improves the reconstruction stability for BLT. In particular the anisotropic Gaussian representation is an effective mathematical approximation in the sense that it parameterizes the major quantitative information of interest, such as the intensity, the center and the size of the source.

Physiologically, the Gaussian shape is indeed a natural choice. According to http://www.humpath.com/tumorigenesis, tumorigenesis is a collection of complex genetic diseases characterized by multiple defects in the homeostatic mechanisms that regulate cell growth, proliferation and differentiation. Cancer is caused by uncontrolled proliferation and the inappropriate survival of damaged cells, which results in tumor formation. Cells have developed several safeguards to ensure that cell division, differentiation and death occur correctly and in a coordinated fashion, both during development and in the adult body. Many regulatory factors switch on or off genes that direct cellular proliferation and differentiation. Damage to these genes, which are referred as tumor-suppressor genes and oncogenes, is selected for in cancer. Most tumor-suppressor genes and oncogenes are first transcribed from DNA into RNA, and are then translated into protein to exert their effects. Recent evidence indicates that small non-protein-coding RNA molecules, called microRNAs (miRNAs), might also function as tumour suppressors and oncogenes. Therefore, in the tumorigenesis stage, cancerous cells are localized, limited in quantity, irregular in shape, and with some decaying behavior away from the center. As a result, they can be modeled effectively as Gaussian sources after they are bioluminescently labeled. Other inflammation foci can be similarly modeled at earlier stages.

Therefore we represent the source *q* as the summation of Gaussians with the intensity *ρ*, the center (*x _{c}, y_{c}, z_{c}*), anisotropic radius (

*r*) and Euler angles (j, f, q), i.e.,

_{x}, r_{y}, r_{z}Computationally, the Gaussian source representation Eq. (5) reduces the burden significantly since we solve BLT with at most *n´(1 + d)(2 + d)/2* unknowns in *d* dimensions, which is in general much smaller than that from voxel representation, e.g., *N* (total number of voxels) unknowns for Eq. (1).

#### 2.3. Gaussian-based BLT

Modeling bioluminescent source in Gaussians, we formulate BLT as the minimization problem with respect to variables *X ^{G}* = {

*ρ*j

_{j}, x_{c,j}, y_{c,j}, z_{c,j}, r_{x,j}, r_{y,j}, r_{z,j},*, f*

_{j}*, q*

_{j}*},*

_{j}, j≤nHere we use the set {*f _{i}, i≤M*} for the measured data, the set {

*P*} of column vectors for the discretized measuring operator, and geometric constraints

_{i}, i≤M*R*on

*X*, which can be regarded as the regularization term. Please note that

^{G}*f*is now nonlinearly dependent on

*X*in this new formulation, although it linearly depends on

^{G}*q*in Eq. (1) in the voxel-based representation.

In Eq. (6), we do not treat *n* as a variable, mainly because there is generally no practical need for that since the tumors are usually localized and limited in quantity. We can just assign a practical estimate to *n* that is larger enough than or comparable with the actual number.

However, in the case that such an estimate is unavailable, we formulate Eq. (6) as a combinatorial optimization problem with respect to the variable *n*. That is we start from an initial guess *n _{0}* of

*n*, which can be smaller than the actual number, and solve Eq. (6) iteratively until the data discrepancy

*d*reaches certain tolerance

*e*. Here the update of

*n*is simply

*n*. For the notation convenience, we will use

_{k + 1}= n_{k}+ 1*X*for

*X*whenever it is appropriate. Then we have the following algorithm.

^{G}**Algorithm 1**: Combinatorial optimization**Given**: initial guess*n*and_{0}*e*= 0.05.- 2.
*d = ||f*_{k}*(**X*_{k}*)*–*f||/||f||;* - 3. Stopping criterion:
**Quit**if*d <e.* - 4.
*n*_{k + 1}= n_{k}+ 1

In Algorithm 1, *e = 0.05* is an empirical number that can identify those *n*’s that are sufficient for the reconstruction; *f _{k}* denotes the measurements with the source

*X*; || × || can be the simple summation of the absolute values.

_{k}Next, we turn to the geometric constraints *R*, which can be imposed naturally from the known geometry of the medium. Please note that we do not assume any *a priori* knowledge on the anatomical structure of the medium, e.g., from other imaging modalities, although they can be easily incorporated into the scheme to improve the reconstruction. All the geometric constraints are either from the shape of the medium or from some commonsense assumptions on the shape of the source to avoid some non-uniqueness in geometric representation.

The first type of constraints is the min-max constraint, the minimum and maximum for each variable of *X ^{G}*. For example, we can use the size of the medium as the min-max constraint for the center, and realistic values for anisotropic radius or intensity such as zeros as the lower bounds. That is

Notice that although we constraint the rotation angles, e.g., *0°*£q<*180°*, the representation of ellipsoid is not unique, e.g., by simultaneously exchanging *r _{x}* and

*r*and considering the supplementary angle of q in 2D. However, this does not affect the reconstruction result.

_{y}Secondly, we restrict the shape variation along different directions so that the source is not too “narrow”, which is physiologically reasonable, i.e.,

Here the parameter *c _{xy}*,

*c*and

_{yz}*c*control the size difference between any two directions for each source.

_{zx}In the case with multiple Gaussians, it is important to impose the following type of constraints, minimal-separation constaints, on centers and radiuses between any two Gaussians so that they will not overlap or get too close to each other. That is for any two Gaussian *i* and *j*,

Here *c _{x}*,

*c*and

_{y}*c*are the control parameters, which can for example be set to

_{z}*ln2*, corresponding to full width at half maximum (FWHM). Without using constraints from Eq. (7.2), multiple inclusions may not be separated. An example is given in the result section.

All the above constraints Eq. (7) can be represented by the following abstract inequalities

Here we use the popular logarithmic penalty functions to enforce those constraints, i.e.,

Combining Eq. (6) and (9), the BLT with Gaussian representation is formulated as the minimization of the sum of the data fidelity and the logarithmic penalty functions of geometric constraints. In the next section we will develop the algorithm for solving Eq. (6) via the popular barrier method.

Please note that when anatomical priors are available from other imaging modalities, such as X-ray CT or MRI, they can be incorporated naturally into the scheme Eq. (6). For example, the min-max constraints Eq. (7.1) can be tightened for centers or radiuses with the measured values from anatomical images; more accurate constraints can be imposed by assigning the appropriate constants in Eq. (7.2) and Eq. (7.3). In contrast, it is difficult to impose those global geometric features locally in the voxel-based BLT. On the other hand, the straightforward co-registration of the images from X-ray CT or MRI may not improve BLT satisfactorily since there is usually no simple one-to-one anatomical correspondence between them due to the fundamentally distinct imaging mechanisms.

#### 2.4. Minimization by Barrier Method

As the standard approach for minimizing nonlinear least squares, we linearize the first term (the data fidelity) in Eq. (6), and solve for the incremental change *dX* iteratively via the following outer loop of iterations

Here *f ^{n}* is the simulated data with

*X*,

^{n}*J*is the Jacobian coming from the linearization, which depends on

*X*, and the detail for the computation of

^{n}*J*is given in the next section for the algorithm implementation.

During each iteration in the outer loop for Eq. (10), we minimize the following with *b = f-f ^{n}* and

*x = dX*

Here we choose the popular barrier method [25] to enforce the constraint term *R*. That is we solve a sequence of the following minimization problems

Minimizing Eq. (12) strictly enforces the geometric constraint *R* since the value of logarithmic penalty functions would otherwise become infinitely large. During each step the solution *x ^{n}* is no more than

*K/t*-suboptimal with

*K*as the total number of constraints. This implies

*x*converges to the optimal point of Eq. (11) such that

^{n}*Jx = b*with strictly enforced geometry constraints as

*t*→∞. The Eq. (12) can be minimized with Newton’s method via the following inner loop

In each iterative step by Eq. (13), we compute the descent direction *dx*, find the moving step *s* through the backtracking line search, and then update *x* and *t* for the next iteration until the stopping criterion is satisfied. That is we solve a sequence of *t*-subproblems via Eq. (13) with the increasing *t*.

#### 2.5. Algorithm Implementation

The flowchart of the algorithm for solving Eq. (6) is as follow: the outer loop comes from the linearization Eq. (10) of the data fidelity term while the inner loop solves each linearization step using barrier method via Eq. (13).

**Outer loop**: Linearization of the data fidelity via Eq. (10)**Given**: initial guess*X*and^{0}*e*= 0.01._{o}**Repeat**: 1. Compute Jacobian*J*from*X*.^{n}**Given**:*t*,^{0}= −2R(X^{n})/||b||^{2}*μ = 2*, initial guess*x*.^{0}= 0**Repeat**: 2.1. Compute the descent direction*dx*;- 2.2. Compute the moving step
*s*via backtracking line search; **Given**:*a = 0.01, s = 1, b = 0.5*.**While**:*L(x*^{n}+ sdx)>L(x^{n}) + α(ÑL_{x})^{T}× (sdx)*s =*b*s*.- 2.3. Update
*x*+^{n + 1}= x^{n}*sdx*and*t*.^{n + 1}= μt^{n} - 2.4. Stopping criterion:
**Quit**if*K/t*.^{n}<e_{i} - 3. Stopping criterion:
**Quit**if*|E*^{n + 1}-E^{n}|/ E^{n}<e_{o}.

In this flowchart, we use the ratio difference of the total source power *E = òqdΩ* as the stopping criterion for the outer loop, i.e., assuming that the total power *E* is stable when the algorithm converges. As the stopping criterion for the inner loop, we use the sub-optimal ratio *K/t*, which naturally measures the difference between the iterative solution of Eq. (12) and the original solution with strictly enforced constraints. We set *e _{i}* empirically to

*0.0001t*. Please note that the formula for

^{0}*t*is for balancing the linearized data fidelity and the penalty function.

^{0}Next, we give the details of computing Jacobian and the descent direction in Algorithm 2 for the completeness.

Let us start with the computation of Jacobian *J = {J _{ij}, i£M, j£n}*, which can be derived from the Jacobian

*J*with respect to the piecewise-constant voxel representation via Eq. (1). The connection between

_{0}= {J_{0,ij}, i£M, j£N}*J*and

*J*is through the following coordinate transformation

_{0}*q*represents the

_{k}*k*th voxel value in Eq. (1),

*x*is the

_{j}*j*th parameter in global representation by Gaussians in Eq. (5).

An efficient way to compute *J _{0}* is through the adjoint method [24]. After computing

*J*,

_{0}*J*follows immediately from (14). However, noticing that there are only a few parameters in the new source representation by Gaussians, i.e.,

*n†M†N*, we directly compute

*J*without computing

*J*via the direct method as follow.

_{0}From Eq. (3), we have *F(∂ϕ/∂q _{k}) = ∂Q/∂q_{k}*, then

Through Eq. (15), we only need to compute the forward solver *n* times, which would be at least *M* times if we compute *J _{0}* first via Eq. (14). Please note that

*∂Q/∂q*is a sparse vector, which is nonzero only at the nodes in the

_{k}*k*th voxel; on the other hand

*∂q*can be analytically derived from Eq. (5). Therefore, we can now compute the Jacobian

_{k}/∂x_{j}*J*merely by computing

*n*forward solvers, which is much more efficient than that in the voxel-based BLT. The trade-off is that we need to compute the Jacobian more than once due to the nonlinear dependence of

*J*on

*X*, which needs to be computed only once in the voxel-based BLT due to the linear dependence of

^{G}*J*on

_{0}*X*. However, since

*n†M*, we still achieve a considerable gain in speed, since the new algorithm usually converges in less than 20 iterations while the ratio

*M/n*is usually at least about a hundred. That is we usually achieve at least 5 gains in speed on the computation of Jacobians.

In the rare cases when *n>M*, the adjoint method is preferred, i.e.,

*M*forward solvers instead.

Next we compute the first and second derivatives of logarithmic penalty functions for each geometric constraint *g _{k}(X^{G})* as in Eq. (7), i.e.,

*R*. From the straightforward computation, we have the following simple formula

_{k}= -ln[-g_{k}(X^{G})]Therefore, we have

*R*is evaluated at

_{k}*X*rather than

^{n}+ x*x*.

The overall algorithm for BLT with Gaussian representation is more efficient than voxel-based BLT, mainly because the minimization is now with respect to only a few parameters of Gaussians instead of at least thousands of voxels in the conventional voxel-based BLT. For the current BLT algorithm with new model by Gaussians, minimization through barrier method in Step 2 is extremely fast although the algorithm seems complicated, which is understandable due to the significantly reduced number of variables; the major computation cost comes from the Jacobian *J* in Step 1, although it takes less computational time than that in voxel-based BLT, which is again due to the reduced number of unknowns.

## 3. Results

In this section, we first compare the performance of Gaussian-based BLT, “G-BLT”, with that of the conventional voxel-based BLT, “V-BLT”, through simulated data in single-inclusion and multiple-inclusion cases, and then validate the proposed method in mouse study.

#### 3.1. Reconstruction with Single Inclusion

We first evaluate the proposed algorithm with single anisotropic Gaussian inclusion, and then perform the single-inclusion reconstructions with various shapes, different noise levels, and different approximation errors in optical background to compare the performance of G-BLT and V-BLT. For the presentation purpose, we only show 2D comparisons. The results and the conclusions are similar for the 3D case. In particular, 3D results will be demonstrated with *in vivo* data.

All the 2D reconstructions are performed on a circular phantom with the diameter of 20 millimeters. 120 data are collected all around the phantom at the boundary with equal angular distance in between. The reconstruction mesh has 2162 triangular elements and 1130 nodes. To avoid the inverse crime, the data are generated with different meshes that are about four times as large as the reconstruction mesh. For example, there are 8664 triangular elements and 4425 nodes for the case with circular inclusion as shown in Fig. 1(a) .

We assume the homogeneous optical background with *μ _{a} = 0.01mm^{−1}* and

*μ*. Without further mentioning, we use

_{s}^{'}= 1mm^{−1}*mm*as the unit for the length,

*mm*as the unit for absorption or scattering coefficients, nano Watts (

^{−1}*nW*) as the unit for the total source power, and

*nW/mm*,

^{2}*nW/mm*as the unit for the source intensity in 2D and 3D respectively. The unit of angle is in degree.

^{3}In G-BLT, as the initial guess, for simplicity we assume there is one Gaussian inclusion with *ρ = 0.1*, *x _{c} = 0*,

*y*,

_{c}= 0*r*. The reconstruction with multiple inclusions as the initial guess will be considered in the next section. Please notice that the nonzero value has to be assigned for

_{x}= 1, r_{y}= 1 and q = 90°*ρ*, otherwise the Jacobian with respect to other parameters will be zero according to Eq. (14). For the constraint variables, we set

*ρ*,

_{max}= 10, ρ_{min}= 0.01, x_{max}= y_{max}= 10, x_{min}= y_{min}= −10, r_{x,max}= r_{y,max}= 5, r_{x,min}= r_{y,min}= 0.1*q*, and

_{max}= 180°,q_{min}= 0°*c*.

_{xy}= 5For the V-BLT, we solve the problem via the standard Levenberg-Marquardt method with L_{2} regularization [24]. That is to minimize

*J*is with respect to the piecewise-constant voxel representation by Eq. (1), and

_{0}*l*is the regularization parameter.

Please notice that one may optimize the choice of regularization term in Eq. (19). However we found out that the results did not improve drastically when using other regularization strategies [14,15], which may be due to the severely illposed nature of the inverse problem in the optical regime modeled by diffusion approximation Eq. (2). In the following, we restrict our discussion based on L2 regularization as in Eq. (19).

### 3.1.1. Inclusions with Gaussian Shapes

Since it is physiologically reasonable to approximate the bioluminescence sources by Gaussians, we first consider reconstructions where we assume the source is a single Gaussian. All the tests have the same Gaussian up to a rotation, with parameters specified in Table 1 . The reconstructed results with both G-BLT and V-BLT are plotted in Fig. 1 and the recovered parameter values are in Table 1 as well.

From Fig. 1, in all tests G-BLT can recover the Gaussian phantom when using the representation of source in Gaussians while V-BLT gives relatively much poorer results. Moreover, G-BLT is quantitatively better than V-BLT in the sense it gives more accurate total power *E* and maximal intensity as documented in Table 1.

### 3.1.2. Inclusions with Various Shapes

We evaluate the algorithm through the simulated data with single inclusions of various shapes. The parameters for various single inclusions are specified in Table 2 . The reconstructed results with both G-BLT and V-BLT are plotted in Fig. 2 and the parameter values are in Table 3 .

Before doing the reconstruction, we simply compare the measurements from Gaussian and non-Gaussian sources as shown in Fig. 3 . Here we embed inclusions deeply at the center of the medium. The plotted measurements are normalized with the absolute sum 1. The figure shows that the measurement profiles from Gaussian sources are almost the same as those from the non-Gaussian ones. It suggests that various shapes can be characterized by anisotropic Gaussians when the sources are deeply embedded into the diffusive medium.

From Fig. 2, visually we can tell that G-BLT is consistently better in localizing the inclusions than V-BLT. Moreover, G-BLT is quantitatively better than V-BLT in the sense it gives more accurate total power *E* and maximal intensity as documented in Table 3.

One minor issue regarding G-BLT is that the reconstructed radiuses along x- and y- axis seem not to match those in true phantoms, whish is possibly due to the Gaussian approximation error of non-Gaussian phantom. For example, with rectangular inclusion, the recovered *r _{x} = 1.752* and

*r*, while the true

_{y}= 1.797*r*and

_{x}= 1*r*. Although

_{y}= 2*r*in the reconstruction, this scale difference is much smaller than the true difference.

_{x}<r_{y}Notice that in this section and the rest of reconstructions for non-Gaussian sources, we did not consider *q* as reconstruction variable. The reason is that G-BLT does not recover the shapes for non-Gaussian sources as exactly as for Gaussian sources even with rotation.

### 3.1.3. Different Noise Level

We evaluate the algorithm with respect to different noise level. 5%, 10%, 20%, 30% percentage of Gaussian noise is added to the measurements *f* respectively, e.g., *f(1 + 5% × Randn)*, where *Randn* represents the random Gaussian distribution with zero mean and unit variance. The reconstructed results with both G-BLT and V-BLT are plotted in Fig. 4
and the parameter values are in Table 4
. Here the simulations are performed on the same circular phantom with single circular inclusion as in Fig. 2(a).

Again, from Fig. 4, visually G-BLT localizes much better than V-BLT. In particular, the result from G-BLT changes a little for the noise up to 20%. On the other hand, it is difficult to tell even the center of object in V-BLT when the noise is increased beyond 10%. Moreover, from Table 4, G-BLT gives more accurate quantitative information than V-BLT.

### 3.1.4. Mismatch of Optical Background

We evaluate the algorithm with different mismatch in optical background. The background is perturbed with ± 10%, ± 30%, ± 50%, ± 70% percentage of error respectively, e.g., *μ _{a}(1 + 10%)* and

*μ*. Notice that we either increase or decrease both absorption and scattering coefficients since empirically we find that the mismatch may cancel if one is increased while the other is decreased. Therefore, we change both in the same direction to simulate the worst-case errors. Please also note that the reconstruction results with V-BLT are not presented since it does not localize the inclusion satisfactorily even in the case without optical background error.

_{s}^{'}(1 + 10%)The results from G-BLT are presented in Fig. 5 and Table 5 , which show that G-BLT is robust to the mismatch error of optical background in the sense that the inclusion can still be localized well despite of the shifting of centers due to the systematic discrepancy between the true values and the values used in the reconstruction.

#### 3.2. Reconstruction with Multiple Inclusions

Now we evaluate the algorithm in the presence of multiple inclusions by considering first Gaussian inclusions and then non-Gaussian shapes.

Despite of the fact that there are two inclusions in the test on Gaussian sources (Fig. 6
) and three inclusions in the test on non-Gaussian sources (Fig. 6), we assume four inclusions as the initial guess in G-BLT, since the exact number of inclusions is usually not known in practice unless the correct anatomical priors are available from other imaging modalities. However, since the bioluminescent source is localized and limited in quantity as discussed previously, we only need to set the number of sources *n* to a reasonable larger value than the actual number, in order to faithfully reconstruct the source distribution without considering *n* as an independent variable. Particularly, as the initial guess, we assume there are four Gaussian inclusions, i.e., *n = 4*, with *ρ _{1} = ρ_{2} = ρ_{3} = ρ_{4} = 0.1*,

*x*,

_{1,c}= 5*y*,

_{1,c}= 5*x*,

_{2,c}= 5*y*,

_{2,c}= −5*x*,

_{3,c}= −5*y*,

_{3,c}= −5*x*,

_{4,c}= −5*y*,

_{4,c}= 5*r*,

_{1,x}= r_{2,x}= r_{3,x}= r_{4,x}= 1*r*and

_{1,y}= r_{2,y}= r_{3,y}= r_{4,y}= 1*q*. Notice that for the same reason mentioned above, we do not consider

_{1}= q_{2}= q_{3}= q_{4}= 90°*q*as reconstruction variable for non-Gaussian sources for simplicity.

The same constants as in the single-inclusion simulations are used. In order to separate the multiple objects, we find out the minimal-separation constraints are usually necessary. This is also intuitively correct since the source representation may otherwise be non-unique, given the severely illposed nature of BLT. Here we set *c _{x}* =

*c*=

_{y}*ln2*, corresponding to that the minimal distance between any two objects is at least the average of FWHM.

### 3.2.1. Inclusions with Gaussian Shapes

The detailed setting is given in Table 6 and Fig. 6(a). The results are plotted in Fig. 6, which clearly shows that G-BLT is able to accurately recover the values of intensity, centers, radiuses and rotation angle. On the other hand, neither shapes nor locations of inclusions can be not be captured from V-BLT.

### 3.2.2. Inclusions with Non-Gaussian Shapes

The detailed setting is given in Table 7 and Fig. 7(a) . The results are plotted in Fig. 7, which clearly shows that G-BLT with proper constraints is able to separate the inclusions while V-BLT fails. The parameter values are documented in Table 7, which are quite consistent with the true values. One minor issue is that the recovered center is not as accurate as one from the single-inclusion cases, which is understandable due to mutual effect of inclusions and the illposedness of the problem.

### 3.2.3. Combinatorial optimization

In this section, we give one example of combinatorial optimization (Algorithm 1) with the number of Gaussians *n* as a variable.

In the phantom [Fig. 8(a)
], there are three Gaussians with the parameters specified in Table 8
. In the reconstruction, we start with the initial guess n_{0} = 1, and continue the combinatorial process until the data discrepancy *d* reaches the tolerance. Please notice that the process should have terminated at *n = 3*. For illustration purpose, we run until *n = 5*.

The detailed setting and the recovered parameters are given in Table 8. Figure 8 shows that this combinatorial process is able to reconstruct the source even the initial guess is smaller than the actual number, e.g., the results from *n = 1* and *n = 2* fail to give a satisfactory reconstruction while those from *n = 3*, *n = 4* and *n = 5* all successfully recover the desired distribution.

It is interesting to notice that although the reconstruction with more than or equal to the actual number of Gaussians is able to recover the actual number of Gaussians, the reconstruction with more Gaussians than the actual number may degrade the image quality due to the apparent non-unique representation of the geometry. For example, two recovered sources in Fig. 8(e) constitute the top Gaussian in Fig. 8(a) when using *n = 4* and two recovered sources in Fig. 8(f) constitute the left Gaussian in Fig. 8(a) when using *n = 5.*

#### 3.3. 3D in vivo validation

After evaluating the propose G-BLT in 2D with various settings, we have established the superiority of G-BLT over V-BLT. Now we validate G-BLT in 3D with *in vivo* experimental data. The details of experimental setups and mouse experiments are given in [4]. We also refer the readers to the above reference for the reconstruction parameters, such as the used mesh and the values for optical parameters.

In the previous voxel-based reconstruction [4], with a permissible region selected according to the high value clusters in the data, two sources were localized with a stronger one of the power 39.8nW/mm^{3} (right) and a smaller one of the power 1.5nW/mm^{3} (left). When the mouse was dissected after the experiment, two tumors were found on both adrenal glands, respectively, as shown in Fig. 9(c)
. The volume of tumor tissues as measured by Vernier calipers was 468 mm^{3} for the tumor (right) 275 mm^{3} for the tumor (left).

Now we use G-BLT, Gaussian-based BLT, to reconstruct instead of V-BLT. Please note that we do not select the permissible region.

As the initial guess, we assume there are four Gaussian inclusions, i.e., *n = 4*, with *ρ _{1} = ρ_{2} = ρ_{3} = ρ_{4} = 10*,

*x*,

_{1,c}= x_{2,c}= x_{3,c}= x_{4,c}= 0*y*,

_{1,c}= y_{2,c}= y_{3,c}= y_{4,c}= 0*z*,

_{1,c}= 5*z*,

_{2,c}= 10*z*,

_{3,c}= 15*z*,

_{4,c}= 20*r*,

_{1,x}= r_{2,x}= r_{3,x}= r_{4,x}= 1*r*and

_{1,y}= r_{2,y}= r_{3,y}= r_{4,y}= 1*r*.

_{1,z}= r_{2,z}= r_{3,z}= r_{4,z}= 1The constants in the geometric constraints are: *ρ _{max} = 1000, ρ_{min} = 0.1, x_{max} = y_{max} = 12, x_{min} = y_{min} = −12, z_{max} = 27, z_{min} = 0, r_{x,max} = r_{y,max} = 6, r_{x,min} = r_{y,min} = 0.5*, and

*c*=

_{xy}= c_{yz}= c_{zx}= 5, c_{x}*c*=

_{y}*c*=

_{z}*ln2*. Here the min-max values on the geometry correspond to the physical size of the domain; the maximum of the intensity is estimated from the data, which is roughly 50 times of the maximal value in the data.

As documented in Table 9
, 2 out of 4 recovered inclusions, with the locations just corresponding to the kidneys, have significantly larger values, i.e., the one with the peak intensity of 16.593nW/mm^{3} (right) and the other one with the peak intensity of 7.692nW/mm^{3} (left). In addition, the peak intensity of the other 2 inclusions (Inclusion 3 and 4) is less than 5% of the maximum (Inclusion 1), which can be from experimental or numerical noise. The reconstructed bioluminescent volumes (*2r _{x} × 2r_{y} × 2r_{y}*) are 506 mm

^{3}(right) and 901 mm

^{3}(left) respectively. The discrepancy between the reconstructed volumes and the measured ones is not surprising since the anatomical volume may not correspond to the bioluminescent volume.

## 4. Conclusions and discussions

In this proof-of-concept study, we have introduced a stable way to recover bioluminescent source through novel source representation by Gaussians. It is mathematically unique and computationally cheap. Besides, there are strong physiological evidences that the actual bioluminescent sources can be modeled by Gaussians, which will be further studied in the future work. The resultant problem with the fixed number of Gaussians is formulated as a nonlinear least-square minimization with appropriate geometric constraints, which is solved in turn by barrier method; on the other hand, the problem with the varying number of Gaussians can be formulated as a combinatorial optimization problem. The superiority of the proposed method has been established through simulations and *in vivo* experimental data. In particular, when the source itself can be approximated by Gaussians, the proposed method is able to accurately recover the intensity, centers, radiuses and rotation angles of the Gaussians.

Due to the significantly reduced number of variables, the proposed method is robust to the measurement noise and the mismatch in optical background; the computational burden is also reduced to the minimization problem with respect to a few parameters instead of at least thousands of parameters in the conventional voxel-based BLT.

In the future work, we shall perform more mouse studies to further develop and optimize the current method into a practical *in vivo* bioluminescence imaging method, for which we may need to incorporate RTE for more accurate modeling, simplified representation for optical heterogeneity, such as piecewise constants, to be simultaneously reconstructed for correcting optical background, and multi-spectral data for better stability.

Although the source is represented as Gaussians in this study, other shape functions may be appropriate for other purposes. An interesting topic to pursuit is the multi-modality imaging since the structural priors can be easily incorporated into the proposed method for more accurate quantitative estimation.

## Acknowledgments

The work is partially supported by National Science Foundation (NSF) grant DMS0811254, National Institutes of Health (NIH)/National Heart, Lung and Blood Institute (NHLBI) grant HL098912 and NIH/National Cancer Institute (NCI) grant CA135151.

## References and links

**1. **C. H. Contag and B. D. Ross, “It’s not just about anatomy: in vivo bioluminescence imaging as an eyepiece into biology,” J. Magn. Reson. Imaging **16**(4), 378–387 (2002). [PubMed]

**2. **G. Wang, E. A. Hoffman, G. McLennan, L. V. Wang, M. Suter, and J. Meinel, “Development of the first bioluminescent CT scanner,” Radiology **299**, 566 (2003).

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

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

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

**6. **G. Wang, X. Qian, W. Cong, H. Shen, Y. Li, W. Han, K. Durairaj, M. Jiang, T. Zhou, and J. Cheng, “Recent development in bioluminescence tomography,” Curr. Med. Imaging Rev. **2**, 453–457 (2006).

**7. **G. Wang, H. Shen, K. Durairaj, X. Qian, and W. Cong, “The first bioluminescence tomography system for simultaneous acquisition of multi-view and multi-spectral data,” Int. J. Biomed. Imaging **2006**, 1–8 (2006).

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

**9. **A. J. Chaudhari, F. Darvas, J. R. Bading, R. A. Moats, P. S. Conti, D. J. Smith, S. R. Cherry, and R. M. Leahy, “Hyperspectral and multispectral bioluminescence optical tomography for small animal imaging,” Phys. Med. Biol. **50**(23), 5421–5441 (2005). [PubMed]

**10. **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). [PubMed]

**11. **C. Kuo, O. Coquoz, T. Troy, D. Zwarg, and B. Rice, “Bioluminescent tomography for in vivo localization and quantification of luminescent sources from a multiple-view imaging system,” Mol. Imaging **4**, 370 (2005).

**12. **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). [PubMed]

**13. **Y. Lv, J. Tian, W. Cong, and G. Wang, “Experimental study on bioluminescence tomography with multimodality fusion,” Int. J. Biomed. Imaging **2007**, 86741 (2007). [PubMed]

**14. **H. Gao and H. K. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation Part 1: l1 regularization,” Opt. Express **18**(3), 1854–1871 (2010). [PubMed]

**15. **H. Gao and H. K. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation part 2: total variation and l1 data fidelity,” Opt. Express **18**(3), 2894–2912 (2010). [PubMed]

**16. **J. Liu, A. Li, A. E. Cerussi, and B. J. Tromberg, “Parametric diffuse optical imaging in reflectance geometry,” IEEE Sel. Top. Quantum Electron. **16**, 555–564 (2010).

**17. **K. M. Case, and P. F. P. F. Zweifel, *Linear Transport Theory* (Addison-Wesley Educational Publishers Inc., 1967).

**18. **A. Ishimaru, *Wave Propagation and Scattering in Random Media* (Academic Press, 1978).

**19. **E. E. Lewis, and W. F. Miller, *Computational Methods of Neutron Transport* (Wiley, 1984).

**20. **H. Gao and H. K. Zhao, “A fast forward solver of radiative transfer equation,” Transp. Theory Stat. Phys. **38**, 149–192 (2009).

**21. **H. Shen and G. Wang, “A tetrahedron-based inhomogeneous Monte Carlo optical simulator,” Phys. Med. Biol. **55**(4), 947–962 (2010). [PubMed]

**22. **W. Cong, H. Shen, A. Cong, Y. Wang, and G. Wang, “Modeling photon propagation in biological tissues using a generalized Delta-Eddington phase function,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **76**(5 Pt 1), 051913 (2007). [PubMed]

**23. **A. D. Klose and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comput. Phys. **220**, 441–470 (2006).

**24. **A. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. **15**, R41–R93 (1999).

**25. **S. Boyd, and L. Vandenberghe, *Convex Optimization* (Cambridge university press, 2004).