1. Introduction

In recent years, in vivo bioluminescence imaging (BLI) has been making a difference in biomedical research, which plays a more and more important role in improving our understanding of biological processes in living subjects and accelerating drug development [1–3]. There is no inherent tissue autofluorescence generated by external excitation light, making BLI sensitive [4, 5]. Furthermore, bioluminescence tomography (BLT) can make use of the information obtained from BLI data measured on the surface of a small animal in reference to a corresponding micro-CT anatomical structure of the same small animal, and vastly improve the image quality and information available in three dimension [6].

Most near-infrared optical imaging studies rely on the development of model-based image reconstruction methods, so the forward models are essential to describe light propagation in turbid media for the inverse problems. Radiative transport equation (RTE) is generally accepted as an accurate model to govern this physical process. However, implementation of RTE is extremely complicated even in very simple settings, and imaging reconstruction with the RTE model is time-consuming [7]. Therefore, in order to make it more tractable, several approximations to RTE have been developed to overcome the difficulties for directly solving it, including the discrete ordinates (S_N) [8–11] and spherical harmonics (P_N) equations [12–15]. However, the computation cost is also very huge: ever S_N approximation needs N(N + 2) coupled equations, where N is the number of direction cosines, or the P_N approximation requires (N + 1)² coupled equations to be solved, where N is the number of Legendre polynomials. Actually, when the scattering interactions dominate over absorption, the first-order diffusion approximation (P₁) is widely utilized for modeling of light propagation in tissue and reconstructing the information of interest [6, 16–28]. However, in some particular regions including low scattering or strong absorbing tissues, turbid media with highly heterogeneous optical background, or in small geometries [9, 29], the P₁ approximation fails to work properly. In order to get over the limitations above and improve the imaging quality with less computational cost, the simplified spherical harmonics (SP_N) method has been developed [30] which requires just (N +1)/2 equations, where N is also the number of Legendre polynomials. In addition, a SP_N forward model extended to the frequency domain has been achieved [31]. SP_N-based inverse problem is also investigated for in vivo homogenous mouse BLT experiments [32]. Nevertheless, its performance in inverse problem still needs to be explored thoroughly, especially when the background of the mouse tissue is highly heterogeneous, the assumption of ${\mu \prime }_s\gg {\mu }_a$ is not applicable any more, or small tissue geometries are encountered and so forth.

Based on the P₁ diffusion approximation, the imaging priors of heterogeneous tissues in the mouse tissue have been investigated extensively in optical imaging [18, 23, 24, 27, 33–35]. In this study, furthermore, the imaging performance of SP_N model-based bioluminescence tomography method in three dimensional heterogeneous tissue is further investigated in a wide optical domain. Since SP₃ can yield transport-like solutions with lowest computational cost in biological tissue compared with SP_N (N > 3), it is actually selected as a case of SP_N here. Additionally, a generalized graph cuts optimization approach [36], which is gradient free, is incorporated with the SP₃ model to improve the reconstruction efficiency on whole body of the mouse. The paper is organized as follows: In the next section, we present the reconstruction methodology for BLT. The detailed implementation in the finite element formulation for the simplified spherical harmonics approximation is presented. In Section 3, firstly, numerical validations of the proposed method with respect to several issues discussed above are performed. Secondly, experimental mouse reconstructions in heterogeneous tissue further demonstrated the predominance of the higher-order method. Finally, in the last section we discuss relative issues and conclude the paper.

2. Methods

Based on the radiative transfer equation [7], the three dimensional simplified spherical harmonics approximation (SP_N) is obtained after a series of inferential reasoning in the planar geometry with the spherical harmonics (P_N) approximation has been carried on [30].

As above mentioned, SP₃ is applied here to evaluate the performance of this higher-order approximation model. The SP₃ equations, the boundary condition equations and the corresponding exiting partial current (J ⁺(r)) on the boundary of objects can be found in Ref. [30, 37]. According to Green formula [38], the SP₃ equations and its boundary condition equations can be deduced into the corresponding variational formulation. And then, the finite element formulation is constructed: the domain of object Ω is discretized as a tetrahedral grid $\mathcal{T}$ , so the composite moments ϕ_i(r) (i = 1, 2) and bioluminescent source distribution $\mathcal{X}$ (r) can be approximated with the piecewise linear bases:

where ϕ_i,k and x_k are the discretized values at node k when using the basis function ν_k(r); N is the total number of nodes on the entire discretized domain $\mathcal{T}$ . Now, the finite element formulation of equation is derived:

where j = 1, …, N, μ_ai(r) = μ_a(r) + μ_s(r)(1 − gⁱ) (i = 1, 2, 3), and

The detailed expression for A_i, B_i, C_i, D_i (i = 1, 2) can be found in Ref. [30]. For convenient representation, the equations above are reformulated into matrix equation:

where the corresponding components in the block matrices denote:

Then, ϕ_i can be represented by the block matrices:

After only the part on the measurable boundary in the exiting partial current equation (Eq. (6) in [32]) is retained, the linear relationship between the unknown source distribution $\mathcal{X}$ inside biological tissue and the exiting partial current J ⁺ on boundary is established, which is very similar to the formulation based on diffusion approximation [17]:

where

The common approach is to use the output-least-squares formulation incorporated with a regularization term. The solution can be determined by minimizing the following energy function:

where J ^+,meas denotes the measured exiting partial current on the boundary, ℳ the simplified spherical harmonics based forward operator and λ the regularization parameter.

Here, a gradient-free inverse solver is employed to optimize the energy function Eq. (9), which is called generalized graph cuts (GGC) algorithm [36]. Since the inverse solver is not dependent on gradient any more, it can perform very efficiently. Moreover, by restricting graph representable, the global optimal solutions can be often found in polynomial time [39], thus, this graph cuts method can provide both fast and exact reconstruction results. The detailed algorithm can be found in Ref. [36].

3. Results

3.1. Numerical validations

The heterogeneous characteristic in tissues has been suggested as a priori information to improve imaging reconstruction quality [18, 23, 24, 27, 33–35]. While the performance of simplified spherical harmonics approximation has been testified in the homogenous tissues [32], the capability of SP_N model still needs to be investigated intensively in heterogeneous media. In this part, various experiments were designed to give an insight into the proposed method from this perspective. The heterogeneous mouse atlas used for simulation validations were based on the data using the micro-CT part of dual-modality optical/micro-CT in vivo imaging system developed in our lab [40, 41] and GPU-accelerated cone-beam reconstruction algorithm [42]. The anatomical structure of mouse was segmented interactively, and several main organs were segmented (Fig. 1).

 

Fig. 1 (a) The mouse phantom scanned by the CT component in our dual-modality in vivo imaging system. (b) The heterogeneous mouse torso used for reconstructions, which contains main organs such as heart, lungs, liver, spleen and muscle. The mouse torso was discreted into volumetric mesh with 5747 nodes and 31245 tetrahedral elements with 940 nodes on the surface.

Download Full Size | PDF

It has been experimentally shown that, both the absorption and reduced scattering properties are the functions of wavelength. Based on the compositive functions and the coefficients [18, 43], the optical properties are estimated in the range of visible and near-infrared (390-1000 nm), and the ratio of ${\mu }_s^{\prime }$ to μ_a for the major organs varies over several orders of magnitude (Fig. 2). In addition, it is also observed that ratio extension is widely divergent for different organs. For instance, liver and heart cover relatively small scope, which is no more than 10, whereas bone is extended in a wide scope. On the other hand, the ratio decreases as the wavelength becomes shorter on the whole. When the condition ${\mu }_s^{\prime }\gg {\mu }_a$ does not hold, the widely used diffusion approximation becomes less valid [7, 30]. In order to investigate the imaging quality of the higher-order approximation, three typical sets of optical properties at 620 nm, 670 nm and 800 nm are considered here, which cover orange, red and near-infrared optical regimes, as listed in Table 1. Moreover, it should be noted that, for both liver and spleen contain abundant blood vessels, they are strong absorbing and low scattering compared with other tissues. Even in 800 nm, the ratio of ${\mu }_s^{\prime }/{\mu }_a$ of the two organs are still less than 5.0, and the value at 620 nm is less than 1.0. Additionally, the anisotropy parameter g is assumed as 0.9. We further stress that, the optical properties here are just applied to validate the proposed method, and in the actual situation, there may be some extent of discrepancy between calculated and the measured data from individual mouse.

 

