## 1.

## Introduction

It is well established that reflectance measurements acquired from within the subdiffusion regime (i.e., source-detector separations smaller than a transport mean free path ${l}_{\mathrm{s}}^{*}$) preserve information about the shape of the scattering phase function $P(\theta )$.^{1}2.3.4.5.^{–}^{6} Along with this increased sensitivity comes the need for improved models of light scattering that are both versatile and contribute improved insights into tissue characterization. However, despite the direct link between scattering and fundamental tissue ultrastructure, many models of light scattering focus on the role of wavelength-dependent empirical parameters, which determine the shape of $P(\theta )$, rather than the more insightful physical properties. Given the wealth of information that may potentially be extracted from within the subdiffusion regime, we argue that the primary mandate of any extended light-scattering model should not be to simply expand the empirical parameter space for the sake of obtaining more optical properties or improved fits but also to forge a more fundamental understanding of the tissue structure, which leads to the observed reflectance signals.

Contrasting with the common treatment of scattering is the admirable way in which quantification of absorption is frequently applied. In absorption analysis, it is well recognized that providing values of the spectrally resolved absorption coefficient ${\mu}_{\mathrm{a}}(\lambda )$ are often less informative than providing the type (e.g., hemoglobin, melanin, fat, etc.), concentration, and organization of absorbing species, which lead to that signal. With regards to hemoglobin absorption fitting, measurement of ${\mu}_{\mathrm{a}}(\lambda )$ provides a vehicle from which to extract more physical information about microvasculature, such as blood volume fraction $v$, hemoglobin oxygen saturation ${S}_{\mathrm{o}}$, and blood vessel radius ${R}_{\text{vessel}}$.^{7}8.9.^{–}^{10} Should the specific value of ${\mu}_{\mathrm{a}}$ at an arbitrary wavelength become required for some computation, the underlying physical properties can be quickly expanded to yield an answer.^{11}

Taking insight from the way in which absorption analysis is applied, we and others have proposed a scattering model, which takes its main inspiration from considering the physical structure of tissue.^{12}13.^{–}^{14} Assuming the validity of the Born approximation (a safe assumption for most solid tissues^{15}), the fundamental physical tissue characteristic, which gives rise to light scattering, is the spatial refractive index (RI) autocorrelation function, ${B}_{n}({r}_{d})$.^{16}17.^{–}^{18} If ${B}_{n}({r}_{d})$ is known, a series of simple mathematical transformations can be used to first derive $P(\theta )$ and subsequently optical properties, such as the scattering coefficient ${\mu}_{\mathrm{s}}$ and anisotropy factor $g$ as well as higher order shape parameters, such as $\gamma $, $D$, etc.^{10}^{,}^{19} Since ${B}_{n}({r}_{d})$ can similarly be defined for discrete particles (e.g., spheres,^{12} red bloods cells,^{20} rods,^{18} etc.) as well as continuous random media, it can serve as a universal mediator between all models of light scattering. The key then becomes to find a model of ${B}_{n}({r}_{d})$ that is not only versatile and mathematically convenient but also accurately describes most tissues.

In order to satisfy both criteria for achieving an appropriate model of ${B}_{n}({r}_{d})$, we employ the three-parameter Whittle–Matérn model.^{12}^{,}^{21} From a mathematical standpoint, this model encompasses a wide range of plausible functional forms for the shape of ${B}_{n}({r}_{d})$, including Gaussian, stretch-exponential, and power-law distributions. As an added benefit for the optics community, one special case of the Whittle–Matérn model includes the commonly used modified Henyey–Greenstein model. Furthermore, a number of recent publications demonstrate the suitability of using the Whittle–Matérn model for characterizing tissue mucosa.^{4}^{,}^{22}23.^{–}^{24}

In this paper, we present a unified model to quantify six physical tissue properties (three ultrastructural and three microvascular) from a single spectral measurement of the spatial reflectance profile $p(r,\lambda )$ at subdiffusion lengthscales. Toward this end, there are four primary goals of the current work: (1) to present a unified model (using two previously developed models) that relates light propagation in biological media to the underlying tissue ultrastructure and microvasculature; (2) to provide a new empirical Monte Carlo-based model that enables rapid generation of $p(r,\lambda )$, specifically at subdiffusion lengthscales; (3) to demonstrate the application of a newly developed inverse algorithm to extract physical properties from a measurement of $p(r,\lambda )$; and (4) to announce MATLAB codes posted online^{25} for use by other researchers, who wish to carry out the methods and analysis presented in this work.

## 2.

## Theory

In the following theory subsections, we present a model of light–tissue interaction, which incorporates scattering from a continuous random medium specified by the Whittle–Matérn model^{21}^{,}^{26}^{,}^{27} and absorption assuming packaging of hemoglobin within blood vessels.^{7}8.^{–}^{9} The full medium model is characterized by six physical sample properties: the variance of RI ${\sigma}_{n}^{2}$, the characteristic lengthscale of RI heterogeneity ${L}_{n}$, the shape parameter of RI distribution $D$, the blood volume fraction $v$, the oxygen saturation ${S}_{\mathrm{o}}$, and the blood vessel radius ${R}_{\text{vessel}}$. As we will demonstrate, these six “physical” properties directly determine the wavelength dependence of four commonly measured “optical” properties: the scattering coefficient ${\mu}_{\mathrm{s}}$, the absorption coefficient ${\mu}_{\mathrm{a}}$, the reduced scattering coefficient ${\mu}_{\mathrm{s}}^{*}$, and the anisotropy factor $g$. Table 1 previews the relationship between the physical and optical properties, which will be explored in depth in the following two subsections.

## Table 1

Relation between optical and physical parameters.

${\mu}_{\mathrm{s}}$ | Function of ${\sigma}_{n}^{2}$, ${L}_{n}$, $D$, and $\lambda $ |

$g$ | Function of ${L}_{n}$, $D$, and $\lambda $ |

${\mu}_{\mathrm{s}}^{*}$ | Function of ${\sigma}_{n}^{2}$, ${L}_{n}$, $D$, and $\lambda $ |

${\mu}_{\mathrm{a}}$ | Function of $v$, ${S}_{\mathrm{o}}$, ${R}_{\text{vessel}}$, and $\lambda $ |

## 2.1.

### Relating Scattering Properties to ${\sigma}_{n}^{2}$, ${L}_{n}$, and $D$

Among the vast number of phase function models that have been proposed, two families dominate within biomedical optics: those are based on (1) the Henyey–Greenstein phase function (HGPF) and (2) the Mie phase function (MPF).

Although originally developed for interstellar light scattering by dust clouds,^{28} the HGPF has made its way into the field of tissue optics due, in large part, to its simple one-parameter (i.e., the anisotropy factor $g$) mathematical form and its reasonable ability to approximate tissue scattering. Unfortunately, the HGPF provides an inherently unphysical form of the phase function due to its lack of a dipole/Raleigh scattering component, a fact that is confirmed by the observed mismatch between the HGPF and goniometric measurements of cells.^{12}^{,}^{29} In order to remedy this shortcoming, the modified HGPF^{30} was developed by adding a ${\mathrm{cos}}^{2}(\theta )$ dipole term and introducing a weighting parameter that scales the relative contribution of the HGPF and dipole term. Nonetheless, while it may provide more satisfying fits with measured data, the modified HGPF remains largely empirically based and therefore provides a limited connection to the underlying physical origin of the phase function.

The MPF provides a more physical understanding of the origin of tissue scattering with a rigorous solution for scattering from spheres of all sizes. Under this model, tissue scattering is represented as an incoherent superposition of scattering from spheres of all sizes. In principle, the number of free parameters (e.g., sphere number, size, and RI) is unlimited. However, for practical reasons, a fractal (i.e., power-law) distribution of sphere sizes as well as internal and external RI values are assumed. Despite the attractiveness of a MPF, we argue that biological specimens are very infrequently organized in nice spherical structures throughout the entire range of structural lengthscales that determine single light scattering (i.e., tens of nanometers to several microns).^{31} Thus, the use of the MPF may lead to an overly simplistic and perhaps misleading conceptual understanding of tissue structure.

Given the limitations of the HGPF and MPF, let us therefore consider the physical composition of tissue and use that as a starting point from which to understand light scattering. Most biological tissues are composed of a variety of arbitrarily shaped structures ranging in size from tens of nanometers (e.g., chromatin fibers or ribosomes) to microns (e.g., collagen fibers) to tens of microns (e.g., cells). Because of this continuous distribution of sizes and widely varying shapes, most biological tissues are best modeled as a continuous random media,^{4}^{,}^{14}^{,}^{21}^{,}^{31}^{,}^{32} such as the examples shown in Fig. 1. Rogers et al.^{12} provide a more thorough argument for modeling light scattering as a continuous random medium. We also note an exception to this argument in biological samples, such as adipose tissue or cell suspensions, in which it can be well argued that the spherical geometry justifies a proper application of Mie theory.

In order to describe the distribution of RI in a continuous random medium, we employ the versatile three-parameter Whittle–Matérn model to define the RI correlation function ${B}_{n}({r}_{d})$:^{4}^{,}^{14}^{,}^{21}^{,}^{27}^{,}^{31}^{,}^{32}

## (1)

$${B}_{n}({r}_{d})={A}_{n}\xb7{\left(\frac{{r}_{d}}{{L}_{n}}\right)}^{\frac{D-3}{2}}\xb7{K}_{\frac{D-3}{2}}\left(\frac{{r}_{d}}{{L}_{n}}\right),$$^{12}where ${r}_{\mathrm{min}}$ is the minimum structural length-scale of fractality. Here, we approximate ${r}_{\mathrm{min}}=2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$, corresponding roughly to the size of a number of fundamental biomolecules, which compose biological tissue (e.g., amino acids, monosaccharides, B-form DNA, etc.).

^{33}

^{,}

^{34}

Figure 2(a) shows examples of ${B}_{n}({r}_{d})$ for $D=2.0$, 4.0, and 6.0. Figure 1 shows one possible realization of the spatial distribution of RI for each of these three distributions. Under the first Born approximation^{16} (also known as the Rayleigh–Gans–Debye theory), the power spectral density ${\mathrm{\Phi}}_{s}({\overrightarrow{k}}_{\mathrm{s}})$ is the three-dimensional Fourier transform of ${B}_{n}({\overrightarrow{r}}_{d})$.^{17} Applying the Whittle–Matérn model, which assumes spherical symmetry, ${\mathrm{\Phi}}_{\mathrm{s}}({k}_{\mathrm{s}})$ can be found as^{31}

## (3)

$${\mathrm{\Phi}}_{\mathrm{s}}({k}_{\mathrm{s}})=\frac{{A}_{n}{L}_{n}^{3}\mathrm{\Gamma}\left(\frac{D}{2}\right)}{{\pi}^{3/2}{2}^{(5-D)/2}}\xb7{(1+{k}_{\mathrm{s}}^{2}{L}_{n}^{2})}^{-D/2},$$Defining scattering for arbitrary polarization, the amplitude scattering matrix per unit volume polarizability $\mathbf{s}(\theta )$ can be found by incorporating the dipole scattering pattern into ${\mathrm{\Phi}}_{\mathrm{s}}({k}_{\mathrm{s}})$:^{35}

## (4)

$$\mathbf{s}(\theta )=-i{(k{n}_{\mathrm{o}})}^{3}\xb7\sqrt{{\mathrm{\Phi}}_{\mathrm{s}}({k}_{\mathrm{s}})}\xb7\left[\begin{array}{cc}\mathrm{cos}\text{\hspace{0.17em}}\theta & 0\\ 0& 1\end{array}\right]=\left[\begin{array}{cc}{s}_{2}& 0\\ 0& {s}_{1}\end{array}\right]\mathrm{.}$$^{18}

## (5)

$$P(\theta ,\phi )=\frac{{|{s}_{2}(\theta )\mathrm{cos}(\phi ){E}_{\parallel}+{s}_{1}(\theta )\mathrm{sin}(\phi ){E}_{\perp}|}^{2}}{{\iint |{s}_{2}(\theta )\mathrm{cos}(\phi ){E}_{\parallel}+{s}_{1}(\theta )\mathrm{sin}(\phi ){E}_{\perp}|}^{2}\mathrm{sin}(\theta )\mathrm{d}\theta \text{\hspace{0.17em}}\mathrm{d}\phi},$$^{18}

^{,}

^{36}

## (6)

$$\sigma (\theta ,\phi )={(k{n}_{\mathrm{o}})}^{4}(1+{\mathrm{cos}}^{2}\text{\hspace{0.17em}}\theta ){\mathrm{\Phi}}_{\mathrm{s}}({k}_{\mathrm{s}}),$$## (7)

$${\mu}_{\mathrm{s}}=2\pi {\int}_{-1}^{1}\sigma (\mathrm{cos}\text{\hspace{0.17em}}\theta )d\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\theta ,$$## (8)

$$g=\frac{2\pi}{{\mu}_{\mathrm{s}}}{\int}_{-1}^{1}\mathrm{cos}\text{\hspace{0.17em}}\theta \xb7\sigma (\mathrm{cos}\text{\hspace{0.17em}}\theta )d\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\theta .$$^{30}These can be calculated as

## (9)

$${g}_{2}=\frac{\pi}{{\mu}_{\mathrm{s}}}{\int}_{-1}^{1}(3\text{\hspace{0.17em}}{\mathrm{cos}}^{2}\text{\hspace{0.17em}}\theta -1)\xb7\sigma (\mathrm{cos}\text{\hspace{0.17em}}\theta )d\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\theta .$$Analytical equations for ${\mu}_{\mathrm{s}}$ and $g$ under the Whittle–Matérn model were calculated by Rogers et al. and can be found in Ref. 21 The reduced scattering coefficient ${\mu}_{\mathrm{s}}^{*}$, which is primarily relevant for multiple scattering samples, is defined as ${\mu}_{\mathrm{s}}^{*}={\mu}_{\mathrm{s}}\xb7(1-g)$. Note that the scattering mean free path ${l}_{\mathrm{s}}$ and the transport scattering mean free path ${l}_{\mathrm{s}}^{*}$ are simply the inverse of ${\mu}_{\mathrm{s}}$ and ${\mu}_{\mathrm{s}}^{*}$, respectively.

Furthermore, it should be noted that the widely recognized power-law decay in ${\mu}_{\mathrm{s}}^{*}(\lambda )$ is not just an empirical observation but has physical origins in the shape of ${B}_{n}({r}_{d})$. In the limit of high $k{L}_{n}$ (satisfied in most tissues), ${\mu}_{\mathrm{s}}^{*}\propto {\lambda}^{D-4}$ for $D\le 4$. For $D>4$, the value of ${\mu}_{\mathrm{s}}^{*}$ will be constant across wavelength. Thus, quantifying the power of the decay in ${\mu}_{\mathrm{s}}^{*}(\lambda )$ can provide insight into the spatial organization of ultrastructures.

The validity of applying the Born approximation in continuous random media to calculate measurable optical parameters has been explored by Çapoğlu et al. using rigorous finite-difference time-domain analysis.^{15} Provided that the relationship ${\sigma}_{n}^{2}{(k{L}_{n})}^{2}\ll 1$ is satisfied, the Born approximation will yield an accurate measure of light scattering. Conceptually speaking, this criterion will be valid for any tissue in which it is possible to define a value for ${l}_{\mathrm{s}}$. In other words, the Born approximation can be accurately applied to all tissues for which the concept of optical properties, such as ${l}_{\mathrm{s}}$, ${\mu}_{\mathrm{s}}$, ${\mu}_{\mathrm{s}}^{*}$, etc., can be applied (i.e., most tissues that are studied in the biomedical optics community).

## 2.2.

### Relating Absorption Properties to $v$, ${S}_{\mathrm{o}}$, ${R}_{\text{vessel}}$

In many models of light–tissue interaction, it is assumed that the absorbing species are uniformly distributed throughout the medium. In some cases such as absorption in the skin due to melanin, this assumption can be fairly accurate. However, when looking at absorption due to the presence of hemoglobin, it becomes necessary to consider the effect of absorption localized within blood vessels (e.g., arteries, veins, etc.). For example, the blood content that supplies tissue mucosa is typically packaged in vessels ranging in size from $\sim 5$ to $50\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ in diameter (e.g., capillaries, arterioles, and venules).

In the typical case where blood content is packaged within vessels, a shading effect occurs in which the “measurable” absorption coefficient appears smaller than the “actual” absorption coefficient located within these vessels.^{7}8.^{–}^{9} This occurs as a result of the surrounding blood cells “shading” the blood cells at the center of a vessel from the outside illumination. Consequently, there is little to no contribution to the effective absorption from these inner blood cells.

In order to correct for this shading effect, the effective absorption coefficient ${\mu}_{\mathrm{a}}^{*}$ can be found as

## (11)

$${\mu}_{\mathrm{a}}^{*}(\lambda )={C}_{\text{pack}}(\lambda )\xb7v\xb7{\mu}_{\mathrm{a}}(\lambda ),$$## (12)

$${\mu}_{\mathrm{a}}(\lambda )=\frac{2.303}{64500}\xb7{C}_{\mathrm{Hb}}\xb7[{S}_{\mathrm{o}}\xb7{\u03f5}_{{\mathrm{HbO}}_{2}}(\lambda )+(1-{S}_{\mathrm{o}})\xb7{\u03f5}_{\mathrm{Hb}}(\lambda )],$$^{37}Finally, using the formalism presented by Van Veen et al., ${C}_{\text{pack}}$ can be calculated as

## (13)

$${C}_{\text{pack}}(\lambda )=\frac{1-\mathrm{exp}(-2{\mu}_{\mathrm{a}}\xb7{R}_{\text{vessel}})}{2{\mu}_{\mathrm{a}}\xb7{R}_{\text{vessel}}},$$Incorporating the effects of hemoglobin absorption into our calculations introduces three physical properties: $v$, ${S}_{\mathrm{o}}$, and the product ${C}_{\mathrm{Hb}}\xb7{R}_{\text{vessel}}$. Note that within this formalism, the effect of ${C}_{\mathrm{Hb}}$ cannot be separated from that of ${R}_{\text{vessel}}$ without additional information. As a result, a common assumption is that ${C}_{\mathrm{Hb}}=150\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{g}/\mathrm{L}$, allowing for the isolation of ${R}_{\text{vessel}}$.^{9}^{,}^{37}^{,}^{38}

## 3.

## Materials and Methods

In this section, we begin by summarizing the aspects of our Monte Carlo simulations, which are relevant to the current study.^{35} We next detail the implementation of three algorithms, which allow us to (1) smooth and compress a large library of Monte Carlo simulated $p(r,\lambda )$ results, (2) quickly and accurately generate $p(r,\lambda )$ using the prerun library of simulations, and (3) invert $p(r,\lambda )$ measurements to determine a particular specimen’s physical properties. Each of the three algorithms are written in MATLAB and are posted on our laboratory website.^{25} We conclude this section by reviewing the measurement of the spatial reflectance profile [denoted as $p({r}_{\mathrm{s}})$] at subdiffusion lengthscales using enhanced backscattering.^{5}

## 3.1.

### Monte Carlo Simulation and $p({r}_{\mathrm{s}},{z}_{\mathrm{max}})$ Library Population

Modeling of $p({r}_{\mathrm{s}})$ is achieved through solution of the radiative transfer equation (RTE).^{39} Given the complexity of solving the RTE, analytical solutions rely on numerous simplifying assumptions, which typically result in loss of accuracy at subdiffusion lengthscales.^{40} Although improved analytical solutions, such as the phase function corrected diffusion approximation by Vitkin et al.^{1} and modified spherical harmonic method by Liemert and Kienle^{41}, have demonstrated excellent accuracy at subdiffusion lengthscales, these results are still limited to a scalar approximation that is only strictly applicable for unpolarized light. For this reason, we choose to take a more brute force approach by solving for $p({r}_{\mathrm{s}})$ with electric field Monte Carlo simulation. Though less elegant than an analytical solution, Monte Carlo nonetheless yields a full solution to the RTE provided enough photon realizations are calculated. In-depth details of the Monte Carlo algorithm and software used to accurately simulate enhanced backscattering spectroscopy (EBS) are discussed in another publication.^{35} The code for this software is open source and can be downloaded on our laboratory website.^{25} Below, we summarize the relevant simulation parameters.

The simulated scattering medium is a slab of statistically homogeneous material with a total thickness of $100{l}_{\mathrm{s}}^{*}$ and RI matching at both boundaries. Each photon is injected into the slab orthogonal to the surface and allowed to scatter according to the Whittle–Matérn model discussed in Sec. 2.1. All “multiply scattered” photons (i.e., those having been scattered two or more times), which are reflected from the medium in the exact “backscattering” direction (i.e., antiparallel to the incident photon direction and with infinitely narrow solid angle), are scored in a two-dimensional (2-D) grid according to their exit distance relative to the entrance point (${r}_{\mathrm{s}}$) and their maximum penetration depth (${z}_{\mathrm{max}}$). For samples with finite thickness, integration over the ${z}_{\mathrm{max}}$ dimension from 0 to the tissue thickness provides a way to arrive at the measurable $p({r}_{\mathrm{s}})$ distribution. Four different 2-D collection grids are used to record the linear copolarization $(xx)$, linear cross-polarization $(xy)$, circular helicity preserving $(++)$, and orthogonal helicity $(+-)$ polarization channels. The radial resolution of these grids is $4\xb7{10}^{-3}{r}_{\mathrm{s}}/{l}_{\mathrm{s}}^{*}$ with total extent of $4{r}_{\mathrm{s}}/{l}_{\mathrm{s}}^{*}$, while the depth resolution is $4\xb7{10}^{-3}{z}_{\mathrm{max}}/{l}_{\mathrm{s}}^{*}$ with total extent of $2{z}_{\mathrm{max}}/{l}_{\mathrm{s}}^{*}$. All photons, which exit outside of the 2-D collection grid, are stored in the final grid element of the corresponding row or column. In order to satisfy conservation of energy, all single scattered intensity is also recorded. Each simulation is terminated when ${10}^{9}$ photon histories have been calculated.

In order to develop a library of $p({r}_{\mathrm{s}},{z}_{\mathrm{max}})$ for tissue relevant optical properties, we performed a series of simulations with varying values of $g$, $D$, and ${\mu}_{\mathrm{a}}/{\mu}_{\mathrm{s}}^{*}$. These values are summarized in Table 2. To achieve different values of $g$, we fixed $\lambda =0.633\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ and varied ${L}_{n}$. The total computation time for the entire library was over 1 million processor-hours performed on the Quest high-performance computing facility at Northwestern University.

## Table 2

Summary of p(rs,zmax) library optical properties.

$g$ | 0.7 to 0.96 in 0.02 steps |

$D$ | 2.0 to 4.0 in 0.1 steps |

${\mu}_{\mathrm{a}}/{\mu}_{\mathrm{s}}^{*}$ | 0, 0.25, 0.5, 0.75, 1, 2, and 3 |

The assumption of RI matching at the boundaries simplifies the parameter space of the model and thereby reduces the required computation time and data storage space. Although this limits the versatility of the model, we argue that this is an acceptable trade-off since it simplifies the methods needed to characterize the more informative internal ultrastructure at the expense of quantifying the boundary mismatch. Moreover, we note that in most situations, it is experimentally feasible to submerge the sample of interest in an index matching material to satisfy the boundary matching condition.

## 3.2.

### Smoothing and Compression of Monte Carlo $p({r}_{\mathrm{s}},{z}_{\mathrm{max}})$ Data

After running simulations with ${10}^{9}$ photons, the calculated $p({r}_{\mathrm{s}},{z}_{\mathrm{max}})$ distributions retained a very low level of residual noise. Still, regardless of how many photon histories are followed, some small amount of noise will necessarily remain. Therefore, further processing of the library of raw Monte Carlo simulations was performed to smooth the data to achieve essentially analytical quality curves. Moreover, this smoothing routine enabled us to compress the data by $\sim 50$ times for much easier portability. The general flow of this processing is shown in Fig. 3 with further details provided later in the paper. All processing of the data is implemented in MATLAB. The raw $p({r}_{\mathrm{s}},{z}_{\mathrm{max}})$ data shown in Fig. 4(a) is smoothed by first calculating the cumulative summation over rows (i.e., the ${r}_{\mathrm{s}}$ dimension). This process reduces the noise level for each subsequent row (enabling improved fitting) and results in a 2-D dataset with the same dimensions as the raw $p({r}_{\mathrm{s}},{z}_{\mathrm{max}})$. Next, each row is fit in a log–log space using a cubic spline with 40 breaks. We found that this number of breaks was large enough to accurately fit the underlying data but small enough to avoid fitting the noise. In order to return the smoothed data to its original form, we then calculate the difference between subsequent rows to get an intermediate smoothed version of $p({r}_{\mathrm{s}},{z}_{\mathrm{max}})$.

The same smoothing process is then repeated over the columns (i.e., the ${r}_{\mathrm{s}}$ dimension) of $p({r}_{\mathrm{s}},{z}_{\mathrm{max}})$. The two intermediate results of smoothing over rows and columns are then averaged. A final smoothing and compression step is performed by fitting each row of the 2-D data with a polynomial of order between 3 and 20. The optimal polynomial order for each row is chosen by minimizing the sum squared error (SSE) between the raw data and the fit. The polynomial coefficients for each row are then stored on the hard drive. Using this algorithm, the full library of data specified by Table 2 can be compressed from $\sim 6$ gigabytes to $\sim 70$ megabytes for each polarization channel, vastly improving data portability.

The polynomial coefficients can be quickly expanded into the final smoothed version of the $p({r}_{\mathrm{s}},{z}_{\mathrm{max}})$ array shown in Fig. 4(b) using matrix multiplication:

## (14)

$$p({r}_{\mathrm{s}},{z}_{\mathrm{max},j})=\sum _{i=0}^{20}{b}_{i}({z}_{\mathrm{max},j})\xb7{r}_{\mathrm{s}}^{i}\phantom{\rule{0ex}{0ex}}\begin{array}{c}=\left(\begin{array}{cccc}{b}_{0}({z}_{\mathrm{max},1})& {b}_{1}({z}_{\mathrm{max},1})& \cdots & {b}_{20}({z}_{\mathrm{max},1})\\ {b}_{0}({z}_{\mathrm{max},2})& {b}_{1}({z}_{\mathrm{max},2})& \cdots & {b}_{20}({z}_{\mathrm{max},2})\\ \vdots & \vdots & \ddots & \vdots \\ {b}_{0}({z}_{\mathrm{max},L})& {b}_{1}({z}_{\mathrm{max},L})& \cdots & {b}_{20}({z}_{\mathrm{max},L})\end{array}\right)\times \left(\begin{array}{c}1\\ {r}_{\mathrm{s}}^{1}\\ \vdots \\ {r}_{\mathrm{s}}^{20}\end{array}\right),\end{array}$$A qualitatively good match between the raw and smoothed $p({r}_{\mathrm{s}})$ for a typical simulation is shown in Fig. 4(c). In order to quantitatively verify the accuracy of our smoothing algorithm, we calculated the coefficient of determination (${R}^{2}$) by comparing each point in the raw data with its corresponding point in the smoothed data. We made this comparison for both the 2-D $p({r}_{\mathrm{s}},{z}_{\mathrm{max}})$ grid as well as the 1-D $p({r}_{\mathrm{s}})$ distribution obtained by integrating over all ${z}_{\mathrm{max}}$ (i.e., $p({r}_{\mathrm{s}})={\int}_{0}^{\infty}p({r}_{\mathrm{s}},{z}_{\mathrm{max}})\mathrm{d}{z}_{\mathrm{max}}$). The mean ${R}^{2}$ values are summarized in Table 3. A mean ${R}^{2}$ value of $>0.94$ for $p({r}_{\mathrm{s}},{z}_{\mathrm{max}})$ for each polarization channel indicates an accurate fit. The slight reduction in ${R}^{2}$ from the ideal value of 1 is due to the small amount of noise present in the raw Monte Carlo simulations.

## Table 3

Summary of the R2 values between raw and smoothed data.

Polarization channel | (xx) | (xy) | (++) | (+−) |
---|---|---|---|---|

Mean ${R}^{2}$ for $p({r}_{\mathrm{s}},{z}_{\mathrm{max}})$ | 0.989 | 0.944 | 0.977 | 0.987 |

Mean ${R}^{2}$ for $p({r}_{\mathrm{s}})$ | 0.999 | 0.993 | 0.996 | 0.999 |

## 3.3.

### Rapid Modeling of $p({r}_{\mathrm{s}},\lambda )$ from Polynomial Library

Having discussed how the polynomial library for $p({r}_{\mathrm{s}})$ is calculated, smoothed, and stored, we now move on to a discussion of the $p({r}_{\mathrm{s}},\lambda )$ model generating algorithm that enables modeling with nearly analytical speed and accuracy. The algorithm inputs are the desired values of ${r}_{\mathrm{s}}$, $\lambda $, sample thickness, six physical model parameters, and the polarization channel.

The process of generating $p({r}_{\mathrm{s}})$ is depicted in Fig. 5 and is as follows: first, we calculate the optical properties from the input physical properties using the equations discussed in Sects. 2.1 and 2.2. We then populate a five-dimensional grid of $p({r}_{\mathrm{s}}\xb7{\mu}_{\mathrm{s}}^{*},{z}_{\mathrm{max}}\xb7{\mu}_{\mathrm{s}}^{*},D,g,{\mu}_{\mathrm{a}}/{\mu}_{\mathrm{s}}^{*})$ using a matrix expansion of the polynomial coefficients previously stored by the smoothing and compression algorithm discussed immediately above. We populate this grid with a subset of curves containing the closest three library values (in order to perform cubic interpolation) for $D$ and $g$, as well as all library values up to and including the maximum desired ratio of ${\mu}_{\mathrm{a}}/{\mu}_{\mathrm{s}}^{*}$. Next, we loop through each of the input wavelengths to compute its specific $p({r}_{\mathrm{s}})$. In this process, we first integrate $p({r}_{\mathrm{s}}\xb7{\mu}_{\mathrm{s}}^{*},{z}_{\mathrm{max}}\xb7{\mu}_{\mathrm{s}}^{*},D,g,{\mu}_{\mathrm{a}}/{\mu}_{\mathrm{s}}^{*})$ over the ${z}_{\mathrm{max}}$ dimension up to the input sample thickness yielding a four-dimensional grid of $p({r}_{\mathrm{s}}\xb7{\mu}_{\mathrm{s}}^{*},{z}_{\mathrm{max}}\xb7{\mu}_{\mathrm{s}}^{*},D,g,{\mu}_{\mathrm{a}}/{\mu}_{\mathrm{s}}^{*})$. Finally, we do a multidimensional cubic interpolation to find $p({r}_{\mathrm{s}})$ for the values of ${\mu}_{\mathrm{s}}^{*}$, $D$, $g$, and ${\mu}_{\mathrm{a}}$ at the specified wavelength. Using a 2.5-GHz Intel Core i5-3210M processor, an entire $p({r}_{\mathrm{s}})$ spectrum for a nonabsorbing sample over the range $\lambda =0.500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ to $0.700\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ in $0.020\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ steps can be calculated in under 150 ms.

In order to validate the accuracy of the $p({r}_{\mathrm{s}})$ generating algorithm, we performed a series of 30 simulations with randomized ${\sigma}_{n}^{2}$, ${L}_{n}$, $D$, and ${\mu}_{\mathrm{a}}$ within the range of values incorporated into our library. We found an excellent fit between the simulations and output of our $p({r}_{\mathrm{s}})$ model generating algorithm, with a median ${R}^{2}$ value of 0.999 or greater for each polarization channel. Table 4 summarizes the ${R}^{2}$ values between the raw and generated data. Figure 6(a) demonstrates the excellent agreement between the simulations and output of our $p({r}_{\mathrm{s}})$ model generating algorithm even for the trial with the lowest ${R}^{2}$ value.

## Table 4

Summary of the R2 values between raw and generated data.

Polarization channel | (xx) | (xy) | (++) | (+−) |
---|---|---|---|---|

Median ${R}^{2}$ for $p({r}_{\mathrm{s}})$ | 0.999 | $>0.999$ | 0.999 | $>0.999$ |

Minimum ${R}^{2}$ for $p({r}_{\mathrm{s}})$ | 0.9865 | 0.996 | 0.9784 | 0.984 |

## 3.4.

