Translator Disclaimer
1 July 2008 White Monte Carlo for time-resolved photon migration
Author Affiliations +
A novel scheme for fully scalable White Monte Carlo (WMC) has been developed and is used as a forward solver in the evaluation of experimental time-resolved spectroscopy. Previously reported scaling problems are avoided by storing detection events individually, turning spatial and temporal binning into post-simulation activities. The approach is suitable for modeling of both interstitial and noninvasive settings (i.e., infinite and semi-infinite geometries). Motivated by an interest in in vivo optical properties of human prostate tissue, we utilize WMC to explore the low albedo regime of time-domain photon migration—a regime where the diffusion approximation of radiative transport theory breaks down, leading to the risk of overestimating both reduced scattering (μs) and absorption (μa). Experimental work supports our findings and establishes the advantages of Monte Carlo–based evaluation.



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 (μa) and the reduced scattering coefficient (μs) 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 10cm1 and absorption above 0.3cm1 .

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

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 μ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, μ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 wa according to the Beer-Lambert law of attenuation. This weight is stated in Eq. 1, where c is the speed of light within the media:

Eq. 1

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 μs and μ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 μa and μs . In practice, however, they based their evaluation on interpolation between MC simulations carried out at different μs , thus utilizing scalability in μa only. On the other hand, their work included evaluation of experimental time-resolved transmittance in the slab geometry (a geometry not allowing μ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 program15). 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 μ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 μs scaling. It also allows an accurate account for finite extension of source and detector areas (e.g., optical fiber diameters).


Materials and Methods


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 μs and μ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 Rf ) 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 μsmax but is always terminated when time reaches tmax . This notation is used because μsmax defines the maximum scattering coefficient for which the resulting database can provide valid data throughout the time interval [0,tmax] . Generally, scaling to a certain scattering μs=μsmaxα results in data valid in the time interval [0,αtmax] .

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 Rf of the involved optical fibers, the source-detector separation ρ , the optical properties μs and μa , and the desired temporal channel width Δt . As this fast forward model is incorporated in the forward solver, the method can be used to evaluate experimental time-domain data.



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 (MCML7). 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 21320491 period, documented excellent equidistribution properties, and fast random number generation. The entire 32-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, θ , representing the emission cone of the source fiber, defined by the NA of the fiber, θmax=arcsin(NAn) . The angular distribution is assumed flat, i.e., cosθ=1ξ(1cosθmax) , where ξ is a random number, ξϵ[0,1) . Although unnecessary in a cylindrically symmetric problem, the azimuthal angle φ is also randomized, φ=2πξ , where ξϵ[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, Δs , is calculated using the scattering coefficient instead of the total attenuation coefficient. It is defined in Eq. 2:

Eq. 2

where ξ is a random variable, ξϵ(0,1] , implying Δsϵ[0,) . 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-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 tmax . 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 tmax .


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[nlog(n)] worst-case time complexity} in-place sorting algorithm Combsort11.18 The O(1) 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: μa and μs being the optical properties corresponding to the resulting time-of-flight distribution, the temporal channel width Δt , source-detector separation ρ , and the fiber radius Rf (experiments typically involve two identical fibers). The first task is to calculate the scaling coefficient α=μsmaxμs (α1) . The r , t pairs within the spatial interval ρ2Rf<αr<ρ+2Rf are extracted and scaled ( r=αr , t=α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 2Rf . Each photon detection event is assigned a weight, wf(r) , according to its spatial position. This weight is based on the overlap between two circles of radius Rf , 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 wf(r) 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, wa(t)=exp(μact) , where c is the local speed of light. The total weight for each detection event is thus wtot=wfwa .

The final step is a simple matter of creating a weighted histogram of the ( t , wtot ) data, where the temporal channel width, Δt , can be specified by the user.


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 ρ2Rf<r<ρ+2Rf ). However, if the spatial distance between two such events is less than 2Rf , 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 2Rf or greater, none of its contribution is erroneous. For event separations between zero and 2Rf , 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 μs=10cm1 , g=0.7 , NA=0.29 , and Rf=0.3mm , we found that 0.5% of the primary events are accompanied by a secondary event less than 2Rf away (and only 0.1% within Rf ). 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.



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-mm -wide spatial bins, centered around 10, 15, and 20mm , were recorded with 10-ps temporal resolution. 108 photons were simulated at g=0.7 and n=1.33 for all combinations of μs=5 , 10cm1 and μa=0.1 , 0.5cm1 .


Time-Resolved Instrumentation

Time-resolved photon migration experiments were conducted using a compact (50×50×30cm3) 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 916nm , respectively) generate pulses having a FWHM of about 70ps , the average power being 1to2mW . Light is injected into the sample and collected using 600-μ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-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 .

Fig. 3

A schematic of the instrumentation in interstitial mode.




In order to compare diffusion-based and WMC-based evaluation of experimental data, measurements were performed on tissue phantoms based on Intralipid (Fresenius Kabi 200mgml ) 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 μ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 577ml water and 22ml Intralipid. Measurements were performed at 2to8ml total volume of added ink solution ( 1-ml increments). The added scatterer series was performed on a phantom based on 577ml water and 4ml ink solution. Measurements were performed at 10to24ml added volume of Intralipid (in 2-ml increments). A cylindrical plastic container, Ø=110-mm and height=70mm , 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-mm depth in the center of the phantom, separated 14.7mm , center to center. Data acquisition was performed for 30s for each measurement, and the IRF was measured before, after, and occasionally between the measurements (monitoring potential temporal drifts in the system).



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 μa , μs , and an amplitude factor, k , are varied in order to minimise a χ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

where ρ is the source-detector separation, and c is the speed of light in the sample. To reach this expression, we use the absorption-independent definition of the diffusion coefficient, D=(3μs)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 μs . That is, for each value in a set of trial values of μs , the optimal choice of μa and k is determined in a suboptimization procedure. The suboptimized error norm, χ̃2(μs) , is given in Eq. 4:

Eq. 4

and is thus to be exhaustively searched on some appropriate μs grid. In this work, χ̃2(μs) is evaluated in steps of 0.05cm1 . An example of this procedure is given in Fig. 4 .

Fig. 4

WMC fitting procedure of the experimental data featuring the 2-ml added ink absorption series measurement at 660nm , ρ=14.7mm .


This article uses an MC database for infinite geometries (interstitial measurements) consisting of about 8×107 detection events, originating from a simulation of 2×108 photons. The database for semi-infinite evaluations is of a similar size, resulting from a simulation of 109 photons. Both these databases are based on simulations where g=0.7 , n=1.33 , NA=0.29 , μsmax=90cm1 , and tmax=2ns . 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 15s .

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.




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 μs=5 , 10cm1 and μa=0.1 , 0.5cm1 , 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±4% in absorption, and 0±2% for reduced scattering (mean± . 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 15mm , μs=5cm1 , and μa=0.1cm1 .



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 650to900nm ).

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-mm source-detector separation, with μs=7.5cm1 , and μa=0.5cm1 .


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-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-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, μ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

where μ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 μa and μs .

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) μs is kept constant, while μa gradually increases (added absorber series), and (2) μa is kept constant, while μs 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 μa is steeper than the true increase, and that the constancy of μs is replaced by a fairly linear increase. For the added scatterer series at μa=0.5cm1 , note the offset-like pattern for μs , as well as the remarkably poor agreement between true and derived μa at low values of μs .

Fig. 9

Derived optical properties for two levels of μs when diffusion modeling is used to interpret WMC data (corresponds to the added absorber series).


Fig. 10

Derived optical properties for two levels of μa when diffusion modeling is used to interpret WMC data (corresponds to the added scatterer series).



Experimental Results

The experimental results from the added absorber series, evaluated using both diffusion theory and WMC, are presented in Fig. 11 . Here, derived μ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 μs , 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

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, μ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 μs , 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 μs 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 WorkVan Staveren 26



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.5mm from the source-detector separation ρ . In contrast, our WMC simulations assumed that the source and detector fibers had a radius of 0.3mm , 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 μa and μs . However, light propagation in the photon migration regime is known to fall under similarity as long as μs=(1g)μ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 μs , 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 (2.25to4.75mm) 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 μa and a constant level of μs . Extrapolating measured μ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 μs together with a constant level of absorption. The small offsets in derived μs 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 (801) . In this work, we mainly use a 5020 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 801 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μt=1(μa+μs) , whereas a WMC approach uses a 1μs average step-size. This means that the photon paths generated by the two methods will differ if μa is nonzero. Since traditional Monte Carlo for light propagation in turbid media assumes that scattering dominates absorption ( μsμ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=μs(μa+μ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 a1 or a0 , but notes that equivalence is not guaranteed in the intermediate region where μsμ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 tmax may impose errors during steady-state analysis of generated WMC databases.



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 μs and μa ). The relevance of this regime is, for example, motivated by the recent interest in the optical properties of human prostate tissue, where μa> 0.3cm1 and μs<10cm1 .


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



M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical-properties,” Appl. Opt., 28 (12), 2331 –2336 (1989). 0003-6935 Google Scholar


S. T. Flock, M. S. Patterson, B. C. Wilson, and D. R. Wyman, “Monte Carlo modeling of light-propagation in highly scattering tissues: I. model predictions and comparison with diffusion-theory,” IEEE Trans. Biomed. Eng., 36 (12), 1162 –1168 (1989). 0018-9294 Google Scholar


A. H. Hielscher, S. L. Jacques, L. H. Wang, and F. K. Tittel, “The influence of boundary-conditions on the accuracy of diffusion-theory in time-resolved reflectance spectroscopy of biological tissues,” Phys. Med. Biol., 40 (11), 1957 –1975 (1995). 0031-9155 Google Scholar


M. L. Pantelides, C. Whitehurst, J. V. Moore, T. A. King, and N. J. Blacklock, “Photodynamic therapy for localized prostatic cancer—light penetration in the human prostate-gland,” J. Urol. (Baltimore), 143 (2), 398 –401 (1990). 0022-5347 Google Scholar


T. Svensson, M. Einarsdóttír, K. Svanberg, and S. Andersson-Engels, “In vivo optical characterization of human prostatic tissue using near-infrared time-resolved spectroscopy,” J. Biomed. Opt., 12 (1), 014022 (2007). 1083-3668 Google Scholar


B. C. Wilson and G. Adam, “A Monte Carlo model for the absorption and flux distributions of light in tissue,” Med. Phys., 10 (6), 824 –830 (1983). 0094-2405 Google Scholar


L. H. Wang, S. L. Jacques, and L. Q. Zheng, “MCML Monte Carlo modeling of light transport in multilayered tissues,” Comput. Methods Programs Biomed., 47 (2), 131 –146 (1995). 0169-2607 Google Scholar


J. Swartling, A. Pifferi, A. M. K. Enejder, and S. Andersson-Engels, “Accelerated Monte Carlo models to simulate fluorescence spectra from layered tissues,” J. Opt. Soc. Am. A, 20 (4), 714 –727 (2003). 0740-3232 Google Scholar


N. Everall, T. Hahn, P. Matousek, A. W. Parker, and M. Towrie, “Photon migration in Raman spectroscopy,” Appl. Spectrosc., 58 (5), 591 –597 (2004). 0003-7028 Google Scholar


A. Kienle and M. S. Patterson, “Determination of the optical properties of turbid media from a single Monte Carlo simulation,” Phys. Med. Biol., 41 (10), 2221 –2227 (1996). 0031-9155 Google Scholar


A. Pifferi, R. Berg, P. Taroni, and S. Andersson-Engels, “Fitting of time-resolved reflectance curves with a Monte Carlo model,” Trends in Optics and Photonics: Advances in Optical Imaging and Photon Migration, 1996). Google Scholar


A. Pifferi, P. Taroni, G. Valentini, and S. Andersson-Engels, “Real-time method for fitting time-resolved reflectance and transmittance measurements with a Monte Carlo model,” Appl. Opt., 37 (13), 2774 –2780 (1998). 0003-6935 Google Scholar


R. Graaff, M. H. Koelink, F. F. M. Demul, W. G. Zijlstra, A. C. M. Dassel, and J. G. Aarnoudse, “Condensed Monte Carlo simulations for the description of light transport,” Appl. Opt., 32 (4), 426 –434 (1993). 0003-6935 Google Scholar


J. Swartling, “Biomedical and atmospheric applications of optical spectroscopy in scattering media,” Lund University, (2002). Google Scholar


H. P. Xu, T. J. Farrell, and M. S. Patterson, “Investigation of light propagation models to determine the optical properties of tissue from interstitial frequency domain fluence measurements,” J. Biomed. Opt., 11 (4), 041104 (2006). 1083-3668 Google Scholar


S. A. Prahl, M. Keijzer, S. L. Jacques, and A. J. Welch, “A Monte Carlo model of light propagation in tissue,” Dosimetry of Laser Radiation in Medicine and Biology, IS 5 102 –111 SPIE, Bellingham, WA (1989). Google Scholar


M. Saito and M. Matsumoto, “SIMD-oriented fast Mersenne twister: a 128-bit pseudorandom number generator,” (2006). Google Scholar


S. Lacey and R. Box, “A fast, easy sort,” BYTE, 16 (4), 315 (1991). 0360-5280 Google Scholar


L. Wang, S. L. Jacques, and L. Zheng, “CONV-convolution for responses to a finite diameter photon beam incident on multi-layered tissues,” Comput. Methods Programs Biomed., 54 141 –150 (1997). 0169-2607 Google Scholar


S. A. Prahl, “Light transport in tissue,” Univ. Texas Austin, (1988). Google Scholar


T. Svensson, J. Swartling, P. Taroni, A. Torricelli, P. Lindblom, C. Ingvar, and S. Andersson-Engels, “Characterization of normal breast tissue heterogeneity using time-resolved near-infrared spectroscopy,” Phys. Med. Biol., 50 (11), 2559 –2571 (2005). 0031-9155 Google Scholar


A. Pifferi, A. Torricelli, “Performance assessment of photon migration instruments: the MEDPHOT protocol,” Appl. Opt., 44 (11), 2104 –2114 (2005). 0003-6935 Google Scholar


F. E. W. Schmidt, M. E. Fry, E. M. C. Hillman, J. C. Hebden, and D. T. Delpy, “A 32-channel time-resolved instrument for medical optical tomography,” Rev. Sci. Instrum., 71 (1), 256 –265 (2000). 0034-6748 Google Scholar


B. W. Pogue and M. S. Patterson, “Review of tissue simulating phantoms for optical spectroscopy, imaging, and dosimetry,” J. Biomed. Opt., 11 (4), 041102 (2006). 1083-3668 Google Scholar


S. J. Madsen, M. S. Patterson, and B. C. Wilson, “The use of india ink as an optical absorber in tissue-simulating phantoms,” Phys. Med. Biol., 37 (4), 985 –993 (1992). 0031-9155 Google Scholar


H. J. van Staveren, C. J. M. Moes, J. van Marle, S. A. Prahl, and M. J. C. Vangemert, “Light scattering in intralipid-10-percent in the wavelength range of 4001100nm,” Appl. Opt., 30 (31), 4507 –4514 (1991). 0003-6935 Google Scholar


S. T. Flock, S. L. Jacques, B. C. Wilson, W. M. Star, and M. J. C. Vangemert, “Optical properties of intralipida phantom medium for light-propagation studies,” Lasers Surg. Med., 12 (5), 510 –519 (1992). 0196-8092 Google Scholar


C. Abrahamsson, T. Svensson, S. Svanberg, S. Andersson-Engels, J. Johansson, and S. Folestad, “Time and wavelength resolved spectroscopy of turbid media using light continuum generated in a crystal fiber,” Opt. Express, 12 (17), 4103 –4112 (2004). 1094-4087 Google Scholar


R. Cubeddu, A. Pifferi, P. Taroni, A. Torricelli, and G. Valentini, “Experimental test of theoretical models for time-resolved reflectance,” Med. Phys., 23 (9), 1625 –1633 (1996). 0094-2405 Google Scholar


R. Graaff, J. G. Aarnoudse, F. F. M. Demul, and H. W. Jentink, “Similarity relations for anisotropic scattering in absorbing media,” Opt. Eng., 32 (2), 244 –252 (1993). 0091-3286 Google Scholar


S. L. Jaques and L. Wang, “Monte Carlo modeling of light transport in tissues,” Optical-Thermal Response of Laser-Irradiated Tissue, 73 –100 Plenum Press, New York (1995). Google Scholar


B. Farina, S. Saponaro, E. Pignoli, S. Tomatis, and R. Marchesini, “Monte Carlo simulation of light fluence in tissue in a cylindrical diffusing fiber geometry,” Phys. Med. Biol., 44 1 –11 (1999). 0031-9155 Google Scholar
©(2008) Society of Photo-Optical Instrumentation Engineers (SPIE)
Erik Alerstam, Stefan Andersson-Engels, and Tomas Svensson "White Monte Carlo for time-resolved photon migration," Journal of Biomedical Optics 13(4), 041304 (1 July 2008).
Published: 1 July 2008

Back to Top