# Subdiffusion reflectance spectroscopy to measure tissue ultrastructure and microvasculature: model and inverse algorithm

**Andrew J. Radosevich**,

**Adam Eshein**,

**The-Quyen Nguyen**,

**Vadim Backman**

Northwestern University, Biomedical Engineering, Tech E310, 2145 Sheridan Road, Evanston, Illinois 60208, United States

*J. Biomed. Opt*. 20(9), 097002 (Sep 28, 2015). doi:10.1117/1.JBO.20.9.097002

#### Open Access

## Abstract

**Abstract.**
Reflectance measurements acquired from within the subdiffusion regime (i.e., lengthscales smaller than a transport mean free path) retain much of the original information about the shape of the scattering phase function. Given this sensitivity, many models of subdiffusion regime light propagation have focused on parametrizing the optical signal through various optical and empirical parameters. We argue, however, that a more useful and universal way to characterize such measurements is to focus instead on the fundamental physical properties, which give rise to the optical signal. This work presents the methodologies that used to model and extract tissue ultrastructural and microvascular properties from spatially resolved subdiffusion reflectance spectroscopy measurements. We demonstrate this approach using *ex-vivo* rat tissue samples measured by enhanced backscattering spectroscopy.

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 $ls*$) preserve information about the shape of the scattering phase function $P(\theta )$.^{1}^{–}^{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 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 a(\lambda )$ provides a vehicle from which to extract more physical information about microvasculature, such as blood volume fraction $v$, hemoglobin oxygen saturation $So$, and blood vessel radius $Rvessel$.^{7}^{–}^{10} Should the specific value of $\mu 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}^{–}^{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, $Bn(rd)$.^{16}^{–}^{18} If $Bn(rd)$ 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 s$ and anisotropy factor $g$ as well as higher order shape parameters, such as $\gamma $, $D$, etc.^{10}^{,}^{19} Since $Bn(rd)$ 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 $Bn(rd)$ 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 $Bn(rd)$, 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 $Bn(rd)$, 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}^{–}^{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.

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}^{–}^{9} The full medium model is characterized by six physical sample properties: the variance of RI $\sigma n2$, the characteristic lengthscale of RI heterogeneity $Ln$, the shape parameter of RI distribution $D$, the blood volume fraction $v$, the oxygen saturation $So$, and the blood vessel radius $Rvessel$. As we will demonstrate, these six “physical” properties directly determine the wavelength dependence of four commonly measured “optical” properties: the scattering coefficient $\mu s$, the absorption coefficient $\mu a$, the reduced scattering coefficient $\mu 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.

$\mu s$ | Function of $\sigma n2$, $Ln$, $D$, and $\lambda $ |

$g$ | Function of $Ln$, $D$, and $\lambda $ |

$\mu s*$ | Function of $\sigma n2$, $Ln$, $D$, and $\lambda $ |

$\mu a$ | Function of $v$, $So$, $Rvessel$, and $\lambda $ |

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 $cos2(\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 $Bn(rd)$:^{4}^{,}^{14}^{,}^{21}^{,}^{27}^{,}^{31}^{,}^{32}

^{12}

^{33}

^{,}

^{34}

Figure 2(a) shows examples of $Bn(rd)$ 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 $\Phi s(k\u2192s)$ is the three-dimensional Fourier transform of $Bn(r\u2192d)$.^{17} Applying the Whittle–Matérn model, which assumes spherical symmetry, $\Phi s(ks)$ can be found as^{31}

Defining scattering for arbitrary polarization, the amplitude scattering matrix per unit volume polarizability $s(\theta )$ can be found by incorporating the dipole scattering pattern into $\Phi s(ks)$:^{35}

^{18}

^{18}

^{,}

^{36}

^{30}These can be calculated as

Analytical equations for $\mu 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 s*$, which is primarily relevant for multiple scattering samples, is defined as $\mu s*=\mu s\xb7(1\u2212g)$. Note that the scattering mean free path $ls$ and the transport scattering mean free path $ls*$ are simply the inverse of $\mu s$ and $\mu s*$, respectively.

Furthermore, it should be noted that the widely recognized power-law decay in $\mu s*(\lambda )$ is not just an empirical observation but has physical origins in the shape of $Bn(rd)$. In the limit of high $kLn$ (satisfied in most tissues), $\mu s*\u221d\lambda D\u22124$ for $D\u22644$. For $D>4$, the value of $\mu s*$ will be constant across wavelength. Thus, quantifying the power of the decay in $\mu 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 n2(kLn)2\u226a1$ 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 $ls$. In other words, the Born approximation can be accurately applied to all tissues for which the concept of optical properties, such as $ls$, $\mu s$, $\mu s*$, etc., can be applied (i.e., most tissues that are studied in the biomedical optics community).

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 $\u223c5$ to $50\u2009\u2009\mu 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}^{–}^{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 a*$ can be found as

^{37}Finally, using the formalism presented by Van Veen et al., $Cpack$ can be calculated as

Incorporating the effects of hemoglobin absorption into our calculations introduces three physical properties: $v$, $So$, and the product $CHb\xb7Rvessel$. Note that within this formalism, the effect of $CHb$ cannot be separated from that of $Rvessel$ without additional information. As a result, a common assumption is that $CHb=150\u2009\u2009g/L$, allowing for the isolation of $Rvessel$.^{9}^{,}^{37}^{,}^{38}

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(rs)$] at subdiffusion lengthscales using enhanced backscattering.^{5}

