System configuration optimization for mesoscopic fluorescence molecular tomography

: Tissue engineering applications demand 3D, non-invasive, and longitudinal assess-ment of bioprinted constructs. Current emphasis is on developing tissue constructs mimicking in vivo conditions; however, these are increasingly challenging to image as they are typically a few millimeters thick and turbid, limiting the usefulness of classical ﬂuorescence microscopic techniques. For such applications, we developed a Mesoscopic Fluorescence Molecular Tomography methodology that collects high information content data to enable high-resolution tomographic reconstruction of ﬂuorescence biomarkers at millimeters depths. This imaging approach is based on an inverse problem; hence, its imaging performances are dependent on critical technical considerations including optode sampling, forward model design and inverse solver parameters. Herein, we investigate the impact of the optical system conﬁguration parameters, including detector layout, number of detectors, combination of detector and source numbers, and scanning mode with uncoupled or coupled source and detector array, on the 3D imaging performances. Our results establish that an MFMT system with a 2D detection chain implemented in a de-scanned mode provides the optimal imaging reconstruction performances.


Introduction
Tissue engineering applications have known an explosive growth over the last decades. Especially, bioprinting technologies have matured such that they can create on-demand living tissue constructs for a large range of biomedical applications spanning from regenerative medicine to drug discovery [1]. To replicate the dynamically-changing biological processes under investigation, such as oncogenesis, it is necessary to develop in vitro bio-constructs that can be perfused via a vascular network to maintain long term tissue viability [1][2][3][4]. Beyond a functional network that can deliver nutrients and oxygen, the perivascular niche drives many biological processes from stem cell differentiation fate [5][6][7][8] to tumorigenesis [9][10][11][12][13]. Indeed, the three-dimensional structure of vascular networks directly affects the formation of many biological functions, such as three-dimensional gene expression, three-dimensional extracellular molecular synthesis, cell three-dimensional scaffold structure, and cell three-dimensional biological activity. Hence, the non-invasive imaging of fluorescence biomarkers that can be associated with the vascular network as well as numerous cellular functions such as cell viability, proliferation and migration for instance, is extremely desirable. However, the resulting bio-constructs are relatively large and thick (a few millimeters) such that they cannot be imaged with conventional fluorescence microscopic techniques in their entirety conversely to simple systems such as tumor spheroids. Recently, Mesoscopic Fluorescence Molecular Tomography (MFMT) has been proposed as a new imaging modality to bridge the gap between the macroscopic and microscopic regimes [14,15]. MFMT provides the unique ability to retrieve either functional or molecular contrasts at millimeter depths with high sensitivity [16] and resolution up to ∼100 µm [17][18][19]. Though, MFMT is based on an optical inverse problem that is intrinsically ill-conditioned [20]. Hence, its imaging performances are highly dependent on the selection of an appropriate forward model, dense spatial measurements and a dedicated inverse solver. To this end, we have developed efficient Monte Carlo tools to model accurately the sub-diffuse and diffuse regime [21,22], compressive sensing based inverse solvers [23], pre-conditioning techniques to reduce the ill-posedness of the sensitivity matrix [24,25], and data reduction techniques to facilitate the computation of the inverse problem [26]. However, it is well known that the intrinsic properties of the Jacobian such as sparsity, condition number, and redundancy have a significant effect on the accuracy of the reconstruction [27]. In this study, we focus on exploring the optimization of the optical system configuration associated with illumination raster scanning sampling and detection array configuration. Especially, we investigate the difference performances associated with coupled or uncoupled detection strategies. First, we set up an in silico model with trunk and branches in different directions to mimic a complex vascular network in bioprinted tissue. Then, the effects of various optical system configurations on performance of reconstruction were evaluated by the proposed metrics. The system configurations to be optimized and reconstruction strategies are described in Section 2 in detail. Section 3 summarizes the in silico experiments with reconstruction results. The measurements resulting from the combination of a high number of scanning positions and a small number of 2D detector array under coupled scanning mode of source and detector array are conducive to high-quality reconstruction. Last, the discussion and conclusion are provided in Section 4.

Principles of MFMT
Conversely to classical fluorescence microscopic techniques, MFMT relies on scattering photons to form a 3D image. To collect such photons, the detectors are positioned away from the illumination spots with increased penetration depths as the distance increases (c.f. Fig. 1(a)). In MFMT, parallel detection of multiple optode distances is performed via detector arrays, either being one dimensional or two dimensional (1D or 2D). Moreover, the source is also raster-scanned over the samples to enable a full view of it. The resulting dataset is a large combination of source detector pair positions over the sample surface that can be used in a classical optical inverse problem to retrieve the 3D bio-distribution of the fluorophore employed. As in Diffuse Optical Tomography (DOT) or Fluorescence Molecular Tomography (FMT), the detector configuration and density are central to MFMT reconstruction fidelity. Historically, MFMT has been performed with 1D arrays (APD or PMT arrays) but recent developments in Electron Multiplying CCD (EMCCD) offer acquisition of dense spatial measurements with high sensitivity at fast speeds. Herein, we investigated the performance of both 1D and 2D arrays based on our previous work.

Detector configuration
The two detector configurations we investigated herein are summarized in Fig. 1. Our original MFMT design was based around a 4×8 avalanche photodiode detector (APD) array (S8550, Hamamatsu Photonics K.K., Japan). However, to ensure fast acquisition, only a portion of the array was used leading to a 1×8 detection module being active. The imaging set up was implemented with a raster scanning of the illumination spot over the sample and de-scanned detection channel. Note that de-scanned detection here means source and detectors coupled scanned mode illustrated by Fig. 4x(a2-d2). A complete description of this platform is provided in [19]. In our previous typical application, the excitation beams scan on the imaging area with 80×200 discrete spots as source positions and emission photons are gathered by 7 detectors for each laser beam position-leading to a total of 112,000 measurements (the first detector not being used due to lack of depth information). Based on this system, we designed a second-generation MFMT system that followed the same optical principles but used an ultra-sensitive, low light Electron Multiplying Charged Coupled Devices (EMCCD, iXon EM+ DU-897 back-illuminated, Andor Technology) for detection, still in a de-scanned configuration. The EMCCD camera offers exceptional Quantum Efficiency (90-95% in visible regime), reduces the system noise due to its industrial grade packaging, and achieves a high level of signal without increasing exposure time via Electron Multiplying (EM) Gain. A Polarizing beam splitter and polarizers minimize specular reflection from the sample surface along the detection path with a fluorescence filter. A scan lens and a tube lens form a conjugate image plane and 4F magnification relay system forms the final image on the EMCCD. The picture and design of the systems are provided in Fig. 2 whereas the technical details of this imaging platform can be found in [26]. In our previous typical application, the excitation beams scan on the imaging area with 31×31 discrete spots as source positions and detectors are obtained from sparsing the EMCCD data.
The number of detectors used in the inverse problem are selected from the EMCCD via binning. The flexible configuration of the EMCCD camera allows us to set the sub-image area to any size and location on the CCD chip. When the camera is running in sub-image mode, only data from the selected pixels will be read out and those data from the remaining pixels will be discarded. In our application, we use the pixel charge aggregation offered by the EMCCD, also known as binning, to create different super-pixels forming different size of detector matrix collecting photons. By the mode, we can define and investigate different 'sub-array' size from within the full image sensor area. Figure 3 illustrates a 7 by 7 detector matrix with a size of 4.064 mm by 4.064 mm binned from the EMCCD camera pixels. In our experiments, we typically selected 48 detectors for each laser beam position, leading to a total of 46,128 measurements through scanning 31 by 31 source positions on the surface of the specimen. We exclude the middle detector for the same reasons mentioned prior that we exclude the first detector in the other detector configuration. Note that if more measurements can lead to a situation in which the number of measurements is greater than the unknown, due to the redundancy of information in diffuse optics, the inverse problem still remains ill-conditioned.

