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



It is well established that reflectance measurements acquired from within the subdiffusion regime (i.e., source-detector separations smaller than a transport mean free path ls*) preserve information about the shape of the scattering phase function P(θ).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 1

Relation between optical and physical parameters.

μsFunction of σn2, Ln, D, and λ
gFunction of Ln, D, and λ
μs*Function of σn2, Ln, D, and λ
μaFunction of v, So, Rvessel, and λ


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.

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

Eq. (1)

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 as12

Eq. (2)

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 as31

Eq. (3)

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.

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 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):35

Eq. (4)

Expanding the electric fields into intensities, the scattering phase function for arbitrary polarization is18

Eq. (5)

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,36

Eq. (6)

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:

Eq. (7)

and the anisotropy factor g by calculating the first moment of σ(cosθ):

Eq. (8)

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

Eq. (9)


Eq. (10)


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

Eq. (11)

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

Eq. (12)

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

Eq. (13)

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


Materials and Methods

In this section, we begin by summarizing the aspects of our Monte Carlo simulations, which are relevant to the current study.35 We next detail the implementation of three algorithms, which allow us to (1) smooth and compress a large library of Monte Carlo simulated p(r,λ) 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 2

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

g0.7 to 0.96 in 0.02 steps
D2.0 to 4.0 in 0.1 steps
μa/μs*0, 0.25, 0.5, 0.75, 1, 2, and 3

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


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

Fig. 3

Flowchart of smoothing and compression algorithm.


Fig. 4

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:

Eq. (14)

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 3

Summary of the R2 values between raw and smoothed data.

Polarization channel(xx)(xy)(++)(+−)
Mean R2 for p(rs,zmax)0.9890.9440.9770.987
Mean R2 for p(rs)0.9990.9930.9960.999


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.

Fig. 5

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 4

Summary of the R2 values between raw and generated data.

Polarization channel(xx)(xy)(++)(+−)
Median R2 for p(rs)0.999>0.9990.999>0.999
Minimum R2 for p(rs)0.98650.9960.97840.984

Fig. 6

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.

Fig. 7

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.

Fig. 8

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.

Eq. (15)

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:

Eq. (16)

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.

Fig. 9

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.

Fig. 10

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.


Fig. 11

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

Fig. 12

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


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

Table 5

Extracted ultrastructural and microvascular properties.

Tissueσn2 (×10−4)Ln (μm)Dv (%)Rvessel (μm)SoR2σnoise2μm−2

Remarking on the physical properties themselves, the reported values of σn2 fall within the expected range for tissue compiled by Çapoğlu et al. in Ref. 15 (i.e., σn20.5×104 to 5×104). The value of Ln varies by more than an order of magnitude between the three tissue types but falls within the expected range of 0.5μm to 10μm.12 For D, we measured ultrastructures with spatial distributions corresponding to fractals (liver and stomach) and stretched exponential (heart). With regards to microvasculature, we measured increasing v from stomach to heart to liver tissues. Since the tissues were washed prior to measurement, these value reflect the blood trapped within the tissue, as opposed to superficial surface blood. For Rvessel, 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 So, we measured values that were much more deoxygenated (i.e., <0.50) than would be expected in live tissue. We attribute this to the ex-vivo nature of the measurements.47

Although σn2, Ln, and D provide useful parametrization of tissue ultrastructure, it is often conceptually enlightening to examine the full shape of Bn(rd) as well as the corresponding random media representations. Figure 13(a) shows the extracted shape of Bn(rd) for measurements of liver, stomach, and heart tissues. The values of Bn at rd=rmin=2nm correspond directly to the values of σn2 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 rd are then determined by Ln 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 Ln=0.49μm. Similarly, for the stomach, the curve is a fractal with D=2.40 for lengthscales up to Ln=8.47μm. For heart tissue, the curve shape is stretched exponential for structural lengthscales smaller than Ln=0.44μ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 σn2 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 Bn in Fig. 13(a).

Fig. 13

