Hybrid reconstruction framework for model-based multispectral bioluminescence tomography based on Alpha-divergence

1. Introduction

Bioluminescence tomography (BLT) is an optical molecular imaging modality based on the principle of localizing sources in organisms by surface measurements, with high sensitivity, low cost, and noninvasive features.^1,2 BLT plays an important role in preclinical researches, such as gene therapy, cancer research, and drug development.^3,4,5 In the imaging process of BLT, the first step is to build a light transmission model to describe the transmission process of bioluminescence in the organism, and then solve the inverse problem of this light transmission model to reconstruct the distribution of bioluminescent sources using the photon density data of the measured surface near-infrared light of the organism.⁶ However, the influence of light scattering and absorption effects, as well as the limitations of boundary measurements, leads to BLT a severely ill-posed inverse problem, which poses a challenge to the reconstruction methods.^7,8

Many effective strategies have been proposed to alleviate the ill-posedness of BLT. For optimal reconstruction based on the characteristics of BLT, it is necessary to combine valid prior knowledge to uniquely reconstruct the sources.⁹ The most common strategies used to alleviate the ill-posedness of BLT are multispectral strategy, feasible region strategy, etc. For example, the multispectral strategy effectively increases the measurements by spectral decomposition using band-pass filters.^10,11 On the other hand, the feasible region strategy alleviates the ill-posedness by limiting the number of unknowns. Generally, the feasible region can be determined using the surface light distribution, anatomical structure, and other prior information.^12,13 Feng et al.¹² suggested that a rough estimate of the feasible region can reduce the number of unknown variables and greatly reduce the reconstruction problem size. Jiao et al.¹³ proposed the strategy of extracting the feasible region with the highest energy at each projection, which has significant advantages in terms of localization error and computation time.

A wide range of reconstruction methods based on regularization have been proposed and verified in BLT. These methods can be roughly divided into two categories: convex and nonconvex. Typical convex regularization methods include $L$1-norm and $L$2-norm regularization.^14,15,16,17 The nonconvex methods are represented by $L_p$-norm ($0<p<1$) and log-sum regularization.^18,19,20 Studies have indicated that the solution of $L$2-norm tends to produce an over-smoothing effect and result in unfavorable loss of details such as boundary.^21,22 On the contrary, the solution of the $L$1-norm and the nonconvex methods has higher sparsity, but the solving process is usually complicated and time-consuming.^15,20,23

As an alternative sparse signal reconstruction method, greedy algorithms have been introduced into BLT,^24,25 and some other optical molecular tomography, such as look ahead orthogonal matching pursuit (LAOMP),²⁶ sparsity adaptive correntropy matching pursuit (SACMP),²⁷ and compressive sampling matching pursuit (CoSaMP).²⁸ It is shown that the greedy algorithms are very efficient for low sparse signals.²⁹ However, the classical greedy strategies require strict adherence to the restricted isometry property, and the solution is easily caught in the local optimum. Thus, BLT based on the traditional greedy strategies suffer from low morphological recovery and poor multi-target resolution. It is necessary and possible to further optimize the solution by improving the support set selection scheme of greedy algorithms.^26,27,28

In this paper, combining the advantages of the greedy strategies as well as regularization methods, we propose a hybrid reconstruction framework (HRF) for multispectral BLT based on Alpha-divergence to improve the reconstruction results. First, the support set of the greedy strategy is used as a feasible region for reconstructing the source. Based on this, the solution is obtained by weighting the $L$1-norm regularization and $L$2-norm regularization using Alpha-divergence. Finally, the global convergence of the solution and the stability of the framework is achieved by the iterative optimization of the weighted solution and the feasible region. The innovations of HRF are as follows: (1) HRF transforms the support set into a feasible region, avoiding the heavy dependence on the support set selection as in the traditional greedy strategies. (2) The results of the $L$1-norm regularization and $L$2-norm regularization are weighted to avoid over-smoothing and over-sparsity, which optimizes the results and ensures the generality of the framework. (3) In addition, weighting factors are defined and the weighting strategy is decided according to Alpha-divergence to ensure the reasonableness of the weighted results.

This paper is organized as follows. In Sec. 2, the multispectral BLT model and the framework for solving BLT reconstruction are explained. In Sec. 3, experimental settings and the evaluation indicators are demonstrated. Then, in Sec. 4, numerical simulations and in vivo experiments validate the reconstruction performance and preclinical feasibility of the proposed framework. Finally, we present the discussion and conclusion in Sec. 5.

