*in silico*studies in a synthetic small animal model for typical illumination patterns. Second, the applicability of this approach in tomographic imaging is validated

*in vitro*using a small animal phantom with two fluorescent capillaries occluded by a highly absorbing inclusion. The simulation study demonstrates an improvement of signal transmitted (∼

**2**orders of magnitude) through the central portion of the small animal model for all patterns considered. A corresponding improvement in the signal at the emission wavelength by 1.6 orders of magnitude demonstrates the applicability of this technique for fluorescence molecular tomography. The successful discrimination and localization (∼1

**mm**error) of the two objects with higher resolution using the optimized patterns compared with nonoptimized illumination establishes the improvement in reconstruction performance when using this technique.

## 1.

## Introduction

Fluorescence molecular tomography (FMT) is a rapidly growing imaging modality that allows the noninvasive assessment of the molecular state of biological tissue *in vivo*.^{1} The advances in transgenic animal models^{2} and functionalized fluorescence molecular markers^{3} have established fluorescence imaging as a pivotal modality in current preclinical biomedical research. A recent development in FMT techniques has been the advent of wide-field excitation schemes where point-sources of light are replaced by spatially extended sources.^{4}5.^{–}^{6} Such tomographic approaches are herein referred to as wide-field optical tomography. Compared with point-source excitation schemes, they allow the collection of dense spatial tomographic information over large surfaces at fast acquisition speed.^{7} Moreover, wide-field illumination patterns can be scaled to image a large volume without additional acquisition time. This makes wide-field molecular tomography a particularly attractive technique for fast whole-body small animal imaging.

The application of wide-field excitation schemes to small animal whole-body FMT was first demonstrated using broad bar-shaped patterns in transmittance mode.^{6} The use of bar patterns (each illuminating half of the imaging volume) had twofold advantages. First, it provided a high signal in the recorded measurements due to the injection of a large number of photons into the medium. Second, the low spatial frequency content in such patterns ensured maximal transmission of spatial information through biological tissue.^{8} Concurrently, a wide-field patterned excitation scheme for fluorescence imaging of thick tissue using spatially sinusoidal patterns (similar to illumination patterns used in structured illumination microscopy) was experimentally demonstrated for optical sectioning studies of shallow tissues in reflectance mode.^{9} The application of similar illumination patterns in a tomographic inverse problem has also been demonstrated in phantoms for the reconstruction of absorptive heterogeneities in thin slab geometry.^{10} Likewise, wavelet-based illumination patterns have been used in FMT of phantoms with faster reconstruction performance.^{4}

Despite developments in wide-field patterned excitation schemes, the highly scattering interaction of photons in the near-infrared spectral range coupled with the reduction in sensitivity due to wide-field illumination presents a highly ill-posed inverse problem in FMT.^{7} An approach towards reducing the ill-posedness in the inverse problem is to maximize the information content in the recorded measurements. As in the case of classical optical tomography techniques, this can be achieved by using advanced measurement data types in the frequency-domain (FD) or time-domain (TD). Alternatively, wide-field FMT provides an approach wherein the recorded tomographic information can be maximized by optimizing the excitation pattern employed. In this regard, a recently developed model-based method employs a computational optimization procedure incorporating the *a priori* knowledge of animal geometry and distribution of optical properties *in vivo*.^{11} However the intraspecies variability of these parameters necessitates a subject-specific optimization protocol limiting its application in experimental settings due to the extreme computational burden.^{12} Moreover, the model-based optimization of patterns does not ensure the quality of measurements obtained in experimental settings rendering such an approach impractical for experimental implementation.

In this work, we present a new wide-field optical tomography technique, referred to as adaptive wide-field tomography (AWFT), wherein the excitation patterns are iteratively updated during acquisition to maximize the information content in the recorded measurements. Specifically, this approach aims at increasing the number of tomographic projections acquired with high signal by reducing the large dynamic range in photons transmitted through the small animal body. This is achieved by locally controlling the amount of light injected on the illumination surface of the sample. In whole-body small animal imaging, the dynamic range of transmitted photons can be attributed to two factors; namely, the animal specific geometry exhibiting variations in thickness across the imaged volume (e.g., neck to tail) and the wide range of optical properties of small animal organs.^{12} The combined effects of these two factors result in significant variations in signal attenuation across the small animal body. The effects of these factors become even more critical for free-space imaging applications in the absence of matching liquid chambers.

The effects of the large dynamic range are particularly critical in time-resolved imaging for two reasons. First, the saturation of the detectors typically used for time-resolved imaging result in nonlinearities in the measured signal, thus limiting the maximum number of photons which can be recorded. A similar limitation on the maximum number of photons detected can also be observed in CW-imaging strategies employing low-cost cameras with smaller bit-depths. Second, the measurements acquired in time-resolved imaging using photon-counting techniques are susceptible to additive noise having a Poisson distribution wherein the noise is dependent upon the photon counts.^{13} In other words, measurements with high photon counts have a higher signal-to-noise ratio (SNR), which drops significantly with decreasing photon counts. As a result, measurements with lower photon counts are typically not used during reconstruction to avoid artifacts arising due to uncertainty in photon counts. The combined effect of these two characteristics of time-resolved techniques, when imaging subjects exhibiting large dynamic range, is the reduced number of measurements with sufficiently high photon counts for model based studies.^{13}^{,}^{14} The optimization of transmitted signal therefore becomes critical in time-resolved optical tomography to ensure the availability of detector readings with sufficiently high SNR.

For instance, detectors positioned in the central part of the animal or positioned below highly attenuating organs acquire weak signals highly affected by noise,^{14} which can lead to poorer reconstruction performances. Figure 1(b) shows an example of the dynamic range observed in an animal model upon illumination by a uniform pattern shown in Fig. 1(a). Here the limited transmission of photons through the thoracic cavity due to highly absorbing organs (e.g., heart and lungs) combined with the high signal transmitted through the edge of the model (due to reduced thickness) results in a high dynamic range leading to a limited number of detectors around the edge of the model with sufficient signal levels. In such scenarios, AWFT iteratively adapts the spatial intensity distribution across the excitation pattern, allowing one to locally increase the laser power, thus increasing the number of detectors acquired with high signal for each individual pattern across the animal body. The iterative correction of the excitation pattern in AWFT is implemented by a measurement-driven pattern optimization scheme, which does not require *a priori* knowledge of model geometry and optical properties. Therefore AWFT provides an experimentally efficient *in situ* pattern optimization scheme in wide-field optical tomography. It is also worth noting that the optimization of the excitation pattern in AWFT based upon the transmitted signal at the excitation wavelength produces a corresponding improvement in the fluorescence signal measurements (using identical patterns), which improves the performance of wide-field FMT.

Herein, AWFT is applied to time-resolved FMT. The improvement in the robustness of excitation and fluorescence measurements in transmittance across the animal body when using the time-gate datatypes is first demonstrated *in silico* using an anatomically accurate synthetic mouse model. Subsequently, the optimization procedure is validated *in vitro* in a solid phantom mimicking a small animal model. The improvement in tomographic information content upon optimization is quantified, and the resulting improvements in the reconstruction performances are demonstrated.

## 2.

## Materials and Methods

## 2.1.

### System Description

Figure 2(a) shows the design of the time-resolved imaging system used in this study and is described in detail in Ref. 15. Briefly, the platform incorporates a tunable Ti-Sapphire laser (Mai Tai HP, Spectra-Physics, USA) used in conjunction with a closed-loop power control system (providing up to 1.8 W power at 800 nm). The time-resolved measurements of the temporal point spread function (TPSF) are made using a gated intensified charge coupled device (ICCD) camera (Picostar HR, LaVision GmbH, Germany). The image intensifier is operated at 300 ps gatewidth at a temporal resolution of 40 ps over a 4 ns time window with images recorded using a 12-bit CCD camera. It should be noted that the noise characteristics of the ICCD are highly dependent on the voltage applied across the image intensifier.^{15} In the studies described herein, the intensifier was operated at 570 V. Thus the gated measurements having less than 200 counts were excluded from reconstruction to avoid noise artifacts.^{16} The fluorescence measurements are acquired using discrete filter sets (Semrock, USA). A noteworthy feature of the system employed in this study is the incorporation of a Pico projector module (Optoma PK101) as the pattern generation mechanism. The in-built light source in the projector is replaced by the output from the laser. The Pico projector allows the easy implementation of grayscale patterns (via a USB connection using PowerPoint presentation software) with a compact light engine as opposed to the Discovery 1100 kit implemented in the previous incarnation of the system.

The spatial and temporal characteristics of the system when using the Pico projector were established using a grayscale gradient pattern shown in Fig. 2(b). The objective of this investigation was to determine the range of gray levels that can be accurately projected and detected by the system. Moreover, the temporal characteristics of the instrument response function (IRF) of the system when using the grayscale pattern were evaluated (as the gray levels are generated on the digital micro-mirror device using modulation frequencies in the kHz range) in order to establish the range of gray levels that can be accurately implemented in the optimization procedure without introducing temporal errors during time-gated acquisition. The temporal response of the pattern was therefore acquired for the grayscale pattern and a uniform full-field pattern using a diffusing white paper. The maximum intensities recorded for each pixel for the gradient pattern were normalized by the corresponding intensity values in the uniform pattern. This was done to mitigate error in grayscale estimates due to nonuniform intensity distribution in the projected pattern owing to the projector optics. The normalized grayscale gradient pattern recorded by the camera is shown in Fig. 2(c), and the profile across the center of the pattern shown in Fig. 2(d) shows an accurate ($<10\%$ error) generation of grayscale intensities for levels greater than 0.25. Figure 2(e) to 2(f) compare the photon arrival times (${t}_{0}$) and the IRF width (at half maximum) across the grayscale pattern and the uniform pattern. It is worth noting that the generation of gray levels has a negligible effect on the system IRF. The spatial variation in the IRF characteristics are, however, nonnegligible for both patterns and must be accounted for in the forward model while setting up the inverse problem. Herein, the IRF recorded during pre-experiment calibration is convolved with the time-resolved forward model to accurately model the light propagation in this system. Note that even though the Pico projector is a consumer-grade product with low-quality optics; it allows one to implement a wide-field optical tomography at very low cost ($<\$200$) and can replace advantageously galvo-scanning based systems even in the case of time-resolved studies.

## 2.2.

### Reconstruction Scheme

As stated previously, in this study we employ time-gates as the datatype for reconstruction. Owing to the complex geometry and wide range of optical properties encountered in small animal models, we implemented a time-resolved Monte Carlo (MC) based method to generate an accurate photon propagation forward model. The use of an MC-based method further ensures the accuracy of the model when using early time-gates as datatypes for solving the inverse problem.^{17} The MC model used in this study has been described in detail in Ref. 16 and is discussed briefly below.

Consider a volume $\mathrm{\Omega}\in {R}^{3}$, with photons incident over an area $A\in {R}^{2}$, then the time-resolved signal detected at the point detector at ${r}_{d}$ is given by Eq. (1).

## (1)

$${U}_{F}({r}_{s},{r}_{d},t)=\underset{\mathrm{\Omega}}{\int}\mathrm{d}{r}^{3}{J}^{\mathrm{fluo}}({r}_{s},{r}_{d},r,t)\eta (r)\phantom{\rule{0ex}{0ex}}\eta (r)=\frac{{\mu}_{\mathrm{af}}^{x}(r)\varphi}{{\mu}_{a}^{x}(r)}.$$Here ${r}_{s}$ is a point in $A$, $\eta (r)$ is the effective quantum yield for a fluorophore with quantum yield ($\varphi $), ${J}^{\mathrm{fluo}}({r}_{s},{r}_{d},r,t)$ is the spatial and temporal sensitivity function with respect to $\eta (r)$ for the time gate $t$ and total absorption coefficient ${\mu}_{a}^{x}$ at the excitation wavelength with a contribution ${\mu}_{\mathrm{af}}^{x}$ from the fluorophore. The absorption coefficient of the fluorophore is linearly related to the fluorophore concentration via the extinction coefficient. However, the uncertainty in the values of quantum yield and extinction coefficient of a fluorophore *in vivo*, which may change due to variations in the microenvironment biochemical properties, introduces uncertainties in the reconstruction of the absolute fluorophore concentration biodistribution. In this work, we therefore report relative effective quantum yield as the quantitative reconstructed parameter. For a fluorophore in the medium that has a mono-exponential decay with lifetime $\tau $, under the assumption of equal absorption and scattering coefficients at the excitation and emission wavelengths, the fluorescence signals can be simulated by convolving the temporal signals generated at the excitation wavelength by the exponential time decay of the fluorophore given by ${e}^{-t/\tau}$. The Jacobians for FMT, ${J}^{\mathrm{fluo}}(t)$ can therefore be calculated using Eq. (2).

## (2)

$${J}^{\mathrm{fluo}}(t)=\underset{0}{\overset{t}{\int}}\mathrm{d}{t}^{\prime}{J}^{\mathrm{exc}}({t}^{\prime})F(t-{t}^{\prime})\phantom{\rule{0ex}{0ex}}F({t}^{\prime})=\underset{0}{\overset{{t}^{\prime}}{\int}}\mathrm{d}{t}^{\prime}\mathrm{IRF}({t}^{\prime \prime}){e}^{-({t}^{\prime}-{t}^{\prime \prime})/\tau}.$$${J}^{\mathrm{exc}}(t)$ is the time-gated weight function at the excitation wavelength and $\mathrm{IRF}(t)$ is the experimentally recorded IRF. It should be noted that the experimentally recorded patterns are used to calculate ${J}^{\mathrm{exc}}(t)$ ensuring that the spatial characteristics of the system are accurately modeled by the MC model. The forward model ${J}^{\mathrm{exc}}$ is computed in a massively parallel computing environment (Blue Gene, Center for Nanotechnology Innovations, RPI) using 1024 nodes allowing the computation of the time-gated Jacobian in less than 7 min for each pattern and all detectors.

To cast the inverse problem, we employed a Born formulation where the experimental time-domain emission measurements are normalized by the continuous wave (CW) excitation flux at the same position. This normalization efficiently mitigates the dependence of the detected fluorescent signal on the optical properties of the examined tissue.^{18} Note also that this formulation employs the CW excitation flux at the same position to alleviate the unavoidable time errors associated with drift and jitter in time-resolved studies. The time-gated Jacobians are computed using the perturbative MC technique, which is the most computationally efficient when large number of detectors are considered.^{19} The inverse problem is solved using a least-squares solver, *lsqr* (MATLAB, Natick, Massachusetts). The solver was stopped after 200 iterations for all studies herein.

## 2.3.

### Optimization Algorithm

The iterative correction of excitation pattern in AWFT modifies the spatial distribution of incident power by modulating the gray levels in the pattern. The optimization procedure therefore employs the measurements from the detectors spanning the area of the pattern on the surface of the model. As a part of pre-optimization calibration, the excitation pattern is measured using a white diffusing paper, and its spatial characteristics are determined. Next a transformation matrix mapping the image used to define the full field pattern to the measured full field is computed. This is done to account for spatial deformations in the projected pattern due to the projector optics.

Consider the animal model shown in Fig. 3(a) where the white border represents the excitation pattern. The optimization of the patterns described here is based on time-resolved data, and the process begins with the acquisition of temporal measurements for each of the $N$ patterns in the user-defined set with the laser operating at power ${\mathrm{Pow}}^{i}$ ($i=1,2,3\dots $) at the $i$th iteration. As the dynamic range limits the maximum photons recorded in each time-resolved curve, for each pattern ${P}_{n}$ ($n=1\dots N$), an excitation map (${E}_{n}$) is constructed from the maximum measured counts (peak counts of the TPSF) at each detector. The underlying principle of the optimization algorithm considers each pixel in the pattern independently and corrects the grayscale level based on the excitation signal in transmittance. Therefore, at each pixel at $(x,y)$ in the $n$th pattern at the $i$th iteration the multiplicative update value $\alpha $ is computed and the corresponding pixel is updated using the formula given in Eq. (3). It should be noted that in Eq. (3), the calculation of $\alpha $ includes the weighted photon counts within a neighborhood of radius $r$, in order to incorporate the effect of diffusion on transmitted photon counts.

## (3)

$$\alpha =\frac{\chi}{\frac{\sum _{k=1}^{K}{E}_{n}({x}_{k},{y}_{k})\times {d}_{k}}{\sum _{k=1}^{K}{d}_{k}}}.$$In Eq. (3), $\chi $ is the maximum desired photon counts, $({x}_{k},{y}_{k})$ is the position of the $k$th neighbor within a radius $r$, ${d}_{k}$ is Euclidean distance of the $k$th neighbor from $(x,y)$, and $K$ is the total number of detectors included in the neighborhood. The selection of detectors within distance $r$ incorporates the effect of photon diffusion when computing the update value for the pattern. Furthermore, we define $\beta $ as the maximum power amplification allowed when acquiring the TPSF. This ensures that the algorithm does not exceed the maximum permissible exposure limit of laser power during optimization. The iteration is completed after each pixel in the pattern has been updated. The update values obtained after this step for pattern ${P}_{n}$ (referred to as ${A}_{n}$) represent the relative change in laser power required at each position on the pattern to achieve desired excitation signal level. The experimental implementation of this updated pattern is achieved by normalizing the pattern to 1 and increasing the laser power in the subsequent step by a factor of the maximum update value, ${A}_{n}^{\mathrm{max}}$.

This process is repeated until ${A}_{n}^{\mathrm{max}}$ is less than or equal to 1 (indicating convergence with all detectors reaching $\chi $ counts). Furthermore, we include a second stopping criterion to ensure the termination of the optimization scheme when ${A}_{n}^{\mathrm{max}}$ converges to a constant value greater than 1, wherein the iterations are terminated when the maximum change in the new pattern (${\parallel {P}_{n}^{i+1}-{P}_{n}^{i}\parallel}_{\infty}$) is less than 10%. The above process is repeated for all patterns in the base pattern set, and the final set of patterns is referred to as the transmittance-optimized patterns. This optimization procedure is outlined in Fig. 3(b).

## 2.4.

### Grayscale Limits During Pattern Optimization

As noted in the system description, the experimental implementation of the adaptive optimization of the patterns was done using the Pico-projector system. The nonlinearity of the projected gray levels as shown in Fig. 2(d) results in inaccurate projection of gray levels below 0.25. In order to mitigate the effect of this hardware characteristic, the optimization of patterns was limited to grayscale values above 0.25 with all pattern positions having values smaller than 0.25 remaining unchanged during optimization. Furthermore, all patterns obtained during the optimization process were acquired post-hoc and used in the Monte Carlo simulations in order to ensure the accuracy of the model during reconstruction. It should also be noted that the *in silico* validation of the algorithm applied no limits on the pattern.

## 2.5.

### Neighborhood Radius

As stated above, the radius of the neighborhood incorporates the effects of diffusion when updating the gray level for a pixel on the pattern. The ideal radius selection for each pixel would typically include the height of the model at the pixel and the local optical properties, leading to subject dependent radius and thus necessitating *a priori* knowledge of the anatomy. In this work, we select the radius based on the distance of the edge of the model from the pattern to ensure that all photons exiting the model are accounted for in the optimization process. It should be noted that a radius of 2.5 mm was selected for the *in silico* studies and 4 mm for the experimental validation. Note that simulations with radii of 3 and 8 mm were also carried out leading to identical convergence rate as reported here when using 2.5 mm (results not shown).

## 2.6.

### Batch Optimization of Patterns

The optimization procedure shown in Fig. 3(b) will result in each pattern converging to the optimality condition at a different rate based upon the region of excitation and will require different power amplification. For instance, photons from a pattern probing the thoracic cavity will be highly attenuated due to highly absorbing and scattering organs (e.g., heart and lungs), while the signal from another pattern in the same set probing the abdominal area will encounter a thicker volume with low absorbing tissue (e.g., stomach, intestines etc.). This leads to variation in the signal attenuation profiles recorded on the surface for each pattern and necessitates the individual optimization of each pattern in the set.

In this implementation of the above algorithm, all patterns were batch optimized to ensure experimental efficiency (acquisition time). In other words, the temporal measurements for all patterns in the set were acquired consecutively and followed by the optimization of each pattern in the set. Specifically, the optimization algorithm was implemented in two steps. First, the pattern update value was unbounded ($\beta =\infty $) during optimization and the optimal patterns were computed. As stated previously, each pattern will have been assigned different power amplification for the next iteration. In the second step, the minimum power amplification among all patterns in the set is put as the upper bound on the value of $\alpha $ [$\beta =\mathrm{min}({A}_{n}^{\mathrm{max}})$] where, $n=1,\dots ,N$. This ensured that the laser power employed during the acquisition step will not saturate the ICCD for any pattern in the set. Moreover, the consistent source power across all patterns allowed a rapid acquisition of the complete pattern set with high experimental efficiency (typically less than 5 min for the optimization of a 36 patterns set). It should be noted that this algorithm can also be applied to CW and FD datatypes without loss of generality.

## 2.7.

### In Silico Validation

The efficacy of the optimization algorithm for whole-body imaging of small animals was validated using a synthetic mouse model (Digimouse).^{20} Figure 4 shows the model employed, where seven major organs in the animal torso were modeled to simulate the propagation of photons in small animals. The optical properties at 730 nm computed using values given in Refs. 21 and 22 were assigned to each organ and have been compiled in Table 1.

## Table 1

Optical properties of major organs in a small animal model at 730 nm.

Organ | μa (cm−1 | μs′ (cm−1 |
---|---|---|

Heart | 0.35 | 23.0 |

Lungs | 0.25 | 30.0 |

Liver | 0.90 | 6.0 |

Spleen | 0.90 | 6.0 |

Stomach | 0.05 | 13.0 |

Kidneys | 0.20 | 20.0 |

Bladder | 0.02 | 8.0 |

Bone | 0.10 | 20.0 |

Muscle | 0.30 | 10.0 |

Furthermore, two $3\times 3\times 3\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ fluorophore inclusions were positioned as shown in Fig. 4(a) to simulate the effect of optical properties and/or geometry on the fluorescence signal intensity. An area of illumination spanning the whole body (as shown in Fig. 4) was selected as the base pattern, and the excitation field was recorded at discrete detectors at 1 mm separation in transmittance mode. Three different illumination patterns were considered for evaluation. First, a uniform bar pattern,^{6} second a wavelet pattern (Battle-Lemarie, Scale 1),^{4} and last a sinusoidal pattern (spatial frequency of $0.125\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$) as shown in Fig. 4(b) to 4(d), respectively.^{9}^{,}^{10} It is worth noting that the Monte Carlo simulations of time-resolved excitation and fluorescence signals also exhibit the Poisson noise characteristics of time-resolved imaging instrumentation. Thus no noise was explicitly added to the simulated measurements.

The optimization scheme was employed to achieve $\chi =4000$ counts at each detector, which is less than the maximum number of photon counts that can be detected by the ICCD camera used in the study (12-bits). After convergence, the fluorescence forward model was computed for the uniform and optimized pattern to determine the corresponding changes in fluorescence signal intensity.^{16}

## 2.8.

### In Vitro Validation

The improvement in performance in tomographic reconstruction was investigated using an agarose phantom mimicking a small animal model. The baseline pattern set employed in this tomographic study are binary bar-shaped patterns described in Ref. 6. The performance of this pattern set in a nonplanar geometry is tested using a truncated hemi-cylindrical phantom shown in Fig. 5 with an absorption coefficient of $0.1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ and reduced scattering coefficient of $13\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$. Furthermore, a cylindrical absorber with six times absorption ($0.6\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$) simulating absorbing organs in the small animal torso (e.g., liver) is placed along the central axis of the phantom. The geometry of the phantom coupled with the absorber results in a high dynamic range across the surface in transmittance.

Two 10 mm long capillary tubes (1.5 mm inner diameter) containing 14 pmol of Cardiogreen (Sigma Aldrich, MO) in 10 *μ*L phosphate-buffered saline were placed at depths of 11 and 9 mm. Moreover, the deeper inclusion was placed below the absorber simulating the occlusion of fluorescent markers by absorbing organs *in vivo*.

The tomographic imaging session employed 36 bar patterns ($34\times 24\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$) described previously. The temporal measurements at excitation wavelength (${\lambda}_{\mathrm{ex}}=780\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$) were acquired at 40 ps interval spanning 2 ns (50 gates) using 300 ps gates ($\mathrm{MCP}\text{\hspace{0.17em}}\text{Voltage}=570\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{V}$; $\text{integration time}=25\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{ms}$). The measurements at excitation wavelength for 36 patterns were acquired in 3 min. The fluorescence measurements (at ${\lambda}_{\mathrm{em}}=832\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$) were acquired for the optimized pattern set ($\mathrm{MCP}\text{\hspace{0.17em}}\text{Voltage}=570\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{V}$; $\text{integration time}=500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{ms}$). The acquisition of measurements spanning 3 ns (75 time-gates) was completed in 18 min. The fluorescence measurements using the uniform pattern set was also acquired for comparison of reconstruction performance.

The tomographic reconstruction of the effective quantum yield of the two inclusions in three-dimensional (3-D) was performed using measurements at 231 point detectors uniformly sampling the surface of the phantom at 2 mm intervals. The time-gate datatype derived at each detector was defined by the number of photons measured relative to the peak photon count and was used to construct the time-gated measurement vector. The effective quantum yield was reconstructed by solving the inverse problem with measurements derived from multiple-gates spanning the rising portion of the TPSF. Furthermore, an x-ray CT image of the phantom was acquired following the optical imaging protocol using small animal CT imaging system (Scanco VivaCT40). The CT data was acquired at an isotropic resolution of 0.038 mm. The images were subsequently resampled to 1 mm isotropic resolution and segmented into three regions—agarose, absorber, and fluorophore—using Amira (VSG US, Burlington, Massachusetts). It should, however, be noted that beside the geometrical boundary information, no *a priori* information on the location of the fluorescence inclusion was enforced in the inverse problem. The 3-D image of the volume (obtained by combining the three regions specified above) was used to generate a voxel-based model geometry with 1 mm voxels. The 3-D model was registered with the optical measurements using fiducial markers positioned outside the phantom. The fully segmented model was subsequently used for validation of the tomographic performance of the pattern optimization scheme.

The dimensions of the 50% isovolume of the reconstructed quantum yield and the position of the centroid of the reconstructed inclusions are used as the resolution metrics to compare the adaptive pattern optimization scheme with the uniform pattern set.

## 3.

## Results

## 3.1.

### Transmittance Optimized Patterns

The optimization of photons transmitted through the synthetic animal model using the full-field pattern shown in Fig. 4 was completed in six iterations for the bar pattern, seven iterations for the wavelet pattern, and eight iterations for the sinusoidal pattern. Figure 6(a) and 6(b) shows the change in the two stopping criterion parameters over the optimization process and for all patterns the process was terminated due to less than 10% change in pattern characteristics. Figure 7(a), 7(c), and 7(e) shows the iterative optimization in the pattern spatial intensity distribution and the corresponding change in transmitted signal intensity on the surface for each of the three patterns. It should be noted that the gray levels on the pattern are adjusted to match the geometry and organ distribution in the animal model with the maximum intensity retained in the thickest portion of the body comprising of the highly absorbing organs like the liver and spleen.

Figure 7(b), 7(d), and 7(f) shows the transmitted excitation field measured on the surface of the animal model. The excitation measurements obtained using the baseline patterns indicate a dynamic range of $\sim 4$ orders of magnitude across the area of illumination due to high fluence near the edge of the model and negligible transmission through the central portion. Following the characteristics of the optimized patterns described above, the optimization process leads to $\sim 2$ orders of magnitude improvement in transmitted light for the central portion of the torso reducing the effective dynamic range in the excitation field to less than 2 orders of magnitude.

## 3.2.

### Fluorescence Signal Optimization

In order to test the hypothesis that the optimization of transmitted excitation field improves the fluorescence measurements, we considered the maximum photon fluence at each inclusion shown in Fig. 4. Figure 8 shows the relative change in the number of photons generated by each of the two inclusions during the optimization process. It is worth noting that the adaptive illumination intensity resulted in more than 600% increase in signal generated at inclusion 1 while retaining the signal level from inclusion 2 at the original level. A comparison of the fluorescence signal measured on the surface at an early time gate (10% of peak) and the maximum gate in Fig. 9(a) to 9(f) demonstrates the effect of nonoptimized measurements in fluorescence tomography with significantly lower signal from deeply embedded fluorescent sources like inclusion 1. Additionally, the optimization process resulted in an improvement in signal from inclusion 1 at both gates suggesting the advantages of this approach for fluorescence tomography. Specifically, the early gate measurements were improved with more than 1.5 order of magnitude increase in signal from inclusion 1 compared with $\sim 1$ order of magnitude increase at the maximum gate. Moreover, the signal from inclusion 2 is not modified for both gates considered and retains the high signal level obtained using the original patterns. The effect of signal improvement for tomographic imaging was quantified by measuring the number of detectors which had a useful signal level (above threshold of 200 counts) at both gates [c.f. Fig. 9(g) to 9(h)]. For all patterns considered, the optimization process resulted in $\sim 60\%$ increase in number of detectors at the early gate and $\sim 75\%$ increase for the maximum gate.

An important characteristic of the measurements shown in Fig. 7 is the near equivalence of the signal obtained using all three types of patterns considered. Specifically, the improvement in signal resulting from pattern optimization (and the corresponding increase in useful detector readings) is the same for all patterns indicating the general applicability of this approach for fluorescence imaging. It suggests also that, due to the highly diffusing and heterogeneous nature of the preclinical models, different patterns will have similar performances when employed in transmittance.

## 3.3.

### Tomographic Reconstruction Using Optimized Patterns

The experimental validation of AWFT using the small animal phantom was performed by batch optimization of the 36 patterns in the baseline set. The iterative scheme was terminated in 4 iterations with the power amplification factor for pattern 31 reaching a value of 0.65 as shown in Fig. 10(a) while the change in pattern characteristics [Fig. 10(b)] was significantly higher than stopping condition, indicating a successful convergence of the optimization scheme.

All patterns were simultaneously optimized in $\sim 3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{min}$ (including $\sim 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{min}$ of pattern setup time per iteration). The optimization procedure of the entire pattern set (4 iterations), including the acquisition of temporal measurements was completed in 24 min (note that the full time-gated data was acquired for each pattern at every iteration). Figure 11(a) and 11(b) compare the structure of a baseline half-field pattern and the corresponding transmittance-optimized grayscale patterns. It can be inferred that the optimization procedure removed the sections of pattern beneath the region with smaller thickness to correct for the dynamic range while retaining the high intensity across the central axis to compensate for the absorber. Figure 11(c) and 11(d) compare the excitation field acquired for the original and optimized patterns for the maximum gate. As expected, the original pattern results in a higher signal towards the edges of the phantom [c.f. Fig. 11(c)] demonstrating the role of model geometry as a source of high dynamic range. Conversely, the transmittance optimized field provides a high signal across the area of the pattern [c.f. Fig. 11(d)] with maximum of 4000 photons across the area of the pattern representing the successful optimization of the transmitted signal.

A similar comparison of the fluorescence field shown in Fig. 11(e) to 11(j) demonstrates that the optimization of the excitation field improves the fluorescence signal from inclusion 1. This demonstrates an improved contrast and higher signal over a larger number of detectors when using AWFT. The results of optimization of all the patterns used in the study are shown in Video 1.

## 3.4.

### Improved Information Content

The impact of the improvement in the overall signal level demonstrated in the measurements on the quality of the tomographic information recorded is quantified using a two-step process. First the percentage increase in number of source-detector pairs at different gates along the temporal measurement which are above the noise threshold level is used to quantify an increase in the number of usable measurements The noise threshold for this system was empirically determined to be 200 photon counts. Figure 12(a) and 12(b) shows $\sim 200\%$ increase in the number of s-d pairs for excitation signal and $\sim 80\%$ in the fluorescence signal at the early gates, with a smaller yet significant improvement for the maximum gate ($\sim 25\%$ for fluorescence field, 15% for excitation field). The improvement in the time-gated signal at the early gates and the late gates establishes the advantages of this imaging technique for high-resolution optical reconstruction^{23}24.^{–}^{25} and lifetime multiplexing applications, respectively.^{14} It is also worth noting that the transmittance optimized patterns result in a higher increase for the late gates when compared with the early gates as the optimization procedure corrects for the effects of the highly absorbing local perturbation in the experimental phantom which reduces signal at the late gates.

Next, the analysis of tomographic information content in the measurements is extended to the quantification of increase in nonredundant information using a model based analysis. Following the work by Freiberger et al., we define the orthogonality of the different rows in the fluorescence weight matrices (obtained by discretizing the fluorescence sensitivity function ${J}^{\mathrm{fluo}}$) as a measure of the information content.^{26}

For a set of $N$ patterns, let $M$ be the total number of measurements (s-d pairs) above the set threshold. Consider a set of measurements ${S}_{n}$ for the $n$th pattern. Then the nonredundant information ${h}_{n}$ for the $n$th pattern is defined as the average orthogonality of each measurement in ${S}_{n}$ with respect to the remaining measurements ($M\setminus {S}_{n}$). ${h}_{n}$ is computed using Eq. (4):

## (4)

$${h}_{i}=1-\frac{1}{|M\backslash {S}_{i}||{S}_{i}|}\sum _{n\in M\setminus {S}_{i}}\sum _{m\in {S}_{i}}\frac{{({j}_{m},{j}_{n})}^{2}}{{\Vert {j}_{m}\Vert}^{2}{\Vert {j}_{n}\Vert}^{2}}.$$Here, ${j}_{m}$ and ${j}_{n}$ are the $m$th and $n$th rows from ${J}^{\mathrm{fluo}}$. The average nonredundant information content for the entire pattern set is therefore given by

The average information content value of 1 represents maximal nonredundancy in the measurement suggesting maximum information content in the measurements. Figure 12(c) shows the comparison of average nonredundant information content for the uniform and the transmittance-optimized pattern set at individual time-gates along the TPSF.

It is observed that the optimization of transmittance excitation field leads to $\sim 6\%$ increase in the information content at the early gates and $\sim 16\%$ increase for gates along the decaying section of the TPSF. This closely follows the trend in increase in number of s-d pairs observed in Fig. 12(a).

## 3.5.

### Reconstruction Results

In order to determine the impact of increased information content on the quality of reconstructions in the image space, the effective quantum yield was reconstructed using three time gates spanning the rising portion of the TPSF (early gates at 20%, 50%, and the maximum gate). These gates were selected on the rising part of the TPSF as it correlates to resolution improvement in the reconstruction, whereas late gates are critical in lifetime multiplexed studies, which we do not consider herein.^{16}^{,}^{27} Figure 13(a) and 13(b) shows the 50% iso-volumes of the effective quantum yield obtained by solving the inverse problem by least-squares minimization. The lower number of source-detector pairs available when using the uniform pattern set (more than 36% reduction) leads to reduced tomographic information in the measurements as discussed previously. This results in poor discrimination of the two objects in the reconstruction and their incorrect localization [c.f. Fig. 13(a)]. Specifically, both objects are reconstructed 5 mm above the expected position along the $z$-axis. Moreover, the lower resolution of the reconstructed volume is evidenced by the connected reconstructed objects. On the contrary, the reconstructed volume obtained using the transmittance optimized patterns allows us to discriminate and accurately localize the two inclusions as shown in Fig. 13(b). The resolution metrics comparing the reconstructed 50% isovolumes with the ground-truth values obtained from the x-ray CT are compiled in Table 2. It is worth noting that the two inclusions are accurately localized with $\sim 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ error.

## Table 2

Comparison of resolution metrics of reconstructed quantum yield with ground truth values derived from x-ray CT.

Inclusion # | x | y | z | ||||
---|---|---|---|---|---|---|---|

Expected | Reconstructed | Expected | Reconstructed | Expected | Reconstructed | ||

Dimensions (mm) | 1 | 9.9 | 9.0 | 1.6 | 4.0 | 1.6 | 8.0 |

2 | 9.8 | 8.0 | 1.6 | 2.0 | 1.6 | 3.0 | |

Position of centroid (mm) | 1 | 20.7 | 21.0 | 16.4 | 16.0 | 8.3 | 9.0 |

2 | 20.7 | 22.0 | 24.7 | 25.0 | 8.3 | 10.0 |

Also, the separation between the two objects was found to be 9 mm (expected—8.3 mm) resulting in less than 1 mm localization error. Furthermore, the reconstructed dimension of the inclusions had $\sim 1$ and $\sim 2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ error along the $x$-axis and $y$-axis, respectively. The resolution of inclusion 2 along the $z$-axis also had less than 2 mm error; however, inclusion 1 had significantly lower z-resolution due to the limited projection angles in the planar imaging configuration. Note that the imaging volume was discretized with $1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ voxels, and hence the reported errors are on the same scale as the smallest element of volume (voxel).

## 4.

## Discussion

The optimization of excitation patterns for more accurate reconstruction in wide-field optical tomography of small animals is an emerging area of investigation owing to the experimental advantages provided by wide-field imaging techniques. Conversely to model-based optimization techniques, in this work we investigated a novel measurement-guided pattern optimization technique with the objective of reducing the dynamic range in intensity of photons transmitted through the small animal model. The iterative pattern correction scheme described in this work was based upon the characteristics of the transmitted excitation field allowing a fast optimization of excitation patterns during acquisition. Furthermore, the associated improvement in fluorescence signal by proxy upon the optimization of the transmitted field established AWFT as a useful technique in FMT. It is worth noting that the acquisition-time optimization of patterns in AWFT introduced minimal time cost to the experimental protocol. Hence it can be readily implemented for *in vivo* studies to acquire optimal data sets by taking into account at each imaging session the geometry of the specimen imaged (animal and posture specific) and its optical parameter characteristics (e.g., spatial distribution of organs, which is dependent on posture, functional state, disease progression).

The *in silico* validation of this approach in a synthetic animal model demonstrated an increase in transmitted photon intensity by $\sim 2.0$ orders of magnitude for all patterns considered. Furthermore, the adaptive correction of excitation pattern was simultaneously able to improve the signal from both fluorescent inclusions located in the animal torso. Specifically, the fluorescent signal from inclusion 2 located in the center of the torso occluded by highly absorbing organs (e.g., liver and spleen) was increased by approximately two orders of magnitude after pattern optimization. More importantly, the optimization scheme provided equivalent results for all three types of patterns investigated here, indicating the applicability of this approach to any type of illumination pattern. The increase in information content with the increase in s-d pairs allows the accurate reconstruction of deep-seated fluorescent inclusions, which are poorly localized when using the nonoptimized patterns.

The application of this technique in FMT was found to especially improve the signal level at the time gates on the early-rising and late-falling portion of the TPSF that correlates with lower photon counts. The corresponding 80% increase in the number of source-detector measurements above the noise threshold acquired in experimental settings demonstrates the advantage of this approach in FMT applications employing early photons for high-resolution optical reconstructions.^{24} Similarly, the improvement in signal level at the late gates implies more robust measurements in fluorescence lifetime based tomographic imaging where multiple fluorophores can be resolved with lower crosstalk owing to improved signal-to-noise ratio.^{28}^{,}^{29} Moreover, the advantages of pattern optimization in FMT also extend beyond the improvement in tomographic information content, in that the improved signal level from deep-seated fluorescent markers will allow the robust estimates of fluorescence lifetime due to the higher SNR of the measurements.

The application of the pattern optimization technique described in this paper is directed toward time-resolved fluorescent tomographic imaging. However, it is worth noting that this approach can be readily translated to other optical tomographic imaging paradigms. For instance, by assuming equivalent characteristics of the continuous-wave datatype and the time-gate at maximum intensity, it can be determined that the pattern optimization scheme also provides added information in tomographic systems employing the CW datatype. Similarly, the optimization scheme can also be an effective imaging approach in systems implementing a cylindrical imaging geometry or a multiview setup based on mirrors where dynamic range is a critical issue in experimental implementation.^{30} We further hypothesize that the approach described here can also be extended to point source excitation schemes, wherein the optimized pattern derived from a full-field illumination pattern can be used as a template to assign the relative power injected into the model at each source location.^{31}

Lastly, the iterative optimization scheme also has applications in related optical imaging applications. For instance, in the case of fluorescence reflectance imaging, which is a commonly used preclinical imaging modality, a measurement guided optimization of excitation field based upon the fluorescence signal detected may improve the depth sensitivity of the technique.

In conclusion, we have for the first time described a measurement-guided pattern optimization scheme, which increases the tomographic information collected in whole-body small animal imaging applications. The preliminary results presented here establish the feasibility and the applicability of the technique. The future work will focus on the impact of specific parameters of the optimization scheme (e.g., variable neighborhood radii) on the convergence of the technique. We next plan to extend the approach to include optimization of patterns based upon the fluorescence field. Pattern optimization based on this criterion will allow the detection of lower concentrations of fluorophores, increasing the sensitivity of the technique. We anticipate the use of measurement-guided excitation optimization schemes to further improve the applicability of FMT in molecular imaging applications.

## Acknowledgments

This work was supported by National Institutes of Health (NIH) grant R21 EB013421 and by the National Science Foundation (NSF) CAREER CBET-1149407. The authors would like to thank Ms. Jin Chen for her assistance with Monte Carlo modeling and reconstruction. We would also like to thank the Center for Nanotechnology Innovations (CCNI) at Rensselaer Polytechnic Institute for their support.

## References

*in vivo*,” Proc. Natl. Acad. Sci. U.S.A. 105(49), 19126–19131 (2008).PNASA60027-8424 http://dx.doi.org/10.1073/pnas.0804798105 Google Scholar

*Ex vivo*fluorescence molecular tomography of the spine,” Int. J. Biomed. Imag. 2012, 1–11 (2012).IJBIBD1687-4196 http://dx.doi.org/10.1155/2012/942326 Google Scholar