Goniometric measurements of thick tissue using Monte Carlo simulations to obtain the single scattering anisotropy coefficient.

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.

1. Introduction

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 [3]. 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 (μ), scattering coefficient (μ) 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 [6]. Moreover, tissues need to be sufficiently thick (typically ~8-10 transport mean free paths) [7] 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. [8]. 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 [10] and the connective tissue disorder osteogenesis imperfecta [11]. 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. [16] 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 [17]. 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 [22] 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 [12] when tissues are greater than one scattering length, i.e. (1/μ), 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, μ, 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.

2. Experimental setup and Monte Carlo simulation model

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 [25], depending on tissue type, where cutting artifacts are larger for thinner samples. The scattering coefficient, μ, of most tissues at visible wavelengths are typically in the range of 200–500 cm−1 [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 μ, 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.

Fig. 1

(a) Experimental setup, (b) Flowchart for obtaining the single scattering anisotropy gsingle from measurements of geff and μ, (c) Geometry for Monte Carlo simulation.

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 [26]. 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 [9]:

Fig. 2

Validation of simulations. (a) Examples of angular distributions for different g values. (b) Snell’s law applied to correct for refracted angle. (c) Dependence of effective anisotropy on scattering length. (d) Wet cell geometry and validation of its minimal effect on the measured effective anisotropy.

The factor α accounts for losses due to refractive index mismatch (such as air-glass interfaces) and can be calculated via Fresnel reflection coefficients [9]. These are determined for complex geometries by Monte Carlo simulation [8] or measured without the tissue present. In our setup, we find this factor to be close to 93%. We calculate µ 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 µ 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 [8] 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 µ (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 μ, 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.

3. Validation of the simulation model

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

We now consider the accuracy of mapping values of (geff, µ) to gsingle by considering ranges of gsingle and corresponding geff values at select µ 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 µ 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 µ. Many tissues lie within this regime, therefore the method can be powerful for characterizing samples of several scattering lengths in thickness with good accuracy.

Table 1

Uncertainty relations between geff and gsingle for ranges of gsingle and μ

gsingle (range)	µsd	Δgsingle /Δgeff
0.00–0.30	0.5	10
0.30–0.50	0.5	2.5
0.50–0.70	0.5	1.3
0.50–0.70	4	5.6
0.70–0.90	0.5	1
0.70–0.90	4	1.2
0.70–0.90	8	2.6
0.90–0.995	0.5	0.9
0.90–0.995	4	0.3
0.90–0.995	8	0.2

We also note that the discrepancy between geff and gsingle even in the case of thin samples with 0.4 <µ <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

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.

4. Experimental validation over a range of mouse muscle tissue thickness

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 [26] and are listed in Table 2 . We note that these values are comparable to other data in the literature for muscle tissue [26] where the refractive index was found to be in the range 1.38–1.46 at 632.8 nm.

Table 2

Refractive index of the muscle tissue measured over a wavelength range 445–1000 nm

λ [nm]	n
445	1.40 ± 0.01
535	1.38 ± 0.02
680	1.38 ± 0.02
1000	1.37 ± 0.02

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 .

Fig. 3

Extraction of the single scattering anisotropy coefficient (gsingle) for mouse muscle tissue of thickness 200 microns.

Table 3

Single scattering anisotropy values extracted from different thickness slices of mouse muscle tissue

	Thickness
	100 microns				150 microns				200 microns
λ [nm]	µs [cm−1]	geff	gsingle		µs [cm−1]	geff	gsingle		µs [cm−1]	geff	gsingle
445	362 ± 50	0.84 ± 0.03	0.96 ± 0.02		310 ± 22	0.83 ± 0.02	0.96 ± 0.01		270 ± 11	0.79 ± 0.04	0.96 ± 0.01
680	327 ± 48	0.89 ± 0.02	0.972 ± 0.004		262 ± 25	0.87 ± 0.02	0.97 ± 0.01		233 ± 12	0.85 ± 0.04	0.97 ± 0.01
890	326 ± 47	0.92 ± 0.02	0.98 ± 0.01		255 ± 24	0.89 ± 0.04	0.98 ± 0.01		219 ± 10	0.91 ± 0.04	0.98 ± 0.01

We attribute the slight thickness dependence of µ 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 µ 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 µ to decrease with increased wavelength. Cheong et al. compiled data on muscle tissue [9], 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.

5. Validation using rat tail tendon

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 [27] 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 [28]. The wavelength dependency will be discussed in the next section.

Table 4

Optical properties of rat tail tendon of diameter ~170 microns

λ [nm]	n	µs [cm−1]	geff	gsingle
445	1.42 ± 0.01	450 ± 30	0.58 ± 0.09	0.91 ± 0.02
494		380 ± 30	0.71 ± 0.12	0.94 ± 0.02
535	1.40 ± 0.01	330 ± 30	0.77 ± 0.10	0.95 ± 0.01
680	1.40 ± 0.01	280 ± 20	0.73 ± 0.05	0.94 ± 0.01
890	1.38 ± 0.01	250 ± 20	0.80 ± 0.05	0.95 ± 0.01
988	1.39 ± 0.01	230 ± 20	0.80 ± 0.03	0.95 ± 0.01
1070	1.37 ± 0.01	220 ± 20	0.78 ± 0.03	0.94 ± 0.01

6. Wavelength dependency of the reduced scattering coefficient

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.

Fig. 4

The wavelength dependency of the reduced scattering coefficient for mouse muscle and rat tail tendon fitted to a power law relation.

7. Simple deconvolution-based model for approximating the single scattering anisotropy

It has been shown analytically [35] 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 µ. Table 5 lists results where simulated data was used to evaluate the error of the method for a range of gsingle and µ. 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 µ 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 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 µ to accurately retrieve gsingle from knowledge of geff.

Table 5

Evaluation of deconvolution-based method based approach for obtaining single via Eq. (6)

gsingle (Used in simulation)	µsd (Range)	Δgsingle (Due to deconv. method)
0.40	0.20–1.00	0.04–0.06
0.40	1.00–2.00	0.06–0.28
0.70	0.20–1.00	0.02–0.06
0.70	2.00–4.00	0.07–0.14
0.80	0.20–1.00	0.01–0.07
0.80	2.00–4.00	0.01–0.07
0.80	4.00–10.00	0.07–0.13
0.90	0.20–1.00	0.00–0.07
0.90	2.00–4.00	0.02–0.01
0.90	4.00–10.00	0.01–0.04
0.95	0.20–1.00	0.02–0.04
0.95	1.00–2.00	0.04–0.02
0.95	2.00–10.00	0.02–0.01
0.98	0.20–10.00	< 0.01

8. Conclusions

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

Table 6

Effect of absorption on effective anisotropy for a range of g values for two scattering lengths, µ = 1 and 4

	geff
	g = 0.20			g = 0.40			g = 0.60			g = 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
0	0.421	0.463		0.459	0.468		0.561	0.488		0.728	0.565		0.834	0.683
5	0.426	0.466		0.472	0.472		0.580	0.495		0.747	0.578		0.851	0.698
10	0.431	0.468		0.481	0.476		0.592	0.502		0.757	0.589		0.859	0.708

Simulation parameters: n = 1.40, d = 0.0100 cm (with wet cell geometry)