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.¹ 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{\mathrm{cm}}^{-1}$ and absorption above $0.3{\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,⁶ Monte Carlo (MC) simulation has become the gold standard for modeling of light propagation in tissue optics.⁷ Besides modeling spatially and temporally resolved light distribution, MC has proven useful in, for example, fluorescence modeling⁸ and Raman spectroscopy.⁹

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 .

Fig. 1

In WMC, the shapes of the photon paths are determined only by the phase scattering function and the sequence of random numbers. Paths generated by simulations in nonabsorbing media can hence be linearly scaled to apply for different values of ${\mu }_s$ . The illustration is adapted from Ref. 12.

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:

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.¹³ 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.¹⁰ 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.⁸ Swartling also raised the important question regarding the equivalency of WMC and traditional MC.¹⁴ 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¹⁵). 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}}∕\alpha$ results in data valid in the time interval $[0,\alpha t_{\mathit{max}}]$ .

Fig. 2

A flowchart of the WMC scheme used. The WMC simulation program performs the simulation in either infinite or semi-infinite homogenously scattering media, using user-supplied WMC input data. The resulting database is sorted and stored to be used later by the forward model. The forward model rescales the database data and generates time-of-flight histograms corresponding to parameters supplied by the user. These parameters include the radius $R_f$ of the involved optical fibers, the source-detector separation $\rho$ , the optical properties ${\mu }_s^{\prime }$ and ${\mu }_a$ , and the desired temporal channel width $\Delta t$ . As this fast forward model is incorporated in the forward solver, the method can be used to evaluate experimental time-domain data.

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⁷). As the photon propagation is similar to traditional MC, we refer to the work of Wang ⁷ and Prahl ¹⁶ 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.¹⁷ 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}}=\arcsin (NA∕n)$ . The angular distribution is assumed flat, i.e., $\cos \theta =1-\xi (1-\cos {\theta }_{\mathit{max}})$ , where $\xi$ is a random number, $\xi ϵ[0,1)$ . Although unnecessary in a cylindrically symmetric problem, the azimuthal angle $\phi$ is also randomized, $\phi =2\pi \xi$ , where $\xi ϵ[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:

Eq. 2

$$\Delta s=\frac{-\ln \left(\xi \right)}{{\mu }_s^{\max }},$$

where $\xi$ is a random variable, $\xi ϵ(0,1]$ , implying $\Delta sϵ[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\log \left(n\right)\right]$ worst-case time complexity} in-place sorting algorithm Combsort11.¹⁸ 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}}∕{\mu }_s$ $(\alpha ⩾1)$ . The $r$ , $t$ pairs within the spatial interval $\rho -2R_f<\alpha r<\rho +2R_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 $2R_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 ¹⁹ and Prahl.²⁰ 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)=\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 -2R_f<r^{\prime }<\rho +2R_f$ ). However, if the spatial distance between two such events is less than $2R_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 $2R_f$ or greater, none of its contribution is erroneous. For event separations between zero and $2R_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{\mathrm{cm}}^{-1}$ , $g=0.7$ , $NA=0.29$ , and $R_f=0.3\mathrm{mm}$ , we found that 0.5% of the primary events are accompanied by a secondary event less than $2R_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\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{\mathrm{cm}}^{-1}$ and ${\mu }_a=0.1$ , $0.5{\mathrm{cm}}^{-1}$ .

2.2.

Time-Resolved Instrumentation

Time-resolved photon migration experiments were conducted using a compact $(50\times 50\times 30{\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\mathrm{nm}$ , respectively) generate pulses having a FWHM of about $70\mathrm{ps}$ , the average power being $1\text{to}2\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 ²³ A schematic illustration of the setup is given in Fig. 3 .

Fig. 3

A schematic of the instrumentation in interstitial mode.

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∕2}\exp (-\frac{3{\mu }_s^{\prime }{\rho }^2}{4c^{\prime }t}-{\mu }_a{c}^{\prime }t),$$

where $\rho$ is the source-detector separation, and $c^{\prime }$ is the speed of light in the sample. To reach this expression, we use the absorption-independent definition of the diffusion coefficient, $D={\left(3{\mu }_s^{\prime }\right)}^{-1}$ . During evaluation of experimental data, the procedure described earlier involves a convolution with the system IRF.

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, ${\tilde{\chi }}^2\left({\mu }_s^{\prime }\right)$ , is given in Eq. 4:

Eq. 4

$${\tilde{\chi }}^2\left({\mu }_s^{\prime }\right)=\underset{k,{\mu }_a}{\min }\left\{{\chi }^2(k,{\mu }_a,{\mu }_s^{\prime })\right\},$$

