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,²⁸ 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³⁰ was developed by adding a ${\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).³¹ 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.¹² 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.

Fig. 1

Examples of the continuous distributions of refractive index (RI) specified by the Whittle–Matérn model for (a) $D=2.0$, (b) $D=4.0$, and (c) $D=6.0$.

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}

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¹⁶ (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)$.¹⁷ Applying the Whittle–Matérn model, which assumes spherical symmetry, ${\mathrm{\Phi }}_{\mathrm{s}}(k_{\mathrm{s}})$ can be found as³¹

Fig. 2

Examples of the functional forms specified by the Whittle–Matérn model for $D=2.0$, $D=4.0$, and $D=6.0$: (a) the normalized $B_n(r_d)$; (b) the normalized ${\mathrm{\Phi }}_{\mathrm{s}}(k_{\mathrm{s}})$; (c) example of $\sigma (\theta )$ for $\lambda =0.633\mu \mathrm{m}$, $A_n=1$, and $L_n=1\mu \mathrm{m}$.

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}})$:³⁵

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}}·(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.¹⁵ 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\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–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

Incorporating the effects of hemoglobin absorption into our calculations introduces three physical properties: $v$, $S_{\mathrm{o}}$, and the product $C_{\mathrm{Hb}}·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\mathrm{g}/\mathrm{L}$, allowing for the isolation of $R_{\text{vessel}}$.^9,37,38

3.1.

Monte Carlo Simulation and $p(r_{\mathrm{s}},z_{\max })$ Library Population

Modeling of $p(r_{\mathrm{s}})$ is achieved through solution of the radiative transfer equation (RTE).³⁹ Given the complexity of solving the RTE, analytical solutions rely on numerous simplifying assumptions, which typically result in loss of accuracy at subdiffusion lengthscales.⁴⁰ Although improved analytical solutions, such as the phase function corrected diffusion approximation by Vitkin et al.¹ and modified spherical harmonic method by Liemert and Kienle⁴¹, 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.³⁵ The code for this software is open source and can be downloaded on our laboratory website.²⁵ Below, we summarize the relevant simulation parameters.

The simulated scattering medium is a slab of statistically homogeneous material with a total thickness of $100l_{\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_{\max }$). For samples with finite thickness, integration over the $z_{\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·{10}^{-3}{r}_{\mathrm{s}}/l_{\mathrm{s}}^{*}$ with total extent of $4r_{\mathrm{s}}/l_{\mathrm{s}}^{*}$, while the depth resolution is $4·{10}^{-3}{z}_{\max }/l_{\mathrm{s}}^{*}$ with total extent of $2z_{\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_{\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\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.

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_{\max })$ Data

After running simulations with ${10}^9$ photons, the calculated $p(r_{\mathrm{s}},z_{\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_{\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_{\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_{\max })$.

Fig. 3

Flowchart of smoothing and compression algorithm.

Fig. 4

Comparison between the raw and smoothed $p(r_{\mathrm{s}},z_{\max })$ data: (a) raw $p(r_{\mathrm{s}},z_{\max })$ in ${\log }_{10}$ scale; (b) smoothed $p(r_{\mathrm{s}},z_{\max })$ in ${\log }_{10}$ scale; (c) $p(r_{\mathrm{s}})$ distribution achieved by integrating $p(r_{\mathrm{s}},z_{\max })$ over the $z_{\max }$ dimension. The legend indicates integration up to $z_{\max }=0.5l_{\mathrm{s}}^{*}$, $1.0l_{\mathrm{s}}^{*}$, and $\infty$. 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 $r_{\mathrm{s}}$ dimension) of $p(r_{\mathrm{s}},z_{\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_{\max })$ array shown in Fig. 4(b) using matrix multiplication:

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_{\max })$ grid as well as the 1-D $p(r_{\mathrm{s}})$ distribution obtained by integrating over all $z_{\max }$ (i.e., $p(r_{\mathrm{s}})={\int }_0^{\infty }p(r_{\mathrm{s}},z_{\max })\mathrm{d}{z}_{\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_{\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.

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}}·{\mu }_{\mathrm{s}}^{*},z_{\max }·{\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}}·{\mu }_{\mathrm{s}}^{*},z_{\max }·{\mu }_{\mathrm{s}}^{*},D,g,{\mu }_{\mathrm{a}}/{\mu }_{\mathrm{s}}^{*})$ over the $z_{\max }$ dimension up to the input sample thickness yielding a four-dimensional grid of $p(r_{\mathrm{s}}·{\mu }_{\mathrm{s}}^{*},z_{\max }·{\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\mu \mathrm{m}$ to $0.700\mu \mathrm{m}$ in $0.020\mu \mathrm{m}$ steps can be calculated in under 150 ms.

Fig. 5

Flowchart of the $p(r_{\mathrm{s}},\lambda )$ model generating algorithm.

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.

Fig. 6

Comparison between raw and generated data with the lowest $R^2$ values: (a) normalized $p(r_{\mathrm{s}})$, Symbols indicate points sampled from the raw data and the solid lines show the generated data; (b) residuals between raw and generated data.

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 )\mathrm{d}{r}_{\mathrm{s}}\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}]\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}]\mu \mathrm{m}$, $\mathrm{d}{r}_{\mathrm{s}}=20\mu \mathrm{m}$, $\lambda =[0.500,0.700]\mu \mathrm{m}$, $\mathrm{d}\lambda =0.020\mu \mathrm{m}$, illumination $\text{spot}\text{diameter}=4\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={10}^{-13}\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_{\max }=0.5l_{\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$.

Fig. 7

Inverse fitting algorithm performance for semi-infinite medium. 1200 $p(r_{\mathrm{s}},\lambda )$ curves were generated with random values for (a) $A_n$, (b) $L_n$, (c) $D$, (d) $v$, (e) $S_{\mathrm{o}}$, and (f) $R_{\text{vessel}}$. The noise variance for panels a to f is ${10}^{-13}$. (g) The degradation of the inverse fitting algorithm with increasing levels of noise variance ${\sigma }_{\text{noise}}^2$. The shaded region indicates the range of noise levels expected for a typical enhanced backscattering spectroscopy experiment. (f) Example $p(r_s)$ curves for ${\sigma }_{\text{noise}}^2$ equal to ${10}^{-12}\mu {\mathrm{m}}^{-2}$ and ${10}^{-14}\mu {\mathrm{m}}^{-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{value}-\text{measured}\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}\mu {\mathrm{m}}^{-2}$, the standard deviation in extracted values is $1.15\times {10}^{-6}$ for $A_n$, $0.28\mu \mathrm{m}$ for $L_n$, 0.05 for $D$, 0.004 for $v$, 0.05 for $S_{\mathrm{o}}$, and $6.08\mu \mathrm{m}$ for $R_{\text{vessel}}$.

Fig. 8

Inverse fitting algorithm variability across different measurement noise levels. Each panel shows the standard deviation (STD) of the measurement error for (a) $A_n$, (b) $L_n$, (c) $D$, (d) $v$, (e) $S_{\mathrm{o}}$, and (f) $R_{\text{vessel}}$; 1200 realizations were performed at each noise level.

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

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 $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\mu \mathrm{m}$ in steps of $0.020\mu \mathrm{m}$ with $0.010\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°\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.¹ 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{rotation}/\min$ using an automated circular rotation stage (Zaber, T-RSW).

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\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\mu \mathrm{m}$ and $g=0.9$ at $\lambda =0.600\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}}^{*}$).¹⁹ Second, since ${\mu }_{\mathrm{s}}^{*}\propto {\lambda }^{D-4}$,²¹ the height and width of $p(r_{\mathrm{s}})$ within the diffusion regime varies more strongly for $D=2$ than for $D=4$.

Fig. 9

Purely scattering $p(r_{\mathrm{s}},\lambda )$ spectrum for varying $D$ with the values of ${\sigma }_n^2$ and $L_n$ fixed such that $l_{\mathrm{s}}^{*}=1000\mu \mathrm{m}$ and $g=0.9$ at $\lambda =0.600\mu \mathrm{m}$. The $p_{xx}$ notation on the $y$-axis indicates curves generated for the linear copolarized channel: (a) $p(r_{\mathrm{s}})$ for $D=2.0$. (b) $p(r_{\mathrm{s}})$ for $D=3.0$. (c) $p(r_{\mathrm{s}})$ 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 $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.²⁵ 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\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\mu \mathrm{m}$ and $g=0.9$ at $\lambda =0.600\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\mu \mathrm{m}$, a value where both oxy and deoxyhemoglobin are strongly absorbing.

Fig. 10

$p(r_{\mathrm{s}})$ spectrum for varying $v$ with the values of ${\sigma }_n^2$ and $L_n$ fixed such that $l_{\mathrm{s}}^{*}=1000\mu \mathrm{m}$ and $g=0.9$ at $\lambda =0.600\mu \mathrm{m}$: (a) $p(r_{\mathrm{s}})$ for $v=1·0.006$ and maximum ${\mu }_{\mathrm{a}}/{\mu }_{\mathrm{s}}^{*}=0.1$; (b) $p(r_{\mathrm{s}})$ for $v=5·0.006$ and maximum ${\mu }_{\mathrm{a}}/{\mu }_{\mathrm{s}}^{*}=0.5$; (c) $p(r_{\mathrm{s}})$ for $v=10·0.006$ and maximum ${\mu }_{\mathrm{a}}/{\mu }_{\mathrm{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.

Fig. 11

Integrated intensity versus wavelength (1-nm resolution) for varying $S_{\mathrm{o}}$ and $R_{\text{vessel}}$. These plots hold constant $D=3.0$ and $v=0.02$ with the values of ${\sigma }_n^2$ and $L_n$ chosen such that $l_{\mathrm{s}}^{*}=1000\mu \mathrm{m}$ and $g=0.9$ at $\lambda =0.600\mu \mathrm{m}$. The intensity $I$ on the $y$-axis is calculated by integrating over the $r_{\mathrm{s}}$ dimension (i.e., $I={\int }_0^{2000}p(r_{\mathrm{s}})\mathrm{d}{r}_{\mathrm{s}}$): $I$ for (a) $S_{\mathrm{o}}=0.0$; (b) $S_{\mathrm{o}}=0.5$; and (c) $S_{\mathrm{o}}=1.0$. In each panel, we calculated $I$ for $R_{\text{vessel}}=0$, 100, and $2000\mu \mathrm{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(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\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\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\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\mu \mathrm{m}$. For completely deoxygenated hemoglobin, the absorption level is highest at $0.555\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\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\mathrm{cm}$ (i.e., thickness $>10l_{\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\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).

Fig. 12

Comparison of experimental $p(r_{\mathrm{s}},\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(r_{\mathrm{s}})$ 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(r_{\mathrm{s}},\lambda )$ over the $r_{\mathrm{s}}$ 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(r_{\mathrm{s}},\lambda )$ between experiment and the model fit. Values on the $x$- and $y$-axes have been multiplied by ${10}^5$. The final column (panels d, h, and l) provides image depictions of $p(r_{\mathrm{s}},\lambda )$ measured experimentally (top) as well as the model fit (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.

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\mu \mathrm{m}$ to $\sim 10\mu \mathrm{m}$.¹² 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.⁴⁷

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_{\min }=2\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\mu \mathrm{m}$. Similarly, for the stomach, the curve is a fractal with $D=2.40$ for lengthscales up to $L_n=8.47\mu \mathrm{m}$. For heart tissue, the curve shape is stretched exponential for structural lengthscales smaller than $L_n=0.44\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).

Fig. 13

Tissue ultrastructure representations. (a) Extracted $B_n(r_d)$ curves for liver, stomach, and heart tissues. The thickness of each curve corresponds to the standard error of 20 measurements. The corresponding random media representations of (b) liver, (c) stomach, and (d) heart. The color map is scaled the same for each image and represents the excess refractive index value.

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$.²¹ 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\mu \mathrm{m}$.

Fig. 14

Expanded optical properties corresponding to the physical properties in Table 5. The first column (panels a, d, and g) shows the spectral dependence of ${\mu }_{\mathrm{s}}^{*}$ and ${\mu }_{\mathrm{a}}$ from 0.500 to $0.700\mu \mathrm{m}$. The second column (panels b, e, and h) shows the spectral dependence of $g$. The third column (panels c, f, and i) shows the spectral dependence of the $\gamma$ parameter. The thickness of each line indicates the standard error over 20 measurements.

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\mu \mathrm{m}$ to $7.52\pm 0.53{\mathrm{cm}}^{-1}$ at $0.700\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\mu \mathrm{m}$. Comparing the effects of absorption and scattering, we note that absorption dominates from 0.500 to $0.600\mu \mathrm{m}$, whereas scattering dominates from 0.600 to $0.700\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\mu \mathrm{m}$ to $0.79\pm 0.02$ at $0.700\mu \mathrm{m}$ and Fig. 14(c) shows a decaying $\gamma (\lambda )$ from a value of $1.82\pm 0.02$ at $0.500\mu \mathrm{m}$ to $1.75\pm 0.02$ at $0.700\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{\mathrm{cm}}^{-1}$ at $0.500\mu \mathrm{m}$ to $5.77\pm 0.37{\mathrm{cm}}^{-1}$ at $0.700\mu \mathrm{m}$. Figure 14(e) shows a very slightly decaying $g(\lambda )$ from $0.87\pm 0.01$ at $0.500\mu \mathrm{m}$ to $0.84\pm 0.02$ at $0.700\mu \mathrm{m}$. Similarly, Fig. 14(f) shows a very slowly decaying $\gamma (\lambda )$ from $1.67\pm 0.03$ at $0.500\mu \mathrm{m}$ to $1.64\pm 0.03$ at $0.700\mu \mathrm{m}$.

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