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

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, μ_s^' for the reduced scattering coefficient, D for the diffusion coefficient, i.e., D = 1/ [3(μ_a + μ_s^')], $\widehat{n}$ for the outgoing normal at the medium boundary, 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_x, r_y, r_z) and Euler angles (j, f, q), i.e.,

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, x_c,j, y_c,j, z_c,j, r_x,j, r_y,j, r_z,j, j_j, f_j, q_j, j≤n},

Here we use the set {f_i, i≤M} for the measured data, the set {P_i, i≤M} of column vectors for the discretized measuring operator, and geometric constraints R on X^G, which can be regarded as the regularization term. Please note that f is now nonlinearly dependent on X^G in this new formulation, although it linearly depends on 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₀ of n, which can be smaller than the actual number, and solve Eq. (6) iteratively until the data discrepancy d reaches certain tolerancee. Here the update of n is simply n_{k + 1} = n_k + 1. For the notation convenience, we will use X for X^G whenever it is appropriate. Then we have the following algorithm.

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_k; || × || can be the simple summation of the absolute values.

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_y and considering the supplementary angle of q in 2D. However, this does not affect the reconstruction result.

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_yz and c_zx control the size difference between any two directions for each source.

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_y and c_z are the control parameters, which can for example be set to 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ⁿ is the simulated data with Xⁿ, J is the Jacobian coming from the linearization, which depends on Xⁿ, and the detail for the computation of 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ⁿ 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ⁿ 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 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).

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 t⁰ is for balancing the linearized data fidelity and the penalty function.

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₀ = {J_0,ij, i£M, j£N} with respect to the piecewise-constant voxel representation via Eq. (1). The connection between J and J₀ is through the following coordinate transformation

An efficient way to compute J₀ is through the adjoint method [24]. After computing J₀, 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.

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₀ first via Eq. (14). Please note that ∂Q/∂q_k is a sparse vector, which is nonzero only at the nodes in the kth voxel; on the other hand ∂q_k/∂x_j can be analytically derived from Eq. (5). Therefore, we can now compute the Jacobian 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^G, which needs to be computed only once in the voxel-based BLT due to the linear dependence of J₀ on 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.,

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_k = -ln[-g_k(X^G)]. From the straightforward computation, we have the following simple formula

Therefore, we have

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) .

 

Fig. 1 Reconstructions with single Gaussian inclusion. (a)-(d) are the phantoms with 0°, 45°, 90° and 135° rotations respectively; (e)-(h) are from the voxel-based BLT; (i)-(l) are from Gaussian-based BLT. Please see Table 1 for the true and recovered parameters. Please note that the maximum intensity recovered via the voxel-based BLT is only up to 20% of the true value. The results show that Gaussian-based BLT not only localizes the inclusion better than, but also quantitatively more accurate than the voxel-based BLT.

Download Full Size | PDF

We assume the homogeneous optical background with μ_a = 0.01mm⁻¹ and μ_s^' = 1mm⁻¹. Without further mentioning, we use mm as the unit for the length, mm⁻¹ as the unit for absorption or scattering coefficients, nano Watts (nW) as the unit for the total source power, and nW/mm², nW/mm³ as the unit for the source intensity in 2D and 3D respectively. The unit of angle is in degree.

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_x = 1, r_y = 1 and q = 90°. 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 ρ, 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_max = 180°,q_min = 0°, and c_xy = 5.

For the V-BLT, we solve the problem via the standard Levenberg-Marquardt method with L₂ regularization [24]. That is to minimize

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.

Table 1. Tests with single Gaussian source. The Gaussian phantom has ρ = 1, x_c = 5, y_c = 0, r_x = 2, r_y = 1. The rotation angles are 0°, 45°, 90°, 135°respectively for Test 1-4. In all tests, the total energy of the phantom E is approximately 6.285. Please notice that the result from Case 2 can be understood according to the non-unique representation of the Gaussian

View Table | View all tables in this article

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 .

Table 2. The parameters for single-object simulations with inclusions of various shapes. E is the total source energy; ρ is the maximal source intensity; x_c and y_c are x- and y-coordinate of the center of the inclusion; r_x and r_y are maximum distances from the center to the boundary along x and y direction respectively

View Table | View all tables in this article

 

Fig. 2 Single-inclusion reconstructions with various shapes. (a)-(d) are the phantoms with circular, square, triangular and rectangular inclusion respectively; (e)-(h) are from the voxel-based BLT; (i)-(l) are from Gaussian-based BLT. Please see Table 2 and 3 for the true and recovered parameters. Please note that the maximum intensity recovered via the voxel-based BLT is only up to 20% of the true value. The results show that Gaussian-based BLT not only localizes the inclusion better than, but also quantitatively more accurate than the voxel-based BLT.

Download Full Size | PDF

Table 3. The reconstructed parameters for single-object simulations with inclusions of various shapes. Please see Table 2 for the corresponding true values

View Table | View all tables in this article

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.

 

Fig. 3 Comparison of boundary measurements from Gaussian sources and non-Gaussian sources. (a)-(c) are from non-Gaussian sources with side length ratio of 1:1, 1:2 and 2:1. (d)-(f) are from Gaussian sources with the radius ratio of 1:1, 1:2 and 2:1. (g)-(i) plot the boundary measurements originating counterclockwise from x-axis, in which the red curve are from Gaussian sources and the blue curves are from non-Gaussian sources. Please note that the data are normalized so that the absolute sum is 1.

Download Full Size | PDF

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_y = 1.797, while the true r_x = 1 and r_y = 2. Although r_x<r_y in the reconstruction, this scale difference is much smaller than the true difference.

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).

 

Fig. 4 Single-inclusion reconstructions with different Gaussian noise levels. The phantom is as the same as (a) in Fig. 2 with single circular inclusion; (a)-(d) are from the voxel-based BLT with 5%, 10%, 20% and 30% noise level respectively; (e)-(h) are correspondingly from Gaussian-based BLT. Please see Table 2 and 4 for the true and recovered parameters. The results show that Gaussian-based BLT is more robust to the noise than the voxel-based BLT.

Download Full Size | PDF

Table 4. The reconstructed parameters for single-circular-object simulations with different Gaussian noise levels. Please see Table 2 for the corresponding true values

View Table | View all tables in this article

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 μ_s^'(1 + 10%). 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.

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.

 

Fig. 5 Single-inclusion reconstructions with different mismatch in optical background. The phantom is as the same as (a) in Fig. 2 with single circular inclusion; (a)-(d) are from Gaussian-based BLT with −10%, −30%, −50%, −70% error in optical background; (e)-(h) are from Gaussian-based BLT with + 10%, + 30%, + 50%, + 70% error in optical background. Please see Table 2 and 5 for the true and recovered parameters. Please also note that the reconstruction results with the voxel-based BLT are not presented since it does not localize the inclusion satisfactorily even in the case without optical background error. The results show that Gaussian-based 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.

Download Full Size | PDF

Table 5. Gaussian-based BLT for single-circular-object simulations with different mismatch in optical background. Please see Table 2 for the corresponding true values

View Table | View all tables in this article

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 ρ₁ = ρ₂ = ρ₃ = ρ₄ = 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_1,y = r_2,y = r_3,y = r_4,y = 1 and q₁ = q₂ = q₃ = q₄ = 90°. Notice that for the same reason mentioned above, we do not consider q as reconstruction variable for non-Gaussian sources for simplicity.

 

Fig. 6 Multiple-inclusion reconstructions with 5% Gaussian noise level. (a) is the phantom with two elliptical inclusions.; (b) is from the voxel-based BLT; (c) is from Gaussian-based BLT. Please see Table 6 for the true and recovered parameters. The results show again that Gaussian-based BLT not only localizes the inclusion better than, but also quantitatively more accurate than the voxel-based BLT.

