Incorporating MRI structural information into bioluminescence tomography: system, heterogeneous reconstruction and in vivo quantiﬁcation

: Co mbining two or more imaging modalities to provide complementary information has become commonplace in clinical practice and in preclinical and basic biomedical research. By incorporating the structural information provided by computed tomography (CT) or magnetic resonance imaging (MRI), the ill poseness nature of bioluminescence tomography (BLT) can be reduced signiﬁcantly, thus improve the accuracies of reconstruction and in vivo quantiﬁcation. In this paper, we present a small animal imaging system combining multi-view and multi-spectral BLT with MRI. The independent MRI-compatible optical device is placed at the end of the clinical MRI scanner. The small animal is transferred between the light tight chamber of the optical device and the animal coil of MRI via a guide rail during the experiment. After the optical imaging and MRI scanning procedures are ﬁnished, the optical images are mapped onto the MRI surface by interactive registration between boundary of optical images and silhouette of MRI. Then, incorporating the MRI structural information, a heterogeneous reconstruction algorithm based on ﬁnite element method (FEM) with L 1 normalization is used to reconstruct the position, power and region of the light source. In order to validate the feasibility of the system, we conducted experiments of nude mice model implanted with artiﬁcial light source and quantitative analysis of tumor inoculation model with MDA-231-GFP-luc. Preliminary results suggest the feasibility and effectiveness


