Comparison of time- and angular-domain scatter rejection in mesoscopic optical projection tomography: a simulation study.

Optical imaging offers exquisite sensitivity and resolution for assessing biological tissue in microscopy applications; however, for samples that are greater than a few hundred microns in thickness (such as whole tissue biopsies), spatial resolution is substantially limited by the effects of light scattering. To improve resolution, time- and angular-domain methods have been developed to reject detection of highly scattered light. This work utilizes a modified version of a commonly used Monte Carlo light propagation software package (MCML) to present the first comparison of time- and angular-domain improvements in spatial resolution with respect to varying sample thickness and optical properties (absorption and scattering). Specific comparisons were made at various tissue thicknesses (1-6 mm) assuming either typical (average) soft tissue scattering properties, μs ' = 10 cm-1, or low scattering properties, μs ' = 3.4 cm-1, as measured in lymph nodes.


Introduction
Optical imaging is used extensively for assessing biological tissue specimens (biopsies) through tissue staining and microscopy of thin tissue slices taken from selected locations of the specimens. However, the time intensive nature of these procedures limits more complete assessment of specimen volumes. Lymph node biopsy, used to stage metastatic cancer in a number of tumor types, is one example of a clinical application that would benefit from imaging strategies capable of evaluating larger volumes of tissue, as clinically important micrometastases can be missed with current protocols [1]. Optical imaging with photons in the near-infrared window (600-1000 nm)-where the absorption is low and scattering dominates-offers the potential to carry out whole specimen analyses for < 1 mm diameter, low-scattering tissues [2], and for larger tissues if scattering can be accounted for [3][4][5]. Two of the more common methods of scatter rejection to improve resolution in imaging of scattering biological tissue include (1) time-domain imaging (also called early-photon imaging) [6][7][8][9], where pulsed light sources are used and images are reconstructed only from the photons that arrived earliest at the detector (having taken the most direct path through tissue); and (2) angular-domain imaging [10][11][12], where only the straightest photons exiting the tissue are collected (presumably limiting collection of photons that have deviated far from the most direct path through the tissue).
In both approaches, there are a number of instrument design considerations that can affect level of scatter rejection (i.e., time resolution in time-domain applications and numerical aperture of the collection optics in angular-domain applications). Yet, while greater levels of rejection (smaller the time window or smaller the aperture) can theoretically provide better spatial resolution, there is a tradeoff in resolution and collected number of photons. As such, the question of how to reject scattered photons and to what extent is complex and will depend on the conditions of the application (size and optical properties of the specimen). This work presents the development of a simulation tool to directly compare time-and angular-domain imaging at different levels of time-and angular-scatter rejection. Specific comparisons were made at various tissue thicknesses assuming either typical (average) soft tissue scattering properties (μ s ' = 10 cm −1 [13]) or low scattering properties (lymph node: μ s ' = 3.4 cm −1 , measured in this work, see Section 2.4 & 3.1). Simulations were carried out with an augmented version of a commonly used Monte Carlo simulation subroutine, mcsub.c [14,15], which was optimized for GPU parallelization and modified such that the path and exit angles of the light that "hit" the designated detector was saved separately from the overall fluence map. This was critical for comparing time-and angular-domain resolution and contrast tradeoffs. The augmented Monte Carlo code is available to anyone upon a request.

Monte Carlo simulation
To compare potential resolution improvements using time-and angular-domain imaging in biological specimens (1-6 mm in diameter), the Monte Carlo (MC) subroutine, mcsub.c, developed and distributed by Jacques et. al [14,15] was augmented in a few ways. Briefly, mcsub.c was designed to mimic the light propagation in finite sized, turbid (scattering) media with homogeneous optical properties, assuming Fresnel reflection at tissue/air interfaces. The code requires inputs of the size of the medium (s), the absorption coefficient of the medium (μ a ), the scattering coefficient of the medium (μ s ), the anisotropy of scattering (g), the refractive index of the medium (n 1 ), the refractive index of the external medium (n 2 ), and characteristics of the incident light beam (e.g., light shape and position of illumination). Briefly, in MC simulation, photon "packets" (which are commonly called photons in MC literature) are launched at a location and in a direction and distribution matching the desired illumination setup. Each packet is made to travel a random distance (step) before scattering (changing direction). The distribution of random distances is linked to the probability of scattering (μ s ) and the distribution of random direction changes is linked to the anisotropy of the scattering (g). At the end of each step, the weight of the photon packet (which begins as 1) is adjusted, according to the absorption properties (μ a ). The process is then repeated until the photon packet leaves the medium, is detected, or reaches a selected lower limit of weight. The mcsub.c program then outputs a vector of escaping flux density versus radial position, ( ) i J r , and a matrix of fractional density map of incident light transported ( ) , i k T r z , recorded on cylindrical coordinate system. In order to make mcsub.c amenable to compare time-and angular-domain imaging in "small" (mesoscopic) tissue samples, and improve execution time, the following adjustments were made: (1) Detector subroutine was added such that one can define a detector's acceptance in space (size, location), time, and angle. (comparing time-vs. angular-domain).
(2) The complete path of each photon package was retained through the "life" of the photon packet (until it exited the medium or reached a weight that would trigger a drop). If the photon packet was accepted by the detector (see (1)), the path of that photon packet was added to a separate detected fractional density map of incident (3) The incident Gaussian beam simulation was modified to propagate the beam as a solution of the paraxial Helmhotlz equation in the tissue until the first interaction point, as described by Hokr et al [16]. (improving approximation of a Gaussian beam propagation in turbid medium).
(4) A variance reduction technique defined by Kawrakow and Fippel [17] was employed to reduce the number of photon packets passing through the medium without having any interaction. Specifically, a factor DIV (0<DIV≤1) was introduced, where each photon package was then propagated a distance, Δs*DIV (where Δs is the standard step size derived by mcsub.c). Roulette was then used to determine if the photon packet was scattered. (improving execution time).
(5) The code was optimized for parallelization on a GPU cluster. Specifically, a new random number generator [18] was used to avoid repeats of the same random seed numbers, when photons were initiated on different threads of the GPU at the same CPU clock time (used in the mcsub.c code as a randomizer for the random number generator). (improving execution time).