Fig. 2 The ratio of ${\mu }_s^{\prime }/{\mu }_a$ between 390-1000 nm for the main organs in a mouse.

Download Full Size | PDF

Table 1. The optical properties for 800, 670 and 620 nm absorption of the major organs. The units for μ_a and ${\mu }_s^{\prime }$ are mm⁻¹.

View Table | View all tables in this article

Monte Carlo (MC) methods are accurate to simulate photon propagating through biological tissues. A MC-based molecular optical simulation environment (MOSE) [44] is used here to statistically estimate the bioluminescence signal on the mouse surface. In the simulation settings, a spherical solid source with the radius of 1.0 mm is located at about the centre of torso with coordinate of (21.45, 33.65, 14.50), abdomen with (18.50, 33.00, 14.50), and dorsum with (23.50, 36.50, 14.50), respectively. All three locations are in the liver. A total of 10⁷ photons at each wavelength for each source location are tracked due to the strong absorption. In all the reconstruction experiments below, we just focus on the reconstructed source distribution, and the information of quantitative value is not involved here. All the reconstructions are performed on our desktop computer with Intel Core 2 Duo 1.86GHz CPU and 3GB RAM.

3.1.1. The cases at 800 nm

Firstly, the performance of SP₃ is evaluated in relatively low absorption optical regime. Figure 3 shows the reconstructed results based on both SP₃ and diffusion approximation (DA) models. When the source is placed at the centre of the torso, the reconstruction results are shown in Figs. 3(a) and (d). It is apparent that, the results of SP₃ more accurately reflect the information for actual source not only in the position but also in the distribution than that of DA. Furthermore, as the DA-based results in Fig. 3(d), there is a large offset (3.17 mm) along the sagittal orientation, and artifacts emerge in the results, which also severely affect the reconstruction quality. One possible explanation is that, the source is adjacent to the common boundary of the liver and lung, and distance is only about 2 mm away from the boundary. Figures 3(b) and (e) plot the reconstructed results when the source is placed at the abdomen, where the tissue background in this region can be almost regarded to be homogenous. While the reconstructed location of SP₃ model deviates more than that of DA, there is any artifact around the reconstructed source, and it is not the case of DA. Furthermore, the reconstruction results for the source in the mouse dorsum are shown in Figs. 3(c) and (f), and the source is close to the vertebra. The reconstructed position offsets are 0.65 mm and 1.53 mm corresponding to SP₃ and DA, respectively (as summarized in Table 2). The reason leading to larger errors of DA than its counterpart most likely is that the diffusion approximation cannot accurately describe photon propagation when small tissue geometries is encountered, and the higher-order approximation has the potential to alleviate this effect. In this part, it is seen that, although both SP₃ and DA forward methods may model light transport in biological tissue in high scattering and low absorption optical regime, the performance of SP₃ is better than that of DA.

 

Fig. 3 The reconstruction comparisons of 800 nm cases between SP₃ and DA approximations in heterogeneous background. Figures (a)-(c) are the SP₃-based reconstruction results when the source is located at (21.45, 33.65, 14.50), (18.50, 33.00, 14.50) and (23.50, 36.50, 14.50), respectively. Figures (d)-(f) are the counterparts with DA-based reconstruction. Cross-sections with red and green boundaries are the centre position of actual and reconstructed sources, respectively. All the reconstructed values above zero are displayed in the results.

Download Full Size | PDF

Table 2. The result summary of both SP₃ and DA models for the 800, 670 and 620 nm cases. Offset denotes the distance between reconstructed and actual centre of the bioluminescent sources.

View Table | View all tables in this article

3.1.2. The cases at 670 nm

In medium absorbing optical regime, the performance of SP₃ compared with DA is also presented. Figure 4 shows the reconstruction results. Figures 4(a) and (d) display the results when the source is placed at the centre of the torso. It is very interesting that the source is reconstructed well based on two forward models. In Figs. 4(b) and (e), the results for the sources placed at the abdomen are depicted. The offset of the location centre is worse for DA compared with SP3, and reconstruction artifacts exist based on DA model again. Actually, the absorption property in liver is 0.345 mm⁻¹ and the ratio of ${\mu }_s^{\prime }/{\mu }_a$ is only 1.97, so the small mean free path should contribute to this problem. The reason can also explain the reconstruction results for the source in the mouse dorsum, as shown in Figs. 4(c) and (f). As the source is located adjacent to vertebra, the reconstructed source based on DA is leaned to be interfered by vertebra, whereas, it does not affect the results based on SP₃ much. The difference between the case in Fig. 4(f) compared with (e) is that there is no artifact around the reconstructed source. The reconstruction information is also summarized in this part in Table 2.

 

Fig. 4 The reconstruction comparisons of 670 nm cases between SP₃ and DA approximations in heterogeneous background. Figures (a)-(c) are the SP₃-based reconstruction results when the source is located at (21.45, 33.65, 14.50), (18.50, 33.00, 14.50) and (23.50, 36.50, 14.50), respectively. Figures (d)-(f) are the counterparts with DA-based reconstruction. Cross-sections with red and green boundaries are the centre position of actual and reconstructed sources, respectively. All the reconstructed values above zero are displayed in the results.

Download Full Size | PDF

3.1.3. The cases at 620 nm

As listed in Table 1, the ratios ${\mu }_s^{\prime }/{\mu }_a$ for tissues is less than 10 except the bone, so the condition ${\mu }_s^{\prime }\gg {\mu }_a$ does not hold for most of the organs. Now, we demonstrate the performance of SP₃ in this situation. Figure 5 shows the reconstructed results based on SP₃ and DA forward models. It is well known that the assumption for the diffusion approximation is somewhat not applicable in the low ratio regime [7], and the SP₃ can describe photon propagation for high absorption using higher order approximation [30, 31]. On the whole, it is very distinct that the reconstruction quality of SP₃ precedes that of DA. The reconstruction information can be found in Table 2.

 

Fig. 5 The reconstruction comparisons of 620 nm cases between SP₃ and DA approximations in heterogeneous background. Figures (a)-(c) are the SP₃-based reconstruction results when the source is located at (21.45, 33.65, 14.50), (18.50, 33.00, 14.50) and (23.50, 36.50, 14.50), respectively. Figures (d)-(f) are the counterparts with DA-based reconstruction. Cross-sections with red and green boundaries are the centre position of actual and reconstructed sources, respectively. All the reconstructed values above zero are displayed in the results.

Download Full Size | PDF

Based on the SP₃, the source in the centre can be reconstructed accurately, whereas its counterpart is affected by artifacts again, as shown in Figs. 5(a) and (d). In Figs. 5(c) and (f), the similar results exist when the source is located in the dorsum. The discrepancy of reconstruction quality between DA and SP₃ is worse than the case above. Actually, the absorption property μ_a used here is 0.542 mm⁻¹ and scattering property ${\mu }_s^{\prime }$ 0.711 mm⁻¹, so two types of properties are very close to each other. Therefore, the non-diffusive effect in the boundary of two organs gets much more conspicuous. While the diffusion approximation fails to recover the real physical process, the higher-order model can still work well. In the dorsum, the source is also reconstructed more reliable based on SP₃ than DA. Be different with the other cases above, the artifacts do not appear any more based on the DA model. The possible reason is that the region in the dorsum is less heterogeneous. It is worth highlighting that, in all the cases for the SP₃ forward model, the artefactual sources do not exist.

3.1.4. Cases with sources in other organs

