Monte Carlo simulations are probabilistic simulations of complex dynamics, usually with simplifying assumptions to make the computations more efficient. The Monte Carlo approach was introduced to physics by Ulam and von Neumann in the late 1940’s to explore neutron diffusion for the Manhattan Project. From that first use, Monte Carlo techniques spread to other areas of physics, other areas of science, and domains as diverse as telecommunications network design and finance.

In a bio-optical context, you can avoid an intractable calculation of wave propagation and interaction with multitudinous cells by simulating a large ensemble of photon-bunch trajectories. Each trajectory is broken into a number of segments. For each segment, microscopic calculation of the photons’ interaction with all the molecules in a region is replaced by a single random event, the probability of which is characterized by the average absorption and scattering probabilities for the tissue. At the end of the simulation, you can calculate the statistical properties of the ensemble of photon bunches, such as mean path length or probability that the trajectory passes through a detector.

In such an approach, photon interactions are classified into two discrete types: absorption and scattering. Absorption is characterized by the absorption coefficient m_a representing Beer’s-law-type attenuation as the light propagates. Elastic scattering is quantified by the scattering coefficient m_s, which represents the probability that the photon scatters.

A Monte-Carlo simulation keeps track of photon trajectory only, and not photon phase. A photon that scatters back into its pre-scatter direction of travel is indistinguishable from a photon that didn’t scatter at all. And most scattering off biological tissue is preferentially in the forward direction. Thus, in performing a Monte-Carlo simulation the “bare” scattering coefficient µ_s, is reduced to the effective scattering coefficient µ_s’ = µ_s(1 – g), in which the “geometry factor” g quantifies the fraction of scattering that occurs in the forward direction. The effective scattering coefficient tracks only those photons that have scattered out of the initial trajectory direction, and are thus distinguishable from unscattered photons.

Consider the scattering of light in human breast tissue (see Figure 1): a situation relevant to mammography using diffuse optical tomography. To obtain a first estimate of probability of photon detection at a distant location, you might assume homogeneous breast tissue (aside from the possibility of an optically distinct tumor region). A Monte Carlo simulation follows the launch, propagation, and interactions (absorption/scattering) experienced by a photon bunch. The bunch is treated collectively: at each computational step the bunch either collectively survives, collectively exits, or – due to absorption – drops below the threshold for further tracking, and is collectively terminated. The simulation contains the following steps:

1. Launch: A photon bunch with initial weight W₀ is launched from the origin, either in a given direction (for a collimated source) or in a randomly chosen direction consistent with the numerical aperture of the source. The weight represents the probability that the bunch remains “alive,” and has not yet been terminated. The initial weight W₀=1.0.
2. Propagate: The bunch travels along its current trajectory a distance D drawn from an exponential probability distribution. The exponent is taken to be µ_tot = (µ_a + µ_s’): that is, the attenuation accounts both for photon absorption and scattering out of its current direction of propagation.Think of D as the distance the photon bunch travels before its next interaction. Using a randomly chosen propagation distance selected in this manner is a standard Monte Carlo “trick” used to ensure that the photon bunch is either collectively absorbed or collectively scattered on each step – no computational resources are used simulating “uninteresting” interaction-free propagation segments.
3. Drop: After propagating the distance D, the photon-bunch weight W is reduced to W·(1 – P_a), where P_a is the probability that the interaction is an absorption, rather than elastic scattering:
   P_a = µ_a /(µ_a + µ_s’).
4. Check exit: The trajectory is checked against the breast-boundary geometry to determine if the photons remaining in the bunch exit the breast on that step. If so, then the bunch is terminated. If the trajectory passes through the detector as it exits, and if the propagation direction falls within the acceptance angle of the detector, the bunch is flagged as “detected.” The final “detected weight” of the photon bunch at the detector is set to the current value of W, rather than to 1.The “number of detection events” may be non-integer, and is always less than the bare number of bunch trajectories that pass through the detector.
5. Roulette: At this juncture, the “fate” of the photon bunch is assessed. The bunch’s current weight is compared to a threshold value (taken in the current simulations to be 0.01). If the weight falls below this threshold value, then: The bunch trajectory is terminated with some probability P_t.
   Otherwise, the weight is multiplied (increased) by 1/P_t, and the bunch continues to the next iteration.
   This methodology is a “variance reduction technique,” standard in Monte Carlo simulations (Kahn, “Random sampling Monte Carlo techniques in neutron attenuation problems I,” Nucleonics, 6, 27, 1950 and Kahn, “Random sampling Monte Carlo techniques in neutron attenuation problems II,” Nucleonics, 6, 61, 1950).
6. Spin: Before going to the next iteration, the photon bunch’s propagation direction is randomly updated, according to the geometry factor g, using the Henyey-Greenstein scattering function probability distribution (Henyey, Greenstein, “Diffuse radiation in the galaxy,” Astrophysical Journal 93, 70, 1941):
   
   The Henyey Greenstein scattering function describes the average value of the cosine of the angular change in propagation direction upon scattering, and is used as a simple parametrization of the complete and complex microscopic scattering distribution. Its use is standard in the simulation of scattering in biological tissue[1].

[1] If ξ is a uniform probability distribution on the interval (0,1], you may obtain the above p(cosθ) by drawing from ξ according to: