## 1.

## Introduction

Bioluminescence imaging (BLI) has attracted much attention and obtained wide applications in recent years because of its significant advantages in specificity, sensitivity, cost-effectiveness, safety and high throughput.^{1}2.3.4.^{–}^{5} Its three-dimensional (3D) tomographic partner, bioluminescence tomography (BLT) has become an invaluable tool for *in vivo* preclinical studies because of its ability to recover the localization and concentration of bioluminescent probes inside small living animals in 3D.^{6} Mathematically, BLT is an inverse source problem with the intrinsic nature of severe ill-posedness.^{7} Based on the accurate light transport model, the goal of BLT is to three-dimensionally reconstruct the distribution of the targeted probes by integrating the multiple two-dimensional bioluminescent images, geometric structures and optical properties of tissues.^{8}9.10.11.12.^{–}^{13} Of the potential applications, BLT has been successfully applied to detection and therapeutic evaluation of solid cancers,^{14}15.^{–}^{16} such as the detection of liver cancer and prostate cancer, and the prediction of prostate cancer metastasis.^{14}^{,}^{15} However, due to the inherent characteristics of cavity cancer in which the tumor cells commonly originate from the mucosa and grow outwards, the bioluminescent light will pass through a non-scattering region constructed by the pouch of the cavity organ when measured on the body surface. The existing BLT reconstruction algorithms based on the approximation model of the radiative transfer equation (RTE)^{8}9.10.11.12.^{–}^{13} are not accurate enough for such a problem because of the poor ability of the approximation models in describing light propagation in cavity organs.^{17}18.19.20.^{–}^{21}

Cavity cancers, such as gastric cancer and cervical cancer, have a relatively lower survival rate of cancer-related death. They remain difficult to cure because when using the existing techniques, they are predominantly diagnosed in the late stage. For example, gastric cancer is the second-leading cause of cancer mortality worldwide with a five-year survival rate of 20 to 30%.^{22}^{,}^{23} Accurate detection of cavity cancers becomes one of the effective approaches to decrease cancer mortality. Motivated by successful applications of the BLT technique in the detection of liver and prostate cancers, BLT may also be used for cavity cancer detection. However, a non-scattering region that commonly exists in biological tissues, such as the stomach, esophagus, intestines and bladder, has great influence on the reconstructed images.^{17} Several solutions have been proposed to handle the light propagation in the non-scattering region, such as the hybrid RTE-diffusion model (HTDM), the hybrid Monte-Carlo-diffusion model (HMDM) and the hybrid radiosity-diffusion model (HRDM).^{18}19.20.^{–}^{21} The computational complexity and expensive time-cost of the RTE and Monte Carlo method (MCM) limit the applications of the HTDM and HMDM in the inverse problem of optical imaging. In view of the significant advantage in computational efficiency, the HRDM has been acceptably applied to optical imaging.^{20}^{,}^{21} However, the application of the HRDM in existing studies is limited to the two-dimensional (2D) and 3D simple and regular geometries, such as concentric circles or homocentric spheres, to the best of our knowledge. For such concentric circles or homocentric spheres, the boundary coupling between the scattering and non-scattering region can be easily implemented. However, it is difficult for complicated geometries because of the irregularity of the boundary.

In this paper, an HRDM-based BLT reconstruction algorithm was applied to the specific problem of cavity cancer detection. In the provided algorithm, the HRDM was first applied to 3D complicated and irregular geometries and then employed as the forward light transport model to characterize the bioluminescent light propagation in both the non-scattering and high scattering region. In the HRDM, the mismatched boundary condition of the refractive index was employed to form a Neumann source at the interface of the non-scattering and scattering region, which was used to couple the radiosity theory and the diffusion equation (DE). Using the finite element technique, the HRDM was discretized and converted into a linear matrix equation that established the linear relationship between the measured bioluminescent signals and the internal bioluminescent probes. Considering the sparse nature of the bioluminescent probes and the insufficiency of the measurements, an ${l}_{1}$ norm regularization term was integrated into the objective function as the penalty term. Finally, the primal-dual interior-point method was used to solve the objective function and reconstruct the internal bioluminescent probe distribution.^{24} The performance of the algorithm was first evaluated with experiments of a concentric cylinder and a digital mouse to simulate two kinds of cavity cancers, cervical and gastric cancers respectively. To illustrate the essentiality and superiority of the HRDM in dealing with the non-scattering problem, the reconstructed results of the HRDM-based algorithm were compared with those of the DE-based one. Furthermore, a physical phantom was designed and used to illustrate the effectiveness of the HRDM-based algorithm in the real optical imaging experiment. Finally, an *in vivo* gastric cancer-bearing nude mouse experiment was conducted. The primary results revealed the ability and feasibility of the HRDM-based algorithm in the biomedical application of gastric cancer detection.

## 2.

## Methodologies

## 2.1.

### Forward Light Transport Model: Hybrid Radiosity-Diffusion Model

Because tumor cells of cavity cancers commonly originate in the mucosa of the cavity organ and grow outward, the emitted bioluminescent light goes through a non-scattering region formed by the pouch of the cavity organ when it was measured at the body surface. Thus, the forward light transport model must possess the ability to describe light propagation in both scattering and non-scattering regions. In the scattering region, the DE and Robin boundary condition are utilized to describe the bioluminescent light propagation:^{25}

## (1)

$$-\nabla \xb7[D(r)\nabla \varphi (r)]+{\mu}_{a}(r)\varphi (r)=q(r),\phantom{\rule[-0.0ex]{2em}{0.0ex}}r\in {\mathrm{\Omega}}_{\mathrm{dt}},\phantom{\rule{0ex}{0ex}}\varphi (r)+\alpha (r)D(r)[\upsilon (r)\xb7\nabla \varphi (r)]=0,\phantom{\rule[-0.0ex]{2em}{0.0ex}}r\in \partial {\mathrm{\Omega}}_{\mathrm{dt}},$$Within the non-scattering region, the propagation of bioluminescent light can be characterized by the radiosity theory as:^{26}

## (2)

$$G({r}^{\prime},r)=\zeta ({r}^{\prime},r)\frac{\mathrm{cos}\text{\hspace{0.17em}}\theta \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}{\theta}^{\prime}}{{|{r}^{\prime}-r|}^{2}}\text{\hspace{0.17em}}\mathrm{exp}(-{\mu}_{a}^{\prime}|{r}^{\prime}-r|),r,{r}^{\prime}\in B,$$When the bioluminescent light propagates across interface $B$ along the direction from the scattering to the non-scattering region, diffuse bioluminescent light transforms into non-diffuse light. Thus, the following refractive-index mismatched boundary condition is established for the transformation:

where ${r}^{\prime}$ is a point on the interface $B$; ${J}_{n}({r}^{\prime})$ is the light flux rate that points to the inside of the non-scattering region; and $\varphi ({r}^{\prime})$ is the light flux constructed by the bioluminescent light propagating in the scattering region.On the other hand, when the bioluminescent light propagates across interface $B$ along the direction from the non-scattering to the scattering region, non-diffuse bioluminescent light transforms into diffuse light. In such a case, a Neumann source is assumed to be formed at interface $B$ because of the light propagation in the non-scattering region:^{26}

## (4)

$${q}_{0}(r)={\int}_{B}\frac{1}{\pi}{J}_{n}({r}^{\prime})G({r}^{\prime},r)\mathrm{d}B,\phantom{\rule[-0.0ex]{2em}{0.0ex}}r,{r}^{\prime}\in B,$$Integrating the bioluminescent light propagation in both the scattering and non-scattering regions, a synthetic light source was constructed by adding the original bioluminescent source and the Neumann source. Substituting the original bioluminescent source $q(r)$ in Eq. (1) with the synthetic light source and putting Eqs. (2)–(4) into Eq. (1), the final integrated hybrid radiosity-diffusion model (HRDM) can be obtained as:^{27}

## (5)

$$-\nabla \xb7[D(r)\nabla \varphi (r)]+{\mu}_{a}(r)\varphi (r)\phantom{\rule{0ex}{0ex}}=q(r)+{\int}_{B}\frac{1}{\pi}\frac{\varphi ({r}^{\prime})}{\alpha ({r}^{\prime})}\zeta ({r}^{\prime},r)\frac{\mathrm{cos}\text{\hspace{0.17em}}\theta \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}{\theta}^{\prime}}{{|{r}^{\prime}-r|}^{2}}\phantom{\rule{0ex}{0ex}}\mathrm{exp}(-{\mu}_{a}^{\prime}|{r}^{\prime}-r|)\mathrm{d}B.$$## 2.2.

### Inverse Reconstruction: Finite Element Discretization and Sparse Reconstruction

