Elastic net-based non-negative iterative three-operator splitting strategy for Cerenkov luminescence tomography

: Cerenkov luminescence tomography (CLT) provides a powerful optical molecular imaging technique for non-invasive detection and visualization of radiopharmaceuticals in living objects. However, the severe photon scattering effect causes ill-posedness of the inverse problem, and the location accuracy and shape recovery of CLT reconstruction results are unsatisfactory for clinical application. Here, to improve the reconstruction spatial location accuracy and shape recovery ability, a non-negative iterative three operator splitting (NNITOS) strategy based on elastic net (EN) regularization was proposed. NNITOS formalizes the CLT reconstruction as a non-convex optimization problem and splits it into three operators, the least square, L 1 / 2 -norm regularization, and adaptive grouping manifold learning, then iteratively solved them. After stepwise iterations, the result of NNITOS converged progressively. Meanwhile, to speed up the convergence and ensure the sparsity of the solution, shrinking the region of interest was utilized in this strategy. To verify the effectiveness of the method, numerical simulations and in vivo experiments were performed. The result of these experiments demonstrated that, compared to several methods, NNITOS can achieve superior performance in terms of location accuracy, shape recovery capability, and robustness. We hope this work can accelerate the clinical application of CLT in the future.

spatial positioning accuracy and shape recovery [16], which is a great challenge to the clinical application of CLT. To improve the performance of CLT reconstruction, some anatomical prior information and regularization methods can be added to CLT reconstruction. In the first CLT reconstruction experiment of small animals, it is assumed that the optical properties of mouse tissues are consistent and uniform [12], but this assumption is not consistent with the real situation. Later, Hu et al. proposed a heterogeneous mouse model. By using mouse anatomical information and different optical properties of organs, this model can reduce the systematic error of CLT and significantly improve the accuracy of tumor localization [17]. Moreover, some regularization theories are gradually applied to CLT reconstruction to ensure the accuracy of Cerenkov source, such as L 2 -norm regularization (Tikhonov method) [12], L 1 -norm regularization (Lasso method) [18,19] and L p -norm (0 < p < 1) regularization (non-convex method) [20]. However, these regularization methods have their defects. L 2 -norm regularization will introduce many artifacts of reconstruction, resulting in excessive smoothing of the results [21]. L 1 -norm regularization and L p -norm regularization cause other problems in the reconstruction results, such as over-sparseness and lack of fluorescence probe boundary information [22]. To overcome these shortcomings, Some methods based on joint regularization were proposed, such as difference of convex algorithm (DCA) based on L 1−2 -norm [23], Sparse-graph manifold learning (SGML) method [24] and so on [25][26][27]. In addition, many researchers found that shrinking the region of interest (ROI) can also improve the reconstruction accuracy of optical tomography reconstruction [28,29]. These works provide some new ideas for alleviating the ill-posedness of CLT reconstruction. In addition to traditional mathematical methods, deep learning strategies are also introduced in CLT [30][31][32]. In 2019, a CLT approach based on multi-layer fully connected networks was proposed [33]. Subsequently, the attention mechanism of deep learning and decoder-encoder and other structures are also used to reconstruct the Cerenkov source [30,32,34]. Although the neural network method improves the performance of CLT reconstruction, it not only needs a huge amount of training data but also has weak generalization ability, that is, the trained artificial neural network can only be applied to specific imaging objects [35].
In this work, a non-negative iterative three-operator splitting (NNITOS) approach based on elastic net (EN) was proposed to improve the accuracy and robustness of CLT reconstruction. EN regularization combined the advantages of L 1/2 -norm regularization and manifold learning regularization. The inverse problem of CLT source reconstruction was transformed by NNITOS into a nonconvex, nonnegative iterative, and mixed regularization minimization problem. Based on EN regularization, the objective minimum optimization function of the inverse problem was split into the least square operator, the L 1/2 -norm regularization operator, and the manifold learning regularization operator using NNITOS. The optimal solution of this objective optimization function was approximated by solving the indeterminate point equation, and because the Cerenkov source was non-negativity, a non-negative constraint was added to the solution of this method. To further speed up the convergence of the solution of the indeterminate point equation, we used a shrinking ROI strategy to constrain the region of the next step solution based on the results of the previous step during the iterative process of the solution.
To evaluate the performance of NNITOS based on EN regularization in CLT reconstruction, numerical simulation and in vivo experiments were carried out. Compared with fast iterative shrinkage thresholding algorithm (FISTA) based on L 1 -norm [36], incomplete variables truncated conjugate gradient (IVTCG) algorithm based on L 1 -norm [37], Gaussian weighted Laplace prior (GWLP) based on L 2 -norm [38]. The experimental results show that the reconstruction accuracy and shape recovery ability are greatly improved by NNITOS.
The structure of the rest of this paper is as follows. Section 2 introduces the CLT forward model, CLT inverse problem and NNITOS strategy. Section 3 introduces the process and results of numerical simulation and in vivo experiments. Section 4 makes some discussion and summary of this paper.

Forward model of CLT
The transmission of Cerenkov photons in biological tissue accords with the Radiative Transfer Equation (RTE), but it has high computational complexity. In the process of biological tissue transmission, the scattering rate of photons is greater than the absorption rate [39], so RTE can be simplified as a Diffusion Approximation (DA) model. The DA model combined with the Robin boundary condition can be expressed as [40,41]: where r represents the nodes in the solution domain Ω. ∂Ω is boundary of the domain. Φ(r) denotes the photon flow rate at the node r. S(r) represents Cerenkov source distribution. D(r) = 1/3(µ α (r) + (1 − g)µ s (r)), where µ α (r) and µ s (r) are optical absorption coefficient and optical scattering coefficient of photons at node r, respectively, and g is the medium anisotropy factor. A(r, n, n ′ ) ≈ 1+R(r) 1−R(r) , where R(r) is the inner reflection coefficient, n and n ′ are the refractive index within Ω and the refractive index of external medium, respectively. v(r) denotes the unit outer normal on ∂Ω. Based on the finite element framework [41], the light propagation model can be discretized, and the relationship between the unknown Cerenkov source distribution and the luminous flux on the surface of biological tissue can be described as the linear equation as: where Y is the luminous flux measured on the surface, U denotes the system weight matrix and X represents the unknown internal Cerenkov source distribution. Therefore, the purpose of the CLT inverse problem is to restore the radioactive probe distribution X in the above linear equation, and the detailed information of the forward model can be found in [42].

Inverse problem
As mentioned above, the linear relationship between the Cerenkov luminescence source and the surface luminous flux in biological tissues has been established. As the result of the non-negativity of the Cerenkov source distribution, and inspired by Three-Operator Splitting (TOS) algorithm [43,44], the NNITOS approach for CLT was proposed to alleviate the ill-posedness of the inverse problem. The TOS algorithm is described by the most straightforward mathematical formula as: where Φ is expressed as our objective function, f is potentially nonconvex and Lipschitz differentiable, and g and h are convex functions and potentially non-smooth. Combining the TOS algorithm and the objective function of CLT reconstruction, our proposed objective function can be expressed as follows: where λ 1 and λ 2 are the regularization parameters, ∥X∥ is the Laplacian manifold regularization term and L is the dynamic graph Laplacian matrix. This combination of different norms is called EN regularization [45,46]. The details of the EN-based NNITOS method are described below.

EN regularization
The EN regularization proposed by us included L 1/2 -norm regularization and manifold regularization. The EN combines two different regularizations, which can combine the advantages of each regularization, so it finds a balance between the sparsity and stationarity of the results. L p -norm regularization can produce more sparse solutions than L 1 -norm regularization, and L 1/2 -norm regularization is representative among all L p -norm regularizations with p in (0, 1) [47]. Thus, we choose L 1/2 -norm regularization instead of L 1 -norm regularization. As a popular learning method, manifold learning is often used to constrain local linear relations [48,49]. We can use manifold learning techniques to explore the geometric characteristics of data distribution [50]. To improve the accuracy of the solution space, in the manifold regularization term, we propose an adaptive grouping that adjusts the grouping according to the intensity of the grid nodes, as shown in Fig. 1. In the finite element framework, the strength of each mesh node is closely related to the tetrahedron, and the central nodex i is shared by multiple tetrahedrons in Fig. 1. Therefore, the strength of the i-th node is transformed into the strength of the i-th group, the n i represents the initial length of the i-th group, and the threshold is E i = X i /n i , where X i is the sum of the energy values of the i-th group. Nodes whose intensity is less than the threshold will be removed from the group. The initial length of each group may be different, but in this paper, n i = 9, refer to related work [51]. In the adaptive grouping, we construct the adjacency matrix g. If the i-th node and the j-th node were in the same group, then g ij = 1, otherwise, g ij = 0. represents the initial group, with the red node as the center nodex i , the yellow node surrounded by the green oval region is the neighbors ofx i , the g ij = 0 of the blue node. (B) indicates that after the adjustment of the group, the nodes whose energy value is less than the threshold E i are excluded from the green ellipse region.
The adjacency matrix g is needed in solving the dynamic graph Laplacian matrix L in EN regularization. The Laplacian manifold regularization term is defined as: where x i and x j denote i-th node and the j-th node, N is the number of finite elements meshes, and w ij is the weight of the edge between node N i and node N j . The w ij is defined as: where ∥︁ ∥︁ N i − N j ∥︁ ∥︁ 2 2 is the Euclidean distance between the i-th node and the j-th node, R is the radius of the Gaussian kernel function which adjusts the probability density function of adjacent nodes, g ij = 1 means that the i-th node and the j-th node are in the same group, whereas g ij = 0 means that the i-th node and the j-th node are not in the same group, as shown in Fig. 1. Obviously, w ij ∈ [0, 1]. After a simple mathematical derivation, Eq. (4) can be simplified to 2X T (D − W)X, where W is an asymmetrical weight matrix, which equivalents to (w ij ) NN . D is the entry diagonal matrix of the vertex, which is defined as: Let L = D − W, we will get the dynamic graph Laplacian matrix proposed by Eq. (4).

NNITOS strategy based on EN
In this section, NNITOS will be explained in detail. When we combine Eq. (3) and Eq. (4), it can be expressed as follows: where λ 1 and λ 2 are selected empirically, which are given in related works [52,53]. It is noted that NNITOS solved our objective function by evaluating the terms h and (f + g) at two different points, x and z.
can be equal to the optimal object value, even if x and z are not close to the solution. We can introduce the distance condition between x and z to address this issue. The following definition of an α-close and β-stationary pair of points is crucial for solving the Eq. (3): where ⟨·, ·⟩ denote the standard Euclidean inner product associated with the norm ∥·∥. α-close β-stationary points (x,z) yield approximate solution to Eq. (3) under appropriate assumptions. For specific assumptions, please refer to the Ref. [54]. By Eq. (9) and Eq. (10) we can know that TOS algorithm is approximately stationary. The TOS algorithm is represented by a fixed-point equation as follows: where the proximal operator of functionġ is defined by prox γ kġ (y k ) := argmin x k is the intermediate solution in the iterative process, γ k is the step-size, k denotes the number of iterative rounds. Therefore, Eq. (11) is denoted as: By designing a well-behaved fixed-point equation, the TOS method decomposes the minimization of the objective function into solving the nearest neighbor mapping of two individuals. When we combine Eq. (8) and Eq. (12), by solving the partial differential, it can be expressed as follows: where the I denotes the unit matrix of N × N. The value of γ 1 is selected refer to related work [44]. We add non-negative constraints to the TOS algorithm. In addition, NNITOS method uses ROI framework. In each iteration, we shrink the ROI based on the results of the iteration. Therefore, the NNITOS method for solving the inverse CLT problem can be summarized as Algorithm 1.

Input:
The system matrix U and the measured luminous flux vector Y.
Step1: Find the nine nearest neighbors of each node and assign them to the same group, construct the adjacency matrix g. calculate the threshold E i for each group.g ij = 0 denote the j-th node is not in the i-th group.

Repeat
Step3: Utilize the Eq. (13) to update x k , z k , y k+1 , The point less than 0 in x k is set to 0.
Step4: Calculate the cosine similarity EC k+1 and Euclidean distance EL k+1 between Ux k and Y.
Step5: Shrinking the ROI according to the reconstruction result x k , set the node with lower energy value in x k to 0.
Step6: Increase the iteration index k = k + 1, and then update step-size Output: x k .

Experiments and results
To verify and systematically evaluate the performance of NNITOS in CLT reconstruction, four groups of simulation experiments and one group of in vivo experiments were introduced in this paper. Three types of algorithms were compared with our proposed methods, namely, FISTA algorithm based on L 1 -norm, IVTCG algorithm based on L 1 -norm, and GWLP based on L 2 -norm. The regularization parameters of these algorithms come from related work [26,18]. In addition, to assess the stability, we also carried out a group of anti-Gaussian noise experiments with dual-source. All the programs in this work were executed on desktop computers with an Intel Core i3-10100 CPU (3.60 GHz) and 8GB RAM.

Evaluation metrics
To evaluate the location accuracy, shape recovery ability, and source intensity accuracy of CLT reconstruction in ROI, the location error (LE), DICE similarity coefficient (DICE), relative intensity error (RIE), and contrast to noise ratio (CNR) were taken as the quantitative evaluation indexes. LE represents the error between the center position of the reconstruction source and the center position of the real source, which is defined as: where X r denotes the vector of the reconstructed source and X t denotes the vector of the true source. DICE is the similarity between the reconstructed source region and the real source region: where R r denotes the region of the reconstructed source and R t denotes the region of the true Cerenkov source. DICE values range from 0 to 1, and the closer the DICE is to 1, the higher the shape recovery of the reconstruction source. RIE is used to evaluate the relative intensity error between the reconstructed Cerenkov source intensity and the actual Cerenkov source intensity: where I r denotes the intensity of the reconstructed source and I t denotes the intensity of the true Cerenkov source. Obviously, the closer the RIE is to 0, the better the intensity recovery of the reconstructed source. The contrast to noise ratio (CNR) is used to evaluate the quality of reconstruction results, and CNR is defined as follows: where µ ROI and µ BCK represent the mean values in the ROI and background, respectively; σ 2 ROI and σ 2 BCK are represent the variances; and m ROI and m BCK represent the number of nodes included in the ROI and background, respectively. The higher value of CNR indicates better quality of the reconstruction.

Numerical simulation
A series of numerical simulations are carried out based on the cylinder model, the main components of the cylinder model are shown in Fig. 2(A). It consists of five organs: muscle, heart, lung, liver, and bone. Moreover, all experiments in this work are based on single spectrum (650 nm). Here, the optical parameters of these organs at 650 nm are consistent with those in Ref. [35,55], as shown in Table 1. In all simulations, for the inverse reconstruction, the cylinder model was discretized into a uniform tetrahedral mesh, including 4626 nodes and 25840 tetrahedral elements using Comsol Multiphysics software [56], and the mean element edge size is 1.37mm, as shown in Fig. 2(B). the molecular optical simulation environment (MOSE) platform [57] based on the Monte Carlo method was used to simulate the surface flux distribution, as shown in Fig. 2(C).  We simulated Cerenkov sources of different sizes, shapes, and positions in organisms. The specific shapes, sizes, and positions of these light sources were shown in Table 2. Four groups of simulation experiments were designed. In one group of experiments, a spherical single light source with a radius of 1mm and 1.5mm was placed at the central coordinate (3, 3, 5) mm, while in the other group, a spherical single light source and a cube single light source with the same radius were placed at the central coordinate (0, -1, 25) mm. These experiments evaluate the shape recovery accuracy of NNITOS by simulating a single source with different positions and different shapes. We also carried out a dual-source simulation experiment to assess the spatial positioning accuracy. The central coordinates of the dual-source were (3, 3, 8) mm and (3, 3, 16) mm, respectively. Finally, to assess the robustness of NNITOS, we designed a set of anti-noise simulation experiment. Specifically, 5%, 10%, 15%, 20%, 25% Gaussian noise was added to the 1.0 mm spherical dual-source simulation experiment with central coordinates of (0, -2, 15) mm and (0, -2, 25) mm.

Implanted animal experiments
To further access the performance of our proposed method, we used the experimental data of implanted adult nude mouse (about 6-8 weeks old/female) and used a CLT/micro-CT dualmodality imaging system to acquire experimental data. Implanted animal experiments were carried out under 2% isoflurane-oxygen mixture gas anesthesia to minimize the suffering of the mouse. All procedures were conducted under the protocols approved by the Animal Ethics Committee of the Northwest University of China and Use Committee. In the experiment, firstly, we used a spherical glass vessel with a radius of 1 mm containing about 800 ± 50µCi 18 F-FDG as the artificial Cerenkov luminescence source, and then implanted it into mouse at the designed locations. Secondly, we began to collect the dual-modality raw data. The mouse was placed on the automatic rotating stage, which was rotating 360 degrees with 1°intervals to capture the images of CLI and X-ray projection. By using band-pass filters (FF01-650/13-25, Semrock, USA). The CLI signal with wavelength of 650nm and white light data were collected by a cooled high-sensitivity EMCCD (-80°C, iXonEM Ultra 888, UK). Before acquisition process of CLT, it is set that long exposure time (5min), high given value (300), high shift speed (13µs), and low speed readout rate (1MHz at 16bits). Then CT volume data was acquired by a micro-CT system (tube voltage of 40kV p, tube current of 300mA).
After in vivo imaging, we proceeded with data processing, as shown in Fig. 3(C), the trunk sections of the mouse with a height of 41mm was obtained. We segmented the main organs of the mouse, including muscle, lung, heart, stomach, liver, and kidney, and integrated these organs into the mouse model. The optical parameters of each organ are consistent with those in the literature [58]. It should be noted that, the mouse was used for three times, and each time the location of the implanted radionuclide probe is different ((17.5, 23.0, 21.0) mm, (14.5, 26.0, 33.0) mm, and (14.5, 20.0, 36.0) mm)). Moreover, the CLI images obtained from each implanted model was mapped and registered to the surface of the mouse torso. To assess the feasibility of this method in vivo, NNITOS was also compared with the other three algorithms mentioned above.

Multi-sized and multi-shape single-source reconstruction
In terms of shape restoration ability, we compare the performance of NNITOS and the above three algorithms in a single source of different sizes and different shapes. The results of 3D reconstruction and the shape of the axial plane prove the effectiveness of our proposed method. As shown in Fig. 4 and Fig. 5, the reconstructed Cerenkov source was represented by a red region under 3D-view. Figure 4 showed the comparative results of using different methods to reconstruct a single source of different sizes, and Fig. 5 showed the comparative results of using different methods to reconstruct a single source of different shapes. These two experiments showed that the reconstruction result of the NNITOS method is closer to the real source. According to Table 3 and Table 4, we found that the NNITOS had a smaller LE, so the reconstructed Cerenkov source center is closer to the real value, and the DICE of NNITOS is smaller, indicating that the reconstructed Cerenkov source region and the actual Cerenkov source region have a good overlap, and the RIE of NNITOS is lower, means that the energy intensity recovery rate of the reconstructed source is very high.

Dual-source reconstruction
To further verify the positioning accuracy of NNITOS reconstruction, we have established a dual-source experiment. Different methods of reconstruction of Cerenkov dual spherical sources are shown in 3D, axial, sagittal and coronal views as shown in Fig. 6. FISTA and IVTCG produce excessive smoothing artifacts around a single source. Although the reconstruction performance of GWLP is better than that of the former two methods, compared with NNITOS, the center position of NNITOS reconstruction is more accurate, the image artifacts are less, and the light source intensity is closer to the real value. Specifically, the quantitative analysis of the reconstruction results of dual-source simulation is shown in Table 5.

Anti-noise performance
To test the robustness of NNITOS, we established an anti-noise experiment. The dual-source simulation as an example, 5% 10% 15% 20% 25% Gaussian noise was added to detect the impact on NNITOS reconstruction results. As shown in Fig. 7, the values of LE, DICE, RIE and CNR fluctuate little, indicating that the method is robust to noise.