and is thus to be exhaustively searched on some appropriate ${\mu }_s^{\prime }$ grid. In this work, ${\tilde{\chi }}^2\left({\mu }_s^{\prime }\right)$ is evaluated in steps of $0.05{\mathrm{cm}}^{-1}$ . An example of this procedure is given in Fig. 4 .

Fig. 4

WMC fitting procedure of the experimental data featuring the $2\text{-}\mathrm{ml}$ added ink absorption series measurement at $660\mathrm{nm}$ , $\rho =14.7\mathrm{mm}$ .

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{\mathrm{cm}}^{-1}$ , and $t_{\mathit{max}}=2\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\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{\mathrm{cm}}^{-1}$ and ${\mu }_a=0.1$ , $0.5{\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 .

Fig. 5

Illustration of the agreement between MCML and the proposed WMC. In this particular case, the source-detector separation is $15\mathrm{mm}$ , ${\mu }_s^{\prime }=5{\mathrm{cm}}^{-1}$ , and ${\mu }_a=0.1{\mathrm{cm}}^{-1}$ .

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\text{to}900\mathrm{nm}$ ).

Fig. 6

Illustration of the breakdown of the diffusion approximation in the low albedo regime of photon migration. In this particular case, the data correspond to an infinite (interstitial) geometry using a $15\text{-}\mathrm{mm}$ source-detector separation, with ${\mu }_s^{\prime }=7.5{\mathrm{cm}}^{-1}$ , and ${\mu }_a=0.5{\mathrm{cm}}^{-1}$ .

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:

$$0.01⩽{\mu }_a⩽1,$$

$$5⩽{\mu }_s^{\prime }⩽15\left[{\mathrm{cm}}^{-1}\right].$$

In order to allow direct comparison with the diffusion modeling used in our previous work on in vivo characterization of human prostate tissue,⁵ 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:

Eq. 5

$$\Delta \mu =\frac{{\mu }_D-{\mu }_{\mathit{WMC}}}{{\mu }_{\mathit{WMC}}},$$

where ${\mu }_D$ denotes the optical properties as derived by diffusion modeling. The result is shown in Fig. 7 and suggests that diffusion modeling will result in fairly large overestimations of both ${\mu }_a$ and ${\mu }_s^{\prime }$ .

Fig. 7

Relative errors due to diffusion modeling. Since the relative errors are positive throughout the examined range, diffusion modeling results in overestimation of both scattering and absorption. Fitting is performed between the 50% level (of peak maximum) on the rising flank and the 20% level on the tail. The gray scale applies to both graphs.

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.

Fig. 8

Relative errors due to diffusion modeling when fitting is performed between the 80% level (of peak maximum) on the rising flank and the 1% level on the tail. As in Fig. 7, diffusion modeling results in overestimation of both scattering and absorption. The gray scale applies to both graphs.

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{\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 }$ .

Fig. 9

Derived optical properties for two levels of ${\mu }_s^{\prime }$ when diffusion modeling is used to interpret WMC data (corresponds to the added absorber series).

Fig. 10

Derived optical properties for two levels of ${\mu }_a$ when diffusion modeling is used to interpret WMC data (corresponds to the added scatterer series).

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.

Fig. 11

The results from the added absorber series (adding ink to a liquid Intralipid phantom). Derived optical properties are given both for diffusion-based evaluation (crosses) and WMC-based evaluation (circles). The results are to be compared with the simulation results presented in Fig. 9.

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.

Fig. 12

The results from the added scatterer series (adding further Intralipid to a liquid phantom based on ink and Intralipid). Derived optical properties are given both for diffusion-based evaluation (crosses) and WMC-based evaluation (circles). The results are to be compared with the simulation results presented in Fig. 10.

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, ²⁶ 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\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\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.³⁰ 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 ¹² 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,¹⁰ who concluded similarity in diffuse reflectance for short fiber separations $\left(2.25\text{to}4.75\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.²⁵

Regarding the fit range, it has been shown that early photons should be excluded when using diffusion as the forward model. Cubeddu ²⁹ 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∕1)$ . In this work, we mainly use a $50∕20$ range. This is to allow direct comparison with previously published modeling of in vivo time-resolved spectroscopy of the human prostate.⁵ Although such a choice might be questioned, the $80∕1$ 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.¹⁴ Conventional Monte Carlo utilizes an average step-size of $1∕{\mu }_t=1∕({\mu }_a+{\mu }_s)$ , whereas a WMC approach uses a $1∕{\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⪢{\mu }_a$ , high albedo), this has not been considered as an important source of discrepancy.³¹ 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∕({\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 ³²). 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{\mathrm{cm}}^{-1}$ and ${\mu }_s^{\prime }<10{\mathrm{cm}}^{-1}$ .