Using the finite element discretization technique, the hybrid light transport model was discretized and converted into the following linear equation that links the internal bioluminescent source distribution and the nodal flux density:

where the component of $M$ can be described as:^{25}

^{,}

^{26}

## (7)

$${m}_{ij}={\int}_{{\mathbf{\Omega}}_{\mathrm{dt}}}D(r)[\nabla {\phi}_{i}(r)]\xb7[\nabla {\phi}_{j}(r)]\mathrm{d}r\phantom{\rule{0ex}{0ex}}+{\int}_{{\mathbf{\Omega}}_{\mathrm{dt}}}{\mu}_{a}(r){\phi}_{i}(r){\phi}_{j}(r)\mathrm{d}r+{\int}_{\partial {\mathbf{\Omega}}_{\mathrm{dt}}\cup B}\frac{1}{\alpha (r)}{\phi}_{i}(r){\phi}_{j}(r)\mathrm{d}r\phantom{\rule{0ex}{0ex}}-{\int}_{B}{\int}_{B}\frac{1}{\pi \alpha ({r}^{\prime})}G(r,{r}^{\prime}){\phi}_{i}(r){\phi}_{j}({r}^{\prime})\mathrm{d}B\mathrm{d}{B}^{\prime},$$Removing the rows of the inverse matrix of $M$ that corresponds to the light flux density at the internal nodes, Eq. (6) was converted into the following simplified formula:

where $A$ is the system matrix obtained by the finite element discretization.Considering the sparse assumption of the internal bioluminescent probe and the insufficient measurements, Eq. (8) was rewritten as a sparse regularization problem, in which an ${l}_{1}$ norm regularization term was utilized and integrated into the objective function:

## (9)

$$\widehat{Q}=\mathrm{arg}\underset{Q}{\mathrm{min}}\text{\hspace{0.17em}}\frac{1}{2}{\Vert AQ-{\mathrm{\Phi}}^{\text{meas}}\Vert}_{2}^{2}+\lambda {\Vert Q\Vert}_{1},$$^{24}

## 3.

## Experiments and Results

In this section, four experiments were conducted to validate the performance of the HRDM-based algorithm, including two heterogeneous simulations, a nylon phantom experiment and an *in vivo* gastric cancer-bearing nude mouse experiment. Furthermore, the DE model-based reconstruction algorithm, in which the DE model instead of the HRDM model was used to describe light propagation in tissues, was employed to illustrate the superiority of the HRDM-based algorithm. To quantitatively evaluate the reconstructed results, the distance error $d$ and the quantification error $E$ were employed:

## 3.1.

### Heterogeneous Models-Based Simulation Demonstrations

In this subsection, two heterogeneous models were designed and employed to simulate the cavity cancer models, including a homocentric cylinder and a digital mouse atlas. To avoid the inverse crime, the simulated boundary measurements were obtained by the Molecular Optical Simulation Environment (MOSE) software,^{28}29.^{–}^{30} which is based on the Monte Carlo method to simulate light propagation in tissues.

## 3.1.1.

#### Cervical cancer-mimic simulation

In this simulation, a cylindrical phantom with a radius of 15 mm and a height of 30 mm was employed, which consisted of a concentric cylinder with an inner radius of 4 mm and an inner height of 18 mm to simulate the cervix. The optical properties of the tissue around the cervix were specified as: the absorption coefficient was $0.024\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and the reduced scattering coefficient was $1.5708\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ (Ref. 31), and the cervix was considered as the non-scattering region. To simulate the internal bioluminescent probe that was marked on the cervical cancer cells, a sphere with a 1 mm radius was used and located at (8, 0, 0) mm, 3 mm away from the cervix. For simplicity, a total bioluminescent optical power of 1 nW was assumed in all simulations in this study.^{32} Thus, a power density of $0.238\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nW}/{\mathrm{mm}}^{3}$ could be obtained through dividing the unity power by the volume of the source. Figure 1 presents the reconstructed results for both the HRDM- and DE-based algorithms, where Fig. 1(a) and 1(b) shows the results of the HRDM-based algorithm, and Fig. 1(c) and 1(d) is those of the DE-based algorithm. Figure 1(a) and 1(c) is the 3D view of the results, and the corresponding axial views are shown in Fig. 1(b) and 1(d). The detailed comparison results are given in Table 1. From Fig. 1 and Table 1, we can find that both the localization and quantification results of the HRDM-based algorithm were encouraging compared with the DE model-based method, with the distance error being 0.77 mm and the quantification error being 5.57%. On the contrary, the distance error of the DE-based algorithm was 2.41 mm and the quantification error reached 50.32%. The overestimated source density was induced by neglecting the non-scattering region. The cylindrical phantom simulation results demonstrated that the HRDM-based algorithm well resolved the internal bioluminescent source in regular geometries.

## Table 1

Results of the HRDM- and DE-based algorithms in the cervical cancer-mimic simulation.

Methods | Reconstructed position | Distance error | Reconstructed density | Quantification error |
---|---|---|---|---|

HRDM | (7.60, 0.05, −0.66) | 0.77 mm | 0.225 nW/mm3 | 5.57% |

DE | (9.25, −0.73, 1.93) | 2.41 mm | 0.358 nW/mm3 | 50.32% |

## 3.1.2.

#### Gastric cancer-mimic simulation

A digital mouse atlas, extracted from CT and cryosection data,^{33} was utilized to perform the gastric cancer-mimic simulation to validate the ability of the HRDM-based algorithm in an irregular shape phantom simulation. In the simulation, only the torso section of the digital mouse was selected, and the detailed optical properties of each organ are listed in Table 2.^{34} In this study, the stomach was modeled as an air bubble, termed as the gastric pouch, which was used to simulate the cavity organ. Because it is difficult to localize the stomach wall during organ segmentation and impossible to discretize it during finite element discretization, the stomach wall was not modeled. An ellipsoid with the size of (1, 1.5, 1) mm was positioned at (26, 12, 23) mm to mimic the inner bioluminescent probe that has been marked on the gastric cancer cells, as shown in Fig. 2(a). Similar to the cervical cancer-mimic simulation, the power of the source was set to be the 1 nW for simplicity, and the power density was calculated as $0.159\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nW}/{\mathrm{mm}}^{3}$. The measured light flux distribution on the surface of the digital mouse was simulated by MOSE, as shown in Fig. 2(b). The reconstructed results for both the HRDM-based and the DE-based algorithms were obtained and are presented in Fig. 3, where Fig. 3(a)–3(c) presents the coronal, axial and sagittal view results of the HRDM based algorithm, and Fig. 3(d)–3(f) shows those for the DE-based algorithm. The pink sphere represents the actual source, and the colored domains around it are the reconstructed sources. The detailed quantitative comparison results are listed in Table 3. Similarly, more accurate and acceptable results were obtained by the HRDM-based algorithm than the DE method. The distance error of the HRDM was 0.63 mm, which was smaller than the size of the actual source and much better than that of the DE (3.05 mm). On the other hand, a more accurate quantification result was obtained by the HRDM, in which the reconstructed source power density was $0.149\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nW}/{\mathrm{mm}}^{3}$ with the quantification error of 5.88%. The reconstructed source power density of the DE was overestimated (greater than 150%) because of neglecting the non-scattering region. The comparison results demonstrated the accuracy of the HRDM-based algorithm in the irregular shape phantom, and also revealed its superiority to the DE-based algorithm when the non-scattering region cannot be neglected.

## Table 2

Optical properties of each organ for the gastric cancer-mimic simulation (in units of mm−1).

Muscle | Heart | Gastric pouch | Liver | Kidneys | Lungs | |
---|---|---|---|---|---|---|

μa | 0.23 | 0.11 | 0(0.21) | 0.45 | 0.12 | 0.35 |

μs′ | 1.00 | 1.10 | 0(1.70) | 2.00 | 1.20 | 2.30 |

Note: The values listed in parentheses are the parameters used in the DE model-based algorithm.

## Table 3

Results of the HRDM- and DE-based algorithms in the gastric cancer-mimic simulation.

Methods | Reconstructed position | Distance error | Reconstructed density | Quantification error |
---|---|---|---|---|

HRDM | (25.97, 11.43, 22.74) | 0.63 mm | 0.149 nW/mm3 | 5.88% |

DE | (28.46, 10.71, 21.76) | 3.05 mm | 0.409 nW/mm3 | 157.05% |

## 3.2.

### Physical Phantom-Based Experimental Demonstration