### Inverse Algorithm to Extract ${\sigma}_{n}^{2}$, ${L}_{n}$, $D$, $v$, ${S}_{\mathrm{o}}$, and R_{vessel}

Our inverse scattering algorithm is performed by minimizing the SSE between an experimental measurement of $p({r}_{\mathrm{s}},\lambda )$ and our model for $p({r}_{\mathrm{s}},\lambda )$ to obtain ${\sigma}_{n}^{2}$, ${L}_{n}$, $D$, $v$, ${S}_{\mathrm{o}}$, and ${R}_{\text{vessel}}$. Prior to calculating the SSE, e first normalize both the experiment and model by their respective total areas [i.e., $\iint p({r}_{\mathrm{s}},\lambda )\text{\hspace{0.17em}}\mathrm{d}{r}_{\mathrm{s}}\text{\hspace{0.17em}}\mathrm{d}\lambda $]. In this way, the fitting is performed by comparing the relative shapes of $p({r}_{\mathrm{s}},\lambda )$ across different wavelengths and is therefore less sensitive to temporal intensity fluctuations.

In order to determine the location of the minimum SSE, we use the built-in MATLAB function “lsqnonlin” to perform a six-dimensional (i.e., one dimension for each physical property) bounded nonlinear least-squares minimization. The bounds are chosen such that each step of the optimization stays within the physiologically relevant range of optical/physical properties and within the bounds of our model library. The bounds are as follows: $D\in [\mathrm{2,4}]$, $v\in [\mathrm{0,1}]$, ${S}_{\mathrm{o}}\in [\mathrm{0,1}]$, and ${R}_{\text{vessel}}\in [\mathrm{0,200}]\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, with ${\sigma}_{n}^{2}$ and ${L}_{n}$ chosen such that ${l}_{\mathrm{s}}^{*}\in [0.5,2.5]$ mm and $g\in [0.7,0.96]$.

One of the main problems associated with optimization problems is discriminating between local and global minima. To combat this difficulty, we implement a two-part algorithm to first quickly find appropriate seed values and subsequently thoroughly fit $p({r}_{\mathrm{s}},\lambda )$ to arrive at the global minimum. In the first step, we use sets of preselected seed physical parameter values spread throughout the bounds of our $p({r}_{\mathrm{s}},\lambda )$ library and perform a coarse minimization with a total of 25 iterations. We repeat this process with six unique sets of physical properties and establish the trial with the lowest SSE as the most likely neighborhood in which to find the global minimum. In the second step, we use the preliminary endpoint values as seeds for a more thorough minimization that performs a maximum of 500 iterations to arrive at the global minimum.

To test the accuracy and stability of our inverse algorithm, we first generated 1200 sets of physical properties within the bounds of our model parameter space and generated the corresponding $p({r}_{\mathrm{s}},\lambda )$ curves. We subsequently added Gaussian noise with zero mean to the $p({r}_{\mathrm{s}},\lambda )$ curves and assessed the accuracy of the inverse algorithm with different noise variance ${\sigma}_{\text{noise}}^{2}$ levels. In performing this test, we used system parameters, which we have the capabilities to measure using EBS: ${r}_{\mathrm{s}}=[\mathrm{40,2000}]\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, $\mathrm{d}{r}_{\mathrm{s}}\text{\hspace{0.17em}}=20\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, $\lambda \text{\hspace{0.17em}}=[0.500,0.700]\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, $\mathrm{d}\lambda \text{\hspace{0.17em}}\text{\hspace{0.17em}}=0.020\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, illumination $\text{spot}\text{\hspace{0.17em}}\text{diameter}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}4\text{\hspace{0.17em}}\mathrm{mm}$, and using the linear copolarization channel. Figure 7 demonstrates the accuracy and stability of our inverse algorithm. Figures 7(a)–7(f) show the excellent correlations between the actual and calculated values for each of the six physical properties using a noise level corresponding to that expected in a typical EBS measurement (i.e., ${\sigma}_{\text{noise}}^{2}\text{\hspace{0.17em}}\text{\hspace{0.17em}}={10}^{-13}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu {\mathrm{m}}^{-2}$). The calculated ${R}^{2}$ values are $>0.99$ for ${A}_{n}$, 0.97 for ${L}_{n}$, 0.99 for $D$, $>0.99$ for $v$, 0.97 for ${S}_{\mathrm{o}}$, and 0.95 for ${R}_{\text{vessel}}$. Figure 7(g) shows the degradation of algorithm accuracy for increasing levels of additive noise variance. Within the range of typical EBS system noise levels, all ${R}^{2}$ values are greater than 0.53, with this minimum value occurring in the ${L}_{n}$ parameter. The reduced accuracy in the ${L}_{n}$ parameter is attributed to the fact that the shape of changes less dramatically across ${z}_{\mathrm{max}}=0.5{l}_{\mathrm{s}}^{*}$ than it does for the other five physical properties. We note that one output of our fitting routine is ${\sigma}_{\text{noise}}^{2}$. Thus, for an experimental measurement, the investigator can gauge the stability of their results based on ${\sigma}_{\text{noise}}^{2}$ and can take measures to decrease the noise (i.e., through longer integration time or higher laser power) to obtain more reliable results. Figure 7(h) shows a demonstration of how $p({r}_{\mathrm{s}})$ would look for differing levels of ${\sigma}_{\mathrm{noise}}^{2}$.

Closing out the validation of our inverse algorithm, Fig. 8 provides estimates of the extracted value variability for varying levels of system noise. In this case, we quantify the variability as the standard deviation of the error (i.e., $\text{error}=\text{actual}\text{\hspace{0.17em}}\text{value}-\text{measured}\text{\hspace{0.17em}}\text{value}$) over multiple different realizations. As expected, for increasing levels of noise, the variability of extracted physical properties also increases. The level of variability is essentially independent of physical property value (data not shown). For a typical experimental noise level with ${\sigma}_{\text{noise}}^{2}={10}^{-13}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu {\mathrm{m}}^{-2}$, the standard deviation in extracted values is $1.15\times {10}^{-6}$ for ${A}_{n}$, $0.28\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ for ${L}_{n}$, 0.05 for $D$, 0.004 for $v$, 0.05 for ${S}_{\mathrm{o}}$, and $6.08\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ for ${R}_{\text{vessel}}$.

## 3.5.

### Measuring $p({x}_{\mathrm{s}},{y}_{\mathrm{s}})$ Using Enhanced Backscattering Spectroscopy

Detailed discussion of the nature of the EBS phenomenon can be found in a number of seminal publications.^{42}43.44.45.^{–}^{46} Additionally, details about the application of EBS to biological tissue can be found in Ref. 5, while discussion about the measurement of $p({x}_{\mathrm{s}},{y}_{\mathrm{s}})$ using EBS can be found in Ref. 31 In brief, the measured angular EBS peak ${I}_{\mathrm{EBS}}({\theta}_{x},{\theta}_{y})$ is the Fourier transform of the product of five functions.

## (15)

$$\begin{array}{c}{I}_{\mathrm{EBS}}({\theta}_{x},{\theta}_{y})\end{array}={\iint}_{-\infty}^{\infty}p({x}_{\mathrm{s}},{y}_{\mathrm{s}})\xb7pc({x}_{\mathrm{s}},{y}_{\mathrm{s}})\xb7s({x}_{\mathrm{s}},{y}_{\mathrm{s}})\xb7c({x}_{\mathrm{s}},{y}_{\mathrm{s}})\xb7mtf({x}_{\mathrm{s}},{y}_{\mathrm{s}}){e}^{ik({x}_{\mathrm{s}}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\theta}_{x}+{y}_{\mathrm{s}}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\theta}_{y})}\text{\hspace{0.17em}}\mathrm{d}{x}_{\mathrm{s}}\text{\hspace{0.17em}}\mathrm{d}{y}_{\mathrm{s}},\phantom{\rule{0ex}{0ex}}=\mathcal{F}\{p({x}_{\mathrm{s}},{y}_{\mathrm{s}})\xb7pc({x}_{\mathrm{s}},{y}_{\mathrm{s}})\xb7s({x}_{\mathrm{s}},{y}_{\mathrm{s}})\xb7c({x}_{\mathrm{s}},{y}_{\mathrm{s}})\xb7mtf({x}_{\mathrm{s}},{y}_{\mathrm{s}})\},$$Although EBS data are measured in angular space, we return the data to spatial coordinates by simply computing an inverse Fourier transform:

## (16)

$$\begin{array}{c}{p}_{\mathrm{eff}}\end{array}=p\xb7pc\xb7s\xb7c\xb7mtf\phantom{\rule{0ex}{0ex}}={\mathcal{F}}^{-1}\{{I}_{\mathrm{EBS}}\},$$It is possible to take the analysis one step further and directly extract function $p$ by dividing function ${p}_{\mathrm{eff}}$ by the four remaining modulating functions. Indeed, we have previously performed such an analysis in Refs. 5 and 47. However, such a procedure amounts to a deconvolution process that only serves to amplify noise. In this work, we therefore limit our solution of the inverse problem to fitting of function ${p}_{\mathrm{eff}}$.

## 3.6.

### Experimental Enhanced Backscattering Spectroscopy Tissue Instrumentation and Measurements

All measurements were taken on the EBS system detailed in Ref. 48. In brief, broadband illumination from a supercontinuum laser (NKT Photonics: SuperK EXW-6) passes into a variable bandpass filter (NKT photonics: SuperK Varia) to select wavelengths between 0.500 and $0.700\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ in steps of $0.020\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ with $0.010\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ bandwidth. The spectrally filtered light is then coupled into a single mode fiber, recollimated using a 100-mm focal length lens, and truncated to a 4-mm-diameter spot size using an iris aperture. The collimated light passes through a linear polarizer and is directed onto the sample. Backscattered light is detected using a $50/50$ plate beam splitter. The backscattered light then passes through a copolarized linear analyzer and is focused onto a $-70\xb0\mathrm{C}$ cooled CCD camera (Princeton Instruments, PIXIS 1024B eXcelon) using a 200-mm Fourier lens. A full set of spectral measurements, including calibrations, can currently be performed in under 20 min with further potential speed optimizations possible.

*Ex-vivo* measurements were taken from a single sacrificed rat within 8 h of organ extraction. All animal procedures were reviewed and approved by the Institutional Animal Care and Use Committee at Northwestern University. Since our model is only strictly valid for the index matching case, all samples were submerged in an index matching liquid during measurement. Assuming a mean tissue RI of 1.38, our index matching liquid used a mixture of 33% glycerol to 67% deionized water by volume.^{1} The resulting RI of the mixture was confirmed using a calibrated refractometer (Reichert, Abbe Mark III). While the actual mean RI for an arbitrary tissue could differ slightly from this value, the difference in the reflection coefficient for any range of values expected in biological tissue would be negligible. In order to reduce speckle and obtain a smooth ensemble average measurement, samples were slowly rotated at $\sim 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{rotation}/\mathrm{min}$ using an automated circular rotation stage (Zaber, T-RSW).

## 4.

## Results

## 4.1.

### Demonstration of the Shape of $p({r}_{\mathrm{s}})$ for Varying Physical Properties

We now use the $p({r}_{\mathrm{s}})$ model generating algorithm described in Sec. 3.3 to demonstrate the effect of the six physical parameters on the shape of $p({r}_{\mathrm{s}})$. Furthermore, these demonstrations serve to illustrate the quality and speed of the $p({r}_{\mathrm{s}})$ model generating algorithm. Each figure in this section, therefore, includes an estimate of the time needed to create all of the displayed data.

We begin by analyzing the effect on $p({r}_{\mathrm{s}},\lambda )$ for varying values of the $D$ parameter. Figure 9 shows a spectrum of $p({r}_{\mathrm{s}})$ for $\lambda =0.400$ to 0.700 in $0.100\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ steps in a purely scattering medium (i.e., $v=0$) with varying $D=2.0$, 3.0, and 4.0. For each panel, the values of ${\sigma}_{n}^{2}$ and ${L}_{n}$ are fixed such that ${l}_{\mathrm{s}}^{*}=1000\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ and $g=0.9$ at $\lambda =0.600\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. The observed changes in the shape of $p({r}_{\mathrm{s}})$ for varying values of $D$ are due to two factors. First, the value of $D$ determines, in part, the shape of the scattering phase function and therefore manifests itself as a change in reflectance at subdiffusion lengthscales (i.e., ${r}_{\mathrm{s}}<{l}_{\mathrm{s}}^{*}$).^{19} Second, since ${\mu}_{\mathrm{s}}^{*}\propto {\lambda}^{D-4}$,^{21} the height and width of $p({r}_{\mathrm{s}})$ within the diffusion regime varies more strongly for $D=2$ than for $D=4$.

Next, we briefly discuss how varying ${L}_{n}$ and ${\sigma}_{n}^{2}$ would be represented in the $p({r}_{\mathrm{s}})$ curves. Corresponding figures are not included in this work for the sake of brevity but can easily be studied using the code posted on our laboratory website.^{25} The primary effect of changing ${L}_{n}$ is to alter the shape of the scattering phase function and thereby change $g$ and ${\mu}_{\mathrm{s}}^{*}$. With increasing ${L}_{n}$, the scattering of light becomes more forward directed, resulting in a higher value of $g$. This change in the shape of the phase function is once again represented as modifications in the shape of $p({r}_{\mathrm{s}})$ at subdiffusion lengthscales (data not shown). The parameter ${\sigma}_{n}^{2}$ is directly proportional to ${\mu}_{\mathrm{s}}$ and therefore ${\mu}_{\mathrm{s}}^{*}$. As a result, its only effect is to scale $p({r}_{\mathrm{s}})$ according to the change in ${\mu}_{\mathrm{s}}^{*}$ without changing its general shape or spectral dependence (data not shown).

We study the effects of the three microvascular properties on the shape of $p({r}_{\mathrm{s}},\lambda )$ in Figs. 10 and 11. Figure 10 shows spectra of $p({r}_{\mathrm{s}},\lambda )$ for $\lambda =0.400$ to 0.700 in $0.100\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ steps for varying $v$. The values of $v$ shown in the title of Figs. 10(a)–10(c) are chosen such that the maximum value of ${\mu}_{\mathrm{a}}/{\mu}_{\mathrm{s}}^{*}=0.1$, 0.5, and 1.0, respectively. The other two microvascular properties are fixed such that ${S}_{\mathrm{o}}=0.5$ and ${R}_{\text{vessel}}=0$. The ultrastructural properties have a fixed $D=3.0$ and values of ${\sigma}_{n}^{2}$ and ${L}_{n}$ chosen such that ${l}_{\mathrm{s}}^{*}=1000\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ and $g=0.9$ at $\lambda =0.600\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. With increasing values of $v$, $p({r}_{\mathrm{s}})$ is increasingly attenuated at larger values ${r}_{\mathrm{s}}$. This occurs because $v$ is directly proportional to ${\mu}_{\mathrm{a}}$ and photons exiting at larger ${r}_{\mathrm{s}}$ travel through larger pathlengths, which are more likely to encounter absorbers that attenuate reflected light. The largest attenuation of $p({r}_{\mathrm{s}},\lambda )$ occurs at $\lambda =0.550\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, a value where both oxy and deoxyhemoglobin are strongly absorbing.

The final demonstration of the $p({r}_{\mathrm{s}},\lambda )$ model generating algorithm is shown in Fig. 11. In this figure, we show the integrated intensity versus wavelength $I(\lambda )={\int}_{0}^{2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}}p({r}_{\mathrm{s}},\lambda )\mathrm{d}{r}_{\mathrm{s}}$ for varying values of ${S}_{\mathrm{o}}$ and ${R}_{\text{vessel}}$ with 1-nm wavelength resolution. The effect of ${R}_{\text{vessel}}$ can be seen by looking at the shading of the curves in each of the three panels of Fig. 11. With increasing values of ${R}_{\text{vessel}}$, the blood vessel shading effect becomes more profound and manifests itself as more muted changes in the hemoglobin absorption band between $0.550-0.600\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ of illumination. Moving from panel a to panel c, the value of ${S}_{\mathrm{o}}$ increases from 0 to 1, representing a shift toward more oxygenated hemoglobin. Within the wavelength range from 0.500 to $0.700\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, the difference between oxy- and deoxyhemoglobin can be identified by observing the absorption peaks occurring at 0.540, 0.555, and $0.576\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. For completely deoxygenated hemoglobin, the absorption level is highest at $0.555\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, which results in a corresponding minimum of $I(\lambda )$. For completely oxygenated hemoglobin, the absorption level is highest for two peaks occurring at 0.540 and $0.576\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, which again result in a corresponding minimum of $I(\lambda )$.

## 4.2.

### Application of Inverse Algorithm to Experimental Tissue Measurements

As a demonstration of the applicability of our inverse model to various biological tissues, we acquired $p({r}_{\mathrm{s}},\lambda )$ measurements from *ex-vivo* rat liver, stomach, and heart tissues using EBS. All measurements were completed within 8 h of sacrificing the animal. Given that the thickness of each tissue sample was $>1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{cm}$ (i.e., thickness $>10\text{\hspace{0.17em}}{l}_{\mathrm{s}}^{*}$), we assumed a semi-infinite tissue thickness when running the inverse algorithm.

Figure 12 shows the resulting fits of the inverse model for measurements taken from the outside surface of an intact rat liver (first row), stomach (second row), and heart (third row). Figures 12(a), 12(e), and 12(i) show the experimentally measured $p({r}_{\mathrm{s}},\lambda )$ with the thickness of each line indicating the standard error of 20 measurements and the color corresponding to the illumination $\lambda $. As expected, the green and yellow wavelengths, which are much more susceptible to hemoglobin absorption, show a quicker decay in reflectance with increasing ${r}_{\mathrm{s}}$. In order to better observe this spectral shape, Figs. 12(b), 12(f), and 12(j) show the $I(\lambda )$ spectrum obtained by integrating $p({r}_{\mathrm{s}},\lambda )$ over ${r}_{\mathrm{s}}$. The characteristic hemoglobin absorption dip is visible from $\sim 0.500$ to $0.600\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. To demonstrate the match between experiment and model, Figs. 12(c), 12(g), and 12(k) show the linear fit correlation between corresponding points on $p({r}_{\mathrm{s}},\lambda )$. In addition to a qualitative match between the experiment and fit, we note a very high ${R}^{2}$ value of 0.99 for each tissue sample. Finally, Figs. 12(d), 12(h), and 12(l) provide average image depictions of $p({r}_{\mathrm{s}},\lambda )$ with the comparison between experimental (top) and model (bottom).

Summarizing the results for the three measured tissues, Table 5 lists the extracted values of the six physical parameters as well as the ${R}^{2}$ value of the model fit and the experimental noise variance ${\sigma}_{\text{noise}}^{2}$. The values are given as the $\text{mean}\pm \text{standard error}$ over 20 measurements. We note that the ${R}^{2}$ values are excellent at $>0.99$ for each fit. Furthermore, the relatively low levels of ${\sigma}_{\text{noise}}^{2}$ indicate that our experimental measurements fall within a range where the inverse model performs admirably. The measured standard deviations are roughly consistent with those reported in Sec. 3.4. Thus, the standard deviations are largely representative of the system noise, as opposed to tissue heterogeneity.

## Table 5

Extracted ultrastructural and microvascular properties.

Tissue | σn2 (×10−4) | Ln (μm) | D | v (%) | Rvessel (μm) | So | R2 | σnoise2μm−2 |
---|---|---|---|---|---|---|---|---|

Liver | $1.36\pm 0.22$ | $0.49\pm 0.12$ | $2.91\pm 0.07$ | $15.99\pm 0.51$ | $23.01\pm 1.43$ | $13.84\pm 1.36$ | 0.99 | $8.86\times {10}^{-14}$ |

Stomach | $3.32\pm 0.61$ | $8.47\pm 3.28$ | $2.40\pm 0.07$ | $4.30\pm 0.48$ | $18.16\pm 3.64$ | $51.44\pm 4.16$ | 0.99 | $2.48\times {10}^{-13}$ |

Heart | $1.05\pm 0.03$ | $0.44\pm 0.02$ | $3.64\pm 0.04$ | $11.48\pm 0.39$ | $6.22\pm 1.41$ | $7.71\pm 2.06$ | 0.99 | $1.82\times {10}^{-13}$ |

Remarking on the physical properties themselves, the reported values of ${\sigma}_{n}^{2}$ fall within the expected range for tissue compiled by Çapoğlu et al. in Ref. 15 (i.e., ${\sigma}_{n}^{2}\sim 0.5\times {10}^{-4}$ to $5\times {10}^{-4}$). The value of ${L}_{n}$ varies by more than an order of magnitude between the three tissue types but falls within the expected range of $\sim 0.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ to $\sim 10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$.^{12} For $D$, we measured ultrastructures with spatial distributions corresponding to fractals (liver and stomach) and stretched exponential (heart). With regards to microvasculature, we measured increasing $v$ from stomach to heart to liver tissues. Since the tissues were washed prior to measurement, these value reflect the blood trapped within the tissue, as opposed to superficial surface blood. For ${R}_{\text{vessel}}$, the smallest sizes corresponding to a mix of capillaries and small arterioles/venules were found in the heart, whereas the liver and stomach showed larger sizes corresponding to medium-sized arterioles/veins. Finally, for ${S}_{\mathrm{o}}$, we measured values that were much more deoxygenated (i.e., $<\sim 0.50$) than would be expected in live tissue. We attribute this to the *ex-vivo* nature of the measurements.^{47}

Although ${\sigma}_{n}^{2}$, ${L}_{n}$, and $D$ provide useful parametrization of tissue ultrastructure, it is often conceptually enlightening to examine the full shape of ${B}_{n}({r}_{d})$ as well as the corresponding random media representations. Figure 13(a) shows the extracted shape of ${B}_{n}({r}_{d})$ for measurements of liver, stomach, and heart tissues. The values of ${B}_{n}$ at ${r}_{d}={r}_{\mathrm{min}}=2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ correspond directly to the values of ${\sigma}_{n}^{2}$ are shown in Table 5. Thus, the curves for liver and heart tissues begin at roughly the same value, while the curve for stomach is about three times higher. The curve shapes at larger ${r}_{d}$ are then determined by ${L}_{n}$ and $D$. For the liver, the curve shape is a fractal with dimension $D=2.91$ (corresponding roughly to the Henyey–Greenstein case at $D=3.0$) until the upper lengthscale of fractality reaches ${L}_{n}=0.49\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. Similarly, for the stomach, the curve is a fractal with $D=2.40$ for lengthscales up to ${L}_{n}=8.47\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. For heart tissue, the curve shape is stretched exponential for structural lengthscales smaller than ${L}_{n}=0.44\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. Figures 13(b)–13(d) show the random media representations of liver, stomach, and heart tissues, respectively. Because of the relatively higher value of ${\sigma}_{n}^{2}$ for stomach tissue, the stomach realization exhibits a higher contrast than that of either liver or heart tissue. Comparing the liver and heart tissue realizations, we see that although the contrast is similar, the heart tissue possesses larger structural features than that of liver tissue, directly mirroring the shape of ${B}_{n}$ in Fig. 13(a).

Having established the physical properties, which give rise to the optical signal, we can now analytically determine the optical properties for arbitrary wavelengths. Using the equations derived by Rogers et al., we convert ${\sigma}_{n}^{2}$, ${L}_{n}$, and $D$ into the spectral dependence of ${\mu}_{\mathrm{s}}^{*}$, $g$, and $\gamma $.^{21} Additionally, using Eq. (12), we convert $v$, ${L}_{n}$, and ${R}_{\text{vessel}}$ into the spectral dependence of ${\mu}_{\mathrm{a}}$. Figure 14 summarizes the spectral dependencies of various optical coefficients for liver, stomach, and heart tissues over the fitted range of wavelengths between 0.500 and $0.700\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$.

For liver tissue, ${\mu}_{\mathrm{s}}^{*}(\lambda )$ exhibits a power-law decay from $11.92\pm 0.73$ ($\text{mean}\pm \text{standard error}$) ${\mathrm{cm}}^{-1}$ at $0.500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ to $7.52\pm 0.53\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ at $0.700\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, while ${\mu}_{\mathrm{a}}$ exhibits the characteristic hemoglobin absorption spectrum with large absorption occurring between wavelengths from 0.500 to $0.600\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. Comparing the effects of absorption and scattering, we note that absorption dominates from 0.500 to $0.600\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, whereas scattering dominates from 0.600 to $0.700\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. Moving to the phase function shape parameters, Fig. 14(b) shows a slowly decaying $g(\lambda )$ from a value of $0.84\pm 0.01$ at $0.500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ to $0.79\pm 0.02$ at $0.700\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ and Fig. 14(c) shows a decaying $\gamma (\lambda )$ from a value of $1.82\pm 0.02$ at $0.500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ to $1.75\pm 0.02$ at $0.700\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. The decaying trend in $g$ is expected since as the ratio of $\lambda /{L}_{n}$ becomes larger, the scattering phase function will become more isotropic, eventually reaching the Rayleigh scattering regime with $g=0$.