Tissue ultrastructure representations. (a) Extracted Bn(rd) 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 σn2, Ln, and D into the spectral dependence of μs*, g, and γ.21 Additionally, using Eq. (12), we convert v, Ln, and Rvessel into the spectral dependence of μ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μ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 μs* and μa from 0.500 to 0.700μ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 γ parameter. The thickness of each line indicates the standard error over 20 measurements.


For liver tissue, μs*(λ) exhibits a power-law decay from 11.92±0.73 (mean±standard error) cm1 at 0.500μm to 7.52±0.53cm1 at 0.700μm, while μa exhibits the characteristic hemoglobin absorption spectrum with large absorption occurring between wavelengths from 0.500 to 0.600μm. Comparing the effects of absorption and scattering, we note that absorption dominates from 0.500 to 0.600μm, whereas scattering dominates from 0.600 to 0.700μm. Moving to the phase function shape parameters, Fig. 14(b) shows a slowly decaying g(λ) from a value of 0.84±0.01 at 0.500μm to 0.79±0.02 at 0.700μm and Fig. 14(c) shows a decaying γ(λ) from a value of 1.82±0.02 at 0.500μm to 1.75±0.02 at 0.700μm. The decaying trend in g is expected since as the ratio of λ/Ln 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 μs*(λ) decaying from 10.20±0.68cm1 at 0.500μm to 5.77±0.37cm1 at 0.700μm. Figure 14(e) shows a very slightly decaying g(λ) from 0.87±0.01 at 0.500μm to 0.84±0.02 at 0.700μm. Similarly, Fig. 14(f) shows a very slowly decaying γ(λ) from 1.67±0.03 at 0.500μm to 1.64±0.03 at 0.700μm.

For heart tissue, Fig. 14(g) shows (λ) decaying from 17.84±0.30cm1 at 0.500μm to 13.70±0.20cm1 at 0.700μm. Figure 14(h) shows a decaying g(λ) from 0.95±0.002 at 0.500μm to 0.92±0.004 at 0.700μm. Finally, Fig. 14(f) shows a very slowly decaying γ(λ) from 2.28±0.01 at 0.500μm to 2.19±0.01 at 0.700μm.


Discussion and Conclusions

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

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

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

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

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

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


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



E. Vitkin et al., “Photon diffusion near the point-of-entry in anisotropically scattering turbid media,” Nat. Commun., 2 587 (2011). Google Scholar


A. Kienle, F. K. Forster and R. Hibst, “Influence of the phase function on determination of the optical properties of biological tissue by spatially resolved reflectance,” Opt. Lett., 26 (20), 1571 –1573 (2001). Google Scholar


A. Amelink and H. J. C. M. Sterenborg, “Measurement of the local optical properties of turbid media by differential path-length spectroscopy,” Appl. Opt., 43 (15), 3048 –3054 (2004). APOPAI 0003-6935 Google Scholar


A. J. Gomes et al., “In vivo measurement of the shape of the tissue-refractive-index correlation function and its application to detection of colorectal field carcinogenesis,” J. Biomed. Opt., 17 (4), 047005 (2012). Google Scholar


A. J. Radosevich et al., “Polarized enhanced backscattering spectroscopy for characterization of biological tissues at subdiffusion length scales,” IEEE J. Sel. Top. Quantum Electron., 18 (4), 1313 –1325 (2012). IJSQEN 1077-260X Google Scholar


S. C. Kanick et al., “Sub-diffusive scattering parameter maps recovered using wide-field high-frequency structured light imaging,” Biomed. Opt. Express, 5 (10), 3376 –3390 (2014). Google Scholar


W. Verkruysse et al., “Modelling light distributions of homogeneous versus discrete absorbers in light irradiated turbid media,” Phys. Med. Biol., 42 (1), 51 (1997). PHMBA7 0031-9155 Google Scholar


R. L. P. van Veen, W. Verkruysse and H. J. C. M. Sterenborg, “Diffuse-reflectance spectroscopy from 500 to 1060 nm by correction for inhomogeneously distributed absorbers,” Opt. Lett., 27 (4), 246 –248 (2002). OPLEDP 0146-9592 Google Scholar


J. C. Finlay and T. H. Foster, “Effect of pigment packaging on diffuse reflectance spectroscopy of samples containing red blood cells,” Opt. Lett., 29 (9), 965 –967 (2004). OPLEDP 0146-9592 Google Scholar


K. W. Calabro and I. J. Bigio, “Influence of the phase function in generalized diffuse reflectance models: review of current formalisms and novel observations,” J. Biomed. Opt., 19 (7), 075005 (2014). Google Scholar


S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol., 58 (11), R37 (2013). PHMBA7 0031-9155 Google Scholar


J. D. Rogers et al., “Modeling light scattering in tissue as continuous random media using a versatile refractive index correlation function,” IEEE J. Sel. Top. Quantum Electron., 20 (2), 173 –186 (2014). IJSQEN 1077-260X Google Scholar


M. Xu and R. R. Alfano, “Fractal mechanisms of light scattering in biological tissue and cells,” Opt. Lett., 30 (22), 3051 –3053 (2005). OPLEDP 0146-9592 Google Scholar


J. M. Schmitt and G. Kumar, “Turbulent nature of refractive-index variations in biological tissue,” Opt. Lett., 21 (16), 1310 –1312 (1996). OPLEDP 0146-9592 Google Scholar


İ. R. Çapoğlu et al., “Accuracy of the born approximation in calculating the scattering coefficient of biological continuous random media,” Opt. Lett., 34 (17), 2679 –2681 (2009). OPLEDP 0146-9592 Google Scholar


M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light, 7th ed.Cambridge University Press, Cambridge, England (1999). Google Scholar


A. Ishimaru, Wave Propagation and Scattering in Random Media. Volume I—Single Scattering and Transport Theory, 1st ed.Academic Press, New York (1978). Google Scholar


C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles, John Wiley and Sons, Hoboken, New Jersey (1983). Google Scholar


V. Turzhitsky et al., “A predictive model of backscattering at subdiffusion length scales,” Biomed. Opt. Express, 1 (3), 1034 –1046 (2010). BOEICL 2156-7085 Google Scholar


J. Lim et al., “Born approximation model for light scattering by red blood cells,” Biomed. Opt. Express, 2 (10), 2784 –2791 (2011). Google Scholar


J. D. Rogers, İ. R. Çapoğlu and V. Backman, “Nonscalar elastic light scattering from continuous random media in the born approximation,” Opt. Lett., 34 (12), 1891 –1893 (2009). OPLEDP 0146-9592 Google Scholar


A. J. Radosevich et al., “Ultrastructural alterations in field carcinogenesis measured by enhanced backscattering spectroscopy,” J. Biomed. Opt., 18 (9), 097002 (2013). Google Scholar


A. J. Radosevich et al., “Buccal spectral markers for lung cancer risk stratification,” PLoS One, 9 (10), e110157 (2014). POLNCL 1932-6203 Google Scholar


N. N. Mutyal et al., “In vivo risk analysis of pancreatic cancer through optical characterization of duodenal mucosa,” Pancreas, 44 (5), 735 –741 (2015). PANCE4 0885-3177 Google Scholar


A. J. Radosevich and J. D. Rogers, (2015) September ). 2015). Google Scholar


C. J. R. Sheppard, “Fractal model of light scattering in biological tissue and cells,” Opt. Lett., 32 (2), 142 –144 (2007). OPLEDP 0146-9592 Google Scholar


P. Guttorp and T. Gneiting, “Studies in the history of probability and statistics XLIX on the Matérn correlation family,” Biometrika, 93 (4), 989 –995 (2006). BIOKAX 0006-3444 Google Scholar


L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J., 93 70 –83 (1941). ASJOAB 0004-637X Google Scholar


J. R. Mourant et al., “Mechanisms of light scattering from biological cells relevant to noninvasive optical-tissue diagnostics,” Appl. Opt., 37 (16), 3586 –3593 (1998). Google Scholar


F. Bevilacqua and C. Depeursinge, “Monte Carlo study of diffuse reflectance at source-detector separations close to one transport mean free path,” J. Opt. Soc. Am. A, 16 (12), 2935 –2945 (1999). Google Scholar