A cylindrical phantom made of nylon was designed and employed to perform the imaging experiment. The dimensions of the phantom were 30 mm in diameter and 30 mm in height. Two different size holes were drilled into the phantom to hold the light source and to simulate the cavity organ. The No. 1 hole was employed to hold the luminescent liquid extracted from a luminescent light stick (Glowproducts), with a radius of 1 mm and a depth of 16 mm; the No. 2 hole had a 3 mm radius and 18 mm depth and was used to simulate the cavity organ. The detailed schematic diagram is shown in Fig. 4. A luminescent light source with a 2 mm height and a 1 mm radius was formed by injecting $6.28\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ luminescent liquid into the No. 1 hole, with its center at (3, 0, 0) mm and a power density of $0.224\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nW}/{\mathrm{mm}}^{3}$.^{25} At the same time, a cavity organ with a 13 mm height and a 3 mm radius was constructed by blocking the top of the No. 2 hole with a nylon rod of 5 mm in length. The central wavelength of the luminescent light source was about 650 nm, so the measured optical properties of the nylon phantom at 660 nm were used: the absorption coefficient was $0.0138\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and the reduced scattering coefficient was $0.91\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$.^{35}

After the phantom was mounted on the rotation stage and rotated a complete 360 deg, four equally spaced bioluminescent images and the reference white-light images were acquired using a dual-modality molecular imaging prototype system of ZKKS-Direct 3D (jointly developed by Guangzhou Zhongke Kaisheng Medical Technology Co., Ltd. and Xidian University), and mapped onto the surface of the phantom using the surface light flux reconstruction method.^{36} Based on the mapped surface light flux distribution, the luminescent light source was reconstructed using the HRDM-based and the DE-based algorithms, as shown in Fig. 5. Therein, Fig. 5(a)–5(c) presents the coronal, axial and sagittal view results of the HRDM-based algorithm, and Fig. 5(d)–5(f) shows those of the DE-based algorithm. The detailed quantitative comparison results are listed in Table 4. From Fig. 5 and Table 4, we find that a similar conclusion can be addressed as the one obtained for the simulation demonstrations. More accurate reconstructed results were obtained by the HRDM-based algorithm than those of the DE-based one, as shown in Table 4. However, the superiority of the HRDM-based algorithm over the DE one was not as obvious as the one presented in the simulation demonstrations. This may be caused by the unpredictable factors existing in the imaging experiment, such as the obliquity of the phantom, the errors caused by the registration and mapping of the 2D bioluminescent images onto the 3D phantom surface, and the inaccuracy of the optical parameters. Overall, the cylindrical phantom experimental results demonstrated the acceptable accuracy of the HRDM-based algorithm and its superiority to the DE-based method.

## Table 4

Results of the HRDM- and DE-based algorithms in the physical phantom experiment.

Methods | Reconstructed position | Distance error | Reconstructed density | Quantification error |
---|---|---|---|---|

HRDM | (3.82, −1.01, −0.99) | 1.63 mm | 0.218 nW/mm3 | 2.61% |

DE | (1.33, 1.18, −1.51) | 2.54 mm | 0.261 nW/mm3 | 16.34% |

## 3.3.

### In vivo Gastric Cancer-Bearing Mouse Experiment

An *in vivo* gastric cancer-bearing mouse experiment was conducted to evaluate the performance of the HRDM-based algorithm. A female athymic BALB/c nude mouse about six weeks old was employed to perform the gastric cancer detection. Prior to the experiment, the mouse was anesthetized and the left upper abdomen of the mouse was carefully opened with a 2 to 3 cm incision to expose the stomach of the mouse. $5\times {10}^{6}$ SGC7901-Luc-GFP cells were then inoculated into the subserosa layer of the gastric wall. Finally, the abdominal wall and skin were sutured with prolene.