2. Method and Material

2.1. The multispectral BLT model

The diffusion equation (DE) is widely used to describe the light propagation model in BLT studies. In the continuous wave mode, the DE and Robin-type boundary conditions³⁰ can be specified as

where $\mathrm{\Phi }(r)$ describes the source density at the location $r$. $\mathrm{\Omega }$ and $\partial \mathrm{\Omega }$ are the domain and its boundary, respectively. $S(r)$ represents the spatial distribution of the internal sources and $A(r)$ denotes the mismatch factor of the boundary. $v(r)$ is the normal vector pointing out of the surface, and the diffusion coefficient is $D(r)=1∕3({\mu }_a(r)+{\mu }_s^{\prime }(r))$, in which the absorption coefficient and the reduced scattering coefficient are represented by ${\mu }_a(r)$ and ${\mu }_s^{{}^{\prime }}(r)$, respectively.

Equation (1) can be expressed as a Green’s function further. Green’s function relates the fluence rate at the object boundary to the distribution of bioluminescence sources.³¹ In this study, multispectral measurements are used to improve reconstruction accuracy. Given the structural information and optical parameters at different wavelengths, the following linear model of the multispectral BLT can be derived by solving Eq. (1) with the finite element method (FEM)² :

where $b(\lambda )$ is the fluence rate measured at the boundary at wavelength $\lambda$. $\eta (\lambda )$ and $G(\lambda )$ are the weight and the Green’s function for $\lambda$, respectively and $s$ is the source distribution.

Simplifying Eq. (2), we can obtain

where $G$ is a $n\ast m(n\ll m)$ sensitive matrix reflecting the connection between the measured boundary signal $b$ and the inner bioluminescence sources $s$.

2.2. The HRF of BLT reconstruction

For the reconstruction process based on Eq. (3), regularization techniques or greedy strategies are indispensable to reconstruct the unknown source due to the inherent ill-posedness. In this section, we propose an HRF by combining the advantages of greedy strategies as well as regularization methods. In this framework, the support set obtained by the greedy strategy is regarded as the feasible region to reduce the problem size and alleviate the ill-posedness. In this feasible region, the $L$1-norm and $L$2-norm methods are adopted, respectively, to obtain two solutions with different levels of sparsity, smoothness, and morphological information. Then making use of the Alpha-divergence, two solutions are combined to get a weighted solution in the current feasible region. Then, a new feasible region is obtained by calculating the residuals with the above weighted solution. The weighted solution and feasible region are iteratively optimized until the stopping condition is satisfied.

2.2.1. Determine feasible region with greedy strategy

In this research, we follow the greedy strategy of orthogonal matching pursuit (OMP)³² to determine the feasible region. Under a hypothetical assumption that the columns of $G$ are approximated to be orthonormal, the bioluminescence sources $s$ can be exactly recovered from its measurements. Greedy strategy abandons exhaustive search in favor of a series of locally optimal single-term updates. On the $k$th iteration, we can identify these active columns of $G$ as a feasible region by minimizing the residual error $r_k=b-G_{I_k}{s}_{I_k}$. Then, we find the column of $G$ that is maximally correlated between the sensitivity matrix and the residual error, which updates the feasible region, and obtain the temporary solution $s_{\mathrm{\text{temp}}}$ through the pseudo-inverse of $G$. In addition, a greedy stopping condition is needed, i.e., the output is available when the residuals are small enough or when the number of iterations exceeds the maximum. The above process is expressed by the following equation :

where $k$ represents the number of iterations and $m$ represents the number of columns in $G$. $V$ denotes the set of all the column vectors of $G$, and $I_k$ means the feasible region determined in the $k$th iteration, which is composed by adding the most relevant column to the feasible region of the ($k-1$)th iteration. $j\in (V-I_{k-1})$ represents the subscript index not included in the current feasible region. $\mathrm{\text{dot}}(a,b)$ represents the inner product of $a$ and $b$, and $i_k$ indicates the most relevant column. $G_{I_k}^{\perp }$ denotes the pseudo-inverse matrix of $G_{I_k}$, where $G_{I_k}$ is a sensitive matrix corresponding to $I_k$. $s_{\mathrm{\text{temp}}}$ is the temporary solution corresponding to $I_k$ and used as the initial input for subsequent calculations. $\epsilon$ represents the operation precision, which is set to $10^{-5}$.

2.2.2. Obtaining two regularized solutions

Based on the feasible region determined by $I_k$, a new round of BLT reconstruction can be done by solving the following optimization problem :

To deal with the ill-posedness of BLT, proper regularization is necessary to obtain a stable reconstruction. In order to avoid the defects of too-sparse or too-smooth caused by using a single model, we adopt a hybrid model to balance the results by weighting solutions generated from two different regularization model. The two regularized solutions $y_k$ and $z_k$ are as follows :

where ${\lambda }_1$ and ${\lambda }_2$ are the regularization parameters of the $L$1 method and the $L$2 method, respectively. It is important to notice that the problem to be solved in Eq. (6) is a scaled-down inverse problem due to the use of feasible region and the $s_{\mathrm{\text{temp}}}$ obtained by Eq. (4) can be input into the inverse algorithm as the initial solution to accelerate the convergence of the solution. There are plenty of algorithms to solve Eq. (6), we chose LSQR³³ and FISTA³⁴ for solving the $L$2-norm and $L$1-norm regularization problems in our research. And the regularization parameters are empirically chosen, which are adjusted in the range of $10^{-3}$–$10^{-1}$, respectively.

2.2.3. Alpha-divergence and solution update

To better combine the reconstruction results, the Alpha-divergence³⁵ is used to measure the consistency between two solutions, which is expressed as

where $\alpha$ is a real parameter, which is set to 0.5 by experience, and $n$ is the number of rows in $G$. It can be seen that the Alpha-divergence is a convex function. Consequently, if the distributions match exactly, it is zero. In other cases, it is always positive. To ensure the direction of convergence, the weight $w_k=1∕{(y_k+z_k)}^2$ is introduced. Therefore, the weighted solution $s_k$ in the current feasible region is determined as follows : where $\mathrm{\text{tol}}$ means divergence threshold which is selected to update the weighted solution, and for the experiment in this paper, the value of tol was taken in the range of 0.2–0.6.

The detailed steps of the HRF based BLT are described in Table 1.

3. Experiment Setting

3.1. Numerical simulations

Numerical simulations were designed to verify the feasibility and performance of the proposed method, which include single-source reconstruction and dual-source resolution. To provide anatomical information of the biological tissues, a 34mm height torso was extracted from a 3D mouse atlas in the numerical simulations.³⁶ The optical parameters for different organs at 610nm, 630nm, and 650nm wavelengths are shown in Table 2.³⁷

For single-source reconstructions, a spherical source with different radii of 0.8mm, 1mm, and 1.5mm was placed in the liver. The true source center was fixed at (20.1, 7, 16.8) mm. To generate simulated measurements, the torso part including the light source was discretized into a tetrahedral mesh, and the DE combined with Robin-type boundary conditions was solved using FEM. The forward meshes for three groups of single-source simulations included about 18,000 nodes and 100,900 elements.

To verify the multiple-targets resolving capability of HRF, three groups of simulations for dual-source with different edge-to-edge distances (EEDs) were conducted. The sources in this section were spherical sources with a radius of 1mm and the EEDs included 3mm, 4mm, and 5mm, and the corresponding coordinates of the dual-source were (16, 7.5, 17) mm and (21, 7.5, 17) mm, (15, 7.5, 17) mm and (21, 7.5, 17) mm, (14, 7.5, 17) mm and (21, 7.5, 17) mm, respectively. The forward mesh for double-source simulation contained about 20,000 nodes and 111,500 elements.

In all simulation experiments, the true intensity of light sources was set to 1nW/mm³, and 5$\%$ random Gaussian noise was added into the generated measurements to mimic the inevitable noise and error in data acquisition. The mesh used for reconstruction was coarser than the forward meshes, containing 9249 nodes and 50,548 elements. And the corresponding weights for measurements at different wavelengths were 0.29, 0.48 and 0.23, respectively.³⁸

3.2. In vivo experiments

