Origins of subdiffractional contrast in optical coherence tomography

: We demonstrate that OCT images quantify subdiffractional tissue structure. Optical coherence tomography (OCT) measures stratified tissue morphology with spatial resolution limited by the temporal coherence length. Spectroscopic OCT processing, on the other hand, has enabled nanoscale sensitive analysis, presenting an unexplored question: how does subdiffractional information get folded into the OCT image and how does one best analyze to allow for unambiguous quantification of ultrastructure? We first develop an FDTD simulation to model spectral domain OCT with nanometer resolution. Using this, we validate an analytical relationship between the sample statistics through the power spectral density (PSD) of refractive index fluctuations and three measurable quantities (image mean, image variance, and spectral slope), and have found that each probes different aspects of the PSD (amplitude, integral and slope, respectively). Finally, we found that only the spectral slope, quantifying mass scaling, is monotonic with the sample autocorrelation shape.


Introduction
Spectroscopic interrogation of biological tissue can reveal valuable information on the structural organization due to scattering [1][2][3][4] as well as the molecular composition from absorption [5][6][7][8]. Microarchitectural alterations can help quantify cellular dynamics and understand the progression of pathology. For example, many of the early changes associated with neoplasia, eg. alterations in cytoskeleton structure [9], remodeling of higher order chromatin [10] and stromal stiffness [11], occur in the nanometer length scale, and thus appear histologically normal.
Quantifying ultrastructure (eg. morphology below the imaging resolution of standard light microscopes) in-vivo, can provide vital insight for scientific discovery and potential targets for therapeutic intervention [12]. Conventional widefield imaging techniques are limited by the physical constraints of the Abbe diffraction limit (∼200 nm for visible light systems). While super resolution optical techniques like Photon Localization Microscopy (PLM) [13][14][15] and Structured Illumination Microscopy (SIM) [16] have allowed visualization of biological structures as small as ∼10 nm, these techniques are most often unavailable for in-vivo imaging. Sample preparation can include tissue excision, fixation, and staining, and these processes combined with high irradiation necessary for imaging may alter or even destroy the native tissue structures of interest [17,18]. Furthermore, while there are a variety of nanoscale imaging modalities, few can evade the trade-off between resolution and penetration depth.
Optical Coherence Tomography (OCT) is an interferometric imaging technique to reconstruct tissue morphology in 3D with up to several millimeters penetration depth and down to a micron resolution [19]. Inverse Spectroscopic Optical Coherence Tomography (ISOCT) has been developed to characterize the refractive index autocorrelation function (ACF) through spectral domain analysis [20]. Under the Born approximation for which weakly scattering tissue applies, the optical scattering spectra can suitably be described by a single power law decay along wavelength [3,21,22]. Specifically, by measuring the power law decay of the backscattered spectrum, µ b (k), one can quantify the distribution of mass density in tissue [22,23]. Tissues with a high density of small scattering particles with respect to the illumination wavelength such as mitochondria or lysosomes will exhibit Rayleigh scattering spectral decay (eg. µ b ∝ k 4 ). While tissues with larger structures, resulting in more forward scattering behavior, leads to a decrease in the power law measurement. Furthermore, this metric has been experimentally shown to be sensitive to tissue structures at length scales from 35-450 nm [20], well below the conventional optical imaging diffraction limit. This development, coined as inverse spectroscopic optical coherence tomography, or ISOCT, enables a standard OCT device to yield high tissue penetration with nanosensing capabilities by quantifying the backscattered spectral response.
However, this has surfaced a point of confusion: if the backscattered spectra is sensitive to changes in ultrastructure from length scales below the diffraction limit, how is this information obscured in the OCT image? Thus, it is not entirely understood if and how subdiffractional information is folded into the image space and whether nanoscale sensitive descriptive metrics can be retrieved by image-based (as opposed to spectroscopic) analyses. It is commonly assumed that OCT images were strictly diffraction (lateral) and bandwidth (axial) limited, and any structural information probed by frequencies below the minimum illuminated wavelength were not resolvable. Additionally, several analyses have been developed to utilize the spectral features (eg. spectral center of mass [24], or spectral maximum [25]) of spectroscopic OCT in order to differentiate tissue types, ultrastructure or characterize pathology without a unified theory as to how. Thus, there is a fundamental gap in the connection between sample statistics (e.g., its power spectral density), and the resulting OCT image.
The goal of this work is to establish a theoretical framework to understand and quantify how ultrastructure is embedded in the spectral response and incorporated into the OCT image. First, we establish an FDTD based spectral OCT simulation to model the OCT image and spectral response from a given sample. We validate this simulation by cross referencing the simulated image and spectrum with that of experimentally measured beads. Then, using random media samples with known statistical distributions, we show through a series of simulations that (1) the backscattered spectral decay is indeed sensitive to structural changes to mass-density distribution in the range of 30-300nm, despite the minimum resolvability of the image being limited by the temporal coherence length (∼1µm). However, (2) the image mean and variance probe the sample's power spectral density amplitude and integral, respectively, which have non-monotonic relationships with the mass-density distribution.

Experimental setup
The experimental setup can be simplified to a basic Michelson Interferometer with an open-air source spectrum spanning 500-700nm, as shown in Fig. 1. An in-depth description of the system is described elsewhere [26]. Each scanned volume spanned a lateral extent of 1.7mm x 1.7mm with 512 x 512 pixels. The system was calibrated with a 1% by volume solution of 80nm polystyrene beads mixed in deionized water to ensure that the effects of spectral dependent system roll-off were minimized.
In order to generate the OCT image, each interferogram signal was first normalized by the reference spectra and its DC component was subtracted. Computing the Fourier transform (FT) converts the wavenumber k-domain signal to the spatial z-domain, or image space. The wavenumber is related to the wavelength as k = 2π/λ. Spectral analysis requires applying a shorttime Fourier transform (STFT) on the normalized interferogram. To do this, after normalization, the input signal I N is convolved with a series of shifted Gaussian windows in the k-domain, G(k, k 0 ) centered at k 0 , which converts the k-space signal into a k 0 − z dependent extraction. With this, one can isolate specific regions in depth for which to compute the backscattered spectra. Furthermore, integrating over z, and squaring the result yields the spectral backscattering Fig. 1. Schematic of benchtop spectral domain OCT system utilizing a visible light source (VIS) coupled to single mode (SM) optical fibers, passing to a beam splitter (BS) and galvanometric scanning mirrors (GV). In the reference arm, the beam is modulated by a neutral density filter (ND) and dispersion compensation quartz plates (DC). Inset shows how the media is integrated into FDTD: a perfectly matched layer (PML) pads the media border, and the far field scattered E-field is collected within the system specific numerical aperture. coefficient µ b (k), as follows: Here, n 0 is the mean refractive index in the sample. The FT exponential argument leading factor of 2 accounts for the forward and reverse paths of light. In experiments, the STFT Gaussian filters have a full width at half max (FWHM) of 0.12µm −1 , resulting in a spatial filter with mean FWHM of 8nm. This process has been outlined in more depth in a previous publication [20].

Simulation
A simulation method to calculate the propagation of the optical electromagnetic field through an inhomogeneous dielectric media with nanometer-level resolution was required in order to model the nanoscale sensitivity of ISOCT. With only the illumination beam and complex refractive index (RI) of the media specified within a three-dimensional cartesian grid, the finite-difference time-domain (FDTD) computational solution of Maxwell's equations [27] was employed to obtain the wavelength-dependent electric and magnetic fields throughout the sample and its near vicinity. To ensure numerical convergence, the FDTD voxel size was progressively reduced below λ/20, ultimately to as fine as 5nm. The FDTD software employed, Angora, was developed in the Backman lab over a period of 13 years [28,29]. Angora is an open-source software that embeds FDTD and subsequent rigorous near-to-far field transformations as the physics kernel within a highly configurable optical system that functions as a "microscope in a computer". Angora has been validated for accurately computing the penetrating and scattered optical electromagnetic fields and synthesizing the microscope images of structures comprised of both nondispersive [29] and dispersive media [30].
Thus, integrating the output of Angora into an ISOCT based simulation, we establish our modeled system. This included illuminating the sample with an appropriate source, computing the scattered electromagnetic near-field adjacent to the sample, collecting the vectoral field within the system-specific exit pupil, refocusing the field through the objective lens, constructively adding the reference, and finally analyzing the signal as consistent with OCT (image) and ISOCT (spectral) based analysis. To allow for flexibility to account for varying system configurations, the simulation can be decomposed into these primary components: illumination, scattering, collection, refocusing, and processing.

Illumination
Although Angora is capable of modeling a variety of spatially coherent incident beams using plane-wave decomposition, eg. Gaussian, Bessel, etc [28], the simplest and most relevant for OCT imaging is the linearly polarized, spectrally Gaussian, plane-wave propagating in the +z direction: Here, E 0 is the real amplitude of the incident electric field (E-field), ω is the angular frequency (radians per second), k is the wavenumber, k 0 is the source spectra's center wavenumber, and σ k is its standard deviation. All simulations were illuminated with a stationary, broad-band source (500-700nm), with a temporally Gaussian profile: center wavelength of 583nm (k 0 = 10.8µm −1 ) and standard deviation of 48nm (σ k = 0.6µm −1 ).

Scattering
The incident beam is introduced into the simulation space, as visualized in the inset of Fig. 1. The media, generated in three dimensions, has an upper layer of air, and is padded on its remaining boundaries by a perfectly matched layer, eg. n 0 = n mean for a continuous media sample, or n 0 = n background for a bead sample. This border is necessary to emulate a sample with infinite extent by absorbing the propagating wave at the boundaries and reducing spurious boundary reflections. The FDTD output is a vector field consisting of the complex, far-field scattering amplitude and direction for every sampled wavelength and observation position: Here, s x and s y are the directional cosine components of the vector fields as projected onto the x and y axis of the entrance pupil coordinate plane.

Collection and refocusing
Modeling the objective lens, the scattered vector field is masked by the collecting numerical aperture (NA) and focused onto the imaging plane by applying the Debye-Wolf Integral to the scattered field, E s [31].
Here, Ω img = ds x ds y /s z is the solid angle bounded by collection NA. The mask is applied by limiting the collection angles to those within the region defined by √︂ s 2 x + s 2 y ≤ NA. The collection NA was chosen to be 0.1 for continuous media and 0.05 for bead simulations to best imitate the physical OCT system.
Next, the reference field was generated to replicate a field incident to a perfectly reflective dielectric media (mirror). This reference field was chosen to be a plane wave with a phase offset of 2z ref . The variable z ref defines the path difference between the reference and scattered fields, and the multiple of 2 accounts for the forward and reverse paths. Its value was chosen to avoid the complex conjugate artifact from occluding the sample image [32].
Here, A is the amplification of the reference with respect to the incident source. In practice, this is dictated by a neutral density filter placed in the reference arm to ensure the reference is sufficiently greater than the sample arm, and the cross-interference term dominates the signal. In these simulations, A was chosen to be 100. Finally, the focused reference and scattered fields are summed together and squared to produce the OCT interferogram.
Here ρ is the detector's responsivity (units Amperes/Watt), and the angular brackets represent the detector's inherent temporal integration [32]. For simplicity, we set ρ = 2, although to obtain absolute measurements in comparison with an OCT system, this value must be explicitly measured.

Processing
The simulated interferogram was processed as done with experimental data. First, each interferogram was normalized by the reference intensity and its DC component was subtracted.
To obtain volumetric structural images, a full FT along the axial dimension converted each interferogram into z-space. Additionally, as described above, after normalization, STFT processing generated the simulated spectra. This processed was repeated for every x-y position to compute the 3D image and averaged together to accumulate the backscatttered spectra.