Introduction
No single modality has the lock on being the best molecular strategy for whole animal imaging under all circumstances [1]. The combination of two or more imaging modalities, especially providing complementary information on tissue structure and tissue-specific molecular processes, is an attractive concept for improving diagnostic specificity and patient care [2]. The combined systems have the capacity to unveil morphological, functional and molecular variations in disease and therefore enhance the diagnostic and therapeutic accuracy. The combinations of positron emission tomography (PET)/ X-ray computed tomography (CT) [3] and single photon emission computed tomography(SPECT)/CT [4] have become commonplace in clinical practice and in preclinical and basic biomedical research [5]. The combination of PET/magnetic resonance imaging (MRI) [6][7][8] has appeared in clinical trials these years and its effects are still to be tested.
In preclinical and basic research, optical techniques are prevalent due to the low cost, high throughput and high sensitivity at the level of molecular. The optical techniques are extensively used to diagnose diseases [9,10], evaluate therapies [11,12], and facilitate drug development [9,13]. However, due to the high scattering of biological tissue, it is difficult to estimate the depth of the optical signal and perform accurate quantitative analysis by using the optical signals alone. Therefore, the optical tomography techniques incorporating with a prior structural information provided by CT have been developed, such as fluorescence molecular tomography (FMT) [14], bioluminescence tomography (BLT) [15] and Cerenkov luminescence tomography (CLT) [16]. For the light source reconstruction problem, the advantage of heterogeneous reconstruction over homogeneous reconstruction has been well demonstrated [17][18][19][20][21]. In order to construct more accurate heterogeneous model, incorporating structural information from MRI is more appealing because MRI gives higher details in soft tissues than CT does. Moreover, since MRI can perform versatile analysis at the structural, functional and molecular levels because of the plenty of availability on MRI contrast agents, it is better suited than CT to molecular imaging of tissue function [22]. Stuker et al. designed a hybrid system of FMT with MRI [2] to obtain the signal of optical and MRI simultaneously. Allard et al. developed an MRcompatible platform to fix the small animal for MRI and bioluminescence imaging (BLI) using two commercial MRI and BLI devices [23]. Unlu et al. developed a combined MR-diffuse optical tomographic imaging system to report the simultaneous measurement of optical and MR contrast agent kinetics [24].
Among the various optical imaging approaches, BLT is a promising area for quantitative reconstruction of bioluminescent source distribution within a small animal from optical signals on the animals body surface [25], which can be applied to study many physiological and pathological processes in various small animal models. Since the biological tissue does not emit photons and no external light source is required for excitation, the background noise in bioluminescent signal is very low [26]. Hence, BLT has become an important tool for biomedical researchers in molecular imaging and the combination of BLT and MRI is attracting increased interest.
Generally, in optical imaging, photons emitting from the small animal are commonly detected using cooled charge-coupled detectors (CCD), which does not operate at the high magnetic field strengths of MRI scanners. In order to acquire the synchronous optical signal, fiber array and photo-multiplier tube (PMT) [2, [27][28][29][30][31][32] are often used to get deep into the magnetic field. However, they can not offer the high space resolution compared with CCD camera. Actually, is it necessary to capture the synchronous bioluminescent signal when scanning the MRI? The answer is not always positive. In most cases, the bioluminescent signal decays over about 1 hour [33]. When it has reached the steady state, it will almost has little change in the next tens of minutes. Moreover, the procedure of bioluminescence imaging takes only a few minutes, and the scanning of MRI is much longer. The two kinds of modalities is inherently not synchronous, and the biggest benefit of synchronous scan is the accurate registration between two molalities. Therefore, if the posture of the small animal keeps the same and the scanning interval between the two modalities is as short as possible, the result of the fusion of the two modalities will be meaningful.
In this paper, we present a small animal imaging system combing multi-view multi-spectral BLT with MRI.The optical device is placed at the end of the clinical MRI scanner to shorten the scanning interval with the maximum possible. The small animal is transferred between the light tight chamber of the optical device and the animal coil of MRI via a guide rail during the experiment, to ensure the same posture of the small animal between two consecutive imaging procedures. In the BLT device, two mirrors are placed beside the animal to obtain the top, left and right views simultaneously. Furthermore, multi-spectral measurements are captured by using band pass filters. After mapping the bioluminescent images onto the structural surface of MRI, the BLT reconstruction and quantitative analysis are performed using our previously developed heterogeneous reconstruction algorithm. In order to validate the feasibility of the system, we first conducted experiment of nude mice model implanted with artificial light source, and then applied the system to quantitative analysis of tumor inoculation model.
The remaining sections are organized as follows. Section 2 gives the description of the system. In Section 3, the reconstruction method of BLT is introduced as well as the evaluation benchmarks. In Section 4, numerical simulation on a heterogeneous cylindrical phantom is analyzed. Section 5 evaluates the reconstruction performance of BLT by the experiment of nude mice implanted with artificial light source. Section 6 gives the in vivo quantitative BLT experiments of tumor inoculation model. Finally, the conclusion and discussion are given in Section 7.

System design
In order to reduce the impact of the optical device to the MRI scanner, the optical device was basically constructed by non-magnetic materials. On the other hand, to protect the CCD camera from strong magnetic field, the optical device should be placed at a safe zone in which the CCD camera could work normally. Figure 1 shows the static magnetic field distribution of the MRI scanning room (GE 750w 3.0T) in the Peking Union Medical College Hospital, Beijing, China. From the field intensity map, the field intensity at the end of the patient bed is about 30 G (red region in Fig. 1). We have tested the CCD camera of our device at about 100 G, there was little effect from the magnetic field. Therefore, considering the magnetic field distribution and the convenience of operation, the optical device is placed at the end of the clinical MRI bed as shown in Fig. 2(a). After the optical imaging, the small animal is transferred from the light tight chamber of the optical device to the animal coil of MRI via a guide rail( Fig. 2(b) and (c)). Finally, the small animal coil is moved into the MRI scanner with the bed.
In order to achieve multi-view optical images simultaneously, mirrors were often adopted to acquire the photon distributions from different views [33][34][35]. Our BLT system uses two mirrors (left and right) which are canted 45 degrees and are placed beside the animal as shown in Fig. 3. Furthermore, the multi-spectral photon distributions are obtained via the filter wheel. The surface of small animal is obtained by MRI and the multi-view mirror is actually used to capture the bioluminescent information from different view angles.
The system uses a highly sensitive CCD camera (Andor iXon3 888 large field of view and megapixel back-illuminated EMCCD) which offers 1024 × 1024 pixels, 13µm pixel size, with single photon detection capability and larger than 90% quantum efficiency. The camera can be cooled down to −80 • C with air cooling (even lower temperature down to −95 • C with chiller liquid cooling). The typical CCD read noise is less than 1 e rms and the dark current is 5 ×   e/pixel/sec. The camera is coupled with a Nikon Micro-NIKKOR manual focus lens with about 30 × 30cm field of view at 60cm distance.

BLT reconstruction
The process of photon propagation in biological tissue can be described accurately by the radiative transfer equation (RTE) [36][37][38]. However, the high computational cost of RTE makes it inappropriate for BLT reconstruction directly. Various approximations of RTE, such as discrete ordinates [39,40], spherical harmonics [41], and the diffusion approximation (DA) [42] are used in practice. In the near infrared light spectrum, photon scattering usually predominates over absorption in the biological tissue. Thus, the photon propagation can be described by the following DA model [43] depending on the light wavelength λ : where x is the position vector, Ω ∈ R 3 is the problem domain, Since we perform the heterogeneous reconstruction, the small animal should be segmented into different organs and then specific optical parameters are assigned to each organs. MRI can provide high image quality of soft tissues which give more accurate segmentation results that enhances the reconstruction performance.
Because the bioluminescence imaging experiment is performed in a totally dark environment, there is no photon travels across the boundary ∂ Ω into the tissue domain Ω. The equation is subject to a Robin boundary condition [26,44]: where v(x) represents the unit outer normal on the boundary ∂ Ω and A is the boundary mismatch factor, A(x; n, n ′ ) [45], where n is the ratio of optical reflective index of the inner tissue to that outside the boundary. Then, the measured quantity is the outgoing photon density on ∂ Ω [46]: The finite element method (FEM) is a powerful tool for solving the DA equation, and has been extensively applied in BLT reconstruction in the past years [15,26,[47][48][49]. After the discretization using FEM, the diffusion equation (1) and the boundary condition (2) can be expressed as the following matrix form: where, M λ = K λ + C λ + B λ is a positive-definite matrix, and K λ , C λ , and B λ are mass, stiff and boundary matrix respectively. F λ is the source weight matrix, Φ λ and S λ are the collection of all the nodal values of the photon desity Φ λ (x) and source density S λ (x) respectively. Then, the photon density Φ λ can be obtained from equation (4): In fact, in our system, only partial photon on the boundary can be acquired for the reconstruction. Therefore, Φ meas λ is used to denote the measurable boundary data. According to the uniqueness theorem, the BLT solution is not unique in the general case [50], a strong prior knowledge must be incorporated into the reconstruction to overcome the ill-posedness effectively. The prior information such as permissible region of source should be imposed on the unknown variables to obtain a meaningful reconstruction performance. Moreover, a multi-spectral BLT reconstruction procedure is adopted to refine the ill-posedness [51]. Considering the permissible source region S p , we obtain the linear relation between the measurable photon density Φ meas λ and source density S p .
where coefficient matrix A λ = M λ −1 F λ , while its columns corresponds to the forbidden region (not in the permissible region) is removed to fit the equation.
By performing a spectral analysis, the coefficient matrix A λ and source density S p are concatenated together by using specific wights ω λ as follow: where, S t = S p λ ω λ , denotes the total photon density, and A and Φ meas is defined as: where N is the total number of spectrum involved, and the wight ω λ > 0 and meets the formula In order to obtain a stable solution, the optimization objective function is defined using L 1 regularization method: where S t sup is the upper bound of the source density, α is the regularization parameter.

Benchmarks
Mostly, the position error of the distance between the actual source position and the reconstructed source position is used to evaluate the performance of reconstruction and position error P is defined as [52]: where (x r , y r , z r ) is the center coordinate of the reconstruction source, and (x 0 , y 0 , z 0 ) is that of the actual source center. The relative error of the power density between the actual source and the reconstructed source plays an important role in quantitative analysis, since the reconstructed power reflects the concentration of luciferase molecules and luciferin. It is defined as [53]: where S r is the power of reconstruction and S 0 is the power of real source. However, just using the above two benchmarks may not reflect the real reconstruction performance. In order to give more comprehensive evaluation, we adopt the measurement of overlap which is often used for the segmentation evaluation [54]. The overlap between the reconstructed region R r and real source region R 0 is defined as: Although it is difficult to reconstruct the accurate structure of the light source, the degree of overlap still demonstrates the performance of BLT reconstruction.

Numerical simulation
The numerical simulation was performed on a heterogeneous cylindrical phantom that contains different ellipsoidal organs. The phantom was generated and decomposed into mesh by the software of COMSOL. Specifically, the cylindrical phantom had a diameter of 30mm and height of 30mm shown in Fig. 4. The light source was a cylinder with a diameter of 2mm and height of 2mm and was located at (8mm, 15mm, 0mm). The Phantom was discredited into 9490 vertex nodes and 52801 tetrahedral elements. Similar to [51], we divided the spectrum into three regions: [400nm, 530nm], [530nm, 630nm] and [630nm, 750nm]. The weights of the spectral regions were ω λ 1 = 0.29, ω λ 2 = 0.48 and ω λ 3 = 0.23. The optical parameters for the component in the phantom are shown in Table 1. In order to avoid the inverse crime problem, the synthetic data was produced by Molecular Optical Simulation Environment (MOSE) which is developed based on Monte Carlo method [55]. The source had a photon density of 300pW /mm 3 . In order to make a fair estimate of the reconstruction results with different view angles, the same permissible region of source was applied for the reconstruction. The permissible region is Ω p s = {(x, y, z)|0 < x < 15, 0 < y < 30, −7.5 < z < 7.5}.   Table 2. Obviously, the reconstruction failed at the view angles of 0 • , 60 • and 300 • , while results at other angles are somewhat acceptable. This is because the data from the view angles of 0 • , 60 • and 300 • lost too much important photon distribution information, on the other hand, the view angles of 120 • , 180 • and 240 • retain most of the photon distribution information. Therefore, when the most of photon distribution from the surface is captured, our method could achieve satisfactory reconstruction performance. Fortunately, most of the optical information could be obtained from the three view angles (top and two sides). In general, if the optical information is too weak from the top view, it is better to consider to turn the small animal over.

System calibration
In the literatures, in order to quantify the absolute energy distribution, the integrating sphere was often used for the camera calibration [18,26]. Since the power of the reconstructed source is directly proportional to the concentration of tumor cells, in this paper, we used the quantitative method in [18] (see Section 6 for details), which does not require the absolute power of the reconstructed source. Beside, the distances between the lens and animal surface should be considered. Although the distance variation is not very large in single view surface, the distance between the lens and top view is shorter than that between the lens and mirror views, which cannot be ignored. Moreover, the reflectivity of the mirrors is not 100 percent. Therefore, we measured the intensity values of a same light source from different views, and a factor of 1.3158 was calculated. The intensities of the two mirror views were multiplied by the factor for calibration purpose. Fig. 7. Overlay images of photographs and bioluminescence images.

Experiment on implanted mouse model
In order to make an accurate evaluation of BLT imaging system and its reconstruction algorithm, three nude mice were used as research objects which were labeled as LS 1 LS 2 and LS 3 . The nude mice were injected 0.2ml of anesthetic at a 0.2g/ml concentration via intraperitoneal injection. Three glass tubes filled with bioluminescent liquid were implanted into the abdomen of the three mice respectively. The bandpass filters involved were [560, 600]nm,  of intensity from the filters were less than the original intensities without filter, we normalized the weights to meet the formula ∑ N i=1 ω λ i . = 1. Thus, the weights for the filters were ω λ 1 = 0.25, ω λ 2 = 0.41 and ω λ 3 = 0.34 respectively. The optical properties of the mouse organs were calculated according to the literature [56] and the relative parameters from website http://omlc.ogi.edu/spectra/index.html, which are shown in Table 3.
The mice were placed in the the optical device for imaging (the exposure time was 20 s), and the overlay images of photographs and bioluminescent images are shown in Fig. 7. Then, the mice were transferred to the MRI scanner by the guide rail for the MRI scanning. The MRI data was T1 weighted image with TR = 9 ms and TE = 4 ms. After these data were obtained, the organs (muscle, heart, lungs, liver, stomach and kidneys) of the MRI data were segmented interactively, and were subdivided into tetrahedrons by the software of AMIRA. In order to map the bioluminescent intensities to the MRI surface, the registration between the 2D optical images and the 3D MRI surface should be conducted first. The optical boundaries of small animal were extracted by simple segmentation method. Since the posture change between the optical and MR scans is minimal, we conducted the rigid registration by calculating the differences between boundaries of optical images and silhouette of MRI. After the optical boundaries coincide with the MRI silhouette, the bioluminescent intensities were mapped to the MRI surface by the principle of small hole imaging. The permissible regions of the three mice were Ω p LS 1 = {(x, y, z)|39 < x < 43, 49 < y < 53, 11 < z < 16}, Ω p LS 2 = {(x, y, z)|43.5 < x < 47, 52 < y < 56, 11 < z < 15} and Ω p LS 3 = {(x, y, z)|44 < x < 49, 43 < y < 49.5, 11 < z < 16}. The real positions of the light source tubes were (42.6, 52.1, 12.3), (45.0, 53.5, 13.0) and (46.6, 46.5, 14.8), which were determined manually in the MRI images.
The reconstructions are shown in Fig. 8 and the benchmarks are shown in Table 4. The position error are about 1 mm for the three mice which are allowable compared with the size of the implanted source. The relative error and overlap are also acceptable for the BLT reconstruction. The results suggest the reliability of the system and reconstruction algorithm.

In vivo experiment on tumor inoculation model
In order to quantitative analyze the BLT reconstruction in real application, it is better to calculate the number of tumor cells in the lesion. Generally, the relationship between total reconstructed power and number of cells varies with different category of cells. However, for one specific tumor cell, the resultant light intensity is directly proportional to the number of luciferase molecules and the concentration of the luciferin [26]. Therefore, Liu et al. performed the linear relationship between total power and the luciferase molecular cells to quantitatively evaluate the BLT reconstruction [18]. Similarly, in this paper, we first conducted pre-experiment to find the relationship between the total reconstructed power and the number of cell MDA-231-GFP-luc. Then, the experiment on tumor inoculation model with same cell was given.

Linear relation between reconstructed power and number of cells
In this experiment, the nude mice were injected the luciferin coupled with different numbers of MDA-231-GFP-luc cells (0.55 × 10 5 , 1.10 × 10 5 , 2.19 × 10 5 , 4.38 × 10 5 ) into the hepatic lobes. Meanwhile, the control group of the same cells mixed with luciferin were injected into a 96well plate. It should be emphasized that the control group was used to evaluate the power error, which was not used for the quantification. Then, the multi-spectral images were acquired 30 minutes after inoculation of the tumor cells. The exposure time was 120 s for both the in vitro and in vivo experiments. Figure 9 shows the relationship between the total power and tumor cell of both in vitro and in vivo experiments. The total power is the reconstructed power multiplied by the reconstructed region for BLT. As expected, the intensity value of in vitro well fit a linear function (correlation coefficient R 2 = 0.9930) as well as the result of BLT (R 2 = 0.9885), and the average of relative error R is 0.3661. However, the BLI can not fit a linear function (R 2 = 0.4914) which also suggests the necessity of BLT. As shown in Fig. 9, the linear function of BLT is x = 370.3704y − 9903.3612 which will be used to evaluate the number of tumor cells in the following experiment.  the the two mice were Ω p TC 1 = {(x, y, z)|44 < x < 47, 45 < y < 48, 12 < z < 17} and Ω p TC 2 = {(x, y, z)|45 < x < 48, 52 < y < 56, 13 < z < 17}.
The overlay optical images are shown in Fig. 10. It is difficult to find the accurate locations and regions of tumor from the MRI data. The reconstructions are shown in Fig. 11, and the real location and region of the tumor might roughly be judged by the convex of the surface. The reconstructed powers and number of cells are shown in Table 5, where the number of tumor cells was calculated by the linear function x = 370.3704y − 9903.3612. The rough order of magnitudes coincide with literature report of tumor inoculation during the first week and two weeks [57].

Conclusion
We have developed a MRI-compatible optical device for BLT reconstruction. Both the simulation and in vivo experiments suggest the reliability of our system and algorithm. However, it is difficult to give an accurate evaluation of the performance of quantitative reconstruction of tumor inoculation by MRI. In our future work, PET will be involved for the evaluation. Accurate quantitative analysis is a primary advantage for bioluminescence tomography in the applications of cancer detection and monitoring, which can be facilitated by fusing MRI modality because of its better imaging quality of soft tissue. The MRI-compatible optical imaging system further promotes the development of fusing MRI and optical imaging. In addition, together with the plenty of availability on MRI contrast agents, we believe that the fusion of MRI and optical imaging has great potentials in the vast biomedical applications.