The sources located in other organs are also considered to make the SP₃ more convincing for applications. In this part, for generality, the medium optical properties of 670 nm are selected. Two cases are considered: in the first case the source is located in the interspace of liver, and in the second, the source is located in the lung. These two locations represent heterogeneous effect of optical background and the anatomical structure.

Figures 6(a) and (c) display the reconstructed results for the first case. The centre offset of reconstructed source for the SP₃ is (0.70, 0.40, 0.50), and (0.90, 0.40, 0.70) for DA. The recovered source based on DA is inclined to the liver more than SP₃. Since the interspace where source is located is only about 2-3 mm wide, the accuracy of the diffusion approximation is limited in such small geometry. However, SP₃ perform better in this situation. Moreover, when the source is placed in the lung, the reconstruction results are shown in Figs. 6(b) and (d). For the relative big space and high ratio of ${\mu }_s^{\prime }/{\mu }_a$ (=11.14) in lung, two results are very similar for SP₃ and DA based models. Additionally, the reconstructed images are not affected by artifacts. The detailed information are also summarized in Table 2.

 

Fig. 6 The reconstruction comparisons of 670 nm cases between SP₃ and DA approximations in heterogeneous background. Figures (a) and (b) are the SP₃-based reconstruction results when the source is located at (19.00, 33.20, 20.70) in muscle below liver, and (22.50, 35.80, 12.00) in lung, respectively. Figures (c) and (d) are the counterparts with DA-based reconstruction. Cross-sections with red and green boundaries are the centre position of actual and reconstructed sources, respectively. All the reconstructed values above zero are displayed in the results.

Download Full Size | PDF

3.2. Mouse experiment

In this section, an in vivo heterogeneous mouse reconstruction experiment was implemented for further demonstrating the feasibility of the proposed method. The experiment was also performed on the dual-modality optical/micro-CT in vivo imaging system developed in our lab [40–42]. The optical detector is a highly sensitive CCD camera (VersArray, Princeton Instruments, Trenton, NJ) coupled with a lens (Nikkor, Nikon, Japan). The bioluminescent source was simulated by a home-made luminescent bead, which had an emission spectrum similar to that of a firefly luciferase-based source. Its dimension is about 1.5 mm in diameter and 2.5 mm long. A nude hairless mouse (Nu/Nu, Laboratory Animal Center, Peking University, China) was used in this experiment.

Before optical and X-Ray data acquisition, the CCD was cooled to −110°C using liquid nitrogen to reduce dark current noise, and the mouse was anesthetized and the bead was surgically inserted into the interspace of liver. The optical data was collected firstly. Complete angle 360°-projected optical data was acquired from the mouse surface, with 60s integration time for each image. After finishing optical acquisition, the mouse was scanned using the micro-CT to obtain the surface and anatomical structure. The CT data was segmented into heterogeneous volumetric mesh for image reconstruction, as shown in Fig. 7(a). This mesh contains 23752 tetrahedral elements and 4560 discretized nodes with 1092 nodes on the surface. The optical properties for each organ were determined with the inverse adding doubling method [45], as listed in Table 3. The estimated anisotropy parameter g for each organ is about 0.9. Since the source was located near the liver, the torso was used for imaging reconstructions, which contained over 60% volume of the mouse body. And then, the optical data was registered with the mesh for complete-angle measured data projecting from the data on CCD in two dimension to the photon distribution on the surface in three dimension (Fig. 7(b)). The source was easily distinguished in CT images and the actual position of the source could be confirmed at (24.89, 20.40, 7.88).

 

Fig. 7 The volumetric mesh of mouse in the experiment. (a) The heterogeneous mouse torso used for imaging reconstructions, including heart, lungs, liver, muscle, and bone. (b) The mapped photon distribution on the mouse surface from four views.

Download Full Size | PDF

Table 3. Optical properties for each organ in the mouse. The units are mm⁻¹.

View Table | View all tables in this article