Twenty-one days after the SGC7901-Luc-GFP cells were inoculated, the D-luciferin solution ($150\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mg}/\mathrm{kg}$ body weight) was injected into the peritoneal cavity of the mouse. Four equally spaced bioluminescent images and the reference white-light images were acquired using the dual-modality molecular imaging prototype system and were mapped onto the surface of the mouse using the surface light flux reconstruction method.^{35} The anatomical structure of the mouse was also obtained using the dual-modality molecular imaging prototype system and the corresponding optical properties are calculated in Table 5.^{37} Based on the mapped surface light flux distribution and the acquired anatomical structure, the distribution of the inoculated SGC7901-Luc-GFP cells was reconstructed using the HRDM- and DE model-based algorithms respectively. Figure 6 shows the detection results, where Fig. 6(a) and 6(b) presents results of the HRDM-based algorithm, Fig. 6(c) and 6(d) shows those of the DE model-based algorithm, Fig. 6(a) and 6(c) is the 3D views of the results, Fig. 6(b) and 6(d) is the enlarged views of the stomach, and Fig. 6(e) presents the necropsy observation of the ill-stomach.^{38} Considering the gold-standard features of the necropsy observation, we qualitatively compared the reconstructed distribution of the SGC7901-Luc-GFP cells with the necropsy observation results. From Fig. 6, we intuitively find that the reconstructed distribution of the SGC7901-Luc-GFP cells using the HRDM-based algorithm was more consistent with the necropsy observation than that obtained by the DE-based algorithm. The *in vivo* gastric cancer-bearing mouse experiment illustrated that the HRDM-based algorithm has great potential for the application of gastric cancer detection.

## Table 5

Optical properties of each organ for the in vivo gastric cancer-bearing mouse experiment (in units of mm−1).

Muscle | Heart | Gastric pouch | Liver | Kidneys | Lungs | |
---|---|---|---|---|---|---|

μa | 0.0086 | 0.1382 | 0(0.0263) | 0.8291 | 0.1550 | 0.4596 |

μs′ | 1.2584 | 1.0769 | 0(1.5492) | 0.7356 | 2.5329 | 2.2651 |

Note: The values listed in parentheses are the parameters used in the DE model-based algorithm.

## 4.

## Discussions and Conclusion

We have provided an HRDM-based reconstruction algorithm for the application of BLT in the detection of cavity cancer, in which a hybrid light transport model and the sparse regularization technique were incorporated. In the algorithm, the HRDM was first applied to 3D complicated and irregular geometries and then utilized as the forward light transport model to describe the bioluminescent light propagation in tissues. Reconstructed results for both the different shape geometries-based simulations and the physical phantom-based experiment have been illustrated. Both the localization and quantification results demonstrated the accuracy and effectiveness of the HRDM-based algorithm and revealed its essentiality and superiority by comparing it with the DE model-based algorithm. Furthermore, an *in vivo* gastric cancer-bearing mouse-based experiment was conducted, and the preliminary results demonstrated the feasibility and potential of the HRDM-based algorithm in the application of gastric cancer detection. However, further applications of the HRDM-based algorithm in the whole-body small animal imaging are restricted because the DE is inaccurate or invalid in the low-scattering region or when the light source is close to the boundary.

To quantitatively illustrate the influence of the position of the light source on the results of the HRDM-based algorithm, a heterogeneous geometry-based simulation was performed by changing the position of the light source. In the simulation, the same size of the geometry as that used in the cervical cancer-mimic simulation was employed. At the same time, the same bioluminescent source was set at (5.5, 0, 0) mm, with a boundary-to-boundary distance between the source and cavity organ being 0.5 mm, which was smaller than one mean free path (0.627 mm). The reconstructed results of the HRDM-based algorithm for the simulation are shown in Fig. 7, where Fig. 7(a) presents the coronal view of the results and Fig. 7(b) is the axial view. We found that poorly reconstructed results were obtained, with the localization error being 2.21 mm and the quantification error being 70.88%. The results revealed the limitations of the HRDM-based reconstruction algorithm, which would be inaccurate when the light source is close to the boundary. To overcome the limitations, high order approximation models of the RTE, such as the simplified spherical harmonics approximation (${\mathrm{SP}}_{N}$) model, should be utilized to replace the DE in the forward hybrid light transport model. The incorporation of the high order approximation model can also clear up other limitations of the HRDM, such as the dependency of the DE on the high scattering property.

In conclusion, the numerical simulation, phantom and *in vivo* experimental results demonstrated that the HRDM-based reconstruction algorithm would broaden the application of BLT technology and may provide great potential for the applications of detecting cavity cancer. Our further study will focus on the improvement of the HRDM-based reconstruction algorithm for the application of whole-body imaging.

## Acknowledgments

This work is supported by the Program of the National Basic Research and Development Program of China (973) under Grant No. 2011CB707702, the National Natural Science Foundation of China under Grant Nos. 81090272, 81101083, 81101084, 81101100, 81000632, 30900334, the National Key Technology Support Program under Grant No. 2012BAI23B06, and the Fundamental Research Funds for the Central Universities.