^{′}

_{s}) and absorption (μ

_{a}). Experimental work supports our findings and establishes the advantages of Monte Carlo–based evaluation.

## 1.

## Introduction

Many applications within the field of biomedical optics rely either on the capability of performing, or the availability of, accurate measurements of optical properties of highly scattering materials. This is reflected by the massive development of theory related to light propagation in turbid media (photon migration), and the numerous techniques available for characterisation of such materials.

Time-resolved spectroscopy (TRS) is one of several techniques available for assessing optical properties of turbid media. By studying the broadening of short (picosecond) light pulses, it allows determination of both the absorption coefficient
$\left({\mu}_{a}\right)$
and the reduced scattering coefficient
$\left({\mu}_{s}^{\prime}\right)$
without the need of absolute measurements of light intensities. At first, time-resolved photon migration was investigated and evaluated using the diffusion approximation of radiative transport theory.^{1} Although proven useful in numerous cases, the diffusion modeling was found erroneous at low albedos or close to radiation sources.^{2, 3} Unfortunately, several tissue types exhibit optical properties in the range where the use of the diffusion approximation may be questioned. The human prostate is one example,^{4, 5} exhibiting reduced scattering below
$10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
and absorption above
$0.3\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
.

In order to extend the range of optical properties and source-detector separations over which TRS can provide accurate data, methods for Monte Carlo-based modeling were developed. Introduced to the field of biomedical optics and photon migration by Wilson and Adam,^{6} Monte Carlo (MC) simulation has become the gold standard for modeling of light propagation in tissue optics.^{7} Besides modeling spatially and temporally resolved light distribution, MC has proven useful in, for example, fluorescence modeling^{8} and Raman spectroscopy.^{9}

Traditionally, MC was performed for a particular set of optical properties at a time. Since the simulation is time-consuming, it is thus not useful as a forward solver in reverese problems, for instance, during evaluation of experimental data (iterative curve-fitting). This obstacle lead to the development of White Monte Carlo (WMC),^{10, 11, 12} in which a single simulation in combination with proper rescaling ensures coverage of a wide range of optical properties. The main feature in WMC, making it feasible for use as a forward model in an iterative solver, is illustrated in Fig. 1
.

During MC simulations of light propagation in homogeneously scattering and nonabsorbing media, photon paths are determined only by the scattering phase-function, the scattering coefficient, and the sequence of random numbers generated by the simulation program. For a given phase function, and considering a particular sequence of random numbers, the photon path will scale linearly with the scattering coefficient, ${\mu}_{s}$ . The recording of photon paths and corresponding time-of-flights thus allows post-simulation transformation, from the scattering coefficient used during simulation, to an arbitrary. Furthermore, the impact of nonzero absorption can be imposed post-simulation simply by giving photons different weights ${w}_{a}$ according to the Beer-Lambert law of attenuation. This weight is stated in Eq. 1, where ${c}^{\prime}$ is the speed of light within the media:

This means that a single Monte Carlo simulation can be used to extract photon time-of-flight distribution not only for different source-detector separations, but also for different optical properties ${\mu}_{s}^{\prime}$ and ${\mu}_{a}$ .The preceding theory has been discussed in several papers.
^{8, 10, 11, 12, 13, 14, 15} Graaff suggested a limited scalable Monte Carlo technique where the optical properties of simulations in slab geometries could be scaled as long as the total attenuation coefficient was held constant.^{13} Two groups simultaneously and independently extended the theory of Graaff Kienle and Patterson suggested a scalable Monte Carlo technique for infinite and semi-infinite homogeneously scattering media and used independent (reference) MC simulations for verification of its performance in the semi-infinite case.^{10} Pifferi suggested a similar approach,^{11, 12} proposing scalability in both
${\mu}_{a}$
and
${\mu}_{s}^{\prime}$
. In practice, however, they based their evaluation on interpolation between MC simulations carried out at different
${\mu}_{s}$
, thus utilizing scalability in
${\mu}_{a}$
only. On the other hand, their work included evaluation of experimental time-resolved transmittance in the slab geometry (a geometry not allowing
${\mu}_{s}$
rescaling). Swartling showed the usefulness of the WMC approach in fluorescence emission spectra modeling.^{8} Swartling also raised the important question regarding the equivalency of WMC and traditional MC.^{14} Xu demonstrated the superior performance of WMC over different light propagation models as a forward model for evaluation of frequency domain data (generated by a traditional Monte Carlo program^{15}). In the same paper, Xu demonstrated scaling of data from an absorbing media using the weighting relationships from traditional Monte Carlo. Xu also reported the so-called *scaling effect* as being the major inconvenience in WMC-based data evaluation. This inconvenience originates from the fact that scaling in
${\mu}_{s}$
is accompanied by a scaling of temporal and spatial bin size, resulting in a need for data resampling and a limited range in scalability.

Motivated by an interest in *in vivo* optical characterization of human prostate tissue, this work is aimed at providing a scheme for fully scalable WMC for time-domain photon migration and demonstrating its value in evaluation of experimental data in the low albedo regime of photon migration. The approach is useful in both infinite and semi-infinite geometries, featuring individual storage of the spatial location and the time-of-flight for all potential detection events. Since this allows post-scaling binning (temporal, as well as spatial), it eliminates the need for the data resampling otherwise accompanying
${\mu}_{s}$
scaling. It also allows an accurate account for finite extension of source and detector areas (e.g., optical fiber diameters).