Modeling of $p(rs)$ 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(rs)$ 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 $100ls*$ 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 ($rs$) and their maximum penetration depth ($zmax$). For samples with finite thickness, integration over the $zmax$ dimension from 0 to the tissue thickness provides a way to arrive at the measurable $p(rs)$ 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 $(+\u2212)$ polarization channels. The radial resolution of these grids is $4\xb710\u22123rs/ls*$ with total extent of $4rs/ls*$, while the depth resolution is $4\xb710\u22123zmax/ls*$ with total extent of $2zmax/ls*$. 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 $109$ photon histories have been calculated.

In order to develop a library of $p(rs,zmax)$ for tissue relevant optical properties, we performed a series of simulations with varying values of $g$, $D$, and $\mu a/\mu s*$. These values are summarized in Table 2. To achieve different values of $g$, we fixed $\lambda =0.633\u2009\u2009\mu m$ and varied $Ln$. 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.

$g$ | 0.7 to 0.96 in 0.02 steps |

$D$ | 2.0 to 4.0 in 0.1 steps |

$\mu a/\mu 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.

After running simulations with $109$ photons, the calculated $p(rs,zmax)$ 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 $\u223c50$ 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(rs,zmax)$ data shown in Fig. 4(a) is smoothed by first calculating the cumulative summation over rows (i.e., the $rs$ 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(rs,zmax)$. 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(rs,zmax)$.

**F4 :**

Comparison between the raw and smoothed $p(rs,zmax)$ data: (a) raw $p(rs,zmax)$ in $log10$ scale; (b) smoothed $p(rs,zmax)$ in $log10$ scale; (c) $p(rs)$ distribution achieved by integrating $p(rs,zmax)$ over the $zmax$ dimension. The legend indicates integration up to $zmax=0.5ls*$, $1.0ls*$, and $\u221e$. Symbols indicate points sampled from the raw data and the solid lines show the smoothed data.

The same smoothing process is then repeated over the columns (i.e., the $rs$ dimension) of $p(rs,zmax)$. 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 $\u223c6$ gigabytes to $\u223c70$ megabytes for each polarization channel, vastly improving data portability.

The polynomial coefficients can be quickly expanded into the final smoothed version of the $p(rs,zmax)$ array shown in Fig. 4(b) using matrix multiplication:

A qualitatively good match between the raw and smoothed $p(rs)$ 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 ($R2$) 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(rs,zmax)$ grid as well as the 1-D $p(rs)$ distribution obtained by integrating over all $zmax$ (i.e., $p(rs)=\u222b0\u221ep(rs,zmax)dzmax$). The mean $R2$ values are summarized in Table 3. A mean $R2$ value of $>0.94$ for $p(rs,zmax)$ for each polarization channel indicates an accurate fit. The slight reduction in $R2$ from the ideal value of 1 is due to the small amount of noise present in the raw Monte Carlo simulations.

Polarization channel | $(xx)$ | $(xy)$ | $(++)$ | $(+\u2212)$ |
---|---|---|---|---|

Mean $R2$ for $p(rs,zmax)$ | 0.989 | 0.944 | 0.977 | 0.987 |

Mean $R2$ for $p(rs)$ | 0.999 | 0.993 | 0.996 | 0.999 |

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

The process of generating $p(rs)$ 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(rs\xb7\mu s*,zmax\xb7\mu s*,D,g,\mu a/\mu 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 a/\mu s*$. Next, we loop through each of the input wavelengths to compute its specific $p(rs)$. In this process, we first integrate $p(rs\xb7\mu s*,zmax\xb7\mu s*,D,g,\mu a/\mu s*)$ over the $zmax$ dimension up to the input sample thickness yielding a four-dimensional grid of $p(rs\xb7\mu s*,zmax\xb7\mu s*,D,g,\mu a/\mu s*)$. Finally, we do a multidimensional cubic interpolation to find $p(rs)$ for the values of $\mu s*$, $D$, $g$, and $\mu a$ at the specified wavelength. Using a 2.5-GHz Intel Core i5-3210M processor, an entire $p(rs)$ spectrum for a nonabsorbing sample over the range $\lambda =0.500\u2009\u2009\mu m$ to $0.700\u2009\u2009\mu m$ in $0.020\u2009\u2009\mu m$ steps can be calculated in under 150 ms.

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

Polarization channel | $(xx)$ | $(xy)$ | $(++)$ | $(+\u2212)$ |
---|---|---|---|---|

Median $R2$ for $p(rs)$ | 0.999 | $>0.999$ | 0.999 | $>0.999$ |

Minimum $R2$ for $p(rs)$ | 0.9865 | 0.996 | 0.9784 | 0.984 |

Our inverse scattering algorithm is performed by minimizing the SSE between an experimental measurement of $p(rs,\lambda )$ and our model for $p(rs,\lambda )$ to obtain $\sigma n2$, $Ln$, $D$, $v$, $So$, and $Rvessel$. Prior to calculating the SSE, e first normalize both the experiment and model by their respective total areas [i.e., $\u222cp(rs,\lambda )\u2009drs\u2009d\lambda $]. In this way, the fitting is performed by comparing the relative shapes of $p(rs,\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\u2208[2,4]$, $v\u2208[0,1]$, $So\u2208[0,1]$, and $Rvessel\u2208[0,200]\u2009\u2009\mu m$, with $\sigma n2$ and $Ln$ chosen such that $ls*\u2208[0.5,2.5]$ mm and $g\u2208[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(rs,\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(rs,\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(rs,\lambda )$ curves. We subsequently added Gaussian noise with zero mean to the $p(rs,\lambda )$ curves and assessed the accuracy of the inverse algorithm with different noise variance $\sigma noise2$ levels. In performing this test, we used system parameters, which we have the capabilities to measure using EBS: $rs=[40,2000]\u2009\u2009\mu m$, $drs\u2009=20\u2009\u2009\mu m$, $\lambda \u2009=[0.500,0.700]\u2009\u2009\mu m$, $d\lambda \u2009\u2009=0.020\u2009\u2009\mu m$, illumination $spot\u2009diameter\u2009=\u20094\u2009mm$, 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 noise2\u2009\u2009=10\u221213\u2009\u2009\mu m\u22122$). The calculated $R2$ values are $>0.99$ for $An$, 0.97 for $Ln$, 0.99 for $D$, $>0.99$ for $v$, 0.97 for $So$, and 0.95 for $Rvessel$. 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 $R2$ values are greater than 0.53, with this minimum value occurring in the $Ln$ parameter. The reduced accuracy in the $Ln$ parameter is attributed to the fact that the shape of changes less dramatically across $zmax=0.5ls*$ than it does for the other five physical properties. We note that one output of our fitting routine is $\sigma noise2$. Thus, for an experimental measurement, the investigator can gauge the stability of their results based on $\sigma noise2$ 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(rs)$ would look for differing levels of $\sigma noise2$.

**F7 :**

Inverse fitting algorithm performance for semi-infinite medium. 1200 $p(rs,\lambda )$ curves were generated with random values for (a) $An$, (b) $Ln$, (c) $D$, (d) $v$, (e) $So$, and (f) $Rvessel$. The noise variance for panels a to f is $10\u221213$. (g) The degradation of the inverse fitting algorithm with increasing levels of noise variance $\sigma noise2$. The shaded region indicates the range of noise levels expected for a typical enhanced backscattering spectroscopy experiment. (f) Example $p(rs)$ curves for $\sigma noise2$ equal to $10\u221212\u2009\u2009\mu m\u22122$ and $10\u221214\u2009\u2009\mu m\u22122$.

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., $error=actual\u2009value\u2212measured\u2009value$) 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 noise2=10\u221213\u2009\u2009\mu m\u22122$, the standard deviation in extracted values is $1.15\xd710\u22126$ for $An$, $0.28\u2009\u2009\mu m$ for $Ln$, 0.05 for $D$, 0.004 for $v$, 0.05 for $So$, and $6.08\u2009\u2009\mu m$ for $Rvessel$.

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

Although EBS data are measured in angular space, we return the data to spatial coordinates by simply computing an inverse Fourier transform:

It is possible to take the analysis one step further and directly extract function $p$ by dividing function $peff$ 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 $peff$.

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\u2009\u2009\mu m$ in steps of $0.020\u2009\u2009\mu m$ with $0.010\u2009\u2009\mu 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 $\u221270\xb0C$ 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 $\u223c1\u2009\u2009rotation/min$ using an automated circular rotation stage (Zaber, T-RSW).

We now use the $p(rs)$ model generating algorithm described in Sec. 3.3 to demonstrate the effect of the six physical parameters on the shape of $p(rs)$. Furthermore, these demonstrations serve to illustrate the quality and speed of the $p(rs)$ 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(rs,\lambda )$ for varying values of the $D$ parameter. Figure 9 shows a spectrum of $p(rs)$ for $\lambda =0.400$ to 0.700 in $0.100\u2009\u2009\mu 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 n2$ and $Ln$ are fixed such that $ls*=1000\u2009\u2009\mu m$ and $g=0.9$ at $\lambda =0.600\u2009\u2009\mu m$. The observed changes in the shape of $p(rs)$ 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., $rs<ls*$).^{19} Second, since $\mu s*\u221d\lambda D\u22124$,^{21} the height and width of $p(rs)$ within the diffusion regime varies more strongly for $D=2$ than for $D=4$.

**F9 :**

Purely scattering $p(rs,\lambda )$ spectrum for varying $D$ with the values of $\sigma n2$ and $Ln$ fixed such that $ls*=1000\u2009\u2009\mu m$ and $g=0.9$ at $\lambda =0.600\u2009\u2009\mu m$. The $pxx$ notation on the $y$-axis indicates curves generated for the linear copolarized channel: (a) $p(rs)$ for $D=2.0$. (b) $p(rs)$ for $D=3.0$. (c) $p(rs)$ for $D=4.0$. In each panel, the arrow indicates the direction of increasing $\lambda $. Total calculation time for all data in this figure is under 5 s on a 2.5-GHz processor with 8 GB of RAM.

Next, we briefly discuss how varying $Ln$ and $\sigma n2$ would be represented in the $p(rs)$ 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 $Ln$ is to alter the shape of the scattering phase function and thereby change $g$ and $\mu s*$. With increasing $Ln$, 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(rs)$ at subdiffusion lengthscales (data not shown). The parameter $\sigma n2$ is directly proportional to $\mu s$ and therefore $\mu s*$. As a result, its only effect is to scale $p(rs)$ according to the change in $\mu 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(rs,\lambda )$ in Figs. 10 and 11. Figure 10 shows spectra of $p(rs,\lambda )$ for $\lambda =0.400$ to 0.700 in $0.100\u2009\u2009\mu 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 a/\mu s*=0.1$, 0.5, and 1.0, respectively. The other two microvascular properties are fixed such that $So=0.5$ and $Rvessel=0$. The ultrastructural properties have a fixed $D=3.0$ and values of $\sigma n2$ and $Ln$ chosen such that $ls*=1000\u2009\u2009\mu m$ and $g=0.9$ at $\lambda =0.600\u2009\u2009\mu m$. With increasing values of $v$, $p(rs)$ is increasingly attenuated at larger values $rs$. This occurs because $v$ is directly proportional to $\mu a$ and photons exiting at larger $rs$ travel through larger pathlengths, which are more likely to encounter absorbers that attenuate reflected light. The largest attenuation of $p(rs,\lambda )$ occurs at $\lambda =0.550\u2009\u2009\mu m$, a value where both oxy and deoxyhemoglobin are strongly absorbing.

**F10 :**

$p(rs)$ spectrum for varying $v$ with the values of $\sigma n2$ and $Ln$ fixed such that $ls*=1000\u2009\u2009\mu m$ and $g=0.9$ at $\lambda =0.600\u2009\u2009\mu m$: (a) $p(rs)$ for $v=1\xb70.006$ and maximum $\mu a/\mu s*=0.1$; (b) $p(rs)$ for $v=5\xb70.006$ and maximum $\mu a/\mu s*=0.5$; (c) $p(rs)$ for $v=10\xb70.006$ and maximum $\mu a/\mu s*=1.0$. Total calculation time for all data in this figure is under 8 s on a 2.5-GHz processor with 8 GB of RAM.

**F11 :**

Integrated intensity versus wavelength (1-nm resolution) for varying $So$ and $Rvessel$. These plots hold constant $D=3.0$ and $v=0.02$ with the values of $\sigma n2$ and $Ln$ chosen such that $ls*=1000\u2009\u2009\mu m$ and $g=0.9$ at $\lambda =0.600\u2009\u2009\mu m$. The intensity $I$ on the $y$-axis is calculated by integrating over the $rs$ dimension (i.e., $I=\u222b02000p(rs)\u2009drs$): $I$ for (a) $So=0.0$; (b) $So=0.5$; and (c) $So=1.0$. In each panel, we calculated $I$ for $Rvessel=0$, 100, and $2000\u2009\u2009\mu m$. Total calculation time for all data in this figure is under 4 min on a 2.5-GHz processor with 8 GB of RAM.

The final demonstration of the $p(rs,\lambda )$ model generating algorithm is shown in Fig. 11. In this figure, we show the integrated intensity versus wavelength $I(\lambda )=\u222b02\u2009\u2009mmp(rs,\lambda )drs$ for varying values of $So$ and $Rvessel$ with 1-nm wavelength resolution. The effect of $Rvessel$ can be seen by looking at the shading of the curves in each of the three panels of Fig. 11. With increasing values of $Rvessel$, the blood vessel shading effect becomes more profound and manifests itself as more muted changes in the hemoglobin absorption band between $0.550\u22120.600\u2009\u2009\mu m$ of illumination. Moving from panel a to panel c, the value of $So$ increases from 0 to 1, representing a shift toward more oxygenated hemoglobin. Within the wavelength range from 0.500 to $0.700\u2009\u2009\mu m$, the difference between oxy- and deoxyhemoglobin can be identified by observing the absorption peaks occurring at 0.540, 0.555, and $0.576\u2009\u2009\mu m$. For completely deoxygenated hemoglobin, the absorption level is highest at $0.555\u2009\u2009\mu 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\u2009\u2009\mu m$, which again result in a corresponding minimum of $I(\lambda )$.

As a demonstration of the applicability of our inverse model to various biological tissues, we acquired $p(rs,\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\u2009\u2009cm$ (i.e., thickness $>10\u2009ls*$), 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(rs,\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 $rs$. 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(rs,\lambda )$ over $rs$. The characteristic hemoglobin absorption dip is visible from $\u223c0.500$ to $0.600\u2009\u2009\mu 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(rs,\lambda )$. In addition to a qualitative match between the experiment and fit, we note a very high $R2$ value of 0.99 for each tissue sample. Finally, Figs. 12(d), 12(h), and 12(l) provide average image depictions of $p(rs,\lambda )$ with the comparison between experimental (top) and model (bottom).

**F12 :**

Comparison of experimental $p(rs,\lambda )$ with inverse model fit for *ex-vivo* rat (a–d) liver, (e–h) stomach, and (i–l) heart tissues. The first column (panels a, e, and i) shows the experimental $p(rs)$ for 11 different wavelengths. The color of each line indicates the visible color of the corresponding illumination $\lambda $. The thickness of each line indicates the standard error over 20 measurements. The second column (panels b, f, and j) shows the $I(\lambda )$ intensity spectrum obtained by integrating $p(rs,\lambda )$ over the $rs$ dimension. Symbols show the experiment with standard error and the solid line is the fit. The third column (panels c, g, and k) shows the linear correlation for corresponding points on $p(rs,\lambda )$ between experiment and the model fit. Values on the $x$- and $y$-axes have been multiplied by $105$. The final column (panels d, h, and l) provides image depictions of $p(rs,\lambda )$ measured experimentally (top) as well as the model fit (bottom).