A Monte Carlo model of light propagation in tissue

1. A Monte Carlo model of light propagation in tissue An introduction to Monte Carlo techniques ENGS168 Ashley Laughney November 13th, 2009

2. Overview of Lecture • Introduction to the Monte Carlo Technique • Stochastic modeling • Applications (with a focus on Radiation Transport) • Random sampling and Probability Distribution Functions • Monte Carlo Treatment of the Radiation Transport Problem • Photon initialization • Generating the propagation distance • Internal reflection • Photon absorption • Photon termination • Photon scattering • Calculation of observable quantities • Implementing Monte Carlo - Sample Code • Radiation treatment planning

3. Stochastic Modeling • Launch a photon (or particle) • Sample physical properties using random variables (random sampling) • Let the photon evolve through the system • Keep track of all important parameters until the photon dies or exits the problem space • Summarize and infer results after ENOUGH photons interact. Ref: VenkatKrishnaswamy

4. Applications (a few) • Modeling photon transport in tissue • Calculating dose distributions for radiation therapy • Solving the problem of neutron diffusion in a fissionable material • Calculating financial derivatives and evaluating investment value and market behavior • Wireless network design • Numerical integration Ref: Wikipedia

5. The radiation transport problem Radiative Transfer Equation (RTE): divergence extinction scattering source Diffusion Approximation: with diffusion coefficient, Ref: Ambrocio, Master Thesis, http://appliedmath.ucmerced.edu/theses/ambrocio_2009.pdf

6. Tissue optical properties

7. Probability density functions (PDF) • A probability distribution describes the range of possible values that a random variable can attain and the probability that the value of the random variable is within any subset of that range. • i.e., what is the chance of getting a value for every possible outcome of the random variable • Coin toss • A PDF is the functional form of a probability distribution Ref: VenkatKrishnaswamy

8. Beer’s Law - a probability distribution • The fraction of photons that will survive after a distance d’<d can be seen as a probability distribution over all d’. • How far will a photon travel in an absorbing medium without an interaction? • 100% chance it will travel 0cm • ~40% chance it will travel 1cm Ref: VenkatKrishnaswamy

9. Cumulative distribution function (CDF) • The probability that a measurement yielding a value of x will lie in the interval [0,x1] is given by the cumulative distribution function. • Where, , represents the probability density distribution of a set size, x є [0, Inf] , that a photon takes between any two scattering events. Beer’s law as a cumulative distribution

10. * CDFs and Random Sampling* • The CDF is always uniformly distributed on the interval [0,1]. • Sampling by inversion of the CDF • Sample a random number ξ from U[0,1] • Equate ξ with the CDF, F(x) = ξ 3) Invert the CDF and solve for x Cumulative distribution functions associated with physical processes can be sampled using random numbers via direct inversion. Ref: http://www.phy.ornl.gov/csep/CSEP/GIFFIGS/MCF12.GIF

11. Pseudo random numbers • Computers can generate random numbers, , with a uniform PDF • The associated CDF is given by Ref: Jacques, Prahl, http://omlc.ogi.edu/software/mc/

12. Flowchart for variable stepsize MC • Implicit Capture ~ each photon launched into the medium is thought to represent a photon packet, where each packet enters the medium carrying the photon weight (i.e. 1J) Ref: Prahl,  A Monte Carlo Model of Light Propagation in Tissue, SPIE (1989)

13. 1. Photon initialization • N photons are launched, each with a "photon weight" initially set to 1 (computationally efficient) • Start coordinates for each photon are identical • Photon’s initial direction chosen via convolution with the beam shape Example of Convolution Ref: Prahl,  A Monte Carlo Model of Light Propagation in Tissue, SPIE (1989) Image: http://support.svi.nl/wikiimg/ft_1.png

14. 2. Propagation distance • A fixed stepsize, Δs, must be small relative to the average mean free path of a photon in tissue. • It is more efficient to choose a different stepsize for each photon step; the PDF for the Δs follows Beer’s law. • A function of a random variable, ξ : [0,1], that is distributed uniformly and yields a random variable with this distribution: Ref: Prahl,  A Monte Carlo Model of Light Propagation in Tissue, SPIE (1989)

15. 3. Internal reflection • The probability that the photon is internally reflected is determined by the Fresnel reflection coefficient • A random number uniformly distributed between 0 and 1 is used to determine if the photon is reflected or transmitted.  Internal reflection • i.e. For a semi-infinite slab, the internally reflected position is updated by only changing the z-component of photon coordinates // angle of incidence on the boundary // angle of transmission given by Snell’s law Ref: Prahl,  A Monte Carlo Model of Light Propagation in Tissue, SPIE (1989)

16. 4. Photon absorption • After each propagation step, a fraction of the photon packet is absorbed and the remainder is scattered. • Or generate photon absorption (weight) according to randomly generated step size and Beer’s law a = the single particle albedo (fraction scattered) // New weight assigned to surviving photon packet Ref: Prahl,  A Monte Carlo Model of Light Propagation in Tissue, SPIE (1989)

17. 5. Photon termination • Propagating a photon packet with minimal weight yields little information. How is the packet terminated? • Roulette is used to terminate a photon when its weight falls below some minimum • The roulette gives a photon of weight w, one chance in m, of surviving with weight mw. Otherwise, w  0 • Unbiased elimination, conservation of energy Ref: Prahl,  A Monte Carlo Model of Light Propagation in Tissue, SPIE (1989)