## 2.

## Materials and Methods

## 2.1.

### White Monte Carlo

A WMC model was developed to serve as the forward model for evaluation of fiber-based time-resolved spectroscopy under interstitial as well as under noninvasive conditions (i.e., infinite and semi-infinite geometry). The model consists of a simulation program written in C and a set of MATLAB scripts performing post-simulation processing. The main objectives were to retain full scalability in both ${\mu}_{s}^{\prime}$ and ${\mu}_{a}$ , while avoiding scaling inconveniences by moving spatial and temporal binning to post-simulation. This eliminates the need for any temporal resampling and allows accurate convolutions to properly account for the finite size (radius ${R}_{f}$ ) and light distribution of optical fibers. A schematic illustration of the WMC model is given in Fig. 2 . Apart from selecting simulation geometry, the user has to provide the WMC simulation program with a few input parameters. The refractive index $n$ , the anisotropy factor $g$ , and the numerical aperture $NA$ of the involved optical fibers are material parameters specific to a particular simulation. The simulation is run at the specified scattering level ${\mu}_{s}^{\mathit{max}}$ but is always terminated when time reaches ${t}_{\mathit{max}}$ . This notation is used because ${\mu}_{s}^{\mathit{max}}$ defines the maximum scattering coefficient for which the resulting database can provide valid data throughout the time interval $[0,{t}_{\mathit{max}}]$ . Generally, scaling to a certain scattering ${\mu}_{s}={\mu}_{s}^{\mathit{max}}\u2215\alpha $ results in data valid in the time interval $[0,\alpha {t}_{\mathit{max}}]$ .

## 2.1.1.

#### Simulation program

The WMC simulation program was written in ANSI C, adapting some parts from the *de facto* standard MC simulation program multi-layer Monte Carlo (MCML^{7}). As the photon propagation is similar to traditional MC, we refer to the work of Wang
^{7} and Prahl
^{16} for basic Monte Carlo theory. Here, we present features differing from their work.

An important change is the adaptation of a state-of-the-art pseudo-random number generator, the Mersenne twister by Saito and Matsumoto.^{17} The implementation used was the double-precision SIMD-oriented fast Mersenne twister (dSIMD) version 1.2.1, featuring a
${2}^{132049}-1$
period, documented excellent equidistribution properties, and fast random number generation. The entire
$32\text{-}\text{bit}$
output of the *time ()* function in C was used to seed the generator.

Photons are launched from the origin of a Cartesian coordinate system (as in MCML). The launch direction is in the positive $z$ direction (downward) with the addition of a deflection angle, $\theta $ , representing the emission cone of the source fiber, defined by the NA of the fiber, ${\theta}_{\mathit{max}}=\mathrm{arcsin}(NA\u2215n)$ . The angular distribution is assumed flat, i.e., $\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\theta =1-\xi (1-\mathrm{cos}\phantom{\rule{0.2em}{0ex}}{\theta}_{\mathit{max}})$ , where $\xi $ is a random number, $\xi \u03f5[0,1)$ . Although unnecessary in a cylindrically symmetric problem, the azimuthal angle $\phi $ is also randomized, $\phi =2\pi \xi $ , where $\xi \u03f5[0,1)$ . Note that taking the angular distribution of the source fiber into account prevents post-simulation scaling of the refractive index. Thus, different WMC simulations are required in order to handle different fiber $NA$ , as well as different refractive indices.

The step size, $\Delta s$ , is calculated using the scattering coefficient instead of the total attenuation coefficient. It is defined in Eq. 2:

where $\xi $ is a random variable, $\xi \u03f5(0,1]$ , implying $\Delta s\u03f5[0,\infty )$ . As in MCML, the actual photon scattering is simulated using the Henyey-Greenstein phase function.The photon detection scheme assumes that the source and detector optical fibers are parallel, having their tips in the same plane (a plane to which the fibers are also assumed to be orthogonal). This corresponds to the settings followed in the interstitial clinical measurements on human prostate tissue presented in Ref. 5, where optical fibers are inserted to the same depth. Regardless of whether simulations are performed in the infinite or semi-infinite geometry, a potential photon detection event requires that a photon path crosses the source-detector plane.

In the case of an infinite medium, photons passing upward through the source-detector plane, being within the acceptance cone of the detector fiber, may be detected. Such a crossing is referred to as a *detection event* and is always registered by storing its radial distance from the source
$r$
, as well as its total time-of-flight
$t$
. For memory conservation, storing is done using
$32\text{-}\text{bit}$
floating point variables. Since the photon is propagating in an infinite medium, and since the position of the detector fiber is undefined, the photon is, however, not terminated at this point. Instead, photon termination occurs only when the time reaches
${t}_{\mathit{max}}$
. Note that this implies that a single photon may generate multiple detection events. These events are considered independent by the proposed WMC scheme, inducing a small error in generated time-of-flight histograms. This issue is further discussed in Sec.
2.1.3.

For the case of a semi-infinite geometry, the source-detector plane equals the medium boundary. Photons escaping the medium at an angle within the acceptance cone of the detector fiber generate detection events, i.e., $r$ and $t$ are recorded. Here, a notable difference from traditional Monte Carlo in semi-infinite geometries is that partial reflections are not allowed by WMC (the weight of the photons are not monitored). Instead, the reflection coefficient is compared to a random number, and the photons are either reflected or transmitted and terminated. Note that this also means that the semi-infinite detection geometry does not suffer from the risk of double-detecting photons. Note also that photons are terminated not only upon boundary crossing, but also when the time exceeds the predefined maximum time-of-flight ${t}_{\mathit{max}}$ .

## 2.1.2.

#### Post-simulation processing

The database of detection events (
$r$
,
$t$
pairs) from the simulation are sorted with respect to
$r$
. This is done using the fast {
$O\left[n\phantom{\rule{0.2em}{0ex}}\mathrm{log}\left(n\right)\right]$
worst-case time complexity} in-place sorting algorithm Combsort11.^{18} The
$O\left(1\right)$
memory usage of the Combsort11 algorithm enables sorting and usage of databases of sizes close to the physical memory of the computer used.

The sorted database is used by a MATLAB function generating photon time-of-flight distributions. This function is employed as the forward model in an iterative solver and is called with the following input parameters:
${\mu}_{a}$
and
${\mu}_{s}^{\prime}$
being the optical properties corresponding to the resulting time-of-flight distribution, the temporal channel width
$\Delta t$
, source-detector separation
$\rho $
, and the fiber radius
${R}_{f}$
(experiments typically involve two identical fibers). The first task is to calculate the scaling coefficient
$\alpha ={\mu}_{s}^{\mathit{max}}\u2215{\mu}_{s}$
$(\alpha \u2a7e1)$
. The
$r$
,
$t$
pairs within the spatial interval
$\rho -2{R}_{f}<\alpha r<\rho +2{R}_{f}$
are extracted and scaled (
${r}^{\prime}=\alpha r$
,
${t}^{\prime}=\alpha t$
). As the
$r$
,
$t$
pairs are sorted with respect to
$r$
, the selecting process turns into a simple task of finding the borders of the spatial interval, which is done quickly using a standard binary search algorithm. The interval itself is motivated by the convolution of each detection event with the size of the source fiber combined with the size of the detector fiber. A photon will hence contribute to the time dispersion histogram as soon as its distance to the center of the detector fiber is less than
$2{R}_{f}$
. Each photon detection event is assigned a weight,
${w}_{f}\left({r}^{\prime}\right)$
, according to its spatial position. This weight is based on the overlap between two circles of radius
${R}_{f}$
, one centered at the detection event and one centered at the detector location. More precisely, the weight is reached after integration over all origin-to-detector angles (rotating the detector circle around the
$z$
axis). This is similar to work by Wang
^{19} and Prahl.^{20} As the fiber radius and the source-detector separation is known when evaluating data, the weight
${w}_{f}\left({r}^{\prime}\right)$
can be precalculated, and the computationally costly integrals do not have to be evaluated for each photon. A uniform near-field irradiance distribution of the fibers was assumed in this work, although any distribution could be utilized with minimal modifications.

The photon detection events are also assigned a weight due to absorption, ${w}_{a}\left({t}^{\prime}\right)=\mathrm{exp}(-{\mu}_{a}{c}^{\prime}{t}^{\prime})$ , where ${c}^{\prime}$ is the local speed of light. The total weight for each detection event is thus ${w}_{\mathit{tot}}={w}_{f}\cdot {w}_{a}$ .

The final step is a simple matter of creating a weighted histogram of the ( ${t}^{\prime}$ , ${w}_{\mathit{tot}}$ ) data, where the temporal channel width, $\Delta t$ , can be specified by the user.

## 2.1.3.

#### Entangled detection events

As mentioned earlier, a single photon may generate multiple detection events (a primary event, a secondary event, and so on). All such events should contribute to the final time-of-flight histogram (given that they are located within the detection range $\rho -2{R}_{f}<{r}^{\prime}<\rho +2{R}_{f}$ ). However, if the spatial distance between two such events is less than $2{R}_{f}$ , their weights cannot be considered independent, and there is a risk for erroneous additions of weight to the final histogram. The worst case is if the two events coincide spatially (zero event distance), a case in which the entire weight contributed by the secondary event is erroneous. In contrast, when the spatial event separation is $2{R}_{f}$ or greater, none of its contribution is erroneous. For event separations between zero and $2{R}_{f}$ , the matter is nontrivial. In general, the erroneously added weight decreases with event separation. An estimate of the total amount of erroneous weight is reached by studying the multiple detection events occurring during the WMC simulations. For the case of ${\mu}_{s}^{\prime}=10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , $g=0.7$ , $NA=0.29$ , and ${R}_{f}=0.3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , we found that 0.5% of the primary events are accompanied by a secondary event less than $2{R}_{f}$ away (and only 0.1% within ${R}_{f}$ ). Assuming that all corresponding weight is erroneous (clearly an exaggeration of the problem), we can note that the final photon time-of-flight histogram contains less than 0.5% erroneous weight. Although this suggests that the problem of entangled detection events is insignificant, it may deserve further attention.

## 2.1.4.

#### Verification

In an effort to validate the WMC simulation program code and the implementation of the scaling relationships, the open source program MCML was carefully modified to monitor the time-of-flight histograms of photons leaving a homogenously scattering, semi-infinite medium. The photon weights leaving the medium within the $1\text{-}\mathrm{mm}$ -wide spatial bins, centered around 10, 15, and $20\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , were recorded with $10\text{-}\mathrm{ps}$ temporal resolution. ${10}^{8}$ photons were simulated at $g=0.7$ and $n=1.33$ for all combinations of ${\mu}_{s}^{\prime}=5$ , $10\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ and ${\mu}_{a}=0.1$ , $0.5\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ .

## 2.2.

### Time-Resolved Instrumentation

Time-resolved photon migration experiments were conducted using a compact
$(50\times 50\times 30\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{3})$
and portable time-domain photon migration instrument primarily intended for spectroscopy of biological tissues in clinical environments. Detailed information on the instrumentation can be found in previous publications.^{5, 21, 22} Briefly, the system is based on pulsed diode laser technology and time correlated single photon counting (TCSPC). Four pulsed diode lasers (at 660, 786, 830, and
$916\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
, respectively) generate pulses having a FWHM of about
$70\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$
, the average power being
$1\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}2\phantom{\rule{0.3em}{0ex}}\mathrm{mW}$
. Light is injected into the sample and collected using
$600\text{-}\mu \mathrm{m}$
GRIN optical fibers
$(NA=0.29)$
. A fast MCP-PMT together with a TCSPC computer card is used to obtain photon time-of-flight histograms. Broadening in fibers and detector yields an instrument response function (IRF) of about
$100\text{-}\mathrm{ps}$
FWHM. The IRF is measured by putting the fiber ends face to face with a thin paper coated on both sides with black toner, similar to the IRF measurements proposed by Schmidt
^{23} A schematic illustration of the setup is given in Fig. 3
.

## 2.3.

### Experiment

In order to compare diffusion-based and WMC-based evaluation of experimental data, measurements were performed on tissue phantoms based on Intralipid (Fresenius Kabi
$200\phantom{\rule{0.3em}{0ex}}\mathrm{mg}\u2215\mathrm{ml}$
) and ink (1:100 stock solution of Pelikan Fount India ink).^{24}

Two measurement series were performed: one added absorber series using ink,^{25} and one added scatterer series using Intralipid.^{26, 27} Neglecting the minor volume change occurring during these measurement series, we can qualitatively expect a constant scattering and linearly increasing absorption in the added absorber series. In addition, when extrapolating of derived
${\mu}_{a}$
toward the zero ink level, the absorption should be close to the absorption exhibited by pure water.

Similarly, the added scatterer series should yield a constant absorption, linearly increasing scattering. Extrapolation of the derived scattering coefficient toward the zero Intralipid level should yield a zero scattering (at least if ink scattering is negligible).

The added absorber series was performed on a phantom based on $577\phantom{\rule{0.3em}{0ex}}\mathrm{ml}$ water and $22\phantom{\rule{0.3em}{0ex}}\mathrm{ml}$ Intralipid. Measurements were performed at $2\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}8\phantom{\rule{0.3em}{0ex}}\mathrm{ml}$ total volume of added ink solution ( $1\text{-}\mathrm{ml}$ increments). The added scatterer series was performed on a phantom based on $577\phantom{\rule{0.3em}{0ex}}\mathrm{ml}$ water and $4\phantom{\rule{0.3em}{0ex}}\mathrm{ml}$ ink solution. Measurements were performed at $10\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}24\phantom{\rule{0.3em}{0ex}}\mathrm{ml}$ added volume of Intralipid (in $2\text{-}\mathrm{ml}$ increments). A cylindrical plastic container, $\xd8=110\text{-}\mathrm{mm}$ and $\text{height}=70\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , was used to hold the phantoms during mixing and measurements. The container was placed on a magnetic stirrer, utilized for mixing of the phantoms after each addition of ink or Intralipid. Fibers were placed in parallel at $30\text{-}\mathrm{mm}$ depth in the center of the phantom, separated $14.7\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , center to center. Data acquisition was performed for $30\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ for each measurement, and the IRF was measured before, after, and occasionally between the measurements (monitoring potential temporal drifts in the system).

## 2.4.

### Modeling

The agreement between the diffusion approximation and MC is investigated by fitting analytical solutions of the diffusion equation to the impulse responses as delivered by WMC (along the lines of, e.g., Refs. 3, 12). Such diffusion modeling is based on iterative nonlinear Levenberg-Marquardt curve fitting, in which
${\mu}_{a}$
,
${\mu}_{s}^{\prime}$
, and an amplitude factor,
$k$
, are varied in order to minimise a
${\chi}^{2}$
merit function.^{3, 28} For the interstitial case (infinite geometry, see, e.g., Ref. 1), the form of the diffusion-based impulse response is stated in Eq. 3:

## Eq. 3

$$y({\mu}_{a},{\mu}_{s}^{\prime},k,t)=k{t}^{-3\u22152}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-\frac{3{\mu}_{s}^{\prime}{\rho}^{2}}{4{c}^{\prime}t}-{\mu}_{a}{c}^{\prime}t),$$When WMC is used for modeling, the fitting procedure is based on an exhaustive search over ${\mu}_{s}^{\prime}$ . That is, for each value in a set of trial values of ${\mu}_{s}^{\prime}$ , the optimal choice of ${\mu}_{a}$ and $k$ is determined in a suboptimization procedure. The suboptimized error norm, ${\stackrel{\u0303}{\chi}}^{2}\left({\mu}_{s}^{\prime}\right)$ , is given in Eq. 4:

## Eq. 4

$${\stackrel{\u0303}{\chi}}^{2}\left({\mu}_{s}^{\prime}\right)=\underset{k,{\mu}_{a}}{\mathrm{min}}\phantom{\rule{0.2em}{0ex}}\left\{{\chi}^{2}(k,{\mu}_{a},{\mu}_{s}^{\prime})\right\},$$This article uses an MC database for infinite geometries (interstitial measurements) consisting of about $8\times {10}^{7}$ detection events, originating from a simulation of $2\times {10}^{8}$ photons. The database for semi-infinite evaluations is of a similar size, resulting from a simulation of ${10}^{9}$ photons. Both these databases are based on simulations where $g=0.7$ , $n=1.33$ , $NA=0.29$ , ${\mu}_{s}^{\mathit{max}}=90\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , and ${t}_{\mathit{max}}=2\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ . A single forward modeling (that is, a generation of a time-of-flight histogram from the previously stored database) can be performed in less than one second. Using the exhaustive search described earlier, the best fit is found in approximately $15\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ .

It should also be noted that, unless otherwise stated, all evaluations are based on fittings in the temporal range defined by the 50% level (of peak maximum) on the rising flank and the 20% level on the tail.

## 3.

## Results

## 3.1.

### White Monte Carlo Performance

The performance of our WMC approach was investigated and verified for the case of diffuse reflectance by comparing its output with that of MCML. The optical properties examined were combinations of ${\mu}_{s}^{\prime}=5$ , $10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ and ${\mu}_{a}=0.1$ , $0.5\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , as described in Sec. 2.1.4. Evaluating the time-resolved MCML data with WMC shows only minor differences, despite major differences in for example spatial binning. The relative difference between derived optical properties and MCML input parameters was only $1\pm 4\%$ in absorption, and $0\pm 2\%$ for reduced scattering $(\text{mean}\pm \mathrm{st.}\mathrm{dev.})$ . An illustration of this agreement is shown in Fig. 5 .

## 3.2.

### WMC Versus Diffusion Modeling

That the use of diffusion modeling is questionable in the low albedo regime is easily shown by simple visual comparison of the impulse responses as delivered by the diffusion approximation and WMC, respectively. An example of this is given in Fig. 6 , using optical properties encountered in the human prostate (at wavelengths in the range $650\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}900\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ ).

In order to quantitatively investigate the performance of diffusion modeling for the range of optical properties explored in this work, the analytical expression for the diffusion-based impulse response was fitted to data delivered by WMC (at $10\text{-}\mathrm{ps}$ time resolution). In order to mimic the experimental conditions, where the system IRF must be accounted for, data was convoluted with a synthetic, $70\text{-}\mathrm{ps}$ FWHM Gaussian prior to curve fitting. The parameter space covered in this investigation was:

In order to allow direct comparison with the diffusion modeling used in our previous work on *in vivo* characterization of human prostate tissue,^{5} fitting was performed in the range defined by the 50% level (of peak maximum) on the rising flank, down to 20% on the tail. The WMC parameters,
${\mu}_{\mathit{WMC}}$
are considered as true optical properties, and relative errors in derived optical properties due to diffusion modeling are calculated as given in:

The influence of the selection of fit range has been discussed in several articles.^{3, 12, 29} To address this issue, we present the performance of time-domain diffusion modeling also for the case of fitting in the range defined by the 80% level (of peak maximum) on the rising flank, down to the 1% level on the tail. The results are shown in Fig. 8
. Although the use of a reduced fit range slightly improves the performance of diffusion modeling, it is clear that significant overestimation remains.

In order to predict the outcome of the experimental procedures described in Sec. 2.3, we employ diffusion modeling to interpret sets of WMC data where (1) ${\mu}_{s}^{\prime}$ is kept constant, while ${\mu}_{a}$ gradually increases (added absorber series), and (2) ${\mu}_{a}$ is kept constant, while ${\mu}_{s}^{\prime}$ gradually increases (added scatterer series). These two simulation series correspond to two-dimensional (2-D) sections in the three-dimensional (3-D) maps shown in Figs. 7 and 8. The outcome is shown in Figs. 9 and 10 , clearly revealing the increasing model errors as we move further and further into the low albedo regime. For the added absorber series, note especially that the slope encountered in derived ${\mu}_{a}$ is steeper than the true increase, and that the constancy of ${\mu}_{s}^{\prime}$ is replaced by a fairly linear increase. For the added scatterer series at ${\mu}_{a}=0.5\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , note the offset-like pattern for ${\mu}_{s}^{\prime}$ , as well as the remarkably poor agreement between true and derived ${\mu}_{a}$ at low values of ${\mu}_{s}^{\prime}$ .

## 3.3.

### Experimental Results