To further assess the performance of HRF, we further performed an $in$$vivo$ experiment. The animal study was reviewed and approved by the Animal Ethics Committee of the Northwest University of China. A 6-week-old female BALB/c nude mouse was used for the abdominal source implantation experiment. In the experiment, the mouse was anesthetized with 3$\%$ isoflurane air mixture during the surgical and imaging procedures. A catheter approximately 5mm in length and 1.5mm in diameter was implanted in the abdomen of the mouse. It was filled with about 10$\mu$L luminescent solution to serve as the light source for the subsequent reconstruction.

The experimental data was acquired by a homemade Micro-CT/BLT system.³⁹ To acquire multispectral surface measurements at 610nm, 630nm, and 650nm, three band-pass filters were used and the EMCCD camera (iXon Ultra, Andor, Northern Ireland, UK) was operated at $-80^{\circ }$C to reduce the dark current and thermal noise. Only single-view BLI was acquired for each wavelength. To obtain the anatomical structure, the X-ray flat panel detector (C7942CA-22, Hamamatsu Photonics, Hamamatsu Japan) was used for CT imaging with an X-ray source voltage of 90kV and power of 27W to obtain $360^{\circ }$ projection data, and the Feldkamp–Davis–Kress (FDK) algorithm was used for CT reconstruction.

Before BLT reconstruction, some pre-processing operations were required. The CT data of the mouse model were discretized into a tetrahedral mesh containing 9591 nodes and 48,111 elements, which is used for the subsequent 2D–3D surface energy mapping and BLT reconstruction. Figure 1(a) shows the bioluminescence image (at 650nm) overlapped on the white light images. The 2D optical image was registered to 3D surface the obtain the measurement of surface light at each wavelength according to the pinhole mapping model. Figure 1(b) shows part of the tetrahedral mesh of the mouse model with the implanted source. Because the implanted catheter can be captured by CT, the position of the true source center can be determined according to the CT coordinates, which was at (23.6, 22.4, 20) mm.

3.3. Framework comparison and evaluation index

To investigate the potential and feasibility of HRF, LSQR as an $L$2-norm regularization method, FISTA as an $L$1-norm regularization method, and OMP as a greedy method are chosen as comparison methods. To quantify the reconstruction performance, the location error (LE), the Dice coefficient (DICE), and the time cost (TIME) are used to evaluate the methods.

$\mathrm{\text{LE}}$ is used to measure the distance variation between the reconstructed center and the ground truth center, whose definition is given by

where $(x,y,z)$ and $(x_0,y_0,z_0)$ are the center coordinates in the reconstructed and true source regions, respectively.

$\mathrm{\text{DICE}}$ describes the shape similarity between two objects :

where $X$ and $Y$ are the reconstruction and true volume of the target area.

4. Experiment Results

4.1. Reconstruction of single-source

Figure 2 shows three-dimensional views of the single-source reconstruction results and transverse views at $z=16.8$mm. The results of OMP show a large localization deviation, and that of LSQR have significant artifacts and large position bias, which indicates the over-smoothing effect consistent with the drawback of the $L$2-norm method. Table 3 presents the average $\mathrm{\text{LE}}$, $\mathrm{\text{DICE}}$, and $\mathrm{\text{TIME}}$ of 10 independent runs under 5$\%$ noise. The quantitative results show that the HRF has the best performance compared to the other three methods, except for the time cost. Although the computation time of HRF is longer than that of OMP and LSQR, it is shorter than that of FISTA.

4.2. Reconstruction of dual-source

The reconstructed images and quantitative results of dual-source are given in Fig. 3 and Table 4, respectively. Among these methods, the performance of OMP is the worst. The reconstructed sources have a significant deviation from their real position at EED of 3mm and 4mm. The reconstructed images of LSQR and FISTA also have some artifacts. From these results, it can be seen that the HRF can clearly distinguish two sources even with 3mm EED. HRF has the best performance in terms of both localization accuracy and shape recovery, although the reconstruction time is not optimal.

Table 4 summarizes the quantitative results of localization error, and shape similarity. As shown in Table 4, from a series of comparison results, the HRF has the largest DICE and the smallest LE. In addition, for the EED of 3mm, as shown in Table 4, the average LE and DICE of HRF reach 0.30mm and 0.62mm, respectively, which are satisfactory results in terms of localization accuracy, and shape recovery ability. For the EEDs of 3mm and 4mm, the reconstructed result of OMP is hardly distinguishable and completely deviates from the true source, resulting in the worst results for DICE and LE in Table 4.

4.3. In vivo experiment