Figure 8 shows the reconstructed results of based on both of the SP₃ and DA models. The reconstructed centre position of the source is (24.50, 21.10, 7.90) for SP₃ model, (23.90, 19.80, 8.40) for DA model. The corresponding location offsets are 0.80 mm and 1.27 mm, respectively. Although both the reconstructed sources only has a little difference from the actual position, artifacts appear near the reconstructed source in the DA-based results, whereas its counterpart not. It is similar with the results in the simulations (as in Section 3.1), and artifacts often exist in the DA-based reconstructions compared with SP₃. Moreover, in the mouse experiment, as one of the simulations with the source located in the narrow interspace of liver, the assumptions for diffusion approximation is not completely sufficed in such small geometry any more. The results of experimental reconstruction further confirm the explanation for the artifacts. In addition, the reconstruction quality in this experiment is better compared with that of the homogenous experiments [32]. The major reason leading the improvements is the incorporation of the heterogeneous anatomical structure and distribution of optical parameters. Other factors maybe complete-angle optical data registration and projection, and non-existence of bladder-like regions in the mouse body. The reconstruction time is 345.74s. The best regularization parameters are selected for the reconstruction and the parameters are in the range of 10⁻⁵ −10⁻³. Briefly, in the experimental reconstructions, it is illustrated that the proposed reconstruction method is potential for experimental reconstructions and is efficient on whole region.

 

Fig. 8 Comparisons of reconstruction results between the SP₃ and DA models. The left-bottom and right-top figures are the results of SP₃ and DA models in lateral cross sectional views, respectively, compared with the source location in the corresponding CT slices. It is noted that all the reconstructed values above zero in the slice are displayed for the results.

Download Full Size | PDF

4. Discussions and conclusion

In this study, we have investigated the imaging performance of SP₃ model-based bioluminescence tomography method in heterogeneous tissue. A gradient-free optimization approach has been applied here as the optimization approach. Compared with the diffusion approximation, the reconstruction performance of this higher-order approximation model was investigated with multiple locations in anatomically and optically heterogonous tissues, especially when the sources were located in small tissue volumes such as the interspace of liver or close to the interface of two different tissues, the higher-order model can effectively improve three dimension BLT reconstruction quality. The numerical validations demonstrated the predominance of SP₃ based model over its counterpart on the whole. In the mouse experiment, the in vivo reconstructions further illustrated the potential and computationally efficient of the proposed method in biomedical studies.

The selection of models to simulate photon transport in tissue is actually a trade-off between the precision and the computation burden. On one hand, the diffusion approximation tends to be inaccurate in low scattering and strong absorbing regions. On the other hand, S_N-based and P_N-based higher-order approximations composing too many coupled equations introduce high computational burden. Fortunately, the SP_N approximation provides a favorable opportunity to increase precision with lower computation burden. Moreover, it is observed that the SP₃ can yield transport-like solutions with lowest computational cost, especially in high-absorbing regimes for modeling light transport in tissues [30,31], and it can also provide similarly accurate imaging quality compared with higher-order (N > 3) approximations in the inverse problem [37].

The heterogeneous optical property distribution of mouse tissues is well known to be one of the most important factors that affect the quality of BLT reconstructions. The absorption properties in all tissues (Fig. 2) are calculated based on the resident oxy-haemoglobin, deoxyhaemoglobin and water concentrations [43], and due to the content tissue diversity, the corresponding parameters are very different from each other. For instance, in the same optical spectra, the absorption property of liver is much higher than the others for it is filled with massive blood vessels. Therefore, the nature of highly heterogenized coefficients makes the heterogeneous information necessary for the image reconstruction problems. Furthermore, the proposed SP₃ based reconstruction algorithm indeed demonstrates the effectiveness in this case, compared with the homogenous ones [37].

In conclusion, we have evaluated the imaging performance of SP_N (SP₃ as a typical case) model-based bioluminescence tomography in heterogeneous tissue. In various (especially high absorbing) optical regimes, numerical and experimental reconstructions have validated that the SP_N model can be beneficial to improve the tomographic bioluminescence imaging quality. In future, we will focus on real mouse experiments with tumor models and the relevant bioluminescent probes for further research.