The experimental results from the added absorber series, evaluated using both diffusion theory and WMC, are presented in Fig. 11 . Here, derived ${\mu}_{a}$ appears to increase linearly when employing WMC evaluation, while diffusion modeling appears to introduce a slight nonlinear increase. In addition, and in good agreement with the predictions presented in Sec. 3.2, diffusion evaluation gives rise to a steeper slope. Linear extrapolation of the WMC-derived absorption coefficient to the zero ink level produces estimations of background absorption as stated in Table 1 . Turning to derived scattering, WMC produces almost constant estimations of ${\mu}_{s}^{\prime}$ , whereas diffusion theory undergoes a significant increase as more ink is added. Also this phenomenon agrees with our predictions.

## Table 1

Background absorption extrapolated from added absorber measurements (evaluated using WMC). As water is the principal absorber in the ink-free phantom, water absorption coefficients are stated for reference.

λ(nm) | μa(cm−1) | |
---|---|---|

μa(0) | Water | |

660 | 0.002 | 0.004 |

786 | 0.016 | 0.022 |

830 | 0.024 | 0.029 |

916 | 0.077 | 0.091 |

The results from the added scatterer series are shown in Fig. 12 . As predicted by our simulations, absorption coefficients derived with diffusion theory decrease heavily with the volume of added Intralipid. In contrast, ${\mu}_{a}$ remains constant when derived using WMC evaluation. Also, when turning the attention to scattering, the predictions from Sec. 3.2 hold. Although the response is linear, diffusion modeling produces significantly larger values of ${\mu}_{s}^{\prime}$ , in an offset-like manner (clearly not approaching zero scattering when extrapolating the Intralipid content to the zero level). A small offset is present also for WMC evaluation and may be related to ink scattering or temporal drifts.

When using the
${\mu}_{s}^{\prime}$
slope (as provided by linear regression) to estimate the scattering of Intralipid, WMC evaluation yields values in good agreement with the levels reported by van Staveren, ^{26} as shown in Table 2
.

## Table 2

Scattering of 200mg∕l Intralipid.

λ(nm) | μs′[cm−1∕(ml∕l)] | |
---|---|---|

This Work | Van Staveren 26 | |

660 | 0.227 | 0.25 |

786 | 0.171 | 0.20 |

830 | 0.166 | 0.19 |

916 | 0.138 | 0.17 |

## 4.

## Discussion

This work establishes White Monte Carlo as a tool for evaluation of experimental time-resolved photon migration. At the same time, we clearly show that the diffusion approximation of light propagation should be used with great care when evaluating time-resolved spectroscopy in highly absorbing turbid media. For example, regardless of the fitting range used, the errors in derived optical properties due to diffusion modeling are about 10% or more in the regime of interest when modeling light propagation in the human prostate. Our hope is that this development will promote accurate assessment of, for example, the *in vivo* optical properties of human prostate tissue.

The minor differences observed when comparing MCML and WMC indicate that our approach, including both simulation and scaling procedures, is valid. The small discrepancies can be explained either by differences in spatial sampling or by limited photon statistics (signal to noise). It should be noted that the generated MCML data sets correspond to an infinitely small source and equal weighting of all photons appearing at a radial distance within $0.5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ from the source-detector separation $\rho $ . In contrast, our WMC simulations assumed that the source and detector fibers had a radius of $0.3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , performing individual photon weighting based on the spatial (radial) location of the detection event. Furthermore, the numerical apertures of the involved optical fibers were also accounted for. In fact, the limited differences suggest that the elaborate procedures used in order to take source and detector characteristics into account might be unnecessary. However, there might exist regimes of optical properties, or source-detector separations, where this still is a valuable approach. In this context, it should also be mentioned that the use of fixed acceptance cones (rather than $NA$ -based) would allow scalability also with respect to the refractive index.

The inability of scaling Monte Carlo simulations with respect to the anisotropy factor
$g$
has been discussed in earlier work on White Monte Carlo.^{10, 12} The question is how a discrepancy between the actual
$g$
and the value used in simulations will influence derived values of
${\mu}_{a}$
and
${\mu}_{s}^{\prime}$
. However, light propagation in the photon migration regime is known to fall under similarity as long as
${\mu}_{s}^{\prime}=(1-g){\mu}_{s}$
is conserved.^{30} This means that an exact knowledge of
$g$
is not crucial and that a single WMC database may be used to model materials exhibiting different anisotropies. Pifferi
^{12} argue that when using diffuse reflectance to solve for
${\mu}_{s}^{\prime}$
, the influence of the exact value of
$g$
is small at least as long as
$0.7<g<0.9$
. Similar findings are reported by Kienle and Patterson,^{10} who concluded similarity in diffuse reflectance for short fiber separations
$\left(2.25\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}4.75\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\right)$
as long as
$g>0.8$
.

The WMC-based data evaluation of the added absorber and added scatterer series removes the expected, but erroneous, patterns in derived optical properties exhibited when basing modeling on diffusion theory. In the added absorber series, WMC produces a highly linear increase in derived
${\mu}_{a}$
and a constant level of
${\mu}_{s}^{\prime}$
. Extrapolating measured
${\mu}_{a}$
indicates a background absorption close to that of pure water. For the case of the added scatterer series, WMC evaluation results in a linearly increasing
${\mu}_{s}^{\prime}$
together with a constant level of absorption. The small offsets in derived
${\mu}_{s}^{\prime}$
when extrapolating toward the zero Intralipid level remains unexplained and may originate from temporal drifts and/or a nonzero scattering contribution from ink. In total, this forms a strong argument that the proposed WMC approach provides accurate estimations of optical characteristics. This is further supported by the fact that generated estimations of Intralipid scattering are in good agreement with previously published values (as shown in Table 2). However, a more detailed and quantitative discussion on accuracy is aggravated by the uncertainties in the optical properties of the utilized optical phantoms. The scattering of Intralipid is somewhat debated and may depend on batch-to-batch variations.^{26, 27} Moreover, ink is known to scatter light, making it difficult to accurately assess its absorption.^{25}

