High-resolution reconstruction of fluorescent inclusions in mouse thorax using anatomically guided sampling and parallel Monte Carlo computing.

We present a method for high-resolution reconstruction of fluorescent images of the mouse thorax. It features an anatomically guided sampling method to retrospectively eliminate problematic data and a parallel Monte Carlo software package to compute the Jacobian matrix for the inverse problem. The proposed method was capable of resolving microliter-sized femtomole amount of quantum dot inclusions closely located in the middle of the mouse thorax. The reconstruction was verified against co-registered micro-CT data. Using the proposed method, the new system achieved significantly higher resolution and sensitivity compared to our previous system consisting of the same hardware. This method can be applied to any system utilizing similar imaging principles to improve imaging performance.

1. Introduction

Fluorescence diffuse optical tomography (FDOT) is a noninvasive method for imaging fluorescent agents in living organisms using far-red or near infrared light. With the rapid development of fluorescently labeled molecules as the targeted imaging probes, FDOT promises great potential for preclinical studies [1-3]. Compared to other imaging modalities, such as magnetic resonance imaging (MRI), x-ray computed tomography (CT), and ultrasound, a major disadvantage of FDOT for preclinical imaging is the difficulty in modeling light propagation in highly scattering and heterogeneous media. Many techniques have been proposed to overcome this limitation by addressing one or more of the difficulties in modeling, e.g., the use of Born normalization [4], higher-order harmonics for the source term [5], mathematical approximation using simplified geometry [6], altered geometry of the animal using a compressing imaging chamber [7], and the time-gated early-photons [8]. Nonetheless, it is still a challenge for FDOT to resolve closely located small fluorescent heterogeneities in the torso of small animals, which is critical for effective preclinical imaging. The thorax of the small animal, in particular, is of great interest for preclinical applications because it contains the heart and lungs that are very important target organs for a large number of cardiovascular, pulmonary, and cancer disease models, which play critical roles in numerous preclinical and translational efforts.

Light propagation in the turbid media, also known as the “forward problem,” has been reported extensively in the literature [9-11]. The main category of the forward modeling methods is based on the radiative transfer equation (RTE). Due to the complexity of the RTE, the diffusion approximation is often applied in the context of diffuse optical imaging, including diffuse optical tomography (DOT) and FDOT. In practice, further simplification of the diffusion equation, most notably using the Born approximations, are frequently used to solve the forward problem in diffuse optical imaging [4]. Another important category of the forward models is based on the Monte Carlo method, which is known for its accuracy but also for its high computational cost [12-14]. Reconstruction of the FDOT image from a series of optical measurements and the knowledge of the optical medium, as well as the boundary conditions, is referred to as the “inverse problem.”

Successfully solving the inverse problem relies on solving the forward problem correctly, because the two problems share one component—light propagation in the medium. The specific choices of the forward and the inverse problem solvers define a particular reconstruction method for diffuse optical imaging. In preclinical imaging, the highly complex geometry typically requires the use of either the finite-element or Monte Carlo methods. On one hand, although the Monte Carlo method is computationally expensive compared to the finite-element method, it relies on fewer assumptions. For example, the diffusion approximation, on which most of the finite-element method is typically based, becomes invalid when the ratio of the absorption to the scattering coefficients is low (e.g., in the liver, blood, or bladder) or, when the source/detector distance is comparable to the transport mean free pathlength. In addition, the simplicity of the Monte Carlo method makes it attractive over the finite-element method, which involves complex boundary conditions. On the other hand, as the computer technology continuously improves, particularly the wide availability of parallel computing, the Monte Carlo method is becoming more and more attractive. In this study, the forward problem solver used a supercomputer-based parallel implementation of our Monte Carlo method [15]; and the inverse problem solver was based on the Born method.

Our goal was to develop new data acquisition and image reconstruction techniques that can provide high spatial resolution and high sensitivity to satisfy the needs of preclinical applications. These techniques included (1) parallel implementation of our Monte Carlo method that constructed the Jacobian matrix directly based on photon paths from the source to the detector, rather than taking the product of the two fields of photon density waves originated from the source and the detector separately; and (2) use of anatomy of the animal to automatically identify and eliminate problematic data points within a large data set.

In diffuse optical imaging, the inverse problem is typically underdetermined in that there are usually fewer measurements than the number of unknowns. Therefore, for high-resolution reconstruction, denser spatial sampling is necessary. Using the free-space full-angle acquisition method [16], the inverse problem can become seemingly overdetermined using a scanning laser and a charge-coupled device (CCD) camera. Although studies show that higher sampling density can improve spatial resolution, e.g., in [17,18], sampling strategy plays an important role for efficient and high-resolution reconstruction [19-21]. Particularly in free-space animal imaging, where the exact surface of the animal is highly irregular and unpredictable, a mechanism is needed to accurately determine source/detector positions using an optical method during the experiment. We developed a geometrical calibration method to optically determine the source and detector positions and to subsequently co-register the optical and x-ray data.

For diffuse optical imaging, the highly heterogeneous macroscopic anatomical structures of the animal significantly alter the propagation of light and thus heavily influence the detected signal. Consequently, incorporation of structural information is critical for correctly modeling light propagation and improving reconstruction in animal experiments, as concluded in many studies such as [22-26]. In this study, we not only used the x-ray micro-CT data registered to the optical data as the source of structural priors for FDOT reconstruction, but also used this information for optical data sampling. Specifically, instead of sampling uniformly, or using any predetermined pattern, we used an anatomically guided sampling strategy to exclude data that would otherwise inversely affect the reconstruction.

2. Methods

2.1. Animal preparation

We chose quantum dots (QD) as the fluorescent contrast agent because of their high quantum yield, excellent photostability, and broad selection of functional groups for targeted imaging in the future. Although the cytotoxicity of cadmium-based QD to living organisms is well-known, the use of long molecular chain polymer coatings, particularly polyethylene glycol, can effectively eliminate the toxicity concerns of the QD for in vivo studies, e.g., as reported in [27-29]. We used commercially available non-conjugated QD branded Qtracker 800 (Invitrogen, Carlsbad, CA). At our excitation wavelength of 642 nm, the extinction coefficient of the QD is 750,000 M−1cm−1; and the quantum yield is 3%, according to the product specifications.

The animal experiment was conducted on a nude mouse (Nu/Nu, male, 25.5 g) cadaver to limit the complexity of animal support, while still providing a biologically realistic model. The QD inclusions were constructed by sealing polyethylene capillaries with paraffin and filled with 1 μL of 200 nmol/L QD saline solution, which were inserted into the thorax via the esophagus and the trachea from an incision on its neck. Each of these inclusions contained 200 femtomoles of QD. The dimensions of the inclusions were 0.9 mm in a diameter and 1.6 mm in length. The distance measured between their closest edges was 2.5 mm. The dimensions and the positions of the inclusions were verified by x-ray micro-CT. These inclusions were enclosed by the heart, the lungs, and the spine, representing a very challenge but nonetheless realistic situation for FDOT. The purpose of such a design was to emulate a reasonably favorable scenario for high-resolution FDOT study of the mouse with lung tumors: the dimensions of the inclusion are equivalent to a millimeter-sized tumor; and the same concentration of the inclusion is achievable in the blood of the mouse (a typical total volume of 2 mL) by injecting a reasonably large dose of QD (200 μL of 2 μM) via the tail vein.

2.2. Experimental setup

Our optical imaging platform consists of a continuous-wave laser diode as the light source, a galvo-scanner that enables raster scanning of the laser on the animal, a rotation stage that vertically suspends the animal and provides unobstructed panoramic view of the animal, a CCD camera (Imager QE, LaVision, Germany) as a high-density photodetector array to measure the transillumination, and a filter wheel to select the appropriate detection wavelength for the camera. The control and synchronization of these devices were implemented using lab-developed LabVIEW (National Instruments, Austin, TX) applications. This experimental setup is similar to the one previously described in [25], except that the diode laser and the filters were replaced to match the excitation and emission spectra of the QD. The laser diode (HL6388, Opnext, Fremont, CA) radiates at 642 nm with an optical output power of 250 mW. The absorption spectrum of the QD monotonically decreases from ~300 nm in the ultraviolet band and extends beyond 700 nm in the near infrared band. The fluorescent emission spectrum of the QD peaks at ~785 nm with a full-width at half-maximum bandwidth of ~75 nm. Our choice of this particular excitation wavelength was based on considerations of optical penetration, the excitation efficiency, and the specifications of commercially available laser diodes.

The experiment started with fluorescence acquisition, using two stacked 750 nm long-pass filters (CVI Melles Griot, Albuquerque, NM) to collect the fluorescent light and a 647 nm long-pass edge filter (Semrock, Rochester, NY) to block the laser. The CCD camera acquired one image, with integration time of 4 s, for every source position on the animal until the galvo-scanner completed a scan. The rotation stage subsequently changed the angular position of the animal, before the next laser scan and corresponding data acquisitions began. This was repeated until the rotation stage resumed its initial angular position. The same procedure was subsequently applied to acquire the excitation images using a neutral density filter to avoid CCD saturation (integration time 0.1 s). We used 6 vertically arranged source positions for each angular position, each separated by ~2 mm and 40 angular positions, which resulted in 240 images for both fluorescence and excitation. We applied time delays of 0.1 s between consecutive source points and 1 s between angular positions. Overall, fluorescence acquisition required 1024 s; acquisition of excitation images needed 88 s; and the geometrical calibration took 1.2 s.

After the acquisition of fluorescence and excitation images, a set of reflectance images were taken for the same 40 angular positions using a diffusive LED white-light source (Thorlabs, Newton, NJ). These reflectance images were used to extract the 3-D surface of the animal using a parallel-beam backprojection algorithm. Details of this method were previously described in [25].

The last step in optical imaging was to acquire calibration images for precise spatial registration. A calibration target was constructed using a rigid translucent sheet material, on which a rectangular grid with known intervals was printed. The calibration target replaced the rotation stage and was positioned coplanar with the axis of rotation. The laser impinging on the target produced a bright spot. Under front illumination, both the printed grids and the laser spot were visible to the CCD camera. One CCD image was taken for each of the 6 laser/galvo-scanner positions.

After optical acquisition, the rotation stage (with the animal undisturbed) was transferred to the adjacent x-ray system to acquire micro-CT data as structural information, which was then co-registered with the optical data using our calibration method. Details of our dual x-ray tube/detector micro-CT system were previously described in [30]. A set of 200 projections were acquired in less than 3 minutes using a single tube/detector with an angular step size of 1.8° per projection. The x-ray tube operated at 80 kVp, 160 mA, and 10 ms per exposure. The radiation dose associated with the micro-CT scan was 7.2 cGy. The projection set was used to reconstruct the tomograms with a Feldkamp algorithm [31] implemented in the Cobra EXXIM software package (EXXIM Computing, Livermore, CA). The micro-CT image was reconstructed as a (512)3 matrix with an isotropic voxel size of 88 μm.

2.3. Data processing

Based on the white-light reflectance CCD images, a sampling region of 9 × 18 mm was defined in the chest area of the animal in the emission and excitation CCD images. Within the sampling region, 4 × 10 sampling points were evenly spaced. The resulting sampling density was 3 × 2 mm horizontally and vertically, respectively. The mean value of the pixels within a 1 × 1 mm area centered at the sampling point was used as the optical signal intensity for the sampling point. The spatial sampling points were identical for the raw CCD images at the emission and the excitation wavelengths. The resulting excitation and emission data sets both had 9600 elements.

The aforementioned optical calibration target provided a means to map the physical dimensions of the animal to the digital structure obtained from surface extraction. Specifically, the unit length in the CCD images was obtained by dividing the distance of the grids by the number of corresponding pixels. The coordinate system of the optical data was subsequently defined. Using the calibration procedure on the x-ray system described in [32], the unit length in the x-ray micro-CT data was determined.

The x-ray micro-CT images were manually segmented into seven different tissue types: muscle, bones, lungs, heart, liver, stomach, and spleen using Amira software (Visage Imaging, San Diego, CA). Due to the way the animal was mounted on the rotation stage, large skin flaps were created, particularly between the forelimbs and the hindlimbs. These skin flaps were also identified in the segmented micro-CT data, and later excluded from optical sampling to avoid sampling and modeling problems. The segmented micro-CT image, which contains structural a priori information of the animal, was subsequently registered to the optically extracted surface using a lab-developed MATLAB (The MathWorks, Natick, MA) program by firstly scaling the data matrices to the same spatial dimensions, and secondly aligning external structural features of the animal visible in both imaging modalities. The segmented and registered micro-CT data served as the structural basis of optical reconstruction. The optical properties, i.e., the absorption coefficient μ, the reduced scattering coefficient μ′, and the index of refraction n of different tissue types were adopted from the literature [33] and used in the forward problem solver (see Table 1 ). Note that the average values of the optical properties at the excitation and emission wavelengths were used.

Table 1

Optical properties used in the Monte Carlo simulations: absorption coefficient (μ), reduced scattering coefficient (μ′), and index of refraction (n)

	Muscle	Bone	Lung	Heart	Liver	Spleen	Kidney	Stomach
μa (mm−1)	0.0349	0.0242	0.0672	0.0275	0.1623	0.1623	0.0311	0.0069
μ′s (mm−1)	0.3709	2.2929	2.1040	0.8875	0.6371	0.6371	2.0661	1.3560
n	1.37	1.37	1.37	1.37	1.37	1.37	1.37	1.37

Spatial registration of the x-ray data to the optical data enabled transformation of the x-ray coordinate system to the optical system, in which FDOT reconstruction was realized. Because the center of the field of view in the optical coordinate system lies on the rotation axis of the animal, the origin of the optical coordinate system can be related to the physical dimensions of the optical imaging platform. The origin of the source rays was defined as the spot at the center of the second mirror of the galvo-scanner; and the origin of the detector rays as the center of the CCD chip of the camera. These two coordinates were known by physically measuring the dimensions of the optical system (tolerance <1 mm). Knowing the origin of the source ray, the point where the source ray intersects the calibration target with known location, and the 3-D animal surface, the coordinates of the source that impinges the animal were computed (Fig. 1 ). Similarly, as soon as the detector positions were determined on the CCD images, the coordinates of the detector positions were also computed.

Fig. 1

Source localization using the calibration target. The rotation axis of the animal, i.e., the blue dashed-line in the 3-D view (left) and the blue dot in the 2-D view (right), is co-planar with the calibration target. The intersection of the laser ray with the calibration target (point a in the 2-D view) is known from the calibration image, which in turn determines the intersection point of the laser with the surface of the animal (point b). Note that the position of the calibration target and the origin of the laser ray are known.

In free-space animal experiments, the air-medium interface is the animal skin in its natural state. As mentioned previously, due to the posture of the animal, large skin flaps typically were formed between the forelimbs and hindlimbs. These thin skin flaps, as well as other irregularities of the skin, caused significant problems for modeling the forward problem correctly, even using our Monte Carlo method. We implemented an anatomically guided sampling strategy to selectively exclude problematic sampling points from image reconstruction, illustrated in Fig. 2 . These problematic sampling points were identified according to the following criteria: (a) a sampling point (either a source or a detector) was intersecting a skin flap; (b) any pixel within 0.5 mm (i.e., 11 pixels in our setup) from a sampling point in the emission or excitation CCD images was saturated; or, (c) the signal-to-noise ratio (SNR) at a sampling point was less than 2. The skin flaps were identified from the micro-CT images manually. The saturation of a pixel was characterized by a value of 4095 in the CCD image resulting from an all-one representation in the 12-bit analog-to-digital converter. The SNR at a sampling point was defined as the ratio of the mean to the standard deviation of the values within the sampling area.

Fig. 2

Anatomically guided spatial sampling strategy based on structural a priori information from x-ray micro-CT. A representative pattern of spatial sampling (a grid of green dots) is superimposed on the surface (left) and the skeletons (middle) of the animal (based on micro-CT). On the right, the constituents of the segmented animal image are: (1) muscle, (2) bones, (3) lungs, (4) heart, (5) liver, and (6) skin flaps, where were used as the a priori information in the forward model; the spatial relations of the excitation laser (red arrow) and the emission fluorescence (green) with respect to the animal is illustrated, in which any emission or excitation position may coincide with the skin flaps and results in difficulty in reconstruction.

2.4. Image reconstruction

2.4.1. Forward problem

Image reconstruction in FDOT involves solving the forward and the inverse problems, both of which play important roles in overall quality of the reconstruction. As the forward problem solver, we used our Monte Carlo method previously reported in [15]. The Monte Carlo methods are a class of computational algorithms that obtain definitive results from repeated random sampling of physical or mathematical systems. The reliability of the Monte Carlo results is strongly dependent on the number of repetitions. From our previous experience of solving the forward problem in FDOT, the minimal acceptable number of photons simulated was on the order of 107, which took ~15 minutes to complete for each illumination source point using a single processor on a desktop personal computer (Dell Optiplex 745, Intel Core 2 Duo 2.4 GHz processor, and 8 GB memory). In this study, we used a total of 240 illumination points, each having 108 simulated photons, which would have required 600 hours to finish on our desktop personal computer. Another limitation of our Monte Carlo method is the computer storage required for intermediate data that pass from the simulations to the construction of the Jacobian matrix. Our previously developed software package consisted of a Monte Carlo simulation module (implemented in C programming language) and a module for constructing the Jacobian matrix (implemented in MATLAB). Using 108 simulated photons in a typical geometry for the mouse, ~1 gigabytes (GB) of data is produced for each illumination source.

Despite the above limitations, the Monte Carlo methods are the most accurate compared to alternatives, such as analytical approximation or numerical finite-element methods. To overcome the limitations of our Monte Carlo software, we unified the two separate computing modules into a single C language implementation which pass data between the two computational phases in main memory and thereby eliminate the need for intermediate disk files. The integrated software package was subsequently re-written to use the message passing interface (MPI) programming library in order to run on a large-scale parallel supercomputer. The newly implemented parallel Monte Carlo software was validated by confirming replication of the same results as the earlier implementation.

Parallel computing was performed on an SGI Altix 4700 shared-memory system, which consists of 36 computer “blades.” Each blade houses a pair of Itanium-2 Montvale 9130M dual-core processors and the four cores on a blade share 8 GB of local memory. A maximum of 132 cores can be allocated for a single run. The processors are connected by a high-speed “NUMAlink” interconnect, by which the local memory on each processor is accessible to all of the other processors. Each processor runs an enhanced version of the SuSE Linux operating system. The reduction of computing time scaled linearly with the number of allocated cores.

2.4.2. Inverse problem

The mathematical formulation modeling FDOT is very similar to that for perturbation DOT [4,9]. Therefore, the same method can be applied to compute the Jacobian matrices for both imaging modalities. In this study, we used our Monte Carlo algorithm previously developed for DOT to compute the Jacobian [34]:

where x is the quantum yield of the fluorophore; U and U are the fluence rate measured at the fluorescence and excitation wavelengths, respectively; K is the complex wave vector; L is the total path length of photon p in medium m; [F] is the value of function F evaluated at location r; N is the total number of simulated photons; and γ is a constant scaling factor that accounts for the hardware-specific characteristics, such as the insertion loss of the filters, gain of the camera, and the optical power of the light source.

The fluorescent image was reconstructed by solving the inverse problem using our previously developed regularized algorithm [21]. Image reconstruction was achieved by iteratively solving an unconstrained minimization problem

where X is the fluorescent image, and H is a cost function regularizing the reconstruction. The cost function H(X) is defined as

where ||v||2 is the L2 norm of vector v, Ω is the region of interest, and r is location. The three terms of the cost function penalize the residual error (i.e., data consistency), intensity, and smoothness of the solution, respectively, whose weights are independently adjusted by the user-specified parameters α and β. This minimization problem was solved using the conjugate gradient method with the reconstruction space constrained within the torso of the animal.

3. Results

The matrix size for Monte Carlo simulation (the forward problem) was 186 × 135 × 260 with an isotropic voxel size of 0.18 mm. The resulting Jacobian matrix for the inverse problem was down-sampled by a factor of 4 in each dimension, producing a matrix size of 46 × 33 × 65 and an isotropic voxel size of 0.72 mm. As a result, the Jacobian matrix occupied 7.6 GB of computer memory for 64-bit data precision. Using 128 cores of the supercomputer, the forward problem of the entire data set can be solved in ~200 minutes. Overall, the computation performance scaled linearly with the number of computing cores. The inverse problem solver was implemented in MATLAB and took ~15 minutes to finish on our desktop computer workstation (Dell Precision 690, dual Intel Xeon E5320 1.86 GHz 64-bit processors, and 40 GB memory). The inverse problem solver typically allocated up to 20 GB of memory. The regularization parameters, α and β in Eq. (3) were both empirically determined to be 0.1.

The primary goal of this study was to demonstrate the capability of the proposed method to resolve two closely positioned small QD inclusions in the mouse thorax. The QD inclusions were visible in the micro-CT image, which measured 0.88 mm in diameter, and in length, 1.6 mm (for the rostral inclusion) and 1.8 mm (for the caudal inclusion). Figure 3 depicts the axial slices of the FDOT reconstruction (normalized to unity) compared to the co-registered micro-CT slices. The slice thicknesses of the slices are 0.72 mm in FDOT and 0.088 mm in micro-CT. The inter-slice distance is 0.72 mm in both imaging modalities. Although an imaging artifact exists in reconstruction, the two closely positioned microliter-sized QD inclusions were successfully resolved in the FDOT reconstruction.

Fig. 3

FDOT reconstruction (the first and the third rows from the top) compared to the corresponding x-ray micro-CT slices (the second and the forth rows). The FDOT data was normalized to unity, in which the small negative values were due to reconstruction artifacts. The actual inclusions were visible in the micro-CT image (rendered in red). The inter-slice distance is 0.72 mm.

The localization error, defined as the distance from the centroid of the reconstructed inclusion to its actual location, measured 0.6 and 2.2 mm, respectively. We theorize the localization error was mainly due to the limitation of the forward model, as a result of the discrepancy between the assumed and actual optical properties of tissue. In principle, it is possible to determine the actual optical properties of tissues using DOT, which in turn can be used as functional a priori information for FDOT reconstruction, e.g., using the techniques described in [35-37]. However, to the best of our knowledge, a high-resolution DOT reconstruction technique has not been previously reported for animal experiments. Instead, we assigned literature values to the high-resolution anatomical structures obtained from micro-CT. In addition, the organs of the animal were assumed to be homogeneous and clearly separated in the forward model, whereas in reality, not only are the organs highly heterogeneous, particularly in the lungs, but also the tissue boundaries were often ambiguous in image segmentation.

We define the resolution in a similar way to the Rayleigh criterion: the maximum of one peak coincides with the half-maximum of the other peak. If these two peaks are Gaussian-shaped with equal amplitude, this definition translates to a dip of 15% between the two maxima. Based on this definition, we used a threshold of 85% to evaluate the shape metrics by measuring the diameter and the length of the reconstructed inclusions, shown in Table 2 . Note that the dimensions of the inclusions are comparable to that of the voxel size, and that each inclusion in perfect reconstruction should occupy only 2 voxels.

Table 2

Shape metrics of reconstructed QD inclusions

Inclusion	FDOT		Micro-CT
	Diameter (mm)	Height (mm)	Diameter (mm)	Height (mm)
#1 (rostral)	0.72 (1 pixel)	1.44 (2 pixels)	0.88	1.6
#2 (caudal)	1.44 (2 pixels)	1.44 (2 pixels)	0.88	1.8

4. Discussion and Conclusions

Although parallel Monte Carlo implementations have been previously reported by several groups, e.g., in [38,39], our parallel implementation is based on a different algorithm to construct the Jacobian matrix. As a result, our implementation is computationally more efficient when the number of detectors is large because it does not compute the photon density field of the detector. In addition by eliminating the intermediate data, applying our Monte Carlo method to preclinical FDOT with high-density optical sampling became practical. In conventional Monte Carlo implementations where photon density fields of both source and detector are computed, it is implicitly assumed that the source and the detector are interchangeable. This assumption is questionable for free-space FDOT because the source is typically a laser, which is highly directional and extends several millimeters into the medium (evident in Monte Carlo simulations, data not shown); and the detector is a surface area on the medium, which is much less directional (with a large acceptance angle). This assumption is particularly problematic for small animal imaging (e.g., in the mouse) where the source/detector distance is only an order of magnitude larger than the transport mean free pathlength (10-20 mm vs. 1-2 mm).

Monte Carlo methods are inherently parallel because of the highly independent nature of simulation processes. In our implementation, the entire simulation consists of multiple acquisition angles, which was specified by the user. From the perspective of simulation software, each angle was a self-contained simulation job, which consisted of a number of sources and detectors. The simulation software divided the job according to the number of sources. For each source, the computer evenly distributed the workload (number of photons) onto the allocated computer nodes. Because the number of photons was orders of magnitude larger than the number of computer nodes (e.g., 100 million vs. ~100, as in our case), the computing workload was almost perfectly balanced, i.e., all the nodes finished their tasks virtually simultaneously. Upon completion of the current source, the results were saved to the corresponding output file; and thereafter simulation for the next source started.

The optical properties of biological tissue reported in the literature are highly inconsistent, due to the differences in experimental conditions and measurement methods. They are particularly scarce for near infrared. One should also recognize that the differences among species, individuals, and physiological conditions can be significant. Although the optical properties that we referred to in this study were for bioluminescence, we found that their model and parameters produced reasonable values in near infrared compared to those available in the literature. We decided to use the referred model and parameters from the same source for consistency.

In this study, we used the same averaged values of optical properties for both the excitation and emission wavelengths using our Monte Carlo software. Note that in conventional Monte Carlo methods, the Jacobian matrix was constructed by taking the product of the fields of photon density waves irradiated from the sources and the detectors at the excitation and emission wavelengths, respectively. In our preliminary comparison, we found that our method produced better image quality (data not shown). We hypothesized that it was because our Monte Carlo algorithm more closely represented the physical process of photon migration from the sources to the detectors within the medium. Further study is underway to better investigate our findings.

Resolution is the ability of an imaging system to resolve two distinct objects, which is highly subjective to the viewer. In this study, resolution was defined as the center of a Gaussian peak coincided with the half-maximum width of another identical Gaussian peak, which translated to a threshold hold values of 85%. Note that different criteria and peak profiles can result in more or less different definition of resolution, e.g., 80% threshold if using the Rayleigh criterion and 95% using Dawes' limit for Airy diffraction patterns.

One limitation of this FDOT study was the reconstruction artifact shown in Fig. 3. This artifact was likely due to the inaccurately known in vivo optical properties, especially in the lung area where the background signal was relatively high. Nonetheless, in preclinical applications, the disease model is typically known, which assists the interpretation of imaging results and avoids misjudgment. Note that under comparatively more extreme experimental conditions, such as the case in this study, imaging artifact tend to manifest more severely. It is certainly a direction for us to continuously work on to improve reconstruction quality in future.

Another limitation was the lack of quantitative results. Although in this study, it was feasible to calibrate the scaling factor γ in Eq. (1) using another QD solution with known concentration, we noticed that this calibration technique could not be readily generalized to in vivo animal experiments, where this calibration technique would become impractical because creating a fluorescent inclusion in the thorax of the animal is difficult, and is also problematic for preclinical studies. Another concern for this calibration technique was the use of the regularization method in reconstruction due to the ill-posed inverse problem. In addition, it is known that regularization in reconstruction affects the quantitative result in diffuse optical imaging [40]. The values of the regularization parameters (0.1 for both) indicated that the reconstruction was non-negligibly influenced by the second and third terms of the cost function defined in Eq. (3). With improved experimental control, signal acquisition, and imaging apparatus, we anticipate a decreased noise level in the optical signal and, as a result, significantly less regularization will be required in the reconstruction. In addition, more accurate optical properties will likely reduce the systematic error in the forward model. Under these circumstances, it may be feasible to use this calibration technique to produce reliable quantitative FDOT reconstruction for in vivo studies.

In conclusion, we demonstrated a new FDOT method consisting of our newly developed parallel Monte Carlo software and anatomically guided sampling strategy for preclinical imaging. The Monte Carlo software fully utilized the power of parallel supercomputing and thus is better suited for handing large data sets in preclinical imaging. The anatomically guided sampling strategy significantly improved fidelity of data acquisition and subsequently the quality of reconstruction. Using the proposed method, the new FDOT system achieved significantly better resolution and sensitivity compared to our previous system, which consisted of essentially the same hardware. This method can be applied to any FDOT system adopting similar imaging principles to improve resolution and sensitivity.