Bioluminescence tomography (BLT) is an effective noninvasive molecular imaging modality for in vivo tumor research in small animals. However, the quality of BLT reconstruction is limited by the simplified linear model of photon propagation. Here, we proposed a multilayer perceptron-based inverse problem simulation (IPS) method to improve the quality of in vivo tumor BLT reconstruction. Instead of solving the inverse problem of the simplified linear model of photon propagation, the IPS method directly fits the nonlinear relationship between an object surface optical density and its internal bioluminescent source. Both simulation and orthotopic glioma BLT reconstruction experiments demonstrated that IPS greatly improved the reconstruction quality compared with the conventional approach.
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
4 January 2019: Typographical corrections were made to the author affiliations.
In vivo bioluminescence tomography (BLT) was developed to realize quantitative three-dimensional (3D) whole body imaging and is frequently used for the noninvasive study of cancer biological behavior [1–3]. Because bioluminescence is only emitted in viable tumor cells [1,2], BLT is particularly important for the study of tumors in small animal experiments. It provides sufficient sensitivity and specificity for identifying tumor position, and the 3D pinpointing of tumor location is very important for tumor behavior research.
The conventional model-based BLT reconstruction depends on the description of photon propagation in biological tissue. It utilizes the first-order approximate model of radiative transfer equation (RTE) [1,4,5] to modeling the photon propagation. With the application of finite element analysis, a linear relationship between the photon density at the surface of an imaged object and the photon density of bioluminescent source inside the object can be transformed from the first-order approximate model as follows:1). However, because system matrix is a simplified linear approximation model of the nonlinear RTE, it inevitably causes errors in BLT reconstructions. Furthermore, the ill-posed inverse problem of the model-based reconstruction also enhances the reconstruction difficulty [6,7].
To reduce the modeling error and overcome the ill-posed problem in BLT reconstruction, different research groups have proposed several approaches. Anatomical structures of organs and tissues, obtained from computerized tomography (CT) or magnetic resonance imaging (MRI) and their corresponding optical coefficients were used as prior information to improve the description of photon propagation in animal bodies [6–14]. High-order approximation models  were utilized to construct the system matrix . Besides these, a great amount of effort was made for optimization using different regularization techniques to cope with the challenge of the ill-posed inverse problem. Sparse prior information was one of the major assumptions used to construct such regularization, and various sparse algorithms were designed to characterize the results under this prior information. They either utilized greedy strategies to enhance the sparsity of results [9,10] or penalized the reconstruction with sparse regularization terms (, , etc.) [6,11–14]. Although the positioning accuracy of the reconstructed bioluminescent source was gradually improved by these methods, it still is essentially influenced by the deviation between the simplified linear model of RTE and the true photon propagation. This fundamental problem has not been properly solved ever since BLT was proposed.
In this Letter, we propose an inverse problem simulation (IPS) method based on the multilayer perceptron (MLP) strategy to reconstruct the density of bioluminescent source. Simulated and in vivo orthotopic glioma mouse models were employed to validate our method. MLP is a machine-learning network designed to construct a nonlinear mapping from input to output . Parameters of a nonlinear mapping are obtained from statistical learning. In contrast to conventional model-based BLT reconstruction, the MLP-based IPS method constructs the inverse of the photon propagation directly by learning the nonlinear mapping relationship between the surface photon density and the bioluminescent source density. Thus, IPS reconstructs the bioluminescent source without involving any forward photon propagation modeling, which fundamentally removed deviations originated from the simplification of RTE and the ill-posed model-based inverse problem.
The IPS method network contains one input layer, one output layer, and four hidden layers. The network input is the density of surface photon , and the network output is reconstructed bioluminescent source . The output of layer is strongly connected to the input of the next layer and added with bias. Connection weights and bias are learned from the network training. The activation function following the connection is called rectified linear unit (ReLU)  and is expressed as16] is utilized behind the activation function to reduce the overfitting problem of the network. In order to reduce the extremely expensive computational cost for establishing the 3D mesh of every single mouse from its CT or MRI data volume, a standardized 3D mesh of a mouse head containing segmented skull, brain, muscles, etc., was adopted as the “standard mesh” in both sample set collection and BLT reconstruction. Because the brain and skull of a mouse can be regarded as rigid objects, which makes multimodality imaging registration relatively easier, the standard mesh can be first registered to the structural imaging (CT/MRI) of different mice by mutual information registration ; then the surface photon acquired from the mouse by bioluminescence imaging (BLI) can be mapped to the standard mesh , and the bioluminescent source can be reconstructed in it. Details of establishing the standard mesh are provided in Supplement 1 (Section S1, Table S1, and Fig. S1). This approach makes the IPS reconstruction only dependent on a constant standard mesh. Furthermore, because glioma is a type of intracranial tumor and only invades inside the brain, the permissible region  is limited inside the brain, and the number of units in the hidden and output layer is equal to the number of nodes in the brain. The number of units in the input layer is equal to the number of nodes at the object surface. The detailed description of the network structure is shown in Supplement 1 (Section S2, Fig. S2).
The IPS network simulates the iterative inverse solution of the iterative shrinkage threshold (IST) method, which is widely used to solve the inverse problem in Eq. (1) with L1-norm [12,19–21]. In the IST method, the estimate of bioluminescent source is iteratively updated as follows:
The sample set, which contains both given surface photon and true bioluminescent source , is necessary for training connection parameters (both weight and bias) and validating IPS. However, it is extremely unpractical to establish a large population sample set by collecting a huge number of tumor-bearing mouse models with both and . This might be the major bottleneck in applying IPS. However, since the Monte Carlo (MC) method has been used as the gold standard to simulate light propagation in objects [22–24], we hypothesized that it could also be utilized to construct the training and test sets for IPS. Finally, we employed MC and collected 3200 cases of BLT simulations with 400 positions of the single bioluminescent source to build the single-source sample set. The collection was performed using MOSE 2.3 . The criteria of the sample set collection are explained in Supplement 1 (Section S1). The standard mesh (13,083 nodes and 70,297 elements) used in MC simulation was obtained from a segmented x-ray CT mouse image, which contained regions of skull, scalp, brain, and cerebrospinal fluid [Fig. 1(a)] . The corresponding optical coefficients of different organs are listed in Supplement 1 (Section S3, Table S2) .
To minimize the network overfitting, a data-assembly method with a simple computation was designed to enlarge the sample set. This method randomly selected samples from the single-source sample set. Then, it constructed the multisource sample set as follows:Supplement 1 (Section S1).
Dense connection weights and bias were learned together by minimizing the mean squared errors between reconstructed source and true source . The optimizer utilized in the network was the Adam algorithm  with learning rate: 0.0001, , and . The network was implemented by Keras 2.1.2 with Tensorflow 1.4.0 backend in Python 2.7. It was trained with , batch . All computer processing was accomplished by a personal computer with a 3.40 GHz Intel Core i7 CPU, 12 GB RAM, and GTX1080 Ti GPU.
To evaluate the performance of IPS in BLT reconstruction, we designed four experiments, including both simulation and in vivo cases. The barycenter error (BCE) between the reconstructed source and the true source was utilized to quantitatively evaluate the results in the simulation. The function of BCE is given by
Experiment 1: Fivefold cross validation was used to evaluate the positioning accuracy of IPS in the single-source BLT samples. The single-source sample set was divided into five parts (640 samples per part), and each part was in turn considered as a test set in each evaluation. The other corresponding parts and multisource sample set were used as a training set. The maximum BCE, mean BCE, and the number of test samples whose BCE is greater than 0.25 mm and 0.4 mm were also calculated to quantitatively evaluate the IPS performance.
Experiment 2: To evaluate the robustness of a network with different organ distributions, a mesh with a thicker skull and no cerebrospinal fluid was constructed [Fig. 1(a)]. Four single-source samples with different source depths were simulated in this mesh and were only used as test samples in this experiment. The single- and multisource sample sets were utilized as the training set.
Experiment 3: To evaluate the IPS performance in the dual-source reconstruction, four dual-source samples with different dual-source barycenter gaps were simulated using the MOSE platform. These samples were only used as test samples in this experiment. The single- and multisource sample sets were utilized as the training set.
Experiment 4: To evaluate the performance of the IPS method in the in vivo tumor mouse model BLT reconstruction, five orthotopic glioma mouse models were scanned by our self-developed small animal multimodality imaging system [28,29]. Their BLI, CT, and MRI images were sequentially acquired; the details are presented in Supplement 1. The reconstructed bioluminescent sources were mapped back to MRI for qualitative analysis . The ex vivo green fluorescent protein (GFP) fluorescent images together with the hematoxylin-eosin (H&E) staining of the cryo-slices were used as the gold standard to evaluate the accuracy of in vivo reconstruction.
The quantitative analysis of Experiment 1 (Table 1) proved that IPS achieved accurate positioning of the bioluminescent source. The mean and maximum BCEs of the single-source reconstruction in fivefold cross validation were only 0.078 and 0.435 mm, respectively. Furthermore, there were only 12 (1.87% in 640 test samples) IPS reconstructions per 640 tests with . These results revealed the feasibility and outstanding accuracy of IPS in single-source reconstructions. The IPS evaluation with different hidden layer numbers and different dropout probabilities is shown in Fig. S3.
Results of Experiment 2 demonstrated that IPS still achieved accurate source reconstruction when the mesh of test samples was different from the standard mesh in the training set [Figs. 1(b)–1(e)]. For different source depths, the mean and maximum BCEs were 0.057 and 0.09 mm, respectively. The detailed quantifications are provided in Supplement 1 (Section S4 and Table S3). These results revealed the good robustness of IPS for applications in nonstandardized meshes.
For dual-source reconstructions in Experiment 3, IPS showed more accurate reconstructions than FIST in all barycenter gap settings (Fig. 2). When the dual-source gap was 4 mm, FIST completely failed to distinguish the two sources, whereas IPS still offered a clear resolution. The quantitative comparisons showed that the dual-source total BCEs of IPS were about 64.6%–91.6% less than those of FIST (Fig. S4). These results demonstrated the superiority of IPS for multisource BLT reconstruction compared with FIST. Besides these results, the antinoise evaluation of IPS is shown in Supplement 1 (Section S5).
The results of in vivo BLT reconstructions in five orthotopic glioma mouse models were also consistent with all simulation experiments. Figure 3 shows the reconstruction of two mouse models as examples, and the results of the other mice are presented in Supplement 1 (Section S6, Fig. S7). The 3D rendering images [Fig. 3(a)] and tomographic slices [Figs. 3(b) and 3(c)] showed that both IPS-based BLT and MRI located the tumor in the same region inside the cranial cavity of each mouse, whereas FIST reconstructed tumor dispersedly distributed and did not overlap with any of them. Compared with ex vivo references (GFP fluorescent images and H&E staining of cryo-slices), both MRI and IPS-based BLT successfully defined the tumor location of each mouse, but FIST-based BLT failed [Fig. 3(d)]. These in vivo results again validated the accuracy of IPS for pinpointing the tumor location in BLT reconstructions.
In conclusion, we proposed a machine-learning network MLP-based IPS method to achieve accurate and quantitative BLT reconstruction. Its superiorities were demonstrated in glioma-bearing mouse models through simulated and in vivo experiments. IPS completely abandoned the forward photon propagation modeling and model-based inverse reconstruction, which makes it different from all the other existing BLT reconstruction strategies. It constructs the inverse of the photon propagation directly by learning the nonlinear mapping relationship between the surface photon density and the bioluminescent source density. Its network parameters were learned from thousands of simulated training samples obtained by the MC method. Our simulated single-source and dual-source experiments proved that IPS achieved better tumor positioning accuracy comparing with the conventional model-based FIST method. Our in vivo orthotopic glioma mouse model experiment with triple-modality imaging and ex vivo pathological references also demonstrated the feasibility and accuracy of IPS for 3D pinpointing tumors inside living animals. However, there are still obvious drawbacks of IPS, such as different training sets needing to be established for different tumor types; it cannot yet reconstruct the 3D tumor morphology beside the location; it needs an extra registration between the standard mesh and structural imaging. More discussion of limitations is presented in Supplement 1 (Section S7). Our future work will focus on solving these challenges. We believe that this novel strategy of employing machine-learning methods instead of using the conventional model-based approach for BLT reconstruction is a generalizable rationale, which holds great potential of opening a new gate for improving in vivo optical diffusion tomography.
National Natural Science Foundation of China (NSFC) (61671449, 81227901, 81527805); Ministry of Science and Technology of the People’s Republic of China (MOST) (2017YFA0205200); Key Research Projects in Frontier Science of Chinese Academy of Sciences (CAS) (QYZDJ-SSW-JSC005).
See Supplement 1 for supporting content.
1. V. Ntziachristos, J. Ripoll, L. H. V. Wang, and R. Weissleder, Nat. Biotechnol. 23, 313 (2005). [CrossRef]
2. T. F. Massoud and S. S. Gambhir, Genes Dev. 17, 545 (2003). [CrossRef]
3. G. Wang, Y. Li, and M. Jiang, Med. Phys. 31, 2289 (2004). [CrossRef]
4. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, Med. Phys. 22, 1779 (1995). [CrossRef]
5. Y. Lv, J. Tian, W. Cong, G. Wang, J. Luo, W. Yang, and H. Li, Opt. Express 14, 8211 (2006). [CrossRef]
6. K. Liu, J. Tian, C. Qin, X. Yang, S. Zhu, D. Han, and P. Wu, J. Biomed. Opt. 16, 046016 (2011). [CrossRef]
7. P. Wu, Y. Hu, K. Wang, and J. Tian, IEEE Trans. Biomed. Eng. 61, 189 (2014). [CrossRef]
8. X. Chen, Q. Zhang, D. Yang, and J. Liang, J. Appl. Phys. 115, 024702 (2014). [CrossRef]
9. M. Chehade, A. K. Srivastava, and J. W. M. Bulte, Tomography 2, 159 (2016). [CrossRef]
10. X. Zhang, Y. Lu, and T. Chan, J. Sci. Comput. 50, 519 (2012). [CrossRef]
11. X. Chen, D. Yang, Q. Zhang, and J. Liang, J. Appl. Phys. 115, 184702 (2014). [CrossRef]
12. C. Qin, S. Zhu, J. Feng, J. Zhong, X. Ma, P. Wu, and J. Tian, J. Biophoton. 4, 824 (2011). [CrossRef]
13. H. Gao and H. Zhao, Opt. Express 18, 1854 (2010). [CrossRef]
14. Y. Lu, X. Zhang, A. Douraghy, D. Stout, J. Tian, T. F. Chan, and A. F. Chatziioannou, Opt. Express 17, 8062 (2009). [CrossRef]
15. Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner, Proc. IEEE 86, 2278 (1998). [CrossRef]
16. A. Krizhevsky, I. Sutskever, and G. E. Hinton, Commun. ACM 60, 84 (2017). [CrossRef]
17. F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, and P. Suetens, IEEE Trans. Med. Imag. 16, 187 (1997). [CrossRef]
18. X. Chen, X. Gao, D. Chen, X. Ma, X. Zhao, M. Shen, X. Li, X. Qu, J. Liang, J. Ripoll, and J. Tian, Opt. Express 18, 19876 (2010). [CrossRef]
19. T. Blumensath and M. E. Davies, J. Fourier Anal. Appl. 14, 629 (2008). [CrossRef]
20. D. Han, J. Tian, S. Zhu, J. Feng, C. Qin, B. Zhang, and X. Yang, Opt. Express 18, 8630 (2010). [CrossRef]
21. S. Ahn, A. J. Chaudhari, F. Darvas, C. A. Bouman, and R. M. Leahy, Phys. Med. Biol. 53, 3921 (2008). [CrossRef]
22. Y. An, J. Liu, G. Zhang, S. Jiang, J. Ye, C. Chi, and J. Tian, IEEE Trans. Med. Imag. 36, 366 (2017). [CrossRef]
23. S. Ren, X. Chen, H. Wang, X. Qu, G. Wang, J. Liang, and J. Tian, PloS ONE 8, e61304 (2013). [CrossRef]
24. L. Wang and S. L. Jacques, Monte Carlo Modeling of Light Transport in Multi-Layered Tissues in Standard C (University of Texas, 1992), pp. 4–11.
25. D. Ancora, A. Zacharopoulos, J. Ripoll, and G. Zacharakis, IEEE Trans. Med. Imag. 36, 1086 (2017). [CrossRef]
26. G. Strangman, M. A. Franceschini, and D. A. Boas, Neuroimage 18, 865 (2003). [CrossRef]
27. D. P. Kingma and J. Ba, “Adam: a method for stochastic optimization,” arXiv:1412.6980 (2014).
28. Y. Gao, K. Wang, S. Jiang, Y. Liu, T. Ai, and J. Tian, IEEE Trans. Med. Imag. 36, 2343 (2017). [CrossRef]
29. K. Wang, C. Chi, Z. Hu, M. Liu, H. Hui, W. Shang, D. Peng, S. Zhang, J. Ye, and H. Liu, Engineering 1, 309 (2015). [CrossRef]