Regarding the fit range, it has been shown that early photons should be excluded when using diffusion as the forward model. Cubeddu
^{29} suggest the use of the range between the 80% level (of peak maximum) on the rising flank down to the 1% level on the tail
$(80\u22151)$
. In this work, we mainly use a
$50\u221520$
range. This is to allow direct comparison with previously published modeling of *in vivo* time-resolved spectroscopy of the human prostate.^{5} Although such a choice might be questioned, the
$80\u22151$
fitting range still suffers from significant model errors. This in combination with the availability of WMC evaluation renders further discussions on the optimal range for diffusion-based modeling somewhat irrelevant.

Furthermore, we would like to call for some attention regarding the step-sizes used in Monte Carlo simulations. Swartling questions whether the White Monte Carlo approach, including all scalable Monte Carlo approaches based on simulations in nonabsorbing media, is equivalent to the *de facto* standard Monte Carlo approach.^{14} Conventional Monte Carlo utilizes an average step-size of
$1\u2215{\mu}_{t}=1\u2215({\mu}_{a}+{\mu}_{s})$
, whereas a WMC approach uses a
$1\u2215{\mu}_{s}$
average step-size. This means that the photon paths generated by the two methods will differ if
${\mu}_{a}$
is nonzero. Since traditional Monte Carlo for light propagation in turbid media assumes that scattering dominates absorption (
${\mu}_{s}\u2aa2{\mu}_{a}$
, high albedo), this has not been considered as an important source of discrepancy.^{31} Furthermore, the methods for reduction of photon weights differ. Traditional MC reduces the photon weight by a factor of
$a$
at each scattering event [
$a={\mu}_{s}\u2215({\mu}_{a}+{\mu}_{s})$
being the albedo], while WMC follows the conventional Beer-Lambert law of absorption. Swartling argues that the two ways of absorbing weight are equivalent at the extremes, that is, when
$a\to 1$
or
$a\to 0$
, but notes that equivalence is not guaranteed in the intermediate region where
${\mu}_{s}\approx {\mu}_{a}$
. Since WMC does not rely on the high albedo assumption, this method should be more reliable in the regime when scattering and absorption coefficients are of comparable magnitude. (For an example of nonwhite MC without being based on the high albedo assumption, we refer to the work of Farina
^{32}). In most types of biological tissue, the difference between traditional MC and WMC can, however, be disregarded.

Finally, it should be noted that if the Monte Carlo approach described in this work is considered for use in steady-state investigations, the procedures for photon termination deserve extra attention. Terminating photons when the time-of-flight exceeds ${t}_{\mathit{max}}$ may impose errors during steady-state analysis of generated WMC databases.

## 5.

## Conclusion

We have developed a quick and accurate scheme for Monte Carlo–based evaluation of experimental time-resolved spectroscopy of highly scattering materials. We show, for the first time, evaluations of experimental interstitial time-domain photon migration using fully scalable White Monte Carlo. The great value of this approach is demonstrated in a regime where diffusion-based modeling results in fairly large errors (overestimations of both ${\mu}_{s}^{\prime}$ and ${\mu}_{a}$ ). The relevance of this regime is, for example, motivated by the recent interest in the optical properties of human prostate tissue, where ${\mu}_{a}>0.3\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ and ${\mu}_{s}^{\prime}<10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ .

## Acknowledgments

We gratefully acknowledge financial support from the Wallenberg Foundation and through the EC Grant Nano-UB Sources (IST-2005-017128).

## References

**,” Appl. Opt., 28 (12), 2331 –2336 (1989). 0003-6935 Google Scholar**

*Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical-properties***,” IEEE Trans. Biomed. Eng., 36 (12), 1162 –1168 (1989). https://doi.org/10.1109/TBME.1989.1173624 0018-9294 Google Scholar**

*Monte Carlo modeling of light-propagation in highly scattering tissues: I. model predictions and comparison with diffusion-theory***,” Phys. Med. Biol., 40 (11), 1957 –1975 (1995). https://doi.org/10.1088/0031-9155/40/11/013 0031-9155 Google Scholar**

*The influence of boundary-conditions on the accuracy of diffusion-theory in time-resolved reflectance spectroscopy of biological tissues***,” J. Urol. (Baltimore), 143 (2), 398 –401 (1990). 0022-5347 Google Scholar**

*Photodynamic therapy for localized prostatic cancer—light penetration in the human prostate-gland***,” J. Biomed. Opt., 12 (1), 014022 (2007). https://doi.org/10.1117/1.2435175 1083-3668 Google Scholar**

*In vivo*optical characterization of human prostatic tissue using near-infrared time-resolved spectroscopy**,” Med. Phys., 10 (6), 824 –830 (1983). https://doi.org/10.1118/1.595361 0094-2405 Google Scholar**

