|Home | About | Journals | Submit | Contact Us | Français|
The scattering anisotropy, g, of tissue can be a powerful metric of tissue structure, and is most directly measured via goniometry and fitting to the Henyey-Greenstein phase function. We present a method based on an independent attenuation measurement of the scattering coefficient along with Monte Carlo simulations to account for multiple scattering, allowing the accurate determination of measurement of g for tissues of thickness within the quasi-ballistic regime. Simulations incorporating the experimental geometry and bulk optical properties show that significant errors occur in extraction of g values, even for tissues of thickness less than one scattering length without modeling corrections. Experimental validation is provided by determination of g in mouse muscle tissues and it is shown that the obtained values are independent of thickness. In addition we present a simple deconvolution-based method and show that it provides excellent estimates for high anisotropy values (above 0.95) when coupled with an independent attenuation measurement.
Quantitative techniques in tissue and biomedical optics have emerged as important tools with many applications in biological research and clinical studies [1,2]. A critical aspect common to many optical methods is determination of the light distribution in tissue upon illumination. For example, this characterization underlies technologies and modalities used for disease diagnostics such as diffuse optical tomography . A commonly used approach uses the radiative transfer equation (RTE) [4,5], through which the principle of energy conservation describes how a beam of light traverses through media characterized by four bulk optical properties: refractive index (n), absorption coefficient (μa), scattering coefficient (μs) and scattering anisotropy (g). In general, approximations are required to analytically and/or numerically solve the RTE for a given sample geometry, where the diffusion approximation is often used . Moreover, tissues need to be sufficiently thick (typically ~8-10 transport mean free paths)  to fall within this regime. Easier and sometimes more accurate solutions are obtained by using Monte Carlo (MC) simulations that track photons through the sample geometry. An additional advantage of the MC method is that tissue can be composed of multiple layers, each having different properties, and arbitrary geometries, as implemented in the Monte Carlo Multi-Layered (MCML) implementation of Wang et al. . In addition, the MC approach is applicable to the quasi-ballistic regime typical for optical microscopy applications.
For successful implementation, both diffuse and MC methods require accurate values of the bulk optical properties. Many methods have been devised for characterizing tissues in terms of the four coefficients and reviews of techniques along with tables of optical properties for various tissue types have been published [1,9]. It is increasingly understood that the angular distribution, , and the anisotropy,, of a single scattering event are effective metrics for characterizing tissue structure. For example, Campagnola and associates showed that the scattering anisotropy was different in diseased tissues relative to normal in both ovarian cancer  and the connective tissue disorder osteogenesis imperfecta . Several methods have been employed to describe the effective angular distribution of multiply scattered light, where goniometry is the most common approach [12–15]. In addition, Chaikovskaya et al.  used an enhanced diffusion approximation of the RTE to predict the angular distribution of scattering and obtained g from very thick tissues (~3mm). In the present work we focus on the quasi-ballistic regime where diffusion approximations generally fail . Other less direct methods for obtaining g such as reflectance confocal microscopy [18,19], low coherence interferometry [20,21] and phase microscopy of thin slices  have been implemented and shown promise for clinical applications.
It has been shown that for tissues the angular distribution of scattering is well described by the Henyey-Greenstein (HG) phase function, which is fully parameterized by the anisotropy g. More recent work, based on deriving optical properties through a refractive index correlation function of the more general Whittle-Matérn family, has shown that the HG phase function is a special case which is applicable for a wide range of biological tissues [20,23,24]. While goniometric data has often been fit to the HG function, a limiting factor in this approach is that sample thickness affects the angular distribution and the fitted scattering anisotropy  when tissues are greater than one scattering length, i.e. (1/μs), in thickness. In this case multiple scattering leads to an apparent lower g value.
In this paper we develop and apply a novel method to accurately determine g over a range of tissue thickness. Rather than exploring new or modified phase functions we assume that the scattered data fits well to the HG function even in the case of multiple scattering but with an effective anisotropy coefficient (validated in Appendix B). The method is based on the combination of an on-axis attenuation measurement for obtaining the scattering coefficient, μs, via the Beer-Lambert law and measuring the angular distribution of the multiply scattered light by goniometry. MC simulations are then used to extract the single scattering anisotropy, g. To validate the approach, we present data and analysis on mouse muscle tissue of different thicknesses and rat tail tendon at excitation wavelengths ranging from 445–890 nm. Lastly we evaluate the results from a simple deconvolution-based method coupled with the attenuation measurement.
Specimens are sliced with a Vibratome (Leica), where either fresh or fixed tissues can be used. The minimal thicknesses are generally in the range of 50-100 microns , depending on tissue type, where cutting artifacts are larger for thinner samples. The scattering coefficient, μs, of most tissues at visible wavelengths are typically in the range of 200–500 cm−1 , and this range corresponds to a few scattering lengths in typical thin sample thicknesses used in microscopy experiments. We shall show that relatively large errors in g can be obtained even for relatively thin tissues. In addition, our analysis accounts for effects such as refractive index, sample geometry and absorption that can cause errors in the obtained g values. We note that for terminology, the single scattering anisotropy will be denoted by gsingle, where effective anisotropy, geff, will be used for anisotropy coefficients obtained from fitting angular scattering data (both simulated and experimental data) to the HG phase function. We also use μsd, to denote the average number of scattering lengths in the sample, where d is the physical tissue thickness.
Figure 1 (a) illustrates the experimental setup. The sample is illuminated with a collimated laser beam and the angular distribution of the scattered component is measured using a photodiode (Thorlabs). The laser source was a femtosecond mode-locked Coherent Chameleon Ultra Ti:sapphire laser tunable from 680–1070 nm, and an externally frequency doubled component therein, yielding a full tunability of about 340–1070 nm. To optimize the signal-to-noise ratio (SNR), the beam was chopped and the reference signal and the photodiode detector were connected to a lock-in amplifier (Signal Recovery). The signal readout from the lock-in was sent to a computer that also synchronizes the rotation of the detector using a motorized rotation stage to measure the scattered intensity as a function of angle.
The flowchart for the full process of obtaining measurements of gsingle is shown in Fig. 1(b). The tissue sample is extracted and mounted on a glass slide in a wet cell (depicted in Fig. 2(d) ). An estimate of refractive index is determined via a total internal reflection measurement of a sample using the method of Li . The second step involves measurement of the on-axis attenuation of the sample, where the goniometer is aligned to measure the transmission with and without the sample. The transmission follows the Beer-Lambert law :
The factor α accounts for losses due to refractive index mismatch (such as air-glass interfaces) and can be calculated via Fresnel reflection coefficients . These are determined for complex geometries by Monte Carlo simulation  or measured without the tissue present. In our setup, we find this factor to be close to 93%. We calculate µs by assuming µa << µs as follows:
The angular distribution of multiple scattering is measured in the forward direction over a 180° range. Here we correct for refraction using Snell’s law and the measured refractive index. The data is fit over the range of 15°–40°to the HG phase function, which is of the form
The upper end is limited as by the critical angle defined by the refractive index, which for typical tissues is within 40–50 degrees, which corresponds to refractive indices in the range 1.31–1.56. The lower limit of 15 degrees is used to avoid the intense, ballistic component of the laser. The fit is performed on both sides of the curve (over 15–40°) and obtains two independent measurements of geff.
These values of geff are used as inputs, along with µsd to Monte Carlo simulations to extract gsingle by determining the best fit to the experimental data. The simulation geometry is illustrated in Fig. 1(c). The major difference from the experimental configuration is that in the MC simulation light is detected on a large planar screen. The size of the screen and the distance D can be controlled in the simulation and controls the angular sampling size and maximum angle. To accurately sample the angle it is important to account for the projection effect that increases with larger angles, which is described by where in solid angle, which needs to be divided out for appropriate sampling. We use the open source MCML framework  to perform the simulations and the fitting and simulation data analysis are run in MATLAB (MathWorks). The simulations are typically performed over a large range of parameters for g (0.00 to 0.995 in 100 steps) and µs (1 to 1000 in 100 steps) with 107 photons for each run. The data is fitted to Eq. (3) to obtain geff and for each simulation, gsingle, geff and μsd, are tabulated. The resulting grid is then stored and used for subsequent fitting to the experimental data to extract gsingle. The simulations were run on the Condor deployed high-throughput computing cluster at the Center for High-Throughput Computing (CHTC) at the University of Wisconsin-Madison. The simulation framework is available upon request.
We provide a validation of the simulation model and then systematically investigate the effect of sample refraction, geometry and absorption on extracted g values. Figure 2(a) shows examples of simulated data for a few different gsingle values for a thin (1 micron) sample with no refraction (n = 1.0), no absorption, and a scattering coefficient of 100 cm1. The angular data results from a single scattering event and geff obtained from fitting matches that of gsingle very closely, providing initial validation of the simulation as for very thin tissue we expect gsingle and geff to be very close in the limit of only single scattering event. Figure 2(b) shows an example of the effect of sample refraction on the angular distribution with and without angular correction via Snell’s law. For this case, the geff decreases to ~0.73 without correction, whereas the proper value of 0.80 is recovered when the data is properly corrected for refraction. The simulation incorporates internal reflection within the sample (which cannot be done analytically) and thus this shows that the effects of this are minimal and the Snell’s law correction provides good results.
Figure 2(c) illustrates the effect of multiple scattering on the effective anisotropy for a range of gsingle values. As evident from the graph, pairs of gsingle and µsd map into unique locations of geff for a range of values. This forms the basis of our method for extracting gsingle from independent experimental measurements of geff and µsd.
We now consider the accuracy of mapping values of (geff, µsd) to gsingle by considering ranges of gsingle and corresponding geff values at select µsd locations. Here we define Δgsingle and Δgeff such that a value of gsingle in the range of Δgsingle corresponds to a value of geff with a range Δgeff following simulations. By quantifying the corresponding ratios we obtain an estimate of errors over a range of parameter space. Table 1 lists the ratios for a range of cases at specified µsd locations. We highlight that for low anisotropies 0 < gsingle < 0.3, the accuracy is low, with the experimental uncertainty typically being magnified 10-fold. On the other end for high anisotropies, 0.90 < gsingle < 0.995, the uncertainty of geff diminishes and decreases with increasing µsd. Many tissues lie within this regime, therefore the method can be powerful for characterizing samples of several scattering lengths in thickness with good accuracy.
We also note that the discrepancy between geff and gsingle even in the case of thin samples with 0.4 <µsd <1.2 is significant. As an example, for gsingle = 0.8 the variation increases from 0.03 to 0.09, and for gsingle = 0.90 it increases from 0.03 to 0.08. This highlights the need for modeling even in the case of relatively thin samples. It is also worth noting the seemingly artifactual behavior for gsingle values in the range 0.00–0.30, where gsingle<geff. This arises because of the model slab geometry, which inherently filters higher angles as these have longer paths to exit the tissue. However, most tissues have anisotropies in the range 0.60–0.98 and this artifact is not practically significant.
Figure 2(d) shows the mounting geometry we refer to as a wet cell. This mounting approach is necessary as measured values of g will increase as tissues dehydrate. To keep the samples hydrated they are mounted on a standard 1 mm glass slide with phosphate buffered saline (PBS) below and on top of the tissue section. A coverglass is placed on top of the sample and sealed with either nail polish or Vaseline. Tissues then stay hydrated for up to a few hours. We investigated the effect of adding the glass cell on geff via simulations as listed in the table on Fig. 2(d). We note that the effect is minimal for tissues thicker than 10 microns.
If our modeling approach is valid, the recovered values of gsingle should be independent of the tissue thickness. In order to validate the method experimentally, fixed mouse (Col1a1 model; Jackson labs) muscle tissues were sliced to three different thicknesses (100, 150 and 200 microns). We used three sections of each thickness and made the bulk optical properties measurements in three locations at three wavelengths: 445 nm, 680 nm and 890 nm. The wavelength dependency provides additional validation of the model, and, as described in section 5, provides information about the tissue structure.
The refractive indices of the muscle tissue samples were measured by the total internal reflection method  and are listed in Table 2 . We note that these values are comparable to other data in the literature for muscle tissue  where the refractive index was found to be in the range 1.38–1.46 at 632.8 nm.
Figure 3 shows plots of the data obtained for 200 micron thick slices of muscle. The scattering coefficient is shown in (a), the effective anisotropy in (b) and the resulting single scattering anisotropy (c) as obtained via best fit of the experimental data to simulations. The data for all three thicknesses is summarized in Table 3 .
We attribute the slight thickness dependence of µs to thickness variations from slicing. We also note that this error does not influence the measurement of gsingle as it is actually the value of µsd that is measured directly with the on-axis attenuation measurements via the Beer-Lambert law. As can be seen from the table, the agreement of gsingle between all three thicknesses is excellent and is well within experimental error, providing a good experimental validation of the technique. We also note the expected trends of both geff and gsingle to increase with wavelength and µs to decrease with increased wavelength. Cheong et al. compiled data on muscle tissue , where the g values for porcine, chicken and bovine tissue were listed to be in the range 0.94–0.97 (with one exception of 0.30). We note that it is difficult to get exact correspondence to that of other studies due to different tissue sources and their respective preparations. However, given our extracted values for the single scattering anisotropy are all in a similar range and in good agreement with the literature, we have confidence in applying our approach to other tissue types.
As a second example, we applied this approach to determine gsingle in the highly ordered rat tail tendon. The tissue was prepared by extracting tendon from Harlan Sprague Dawley rats and subsequent fixation for 1 hour using 4% formalin. The bulk optical properties were measured as before. The thickness of the tissue was obtained by adding a dilute amount of fluorescent beads on both the glass surfaces (such that the beads adhere to the glass), where 3D confocal fluorescence microscopy  provided an estimate of ~170 microns. The experimental data and extracted geff are listed in Table 4 . As in the case of muscle, the simulations extracted consistent values of geff with those of the literature, where g~0.97 has been reported at 650 nm . The wavelength dependency will be discussed in the next section.
It has been shown empirically that the reduced scattering coefficient has a wavelength dependency in the form of a power law: [29,30]. More recently a theoretical foundation has been established for this relationship [23,31,32]. In a fairly general treatment, it has been shown that the scattering power b has a specific relationship with the refractive index correlation function via , where the parameter m defines the shape of the refractive index correlation function and therefore has a direct relationship with the tissue structure [20,23,33].
Figure 4 shows the wavelength dependency of the reduced scattering coefficient for the mouse muscle and rat tail tendon data from the previous sections. The data was fit to a power law and the scattering power was obtained as b = −1.30 ± 0.12 and b = −1.06 ± 0.05 for the muscle and tendon respectively. That is, the reduced scattering of muscle is slightly more dependent on the wavelength than the tendon. In terms of the structure parameter m, the fits yield m = 1.47 ± 0.03 and m = 1.36 ± 0.06 respectively. The range of corresponds to a mass fractal medium, where this regime is characterized by scatterers that are self-similar over a characteristic size scale, rather than similarity in scatterer size. As has been noted in the literature, many tissues fall into this regime [31,34]. Therefore the wavelength behavior of the tissues investigated is quite typical for soft tissues and provides further validation of our method.
It has been shown analytically  and validated computationally [15,36] that, in the case of forward peaked scattering, the effective anisotropy following a convolution of a discrete number N of scattering events can be given by
However, for real tissues, there will be distribution in the number of scattering events before photons exit. Additionally, in goniometry, larger scattering angles have a larger component of multiply scattered photons. It is possible to estimate the average number of scattering events as detected in goniometry (always at least 1 scattering event) by
In order to obtain gsingle from geff and Navg, we invert Eq. (4) and obtain
Here we evaluate Eq. (6) as a technique for obtaining the single scattering anisotropy directly from independent measurements of geff and µsd. Table 5 lists results where simulated data was used to evaluate the error of the method for a range of gsingle and µsd. As can be seen the error is high for gsingle lower than 0.70, in which case the Monte Carlo based method should be used. However, for gsingle above 0.95 the deconvolution-based method generally produces excellent results. This is presumably because for high anisotropy values the scattering is more densely packed around the forward peak and there is less distribution in µsd for different photons and the approximation in Eq. (6) is more valid than for broader angular distributions. Applying this simpler method to the data from previous sections (Tables 3 and and4),4), we found excellent agreement for gsingle with deviations of less than 0.013. This highlights the usefulness of the deconvolution in the case of very highly forward-scattering tissues (gsingle above 0.95) as were used in this study. However, most tissues have somewhat lower anisotropies, thus in the general case utilizing the simulation approach is more accurate and provides a better picture of the uncertainty in the obtained anisotropy values. Lastly we emphasize that both methods still require an independent measurement of µsd to accurately retrieve gsingle from knowledge of geff.
We have presented a method that obtains the single scattering anisotropy from tissues of up to a several scattering lengths in thickness by combining the measurement of the angular distribution with an independent measurement of the scattering coefficient and then using Monte Carlo simulations to provide the best fit to the experimental data. The method has been validated both by simulations over a broad range of model parameters and also experimentally through measurements of a thickness series of muscle tissue. The largest advantage of this method is that it can be used for samples of a several scattering lengths of thickness which fall within the quasi-ballistic regime found in microscopy applications. The method is also important for thin samples as even for samples of less than one scattering length there can be relatively large errors in the g obtained from fitting as well as errors induced refractive index mismatch. Furthermore we have demonstrated that the method is relatively insensitive to slight errors in the absorption and refractive index used in the simulations as long as they are within typical ranges found in tissue (Appendix A). The MC approach is broadly applicable to tissues, and we also showed that the attenuation measurement coupled with a simple deconvolution based approach can be applied to highly anisotropic (g>0.95) tissues. Future work may consider utilizing other phase functions such as those derived based on Mie theory and the Whittle-Matérn family of refractive index correlation functions that are more analytical in nature and can provide additional information about the underlying tissue structure
We gratefully acknowledge support from National Institutes of Health (NIH) National Institute of Biomedical Imaging and Bioengineering (NIBIB) R01-EB000184, National Cancer Institute R01 CA136590-01A1, National Science Foundation (NSF) Chemical, Bioengineering, Environmental, and Transport Systems (CBET) 0959525 and the Wisconsin Institutes for Discovery. G. H. acknowledges support from a Leifur Eiríksson Foundation Scholarship. We thank Joseph Szulczewski and the Keely laboratory for supplying and sectioning the mouse muscle tissues used in this study. We thank the Research Animal Resources Center (RARC) of University of Wisconsin—Madison for providing the rat tissues used in this study. We acknowledge the Center for High Throughput Computing (CHTC) at the University of Wisconsin—Madison and the Open Science Grid (OSG) for use of their computer cluster for performing simulations.
We evaluate the effect of variations in refractive index and absorption on the effective anisotropy. This is important for evaluating whether accurate measurements of n and µa are needed or if approximations are sufficient. Values of µa of tissues are typically in the range 0–10 cm−1 in the visible and extending to near infrared wavelengths . This may be higher for tissues with blood, in which case an estimate of µa is desirable. Refractive indices are generally within 1.35–1.50 in this wavelength range [1,26]. Figure 5
illustrates the effect of absorption on the angular distribution for three cases with typical values of g, µs and n. The effective anisotropy increases with increased absorption, as the larger angles have a longer path to exiting the tissue (slab geometry) and are filtered to a larger degree, making the beam more forward peaked, thereby increasing the obtained geff. The data shown in Table 6
0.80||g = 0.90|
|µa||µs = 100||µs = 400||µs = 100||µs = 400||µs = 100||µs = 400||µs = 100||µs = 400||µs = 100||µs = 400|
aSimulation parameters: n = 1.40, d = 0.0100 cm (with wet cell geometry)
document the effect for a range of gsingle and µs. In summary, for each row, we note that the effect of increasing µa from 0 to 10 cm−1 is contained within about 0.03 in Δgeff. Thus to obtain higher accuracy, it is important to include a measure of µa in the simulations. However, in many situations an error would lead to, which may be comparable to the experimental error.
We also evaluated the effect of introducing a slight error in the refractive index of the sample. For demonstration, we evaluated the effect of using n=1.40 in the analytical Snell law correction of the refracted angle, where this was between 1.35 and 1.45 over a range of μsd from 0.5 to 10 and g between 0.00 and 0.95. We found that this error had a very minimal effect on geff such that Δgeff was less than 0.01 in all cases (data not shown).
The method presented in this paper relies on fitting both experimental and simulated data to the HG phase function as defined in Eq. (3). As the function is generally intended to be applied to single scattered data, it is important to establish that upon multiple scattering that the shape of the HG function is preserved but with an effective anisotropy coefficient in place of the single scattering anisotropy.
As validation, Fig. 6
shows the effect on the simulated phase functions for increasing µsd for a case where gsingle=0.90. The effective phase functions are fit to the HG function and the obtained geff. In addition, we calculate the R2 values to quantify the goodness of fit between the effective phase function and the HG function. As can be seen an excellent fit is obtained in the range where µsd <5 (typical range for the quasi-ballistic regime) with R2 above 0.99. For µsd higher than 10, the effective phase function changes shape and the HG function does not fit quite as well (R2=0.97) but is still reasonable. We note that the curve saturates for µsd higher than 20, e.g. for µsd=30, the simulated data falls exactly on the same curve as for µsd=20. This indicates that the diffuse regime has been reached at that thickness. Furthermore, we have noticed that in this regime the effective phase function generally is invariant with increasing µsd in the forward direction. While we have only considered the forward direction in this study, the results of Chaikovskaya et al.  suggest that by including the backward direction (reflected portion) of the phase function, it is possible to obtain gsingle for thick tissues.