A. J. Radosevich et al., “Structural length-scale sensitivities of reflectance measurements in continuous random media under the Born approximation,” Opt. Lett., 37 (24), 5220 –5222 (2012). OPLEDP 0146-9592 Google Scholar


J. Yi and V. Backman, “Imaging a full set of optical scattering properties of biological tissue by inverse spectroscopic optical coherence tomography,” Opt. Lett., 37 (21), 4443 –4445 (2012). OPLEDP 0146-9592 Google Scholar


A. Leslie et al., “Polymorphism of DNA double helices,” J. Mol. Biol., 143 (1), 49 –72 (1980). JMOBAK 0022-2836 Google Scholar


H. G. Hansma and J. H. Hoh, “Biomolecular imaging with the atomic force microscope,” Annu. Rev. Biophys. Biomol. Struct., 23 (1), 115 –140 (1994). ABBSE4 1056-8700 Google Scholar


A. J. Radosevich et al., “Open source software for electric field Monte Carlo simulation of coherent backscattering in biological media containing birefringence,” J. Biomed. Opt., 17 (11), 115001 (2012). Google Scholar


H. van de Hulst, Light Scattering by Small Particles, Dover, Mineola, New York (1981). Google Scholar


S. Prahl, “Tabulated molar extinction coefficient for hemoglobin in water,” (2012) December ). 2012). Google Scholar


D. Yudovsky and L. Pilon, “Rapid and accurate estimation of blood saturation, melanin content, and epidermis thickness from spectral diffuse reflectance,” Appl. Opt., 49 (10), 1707 –1719 (2010). APOPAI 0003-6935 Google Scholar


S. Chandrasekhar, Radiative Transfer, Courier Corporation, Mineola, New York (2013). Google Scholar


T. J. Farrell, M. S. Patterson and B. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys., 19 (4), 879 –888 (1992). MPHYA6 0094-2405 Google Scholar


A. Liemert and A. Kienle, “Exact and efficient solution of the radiative transport equation for the semi-infinite medium,” Sci. Rep., 3 (2013). SRCEC3 2045-2322 Google Scholar


Y. Kuga and A. Ishimaru, “Retroreflectance from a dense distribution of spherical particles,” J. Opt. Soc. Am. A, 1 (8), 831 –835 (1984). JOAOD6 0740-3232 Google Scholar


L. Tsang and A. Ishimaru, “Backscattering enhancement of random discrete scatterers,” J. Opt. Soc. Am. A, 1 (8), 836 –839 (1984). JOAOD6 0740-3232 Google Scholar


P. E. Wolf and G. Maret, “Weak localization and coherent backscattering of photons in disordered media,” Phys. Rev. Lett., 55 (24), 2696 –2699 (1985). PRLTAO 0031-9007 Google Scholar


M. P. V. Albada and A. Lagendijk, “Observation of weak localization of light in a random medium,” Phys. Rev. Lett., 55 (24), 2692 –2695 (1985). PRLTAO 0031-9007 Google Scholar


E. Akkermans, P. E. Wolf and R. Maynard, “Coherent backscattering of light by disordered media: analysis of the peak line shape,” Phys. Rev. Lett., 56 (14), 1471 –1474 (1986). PRLTAO 0031-9007 Google Scholar


A. J. Radosevich et al., “Depth-resolved measurement of mucosal microvascular blood content using low-coherence enhanced backscattering spectroscopy,” Biomed. Opt. Express, 1 (4), 1196 –1208 (2010). BOEICL 2156-7085 Google Scholar


A. J. Radosevich et al., “Measurement of the spatial backscattering impulse-response at short length scales with polarized enhanced backscattering,” Opt. Lett., 36 4737 –4739 (2011). OPLEDP 0146-9592 Google Scholar

Biographies for the authors are not available.

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Andrew J. Radosevich, Adam Eshein, The-Quyen Nguyen, and Vadim Backman "Subdiffusion reflectance spectroscopy to measure tissue ultrastructure and microvasculature: model and inverse algorithm," Journal of Biomedical Optics 20(9), 097002 (28 September 2015).
Published: 28 September 2015

Back to Top