Abstract

Monte-Carlo simulation is an important tool in the field of biomedical optics, but suffers from significant computational expense. In this paper, we present the multicanonical Monte-Carlo (MMC) method for improving the efficiency of classical Monte Carlo simulations of light propagation in biological media. The MMC is an adaptive importance sampling technique that iteratively equilibrates at the optimal importance distribution with little (if any) a priori knowledge of how to choose and bias the importance proposal distribution. We illustrate the efficiency of this method by evaluating the probability density function (pdf) for the radial distance of photons exiting from a semi-infinite homogeneous tissue as well as the pdf for the maximum penetration depth of photons propagating in an inhomogeneous tissue. The results agree very well with diffusion theory as well as classical Monte-Carlo simulations. A six to sevenfold improvement in computational time is achieved by the MMC algorithm in calculating pdf values as low as 10-8. This result suggests that the MMC method can be useful in efficiently studying numerous applications of light propagation in complex biological media where the remitted photon yield is low.

© 2005 Optical Society of America

1. Introduction

Numerical studies of light propagation in complex optical media such as biological tissue [1]-[5] are often performed through the Monte-Carlo simulation approach due to lack of analytic solutions to the radiative transfer equation [6]. Although Monte-Carlo modeling has a broad generality and can cope with arbitrary media complexity, it suffers from poor efficiency; substantial computational effort is required in order to obtain acceptable statistical precision. The poor efficiency of classical Monte-Carlo (CMC) simulations is particularly detrimental in problems involving back-reflected light from biological media having multiple-scattering characteristics but low backscattering. Approaches to combat this drawback involve the use of variance reduction techniques which reduce the estimator variance of the scored physical quantity, thus allowing a given statistical precision to be achieved using fewer light propagation paths. Among these methods are implicit photon capture [5] and importance sampling (IS) [2].

In the implicit photon capture technique, many photons (i.e., a packet of photons) rather than single photons are propagated simultaneously along the tissue, therefore improving the efficiency of the CMC simulation by eliminating the all-or-none nature of the single photon absorption events [5]. For the importance sampling method, the trajectories of photons propagating in the random medium are sampled from an importance proposal distribution, which tends to coerce the event of interest to occur more frequently. Typical events include the remitted radial distance of the photon from the laser source or the location from which the photon emerges given that it has escaped from the tissue surface. The fundamental difficulty with IS is the design of a “good” importance proposal distribution. Choosing a “good” distribution can result in significant run-time savings; the penalty for selecting a “bad” distribution can be longer run times than for the CMC simulation. Common “good” importance proposal distributions which force incident photons to penetrate deeper into the medium, while simultaneously enhancing backscattering are based on physical intuition and involve biasing the scattering lengths and the polar angles at certain steps throughout the simulation [2]. However, as the complexity of the light transport model in the biological medium increases, it becomes harder to design suitable photon trajectory distributions a priori.

In this paper, we apply the multicanonical Monte-Carlo (MMC) technique [7] to simulations of light propagation in biological media. The MMC method was originally designed to study first-order phase transitions [7], and it has been recently applied to compute also polarization mode dispersion statistics [8], bit error rates in optical communication systems [9] and noise statistics of semiconductor optical amplifiers [10]. Like ordinary IS, the MMC algorithm increases the number of rare events of interest by sampling the photon trajectories from an importance proposal distribution. However, the advantage of MMC is that it approaches the optimal importance distribution iteratively with little (if any) a priori knowledge of how to bias. The sampling from the sub-optimal importance distribution used within a specific MMC iteration is performed through random walks. To demonstrate the power of the MMC simulation method, we highlight simulations for which the CMC techniques are computationally inefficient. For brevity in the remainder of the manuscript, we will designate the acronym CMC to represent the classical Monte-Carlo simulations that propagate photon packets rather than single photons. Our results include estimating the rare event probabilities of the radial distance of photons exiting a semi-infinite homogeneous tissue as well as the rare event probabilities of the maximum penetration depth of photons propagating in an inhomogeneous tissue. The improved efficiency of the MMC algorithm suggests that multicanonical sampling will be useful in more complex optical forward and inverse problems for which a large space of feasible photon trajectories is to be explored and measures of resolution and adequate statistical precision of the solutions are desirable.

This paper is organized as follows. In Section 2, we outline the theoretical framework of the MMC algorithm and in Section 3 we describe its applicability to problems related to light propagation in biological media. Section 4 is devoted to the simulation results and discussion. Finally, conclusions are drawn is Section 5.

2. Multicanonical Monte-Carlo (MMC) algorithm

The MMC method [7] is a numerical technique that like IS algorithms increases the number of samples in the tail region of the computed probability density function (pdf) by sampling more frequently the values that have higher impact on the parameter being estimated. This biased sampling of the important values reduces the variance of the calculated pdf estimator. However, in contrast to standard IS where one has to intuitively choose an appropriate bias that will encourage sampling of the important values, the MMC algorithm uses an adaptive iterative procedure to select this bias with little (if any) a priori knowledge of how to bias. The iterative procedure updates the bias for the subsequent iteration by drawing samples trough a random walk which is performed in the important regions of the state space. As the iterations evolve the bias histogram tends to equilibrate at the optimal bias. Next, we describe the mathematical framework of the MMC method.

Let fX (x) represent a given pdf defined on a high-dimensional space, χ, that describes the set of possible configurations of a system. In light propagation simulations this set corresponds to the possible trajectories of photons traveling in the random medium. Also, let the random variable Y(x) describe the scored physical quantity and let Li represent the event {yi < Y < yi + dy }. Then, the probability of the event Li is simply

fY(yi)dy=χ1Li(x)fX(x)dx

where fY (yi ) is the pdf of Y evaluated at yi and 1Li (x) is the indicator function of the event Li

1Li(x)={1xLi0xLi

Hence, fY (yi ) can be estimated using the CMC estimator which reads as

f̂YCMC(yi)=1dy·NΣn=1N1Li(xn),iN

where {xn}n=1N are N independent identically distributed (iid) samples drawn from fX (x). Eq. (3) indicates that the CMC estimator simply computes the percentage of samples (out of the N drawn samples) corresponding to the event Li . In order to reduce the variance of f^YCMC (yi ), the MMC method employs the importance sampling (IS) estimation principle [11], [12] and introduces an importance proposal distribution, q(x). This allows Eq. (1) to be written as

fY(yi)dy=χ1Li(x)w(x)q(x)dx

where w(x)= fX(x) / q(x) is the so called importance weight [11], [12] and the corresponding IS estimator of fY(yi) is expressed as

f̂YIS(yi)=1dy·NΣn=1N1Li(xn)w(xn),iN

with xn being generated now from q(x). Notice that the IS estimator use the importance weight w(xn) in order to correct for the use of samples xn drawn from the biased distribution q(x) -rather than from the original distribution fX(x) - therefore ensuring that the pdf estimation is unbiased alike the CMC estimator. Also, note that Eqs. (5) and (3) are identical when q(x)= fX(x) and that the optimal importance distribution which minimizes the variance of the IS estimator reads as [11], [12]

qopt(x)=1Li(x)fX(x)χ1Li(x)fX(x)dx=1Li(x)fX(x)fY(yi)dy,iN

However, the optimal importance distribution is not very useful since it depends on fY(yi), i∈ N which is unknown. Hence, while the aim in IS is to design a “good” approximation of q opt (x) based on physical intuition, the idea behind the MMC algorithm is to iterate over sub-optimal importance proposal distributions that converge toward q opt (x). Specifically, the sub-optimal importance proposal distribution used in the k-th iteration of the MMC method, q opt(k) (x) is given by

qopt(k)(x)1Li(x)fX(x)fY(k1)(yi),iN

where fY(k1) (yi) is the estimation of the pdf of Y obtained in the previous iteration. The sub-optimal importance proposal distribution q opt(k) (x) tends to equilibrate at the optimal bias q opt (x) as the iteration number increases. fY(0) (yi), i∈ N can be a uniform distribution or can be calculated by the CMC estimator (Eq. (3)) using a small number of samples. Within the k-th MMC iteration, the Markov chain Monte Carlo (MCMC) strategy [12] and specifically the Metropolis algorithm [13] is employed to generate samples xn from q opt(k) (x). A Metropolis step involves proposing a candidate sample xn*= xn + δx where xn is the current sample and δx is sampled from a symmetric distribution (e.g., uniform distribution with zero mean). The candidate sample is then accepted as the next sample with probability min{1, q opt(k) (xn*) / q opt(k) (xn)=[fX(xn*) fY(k1) (y(xn))] / [fX(xn) fY(k1) (y(xn*))]} otherwise the next sample remains xn. Notice that by using the Metropolis method we only need to know q opt(k) (x) up to a constant of proportionality. Essentially, the Metropolis algorithm generates samples from q opt(k) (x) by exploring the state space χ using a random walk, or equivalently a homogeneous, irreducible aperiodic Markov chain which employs the detailed balance condition [12] to converge toward the invariant distribution q opt(k) (x) independently of the starting state. Consequently, the success or failure of the algorithm hinges on the choice of distribution of δx as well as on the number of Metropolis steps (N) and the number of MMC iterations that are simulated. Note that the selection of these parameters is performed empirically such that each MMC iteration covers a larger values range of the simulated random variable with adequate standard deviation (including those values located in the pdf’s tail) and that the last MMC iteration achieves the desired probability level. Finally, the updated estimation of the pdf of Y evaluated at yi , i∈ N which results from the k-th MMC iteration is computed using Eq. (5) and reads as

f̂Y(k)(yi)={fY(k1)(yi)·Σn=1N1Li(xn)ifatleastonesamplexnLifY(k1)(yi)otherwise

with C being a normalizing constant ensuring that ∫-∞ fY(k) (y) dy = 1.

To conclude, the MMC technique uses the IS principle (Eq. (5)) together with sub-optimal importance proposal distributions that are computed iteratively (Eq. (7)) to approach the optimal importance proposal distribution (Eq. (6)). An N step Metropolis algorithm is constructed to mimic samples drawn from these sub-optimal importance proposal distributions.

3. Application of the MMC method to light propagation in biological media

The set of possible photon packet trajectories in a biological medium is determined by the polar angles m}, the azimuthal angles m}, the scattering lengths {lm}, and the reflection events {rm} at the boundaries of layers having distinct refractive indices. m = 0,1, ⋯,s, where s indicates the index of the last scattering event. Let fΘm),fΦm), fL(lm) and fR(rm) represent the pdf’s of the iid random variables θm , ϕm , lm and rm , hence fX (x) = m=0s f Θ(θm )f Φ(ϕm )fL (lm )fR (rm )where x = (θ1 , ⋯,θs ,ϕ1 , ⋯,ϕs ,l1 , ⋯, ls , r1 , ⋯, rs , and χ=(Θ,Φ,L,R). Also, let y(x) denote the event of interest for which we wish to compute the pdf. Without loss of generality, we concentrate on events which correspond to photons exiting the medium. Therefore, the computed pdf’s are essentially conditional pdf’s of y(x) given that the photon packet exits the medium. Similar principles can be applied for the events related to photons being absorbed along their propagation in the medium. Finally, notice that we denote conditional pdf’s and pdf’s similarly throughout this work in order to simplify notations and the correct interpretation of the pdf is rather understood from the context.

Simulations of photon packet trajectories in biological media are more efficient than single photon trajectories [4], [5]. Hence, we first employ the MMC algorithm to estimate f no_abs Y,Wpacket (y,w packet) (that is,f no_abs Y,Wpacket|EXIT(y,w packet | EXIT)), which is the joint pdf of y(x) without absorption and the corresponding photon packet weight, w packet (x), given that the photon packet exists the medium, and then we use the relationship

fY(y)=01wpacket·fY,Wpacketno_abs(y,wpacket)dwpacket/01wpacket·fY,W,packetno_abs(y,wpacket)dwpacketdy

to compute the pdf of y (x) with absorption. This approach significantly contributes to the overall efficiency of the MMC method and is preferable over MMC simulations of single photon trajectories which directly estimate fy(y).

In view of the fact that fX(x) is a separable function, each Metropolis step within the k-th MMC iteration is implemented as follows. We first perturb separately θm , ϕm , lm and rm , m=0,1, ⋯,s by ~ U(-εθ / 2,εθ /2), ~ U(-εϕ / 2,εϕ /2), dl ~ U(-εl / 2,εl /2) and dr ~ U(-εr /2,εr/2), respectively, with U symbolizing uniform distribution and εθ , εϕ , εl , εr being selected such that most of the estimated pdf modes are visited equivalently. Perturbation coefficients leading to an acceptance ratio - defined as the ratio of the accepted Metropolis steps to the total number of Metropolis steps (N) - of 30%-45% were found to be appropriate in our studies. Next, we accept or reject independently each of the proposed θm* = θm +, ϕm* = ϕm + , lm* = lm + dl and rm* = rm + dr, m=0,1,⋯,s with probability min{1, fΘ(θm*) / fΘm)}, min{1,fΦ(ϕm*) / fΦ (ϕm)}, min{1, fL(lm*) /fL (lm)}, and min{1, fR(rm*) / fR (rm)}, respectively. The fact that ϕm and rm are uniformly distributed and are subjected to periodic boundary conditions results in accepting ϕm* and rm* with probability one. Following the selection of the test trajectory x * , y (x*) is then calculated by propagating the photon packet along the biological media as described by Wang et al. [5]. Finally, the entire Metropolis step from x to x * is accepted with probability min{1,f no_abs(k-1) Y,Wpacket (y(x),w packet (x)) /f no_abs(k-1) Y,Wpacket (y(x*),w packet (x*))}. Note that the Metropolis algorithm is designed to accept only steps which result in photon packets that exit the medium, hence estimating the conditional pdff no_abs Y,Wpacket|EXIT (y,w packet | exit).