Experimental validation of simulation
With the volumetric refractive index of the media defined as an input, the FDTD-ISOCT simulation outputs a spatially dependent interferogram: I(x, y, k). After processing, the backscattered spectrum and OCT images are computed. Naturally, validation of this simulation is required, and polystyrene microspheres make for predictable phantoms. Experimental data was collected to compare the images and spectra of individual beads. A solution of 5.2µm polystyrene beads with a coefficient of variance of 4%, was diluted with deionized water at a concentration of 0.1% to allow individual beads to be isolated using image analysis. An area of 1.7mm x 1.7mm was scanned with 512 x 512 pixels in the lateral dimensions. From the OCT volume, individual beads were automatically detected, and centered within an encompassing window of dimensions 10x10x30µm. In computing the spectra, regions of the OCT volume without beads were masked out before computing the STFT. After processing, the median of 500 individual beads was computed to generate the image (Figs. 2(b) and 2(d)) and spectra (Fig. 2(e)) in order to be robust to outliers, such as coalesced or damaged microspheres.
Then, a single 5.2µm polystyrene bead was simulated with a refractive index of 1.59 in a background of water, n 0 = 1.33. In the visible light range simulated, the dispersive effects of water were ignored. The media was 400 voxels along each dimension, with a voxel length of 20nm, which was sufficiently small to ensure numerical dispersion was insignificant [27]. The image was generated in 3D, and lateral and axial cross sections are shown in Figs. 2(a) and 2(c), while the spectra is shown in Fig. 2(e). It is clear that the simulated image and spectra agree with experimental results within a reasonable margin of variability.

Modeling biological tissue as a continuous random media
Towards the purpose of light scattering simulations, biological tissue, albeit highly complex in structure, can been approximated via parametric models by assuming various degrees of simplification. Some approaches consider discrete spherical scatterers with a fractal size distribution to emulate the fractality of tissue organization [33,34]. Using electron microscopy, however, it has been shown that the mass-density distribution of biological samples can be more accurately described as a continuous random media, rather than a collection of discrete particles [2,10,35]. Rogers et al. thoroughly discussed the use of one such model, the Whitle-Matérn (WM) model and its utility in accurately describing tissue structure [2]. Here, we use the flexible, 3-parameter WM model to parameterize the statistically isotropic spatial auto-correlation function (ACF) of the media's refractive index fluctuations: Here, r is the separation distance, K ν {.} is the modified Bessel function of the second kind with order ν, A n is the fluctuation strength of RI variations, L n is the characteristic length of heterogeneity, and D is the mass scaling parameter, defining the functional form of the distribution [36,37]. When 2<D<3, the tissue is organized as a mass fractal, for which D is related to the fractal dimension. Although there are several interpretations for the power law scaling parameter as it relates to the fractal dimension [38], it must first be derived from the scaling property of fractals, for which its mass within a radius r scales with fractal dimension D f , as M(r) ∼ r D f [39]. This is analogous to a pure power-law ACF with scaling related to the fractal dimension: B n (r) ∝ r D f −3 [23,40]. Although this formulation is not identical to the WM form, within the regime of D between 2 and 3, D is a reasonable approximation for D f . For 3<D<4, however, the correlation function takes on the form of a stretched exponential, D = 4 yields an exponential decay and as D approaches infinity, the ACF approaches that of a Gaussian. Also, for a fixed L n , larger values of D imply larger length scales of tissue heterogeneity [20].
A n was set according to Eq. (7), in order to normalize the standard deviation of RI fluctuations to a constant value σ n at the smallest resolvable spatial separation r min [2]: Based on the Wiener-Khinchin theorem, the sample's power spectral density (PSD) can then be formulated from Eq. (6) as follows [2]: Here, k s is the sampled wavenumber in the media (eg. for a purely backscattered response k s = 2k), and Γ is the gamma function. Furthermore, if k s L n >>1 and the Born approximation applies, as in the case of primarily forward or weakly scattering media, the backscattering coefficient µ b can be written with respect to these physical parameters [20]: Thus, measuring µ b (k), one can readily extract the mass scaling parameter D using a log transform and linear regression.