Download Full Size | PDF

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.

Table 6. Multiple-Gaussian-inclusion simulation with Inclusion 1 [the left ellipse in Fig. 6(a)] and Inclusion 2 [the right ellipse in Fig. 6(a)]. The total energy E of the phantom is 12.570; the recovered total energy from G-BLT is 12.596 while that from V-BLT is 11.582

View Table | View all tables in this article

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.

Table 7. Multiple-inclusion simulations. Three circular inclusions of 1mm diameter are located at (0, 0), (4, 0) and (8, 0) respectively with the unit source intensity. The total energy E is 9.364; the recovered total energy from G-BLT is 9.290 while that from V-BLT is 8.282. Please note that the presented values are from G-BLT with minimal-separation constraints as in Eq. (7.3). Please also notice that three inclusions are assumed in the initial guess for the reconstruction with G-BLT although there are actually three inclusions. After the reconstruction, 3 out of 4 inclusions have distinguished peak values

View Table | View all tables in this article

 

Fig. 7 Multiple-inclusion reconstructions with 5% Gaussian noise level. (a) is the phantom with three circular inclusions.; (b) is from the voxel-based BLT; (c) is from Gaussian-based BLT without minimal-separation constraint; (d) is from Gaussian-based BLT with the minimal-separation constraint. Please see Table 7 for the true and recovered parameters. The results show that Gaussian-based BLT with proper geometric constraints is in general able to separate multiple inclusions better than the voxel-based BLT.

Download Full Size | PDF

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₀ = 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.

 

Fig. 8 Combinatorial optimization with the varying n for G-BLT. (a) is the phantom with three Gaussian inclusions.; (b)-(f) are the reconstructed images from n = 1, 2, 3, 4, 5 respectively.

Download Full Size | PDF

Table 8. Combinatorial optimization with the varying n for G-BLT. Three Gaussian inclusions are located at (0, 5), (4.33, −2.5) and (−4.33, −2.5) respectively with ρ = 1, r_x = 1 and r_y = 2. The number of Gaussians n is treated as a variable in this combinatorial problem. The value of n in the reconstruction is set to 1 as the initial guess and iteratively increased until the solution converges. The following table records the reconstructed values for n = 1, 2, 3, 4 and 5. Please notice that the reconstruction should have terminated at n = 3 according to the given criterion if it were not for illustration purpose. The total energy E is 18.843; the recovered total energy is 21.164, 20.861, 18.879, 18.908 and 18.912 respectively. The measurement difference d is 0.167, 0.132, 0.00315, 0.00322 and 0.00327 respectively

View Table | View all tables in this article

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³ (right) and a smaller one of the power 1.5nW/mm³ (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³ for the tumor (right) 275 mm³ for the tumor (left).

 

Fig. 9 In vivo mouse studies on Gaussian-based BLT. (a) shows two reconstructed bioluminescent sources with dominant power on two kidneys respectively; (b) is the reconstructed bioluminescent isosurface of the intensity value 7; (c) is the histological verification with two tumors at the same locations on the dissected kidneys. Please see Table 8 for the recovered parameters.

Download Full Size | PDF

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 ρ₁ = ρ₂ = ρ₃ = ρ₄ = 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_1,y = r_2,y = r_3,y = r_4,y = 1 and r_1,z = r_2,z = r_3,z = r_4,z = 1.

The 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³ (right) and the other one with the peak intensity of 7.692nW/mm³ (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³ (right) and 901 mm³ (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.

Table 9. In vivo study using Gaussian-based BLT. Please also notice that four inclusions are assumed in the reconstruction although there are actually two tumors from the histological verification. The reconstructed result also shows exactly two inclusions with the distinguished values, i.e. Inclusion 1 and 2 in the following table

View Table | View all tables in this article

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.