## Abstract

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

## Corrections

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:

where $\varphi $ and $x$ are the densities of the surface photon and bioluminescent source, respectively. The linear relationship between $\varphi $ and $x$ is modeled by the system matrix $A$. Bioluminescent source $x$ is reconstructed by solving the inverse problem of Eq. (1). However, because system matrix $A$ 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 [8] were utilized to construct the system matrix $A$. 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 (${L}_{1}$, ${L}_{p(0<p<2)}$, 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 [15]. 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 $\varphi $, and the network output is reconstructed bioluminescent source $x$. The output of layer $j$ is strongly connected to the input of the next layer $k$ and added with bias. Connection weights ${W}_{jk}$ and bias ${B}_{k}$ are learned from the network training. The activation function following the connection is called rectified linear unit (ReLU) [16] and is expressed as

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 $x$ is iteratively updated as follows:

The sample set, which contains both given surface photon $\varphi $ and true bioluminescent source ${x}_{\text{tr}}$, 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 $\varphi $ and ${x}_{\text{tr}}$. 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 [23]. 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)] [25]. The corresponding optical coefficients of different organs are listed in Supplement 1 (Section S3, Table S2) [26].

To minimize the network overfitting, a data-assembly method with a simple computation was designed to enlarge the sample set. This method randomly selected $K$ samples from the single-source sample set. Then, it constructed the multisource sample set as follows:

where ${K}_{\mathrm{set}}$ is the set of the selected samples. ${\varphi}_{i}$ and ${\varphi}_{\text{trmul}}$ are the surface photon of the $i$th single- and multisource samples, respectively. ${x}_{\text{tr}-i}$ and ${x}_{\text{trmul}}$ are the given true bioluminescent sources of the $i$th single- and multisource samples, respectively. We randomly set $K$ in the range of $2\le K\le 4$ in each construction. Then, 4500 cases of multisource samples were constructed to enlarge the sample set. The criteria for the multisource sample set collection are explained in Supplement 1 (Section S1).Dense connection weights ${W}_{jk}$ and bias ${B}_{k}$ were learned together by minimizing the mean squared errors between reconstructed source $x$ and true source ${x}_{\text{tr}}$. The optimizer utilized in the network was the Adam algorithm [27] with learning rate: 0.0001, ${\beta}_{1}\text{:}\text{\hspace{0.17em}}0.90$, and ${\beta}_{2}\text{:}\text{\hspace{0.17em}}0.99$. The network was implemented by Keras 2.1.2 with Tensorflow 1.4.0 backend in Python 2.7. It was trained with $\text{epochs}=50$, batch $\text{size}=128$. 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
$x$ and the true source ${x}_{\text{tr}}$ 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 [17]. 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.

Furthermore, the fast IST (FIST) method [20,21] with diffusion equation (DE) model [12] was utilized for comparisons in Experiment 3 and 4.

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 $\mathrm{BCE}>0.25\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$. 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.

## Funding

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.

## REFERENCES

**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]