Implanted animal experiment
To evaluate the practicability of NNITOS imaging in vivo, we used the experimental data of a mouse with radioactive probes in vivo. According to the mouse trunk boundary registration, the reconstruction results obtained by different methods are shown in Fig. 8. The three-dimensional reconstruction results and axial plane prove the effectiveness of different methods, our proposed method successfully overcame the over-smooth problem in vivo experiment. NNITOS has the smallest LE, the largest source matching area, and the best energy recovery ability, and  the reconstruction performance is better than other methods. The quantitative analysis of the reconstruction results of the real implanted animal experiment is shown in Table 6.

Discussion and conclusion
CLT can be used to visualize the three-dimensional distribution of radioactive probes in living animals, which solves the problem of insufficient depth information in traditional CLI methods. At present, it is considered that CLT is a highly promising imaging technology for clinical application. However, due to the ill-posedness of the inverse problem caused by the strong scattering of Cerenkov light in vivo, the performance of CLT reconstruction is not satisfactory. Further improving the reconstruction performance is of great significance to promote the application of CLT in preclinical and clinical research. Therefore, this paper proposes an NNITOS method to achieve better CLT reconstruction performance. EN regularization combines the L 1/2 -norm and manifold learning regularization to balance the sparsity, smoothness, and edge preservation of Cerenkov sources. Because the regularization term is complex, we use the NNITOS method to split the inverse problem into three operators: the least square operator, the L 1/2 -norm operator, and the manifold learning operator to reduce the computational complexity. On this basis, we use the method of shrinking ROI, which shrinks ROI based on the prior information generated by the previous iteration, and further controls the sparsity of the reconstructed Cerenkov source. The non-negative constraint of each iteration also ensures the sparsity of the reconstructed Cerenkov source, reduces the error, and helps to improve the performance of CLT reconstruction.
To evaluate the performance of NNITOS in CLT reconstruction, we used FISTA, IVTCG and, GWLP methods to compare with NNITOS, and carried out numerical simulation and in vivo experiments. The experimental results show that: 1) multi-sized source and multi-shape source simulation experiments show that NNITOS also has advantages in restoring source shape. 2) Dual-source simulation experiments show that NNITOS guarantees high reconstruction localization accuracy. 3) NNITOS exhibits good robustness in experimental tests against noise. 4) All simulation experiments show that NNITOS has better intensity recovery ability for Cerenkov sources than other methods. 5) In vivo experiments demonstrate the feasibility of NNITOS in biomedical research.
In addition, we also find that the decrease of position error may not cause the increase of the DICE coefficient or the decrease of the RIE coefficient. That is, there is a situation in which the centers of two regions are very close, but there may be great differences in shape, size, or strength between them. For example, in the spherical single source reconstruction in multi-shape simulation, compared with the FISTA method, the IVTCG method has a lower location error, but a higher RIE coefficient. Compared with the GWLP, the IVTCG method has a higher location error, but a higher DICE coefficient. Compared with other methods, NNITOS performs very well in these three evaluation indicators.
Although the NNITOS strategy performs well in CLT reconstruction, it still has some shortcomings. Firstly, in the manifold learning regularization, to group these nodes, we must find out the k-nearest nodes of each node, which is the initial group, and then automatically adjust the grouping according to the threshold. The selection of k will affect the accuracy of manifold learning gradient calculation, k is larger, the calculation accuracy is greater, but it takes more time. Therefore, how to determine the value of k of the initial group in manifold learning needs further research. In addition, the regularization parameters of NNITOS and some other internal parameters are chosen empirically. It is necessary to propose an adaptive parameter selection algorithm to objectively set these parameters. Finally, the clinical application of CLT based on the NNITOS strategy needs to be further studied.
Overall, an EN regularization-based NNITOS strategy is proposed for CLT morphological reconstruction. Compared with other traditional methods mentioned above, the NNITOS strategy achieved more precise results in the localization and morphological recovery of radioactive probes. And it has good robustness in the reconstruction. We believe that this strategy is beneficial to the preclinical and clinical research of CLT and promotes the development of theoretical research of optical tomography. Disclosures. The authors declare no conflicts of interest.
Data availability. Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon request.