Optode scanning mode
To perform MFMT, current approaches still favor raster scanning of an illumination spot over the sample over structured light illumination [28] as the former one is still expected to provide better overall resolution [29]. The scanning of the illumination spot is performed via use of galvanometer mirror pairs or micro-stage controlled scanner. However, one main aspect of  the optical configuration in MFMT is the relationship between the illumination scans and the detection modules. Two implementations can be followed, in the simplest implementation, the detector acquisition is uncoupled to the scanning mechanism leading to a constant detection FOV over the sample and leading to the illumination spot moving over this FOV of 6.2 mm by 6.2 mm. This mode is referred to as "uncoupled scanning mode" henceforth (as illustrated by Fig. 4 a1-d1). The second approach consists of having the optical path of the detection module coupled with the scanning module in such a way that the detector FOV is not fixed but rather scanned over the sample. In this configuration, the source-detector relationship is constant for all scanning points acquired. This mode is referred to as "coupled scanning mode" henceforth (as illustrated by Fig. 4 a2-d2). This mode is a bit more challenging to implement optically and can limit the FOV on the sample due to the physical limits imposed by the required scanning lens. However, it greatly simplifies the computation of the forward model as the number of forward simulations required to form the sensitivity matrix can merely be limited to the number of detectors plus one source MC simulation.

Forward model
The mesoscopic regime is characterized by photons undergoing multiple scattering events. Hence, the diffusion equation which is typically employed for deep imaging is not valid in such shallow tissue. Herein we utilize the GPU-based software Monte Carlo eXtreme (MCX) to generate our sensitivity matrix [22] for all simulations (including the creation of synthetic data and the calculation of the forward model for experimental validation). Overall, the photon propagation is modeled over a 6.2 mm by 6.2 mm field of view (FOV) with 10 8 photon packets simulated for each optode. For simplicity, we assumed a planar boundary condition.
In all cases, herein, the adjoint method was used to compute the different Jacobians due to its computational efficiency [29]. In the case of a de-scanned configuration, we replicated the strategy as laid out in [30]. Briefly, in such cases, thanks to the symmetry of the problem, the forward model is computed only for one source and the number of detectors considered. Then the overall Jacobian is populated using this submatrix based on the position of the source in the image space. Hence, this is a very fast and efficient computational scheme. Conversely, in the case of a scanned configuration, the forward model needs to be computed for every source and detector combination, leading to an increased computational burden. In all cases, the entire scanning area with FOV of 6.2 mm × 6.2 mm was simulated. Additionally, a padding area is incorporated to account for boundary conditions at the border of the scanning area. In the cases, the separation between adjacent detectors binned from 2 by 2 pixels is set to 0.64 mm, which leads to a detector array with a size of 4.064 mm by 4.064 mm (as shown in Fig. 4). Note that some scanning spots and scanning paths are omitted for the sake of simplicity in Fig. 4. For demonstrative purposes, we only showed the extreme optodes that are required to produce the sensitivity matrix with all other source-detector pairs laid between those points. In Fig. 4, black squares represent the scanning area on which the raster scanning was performed and green squares show the active area of EMCCD representing the detector array.

Inverse solver
An l 1 -regularized least square based inverse solution was employed to retrieve the distribution of the fluorescence inclusion in a homogenous media which can be cast as a linear model of the form with a noise perturbations v, where X ∈ R n is the vector of unknowns, Y ∈ R m is the vector of independent measurements, v ∈ R m is the noise, and A ∈ R m×n is the Jacobian matrix. Like many fluorescence optical tomography applications, the number of independent observations m is always not large enough compared to number of unknowns n in our MFMT to lead to an under-conditioned inverse problem. L 1 -regularized least squares based methods have been proved as an effective means to retrieve the solution of under-conditioned problem [23,27]. More importantly, l 1 -regularized least squares typically yields a sparse vector X. Then, the under-conditioned problem can be transformed the following Lagrange multiplier expression, the l 1 -norm of X and λ>0 is the regularization parameter. The l 1 -regularized least square problem (2) is convex but not differentiable. Here we follow the standard convex optimization methods to reformulate Eq. (2) as a convex quadratic problem with linear inequality constraints. The equivalent quadratic program can be solved by interior-point methods [25,27]. For each of the following simulations and experiment cases, we employed the same inverse solution methodology proposed in [25]. The optimal regularization parameter λ can be chosen through L-curve analysis for different Jacobian matrix [31].

Image-based quantitative metrics
To ascertain the merit of each configuration explored herein, it is important to define meaningful quantitative metrics. Such metrics are typically derived either from the characteristics of the forward model (Jacobian) or directly from the reconstructed images. For instance, Singular Value Analysis (SVA) of the Jacobian has been extensively used in the field to investigate the different parameters associated with the measurement space, including the optode configuration [32][33][34] and/or the data set type. Similarly, the Cramér-Rao lower bound method that is often used in statistical signal processing has been employed to investigate the optimal system configuration for diffuse optical tomography [35,36]. Though, in each case, if these metrics provide useful information about the relative merit of the instrument parameter under investigation, they do not always correlate with the observed imaging performances [37,38]. Hence, herein, we decided to employ four volumetric metrics that are typically employed in the image processing field to evaluate similarity between two set of images, in our case, the reconstructed image and the ground truth. To perform comparison fairly for different scenarios, these metrics are normalized so that they produce values between 0 (for two absolute different results) and 1 (for two totally coincident results). The adopted metrics here provide a comprehensive evaluation in term of similarity (normalized sum of square of the differences, nSSD), displacement of mass center (normalized sum of the absolute differences, nSAD), determination of spatial alignment (normalized correlation, nR), and disparity (normalized disparity, nD) [25]: where A(i, j, k), B(i, j, k) are the value of the normalized numerical model and reconstructed result at coordinates (i, j, k), N = R × C × Z is the total voxels in a cubic phantom and R, C, Z are the number of voxels three axes, respectively. A , B are binaries (0 or 1) of A, B through thresholding, and ⊕ represents the exclusive or (XOR) operation. Thanks to the normalized expressions, the values of these four metrics are adjusted to the range between 0 and 1 with a consistent trend. A larger value of the metrics indicates the reconstruction is totally coincident with the ground truth.

Simulation models
In silico models were designed to optimize the configuration of optical system and evaluate the performances of the proposed method in terms of reconstruction accuracy and spatial resolution.
First, an elaborate model representing a vascular tree with main trunk and three groups of off-shoots with 200 µm separation along three coordinate axes was set up. The different diameters of the vascular tree components are 200 µm and 400 µm, respectively, as shown in Fig. 5. The fluorophore concentrations are assumed to be homogeneous over the vessel with effective quantum yield equal to one. The simulated imaging domain is a 6.2 mm × 6.2 mm surface area with a 3 mm depth. In this paper, we discretized the imaging domain into 31×31×15 voxels with 200 µm discretization, leading to a Jacobian of size 46,128 (961 scanning spots, 48 detectors) by 14,415. Note that we neglect the center detector as it provides no depth information but also lead to dynamic range limitations experimentally. We also assume a homogeneous scattering coefficient (µ s =5 mm −1 ) and endogenous absorption coefficient (µ a = 0.002 mm −1 ) over the domain. The Henyey-Greenstein phase function was used for all computation herein with g = 0.9. These parameters derive from the optical properties of the collagen scaffold employed in our bio-printing application at collagen density of ∼9 mg/ml [19]. The vascular tree was first employed to compare the performance of reconstruction between one dimensional (1D) and two dimensional (2D) detectors placement to investigate the merit of using 2D detectors. The same number of detectors is radially spaced (1D) and is placed in a square grid formation (2D), having the same maximum source-detector separation, as shown in Fig. 1., we further observed the effect of detector number on the reconstruction. In this experiment, again we kept the maximum source-detector separation fixed (i.e. 4.064 mm) while changing the density of the detector array: 3 by 3, 4 by 4, 5 by 5 etc. Secondly, the number of source positions (scanning spots) was examined with regards to reconstruction performance while keeping the number of detectors fixed at 49 detectors (7 by 7). Although the 1 st and 2 nd experiments provided us some understanding of how the data size affect the reconstruction, it is not clear that how these two components of the measurement vector would affect one another. To address this, we design another experiment with approximately the same amount of measurement data by changing the number of detectors and number of sources respectively. Lastly, two scanning mode experiments were designed in order to determine which model would lead to optimal reconstructions. For all the above experiments, we added 10% Gaussian noise to the synthetic data -mimicking the real noisy experimental data. In all cases, both qualitative results and four quantitative indexes are given to evaluate the reconstruction performance in different scenarios. Meanwhile, the coupled scanning mode is adopted for all following simulations and experiments except for the simulations for coupled and uncoupled scanning mode in Section 3.4.

Effect from different layouts of detectors on reconstruction performance
We evaluated the effect of different detector layouts (1D vs. 2D, as shown in Fig. 1) when the total detector number is fixed. We tested typical layouts of detector array such as 3×3 vs. 1×9, 4×4 vs. 1×16, and so on. A visual inspection and quantification metrics both showed that more detectors can help to achieve better reconstruction results. Furthermore, 2D detector arrays always delivered better reconstruction performances than 1D detector arrays with the same detectors. A good reconstruction is first obtained by a 6×6 detector array, followed closely by a 7×7 array, which underscores the impact of increasing the number of detectors, as shown in Fig. 6 and Fig. 7. We further investigated the effect by increasing the detector number beyond 7×7 but the gains were minimal whereas the computational burden increased significantly (see Fig. 6). Our comparison study showed that in both cases reconstruction fidelity increased as we increased the number of detectors. The evaluation metrics showed that 2D detector arrays delivered better reconstruction performance than 1D detectors for all our performance metrics Fig. 7. The first row shows the 3D reconstructions of 1D detector arrays from 1 by 9 to 1 by 81 (a1-d1), while the second row shows the reconstructions of 2D detector arrays from 3 by 3 to 9 by 9 (a2-d2).
(see Fig. 6, 7). Although the difference is small, it is technically difficult to employ a line of more than 49 detectors for an optical setup. Thus, 2D detector arrays offer simplicity in optical setup along with better performance.
As it is shown in Fig. 6, evaluation metrics indicate an asymptotic performance with the increase in detector number. However, that brings a computational burden as the measurement data size increases as well as a technical difficulty to place such a high number of detectors in a line. To highlight the optimum performance and complete characteristic of higher detector number, we carried the study one step ahead (to the maximum capacity of our personal computer) to investigate the effect of size of detector matrix such as 3×3, 4×4, 5×5, . . . . . . , and 9×9, on reconstruction performance. This simulation in Fig. 7 showed us increasing the detector number beyond 7×7 doesn't bring a significant benefit. Overall most of the metrics reached a plateau at 8×8 and 9×9 (see Fig. 6). Another restriction for the optical system is having a symmetric layout of the detectors. We reserved the middle detector for the source. The rest of the detectors are required to be symmetric about the source point, which adds the constraint of having an odd number of detectors. Consequently, we employed 7×7 detector array for our optical system.

Optimization of the number of scanning spots
The location and number of scanning spots are other important parameters as they determine the number of the measurements, the properties of sensitivity matrix, and the experimental acquisition time. As mentioned above, due to the advantage of 2D detector layout, we chose to use a 7×7 detector array. Below we demonstrate that the effect of scanning spots is stronger than the effect of the detector number. Similar to Fig. 6, metrics from scanning spots showed an asymptotic behavior, reaching a plateau at 26×26 positions, as shown in Fig. 8(a). Figure 8(b-d) shows that increasing the number of sampling points contribute to higher reconstruction qualities. To counterbalance the loss of a relatively small detector number, we decided to adopt more scanning spots (31×31). The cooperative effects of both number of detectors and number of source positions is further investigated in the next section.

Combined effect of detector and source number on reconstruction
The total number of measurements is the product of number of scan positions and number of detectors. Mathematically, the numbers of equations are approximately the same for all cases. However, it is not clear how the information content changes with different numbers of measurements. Here we demonstrate 6 different measurement numbers and their reconstructions to compare the data quality (see Fig. 9, 10). Results indicated that the number of scanning spots has more impact on the data quality and the reconstruction than the detector number does. This confirmed our selection of number of scanning spots and detector number above. Although the measurement numbers (number of scanning spot multiplied with the detector number) are approximately the same, the reconstruction performances are significantly different. Both factors were shown to be effective, however this simulation shows that the number of sources has a more drastic impact on the reconstruction performance.

Coupled and uncoupled scanning mode on reconstruction performance
Based on the layout of detectors, number of detectors, and scanning spots optimized above, we further investigate the effect of the coupled and uncoupled scanning mode of source and detector array on reconstruction quality. We employed two numerical models: a simple vessel model (M1) and a complex vascular structure (M2) to test the performance of both scanning modes. The visual results of the models plainly demonstrate the reconstruction achieved via coupled scanning mode (see Fig. 11(c) and (f)) outperforms its counterpart's output (see Fig. 11(b) and (e)) under the same circumstances. As also can be seen from Table 1, all the metrics verify that coupled scanning mode contributes to attain high fidelity results. Especially, the normalized disparity (nD) between original phantom and reconstruction via coupled scanning mode even is about twice as much than the uncoupled one. Additionally, the computation times to generate the full Jacobian matrix in the case of a coupled scanning mode is 8∼10 times faster than in the case of the uncoupled configuration.

Experimental phantom reconstruction with the optimal optical configuration
Last, we test the best optical configuration experimentally using our previously reported platform [25]. Hence, the scanning spots and 2D detectors array are respectively set to 31×31 and 7×7 and with coupled scanning mode. An agar phantom with optical properties of µ a = 0.002 mm −1 , µ s ' = 1 mm −1 , n = 1.34, and g = 0.81 was built with a PEGDA fluorophore letter R (GFP 488/509) placed 0.7∼0.8 mm beneath the sample surface. Figure 12(a) shows the size and shape of the PEGDA fluorophore letter R before it was placed into the phantom. The image space was discretized in 200 × 200 × 200 µm 3 voxels as in the in silico study. Several commonly used inverse solvers, such as conjugate gradient (CG) and l 1 norm based inverse solver [23], are employed to reconstruct the letter "R" from phantom to test the optimized  optical configurations. In order to get optimal results for l 1 norm based inverse solver, we employed the L-curve method to find the regularization parameter λ when residual norm tend to least value [19]. As can be seen from Fig. 12 both inverse solvers achieve a good reconstruction performance with the help of the optimized system configuration. Here we employ Optical Coherence Tomography (OCT) to get the ground truth of the phantom of letter R. Compared with the ground truth, l 1 norm based solver achieve higher metrics value than its counterpart as shown in Table 2. Note that this result is expected as the l 1 norm based was employed to identify the optimal configuration settings. Though, we have also reported on similar enhancement when using the l 1 norm over CG in earlier work [23] that did not use the l 1 norm ab initio as the best solver.

Conclusions
We optimized the system configuration for our EMCCD based MFMT system by virtue of the flexible binning mode and single photon detection capability of the EMCCD. According to the evaluation metrics of reconstruction results, we selected a high number of scanning spots (31×31) and a small number of 2D detector arrays (7×7) as the optimized combination. As for scanning mode, we employ coupled source and detector array to utilize the symmetry of FOV to fast get the sensitivity matrix as well as the preferable formation of it. We employed this optimized optical configuration to test the reconstruction performance on in vivo experiments. The results from both inverse solvers produce better reconstructions. Next, we will further explore more robust algorithms to improve the reconstruction quality.