The quantitative results of in vivo experiment are shown in Table 5. Figure 4 shows the 3D views, sagittal ($y=22.4$mm), coronal ($x=23.6$mm), and transverse ($z=20$mm) views, respectively. From the 3D view and 2D view, it can be seen that the HRF is more accurate in source position reconstruction, which is also proved by the LE in Table 5. As can be seen from Table 5, the results of in vivo and simulation experiments are consistent, and the DICE of the HRF is much higher than that of OMP, LSQR, and FISTA, reaching 75$\%$. This is because the reconstructed area of the HRF is more consistent with the actual area as seen from the 3D view, and the reconstructed part almost completely covers the source. The reconstructed positions of the LSQR and FISTA deviate from the actual position and are tilted toward the animal surface, with LEs of 0.89mm and 0.63mm, respectively.

5. Discussion and Conclusion

In this paper, to improve the accuracy of BLT reconstruction, an HRF for model-based multispectral BLT based on Alpha-divergence is proposed. It combines the advantages of the greedy strategy and regularization methods. It embeds a hybrid optimization process into the iteration steps of OMP by taking the support set as a feasible region and using a hybrid approach to optimize the solutions derived by two different methods under the feasible region. With the alternative update of the feasible region and the weighted solution, the HRF overcomes the limitations of using a single regulation method or greedy strategy. It weighs the results of $L$1-norm and $L$2-norm methods to avoid over-smoothing and over-sparsity. Moreover, the weighting factors are defined and the weighting strategy is determined according to the Alpha divergence, which ensures the reasonableness of the weighted solutions.

The results of single-source simulations showed that the HRF performed best in LE and DICE, which reflected the advantages of the HRF in localization accuracy and morphological recovery. For the dual-source experiment with EEDs of 3–5mm, HRF maintained the best performance, which indicated the advantage of HRF in improving spatial resolution. The preliminary result of $in$$vivo$ experiment shows that the reconstructed area of HRF could better approximate the actual target distribution. The slightly higher time cost of HRF is mainly because two regularization algorithms are embedded in the solution process of OMP to obtain better reconstruction. There is no doubt that the inner iteration increases the complexity of the reconstruction. However, the greedy strategy in HRF provides a priori information about the feasible region and the initial solutions for the subsequent regularization-based reconstruction, so the size of the inverse problem is reduced and the inner iteration process is accelerated. In general, although the computation time of HRF is longer than that of OMP and LSQR, it is acceptable and shorter than that of FISTA. This demonstrated that HRF had great potential for practical applications.

Basically HRF belongs to model-based reconstruction framework, and the quality of the reconstructed images is inevitably affected by the forward model and the FEM mesh used for reconstruction. For the mesh, if it is too sparse, the location accuracy and the shape of the reconstructed light source will be affected, and if the mesh is too dense, it will increase the ill-posedness and computational burden. In addition, the selection of parameters is an important influencing factor. The parameters in HRF can be divided into the following categories, (i) tol, (ii) regularization parameters, (iii) stopping conditions for the iterative steps. Among them, the key parameters are tol and the regularization parameters. The tol is a threshold value that relates to the Alpha-divergence of two solutions. Generally, the more similar the solutions obtained by the two methods, the smaller the tol value will be. In this implementation the tol is determined empirically, ranging from 0.2 to 0.6. As is well known, the choice of regularization parameter is important for an ill-posed inverse problem. Currently, the regularization parameters used in this work are manually adjusted. The adaptive method for parameter choice based on the generalized cross-validation criterion or using the L-curve method can potentially improve the results.

In fact, the hybrid framework is quite flexible, the two regularization methods used in HRF are not limited to $L$1-norm regularization methods or $L$2-norm regularization methods. In addition, the HRF is general and extensible to the similar imaging modalities, such as FMT and XLCT, which is worthy of further investigation.

In conclusion, we propose a new HRF to reduce the ill-posedness of BLT and improve the reconstruction quality. Compared with the comparison methods, HRF has the advantages of accurate localization, excellent robustness, and high practicality. It is believed that this new framework will further facilitate various preclinical applications of BLT and promote the development of optical molecular tomography in theoretical research.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This study was funded by the National Natural Science Foundation of China under Grants Nos. 11871321, 61901374, 61906154, and 61971350 and Postdoctoral Innovative Talents Support Program under Grants No. BX20180254.