OSA's Digital Library

Optics Express

Optics Express

  • Editor: Michael Duncan
  • Vol. 13, Iss. 24 — Nov. 28, 2005
  • pp: 9822–9833
« Show journal navigation

Multicanonical Monte-Carlo simulations of light propagation in biological media

A. Bilenca, A. Desjardins, B. E. Bouma, and G. J. Tearney  »View Author Affiliations


Optics Express, Vol. 13, Issue 24, pp. 9822-9833 (2005)
http://dx.doi.org/10.1364/OPEX.13.009822


View Full Text Article

Acrobat PDF (2276 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

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

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]

]-[5

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]

] are often performed through the Monte-Carlo simulation approach due to lack of analytic solutions to the radiative transfer equation [6

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

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

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]

] and importance sampling (IS) [2

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]

].

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

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]

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

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]

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

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]

] to simulations of light propagation in biological media. The MMC method was originally designed to study first-order phase transitions [7

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]

], and it has been recently applied to compute also polarization mode dispersion statistics [8

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

], bit error rates in optical communication systems [9

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]

] and noise statistics of semiconductor optical amplifiers [10

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]

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

2. Multicanonical Monte-Carlo (MMC) algorithm

The MMC method [7

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]

] 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
(1)

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
(2)

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

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

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

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

], [12

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]

] 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
(4)

where w(x)= fX(x) / q(x) is the so called importance weight [11

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

], [12

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]

] and the corresponding IS estimator of fY(yi) is expressed as

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

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

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

], [12

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]

]

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

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
(7)

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

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

Simulations of photon packet trajectories in biological media are more efficient than single photon trajectories [4

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

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]

]. 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
(9)

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

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]

]. 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
(10)

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
(11)

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

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]

].

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

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

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]

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

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

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]

], 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.

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

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]

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

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

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]

] 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]

OCIS Codes
(000.5490) General : Probability theory, stochastic processes, and statistics
(170.3660) Medical optics and biotechnology : Light propagation in tissues
(290.1350) Scattering : Backscattering

ToC Category:
Research Papers

History
Original Manuscript: October 4, 2005
Revised Manuscript: October 3, 2005
Published: November 28, 2005

Citation
A. Bilenca, A. Desjardins, B. Bouma, and G. Tearney, "Multicanonical Monte-Carlo simulations of light propagation in biological media," Opt. Express 13, 9822-9833 (2005)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-13-24-9822


Sort:  Journal  |  Reset  

References

  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), <a href="http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-3-159">http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-3-159</a> [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, 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]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Figures

Fig. 1. Fig. 2. Fig. 3.
 
Fig. 4.
 

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited