Non-convex sparse regularization approach framework for high multiple-source resolution in Cerenkov luminescence tomography

Hongbo Guo; Zhenhua Hu; Xiaowei He; Xiaojun Zhang; Muhan Liu; Zeyu Zhang; Xiaojing Shi; Sheng Zheng; Jie Tian

doi:10.1364/OE.25.028068

1. Introduction

Cerenkov radiation is generated when a high-energy charged particle moves uniformly in the substance at a velocity exceeding the speed of light in that medium [1, 2]. Cerenkov luminescence imaging (CLI) technique that harnessed Cerenkov radiation to image radionuclides with optical imaging instruments, has received significant interests in preclinical research and clinical medical imaging [3, 4]. In preclinical applications, CLI has been used for lymph node mapping [5, 6], Cerenkov luminescence endoscopy [7–9] and tumour resection [10, 11]. In clinical applications, CLI is mainly applied for human thyroid imaging [12], endoscopic examination of rectal cancer [13], real-time imaging of human breast cancer [14], and determining tumour resection margins and lymph node status [15]. The main advantages of CLI include the ability to easily overlay anatomical images, widely available U.S. food and drug administration (FDA) approved molecular probes for clinic use, simple and relatively inexpensive requirements of imaging instrument, and high surface resolution [7, 16, 17]. Cerenkov luminescence tomography (CLT) as one of the most attractive applications of CLI has been explored for obtaining 3D spatial distribution of radiopharmaceutical in living animals [18, 19]. Different from 2D imaging, CLT provides the depth and semiquantitation information of radioactive probes [20] in biomedical applications. Thus, CLT has been considered as a highly promising technique for preclinical diagnosis and clinical translation [21, 22].

With the help of the clinical application of CLI in tumour and lymph node imaging [5, 6, 10, 11, 15], CLT has the potential to be used for cancer staging. In clinic research, cancer staging generally takes into account the tumor size, the degree of invasion and the number of metastases for evaluating the clinical progression of tumor [23]. The most clinically useful staging system is the tumor, node, and metastasis (TNM) staging system [24]. Generally, the radiopharmaceutical can be intaken by the tumor cells. If malignant lesions are present, optical molecular imaging will have a difference with the healthy tissue. Moreover, if using optical molecular technique for cancer staging, tumour and lymph node can be regarded as multiple-source. And the correlation of location (i.e. the spatial coordinates), intensity (i.e. the dose of the uptake radionuclides), and size (i.e. the lesions area) can be used to assess the state of the cancer. Multiple-source resolution is the process of determining the details of source, such as location, intensity, size etc. However, accurately reconstructing multiple-source is a highly challenge issue. On one hand, the multiple-source resolution in CLT is an ill-posed problem because of light scattering and absorption, and highly surface-weighted of CLI signal [25–31]. On the other hand, the interaction effect of multiple-source, very complex surface-weighted CLI signal and the diversity of tumor biological characteristics make multiple-source resolution more complicated. Therefore, how to solve the multiple-source resolution problem and enhance the reconstruction quality and stability is the main focus of this study.

In previous study of CLT, priori knowledge strategy as one typical technique has been widely used to enhance the location accuracy and quantitative precision. Zhong et al. incorporate the CT anatomic structural information into the CLT inverse problem to improve the quality of source reconstruction [32]. Furthermore, Hu et al. have located the Na¹³¹I radioactive source accurately by using the single photon emission computed tomography (SPECT) functional knowledge [33]. Adopting some priori knowledge, the ill-posedness of inverse problem can be released and the stability of CLT reconstruction has made great improvement. Furthermore, some priori knowledge can be used as a standard for evaluating the reconstruction results. Thus it’s necessary to adopt some priori knowledge for reconstructing multiple-source accurately.

In addition to priori knowledge, many effective reconstruction algorithms have been put forward. In literatures at present, many inverse effective reconstruction algorithms have been testified to have a good performance in single source localization, including newton method with an active set strategy [34], Tikhonov regularization [18, 32], orthogonal matching pursuit (OMP) [35], region growing method [36], pre-conditioned conjugate gradient method [26], non-negative least square (NNLS) algorithm [37], incomplete variables truncated conjugate gradient (IVTCG) [27, 38]. These algorithms are all based on ℓ₂-norm regularization or ℓ₁-norm regularization and belong to the category of convex sparse regularization approach (CSRA) framework. In other optical imaging techniques, some effective algorithms based on non-convex sparse regularization approach (nCSRA) framework have been developed and tested [39, 40]. In mathematics, nCSRA framework is commonly proposed to solve the ℓ_p regularization algorithms with (0 < p < 1). Compared with CSRA, the sparseness and stability of nCSRA have a great improvement [41]. And also in reconstruction quality, algorithms based on nCSRA framework have a better performance in restricting reconstruction artifacts comparing to CSRA framework. Therefore nCSRA framework can provide a promising strategy to solve the multiple-source resolution problem in CLT.

In this work, the nCSRA framework was proposed to solve the multiple-source resolving problem in CLT. In nCSRA framework, the non-convex problem was converted into an iterative, weighted and hybrid ℓ₁-norm regularization problem. And two typical algorithms homotopy and iterative shrinkage-thresholding algorithm (ISTA) were exploited to compare the performance of nCSRA framework and CSRA framework. The major contributions of this work are threefold: (1) nCSRA framework with high sparseness and stability was presented specially for the multiple-source resolution problem. (2) in order to test the performance of nCSRA framework, groups of comparison experiments has been designed in numerical simulations and in vivo imaging experiments. It has experimentally shown that the proposed nCSRA framework performance better than CSRA framework in spatial resolution, intensity resolution, and size resolution. (3) an integrated optical/CT system and a high resolution microPET scanner have been adopted to obtain anatomic structural and functional imaging information for improved CLT reconstruction and results evaluation. The rest of the paper is organized as following: Section 2 provides the nCSRA framework, imaging system and some relative experimental designs. Section 3 presents the results of numerical simulations and in vivo imaging experiments. Finally, some discussions and conclusions have been made for this paper in Section 4.

2. Materials and methods

2.1. Photon propagation model

The surface detectable CLI spectrum is in the visible light range that scattering effects are dominant over absorption effects. Thus, diffusion approximation to radiative transfer equation is widely used to describe light propagation in CLT. The diffusion approximation combined with the Robin-type boundary condition can be expressed as [42]:

{\begin{array}{l} - \nabla \cdot (D (r)) + μ_{a} (r) Φ (r) = S (r), (r \in Ω) \\ Φ (r) + 2 A (r) D (r) (v (r) \cdot D (r)) = 0, (r \in \partial Ω) \end{array}

where Φ(r) denotes the light fluence rate at position r ∈ Ω. S(r) is the Cerenkov luminescence source. D(r)=1/3(µ_a(r) + (1 − g) µ_s(r)) is the optical diffusion coefficient with µ_a(r) and µ_s(r) being the optical absorption coefficient and scattering coefficient, and g is the anisotropy coefficient. A(r) is the boundary mismatch factor accounting for different refractive indices across the boundary ∂Ω, ν(r) denotes the unit outer normal. In the finite element method (FEM) framework, a linear relationship between the unknown Cerenkov luminescence source and the detected surface partial luminescent flux fluence rate can be obtained as [43]:

M X = Φ

where M is the symmetric matrix, X is the unknown Cerenkov luminescence source distribution, and Φ is the surface detected photon flux.

2.2. Non-convex sparse regularization algorithm (nCSRA) framework

In the minimization scheme with nCSRA, the CLT reconstruction can be formulated as:

\min_{X} \frac{1}{2} {‖ M X - Φ ‖}_{2}^{2} + ζ {‖ X ‖}_{p}^{p}, 0 < P < 1

where ‖X‖_p represents the ℓ_p quasi-norm, defended by

{‖ X ‖}_{p} = {(\sum_{i = 1}^{n} {| X_{i} |}^{p})}^{1 / p}

; ζ is the regularization parameter, which is used to balance the relative weighing between the objective function term

{‖ M X - Φ ‖}_{2}^{2}

and the penalty term

{‖ X ‖}_{p}^{p}

. Because the ℓ_p-norm regularization is a nonconvex and no-smooth optimization. An intermediate operation was conducted as:

{‖ X ‖}_{p}^{p} = \sum_{i = 1}^{n} {| X_{i} |}^{p} = \sum_{i = 1}^{n} \frac{| X_{i} |}{{| X_{i} |}^{1 - p}}

At t + 1 iteration, if |X_i|^1−p in Eq. (4) was approximate to

{| X_{i}^{t} |}^{1 - p}

. Eq. (3) can be converted into the following form of an iteration solution [41]:

X^{t + 1} = \min_{X} \frac{1}{2} {‖ M X - Φ ‖}_{2}^{2} + ζ \sum_{i = 1}^{n} \frac{| X_{i} |}{{| X_{i}^{t} |}^{1 - p}}

Eq. (5) can be equated to a weighted ℓ₁-norm regularization:

X^{t + 1} = \min_{X} \frac{1}{2} {‖ M X - Φ ‖}_{2}^{2} + \sum_{i = 1}^{n} λ_{i} | X_{i} |

where

λ_{i} = \frac{ζ_{i}}{{| X_{i}^{t} |}^{1 - p}}

, ζ_i is the ith element of regularization parameter ζ.

Eq. (6) can be solve by the some common used algorithm based on ℓ₁-norm regularization (homotopy [44] and iterative shrinkage-thresholding algorithm (ISTA) [45] in this study). And in practice, a fix positive (k = 1) was introduced in $λ_{i} = \frac{ζ_{i}}{{| X_{i}^{t} |}^{1 - p} + κ}$ to avoid any component of X^t being zero.

2.2.1. Homotopy

The homotopy formulation for Eq. (6) can be expressed as:

{\begin{array}{l} \min_{x} {‖ λ X ‖}_{1} + \frac{1}{2} {‖ M X - Φ ‖}_{2}^{2} + (1 - ε) ρ^{T} X \\ ρ = - λ s i g n (\hat{X}) - M^{T} (M \hat{X} - Φ) \end{array}

where ε 2 ∈ [0, 1] is the Homotopy parameter,

\hat{X}

is the given warm-start vector at ε = 0. In every iterations, the update direction:

\partial X = {\begin{matrix} {(M_{Γ}^{T} M_{Γ})}^{- 1} ρ_{Γ} & o n Γ \\ 0 & o t h e r w i s e \end{matrix}

Where Γ is the the support of

\hat{X}

. For the previous solution X^∗, when increasing ε by a small value σ, the step size of update δ can be computed by: δ = min (δ⁺, δ⁻). In which

{\begin{array}{l} δ^{+} = \min_{i \in Γ^{c}} (\frac{λ_{i} - h_{i}}{d_{i}}, \frac{- λ_{i} - h_{i}}{d_{i}}) \\ h_{i} = M_{i}^{T} (M X^{*} - Φ) + (1 - ε) ρ_{i} & , i \in Γ^{C} \\ d_{i} = σ (M_{i}^{T} M \partial X - ρ_{i}) \end{array}

δ^{-} = \min_{i \in Γ} {(\frac{- X^{*}}{\partial X^{*}})}_{+}

So the solution can be updated by X^∗ = X^∗ + δ ∗ ∂X.

2.2.2. Iterative shrinkage-thresholding algorithm (ISTA)

ISTA is based on the shrinkage function: shrink (a, z) = max(a – z, 0) * sign(z). With a sufficiently small step size ξ, the analytical solution of Eq. (6) can be derived as:

{\begin{array}{l} X_{i}^{t + 1} = shrink ((X^{t} + 2 ξ M^{T} (Φ - M X {))}_{i}, ξ λ_{i}) \\ ξ < 1 / {‖ M^{T} M ‖}_{2} \end{array}

2.3. Numerical simulations

A series of numerical simulations based on a digital mouse model were conducted to evaluate the performance of nCSRA framework. The torso section of the model with a height of 25mm was selected and divided into: muscle, heart, liver, lung, kidney and bone, as shown in Fig. 1(a). The empirical optical parameters for different organs were listed in Table 1 [35]. In numerical simulations, the FEM was employed to simulate the surface flux distribution. To mimic measurement uncertainty and inevitable noise, 10% additive Gaussian noise was added to the surface simulated measurements. For inverse reconstruction, the digital mouse model was discretized into a uniform tetrahedral mesh, including 14,465 nodes and 74,903 tetrahedral elements, as shown in Fig. 1(b).

Fig. 1 Numerical simulation settings. (a) The torso section of digital mouse model; (b) The uniform tetrahedral mesh used in the inverse reconstruction problem; (c–g) The real 3D distribution of Cerenkov sources in numerical simulations.

Download Full Size | PDF

Table 1. Optical properties for numerical simulation and in vivo experiments.

View Table | View all tables in this article

In order to simulate the multiple-source resolution problem more practical, including spatial resolution, intensity resolution and size resolution, three groups of simulations were conducted. And spatial resolution, intensity resolution and size resolution were be used as a indicator for testing the performance of nCSRA and CSRA frameworks. Spatial resolution referred to the minimum source separation which can be distinguished by reconstruction algorithm; Intensity resolution was the ability to overcome the interaction effect for accurately discriminating between source intensities; Size resolution reflected the ability to overcoming the interaction effect for accurately reconstructing source volume; Detail source parameters were listed in Table 2.

Table 2. Cerenkov source parameters in numerical simulations.

View Table | View all tables in this article

2.3.1. Spatial resolution

Two spherical sources with same volume and equal intensity were placed in the mouse abdomen. Three center-to-center separations of 4.8, 3.6 and 2.4mm were designed to test the spatial resolution, as shown in Figs. 1(c)–1(e). And the reconstructed source center, location error (LE) were listed to evaluate the performance of the reconstruction results.

2.3.2. Intensity resolution

Two spherical sources with equal center-to-center separation and volume. While the intensity ratios of two sources were set to 1: 2, 1: 4, and 1: 8. In this simulations, reconstructed source intensity, and relative intensity error (RIE) were listed to evaluate the source intensity accuracy of the reconstruction results.

2.3.3. Size resolution

Two spherical sources with equal intensity and center-to-center separation. The volume of No. 1 source was fixed to 3.09mm³, while the volume of No. 2 source was set to 6.03mm³ and 1.30mm³, respectively, as shown in Figs. 1(f)–1(g). And the reconstructed source volume, and the relative volume errors (RVE) were used to evaluate the ability of size resolution.

2.4. In vivo imaging experiments

In vivo imaging experiments were performed to assess the in vivo resolving power for multiple sources. Female BALB/C nude mice (4–6 weeks old, n = 4) were used for establishment of the subcutaneous tumour mouse model. For every mouse model, 4T1 breast tumor cells were subcutaneously injected into the back of the mice, and more detail information of the cell injection as listed in table 3. All animal experimental procedures were carried out strictly in accordance with the guidelines of the Institutional Animal Care and Use Committee of Chinese PLA General Hospital. All animal procedures were performed under isoflurane gas anaesthesia (3% isoflurane-air mixture) to minimize suffering to mice. Five days later, every mouse was injected with about 200 ± 50µCi fluorine-18-fluorodeoxyglucose (¹⁸F-FDG) through the tail vein, and optical and CT volume data were obtain by using an existing integrated optical/CT system [27]. Then, every mouse was subjected to a PET scan of 5 min with the GENISYS4 microPET scanner. Lastly, these 4T1 subcutaneous tumours were resected for measuring tumor size with vernier caliper. For evaluating the reconstruction results, the center coordinate of actual tumor was determined by CT images. Moreover the PET images were used to conduct a semiquantitation study about the actual CLI intensity in tumor.

Table 3. Detailed information on cell injection in tumor model-establishment process.

View Table | View all tables in this article

In the CT volume acquisition process, after some needed preparations (tube warming up, X-ray calibration, and CT attenuation-corrected), the tube voltage and tube current were set as 40kV p and 300mA. Then the rotating stage was rotating 360 degrees with 1° intervals to capture the X-ray projection images. Then the projection data were reconstructed by Feldkamp-Davis-Kress (FDK) reconstruction algorithm. In optical acquisition process, the CLI and white light data were collected by a highly sensitive TE-cooled backilluminated EMCCD camera. Before optical acquisition, the EMCCD camera was cooled to −80°C to reduce the effects of thermal noise. In the CLI acquisition process, long exposure time (5min), high gain value (300), high shift speed (12.9µsecs), low speed readout rate (1MHz at 16-bit), and 4 × 4 binning were set up to increase signal-to-background ratio. While for white light, short exposure time (1s), no gain value, and 1 × 1 binning were set up to improve the image resolution.

2.5. Statistical analysis

Statistical comparisons were made using Student’s t-test and GraphPad Prism 5 software. P values 0.05 were considered to indicate significance. Box plots and linear fitting results were made using Origin software. Statistical data were expressed as mean ± SD (Standard Deviation).

3. Results

3.1. Numerical simulation

3.1.1. Spatial resolution

The reconstructed results with different center-to-center separations are shown in Fig. 2. Reconstructed Cerenkov signal was delineated with red in 3D view, meanwhile in transverse view the white circles indicated the actual positions of the sources in the slice over the centers of the sources at Z = 7mm. From these results, we can found that Cerenkov sources with center-to-center separation 2.4mm could be distinguished clearly by algorithms based on nCSRA framework. However for CSRA framework, the reconstructed sources were blurred with some image artifacts even when source separation distance was larger than 3.6mm, as shown in Figs. 2(a) and 2(b). What’s worse was that algorithms based on CSRA framework could not distinguish two sources when source separation was 2.4mm, as shown in Figs. 2(c). For better evaluating the reconstruction results, some quantitative indicators (reconstructed source center and location error (LE)) are listed in Table 4 and shown in the box plot, as shown in Fig. 3. In the box plot, red dots and blue square dots stood for the location error of each Cerenkov sources and the average location error, respectively. The statistical results demonstrated that the performance of homotopy based on nCSRA framework was the best with LE 0.46 ± 0.22mm, next was ISTA based on nCSRA framework with LE 0.57 ± 0.39mm and homotopy based on CSRA framework with LE 0.63 ± 0.43mm, and the worst was ISTA based on CSRA framework with LE 0.78 ± 0.51mm. The statistical results also indicated that the differences of these two framework was significant (p = 0.0481) with center-to-center separation less than or equal to 3.6mm; These quantitative results demonstrated that nCSRA framework was able to improve the spatial resolution capability of algorithms, compared with CSRA framework.

Fig. 2 Results of the spatial resolution experiments.(a)–(c) The reconstructed Cerenkov sources with different center-to-center separations in 3D view and transverse view. In 3D view, the red areas denote the reconstructed Cerenkov sources. In transverse view, the white circles represent the actual sources in the slice over the centers of the sources at Z = 7mm.

Download Full Size | PDF

Table 4. Quantitative results of the spatial resolution experiments.

View Table | View all tables in this article

Fig. 3 Box plot about the comparison of location error. In the 3D view, the red regions denote the reconstructed sources. Red dots and blue square dot mean the location error of every sources and the average location error, respectively.

Download Full Size | PDF

3.1.2. Intensity resolution

The reconstruction results with different intensity ratios are displayed in Fig. 4. Based on nCSRA framework, the reconstructed intensities exhibited a good consistency with the real intensity ratio. However, the consistency of the reconstruction intensity based on CSRA framework had a great influence by interaction effect of sources. When intensity ratio was 1: 4 and 1: 8, the reconstructed intensity ratio of two sources hardly changed. Some quantitative indicators (reconstructed intensity and relative intensity error (RIE)) were also listed in Table 5. It’s obviously that the reconstructed source intensities obtained by algorithms based on nCSRA framework were very close to the true value with smaller RIE. But, the RIEs of CSRA framework were progressively increasing with enlarging the difference of sources intensity. More accurately, the performance of homotopy based on nCSRA framework was the best with RIE 11.12 ± 5.20%, next was ISTA based on nCSRA framework with RIE 15.91 ± 11.47% and homotopy based on CSRA framework with RIE 34.48 ± 25.43%, and the worst was ISTA based on CSRA framework with RIE 55.02 ± 17.68%. From Fig. 5, fitted curves of algorithms based on nCSRA framework were much closer to the line of Y = X with higher R² than that based on CSRA framework. Linear fitted results also indicated that the reconstructed source intensity of CSRA framework was much smaller than the actual value due to the limitation of the sparseness. In general, these intensity resolution simulations results indicated that nCSRA framework enhanced the intensity resolution capability compared with the CSRA framework.

Fig. 4 Results of the intensity resolution experiments. (a)–(c) Reconstructed sources by algorithms based on nCSRA and CSRA frameworks with different intensity ratios in 3D view and transverse view. In 3D view, the red areas denote the reconstructed Cerenkov sources. In transverse view, the white circles represent the actual sources in the slice over the centers of the sources at Z = 7mm.

Download Full Size | PDF

Table 5. Quantitative results of the intensity resolution simulations.

View Table | View all tables in this article

Fig. 5 Reconstructed intensities of every source and corresponding linear fitted curves. The black square, red circle, blue regular triangle, pink inverted triangle and green rhombus indicate the reconstructed intensities by algorithms based on nCSRA and CSRA frameworks and actual source intensities, respectively; linear fix of intensity values were expressed by indicate color lines and linear functions.

Download Full Size | PDF

3.1.3. Size resolution

Fig. 6 showed the reconstructed Cerenkov signals based on nCSRA and CSRA framework in 3D view and transverse view with different size ratios. It should be noted that the reconstructed source size of CSRA framework were much large than that based on nCSRA framework, epically for the results of ISTA. Compared with the actual size, nCSRA had a better performance in case of small source size. While for CSRA framework, the reconstructed source size was much larger than the actual size. The reconstructed source volumes and relative volume error (RVE) were listed in Table 6. It was easy to find that RVEs obtained by algorithms based on nCSRA framework were much smaller than that based on CSRA framework. More accurately, the performance of homotopy based on nCSRA framework was the best with RVE 14.85 ± 8.00%, next was ISTA based on nCSRA framework with RVE 21.48 ± 12.38% and homotopy based on CSRA framework with RVE 97.48 ± 90.44%, and the worst was ISTA based on CSRA framework with RVE 263.95 ± 190.93%. The linear fitted curves of reconstructed source volume and true volume were plotted as shown in Fig. 7. Compared with the CSRA framework, the linear representations of algorithms based on nCSRA framework were much more close to the line of Y = X with higher R². Thus, nCSRA framework could obtain a reasonably source volume and had a great improvement in size resolution.

Fig. 6 Reconstruction results in size resolution experiments. (a)–(b) Reconstructed sources by algorithms based on nCSRA and CSRA frameworks with different size ratios in the 3D view and the transverse view. In 3D view, the red areas denote the reconstructed Cerenkov sources. In transverse view, the white circles represent the actual sources in the slice over the centers of the sources at Z = 7mm

Download Full Size | PDF

Table 6. Quantitative results of the size resolution experiments.

View Table | View all tables in this article

Fig. 7 Reconstructed source volumes and corresponding linear fitted curves. The black square, red circle, blue regular triangle, pink inverted triangle and green rhombus indicate the reconstructed intensities by algorithms based on nCSRA and CSRA frameworks and actual source volumes, respectively; linear fix of intensity values were expressed by indicate color lines and linear functions.

Download Full Size | PDF

3.2. In vivo imaging experiments

For in vivo imaging of 4T1 breast tumor mouse models, Fig. 8(a) showed the results of the white-light image overlaid with Cerenkov luminescence image. All tumor images were consistent with the previous injection location and dose. Beside normal high glucose metabolism organs and imaging prob clearance tissue (heart, bladder, liver, and brain etc.), actual tumor areas could be clearly distinguished as yellow ellipses delineated in PET images, as shown in Fig. 8(b). The reconstructed Cerenkov luminescence signal fused with CT images had good consistency with the real ¹⁸F-FDG distribution obtained by the commercial small animal PET device. The comparison results showed that: (1) compare with CSRA framework, algorithms based on nCSRA framework was able to restrict reconstructed artifact in some degree; (2) in terms of reconstructed tumor size, the volume reconstructed by algorithms based on nCSRA framework was much smaller than that based on CSRA framework. (3) for some tumors with small separation, algorithms based on CSRA framework could not distinguish them. Fortunately, nCSRA framework had a good performance in spatial resolution. For every mouse model, the photograph of corresponding resected tumors are shown in Fig. 8(c). For quantitative comparison, the box plot listed the location error in Fig. 8(d). Compared with the CSRA framework, nCSRA had a great improvement in average location error (less than 1.0mm), low discrete degree and significant differences (p = 0.009). In order to make a semiquantitation study, the recovery Cerenkov luminescence intensity as the function of the PET pixel value was shown in Fig. 8(e). By comparing R², the reconstructed intensities obtained by algorithms based on nCSRA framework had a great consistency with the real ¹⁸F-FDG intake in tumor area. After measuring resected tumors size with a vernier caliper, reconstructed tumor volume and the actual tumor volume were compared and plotted as shown in Fig. 8(f). The linear fitting results indicated that the size resolution capability of nCSRA framework outperformed than CSRA framework. And the performance of ISTA based on nCSRA framework was best, secondly was homotopy based on nCSRA framework. And the reconstructed volumes of ISTA based on CSRA framework were much larger than the true values. From these comparison results in in vivo experiments, it was further evident that algorithms based on nCSRA framework had a better performance than that based on CSRA framework in spatial resolution, intensity resolution, and size resolution.

Fig. 8 The results of in vivo imaging experiments. (a) White-light image overlaid with Cerenkov luminescence image for every mouse model. (b) Reconstructed Cerenkov luminescence signal fused with CT imaging by nCSRA and CSRA frameworks and real ¹⁸F-FDG distribution; (c) Photograph of corresponding resected tumors; (d) Box plot of the comparison results in location error; (e) Semiquantitation comparison results of reconstructed intensities and corresponding linear fitted curves; (f) Comparison results in reconstructed source volumes and corresponding linear fitted curves.

Download Full Size | PDF

4. Discussion and conclusion

As a promising preclinical imaging technique, CLI has been paid much attention in imaging theories, acquisition equipments, and biomedical applications. This work mainly focused on the inverse reconstruction algorithm for better recovering the 3D distributions of radiopharmaceutical. More specifically, we proposed a nCSRA framework to solve the multiple-source resolving problem in CLT. Homotopy and iterative shrinkage-thresholding algorithm was adopted to test the performance based on nCSRA framework and CSRA framework. The comparison results in numerical simulations and in vivo experiments demonstrated that algorithms based on the proposed nCSRA framework had a better multiple-source resolution capability in terms of position resolution, intensity resolution and size resolution.

In algorithms design aspect, the non-convex regularization problem was converted into a series of weighted ℓ₁ regularization problem and solved by homotopy and ISTA algorithm. In this study, nonconvex ℓ_p methods with p near 0.5 performed the best. If a higher value of p was used, the solution of nCSRA framework was more close to ℓ₁ norm regularization. Homotopy and ISTA algorithm fall under two different categories of mathematical theory, i.e. the greedy and approximate iteration strategies. Homotopy updates by means of determining the support sets of system matrix M_Γ and have a great sparseness and computational efficiency both in CSRA framework and nCSRA framework. The runtime of homotopy was about 14.6s, while the runtime of ISTA was 56.8s. The computational efficiency of homotopy is about 4 times faster than ISTA. ISTA is an approximate iteration method with poor sparseness. Only in nCSRA framework, the result has a great sparseness. Because the sparseness of ISTA has a great changes in CSRA and nCSRA framework. So in the results part, the ISTA-CSRA vs ISTA-nCSRA shows much more difference than Homotopy-CSRA vs Homotopy-nCSRA. In theory, Approximate iteration strategies are need to build a appropriate approximates function (such as, Shrink (a, z) in this study) and more focus on analytical processing. Thus stability and robustness are the superiority of ISTA. Although algorithm characteristics of two algorithms are different, two algorithms all have a great improvement in nCSRA framework. In these points, our proposed framework has a great suitability, which can be generalized to some other inverse reconstruction algorithms and strategies.

In experimental design aspect, numerical simulations and in vivo experiments have been conducted to testify the feasibility of the proposed framework. Numerical simulations results illustrate that: (1) algorithms based on nCSRA can distinguish two sources with center separation of 2.4mm; (2) nCSRA has a great intensity resolution capability to recovery two sources of intensity ratio 1: 8 with RIE less than 30%; (3) compared with CSRA framework, the size resolution capability of nCSRA has a great improvement. It should be noted that the reconstructed source volume shown in Fig. 6 seems smaller than the value listed in Table 6. That is because of the Gaussian-like distribution of reconstructed intensity. In in vivo imaging experiments, considering the tumor separation, source intensity, tumor size, and source number, four tests based on 4T1 breast tumor models were conducted to assess the in vivo resolving power for multiple sources. And all results further verified the improvement of nCSRA framework in spatial resolution, intensity resolution, and size resolution. However limited by the penetration depth of CLI and the mouse model, only the subcutaneous tumor model has been tested. Fortunately, the nCSRA was not only applicable to CLT but also suitable for bioluminescence tomography (BLT), fluorescence molecular tomography (FMT), and some new imaging techniques based on Cerenkov-excited luminescence etc [46]. So one of our further work is distinguish deep tumors in radiopharmaceutical excited fluorescence tomography by using the proposed nCSRA framework. Our further work will also include the analysis and classification process of the reconstructed results and building the correlation between the details of sources and tumor staging.

In conclusion, an effective nCSRA framework was proposed to solve the multiple-source resolution problem in CLT. The results of numerical simulations and in vivo experiments demonstrated that nCSRA framework could significantly enhance the multiple-source resolution capability in aspect of spatial resolution, intensity resolution, and size resolution. It is believed that this framework can be applicable to some other optical tomography, and thus it will benefits various preclinical applications and facilitates the development of optical molecular tomography in theoretical study.

Funding

National Natural Science Foundation of China (NSFC) (81227901, 81527805, 61231004, 11571012, 61622117, 81671759, 61302024, 61372046); National Key Research and Development Program of China (2016YFC0102600); Strategic Priority Research Program from Chinese Academy of Sciences (XDB02060010); International Innovation Team of CAS (20140491524); Beijing Municipal Science & Technology Commission (Z161100002616022); Scientific Instrument Developing Project of the Chinese Academy of Sciences (YZ201672).

References and links

1. R. Robertson, M. S. Germanos, C. Li, G. S. Mitchell, S. R. Cherry, and M. D. Silva, “Optical imaging of cerenkov light generation from positron-emitting radiotracers,” Phys. Med. Biol. 54(16), N355–N365 (2009). [CrossRef] [PubMed]

2. A. Ruggiero, J. P. Holland, J. S. Lewis, and J. Grimm, “Cerenkov luminescence imaging of medical isotopes,” J. Nucl. Med. 51(7), 1123–1130 (2010). [CrossRef] [PubMed]

3. R. S. Dothager, R. J. Goiffon, E. Jackson, S. Harpstrite, and D. Piwnica-Worms, “Cerenkov radiation energy transfer (CRET) imaging: a novel method for optical imaging of pet isotopes in biological systems,” PloS One 5(10), e13300 (2010). [CrossRef] [PubMed]

4. A. E. Spinelli, D. D’Ambrosio, L. Calderan, M. Marengo, A. Sbarbati, and F. Boschi, “Cerenkov radiation allows in vivo optical imaging of positron emitting radiotracers,” Phys. Med. Biol. 55, 483–495 (2010). [CrossRef]

5. D. L. Thorek, D. S. Abou, B. J. Beattie, R. M. Bartlett, R. Huang, P. B. Zanzonico, and J. Grimm, “Positron lymphography: multimodal, high-resolution, dynamic mapping and resection of lymph nodes after intradermal injection of 18F-FDG,” J. Nucl. Med. 53(9), 1438–1445 (2012). [CrossRef] [PubMed]

6. R. Madru, T. A. Tran, J. Axelsson, C. Ingvar, A. Bibic, F. Stahlberg, L. Knutsson, and S.-E. Strand, “68Ga-labeled superparamagnetic iron oxide nanoparticles (spions) for multi-modality pet/mr/cherenkov luminescence imaging of sentinel lymph nodes,” Am. J. Nucl. Med. Mol. Imag. 4(1), 60 (2014).

7. H. Liu, C. M. Carpenter, H. Jiang, G. Pratx, C. Sun, M. P. Buchin, S. S. Gambhir, L. Xing, and Z. Cheng, “Intraoperative imaging of tumors using cerenkov luminescence endoscopy: a feasibility experimental study,” J. Nucl. Med. 53(10), 1579–1584 (2012). [CrossRef] [PubMed]

8. C. M. Carpenter, X. Ma, H. Liu, C. Sun, G. Pratx, J. Wang, S. S. Gambhir, L. Xing, and Z. Cheng, “Cerenkov luminescence endoscopy: Improved molecular sensitivity with β–emitting radiotracers,” J. Nucl. Med. 55(11), 1905–1909 (2014). [CrossRef] [PubMed]

9. T. Song, X. Liu, Y. Qu, H. Liu, C. Bao, C. Leng, Z. Hu, K. Wang, and J. Tian, “A novel endoscopic cerenkov luminescence imaging system for intraoperative surgical navigation,” Mol. Imaging 14(8), 443–449 (2015).

10. J. P. Holland, G. Normand, A. Ruggiero, J. S. Lewis, and J. Grimm, “Intraoperative imaging of positron emission tomographic radiotracers using cerenkov luminescence emissions,” Molec. Imag. 10(3), 177–186 (2011).

11. D. Fan, X. Zhang, L. Zhong, X. Liu, Y. Sun, H. Zhao, B. Jia, Z. Liu, Z. Zhu, J. Shi, and F. Wang, “68Ga-labeled 3PRGD2 for dual pet and cerenkov luminescence imaging of orthotopic human glioblastoma,” Bioconjugate Chem. 26(6), 1054–1060 (2015). [CrossRef]

12. A. E. Spinelli, M. Ferdeghini, C. Cavedon, E. Zivelonghi, R. Calandrino, A. Fenzi, A. Sbarbati, and F. Boschi, “First human cerenkography,” J. Biomed. Opt. 18, 020502 (2013). [CrossRef]

13. H. Hu, X. Cao, F. Kang, M. Wang, Y. Lin, M. Liu, S. Li, L. Yao, J. Liang, J. Liang, Y. Nie, X. Chen, J. Wang, and K. Wu, “Feasibility study of novel endoscopic cerenkov luminescence imaging system in detecting and quantifying gastrointestinal disease: first human results,” Eur. Radiol. 25(6), 1814 (2015). [CrossRef] [PubMed]

14. L. A. Jarvis, R. Zhang, D. J. Gladstone, S. Jiang, W. Hitchcock, O. D. Friedman, A. K. Glaser, M. Jermyn, and B. W. Pogue, “Cherenkov video imaging allows for the first visualization of radiation therapy in real time,” Int. J. Radiat. Oncol. Biol. Phys. 89(3), 615–622 (2014). [CrossRef] [PubMed]

15. M. R. Grootendorst, M. Cariati, S. E. Pinder, A. Kothari, M. Douek, T. Kovacs, H. Hamed, A. Pawa, F. Nimmo, J. Owen, V. Ramalingam, S. Sethi, S Mistry, K. Vyas, D. S. Tuch, A. Britten, M. V. Hemelrijck, G. J. Cook, C. Sibley-Allen, S. Allen, and A. Purushotham, “Intraoperative assessment of tumor resection margins in breast-conserving surgery using 18F-FDG cerenkov luminescence imaging: A first-in-human feasibility study,” J. Nucl. Med. 58(6), 891–898 (2017). [CrossRef]

16. B. J. Beattie, D. L. Thorek, C. R. Schmidtlein, K. S. Pentlow, J. L. Humm, and A. H. Hielscher, “Quantitative modeling of cerenkov light production efficiency from medical radionuclides,” PLoS One 7(2), e31402 (2012). [CrossRef] [PubMed]

17. D. L. Thorek, A. Ogirala, B. J. Beattie, and J. Grimm, “Quantitative imaging of disease signatures through radioactive decay signal conversion,” Nat. Med. 19(10), 1345–1350 (2013). [CrossRef] [PubMed]

18. C. Li, G. S. Mitchell, and S. R. Cherry, “Cerenkov luminescence tomography for small-animal imaging,” Opt. Lett. 35(7), 1109–1111 (2010). [CrossRef] [PubMed]

19. Z. Hu, M. Liu, Z. Zhang, H. Guo, and J. Tian, “A novel radiopharmaceutical-excited fluorescence tomography of the mice bearing hepatocellular carcinoma,” J. Nucl. Med. 57(2), 1421 (2016).

20. Z. Hu, X. Ma, X. Qu, W. Yang, J. Liang, J. Wang, and J. Tian, “Three-dimensional noninvasive monitoring iodine-131 uptake in the thyroid using a modified cerenkov luminescence tomography approach,” PloS One 7(5), e37623 (2012). [CrossRef] [PubMed]

21. F. Boschi and A. E. Spinelli, “Cerenkov luminescence imaging at a glance,” Curr. Mol. Imag. 3(2), 106–117 (2014). [CrossRef]

22. A. Spinelli, M. Marengo, R. Calandrino, A. Sbarbati, and F. Boschi, “a novel multimodal approach to molecular imaging,” Q. J. Nucl. Med. Molec. Imag. 56, 279–289 (2012).

23. C. L. Chaffer and R. A. Weinberg, “A perspective on cancer cell metastasis,” Science 331(6024), 1559–1564 (2011). [CrossRef] [PubMed]

24. F. L. Greene and L. H. Sobin, “The TNM system: our language for cancer care,” J. Surg. Oncol. 80(3), 119–120 (2002). [CrossRef] [PubMed]

25. C. Qin, J. Feng, S. Zhu, X. Ma, J. Zhong, P. Wu, Z. Jin, and J. Tian, “Recent advances in bioluminescence tomography: methodology and system as well as application,” Laser Photon. Rev. 8(1), 94–114 (2014). [CrossRef]

26. Z. Hu, J. Liang, W. Yang, W. Fan, C. Li, X. Ma, X. Chen, X. Ma, X. Li, X. Qu, J. Wang, F. Cao, and J. Tian, “Experimental cerenkov luminescence tomography of the mouse model with spect imaging validation,” Opt. Express 18(24), 24441–24450 (2010). [CrossRef] [PubMed]

27. H. Guo, X. He, M. Liu, Z. Zhang, Z. Hu, and J. Tian, “Weight multispectral reconstruction strategy for enhanced reconstruction accuracy and stability with cerenkov luminescence tomography,” IEEE T. Med. Imag. 36(6), 1337–1346 (2017). [CrossRef]

28. C. Ran, Z. Zhang, J. Hooker, and A. Moore, “In vivo photoactivation without light: use of cherenkov radiation to overcome the penetration limit of light,” Mol. Imaging Biol. 14(2), 156–162 (2012). [CrossRef]

29. K. Tanha, A. M. Pashazadeh, and B. W. Pogue, “Review of biomedical čerenkov luminescence imaging applications,” Biomed. Opt. Express 6(8), 3053–3065 (2015). [CrossRef] [PubMed]

30. R. Weissleder, “A clearer vision for in vivo imaging,” Nat. Biotechnol. 19(4), 346 (2001). [CrossRef]

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

32. J. Zhong, C. Qin, X. Yang, S. Zhu, X. Zhang, and J. Tian, “Cerenkov luminescence tomography for in vivo radiopharmaceutical imaging,” Int. J. Biomed. Imag. 2011, 641618 (2011).

33. Z. Hu, X. Chen, J. Liang, X. Qu, D. Chen, W. Yang, J. Wang, F. Cao, and J. Tian, “Single photon emission computed tomography-guided cerenkov luminescence tomography,” J. Appl. Phys. 112(2), 024703 (2012). [CrossRef]

34. Z. Hu, W. Yang, X. Ma, W. Ma, X. Qu, J. Liang, J. Wang, and J. Tian, “Cerenkov luminescence tomography of aminopeptidase N (APN/CD13) expression in mice bearing HT1080 tumors,” Molec. Imag. 12(3), 173–181 (2013).

35. H. Liu, X. Yang, T. Song, C. Bao, L. Shi, Z. Hu, K. Wang, and J. Tian, “Multispectral hybrid cerenkov luminescence tomography based on the finite element spn method,” J. Biomed. Opt. 20(8), 086007 (2015). [CrossRef]

36. X. Ding, K. Wang, B. Jie, Y. Luo, Z. Hu, and J. Tian, “Probability method for cerenkov luminescence tomography based on conformance error minimization,” Biomed. Opt. Express 5(7), 2091–2112 (2014). [CrossRef] [PubMed]

37. A. E. Spinelli, C. Kuo, B. W. Rice, R. Calandrino, P. Marzola, A. Sbarbati, and F. Boschi, “Multispectral cerenkov luminescence tomography for small animal optical imaging,” Opt. Express 19(13), 12605–12618 (2011). [CrossRef] [PubMed]

38. M. Liu, H. Guo, H. Liu, Z. Zhang, C. Chi, H. Hui, D. Dong, Z. Hu, and J. Tian, “In vivo pentamodal tomographic imaging for small animals,” Biomed. Opt. Express 8(3), 1356–1371 (2017). [CrossRef] [PubMed]

39. X. Chen, D. Yang, Q. Zhang, and J. Liang, “L1/2 regularization based numerical method for effective reconstruction of bioluminescence tomography,” J. Appl. Phys. 115(18), 184702 (2014). [CrossRef]

40. D. Zhu and C. Li, “Nonconvex regularizations in fluorescence molecular tomography for sparsity enhancement,” Phys. Med. Biol. 59(12), 2901 (2014). [CrossRef] [PubMed]

41. H. Guo, J. Yu, X. He, Y. Hou, F. Dong, and S. Zhang, “Improved sparse reconstruction for fluorescence molecular tomography with L1/2 regularization,” Biomed. Opt. Express 6(5), 1648–1664 (2015). [CrossRef] [PubMed]

42. Y. Lv, J. Tian, W. Cong, G. Wang, J. Luo, W. Yang, and H. Li, “A multilevel adaptive finite element algorithm for bioluminescence tomography,” Opt. Express 14(18), 8211–8223 (2006). [CrossRef] [PubMed]

43. C. Leng and J. Tian, “Mathematical method in optical molecular imaging,” Sci. China Inform. Sci. 58(3), 1–13 (2015). [CrossRef]

44. J. He, “Homotopy perturbation technique,” Comput. Method. Appl. M. 178(3), 257–262 (1999). [CrossRef]

45. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci. 2(1), 183–202 (2009). [CrossRef]

46. Z. Hu, Y. Qu, K. Wang, X. Zhang, J. Zha, T. Song, C. Bao, H. Liu, Z. Wang, J. Wang, Z. Liu, H. Liu, and J. Tian, “In vivo nanoparticle-mediated radiopharmaceutical-excited fluorescence molecular imaging,” Nat. Commun. 6, 7560 (2015). [CrossRef]

Material	µ_a(mm⁻¹)	µ_s(mm⁻¹)	g
Muscle	0.12	15.58	0.97
Heart	0.08	10.07	0.90
Lung	0.26	31.56	0.93
Liver	0.47	10.00	0.93
Kidney	0.09	23.59	0.90
Bone	0.08	10.07	0.93

Mouse	Total number of injected cells	Number of injections	Separation of injection sites	Ratio of injection dose
No. 1	1.00 × 10⁶	2	5mm	1:1
No. 2	1.50 × 10⁶	2	5mm	2:1
No. 3	0.75 × 10⁶	2	3mm	2:1
No. 4	2.00 × 10⁶	3	3mm	1:2:3

Separation	Framework	Method	Reconstructed center(mm)		LE(mm)

			No. 1	No. 2	No. 1	No. 2
4.8(mm)	nCSRA	Homotopy	(18.0, 19.8, 6.8)	(13.0, 19.8, 6.7)	0.25	0.36
	nCSRA	ISTA	(17.7, 19.9, 6.8)	(12.7, 19.9, 7.2)	0.44	0.53

	CSRA	Homotopy	(18.0, 19.8, 6.8)	(13.0, 19.8, 6.7)	0.25	0.36
	CSRA	ISTA	(17.4, 19.6, 7.0)	(13.1, 19.5, 7.0)	0.59	0.15

3.6(mm)	nCSRA	Homotopy	(18.0, 19.8, 6.8)	(14.7, 19.8, 7.3)	0.25	0.45
	nCSRA	ISTA	(18.0, 19.8, 6.8)	(14.7, 19.8, 7.3)	0.25	0.45

	CSRA	Homotopy	(18.0, 19.8, 6.8)	(15.0, 19.5, 7.2)	0.25	0.60
	CSRA	ISTA	(18.0, 19.8, 6.8)	(15.4, 19.0, 7.4)	0.25	1.24

2.4(mm)	nCSRA	Homotopy	(16.2, 19.5, 6.8)	(13.6, 16.7, 7.1)	0.62	0.83
	nCSRA	ISTA	(16.6, 20.0, 7.1)	(13.4, 20.3, 6.3)	0.41	1.34

	CSRA	Homotopy	(16.4, 18.8, 6.2)	(15.4, 19.1, 7.0)	1.20	1.11
	CSRA	ISTA	(15.7, 19.2, 6.8)	(15.4, 19.0, 7.4)	1.18	1.22

Intensity Ratio	Framework	Method	Reconstructed Intensity(nw/mm³)		RIE(%)

			No. 1	No. 2	No. 1	No. 2
1:2	nCSRA	Homotopy	11.65	22.26	16.5	11.3
	nCSRA	ISTA	12.30	21.33	23.0	5.2

	CSRA	Homotopy	7.75	27.04	22.5	35.0
	CSRA	ISTA	5.97	13.66	40.3	31.7

1:4	nCSRA	Homotopy	11.81	43.17	18.1	7.9
	nCSRA	ISTA	12.21	41.42	22.1	3.6

	CSRA	Homotopy	4.27	38.62	57.3	3.5
	CSRA	ISTA	4.38	14.50	56.2	63.8

1:8	nCSRA	Homotopy	9.28	75.81	7.7	5.2
	nCSRA	ISTA	13.21	72.42	32.1	9.5

	CSRA	Homotopy	2.91	65.82	70.9	17.7
	CSRA	ISTA	4.38	14.49	56.2	81.9

Volume Ratio	Framework	Method	Reconstructed volume(mm³)		RVE(%)

			No. 1	No. 2	No. 1	No. 2
3.09:6.03	nCSRA	Homotopy	3.88	5.43	25.57	9.95
	nCSRA	ISTA	3.88	3.83	25.57	36.48

	CSRA	Homotopy	8.24	5.43	166.67	9.95
	CSRA	ISTA	11.93	14.72	286.08	144.11

3.09:1.30	nCSRA	Homotopy	2.59	1.40	16.18	7.69
	nCSRA	ISTA	2.59	1.40	16.18	7.69

	CSRA	Homotopy	4.00	3.69	29.45	183.85
	CSRA	ISTA	6.21	8.12	100.97	524.62

Non-convex sparse regularization approach framework for high multiple-source resolution in Cerenkov luminescence tomography

Abstract

1. Introduction

2. Materials and methods

2.1. Photon propagation model

2.2. Non-convex sparse regularization algorithm (nCSRA) framework

2.2.1. Homotopy

2.2.2. Iterative shrinkage-thresholding algorithm (ISTA)

2.3. Numerical simulations

2.3.1. Spatial resolution

2.3.2. Intensity resolution

2.3.3. Size resolution

2.4. In vivo imaging experiments

2.5. Statistical analysis

3. Results

3.1. Numerical simulation

3.1.1. Spatial resolution

3.1.2. Intensity resolution

3.1.3. Size resolution

3.2. In vivo imaging experiments

4. Discussion and conclusion

Funding

References and links

Cited By

Figures (8)

Tables (6)

Equations (11)

Optics Express

Group	ID	source center (mm)	Separation (mm)	Intensity (nw/mm³)	Volume (mm³)
Spatial resolution	No. 1	(18.0, 19.6, 7.0)	4.8	10	3.09
	No. 2	(13.2, 19.6, 7.0)	4.8	10	3.09

	No. 1	(18.0, 19.6, 7.0)	3.6	10	3.09
	No. 2	(14.4, 19.6, 7.0)	3.6	10	3.09

	No. 1	(16.8, 19.6, 7.0)	2.4	10	3.09
	No. 2	(14.4, 19.6, 7.0)	2.4	10	3.09

Intensity resolution	No. 1	(18.0, 19.6, 7.0)	4.8	10	3.09
	No. 2	(13.2, 19.6, 7.0)	4.8	20	3.09

	No. 1	(18.0, 19.6, 7.0)	4.8	10	3.09
	No. 2	(13.2, 19.6, 7.0)	4.8	40	3.09

	No. 1	(18.0, 19.6, 7.0)	4.8	10	3.09
	No. 2	(13.2, 19.6, 7.0)	4.8	80	3.09

Size resolution	No. 1	(18.0, 19.6, 7.0)	4.8	10	3.09
	No. 2	(13.2, 19.6, 7.0)	4.8	10	6.03

	No. 1	(18.0, 19.6, 7.0)	4.8	10	3.09
	No. 2	(13.2, 19.6, 7.0)	4.8	10	1.30