Research Papers: Sensing

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

[+] Author Affiliations
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
History: Received June 23, 2015; Accepted August 31, 2015
Text Size: A A A

Open Access Open Access

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.

Figures in this Article

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(θ).16 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(θ), 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 μa(λ) 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 μa(λ) 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.710 Should the specific value of μ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.1214 Assuming the validity of the Born approximation (a safe assumption for most solid tissues15), the fundamental physical tissue characteristic, which gives rise to light scattering, is the spatial refractive index (RI) autocorrelation function, Bn(rd).1618 If Bn(rd) is known, a series of simple mathematical transformations can be used to first derive P(θ) and subsequently optical properties, such as the scattering coefficient μs and anisotropy factor g as well as higher order shape parameters, such as γ, 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,2224

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,λ) 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,λ), 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,λ); and (4) to announce MATLAB codes posted online25 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 model21,26,27 and absorption assuming packaging of hemoglobin within blood vessels.79 The full medium model is characterized by six physical sample properties: the variance of RI σ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 μs, the absorption coefficient μa, the reduced scattering coefficient μ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 Grahic Jump Location
Table 1Relation between optical and physical parameters.
Relating Scattering Properties to σn2, Ln, 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 HGPF30 was developed by adding a cos2(θ) 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.

Graphic Jump Location
Fig. 1
F1 :

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 Bn(rd):4,14,21,27,31,32Display Formula

where Kν(·) is the modified Bessel function of the second kind with order ν, An is the fluctuation strength of RI, Ln is the characteristic length of heterogeneity, and D determines the shape of the distribution (e.g., power-law/fractal for D<3, decaying exponential for D=4, Gaussian as D, etc.). We further note that when D=3, the resulting phase function equation is identical to the modified HGPF. In purely mathematical terms, An is simply an amplitude factor, which stretches Bn up and down along the y-axis, Ln stretches Bn left and right along the x-axis, and D changes the functional form of the distribution. In large part, An is simply a normalization factor that has limited physical meaning outside of the Whittle–Matérn model functional form. We, therefore, convert to the more intuitive RI variance σn2 as12Display Formula
where rmin is the minimum structural length-scale of fractality. Here, we approximate rmin=2nm, 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 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 approximation16 (also known as the Rayleigh–Gans–Debye theory), the power spectral density Φs(ks) is the three-dimensional Fourier transform of Bn(rd).17 Applying the Whittle–Matérn model, which assumes spherical symmetry, Φs(ks) can be found as31Display Formula

where ks equals 2knosin(θ/2), k is the freespace wavenumber, and no is the mean medium RI. Figure 2(b) shows examples of Φs for varying D.

Graphic Jump Location
Fig. 2
F2 :

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 Bn(rd); (b) the normalized Φs(ks); (c) example of σ(θ) for λ=0.633μm, An=1, and Ln=1μm.

Defining scattering for arbitrary polarization, the amplitude scattering matrix per unit volume polarizability s(θ) can be found by incorporating the dipole scattering pattern into Φs(ks):35Display Formula

Expanding the electric fields into intensities, the scattering phase function for arbitrary polarization is18Display Formula
where E and E are the orthogonal components of the electric field. Next, the differential scattering cross-section per unit volume for unpolarized light σ(θ,φ) can be calculated as18,36Display Formula
where θ is the scattering angle and φ is the azimuthal angle. Conceptually, function σ(θ,φ) specifies the directionality of scattered light intensity for a single scattering event. We note that since unpolarized light is rotationally symmetric about φ, the right-hand side of Eq. (6) is independent of φ. Figure 2(c) shows examples of σ(θ) for varying D. For arbitrary σ(θ), the scattering coefficient μs is defined by integrating over all solid angles: Display Formula
and the anisotropy factor g by calculating the first moment of σ(cosθ): Display Formula
Similarly, if desired, higher order parameters characterizing σ can also be calculated. Two such parameters include the second moment g2 and the shape parameter γ.30 These can be calculated as Display Formula
Display Formula

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

Furthermore, it should be noted that the widely recognized power-law decay in μs*(λ) 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), μs*λD4 for D4. For D>4, the value of μs* will be constant across wavelength. Thus, quantifying the power of the decay in μs*(λ) 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 σn2(kLn)21 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, μs, μs*, etc., can be applied (i.e., most tissues that are studied in the biomedical optics community).

Relating Absorption Properties to v, So, Rvessel

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 5 to 50μ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.79 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 μa* can be found as Display Formula

where Cpack is the blood vessel packaging correction factor, v is the blood volume fraction, and μa is the absorption coefficient of the whole blood contained with the blood vessels. The spectral dependence of μa can be found as Display Formula
where CHb is the concentration of hemoglobin, So is the hemoglobin oxygen saturation, εHbO2 is the molar extinction coefficient of oxygenated hemoglobin, and εHb is the molar extinction coefficient of deoxygenated hemoglobin. In this paper, we use values for εHbO2 and εHb taken from Prahl’s online database.37 Finally, using the formalism presented by Van Veen et al., Cpack can be calculated as Display Formula
where Rvessel is the blood vessel radius.

Incorporating the effects of hemoglobin absorption into our calculations introduces three physical properties: v, So, and the product CHb·Rvessel. 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=150g/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,λ) results, (2) quickly and accurately generate p(r,λ) using the prerun library of simulations, and (3) invert p(r,λ) 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

Monte Carlo Simulation and p(rs,zmax) Library Population

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 Kienle41, 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 (+) polarization channels. The radial resolution of these grids is 4·103rs/ls* with total extent of 4rs/ls*, while the depth resolution is 4·103zmax/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 μa/μs*. These values are summarized in Table 2. To achieve different values of g, we fixed λ=0.633μ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.

Table Grahic Jump Location
Table 2Summary 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.

Smoothing and Compression of Monte Carlo p(rs,zmax) Data

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 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(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).

Graphic Jump Location
Fig. 3
F3 :

Flowchart of smoothing and compression algorithm.

Graphic Jump Location
Fig. 4
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 . 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 6 gigabytes to 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(rs,zmax) array shown in Fig. 4(b) using matrix multiplication: Display Formula

where bi are the coefficients for the i’th polynomial and L indicates the maximum zmax extent of the grid.

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)=0p(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.

Table Grahic Jump Location
Table 3Summary of the R2 values between raw and smoothed data.
Rapid Modeling of p(rs,λ) from Polynomial Library

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,λ) model generating algorithm that enables modeling with nearly analytical speed and accuracy. The algorithm inputs are the desired values of rs, λ, 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·μs*,zmax·μs*,D,g,μa/μ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 μa/μs*. Next, we loop through each of the input wavelengths to compute its specific p(rs). In this process, we first integrate p(rs·μs*,zmax·μs*,D,g,μa/μs*) over the zmax dimension up to the input sample thickness yielding a four-dimensional grid of p(rs·μs*,zmax·μs*,D,g,μa/μs*). Finally, we do a multidimensional cubic interpolation to find p(rs) for the values of μs*, D, g, and μ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 λ=0.500μm to 0.700μm in 0.020μm steps can be calculated in under 150 ms.

Graphic Jump Location
Fig. 5
F5 :

Flowchart of the p(rs,λ) model generating algorithm.

In order to validate the accuracy of the p(rs) generating algorithm, we performed a series of 30 simulations with randomized σn2, Ln, D, and μ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.

Table Grahic Jump Location
Table 4Summary of the R2 values between raw and generated data.
Graphic Jump Location
Fig. 6
F6 :

Comparison between raw and generated data with the lowest R2 values: (a) normalized p(rs), Symbols indicate points sampled from the raw data and the solid lines show the generated data; (b) residuals between raw and generated data.

Inverse Algorithm to Extract σn2, Ln, D, v, So, and Rvessel

Our inverse scattering algorithm is performed by minimizing the SSE between an experimental measurement of p(rs,λ) and our model for p(rs,λ) to obtain σ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., p(rs,λ)drsdλ]. In this way, the fitting is performed by comparing the relative shapes of p(rs,λ) 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[2,4], v[0,1], So[0,1], and Rvessel[0,200]μm, with σn2 and Ln chosen such that ls*[0.5,2.5] mm and g[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,λ) 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,λ) 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,λ) curves. We subsequently added Gaussian noise with zero mean to the p(rs,λ) curves and assessed the accuracy of the inverse algorithm with different noise variance σnoise2 levels. In performing this test, we used system parameters, which we have the capabilities to measure using EBS: rs=[40,2000]μm, drs=20μm, λ=[0.500,0.700]μm, dλ=0.020μm, illumination spotdiameter=4mm, 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., σnoise2=1013μm2). 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 σnoise2. Thus, for an experimental measurement, the investigator can gauge the stability of their results based on σ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 σnoise2.

Graphic Jump Location
Fig. 7
F7 :

Inverse fitting algorithm performance for semi-infinite medium. 1200 p(rs,λ) 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 1013. (g) The degradation of the inverse fitting algorithm with increasing levels of noise variance σnoise2. The shaded region indicates the range of noise levels expected for a typical enhanced backscattering spectroscopy experiment. (f) Example p(rs) curves for σnoise2 equal to 1012μm2 and 1014μm2.

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=actualvaluemeasuredvalue) 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 σnoise2=1013μm2, the standard deviation in extracted values is 1.15×106 for An, 0.28μm for Ln, 0.05 for D, 0.004 for v, 0.05 for So, and 6.08μm for Rvessel.

Graphic Jump Location
Fig. 8
F8 :