*A Monte Carlo model for the absorption and flux distributions of light in tissue***,” Comput. Methods Programs Biomed., 47 (2), 131 –146 (1995). https://doi.org/10.1016/0169-2607(95)01640-F 0169-2607 Google Scholar**

*MCML Monte Carlo modeling of light transport in multilayered tissues***,” J. Opt. Soc. Am. A, 20 (4), 714 –727 (2003). https://doi.org/10.1364/JOSAA.20.000714 0740-3232 Google Scholar**

*Accelerated Monte Carlo models to simulate fluorescence spectra from layered tissues***,” Appl. Spectrosc., 58 (5), 591 –597 (2004). https://doi.org/10.1366/000370204774103426 0003-7028 Google Scholar**

*Photon migration in Raman spectroscopy***,” Phys. Med. Biol., 41 (10), 2221 –2227 (1996). https://doi.org/10.1088/0031-9155/41/10/026 0031-9155 Google Scholar**

*Determination of the optical properties of turbid media from a single Monte Carlo simulation***,” Trends in Optics and Photonics: Advances in Optical Imaging and Photon Migration, 1996). Google Scholar**

*Fitting of time-resolved reflectance curves with a Monte Carlo model***,” Appl. Opt., 37 (13), 2774 –2780 (1998). 0003-6935 Google Scholar**

*Real-time method for fitting time-resolved reflectance and transmittance measurements with a Monte Carlo model***,” Appl. Opt., 32 (4), 426 –434 (1993). 0003-6935 Google Scholar**

*Condensed Monte Carlo simulations for the description of light transport***,” Lund University, (2002). Google Scholar**

*Biomedical and atmospheric applications of optical spectroscopy in scattering media***,” J. Biomed. Opt., 11 (4), 041104 (2006). https://doi.org/10.1117/1.2241609 1083-3668 Google Scholar**

*Investigation of light propagation models to determine the optical properties of tissue from interstitial frequency domain fluence measurements***,” Dosimetry of Laser Radiation in Medicine and Biology, IS 5 102 –111 SPIE, Bellingham, WA (1989). Google Scholar**

*A Monte Carlo model of light propagation in tissue***,” (2006). Google Scholar**

*SIMD-oriented fast Mersenne twister: a $128\text{-}\text{bit}$ pseudorandom number generator***,” Comput. Methods Programs Biomed., 54 141 –150 (1997). https://doi.org/10.1016/S0169-2607(97)00021-7 0169-2607 Google Scholar**

*CONV-convolution for responses to a finite diameter photon beam incident on multi-layered tissues***,” Phys. Med. Biol., 50 (11), 2559 –2571 (2005). https://doi.org/10.1088/0031-9155/50/11/008 0031-9155 Google Scholar**

*Characterization of normal breast tissue heterogeneity using time-resolved near-infrared spectroscopy***,” Appl. Opt., 44 (11), 2104 –2114 (2005). https://doi.org/10.1364/AO.44.002104 0003-6935 Google Scholar**

*Performance assessment of photon migration instruments: the MEDPHOT protocol***,” Rev. Sci. Instrum., 71 (1), 256 –265 (2000). https://doi.org/10.1063/1.1150191 0034-6748 Google Scholar**

*A 32-channel time-resolved instrument for medical optical tomography***,” J. Biomed. Opt., 11 (4), 041102 (2006). https://doi.org/10.1117/1.2335429 1083-3668 Google Scholar**

*Review of tissue simulating phantoms for optical spectroscopy, imaging, and dosimetry***,” Phys. Med. Biol., 37 (4), 985 –993 (1992). https://doi.org/10.1088/0031-9155/37/4/012 0031-9155 Google Scholar**

*The use of india ink as an optical absorber in tissue-simulating phantoms***,” Appl. Opt., 30 (31), 4507 –4514 (1991). 0003-6935 Google Scholar**

*Light scattering in intralipid-10-percent in the wavelength range of $400\u20131100\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$***,” Lasers Surg. Med., 12 (5), 510 –519 (1992). https://doi.org/10.1002/lsm.1900120510 0196-8092 Google Scholar**

*Optical properties of intralipida phantom medium for light-propagation studies***,” Opt. Express, 12 (17), 4103 –4112 (2004). https://doi.org/10.1364/OPEX.12.004103 1094-4087 Google Scholar**

*Time and wavelength resolved spectroscopy of turbid media using light continuum generated in a crystal fiber***,” Med. Phys., 23 (9), 1625 –1633 (1996). https://doi.org/10.1118/1.597739 0094-2405 Google Scholar**

*Experimental test of theoretical models for time-resolved reflectance***,” Opt. Eng., 32 (2), 244 –252 (1993). https://doi.org/10.1117/1.2170988 0091-3286 Google Scholar**

*Similarity relations for anisotropic scattering in absorbing media***,” Optical-Thermal Response of Laser-Irradiated Tissue, 73 –100 Plenum Press, New York (1995). Google Scholar**

*Monte Carlo modeling of light transport in tissues***,” Phys. Med. Biol., 44 1 –11 (1999). https://doi.org/10.1088/0031-9155/44/1/002 0031-9155 Google Scholar**

*Monte Carlo simulation of light fluence in tissue in a cylindrical diffusing fiber geometry*