It is worth mentioning that the described perturbation scheme explores the high dimensional space χ more efficiently than the method presented in section 2 due to the independent perturbations of the pdf‘s comprising fX(x). However, notice that this technique is applicable only in cases for which fX(x) is a separable function.

Following N Metropolis steps, the joint pdf f Y,wPacket(y,w packe t) is updated similarly to Eq. (8):

f̂Y,Wpacketno_abs(k)(yi,wpacketj)={fY,Wpacketno_abs(k)(yiwpacketj)·Σn=1N1Lij(xn)ifatleastonesamplexnLijfY,Wpacketno_abs(k)(yiwpacketj)otherwise

where Lij represents the event {yi < Y < yi +dy, w packetj < W packet < W packetj+dw packet} and C stands for a normalizing constant ensuring that 01 no_abs(k) Y,Wpacket(yi , w packet) dw packet dy = 1. Finally, the estimated fY(y) obtained from the k-th MMC iteration is

f̂Y(k)(y)=01wpacket·f̂Y,Wpacketno_abs(k)(y,wpacket)dwpacket/01wpacket·f̂Y,W,packetno_abs(k)(y,wpacket)dwpacketdy

Figure 1 summarizes the pseudo code of the MMC algorithm applied in computations of pdf’s emerging in simulations of light transport in biological media using the Heyney-Greenstein pdf for the cosine of the deflection angle [5].

 

Fig. 1. Multicanonical Monte-Carlo (MMC) algorithm for light transport in biological media.

Download Full Size | PPT Slide | PDF

4. Simulation results and discussion

Two scenarios of light propagation in a biological medium are examined in this section. These scenarios serve as demonstrations for the potential use of the multicanonical sampling in stochastic simulations of cases for which CMC methods are impractical or computationally inefficient.

The results presented hereon were simulated using Matlab software running on a conventional Pentium 4 desktop computer.

4.1 Radially resolved steady state diffuse reflectance pdf for photons propagating in a semi-infinite homogeneous random medium

To illustrate the efficacy of the MMC method compared with that of CMC simulations, we first used multicanonical sampling to compute the pdf of the radial distance of photons exiting a semi-infinite homogeneous tissue with the following optical properties: μa = 0.25 mm-1, μs = 50 mm-1 and g = 0.9. We assumed that the tissue-ambient medium boundary is refractive-index matched. In the simulations, the photon packets were launched orthogonally onto the tissue at the origin. Tracking a photon packet was terminated either when the photon packet exited the medium, the weight of the photon packet was sufficiently low (< 10-15) or when the number of scattering events exceeded a predefined value. For this study, the maximum number of scattering events was 2000, resulting in a state space χ of dimension (1999 polar angles) × (1999 azimuthal angles) × (2000 scattering lengths) and ensuring unbiased low probability values which were further confirmed by trial simulations using 10000 scattering events. Furthermore, the joint pdf of the radial distance without absorption and the corresponding photon packet weight is scored using a linear scale for the radial distance axis and a semi-logarithmic scale for the packet weight axis.

The MMC simulations employed six iterations which were found to be sufficient in order to estimate diffuse reflectance probability values as low as ~10-8 with standard deviation of one half-decade. The first iteration consisted of 2187 photon packets and used a CMC simulation which introduced the first coarse estimation for the diffuse reflectance pdf. Note that alternatively one can set f no_abs(0) Y,Wpacket (y(x),w packet (x)) = 1/D with D = [0, ymax ] × [0,1] and execute an MMC iteration. These two approaches for computing the first coarse estimation of the diffuse reflectance pdf showed similar results. The remaining iterations comprised five multicanonical sampling procedures with 15000, 22500, 33750, 50625 and 75938 photon packets. Accordingly, a total number of 2·105 photon packets was simulated during 80 minutes. Notice that the number of photon packets used in each iteration was increased in order to encourage the Metropolis random walk to accept samples also in the tail of the estimated conditional pdf and not only at shorter radial exit distances. Furthermore, the number of photon packets in each iteration was determined such that it yielded an improvement of approximately one and a half orders of magnitude in the computed diffuse reflectance pdf values.

Figure 2 shows four (out of six) of the MMC iterations. It includes also the analytic expression of the diffuse reflectance pdf [14] in dashed line. We observe that the estimation of the diffuse reflectance pdf resulting from the first MMC iteration is limited to values higher than 10-3 as indicated by the straight horizontal line appearing at the pdf tail. This line places a threshold of one half-decade for the fluctuations of the estimated pdf. Note that for the tissue optical parameters used in this example, pdf values higher than 10-3 correspond to radial exit distances lower than 2 mm. The subsequent MMC iterations cover a larger radial exit distance range which finally extends over approximately 8 mm. Notice that the last MMC iteration considerably reduces the estimation variance of the preceding MMC iterations and generates a pdf estimator which successfully follows the tail of the theoretical diffuse reflectance pdf down to values of ~10-8. To conclude, Fig. 2 illustrates the usefulness of multicanonical sampling in increasing the dynamic range of the estimated pdf region with adequate statistical accuracy. Moreover, each MMC iteration improves the tail estimation in approximately one order of magnitude while maintaining a similar order of simulated photon packets number.

 

Fig. 2. MMC estimators of the diffuse reflectance pdf of a semi-infinite homogeneous random medium.

Download Full Size | PPT Slide | PDF

To compare the performances of the MMC and CMC techniques, we performed a CMC simulation employing 2·105 photon packets - a number identical to that simulated previously by the MMC method. Figure 3 presents the estimated diffuse reflectance pdf obtained from the CMC simulation (in red color) and MMC simulation (in green color) using an identical number of photon packets. Once more, the dashed line describes the theoretical pdf calculated by diffusion theory [14], and the straight horizontal lines represent the minimum achievable pdf value for which the fluctuations of the estimated pdf have not exceed a threshold of one half-decade. Clearly, the MMC estimator of the diffuse reflectance pdf is significantly superior to the CMC estimator - it covers a larger range of radial exit distances and approaches pdf values that are more than two orders of magnitude smaller while preserving an acceptable estimation variance over the entire radial exit distance range.

It is important to point out that the CMC simulation ran 2.5 times faster than the MMC simulation, even though both simulation codes employed an identical photon transport function and had an identical number of photon packets. This difference in the run time stemmed mainly from the fact that longer photon trajectories should be computed throughout the MMC iterations in order to achieve lower values of the diffuse reflectance pdf. This point is further confirmed by executing a CMC simulation using 3.9·106 photon packets as shown in Fig. 3 in blue color. While both CMC and MMC techniques attained similar pdf values as low as 10-8 with similar statistical accuracy, the CMC simulation required nine hours to complete - 6.75 times slower than the MMC algorithm. These results illustrate the improved performance of the MMC algorithm in estimating adequately the diffuse reflectance pdf over a larger range of radial exit distances.

 

Fig. 3. MMC and CMC estimators of the diffuse reflectance pdf of a semi-infinite homogeneous random medium.

Download Full Size | PPT Slide | PDF

4.2 pdf of the maximum penetration depth of photons that propagate and exit a two-layered random medium

To further explore the performance of multicanonical sampling of photon packets propagating in biological media, we performed MMC simulations for a two-layered random medium. In these simulations, a pencil beam of photons entered the top layer of the medium along the z direction and traveled inside the inhomogeneous media. Then, the maximum penetration depth Zmax of photon packets that successfully exited the top layer of the medium (i.e., Zmax | EXIT) was recorded and its pdf was computed. Tracking a photon packet was terminated either when the weight of the photon packet was sufficiently low (< 10-15) or when the number of scattering events exceeded a predefined value. For this study, the maximum number of scattering events was 3000 (therefore guaranteeing unbiased pdf estimators) and a semi-logarithmic grid was used for the packet weight axis of the joint pdf of Zmax without absorption and the corresponding photon packet weight.

The scattering properties of the two layers were identical, that is μs1 = μs2 = 50 mm-1 and g1 = g2 = 0.9, whereas the absorption coefficients of the first and second layer were μa1 = 0.25 mm-1 and μa2 = 0.025 mm-1, respectively. Furthermore, the top-layer thickness was 3.75 mm and the refractive indices for the two layers as well as for the overlying ambient medium were equal.

First, we confirmed that the MMC estimation for the diffuse reflectance pdf is consistent with the analytic result [15]. Next, six MMC iterations were executed during 80 min using a total number of 2·105 photon packets in order to compute the pdf of the random variable Zmax | EXIT. Figure 4 presents the last iteration of the MMC pdf estimation for Zmax | EXIT (in green color). The straight horizontal line at the pdf tail shows that pdf values lower than 10-7 (for which the fluctuations of the estimated pdf have not exceed a threshold of one half-decade) were obtained throughout the simulation. Furthermore, two distinct slopes are clearly observed, where the slope for 0 < Zmax < 3.75 mm is larger than that for Zmax > 3.75 mm. This result is expected and stems from the lower absorption coefficient of the second layer.

 

Fig. 4. MMC and CMC estimators of the pdf of the maximum penetration depth of light in a two-layered inhomogeneous random medium.

Download Full Size | PPT Slide | PDF

We compared the performance of the multicanonical and Monte-Carlo sampling methods by carrying out two CMC simulations. The first CMC simulation lasted 30 minutes and estimated the pdf of the maximum penetration depth using 2·105 photon packets (Fig. 4, red dashed line). The second CMC simulation employed 3.75·106 photon packets and was completed in eight hours (Fig. 4, in blue line). We note that with the CMC method, simulating 2·105 photon packets is not adequate to predict the correct pdf for the Zmax | EXIT. Moreover, the number of photon packets should be increased by at least one order of magnitude in order to obtain an acceptable estimation of this pdf. However, the improved estimation variance comes at the expense of computational time which in this case is six-fold slower than that achieved by multicanonical sampling for a similar level of estimation variance. This study therefore demonstrates the superior efficacy of the MMC algorithm over the CMC method in sampling back reflected light in order to successfully identify tissue inhomogeneities.

5. Conclusions

In propagation simulations, where importance proposal distributions can be intuitively determined a priori, importance sampling (IS) [2] will be highly efficient, since the photon trajectories that give rise to the tail-states are predetermined. In this paper, we have presented the application of an adaptive importance sampling technique based on multicanonical Monte-Carlo (MMC) simulations to the studies of light propagation in biological media. Like IS, MMC efficiently accesses the low-probability tails of numerous pdf’s of interest in biological light transport simulations. Unlike standard IS algorithms, however, the MMC method is blind to the physics of the light propagation model. As a result, MMC may be considered more general than IS allowing it to be applied to a wider variety of simulation geometries in which the importance proposal distributions can not be intuitively advised. It is important to point out that the efficiency of MMC sampling hinges on the choice of the perturbation coefficients used in the Metropolis random walk as well as on the number of the simulated Metropolis steps and MMC iterations. These parameters are selected empirically such that (1) most of the estimated pdf modes are visited equivalently through the random walk, (2) the desired probability level is obtained with adequate standard deviation, and (3) each MMC iteration successfully increases the dynamic range in the pdf tail region.

We have shown the efficacy of the MMC method by evaluating the pdf for the radial distance of photons exiting from a semi-infinite homogeneous medium as well as the pdf for the maximum penetration depth of photons propagating in an inhomogeneous medium. The computed pdf’s were consistent with those estimated by classical Monte-Carlo (CMC) simulations. Furthermore, in these examples a six- to seven-fold improvement in computational time was achieved using multicanonical sampling. These results illuminate the potential usefulness of the MMC algorithm for an efficient sampling of pdf’s in complex biological simulations. Based on these encouraging results we feel that further investigation of the computational efficiency of MMC for other simulation parameters_as well as extension of MMC to model temporal, polarization, and coherence effects are merited.

Acknowledgments

This research was supported in part by the DoD Medical Free Electron Laser Program F49620-021-1-1-0014.

References and links

1 . S. T. Flock , M. S. Patterson , B. C. Wilson , and D. R. Wyman , “ Monte Carlo modeling of light propagation in highly scattering tissues - 1. Model predictions and comparison with diffusion theory ,” IEEE Trans. Biomed. Eng. 36 , 1162 – 1167 ( 1989 ). [CrossRef]   [PubMed]  

2 . J. M. Schmitt and K. Ben-Letaief , “ Efficient Monte Carlo simulation of confocal microscopy in biological tissue ,” J. Opt. Soc. Am. A 13 , 952 – 961 ( 1996 ). [CrossRef]  

3 . G. Yao and L. Wang , “ Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media ,” Phys. Med. Biol. 44 , 2307 – 2320 ( 1999 ). [CrossRef]   [PubMed]  

4 . D. A. Boas , J. P. Culver , J. J. Stott , and A. K. Dunn , “ Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head ,” Opt. Express 10 , 159 – 170 ( 2002 ), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-3-159 [PubMed]  

5 . L. Wang , S. L. Jacques , and L. Zheng , “ MCML - Monte Carlo modeling of light transport in multi-layered tissues ,” Comput. Methods and Programs in Biomed. 47 , 131 – 146 ( 1995 ). [CrossRef]  

6 . A. Ishimaru Wave Propagation and Scattering in Random Media ( Academic Press, Inc., San Diego 1978 ).

7 . B. A. Berg and T. Neuhaus , “ Multicanonical ensemble: A new approach to simulate first-order phase transitions ”, Phys. Rev. Lett. 68 , 9 – 12 ( 1992 ). [CrossRef]   [PubMed]  

8 . D. Yevick , “ Multicanonical communication system modeling - application to PMD statistics ,” IEEE Photon. Tech. Lett. 14 , 1512 – 1514 ( 2002 ). [CrossRef]  

9 . R. HolzlÖhner and C. R. Menyuk , “ Use of multicanonical Monte Carlo simulations to obtain accurate bit error rates in optical communications systems ,” Opt. Lett. 28 , 1894 – 1896 ( 2003 ). [CrossRef]   [PubMed]  

10 . A. Bilenca and G. Eisenstein , “ Statistical noise properties of an optical pulse propagating in a nonlinear semiconductor optical amplifier ,” J. Quantum Electron. 41 , 36 – 44 ( 2005 ). [CrossRef]  

11 . D. P. Landau and K. Binder , A Guide to Monte Carlo Simulations in Statistical Physics ( Cambridge, MA/New York , 2000 ).

12 . C. Andrieu , N. De Freitas , A. Doucet , and M. I. Jordan , “ An introduction to MCMC for machine learning ,” Machine Learning 50 , 5 – 43 ( 2003 ). [CrossRef]  

13 . N. Metropolis , A. Rosenbluth , M. Rosenbluth , M. Teller , and E. Teller , “ Equation of state calculations by fast computing machines ,” J. Chem. Phys. 21 , 1087 – 1092 ( 1953 ). [CrossRef]  

14 . T. J. Farrell , M. S. Patterson , and B. Wilson , “ A diffusion theory model of spatially resolved, steady-state diffusive reflectance for the noninvasive determination of tissue optical properties in vivo ,” Med. Phys. 19 , 879 – 888 ( 1992 ). [CrossRef]   [PubMed]  

15 . G. Alexandrakis , T. J. Farrell , and M. S. Patterson , “ Accuracy of the Diffusion Approximation in Determining the Optical Properties of a Two-Layer Turbid Medium ,” Appl. Opt. 37 , 7401 – 7409 ( 1998 ). [CrossRef]  

References

  • View by:
  • |
  • |
  • |

  1. S. T. Flock , M. S. Patterson , B. C. Wilson , and D. R. Wyman , “ Monte Carlo modeling of light propagation in highly scattering tissues - 1. Model predictions and comparison with diffusion theory ,” IEEE Trans. Biomed. Eng.   36 , 1162 – 1167 ( 1989 ).
    [Crossref] [PubMed]
  2. J. M. Schmitt and K. Ben-Letaief , “ Efficient Monte Carlo simulation of confocal microscopy in biological tissue ,” J. Opt. Soc. Am. A   13 , 952 – 961 ( 1996 ).
    [Crossref]
  3. G. Yao and L. Wang , “ Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media ,” Phys. Med. Biol.   44 , 2307 – 2320 ( 1999 ).
    [Crossref] [PubMed]
  4. D. A. Boas , J. P. Culver , J. J. Stott , and A. K. Dunn , “ Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head ,” Opt. Express   10 , 159 – 170 ( 2002 ), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-3-159
    [PubMed]
  5. L. Wang , S. L. Jacques , and L. Zheng , “ MCML - Monte Carlo modeling of light transport in multi-layered tissues ,” Comput. Methods and Programs in Biomed.   47 , 131 – 146 ( 1995 ).
    [Crossref]
  6. A. Ishimaru Wave Propagation and Scattering in Random Media ( Academic Press, Inc., San Diego 1978 ).
  7. B. A. Berg and T. Neuhaus , “ Multicanonical ensemble: A new approach to simulate first-order phase transitions ”, Phys. Rev. Lett.   68 , 9 – 12 ( 1992 ).
    [Crossref] [PubMed]
  8. D. Yevick , “ Multicanonical communication system modeling - application to PMD statistics ,” IEEE Photon. Tech. Lett.   14 , 1512 – 1514 ( 2002 ).
    [Crossref]
  9. R. HolzlÖhner and C. R. Menyuk , “ Use of multicanonical Monte Carlo simulations to obtain accurate bit error rates in optical communications systems ,” Opt. Lett.   28 , 1894 – 1896 ( 2003 ).
    [Crossref] [PubMed]
  10. A. Bilenca and G. Eisenstein , “ Statistical noise properties of an optical pulse propagating in a nonlinear semiconductor optical amplifier ,” J. Quantum Electron.   41 , 36 – 44 ( 2005 ).
    [Crossref]
  11. D. P. Landau and K. Binder , A Guide to Monte Carlo Simulations in Statistical Physics ( Cambridge, MA/New York , 2000 ).
  12. C. Andrieu , N. De Freitas , A. Doucet , and M. I. Jordan , “ An introduction to MCMC for machine learning ,” Machine Learning   50 , 5 – 43 ( 2003 ).
    [Crossref]
  13. N. Metropolis , A. Rosenbluth , M. Rosenbluth , M. Teller , and E. Teller , “ Equation of state calculations by fast computing machines ,” J. Chem. Phys.   21 , 1087 – 1092 ( 1953 ).
    [Crossref]
  14. T. J. Farrell , M. S. Patterson , and B. Wilson , “ A diffusion theory model of spatially resolved, steady-state diffusive reflectance for the noninvasive determination of tissue optical properties in vivo ,” Med. Phys.   19 , 879 – 888 ( 1992 ).
    [Crossref] [PubMed]
  15. G. Alexandrakis , T. J. Farrell , and M. S. Patterson , “ Accuracy of the Diffusion Approximation in Determining the Optical Properties of a Two-Layer Turbid Medium ,” Appl. Opt.   37 , 7401 – 7409 ( 1998 ).
    [Crossref]

2005 (1)

A. Bilenca and G. Eisenstein , “ Statistical noise properties of an optical pulse propagating in a nonlinear semiconductor optical amplifier ,” J. Quantum Electron.   41 , 36 – 44 ( 2005 ).
[Crossref]

2003 (2)

C. Andrieu , N. De Freitas , A. Doucet , and M. I. Jordan , “ An introduction to MCMC for machine learning ,” Machine Learning   50 , 5 – 43 ( 2003 ).
[Crossref]

R. HolzlÖhner and C. R. Menyuk , “ Use of multicanonical Monte Carlo simulations to obtain accurate bit error rates in optical communications systems ,” Opt. Lett.   28 , 1894 – 1896 ( 2003 ).
[Crossref] [PubMed]

2002 (2)

1999 (1)

G. Yao and L. Wang , “ Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media ,” Phys. Med. Biol.   44 , 2307 – 2320 ( 1999 ).
[Crossref] [PubMed]

1998 (1)

1996 (1)

1995 (1)

L. Wang , S. L. Jacques , and L. Zheng , “ MCML - Monte Carlo modeling of light transport in multi-layered tissues ,” Comput. Methods and Programs in Biomed.   47 , 131 – 146 ( 1995 ).
[Crossref]

1992 (2)

B. A. Berg and T. Neuhaus , “ Multicanonical ensemble: A new approach to simulate first-order phase transitions ”, Phys. Rev. Lett.   68 , 9 – 12 ( 1992 ).
[Crossref] [PubMed]

T. J. Farrell , M. S. Patterson , and B. Wilson , “ A diffusion theory model of spatially resolved, steady-state diffusive reflectance for the noninvasive determination of tissue optical properties in vivo ,” Med. Phys.   19 , 879 – 888 ( 1992 ).
[Crossref] [PubMed]

1989 (1)

S. T. Flock , M. S. Patterson , B. C. Wilson , and D. R. Wyman , “ Monte Carlo modeling of light propagation in highly scattering tissues - 1. Model predictions and comparison with diffusion theory ,” IEEE Trans. Biomed. Eng.   36 , 1162 – 1167 ( 1989 ).
[Crossref] [PubMed]

1953 (1)

N. Metropolis , A. Rosenbluth , M. Rosenbluth , M. Teller , and E. Teller , “ Equation of state calculations by fast computing machines ,” J. Chem. Phys.   21 , 1087 – 1092 ( 1953 ).
[Crossref]

Alexandrakis, G.

Andrieu, C.

C. Andrieu , N. De Freitas , A. Doucet , and M. I. Jordan , “ An introduction to MCMC for machine learning ,” Machine Learning   50 , 5 – 43 ( 2003 ).
[Crossref]

Ben-Letaief, K.

Berg, B. A.

B. A. Berg and T. Neuhaus , “ Multicanonical ensemble: A new approach to simulate first-order phase transitions ”, Phys. Rev. Lett.   68 , 9 – 12 ( 1992 ).
[Crossref] [PubMed]

Bilenca, A.

A. Bilenca and G. Eisenstein , “ Statistical noise properties of an optical pulse propagating in a nonlinear semiconductor optical amplifier ,” J. Quantum Electron.   41 , 36 – 44 ( 2005 ).
[Crossref]

Binder, K.

D. P. Landau and K. Binder , A Guide to Monte Carlo Simulations in Statistical Physics ( Cambridge, MA/New York , 2000 ).

Boas, D. A.

Culver, J. P.

Doucet, A.

C. Andrieu , N. De Freitas , A. Doucet , and M. I. Jordan , “ An introduction to MCMC for machine learning ,” Machine Learning   50 , 5 – 43 ( 2003 ).
[Crossref]

Dunn, A. K.

Eisenstein, G.

A. Bilenca and G. Eisenstein , “ Statistical noise properties of an optical pulse propagating in a nonlinear semiconductor optical amplifier ,” J. Quantum Electron.   41 , 36 – 44 ( 2005 ).
[Crossref]

Farrell, T. J.

G. Alexandrakis , T. J. Farrell , and M. S. Patterson , “ Accuracy of the Diffusion Approximation in Determining the Optical Properties of a Two-Layer Turbid Medium ,” Appl. Opt.   37 , 7401 – 7409 ( 1998 ).
[Crossref]

T. J. Farrell , M. S. Patterson , and B. Wilson , “ A diffusion theory model of spatially resolved, steady-state diffusive reflectance for the noninvasive determination of tissue optical properties in vivo ,” Med. Phys.   19 , 879 – 888 ( 1992 ).
[Crossref] [PubMed]

Flock, S. T.

S. T. Flock , M. S. Patterson , B. C. Wilson , and D. R. Wyman , “ Monte Carlo modeling of light propagation in highly scattering tissues - 1. Model predictions and comparison with diffusion theory ,” IEEE Trans. Biomed. Eng.   36 , 1162 – 1167 ( 1989 ).
[Crossref] [PubMed]

Freitas, N. De

C. Andrieu , N. De Freitas , A. Doucet , and M. I. Jordan , “ An introduction to MCMC for machine learning ,” Machine Learning   50 , 5 – 43 ( 2003 ).
[Crossref]

HolzlÖhner, R.

Ishimaru, A.

A. Ishimaru Wave Propagation and Scattering in Random Media ( Academic Press, Inc., San Diego 1978 ).

Jacques, S. L.

L. Wang , S. L. Jacques , and L. Zheng , “ MCML - Monte Carlo modeling of light transport in multi-layered tissues ,” Comput. Methods and Programs in Biomed.   47 , 131 – 146 ( 1995 ).
[Crossref]

Jordan, M. I.

C. Andrieu , N. De Freitas , A. Doucet , and M. I. Jordan , “ An introduction to MCMC for machine learning ,” Machine Learning   50 , 5 – 43 ( 2003 ).
[Crossref]

Landau, D. P.

D. P. Landau and K. Binder , A Guide to Monte Carlo Simulations in Statistical Physics ( Cambridge, MA/New York , 2000 ).

Menyuk, C. R.

Metropolis, N.

N. Metropolis , A. Rosenbluth , M. Rosenbluth , M. Teller , and E. Teller , “ Equation of state calculations by fast computing machines ,” J. Chem. Phys.   21 , 1087 – 1092 ( 1953 ).
[Crossref]

Neuhaus, T.

B. A. Berg and T. Neuhaus , “ Multicanonical ensemble: A new approach to simulate first-order phase transitions ”, Phys. Rev. Lett.   68 , 9 – 12 ( 1992 ).
[Crossref] [PubMed]

Patterson, M. S.

G. Alexandrakis , T. J. Farrell , and M. S. Patterson , “ Accuracy of the Diffusion Approximation in Determining the Optical Properties of a Two-Layer Turbid Medium ,” Appl. Opt.   37 , 7401 – 7409 ( 1998 ).
[Crossref]

T. J. Farrell , M. S. Patterson , and B. Wilson , “ A diffusion theory model of spatially resolved, steady-state diffusive reflectance for the noninvasive determination of tissue optical properties in vivo ,” Med. Phys.   19 , 879 – 888 ( 1992 ).
[Crossref] [PubMed]

S. T. Flock , M. S. Patterson , B. C. Wilson , and D. R. Wyman , “ Monte Carlo modeling of light propagation in highly scattering tissues - 1. Model predictions and comparison with diffusion theory ,” IEEE Trans. Biomed. Eng.   36 , 1162 – 1167 ( 1989 ).
[Crossref] [PubMed]

Rosenbluth, A.

N. Metropolis , A. Rosenbluth , M. Rosenbluth , M. Teller , and E. Teller , “ Equation of state calculations by fast computing machines ,” J. Chem. Phys.   21 , 1087 – 1092 ( 1953 ).
[Crossref]

Rosenbluth, M.

N. Metropolis , A. Rosenbluth , M. Rosenbluth , M. Teller , and E. Teller , “ Equation of state calculations by fast computing machines ,” J. Chem. Phys.   21 , 1087 – 1092 ( 1953 ).
[Crossref]

Schmitt, J. M.

Stott, J. J.

Teller, E.

N. Metropolis , A. Rosenbluth , M. Rosenbluth , M. Teller , and E. Teller , “ Equation of state calculations by fast computing machines ,” J. Chem. Phys.   21 , 1087 – 1092 ( 1953 ).
[Crossref]

Teller, M.

N. Metropolis , A. Rosenbluth , M. Rosenbluth , M. Teller , and E. Teller , “ Equation of state calculations by fast computing machines ,” J. Chem. Phys.   21 , 1087 – 1092 ( 1953 ).
[Crossref]

Wang, L.

G. Yao and L. Wang , “ Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media ,” Phys. Med. Biol.   44 , 2307 – 2320 ( 1999 ).
[Crossref] [PubMed]

L. Wang , S. L. Jacques , and L. Zheng , “ MCML - Monte Carlo modeling of light transport in multi-layered tissues ,” Comput. Methods and Programs in Biomed.   47 , 131 – 146 ( 1995 ).
[Crossref]

Wilson, B.

T. J. Farrell , M. S. Patterson , and B. Wilson , “ A diffusion theory model of spatially resolved, steady-state diffusive reflectance for the noninvasive determination of tissue optical properties in vivo ,” Med. Phys.   19 , 879 – 888 ( 1992 ).
[Crossref] [PubMed]

Wilson, B. C.

S. T. Flock , M. S. Patterson , B. C. Wilson , and D. R. Wyman , “ Monte Carlo modeling of light propagation in highly scattering tissues - 1. Model predictions and comparison with diffusion theory ,” IEEE Trans. Biomed. Eng.   36 , 1162 – 1167 ( 1989 ).
[Crossref] [PubMed]

Wyman, D. R.

S. T. Flock , M. S. Patterson , B. C. Wilson , and D. R. Wyman , “ Monte Carlo modeling of light propagation in highly scattering tissues - 1. Model predictions and comparison with diffusion theory ,” IEEE Trans. Biomed. Eng.   36 , 1162 – 1167 ( 1989 ).
[Crossref] [PubMed]

Yao, G.

G. Yao and L. Wang , “ Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media ,” Phys. Med. Biol.   44 , 2307 – 2320 ( 1999 ).
[Crossref] [PubMed]

Yevick, D.

D. Yevick , “ Multicanonical communication system modeling - application to PMD statistics ,” IEEE Photon. Tech. Lett.   14 , 1512 – 1514 ( 2002 ).
[Crossref]

Zheng, L.

L. Wang , S. L. Jacques , and L. Zheng , “ MCML - Monte Carlo modeling of light transport in multi-layered tissues ,” Comput. Methods and Programs in Biomed.   47 , 131 – 146 ( 1995 ).
[Crossref]

Appl. Opt. (1)

Comput. Methods and Programs in Biomed. (1)

L. Wang , S. L. Jacques , and L. Zheng , “ MCML - Monte Carlo modeling of light transport in multi-layered tissues ,” Comput. Methods and Programs in Biomed.   47 , 131 – 146 ( 1995 ).
[Crossref]

IEEE Photon. Tech. Lett. (1)

D. Yevick , “ Multicanonical communication system modeling - application to PMD statistics ,” IEEE Photon. Tech. Lett.   14 , 1512 – 1514 ( 2002 ).
[Crossref]

IEEE Trans. Biomed. Eng. (1)

S. T. Flock , M. S. Patterson , B. C. Wilson , and D. R. Wyman , “ Monte Carlo modeling of light propagation in highly scattering tissues - 1. Model predictions and comparison with diffusion theory ,” IEEE Trans. Biomed. Eng.   36 , 1162 – 1167 ( 1989 ).
[Crossref] [PubMed]

J. Chem. Phys. (1)

N. Metropolis , A. Rosenbluth , M. Rosenbluth , M. Teller , and E. Teller , “ Equation of state calculations by fast computing machines ,” J. Chem. Phys.   21 , 1087 – 1092 ( 1953 ).
[Crossref]

J. Opt. Soc. Am. A (1)

J. Quantum Electron. (1)

A. Bilenca and G. Eisenstein , “ Statistical noise properties of an optical pulse propagating in a nonlinear semiconductor optical amplifier ,” J. Quantum Electron.   41 , 36 – 44 ( 2005 ).
[Crossref]

Machine Learning (1)

C. Andrieu , N. De Freitas , A. Doucet , and M. I. Jordan , “ An introduction to MCMC for machine learning ,” Machine Learning   50 , 5 – 43 ( 2003 ).
[Crossref]

Med. Phys. (1)

T. J. Farrell , M. S. Patterson , and B. Wilson , “ A diffusion theory model of spatially resolved, steady-state diffusive reflectance for the noninvasive determination of tissue optical properties in vivo ,” Med. Phys.   19 , 879 – 888 ( 1992 ).
[Crossref] [PubMed]

Opt. Express (1)

Opt. Lett. (1)

Phys. Med. Biol. (1)

G. Yao and L. Wang , “ Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media ,” Phys. Med. Biol.   44 , 2307 – 2320 ( 1999 ).
[Crossref] [PubMed]

Phys. Rev. Lett. (1)

B. A. Berg and T. Neuhaus , “ Multicanonical ensemble: A new approach to simulate first-order phase transitions ”, Phys. Rev. Lett.   68 , 9 – 12 ( 1992 ).
[Crossref] [PubMed]

Other (2)

A. Ishimaru Wave Propagation and Scattering in Random Media ( Academic Press, Inc., San Diego 1978 ).

D. P. Landau and K. Binder , A Guide to Monte Carlo Simulations in Statistical Physics ( Cambridge, MA/New York , 2000 ).

Cited By

OSA participates in Crossref's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (4)

Fig. 1.
Fig. 1.

Multicanonical Monte-Carlo (MMC) algorithm for light transport in biological media.

Fig. 2.
Fig. 2.

MMC estimators of the diffuse reflectance pdf of a semi-infinite homogeneous random medium.

Fig. 3.
Fig. 3.

MMC and CMC estimators of the diffuse reflectance pdf of a semi-infinite homogeneous random medium.

Fig. 4.
Fig. 4.

MMC and CMC estimators of the pdf of the maximum penetration depth of light in a two-layered inhomogeneous random medium.

Equations (11)

Equations on this page are rendered with MathJax. Learn more.

f Y ( y i ) dy = χ 1 L i ( x ) f X ( x ) d x
1 L i ( x ) = { 1 x L i 0 x L i
f ̂ Y CMC ( y i ) = 1 dy · N Σ n = 1 N 1 L i ( x n ) , i N
f Y ( y i ) dy = χ 1 L i ( x ) w ( x ) q ( x ) dx
f ̂ Y IS ( y i ) = 1 dy · N Σ n = 1 N 1 L i ( x n ) w ( x n ) , i N
q opt ( x ) = 1 L i ( x ) f X ( x ) χ 1 L i ( x ) f X ( x ) d x = 1 L i ( x ) f X ( x ) f Y ( y i ) dy , i N
q opt ( k ) ( x ) 1 L i ( x ) f X ( x ) f Y ( k 1 ) ( y i ) , i N
f ̂ Y ( k ) ( y i ) = { f Y ( k 1 ) ( y i ) · Σ n = 1 N 1 L i ( x n ) if at least one sample x n L i f Y ( k 1 ) ( y i ) otherwise
f Y ( y ) = 0 1 w packet · f Y , W packet no _ abs ( y , w packet ) dw packet / 0 1 w packet · f Y , W , packet no _ abs ( y , w packet ) d w packet dy
f ̂ Y , W packet no _ abs ( k ) ( y i , w packet j ) = { f Y , W packet no _ abs ( k ) ( y i w packet j ) · Σ n = 1 N 1 L ij ( x n ) if at least one sample x n L ij f Y , W packet no _ abs ( k ) ( y i w packet j ) otherwise
f ̂ Y ( k ) ( y ) = 0 1 w packet · f ̂ Y , W packet no _ abs ( k ) ( y , w packet ) dw packet / 0 1 w packet · f ̂ Y , W , packet no _ abs ( k ) ( y , w packet ) d w packet dy

Metrics