Inverse fitting algorithm variability across different measurement noise levels. Each panel shows the standard deviation (STD) of the measurement error for (a) An, (b) Ln, (c) D, (d) v, (e) So, and (f) Rvessel; 1200 realizations were performed at each noise level.

Measuring p(xs,ys) Using Enhanced Backscattering Spectroscopy

Detailed discussion of the nature of the EBS phenomenon can be found in a number of seminal publications.4246 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(θx,θy) is the Fourier transform of the product of five functions. Display Formula

where F specifies the 2-D Fourier transform operation, p is the spatial reflectance profile, pc is the phase correlation function, s is a modulation due to finite illumination spot size, c is a modulation due to finite spatial coherence length, and mtf is the imaging system’s modulation transfer function. We remark that in the case of EBS, function p represents all multiply scattered rays (i.e., those undergoing two or more scattering events) that exit the medium in a direction antiparallel to the incident direction. Furthermore, in a number of special cases, Eq. (15) can be simplified. For the polarization preserving channels (i.e., linear copolarized and helicity preserving), function pc=1 and can therefore be neglected. For samples in which the illumination spot diameter is much greater than ls*, we can approximate s=1. Finally, when spatial coherent laser illumination is used, we can assume that the spatial coherence length is sufficiently large such that c=1.

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

where peff is the “effective” reflectance profile and the function arguments have been removed for simplicity.

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.

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μm in steps of 0.020μm with 0.010μ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°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 1rotation/min using an automated circular rotation stage (Zaber, T-RSW).

Demonstration of the Shape of p(rs) for Varying Physical Properties

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,λ) for varying values of the D parameter. Figure 9 shows a spectrum of p(rs) for λ=0.400 to 0.700 in 0.100μ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 σn2 and Ln are fixed such that ls*=1000μm and g=0.9 at λ=0.600μ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 μs*λD4,21 the height and width of p(rs) within the diffusion regime varies more strongly for D=2 than for D=4.

Graphic Jump Location
Fig. 9
F9 :

Purely scattering p(rs,λ) spectrum for varying D with the values of σn2 and Ln fixed such that ls*=1000μm and g=0.9 at λ=0.600μ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 λ. 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 σ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 μ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 σn2 is directly proportional to μs and therefore μs*. As a result, its only effect is to scale p(rs) according to the change in μ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,λ) in Figs. 10 and 11. Figure 10 shows spectra of p(rs,λ) for λ=0.400 to 0.700 in 0.100μ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 μa/μ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 σn2 and Ln chosen such that ls*=1000μm and g=0.9 at λ=0.600μm. With increasing values of v, p(rs) is increasingly attenuated at larger values rs. This occurs because v is directly proportional to μ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,λ) occurs at λ=0.550μm, a value where both oxy and deoxyhemoglobin are strongly absorbing.

Graphic Jump Location
Fig. 10
F10 :

p(rs) spectrum for varying v with the values of σn2 and Ln fixed such that ls*=1000μm and g=0.9 at λ=0.600μm: (a) p(rs) for v=1·0.006 and maximum μa/μs*=0.1; (b) p(rs) for v=5·0.006 and maximum μa/μs*=0.5; (c) p(rs) for v=10·0.006 and maximum μa/μ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.

Graphic Jump Location
Fig. 11
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 σn2 and Ln chosen such that ls*=1000μm and g=0.9 at λ=0.600μm. The intensity I on the y-axis is calculated by integrating over the rs dimension (i.e., I=02000p(rs)drs): 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μ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,λ) model generating algorithm is shown in Fig. 11. In this figure, we show the integrated intensity versus wavelength I(λ)=02mmp(rs,λ)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.5500.600μ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μm, the difference between oxy- and deoxyhemoglobin can be identified by observing the absorption peaks occurring at 0.540, 0.555, and 0.576μm. For completely deoxygenated hemoglobin, the absorption level is highest at 0.555μm, which results in a corresponding minimum of I(λ). For completely oxygenated hemoglobin, the absorption level is highest for two peaks occurring at 0.540 and 0.576μm, which again result in a corresponding minimum of I(λ).

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(rs,λ) 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 >1cm (i.e., thickness >10ls*), 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,λ) with the thickness of each line indicating the standard error of 20 measurements and the color corresponding to the illumination λ. 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(λ) spectrum obtained by integrating p(rs,λ) over rs. The characteristic hemoglobin absorption dip is visible from 0.500 to 0.600μ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,λ). 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,λ) with the comparison between experimental (top) and model (bottom).

Graphic Jump Location
Fig. 12
F12 :

Comparison of experimental p(rs,λ) 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 λ. The thickness of each line indicates the standard error over 20 measurements. The second column (panels b, f, and j) shows the I(λ) intensity spectrum obtained by integrating p(rs,λ) 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,λ) 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,λ) measured experimentally (top) as well as the model fit (bottom).