Simulated phantom experiment
The augmented Monte Carlo software was used to evaluate how the scatter rejection by timeand angular-domain restrictions can influence spatial resolution and contrast in transmission optical projection tomography (OPT) of average and low (lymph node) scattering coefficient tissues over a range of tissue thicknesses: 1-6 mm. Average soft tissue optical properties for 780 nm wavelength light (absorption coefficient: μ a = 0.2 cm −1 , scattering coefficient: μ s = 100 cm −1 , anisotropy factor: g = 0.9) were taken from Jacques [13], and lower scattering lymph node optical properties (μ a = 0. The illumination source was simulated as a Gaussian beam reaching its waist at the surface of the sample and having a full-width-half-maximum waist of 100 μm. The size of the source was selected to minimize the spot size while maintaining a non-divergent beam through the thickness of the samples (divergence increases with smaller spot sizes) [19]. Detector size (100 μm diameter, circular) was selected to match the size of the illumination point since the resolution will be dominated by the light source beam width. Source and detector were simulated to be directly opposite on either side of the scattering media. The tissue index of refraction was set to 1.4 [13] and the surrounding was simulated as air. Under each imaging condition and assuming a small perturbation of the absorption ( Fig.  6(a)), the detected intensity at x-axis location, x m , and y-axis location, y n was estimated by: x y z is the detected fractional density map of transported incident light obtained in Monte Carlo simulations defined on a 10 μm pixel grid over a 2d x 2d x d mm 3 uniform object area, with d denoting the object thickness, 0 N representing the number of launched photons per source-detector location, and δμ representing the perturbation in the absorption coefficients of the scattering medium [20]. The δμ perturbation map was defined as zero at all locations except for two, 200x200x200-μm 3 -square inclusions separated by 200 μm in the central axial plane of the objects. The level of perturbation was assumed to be 0.1 × (μ a + μ s ) -i.e., 10% of the attenuation coefficient for each case. To simulate tomographic data collection, 2D projections were calculated at 10 μm resolution over a d x d field of view with 600 fJ of photons detected per point, for 72 projections about the objects (assuming projections acquired every 5°). Poisson (shot) noise was added with the poissrnd() built-in function in MATLAB to each projection, and data were reconstructed in 3D by standard filtered backprojection using the iradon() built-in function in MATLAB. This setup would be representative of a spherical tissue object rotated in an optical property-matching solution. All tomographic data were presented as reconstructed 2D slices, perpendicular to the axis of rotation and at the center of mass of the inclusions.

Analyzing spatial resolution
To evaluate achievable resolution the object to background contrast was calculated as: 1 2 where, I 1 and I 2 were the mean intensities measured from 50x50 μm 2 regions about the center of mass of one absorber (I 1 ) and in the center of mass of the space between the absorbers (I 2 ), in the reconstructed images. Mean squared error (MSE) between estimated and simulated values was measured to assess the accuracy of the reconstructed images as: where  denotes 500x500 μm 2 region containing inclusions; N is the number of pixels in this region; ( ) Note that the mean background of each reconstructed image should be zero . To evaluate limits of the spatial resolution (independent of imaging time or laser power) both contrast and MSE were calculated using noise free projections. For a discussion about imaging time and ANSI safety limit on laser power see Section 3.2.1.

Optical properties of lymph nodes
Lymph node optical properties were estimated in two steps: 1) using a time-domain optical imaging system, diffusion approximation of the radiative transfer equation, and thick tissue samples to estimate a μ and s μ′ ; 2) Beer-Lambert law and thin tissue sample to estimate g . 1). Estimation of a μ and s μ′ ; Thirty lymph nodes were surgically removed from the neck tissue of freshly slaughtered pigs provided by a local butcher. The nodes were packed into a 4x4-cm 2 clear square glass container. On the same day, the lymph node filled container was placed in the imaging field of a time-domain optical imaging system [19]. Briefly, the system employed a 780 nm femtosecond pulsed laser (Mendocino, Calmar Laser, Palo Alto, CA) for illumination and a single photon avalanche diode detector (PDM, Micro Photon Devices, Italy) connected to a time correlated single photon counting module (HydraHarp, PicoQuant, Berlin, Germany) to measure the temporal point spread function of light arrival at a time resolution of 4 ± 12 ps. The system was used to collect the transit time distribution of photons passing through the 4-cm-thick lymph node sample. At this thickness, it was assumed that the diffusion approximation of the radiative transfer equation was sufficient to model the collected time-domain signal [21]. Specifically, the following expression was used to estimate the average absorption coefficient, a μ , and reduced scattering coefficient, s μ′ , of lymph nodes through least squares optimization using MATLAB (Mathworks, Natick, MA) code: where φ(t) is the measured temporal point spread function of the laser after passing through the lymph node sample, φ 0 (t) is the measured instrument response function of the system, * represents the convolution operator, a is a scaling factor and fitting parameter, d is the thickness of the sample (4 cm), and v represents the speed of light in the medium (with an assumed index of refraction of 1.4 [13]).
2). Estimation of g ; Lymph nodes were then frozen in TissueTek OCT Compound (Sakura Finetek, Torrance, CA) and cryostat sectioned in d = 100-µm thick slices. The samples were mounted using a wet cell geometry as developed by Hall et. al [22] to prevent dehydration of the tissues. Briefly, samples were placed on 1-mmthick glass slides with phosphate buffered saline above and below the tissue, and topped with a coverslip. The slide was sealed using Vaseline. Thin tissue slices permitted the condition of single scattering, and thus, use of the Beer-Lambert law to model transmission: where P is the measured transmitted power after passing through the lymph node sample, P 0 is the measured incident power, d is the tissue thickness, and µ t is the total attenuation coefficient comprised of the absorption and scattering coefficients (µ t = µ a + µ s ). Each tissue section was illuminated with the 780 nm femtosecond laser in two spots -one in the center of the sample and the other on the edge to account for variability in the biological structure of lymph nodes and corresponding optical properties [23]. Incident and transmitted power were measured with a photodiode power sensor (S120C, Thorlabs, Newton, NJ). The scattering anisotropy factor was calculated using g = 1 -µ s '/µ s . To verify the methods described above, the procedure was carried out on optical tissuemimicking phantoms (Biomimic, INO, Quebec, Canada) where optical properties were known (µ a = 0.05 cm −1 , µ s ' = 10 cm −1 , g = 0.62, n = 1.51). The thick tissue diffuse approximation experiment was conducted with a 3x3x4.8 cm 3 block of the phantom, while the thin tissue attenuation experiment was done on a 3x3x0.4 cm 3 slab.

Lymph node optical property estimation
The absorption and reduced scattering coefficients were obtained from fitting measured temporal pulse spread functions with the diffusion approximation to the radiative transfer equation [Eq. (4)]. The estimated optical property values were µ a = 0.30 cm −1 and µ s ' = 3.37 cm −1 . To further characterize the lymph nodes, transmittance measurements were acquired from thin tissue samples, and Eq. (5) was used to calculate µ t and subsequently obtain µ s so that the anisotropy factor could be attained. The average attenuation coefficients were µ t = 46.4 ± 17.2 cm −1 and µ t = 40.3 ± 19.6 cm −1 at the edge and center of the lymph node sections respectively. Higher values at the edge compared to the bulk are consistent with findings from an in-depth analysis of local attenuation coefficient in human lymph nodes conducted by Scolaro et al. [23]. It is expected that regions of fibrous stroma from the capsule (periphery of the node) have higher attenuation, while the paracortex and medullary sinuses (bulk of the node) will have lower attenuation; and the obtained results follow this trend. The above values are also comparable to the coefficients presented in the literature, which were found to be in the range 45 -153 cm −1 at 1320 nm using optical coherence tomography [23]. Note: lymph nodes primarily consist of lymph fluid, fibrous tissue, and fatty tissue, which all have scattering coefficient values that are relatively constant between 780 and 1320 nm [13]. Furthermore, Scolaro et al. assumed negligible absorption coefficient, yet water is 4 ordersof-magnitude more absorbing at 1320 nm compared to 780 nm, which could have yielded slight overestimates in lymph node scattering coefficient estimates at 1320 nm. A mean of µ t = 43.3 cm −1 was used for calculations and provided a resultant anisotropy factor of 0.92again consistent with other values of g in the literature for biological soft tissue [13].
The optical properties of a phantom with known absorption coefficient, reduced scattering coefficient, and anisotropy factor were characterized to confirm the validity of the experimental methods. The obtained values were µ a = 0.055 cm −1 , µ s ' = 9.377 cm −1 and g = 0.636, providing an error of 10%, 6.2%, and 2.6% for the absorption coefficient, reduced scattering coefficient and anisotropy factor, respectively (target values: µ a = 0.05 cm −1 , µ s ' = 10 cm −1 and g = 0.62-provided from the manufacturer). These coefficient values fall within the expected margin of error that results from a time resolved transmittance characterization method (11.3% for absorption coefficient and 6.8% for reduced scattering coefficient [24]). Although the expected error for g is not reported, the obtained result is in good agreement with the target value, and within 2.6% of the expected value (g = 0.62 ± 0.015). Based on these findings, the validity of the aforementioned techniques to characterize optical properties was confirmed, and the estimated lymph node optical properties could then be applied to simulation studies with confidence.

Simulated phantom results
First, in Fig. 1, the detected fractional density map of transported incident light recorded on a cylindrical coordinate system, normalized to maximum of one, is reported for the studied range of time-and angular-domain restriction for the low-scattering tissue (i.e., lymph node like) and thicknesses of 3-6 mm. In general, angular-domain restriction provided the highest improvement in potential spatial resolution as one can determine from the density map profile at the half thickness of the object, relative to the object thickness, with improvements over no restriction diminishing with increasing object thickness. On the other hand, time-domain restriction tended to provide greater improvements in potential spatial resolution for thicker samples. Though results are not provided in Fig. 1, if one would increase scattering properties, µ s , or lower scattering anisotropy, g, similar trend can be observed as with increased sample thickness.

ANSI safety limits
One limitation to scatter rejection approaches that requires consideration is the potential for loss of signal-to-noise ratio as restrictions become more substantial. Based on the light source and detector sizes simulated in this work, the amount of energy that would be needed to detected 600 fJ of light energy, i.e. ~2.4 × 10 6 photons at 780 nm, (necessary to clearly see 10% perturbations in optical properties for the most restrictive OPT case evaluated in this work: 0.005 NA angular restriction at 6 mm diameter lymph node object-see Fig. 3) at the detector for different lymph node-like tissue thicknesses is presented in Fig. 2. Assuming a dwell time of 1 ms (< 15 min tomographic imaging at 5x5 mm fields-of-view at 50 μm steps about 72 projections), the laser power required to achieve 600 fJ at the detector in the most restrictive case evaluated (angular restriction of NA = 0.005) and at the thickest distance evaluated (6 mm) would be 60 mW. When this power is focused to a 100-μm diameter spot, using for an example a 10 MHz pulsed laser, the energy deposited per pulse would be less than 0.1 mJ/cm 2 , 2 orders-of-magnitude below the ANSI skin safety limit for a single pulsed near-infrared laser (30 mJ/cm 2 ) [25]. While repeated pulsing as required in the proposed setup can amplify the potential for tissue damage, it should be noted that in an experiment evaluating cellular reproduction of exposed living hamster embryos over 24 h of exposure at 8 orders-of-magnitude higher light power density [26] did not produce appreciable effects. It should be noted, that for angular restriction applications, there is no need to use a single pixel detector, but rather a CCD camera could be used for detection, allowing imaging times to reduce significantly (unpublished results). Even for time-domain applications, optical gating allows the use of CCD detectors to acquire many projection simultaneously [27], and there is work to advance time-correlated single-photon counting SPAD arrays [28].