For stomach tissue, Fig. 14(d) shows ${\mu}_{\mathrm{s}}^{*}(\lambda )$ decaying from $10.20\pm 0.68\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ at $0.500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ to $5.77\pm 0.37\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ at $0.700\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. Figure 14(e) shows a very slightly decaying $g(\lambda )$ from $0.87\pm 0.01$ at $0.500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ to $0.84\pm 0.02$ at $0.700\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. Similarly, Fig. 14(f) shows a very slowly decaying $\gamma (\lambda )$ from $1.67\pm 0.03$ at $0.500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ to $1.64\pm 0.03$ at $0.700\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$.

For heart tissue, Fig. 14(g) shows $(\lambda )$ decaying from $17.84\pm 0.30\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ at $0.500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ to $13.70\pm 0.20\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ at $0.700\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. Figure 14(h) shows a decaying $g(\lambda )$ from $0.95\pm 0.002$ at $0.500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ to $0.92\pm 0.004$ at $0.700\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. Finally, Fig. 14(f) shows a very slowly decaying $\gamma (\lambda )$ from $2.28\pm 0.01$ at $0.500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ to $2.19\pm 0.01$ at $0.700\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$.

## 5.

## Discussion and Conclusions

The three main focuses of this work are (1) to summarize a unified model for calculating scattering and absorption properties using tissue ultrastructure and microvasculature as inputs, (2) to present an empirical Monte Carlo model for quickly generating subdiffusion regime tissue spatial reflectance profiles $p({r}_{\mathrm{s}})$, and (3) to apply this model toward an inverse algorithm capable of quickly and accurately extracting tissue ultrastructural and microvascular properties. We note that all MATLAB codes required for accomplishing these three goals are posted on our laboratory website.^{25}

Toward the first and second goals, we employed a model of scattering that represents biological tissue as a continuous random medium described by the Whittle–Matérn family of correlation functions^{21}^{,}^{26}^{,}^{27} and a model of absorption that assumes packaging of hemoglobin within blood vessels.^{7}8.^{–}^{9} Under these assumptions, we conducted a series of Monte Carlo simulations for a total of over 1 million processor-hours using an open-source code previously described in Ref. 35. We subsequently smoothed and compressed this vast dataset from $\sim 6$ gigabytes down to $\sim 70$ megabytes for easier portability and improved loading speed. Finally, we developed an algorithm to quickly expand and interpolate this dataset into experimentally observable $p({r}_{\mathrm{s}},\lambda )$ curves with tissue relevant optical properties in a matter of a few seconds to a few minutes depending on desired sample parameters.

Toward the third goal, we developed an inverse algorithm that uses the rapidly generated $p({r}_{\mathrm{s}},\lambda )$ curves to extract six physical tissue parameters from a given experimental measurement. This algorithm uses a gradient search to minimize the SSE between the model and experimental data. In order to find the global minimum, we perform a cursory search with six different initial seeds values and subsequently use the trial with the lowest error to perform a more rigorous search that determines the final values of the physical properties. Depending on the noise level and computer processor speed, a typical $p({r}_{\mathrm{s}},\lambda )$ measurement can be fit in between 3 and 7 min. Furthermore, under typical noise levels expected while using EBS to measure $p({r}_{\mathrm{s}},\lambda )$, the inverse algorithm provides accurate quantification of physical properties with an average ${R}^{2}>0.90$ between actual and measured values.

Applying the inverse algorithm to experimental data, we measured physical properties from *ex-vivo* rat liver, stomach, and heart tissues. For each of these samples, the model fit was excellent with ${R}^{2}>0.99$. This suggests the applicability of our model to solid tissues such as these examples. Despite these encouraging results, it is recommended that the applicability of the Born approximation and Whittle–Matérn model be considered prior to implementing our inverse algorithm. For instance, tissues with spherical form factors (e.g., adipose tissue) may be better suited to the application of Mie theory. Furthermore, the user should also interpret the validity of the extracted physical parameters based on the ${R}^{2}$ value of the fit and the noise level of their measurement.

A number of limitations and assumptions should be reiterated. First, as stated earlier, the validity of our model relies on assumptions that the Whittle–Matérn and hemoglobin vessel packaging models adequately describe the tissue under investigation. For most tissue types, there is evidence that these models offer a good approximation of underlying tissue structure.^{8}^{,}^{12}^{,}^{22}^{,}^{32} However, there may be limitations for (1) tissues with spherical scattering geometries or (2) anisotropic tissue structures such as bone or muscle. We note that all samples presented in this work are rotated such that structural anisotropy is averaged away. Furthermore, under the current form of the Whittle–Matérn model, dispersion of RI is assumed to be negligible, meaning the shape of ${B}_{n}({r}_{d})$ does not change with wavelength. A hypothesis in favor of this assumption is that while local RI varies with wavelength, on average across many different molecular species, the variance and distribution of excess RI that leads to scattering remains relatively stable. Given the practical difficulty of measuring the dispersion of RI in biological tissue with nanoscale resolution, we use indirect methods such as the excellent model fits demonstrated in this work to corroborate this hypothesis. A second main assumption is that the tissue under investigation is statistically homogeneous and therefore does not change with depth. Although this is not strictly true for most tissue types, we argue that the physical properties of most tissues do not change drastically within the superficial depths interrogated at subdiffusion transport lengthscales. Moreover, although it is possible in principle to incorporate a multilayered inverse model, in practice, it becomes very difficult to fit for the rapidly expanding number of free parameters added with each new layer. A third limitation is that the current model is only valid for media in which there is no RI contrast at the boundary. Therefore, in order to achieve the most accurate results, tissues must be submerged within an index matching liquid at the time of the measurement. A fourth limitation of the current model is that absorption analysis is limited to quantification of hemoglobin. For many tissue types and spectral measurement regimes, it is safe to assume that hemoglobin is the dominant absorber. Still, we note that our model can easily be extended to include the influence of other biological absorbers, such as melanin, bilirubin, fat content, and water content based on their respective molar extinction coefficients.^{11}

Although a primary endpoint of many methods of optical tissue characterization is determining optical coefficients (${\mu}_{\mathrm{s}}$, ${\mu}_{\mathrm{a}}$, $g$, $\gamma $, etc.) or other empirical parameters, we have argued that a more complete understanding of tissue structure can be gained by going one step further and extracting physical tissue properties. We have therefore provided a framework from which to relate empirical reflectance measurements back to their underlying ultrastructural and microvascular properties. Despite the noted limitations, the models and algorithms presented in this work can be extended to other applications to help bring improved quantification of fundamental tissue structure to the field of optical characterization.

## Acknowledgments

This study was supported by National Institutes of Health under Grant Nos. RO1CA128641 and R01EB003682. A.J.R. was supported by a National Science Foundation Graduate Research Fellowship under Grant No. DGE-0824162. The simulations in this paper were made possible by a supercomputing allocation grant from Northwestern University’s Quest high-performance computing resource. We would like to thank Dr. Mitra Hartmann and Chris Bresee for providing rat tissue samples.

## References

*In vivo*measurement of the shape of the tissue-refractive-index correlation function and its application to detection of colorectal field carcinogenesis,” J. Biomed. Opt. 17(4), 047005 (2012). http://dx.doi.org/10.1117/1.JBO.17.4.047005 Google Scholar

*In vivo*risk analysis of pancreatic cancer through optical characterization of duodenal mucosa,” Pancreas 44(5), 735–741 (2015).PANCE40885-3177 http://dx.doi.org/10.1097/MPA.0000000000000340 Google Scholar

*in vivo*,” Med. Phys. 19(4), 879–888 (1992).MPHYA60094-2405 http://dx.doi.org/10.1118/1.596777 Google Scholar