Fluorescence molecular imaging is a very popular method to study biological molecules, pathways, and events in living cells and tissues.¹ With its sensitivity, specificity, noninvasiveness, and cost-effectiveness, fluorescence molecular imaging is widely applied in preclinical and clinical applications. However, after propagating over a few hundred micrometers, fluorescence light would become highly diffusive. The state-of-the-art microscope operates only at depths less than one transport mean free path length in tissue. With the intense scattering signal, it becomes difficult to image a fluorescence probe distribution within a thick biological sample using the current fluorescence microscopy.^{2, 3}

In this letter, we propose a fluorescence tomographic method to reconstruct an optical molecular probe distribution using a standard fluorescence microscope on a thick biological sample. The light propagation through the biological tissue is a complicated process that involves both absorption and scattering. The photon propagation model governs the interaction between light and tissues, and is essential for tomographic imaging with visible and near-IR light.

The popular diffusion approximation model (DA) is fairly accurate and efficient in the cases of high scattering and weak absorption in a large tissue sample. However, the DA model suffers from a substantial discrepancy near either sources or boundaries,⁴ which happens to be the case of microscopic imaging. Monte Carlo simulation (MC) is a statistical technique that accurately traces photon paths in the biological environment. The MC method can give an accurate estimation for light propagation.⁵ However, the MC method takes long computation time. Recently, the massive parallel approach using general-purpose graphic processing units (GPGPU) has been adopted to speed up the MC simulation. Alerstam showed that GPU-based acceleration could be applied in the MC simulation to obtain a massive speedup over the traditional CPU implementation.⁶ Recently, we implemented a GPU-based Monte Carlo simulator on a multi-GPU system with the CUDA platform. Using four NVidia Tesla C1050 GPUs, we are able to archive a speedup of over 800 times compared to the same implementation on a single CPU core, enabling the Monte Carlo method to be used in the inverse process.

Using the GPU-based Monte Carlo simulation, we developed a Monte Carlo algorithm for fluorescence tomography. Usually, microscopic imaging focuses on a region of interest (ROI) in the biological sample and acquires surface fluorescence scattering signals emitted from optical molecular probes in the ROI. The ROI volume can be discretized into finite volumetric elements. The acquired photon fluence rate data on the surface of the ROI is a linear combination of fluorescence source intensities,

where Φ is a vector of the photon fluence rate measured at detector pixels, S is a vector of fluorescence source intensities, and A is the discrete Green function. A column vector of A is constructed from the photon fluence rate at each detector pixel computed from the GPU-based Monte Carlo simulator, assuming that each potential elemental fluorescence source is with a unit intensity.

In practice, the measurable quantity Φ in Eq. 1 can only be collected on a partial external surface of an object, either upside or downside of a biological sample, while the fluorescence sources are distributed in a three-dimensional (3D) domain. Clearly, the source reconstruction is an ill-conditioned problem, which often produces false solutions due to the existence of multiple solutions and the presence of measurement noise. Regularization is often used to stabilize the reconstruction. With the use of regularization, such as Tikhonov regularization, feasible region constraint,⁷ and compressive sensing-based l ₁-norm regularization,^{8, 9} fluorescence tomography shows improvements, but the ill-posedness of the problem remains.

To find a stable solution, here we introduce a well-posed model to regularize the inverse source problem: Given the prior optical parameters of the sample, including the absorption coefficient μ_a, scattering coefficient μ_s, and anisotropic factor g, then the positions and powers of luminescent sources can be uniquely determined from the measurable diffused light signal on the boundary of the ROI under the assumption that the number of isolated luminescent sources n is known in the ROI. Based on this concept, Eq. 1 can be converted to an optimization problem:¹⁰

where the vector S only includes n positive nonzero values [TeX:] $\left\{ {S_{i_1 },S_{i_2 }, \cdots,S_{i_n } } \right\}$ $\left\{S_{i_1},S_{i_2},\cdots ,S_{i_n}\right\}$ , which represents the power of the fluorescence sources, and its index {i ₁, i ₂, ⋅⋅⋅, i _n} corresponds to the source positions. The new reconstruction model based on Eq. 2 has a unique solution,¹¹ and can effectively overcome the ill-posedness of fluorescence tomography.

Equation 2 is an optimization problem subject to the l ₀-regulization, which is non-deterministic polynomial time hard (NP-hard), and cannot be efficiently solved using a commonly adopted gradient-based optimization method. In the work, we developed a differential evolution (DE)-based heuristic reconstruction algorithm to solve the optimization problem.¹⁰ DE is based on a spontaneous self-adaptive vector difference operator, and displayed a superior converging behavior with high precision. Experience with numerous benchmark and real-life problems suggests that the DE often has a fast and reliable converging behavior with a relatively small population size. Each candidate solution in DE is a weighted combination of the columns in A. Due to the rapid convergence of DE, the optimizer can delineate the feasible region of the fluorescence sources in a few generations (i.e., iterations), and many of the columns in the discrete Green function A are never used. Taking advantage of this observation, we replace the full Green function construction with an on-the-fly evaluation scheme: whenever a column in A is needed in a candidate solution, and has not been evaluated, that column will be computed by performing the Monte Carlo simulation. Compared with the full Green function evaluation, the on-the-fly scheme saves about 40 to 70% of Monte Carlo runs, which further speeds up the reconstruction process.

The constraint on the number of light sources, however, cannot be always obtained in advance. We can heuristically determine the number of light sources by gradually increasing the number of sources and solve the optimization problem in Eq. 2. For each reconstructed source distribution, if two or more sources are so close that they cannot be isolated within a user-provided spatial resolution, they are combined into a single source. After the source combination, the source distribution is equivalent to a previous reconstructed counterpart with a smaller number of sources. Such a number can always be iteratively and heuristically determined to interpret measured data effectively and efficiently.¹⁰

We investigated the accuracy and robustness of the proposed fluorescence microtomography approach in experiments on biomimetic tissue scaffolds. Three poly (D,L-lactide) (PDLA) scaffolds were fabricated using an electrospinning technique¹² with thicknesses of 0.65, 0.8, and 1.4 mm, respectively. The optical properties were measured using a spectrophotometer (Varian Cary 5000 UV-Vis-NIR, GMI, Inc.). The absorption coefficient is approximately 0.01 mm⁻¹, the scattering coefficient 15.0 mm⁻¹, and the anisotropic factor 0.9. Fluorescence microbeads of 90 μm diameter (Spherotech Inc.) were laminated and pressed onto each scaffold sample to mimic green fluorescence labeled cells. The actual microbead distribution was captured under a fluorescence microscope, as shown in Fig. 1, and serves as the control image. The scaffold was illuminated using a uniform field of blue laser light (λ = 473 nm) to excite the microbeads, and the emission signal captured using an EMCCD camera with a 1x long working distance lens (Mitutoyo, M Plan 1×), creating a field of view of 8.192 mm × 8.192 mm, which is the same as the size of the EMCCD sensor. From each set of captured emission signals, we reconstructed the microbead distribution, as shown in Fig. 2. The on-fly matrix evaluation and optimization took about 20 minus. The results show no significant difference in localization accuracy over all the scaffold thicknesses. Figures 2b, 2e, 2h present maximum position errors of ∼250, ∼275, and ∼500 μm respectively. The theoretical spatial resolution is finer based on the numerical studies, which is around 100 to 150 μm (∼2 to 3 pixels) in the lateral and axial directions with the current imaging setup. These reconstruction errors were relative to the size of ROI, and induced from the measurement noise rather than the ill-posedness of the reconstruction. This claim was confirmed in our numerical simulation where a reduced measurement noise was correlated to an improved location accuracy.

In summary, we have developed a novel fluorescence microtomography approach to reconstruct fluorescence reporters in thick biological samples. This tomographic imaging method is based on a well-posed inverse model we have established, and has a unique and stable solution. Our experiments have demonstrated that the proposed fluorescence microtomographic method is able to localize and quantify an optical molecular probe distribution accurately, up to a spatial resolution of several hundred micrometers in relatively thick and highly opaque samples. Using a powerful microscope setup, we expect further improvements in spatial resolution for fluorescence microtomography.