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
(

*μ*_{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 [

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/

*μ*_{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.

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,

*μ*_{s},
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

*g*_{single}, where effective anisotropy,

*g*_{eff}, 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

*μ*_{s}d, to denote the average number of
scattering lengths in the sample, where

*d* is the physical tissue
thickness.

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

*g*_{single} is shown in . The tissue sample is extracted and mounted on a glass slide in a wet cell
(depicted in
). 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]:

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

*µ*_{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
*g*_{eff}.

These values of

*g*_{eff} are used as inputs, along with

*µ*_{s}d to Monte Carlo simulations to extract

*g*_{single} by determining the best fit to the experimental data. The
simulation geometry is illustrated in . 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

*µ*_{s} (1
to 1000 in 100 steps) with 10

^{7} photons for each run. The data is fitted to

Eq. (3) to obtain

*g*_{eff} and
for each simulation,

*g*_{single},

*g*_{eff} and

*μ*_{s}d, are tabulated. The resulting grid is then stored and
used for subsequent fitting to the experimental data to extract

*g*_{single}. 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.