1. 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–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–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–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⁻¹ [13]) or low scattering properties (lymph node: μ_s’ = 3.4 cm⁻¹, 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.

2. Method

2.1 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₁), the refractive index of the external medium (n₂), 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, $J\left(r_i\right)$, and a matrix of fractional density map of incident light transported $T\left(r_i,z_k\right)$, 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:

2.2 Simulated phantom experiment

The augmented Monte Carlo software was used to evaluate how the scatter rejection by time- and 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⁻¹, scattering coefficient: μ_s = 100 cm⁻¹, anisotropy factor: g = 0.9) were taken from Jacques [13], and lower scattering lymph node optical properties (μ_a = 0.3 cm⁻¹, μ_s = 43 cm⁻¹, g = 0.92) were estimated experimentally (Sections 2.4 & 3.1). For each optical property group and every thickness, four levels of angular-domain restriction were tested (numerical apertures (NA) of 0.124, 0.059, 0.025, and 0.005), including all photon arrival times, and at NA = 0.124, three time-domain restrictions (also called time gates) were tested (1, 0.3, and 0.1 ps). 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:

2.3 Analyzing spatial resolution

To evaluate achievable resolution the object to background contrast was calculated as:

Mean squared error (MSE) between estimated and simulated values was measured to assess the accuracy of the reconstructed images as:

2.4 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 ${\mu }_a$ and ${{\mu }^{\prime }}_s$; 2) Beer-Lambert law and thin tissue sample to estimate $g$.

To verify the methods described above, the procedure was carried out on optical tissue-mimicking phantoms (Biomimic, INO, Quebec, Canada) where optical properties were known (µ_a = 0.05 cm⁻¹, µ_s’ = 10 cm⁻¹, g = 0.62, n = 1.51). The thick tissue diffuse approximation experiment was conducted with a 3x3x4.8 cm³ block of the phantom, while the thin tissue attenuation experiment was done on a 3x3x0.4 cm³ slab.

3. Results and discussion

3.1 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⁻¹ and µ_s’ = 3.37 cm⁻¹. 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⁻¹ and µ_t = 40.3 ± 19.6 cm⁻¹ 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⁻¹ 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 orders-of-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⁻¹ was used for calculations and provided a resultant anisotropy factor of 0.92—again 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⁻¹, µ_s’ = 9.377 cm⁻¹ 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⁻¹, µ_s’ = 10 cm⁻¹ 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.

3.2 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.

 

Fig. 1 Detected photons’ transported density maps determined by Monte Carlo simulation plotted as a function of axial distance (z-axis; along the direction of the illumination vector) and the radial distance (x-axis; perpendicular to the direction of the illumination vector). This figure provides a subset of results for the low scattering tissue (similar to lymph node tissue; μ_a = 0.3 cm⁻¹, μ_s = 43 cm⁻¹, g = 0.92), for object thickness between 3 and 6 mm, angular-restriction between NA = 0.005-0.124, and time-domain restriction for 0.1-1 ps. All simulations were based on a 100-μm FWHM Gaussian light source incident on the object at its waist and a 100-μm diameter detector.

Download Full Size | PDF

3.2.1 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⁶ 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 orders-of-magnitude below the ANSI skin safety limit for a single pulsed near-infrared laser (30 mJ/cm²) [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].

 

Fig. 2 Monte Carlo simulation results estimating the amount of light energy required at illumination to detect 600 fJ of energy (~2.4 × 10⁶ photons at 780 nm) at various thickness of a lymph node-like tissue (μ_a = 0.3 cm⁻¹, μ_s = 43 cm⁻¹, g = 0.92), assuming a 100-μm FWHM Gaussian light source incident on the object at its waist and a 100-μm diameter detector at three example restrictions along with the no restriction case.

Download Full Size | PDF

3.2.2 Visual comparison of the tomographic images with the imaging noise

Images for a visual comparison are shown in (Section 2.2) in Fig. 3. Here we show 3D reconstructions, from the projections with ~2.4 × 10⁶ detected photons per pixel (equivalent to 600fJ), for a 6 mm thick object with lymph node-like optical properties and a 3 mm thick object with average optical properties for each scatter rejection method. One can observe that without any scatter rejection, it was not possible to resolve the objects in the 3 mm average tissue simulation; however, the low scattering nature of the lymph node simulation did make it possible to “see something” even with no scatter rejection for a 6 mm thick object. In both scattering cases, a 0.1 ps time gate improved observation of the inclusion, further supporting the ability of narrow window (optical gating) for improving optical imaging resolution in the mesoscopic range [27]. Angular restriction with an NA = 0.005 was not sufficient to improve identification of the inclusions for the higher scattering 3 mm thick average tissue-like object; whereas, significant improvements in identification were observed for the lower scattering 6 mm thick lymph node like object.

 

Fig. 3 Transverse plane images of 3D scatter-rejection optical projection tomography reconstructions. Two, 200x200x200 μm³ inclusions with a 10% increase in attenuation coefficient were simulated, separated by 200 μm. Projections were simulated at 10 μm spacing in a d x d grid (d = thickness of object) for 72 evenly spaced projections about each object. The detected energy of light at each detector location was assumed to be 600 fJ. Poisson noise was estimated and results for one noise-realization using standard filtered backprojection reconstruction are shown in under all conditions. The top row of images pertains to results from a 3 mm thick average soft tissue (μ_a = 0.2 cm⁻¹, μ_s = 100 cm⁻¹, g = 0.9) object. The second row pertains to a 6 mm thick lymph node like tissue (μ_a = 0.3 cm⁻¹, μ_s = 43 cm⁻¹, g = 0.92; at 780 nm). The first column of images is for the case with no scatter-restriction, the second column represents images for an angular restriction of NA = 0.059, the third column represents images for an angular-domain restriction of NA = 0.005, the fourth column represents images for a time-domain restriction of 1 ps, and the fifth column represents images for a time-domain restriction of 0.1 ps. Units are in cm⁻¹ (attenuation).

Download Full Size | PDF

3.2.3 Small region contrast

The detected-photons’ transported density maps of various optical projection setups (e.g., Fig. 1) in average and lymph node-like biological tissues were further analyzed in terms of a contrast metric, as defined in Eq. (2). Figure 4 presents the comparison of normalized contrast for a range of time- and angular-domain scatter rejection methods in terms of detecting a 200-μm inclusion in average tissue thicknesses of 1-3 mm diameter and in lymph node tissues of 1-6 mm diameter. These results further supported the conclusions from Fig. 1: that angular-domain restriction is only effective at improving spatial resolution for thinner and/or lower scattering tissues. Improvements diminish as the thickness and/or the level of scattering of the medium increases, until all possible photon paths are equally probable for any acceptance angle (i.e., the diffusion regime, which is typically assumed to be around 1-2 cm for typical biological tissue [29]). More specifically, significant angular restriction (NA = 0.005) only provided improvements in contrast compared to more standard imaging (e.g., NA = 0.124) up to tissue thicknesses of about 2 mm for average tissue (Fig. 4(a)); whereas in lower scattering lymph nodes, it was predicted that improvements in contrast with angular restriction could be observed at least up to 6 mm (Fig. 4(b)). This is significant, as the vast majority of excised lymph nodes have elliptical shape with short-axis diameters of up to 6 mm (results not shown), thus demonstrating the potential to use relatively inexpensive angular domain imaging methods to significantly improve spatial resolution in whole lymph node imaging. Moreover, the ability to detect 200-μm diameter objects is significant as this is the size that is recognized to correlate to the smallest clinically relevant metastases in lymph node biopsy pathology (e.g [30,31].).

 

Fig. 4 Contrast values obtained for different sample types with different restrictions. The values are plotted for average tissue optical properties (μ_a = 0.2 cm⁻¹, μ_s = 100 cm⁻¹, g = 0.9; column 1), and for lymph node-like optical properties (μ_a = 0.3 cm⁻¹, μ_s = 43 cm⁻¹, g = 0.92; column 2). All simulations were based on a 100-μm FWHM Gaussian light source incident on the object at its waist and a 100-μm diameter detector, for 2D plane images of 3D scatter-rejection optical projection tomography reconstructions. Two, 200x200x200 μm³ inclusions with a 10% increase in attenuation coefficient were simulated, separated by 200 μm. Projections were simulated at 10 μm spacing in a d x d grid (d = thickness of object) for 72 evenly spaced projections about each object. Variations in time- and angular-domain restriction are presented in rows 1 and 2, respectively.

Download Full Size | PDF

The contrasts of various time-domain restrictions (0.1, 0.3, and 1 ps duration) are presented for average tissue and lymph node-like tissue in Figs. 4(b) and 4(d), respectively. A 1 ps time-resolution is about the limit for time-correlated single photon counting technology at present [32]. This temporal resolution was shown to be unable to improve resolution over standard optical detection in average tissue until thickness was greater than 1 mm, and greater than 4 mm in lymph node-like tissue. This makes sense conceptually because as the object size increases, the range of pathlengths of photons collected within a set arrival-time interval do not change, thus amplifying the improvement of early photon. It is this principle that makes it possible to improve imaging resolution with time-domain rejection imaging even for whole mouse imaging (~30 mm diameter) [7]. Another limitation of this form of time-domain restriction is that detection rates are restricted by the dynamic range of detectors and most of the dynamic range is used up by later arriving photons; however, gated detectors [33] and dead-time enhanced early photon detection [34] can help overcome this limitation. More restrictive early photon gating is possible with optical gating approaches (width proportional to the duration of each light pulse emitted by the light source: ~0.1 ps for a “femtosecond” laser) [27]. With 0.1 ps time-domain restriction, substantial improvements in contrast were observed at nearly all thicknesses and tissues investigated, with improvements increasing with tissue thickness. These results suggest that optical gating with femtosecond lasers has the potential to significantly improve resolution in mesoscopic imaging; however, it should be noted that, experimentally, optimal selection of the location of the 0.1 ps gate can be difficult, particularly considering the variability in specimen thickness across a field of view. Additionally, it should be noted that to avoid tissue damage (requiring some restriction on the laser power used) while maintaining signal-to-noise in terms of the photon budget, 0.1 ps time gates are likely not possible for tissue imaging past 3 mm diameter in typical tissue or 6 mm diameter in lymph node like tissue, without significantly increasing imaging time.

3.2.4 Mean squared error in 3D reconstructions

Mean squared error (MSE) represents another means of analyzing the accuracy of an OPT imaging and reconstruction protocol. Figure 5 presents MSE for noise-free OPT reconstructions for a range of angular- and time-domain restrictions (NA = 0.124-0.005 for angular domain, and all photons to 0.1 ps early photons for time-domain). These noise-free results present what is ultimately achievable given high enough signal. Similar to the contrast results (Fig. 4), time-domain restrictions provided improvements in reconstructions commensurate with the level of restriction for all tissue thickness and tissue types (average soft tissue and lymph node like). Additionally, angular-domain restriction provided no appreciable improvements in MSE for any thickness tested (1-3 mm) of average soft tissue like objects; however, the most restrictive angular-domain restriction (NA = 0.005) provided significant improvements in MSE compared to all other angular-restrictions, particularly at the 6 mm thick lymph node like object.

 

Fig. 5 Mean squared error values obtained for different sample types with different restrictions. The values are plotted for average tissue optical properties (μ_a = 0.2 cm⁻¹, μ_s = 100 cm⁻¹, g = 0.9; column 1), and for lymph node-like optical properties (μ_a = 0.3 cm⁻¹, μ_s = 43 cm⁻¹, g = 0.92; column 2). All simulations were based on a 100-μm FWHM Gaussian light source incident on the object at its waist and a 100-μm diameter detector, for 2D plane images of 3D scatter-rejection optical projection tomography reconstructions. Two, 200x200x200 μm³ inclusions with a 10% increase in attenuation coefficient were simulated, separated by 200 μm. Projections were simulated at 10 μm spacing in a d x d grid (d = thickness of object) for 72 evenly spaced projections about each object. Variations in time- and angular-domain restriction are presented in rows 1 and 2, respectively.

Download Full Size | PDF

4. Conclusion

The culmination of the simulation and lymph node optical property experimental results presented in this work provides guidance for the development of OPT systems for tissue biopsy and in particular lymph nodes. Specifically, it was demonstrated that low cost angular domain OPT could be used to achieve detection of the smallest clinically relevant micrometastases. Ultimately, this can help develop lymph node imaging systems able to outperform standard pathology protocols in terms of processing time and sensitivity of detecting cancer spread. Future steps will include the development of lymph node staining 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 model-based 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.