Visual
Images for a reconstruction 600fJ), for a object with av without any s tissue simulat it possible to scattering cas the ability of mesoscopic ra identification whereas, sign mm thick lym 2 [30,31].). ane images of 3 200x200x200 μm ted, separated by 2 kness of object) fo t at each detector l for one noisewn in under all co erage soft tissue ( a 6 mm thick lym first column of im nts images for an an angular-doma time-domain restr iction of 0.1 ps. Un st ported density node-like biolo n Eq. (2). Figu gular-domain sc sue thicknesse sults further su effective at im ments diminish l possible phot gime, which is re specifically contrast compa ut 2 mm for av ed that improv mm ( Fig. 4(b)).
al shape with ng the potenti cantly improve ect 200-μm dia the smallest cl 3D scatter-rejectio m 3 inclusions wit 200 μm. Projection or 72 evenly spaced location was assum realization using onditions. The top (μ a = 0.2 cm −1 , μ s mph node like tissu mages is for the c angular restrictio ain restriction of riction of 1 ps, and nits are in cm −1 (at maps of vario ogical tissues ure 4 presents t catter rejection s of 1-3 mm d upported the c mproving spatia as the thickne ton paths are typically assu , significant an ared to more sta verage tissue (F vements in con This is signif short-axis dia al to use rela e spatial resolu ameter objects linically releva on optical projec th a 10% increa ns were simulated d projections abou med to be 600 fJ. P g standard filtere p row of images p s = 100 cm −1 , g = ue (μ a = 0.3 cm −1 , case with no scat on of NA = 0.059, NA = 0.005, th d the fifth column ttenuation). ous optical proj were further a the comparison n methods in te diameter and in conclusions fro al resolution fo ess and/or the l equally proba umed to be aro ngular restricti andard imagin Fig. 4(a)); whe ntrast with angu ficant, as the v ameters of up atively inexpe ution in whole is significant a ant metastases ction tomography ase in attenuation d at 10 μm spacing ut each object. The Poisson noise was ed backprojection pertains to results = 0.9) object. The μ s = 43 cm −1 , g = tter-restriction, the , the third column he fourth column represents images jection setups ( analyzed in te n of normalized erms of detectin n lymph node t om Fig. 1: that or thinner and evel of scatteri able for any ac ound 1-2 cm fo ion (NA = 0.0 ng (e.g., NA = 0 ereas in lower s ular restriction vast majority o to 6 mm (re ensive angular e lymph node as this is the si in lymph nod protocols to enhance cancer cell contrast through targeted fluorescent of other optical imaging agents, and the long term ability to detect micrometastases in excised human nodes will require task-based evaluation metrics to demonstrate improved sensitivity and specificity over existing clinical protocols. Laminar optical tomography (LOT) [35] and mesoscopic fluorescence molecular tomography (MFMT) are competing imaging strategies to the mesoscopic OPT methods explored in this work, permitting 3D visualization of biological tissues at depths of 2-3 mm with high resolution (~100-200 µm). Their application has been demonstrated for in vivo studies of tumor [36,37] and brain activity [38] in mice, and evaluation of engineered tissues [39,40]. These imaging methods have the advantage of being capable of imaging the top surfaces of much thicker tissues compared to transmission OPT. The advantage of OPT, demonstrated in this paper, is that it is possible to achieve similar spatial resolution by simple filtered backprojection imaging reconstruction; hence, it is not necessary to utilize modelbased reconstructions that can be timely, particularly for LOT/MFMT, which usually employs modeling of photon propagation via Monte Carlo simulation (yet it should be noted that data reduction techniques are being explored to account for the computational expense in solving the inverse problem [41]).
A broader impact of this work is the development and dissemination of a Monte Carlo software optimized for mesoscopic imaging. The nuances demonstrated in the work between optical properties and methods of scatter rejection with respect to their effects on spatial resolution in mesoscopic optical imaging reinforce the complexity of optical imaging in this regime. They also highlight the critical role that an accurate, fast, and informative photon propagation simulator can play in further developing and optimizing the future of mesoscopic optical imaging systems. The software utilized and provided by this work offers such ability, as an augmented version of the Monte Carlo subroutine employed by the "MCML" software widely used for visible and near-visible photon propagation modeling in scattering biological tissues [42]. Specifically, the code was (1) optimized for GPU-parallelization, (2) amended to allow the paths of photons arriving at a detector within definable time windows and/or at definable angles of incidence to be separated from the paths of all other photons, and (3) amended to allow modeling of coherent propagation of a Gaussian light source prior to the first interaction point. While neither GPU parallelization of MCML [43], addition of time and angular constraints [11], nor Gaussian propagation [16] in photon propagation are on their own new, the presented code is the first code to combine all of these optimizations, providing unique capability to accurately compare time-and angular-domain systems, which encompasses the immediate future of advances in optical projection tomography in the mesoscopic range.

Disclosures
The authors declare that there are no conflicts of interest related to this article.