18. 6. Photon scattering • The PDF for the scattered cosine of the deflection angle (cosθ) in tissue is characterized by the Henyey-Greenstein phase function. • The azimuth angle is uniformly distributed between [0,2π], and may be generated by multiplying a random number ξ:[0,1] by 2π • The photon is scattered at an angle (θ, ) * For isotropic scattering, *Assumes phase function has no azimuth dependence Ref: Prahl,  A Monte Carlo Model of Light Propagation in Tissue, SPIE (1989)

19. 7. Observable quantities • Analytic Solution to RTE • fluence rate resulting from photons launched at a single point corresponds to the Green’s function for the medium. • Monte Carlo Solution to RTE  defines grid over solution space and scores physical quantities (reflection, transmission, absorption = energy deposited) at each grid element as the program traces N photons. Ref: VenkatKrishnaswamy

20. Abbreviated review

21. Implementing Monte Carlo: A steady state example (“mc321.c”) • The following slides walk through a steady-state Monte Carlo simulation by Steve Jacques and Scott Prahl at the Oregon Medical Laser Center. All code discussed is available online: http://omlc.ogi.edu/software/mc/ • Problem Definition: • Photons are launched from an isotropic point source of unit power, P = 1W into an infinite, homogeneous medium with no boundaries • The medium has the optical properties of absorption, scattering and anisotropy • N Photons are launched, each initialized with weight, W = 1. • Solution space is divided into an array of bins (position is defined by distance r from source), 3 options – spherical shell, cylindrical shell, planar shell • Each bin accumulates photon weights deposited due to absorption by N photons Map of fluence Array containing accumulated weight of absorbed photons Concentration of Photons

22. Implementing Monte Carlo: Definitions of variables and arrays

23. Implementing Monte Carlo: User input

24. Implementing Monte Carlo: Launching photons

25. Implementing Monte Carlo: Moving photons ~ HOP Beer’s Law 

26. Implementing Monte Carlo: Moving photons ~ DROP

27. Implementing Monte Carlo: Moving photons ~ SPIN

28. Implementing Monte Carlo: Moving photons ~ SPIN

29. Implementing Monte Carlo: Moving photons ~ CHECK ROULETTE

30. Implementing Monte Carlo: Output bin arrays as fluence rate

31. Implementing Monte Carlo: Example Output Ref: Jacques, Prahl, http://omlc.ogi.edu/software/mc/

32. Advance Application: Radiation Treatment Planning • Motivation: • Mortality caused by (1) providing too little radiation to the tumor for cure, or (2) providing too much radiation to nearby healthy tissue. • Need: Improved radiation therapy planning • Ionizing and non-ionizing radiation can be used • Photon therapy accounts for 90% of all radiation treatment in US • Photons, electrons, neutrons, heavy charged particles (protons) • Dose distribution is key parameter of interest in treatment planning Diagnosed with life-threatening forms of cancer annually Receive radiation treatment Die anyway Are considered curable Ref: Venkat Krishnaswamy, https://www.llnl.gov/str/Moses.html

33. Current dose estimation techniques • Generate 3D electron-density map of body using stack of CTs • Model the body as a homogeneous bucket of water • One way to estimate dose in tissue is to use a water phantom • use ionization chambers/chamber arrays to detect dose distribution • currently used in clinical treatment planning • complicated experiments • heterogeneities hard to model Ref: VenkatKrishnaswamy

34. Monte Carlo-based treatment planning • Voxelize medium of interest using CT/MR patient images • Compute solution space geometry • Assign material data to each voxel (from atomic and nuclear-interaction databases) • Launch radiation particles one at a time and let them evolve • Store accumulated dose per voxel for N radiation particles • After following many particle histories, an accurate estimation of dose is obtained. Ref:VenkatKrishnaswamy

35. PEREGRINE3D Monte Carlo Treatment Planing Ref:https://www.llnl.gov/str/Moses.html

36. PEREGRINE Defining the radiation source and patient • 3D Transport mesh of patient generated from stack of CT images • Radiation Source • upper portion of accelerator does not vary between treatments, but the lower portion is modified by collimators, blocks and wedges to customize patient treatment. • PEREGRINE library accounts for modification in lower half of accelerator Ref: https://www.llnl.gov/str/Moses.html

37. PEREGRINE Calculating Dose Five-field treatment for a lung tumor; 6 MV photon beam Seven-field conformal boost to the prostate; 18MV photon beam Predicted dose build up for treatment of a brain tumor Ref: https://www.llnl.gov/str/Moses.html, Venkat

38. Suggested Reading Literature: ED Cashwell, CJ Everett, "A Practical Manual on the Monte Carlo Method for Random Walk Problems,“ Pergammon Press, New York, 1959. BC Wilson, G Adams, A Monte Carlo model for the absorption and flux distributions of light in tissue, Med. Phys. 10:824-830, 1983. L Wang, SL Jacques, L Zheng, MCML - Monte Carlo modeling of light transport in multi-layered tissues, Computer Methods and Programs in Biomedicine 47:131-146, 1995. L Wang, SL Jacques, "Monte Carlo Modeling of Light Transport in Multi-layered Tissues in Standard C,” 1992-1998. Download as 177-page manual in pdf format. Source Code: http://mcnp-green.lanl.gov/index.html Radiation MC, state of the art http://omlc.ogi.edu/software/mc/ Photon MC, MCML and others http://omlc.ogi.edu/software/polarization Photon MCML (Vector)