Spectral length scale sensitivity
To explore the sensitivity of mass scaling D to particular length scales, we generate a continuous random media and remove the contribution from specific length scales by applying Gaussian filters with varying filter sizes. The RI distribution of the random media sample was computed with the ACF defined in Eq. (6). The model parameters for this simulation were D = 2.7, L n = 1µm, and r min = dx/2, which normalizes the standard deviation of fluctuations of the numerical sample to σ n = 0.03. These parameters have been shown previously to be consistent with various types of tissue [22,41]. The media was generated with 650 pixels 3 along each dimension and voxel length dx = 20nm, thus spanning 13µm along each dimension. Next, the media was convolved with Gaussian low pass filter (denoted as W LP ) with FWHM ranging between 20 − 100nm and with high-pass filters (W HP ), ranging from 100 − 500nm. Examples of filtered media are shown in Figs The OCT images and spectra were simulated and processed as described previously. Example B-scans are shown in Figs. 3(d)-f. We note here that although the standard deviation of RI fluctuations is consistent for all samples, reflections from top and bottom interfaces are strongest for low pass samples, as visible in subplot (e).
Simulated spectrum are shown in Fig. 4. The spectral shape parameter D was calculated for each simulation and is shown next to each spectrum. Percent deviation from simulated D values are shown in 4(c), with filter sizes that deviate from 5% indicated in yellow (∼30-300nm).
Additionally, the ACFs of the simulated images along the axial dimension were computed, normalized by the unfiltered variance and are plotted in Figs. 5(c) and Figs. 6(c). It is clear from these figures, that although the backscattering spectra are highly sensitive to mass-density perturbations in these subdiffractional length scales, the simulates ACFs are not. With the exception of the variance of the simulated images, eg. B ′ n (0) as shown in the inset, the structure of OCT images do not change appreciably as filter size changes, as evident by the similarity between ACF shapes.   4. Backscattered spectra is shown for low pass (a) and high pass (b) filtered media after processing the spectral OCT simulation. Percent change in D for all filter sizes are shown in (c). Here, 5% change is highlighted with grey, and the filter size regime yielding changes greater than 5% is indicated with yellow.   6. Sample media ACF (a) and PSD (b) for varying high pass filter sizes, with FWHM displayed in legend. The unfiltered media is denoted with filter size W HP ∞nm . After OCT simulation, the ACFs for corresponding simulated images are taken along the axial dimension, normalized by the r = 0 value and are plotted together in (c). The simulated image ACF normalization is again shown in the inset of (c).

Spatial length scale sensitivity
After establishing the spectral response to perturbations of subdiffractional length scales, we attempt to quantify the sensitivity of image space analysis. While there are a plethora of metrics that can be extracted from image space processing, the scope of this paper specifically analyzed the first and second moment of the image intensity, e.g. the mean and variance. The following derivations sought to relate the expected OCT image mean and variance to the sample statistics through its PSD.
Incorporating the assumptions of a weakly scattering, non-absorbing media and a linearly polarized incident beam [41,42], we can formulate the normalized backscattering spectral intensity I N (x, y, k) using the basic interference equation: I N is computed by normalizing I, the measured interferogram, by the reference intensity I R . The variables t 01 and t 10 are the Fresnel transmission coefficients of the forward and reverse paths through the air-sample interface, ρ is the detector gain, and k is the wavenumber within the sample. Whereby, following the notation in [41], n 1D is equivalent to the media's RI fluctuations, n ∆ = n(x,y,z)−n 0 n 0 , smoothed by the windowing functions of the numerical aperture and the source spectra (T kNA and T 2k respectively): Thus, the sample PSD can be written as: Since the sample is only considered in the region of positive depth, z>0, the property of Fourier transforms for real, causal functions will enable further simplification: Substituting the resulting PSD of the media's fluctuations Eq. (12) and Eq. (13) into Eq. (10), and introducing a Gaussian STFT with windowing function G(k 0 ), we obtain the simplified expression for the statistically averaged normalized backscattered spectrum: Relating the DC in the spatial domain to the minimum sampled frequency k min in Fourier domain, and considering the STFT is not applied to image space analysis (eg. in the limit as the G(k 0 ) filter width approaches ∞ ), the mean image intensity can be simplified as: We introduce L to be the thickness of the sample. From Parsevel's theorem, the expected image variance can be represented in terms of the sample PSD as follows: Here, k 0 is the center illumination wavenumber in the sample. Finally, for completeness, we include here a relationship between the statistically averaged backscattered spectrum, shown in Eq. (14), and the backscattered coefficient, after performing a short-time Fourier transform: To verify the theoretical formulation for these image space metrics, we simulate a series of continuous random media with varying size distribution parameter D, ranging from 2.1 − 3.9. Each volume measured 10µm along each dimension, maintained a constant standard deviation (σ n = 0.03), mean index of refraction (n 0 = 1.38), several values for L n (0.5, 1.0, and 1.5µm), and was repeated 5 times. The OCT images were simulated and the image mean and variance were calculated within the bounding box of the sample and are plotted in Fig. 7(a) and 7(b) with the analytical formulations computed from Eq. (15) and Eq. (17) overlaid. Additionally, the mass scaling D was measured after simulation using spectral fitting, and is plotted in Fig. 7(c). As shown in Eq. (17), the log-slope is equivalent to 4 − D, from which D can be measured, and is plotted along the y-axis.

Conclusion
In this paper, we have demonstrated the use of FDTD in the simulation of spectroscopic OCT for the analysis of nanostructural characterization. We established a theoretical framework to relate the sample statistics, through its power spectral density, to the backscattered spectrum and the OCT image. This framework provides valuable insight into how ultrastructural statistical properties of a continuous random media are encoded in the diffraction limited OCT image, specifically the mean and variance. It was found that while the spectral slope is highly sensitive to length scale alterations below the diffraction limit, in the image space, these alterations manifested in changes to the image mean and variance, albeit non-monotonically.
Since the standard deviation of RI fluctuations were normalized across all samples, the contrast measured between samples is solely due to the differences in mass-density size distribution. We observe from Eqs. (15 -(17)), that each metric (image mean, image variance, and D), probe different aspects of power spectral density of the media's refractive index fluctuations (amplitude, integral and log-slope, respectively). Furthermore, while we utilized the Whittle-Matén model to describe the sample's ACF, the relationship between the PSD and the image and spectral based metrics applies to any PSD under the Born approximation, for example a purely power law ACF: B n (r) ∝ r D−3 .
For the WM model specifically, while the spectral slope changes linearly with mass-density scaling parameter D (Fig. 7(c)), the relationship between image mean ( Fig. 7(a)) and variance ( Fig. 7(b)) invert when D ≈ 3, and that non-monotonicity is dependent on the characteristic length L n . Furthermore, while we observe that image mean and variance have smaller inter-simulation variation, the curves depend highly on L n . We also see this inversion when applying the low-pass filters to simulated media, as shown in the inset of Fig. 5(c). There is a particular point at which D and L n generate a sample for which the characteristic scattering size of the media is nearing the wavelength of light illuminated. When L n is reduced to sub-wavelength values, the media appears as a collection of Rayleigh particles, and larger D is required to find that maximal point. As L n increases beyond the size of the wavelength, the media appears more homogeneous from the perspective of the illuminated wavelength, and scattering is reduced. For a fixed L n , however, changing D also changes the average mass density size.
The spectral slope (Fig. 7(c)), on the other hand is independent of L n , and is linear and monotonic with size distribution D. Additionally, from this analysis, we see that within specific regimes (eg. D< ∼ 3 or D> ∼ 3), mass scaling D is monotonic and correlated with image intensity. For example, for D<3, we see a positive correlation between image intensity and D; this is consistent with previous OCT measurements [43]. In past experiments, however, image mean and variance were not found to be reliable metrics to differentiate between tissue types due to inconsistent variations within a tissue layer [43]. Additionally, image mean is not a practical biomarker, as it cannot be compared to samples from other systems without a reliable phantom with known reflectance properties.
Future work will seek to improve simulation limitations for a more experimentally accurate, and tissue specific analysis. For example, optical aberrations, noise or tissue absorption could be considered, which all become more problematic in deeper layers of tissue. Additionally, it would be useful to examine the effect of sample statistics for multilayered tissue structures. Despite its limitations, we believe the incorporation of FDTD-ISOCT simulation for nanoscale resolved simulations will be most invaluable to aid in characterization, analysis and optimization for biological research.