Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head

: We describe a novel Monte Carlo code for photon migration through 3D media with spatially varying optical properties. The code is validated against analytic solutions of the photon diffusion equation for semi-infinite homogeneous media. The code is also cross-validated for photon migration through a slab with an absorbing heterogeneity. A demonstration of the utility of the code is provided by showing time-resolved photon migration through a human head. This code, known as ‘tMCimg’, is available on the web and can serve as a resource for solving the forward problem for complex 3D structural data obtained by MRI or CT.


Introduction
Imaging tissue with diffuse light (i.e.diffuse optical tomography) produces images with poor structural details (for instance non-invasive measurements of the human cortex have a resolution typically no better than 5-10 mm), but the images contain exquisite functional information that complements the information obtained by other imaging modalities such as MRI, X-Ray CT, and PET [1].Because of limited spatial resolution, diffuse optical tomography is increasingly being used in combination with imaging methods that provide high-resolution structural information, such as MRI, CT, X-Ray, and ultrasound [2-4].To incorporate high-resolution structural information into the diffuse optical tomography problem it is necessary to have accurate forward models for highly heterogeneous medium with arbitrary boundaries.We present a computer code for accurate forward solutions of complex 3D media and demonstrate its utility using a 3D MRI segmentation of the human head as the structure for photon migration.
During the last 15 years, the modeling of photon migration within tissue has allowed for the quantitative determination of tissue hemoglobin oxygen saturation for basic physiological research and clinical applications [1,5,6].This progress has motivated the extension of photon migration techniques to tissue imaging applications using diffuse optical tomography (DOT) methods.These imaging methods enable determination of spatially varying functional and structural information.More accurate mapping of tissue will increase the utility of photon migration technology for a broader range of applications including: optical breast imaging [7-9], functional brain imaging [6, 10-12], as well as imaging stroke [6] and head trauma [13,14].Progress in diffuse optical tomography research has been promising.Spatial variations in absorption and scattering coefficients can be experimentally imaged under a wide variety of conditions in simplified tissue phantoms with piece-wise continuous optical properties [15,16].In addition, some groups have demonstrated the possibility of imaging within more realistic head models [2, 17,18].
Improving the accuracy of diffuse optical tomography of complex tissue structures depends on the development of more sophisticated methods for solving the forward photon migration problem.While finite-element and finite-difference solutions of the diffusion equation are likely to provide accurate solutions, they need to be compared against solutions of the radiative transport equation [19,20] to verify accuracy near internal and external boundaries.
Furthermore, there are relevant conditions under which the diffusion approximation introduces significant errors including highly absorbing regions (e.g.brain bleeds) and weakly scattering regions (such as cerebral spinal fluid in the brain) [21][22][23].Hielscher et al [24] have implemented a finite-difference solution of the radiative transport equation and compared it against solutions of the diffusion equation for homogeneous and complex heterogeneous media like the human head.Arridge et al [25,26] have developed a hybrid diffusion -radiosity method for handling nonscattering voids within highly scattering media.Dorn has developed an efficient approach for utilizing the transport equation in the tomography problem [27].While these approaches hold great promise, they are complex and there are several issues involved with establishing their accuracy (e.g. the number of nodes used to discretize the boundaries, spatial and angular coordinates).Monte Carlo solution methods, on the other hand, are conceptually simpler to implement and rely on fewer assumptions, but at the expense of computational speed.
We describe a novel Monte Carlo code that allows us to obtain results, for example, in a complex 3D head model with a signal-to-noise ratio greater than 100 up to distances of 30 mm with a 1 mm 2 detector with 10 8 photons propagated within 5-10 hours of computer time on a Pentium III 1000 MHz CPU.The method is validated against an accepted analytic solution for a semi-infinite medium, and cross-validated for a slab geometry with an absorbing inclusion.We then illustrate the utility of the Monte Carlo method with a novel simulation and movie of time-resolved photon migration through a human head and deduce the depth sensitivity of different measurement types.The human head model was obtained from a 3D segmented anatomical MRI.The ability to perform such simulations in a medium with arbitrary boundaries and spatial variation in the optical properties, in our publicly available code 'tMCimg' [28], is a significant advancement over the capabilities of the widely used and publicly available code known as 'MCML' developed by Wang and Jacques [29] which models photon propagation through layered media with planar boundaries.

Method
The rules used in Monte Carlo programs for photon migration in biological tissue are straightforward as described by Wang and Jacques [29,30].Differences between programs arise in how the photon fluence or flux is recorded during the Monte Carlo simulation.Below, we briefly describe the rules for photon migration and how we record the photon fluence within the medium and the flux of photons exiting the medium.
To begin, the initial position and direction of the photon is defined.In our case, we consider a point source of photons with a defined direction, entering the medium on the surface.The source can be diverging, as defined by the source numerical aperture (NA), in which case the azimuthal angle is determined by a random number uniformly distributed from 0 to 2π and the elevation (angle from the source axis) is determined from another random number distributed uniformly from 0 to sin -1 (NA).
Given the initial position and direction of the photon, the length to the first scattering event is calculated from an exponential distribution.Absorption of the photon is considered by decreasing the photon weight by exp(-µ a L) where µ a is the absorption coefficient and L is the length traveled by the photon [29].The photon is moved this length.A scattering angle is calculated using the probability distribution given by the Henyey-Greenstein phase function [29].A new scattering length is then determined from an exponential distribution.The photon is propagated this new length in the new direction.This process continues until the photon exits the medium or has traveled longer than a predefined period of time.We typically stop the photon after 10 ns since the probability of photon detection in tissue after such a period of time is exceedingly small.When the photon does try to leave the medium, the probability of an internal reflection is calculated using Fresnel's equation [29,31].If a reflection occurs, then the photon is reflected back into the medium the appropriate distance and migration continues.Otherwise, the migration of that particular photon halts and a new photon is launched into the medium at the predefined source location.
We consider photon migration through media with spatially varying optical properties by representing the optical properties within voxels on a cubic grid.We typically consider a grid spacing of 1 mm, but this can be scaled.As the photon is propagated from one scattering event to the next, a check is made every 1 grid spacing to see if the scattering or absorption coefficient has changed.If the scattering coefficient has changed, the remainder of the scattering length is renormalized by µ s old /µ s new where µ s is old and new scattering coefficient, as described in Jacques and Wang [32].If the absorption coefficient has changed, then the photon weight is decreased with the new absorption coefficient [32].The scattering angle is determined by the value of the scattering anisotropy factor, g, within the particular voxel containing the scattering event.
The photon fluence within the medium is determined by accumulating the photon weight every 1 grid spacing in the voxel corresponding to the present position of the photon.When a photon exits the scattering medium and enters the surrounding non-scattering medium, the exiting photon flux is determined by accumulating the photon weight in a bin representing the non-scattering voxel that the photon entered.In addition, if the photon exits in a location that has been predefined as a detector location, then information is recorded in a history file identifying the detector and pathlength of the photon in each tissue type prior to being detected.The tissue types are specified as input to the Monte Carlo code and refer to voxel groups with identical optical properties.The code supports up to 100 different tissue types.This information can be used in post-processing to analyze the effect of absorption changes within different tissue types as described further below.
After all simulated photons have been propagated, it is necessary to normalize the exiting photon flux and the photon fluence within the medium.The exiting photon flux , J out (r), is divided by the number of simulated photons.Normalizing the photon fluence, Φ(r), is more involved.To conserve energy, the exiting photon flux plus the number of photons absorbed in the medium must equal the number of simulated photons, which we normalize to 1.The number of photons absorbed at a given point is Φ(r) µ a (r).Therefore, where r i indicates the position of each surface and volume voxel, V voxel is the voxel volume, and A i is the area of the surface element at position r i .The normalization factor for Φ(r) is determined from this relation.We now described how the detected photons stored in the history file are analyzed.First, it is necessary to realize that the photon fluence can be written in the form [33] ( ) ( ) ( ) where Φ(t) is the measured photon fluence, N photons (t) is the number of photons collected in a time-gate of width ∆t centered at time t, exp(-µ a,j L i,j ) accounts for the effects of absorption in each region where L i,j is the pathlength of photon i through region j, and the photon migration time is related to the photon pathlength by the speed of light in the medium.The history file of detected photons provides all the necessary information for calculation of eq.(2).

Solutions of the Diffusion Equation for Comparison with Monte Carlo
The solution of the photon diffusion equation [19,31,34] for a semi-infinite homogeneous medium for a continuous-wave point source of light is given by [35,36] ( ) ( ) ( ) (3) and for a time-domain measurement with an impulse point source of light [34,36] ) is the photon fluence at the detector position r d at time t, generated by a point source of amplitude S at position r s .D=v/(3µ s ') is the photon diffusion coefficient [37,38], µ s ' is the reduced scattering coefficient, µ a is the absorption coefficient, and v is the speed of light in the medium.The semi-infinite boundary condition is satisfied by the method of images [31,36].The position of the image source is indicated by r s,i .
When the optical properties are spatially varying, there are two standard approaches to finding approximate linear solutions: the first Born approximation [19,39] Φ=Φ o +Φ pert (5) and the Rytov approximation [39] Φ =Φ o exp(Φ pert ).( 6) Each approach decomposes the total fluence Φ into Φ o which only depends on the background optical properties µ ao and µ so ', and Φ pert which is linearly related to spatial variations in the optical properties δµ a and δµ s '.Historically, the Born approximation is more common, but for the diffuse optical tomography problem the Rytov approximation tends to be more accurate as it accounts for some non-linear saturation due to increasing perturbation in the absorption coefficient [10].For the Rytov approximation, assuming absorption variations only, r s and r d are the position of the source and detector respectively, G is the Greens function of the photon diffusion equation for the background optical properties given the boundary conditions.For the Born approximation, Φ pert is not normalized by Φ o (r s ,r d ) If the background is homogeneous then Φ o and G can be expressed analytically in some simple geometries [34,40], where the difference between Φ o and G is amplitude factor S (see eq. ( 4)).

Validation of the Monte Carlo Code in a Homogeneous Semi-Infinite Medium
The first step in validating the Monte Carlo code is to compare the results for the remitted flux of photons from a semi-infinite homogeneous medium with an analytic solution of the diffusion equation as has been done in the past by several others [29,36].For the simulation, a collimated point source was incident on a semi-infinite medium with a scattering coefficient of 1 mm -1 , a scattering anisotropy of g = 0.01, and an absorption coefficient of 0.005 mm -1 .The medium had dimensions of 60 x 60 x 60 mm with the source positioned at (x,y,z) = (30,30,1) mm.All boundaries were treated as index matched.The source was sufficiently far from the edges so that the medium is effectively semi-infinite.A simulation was executed with 10 8 photons which took approximately 6 hours on a 1 GHz Pentium 3.

(B) (A)
The results for the continuous-wave simulation are shown in fig. 1. Fig. 1a compares the radially resolved flux of photons exiting the medium determined by Monte Carlo (eq.( 2)) and eq.(3), while fig.1b compares the fluence field within the medium with contours drawn every 10 dB.In both cases we observe decent agreement.In fig.1a, the systematically smaller Monte Carlo result at larger separations results from the influence of the edge at 30 mm which is not modeled in the solution of the diffusion equation.The discrepancy observed in the fluence field in fig.1b results from the isotropic point source approximation, within the diffusion equation, of the collimated source.This isotropic point source approximation gives rise to decent agreement for the remitted flux of photons on the surface (as shown in fig.1a), but at the expense of decent agreement within the highly scattering medium, especially in the forward direction.Better agreement is found within the medium by extending the isotropic point source slightly deeper into the medium, but at the expense of reduced agreement for the remitted flux.A more accurate solution is to model the exponentially attenuated collimated source within the diffuse equation [31].
The results for the time-domain simulation are shown in fig. 2. In fig.2a the timeresponse profile is shown for the remitted flux at x = 15 mm from the source and for the photon fluence at position (x, z) = (15, 10) mm within the medium, for both the Monte Carlo simulation and the diffusion equation.Fig. 2b compares the spatial distribution of the photon fluence within the medium by plotting the half-maximum contours at time points of 0.1 ns, 0.5, 1.0, 1.5, and 2.0 ns.The temporal and spatial agreement is fairly good.The systematic decrease for the Monte Carlo photon fluence within the medium at longer times results from the edge effects.At 0.1 ns the Monte Carlo spatial distribution is much more compact than that for the diffusion equation as expected from the non-casual behavior of the diffusion equation.We have three approaches for analyzing the effect of an absorbing inclusion with different absorption coefficients on the detection of diffuse photons: 1) run a Monte Carlo simulation for each absorption coefficient, 2) utilize the history file of detected photons which records the pathlength each detected photon spent in each piece-wise constant region, and 3) utilize the Born and/or Rytov approximation which can be produced from the product of the fluence field from a source into the medium, Φ o (r s ,r), and from the detector into the medium, G(r,r d ),

Validating the Monte Carlo code in a slab geometry with an absorbing inclusion
(see eq. ( 7), so called adjoint fields [41,42]).Comparing these three approaches serves to cross-validate the methodology for media with spatially varying optical properties.The comparison is performed in a 40 mm thick slab with a cubic absorbing inclusion 15 x 15 x 15 mm centered 17 mm into the slab from the source.The source and detector are positioned directly across each other centered on the inclusion (see fig. 3a).The scattering coefficient of the medium was set to 1 mm -1 , the scattering anisotropy factor was set to 0.01, the absorption coefficient of the background was set to 0.005 mm -1 , while that of the inclusion was increased from 0.005 mm -1 to 0.065 mm -1 .A full Monte Carlo simulation (which traces out the scattering path of every photon) was executed with the different absorption coefficients to calculate the flux of photons received at the detector.These results normalized by the flux detected with no absorbing inclusion are shown by the square symbols in fig.3b.Notice that as the absorption coefficient initially increases, there is an approximate linear decrease in the intensity which then saturates for higher absorption coefficients.The detected flux of photons for different absorption coefficients can also be calculated from the photon history file produced by a single Monte Carlo simulation, as described above in eq. ( 2).These results are shown by the solid line in fig.3b.As expected, the results produced by the history file (generated by a previous full Monte Carlo simulation but using a different value of the absorption coefficient in the post-processing) agree with full Monte Carlo results.Finally, the results produced by the Born and Rytov approximations are shown by the dashed and dot-dashed lines respectively.These approximations are linear and exponential respectively and therefore show the expected agreement for small changes in the absorption coefficient, but the approximations break down for larger absorption coefficients as evidenced by the poor agreement with the non-linear Monte Carlo results.Notice that the Rytov approximation shows the expected better agreement with the Monte Carlo results.Finally, in fig. 4 we show the photon fluence within the medium without the absorber and the change with the absorbing inclusion demonstrating the shadowing effect of the absorbing inclusion.

Full 3D Head
Figure 5 shows an anatomical MRI of a human head, segmented into five tissue types (air, scalp, skull, cerebral spinal fluid, and gray/white matter; see [43]), with a contour overlay indicating the photon migration spatial sensitivity profile for (a) continuous-wave, (b) 200 MHz modulation, and (c,d) pulsed measurements.The head was contained within 151x171x232 1mm 3 voxels and the Monte Carlo simulation recorded the temporal response and took approximately 10 hours.The continuous-wave and 200 MHz results were obtained from the modulus of Fourier transform of the temporal response.One contour line is shown for each half order of magnitude (10 dB) signal loss, and the contours end after 3 orders of magnitude in loss (60 dB).For the 3D Monte Carlo simulation, we assumed that µ s ' = 1 mm -1 and µ a =0.04 mm -1 for the scalp and skull, µ s ' = 0.01 mm -1 and µ a =0.001 mm -1 for the CSF, and µ s ' = 1.25 mm -1 and µ a =0.025 mm -1 for the gray/white matter.Note how the contours extend several millimeters into the brain tissue, indicating sensitivity to changes in cortical optical properties.The depth penetration difference between the continuous-wave and 200 MHz measurements is difficult to discern.A ratio of the two sensitivity profiles (not shown) shows that the 200 MHz profile is shifted slightly towards the surface.The time-domain Figure 6.A movie of the propagation of a pulse of light through a 3D human head.The color scale is logarithmic and spans 10 orders of magnitude from a peak in the dark red to a minimum in the dark blue.The AVI movie file size is 1.0 Mega-Bytes.sensitivity profiles suggest the possibility of obtaining greater penetration depths in the head from measurements made at longer delay times.This is further supported by the movie of the temporal evolution of a light pulse within the head as shown in fig.6.This movie illustrates the usefulness of the code for visualizing the temporal evolution of photon migration through a heterogeneous medium.
In fig.7 we show an estimate of the signal-to-noise ratio obtained by running 10 8 photons on the 3D human head, as determined from running 9 independent Monte Carlo simulations with different random number seeds.Fig. 7a shows the flux of photons exiting from the head as a function of distance from the source voxel to each voxel into which photons escaped, i.e. voxels describing the air surrounding the head.This flux is normalized by the number of simulated photons and of course depends on the 3D geometry and spatially varying optical properties.The expected exponential decay of the photon flux with distance is observed.The structure is seen in the data because of the cross-sectional area of different air voxels against the head.Some air voxels adjoined the head on 1, 2, or more surfaces, while some voxels only touched at the corner.The needed correction factor is just the effective cross-sectional area of each air voxel (see eq. ( 1)).At present, the code does not correct for this crosssectional area for a curved surface.The noise was determined from the standard deviation of the 9 independent Monte Carlo simulations.This deviation is seen in fig.7a  Note that an SNR greater than 100 is achieved for separations smaller than 2 cm, while an SNR greater than 10 is still found for separations up to 5 cm.The SNR can be improved by the square-root of the area of the detector by averaging over neighboring voxels.The voxel size in these measurements were 1 mm.Typical detector sizes for human head measurements are 3-5 mm.Using detectors of this size will increase the SNR by a factor of 3-5.The SNR will also improve as the square-root of the number of simulated photons.Finally, the SNR will also improve by decreasing the reduced scattering coefficient and/or the absorption coefficient so that a greater number of photons are able to travel greater distances through the tissue (note that reducing the scattering coefficient can actually decreasing the SNR within a few millimeters of the source).

Summary
We have described a computer code for calculating the migration of light through 3D highly scattering media with spatially varying optical properties and arbitrary boundary conditions.As diffuse optical tomography is combined with and guided by other imaging modalities with superior spatial resolution (and thus superior structural information), it becomes increasingly important to have an accurate photon migration forward solver.This Monte Carlo computer code can take structural information provided by MRI, or X-Ray CT for example, and accurately solve the photon migration forward problem in a reasonable amount of time.In addition to allowing validation of more efficient forward solvers, this code should prove useful on its own.As a few examples, this code could be used for investigating the feasibility of quantitative tissue characterization with different data types, cerebral oximetry with time domain data for example, aiding in the development of an optimal source-detector geometry for imaging brain function, as well as serving as the forward solver for MRI guided diffuse optical tomography of the head.DAB acknowledges funding from Advanced Research Technologies, NIH R29-NS38842, NIH P41-RR14075 and from the Center for Innovative Minimally Invasive Therapies.This research was funded in part by the US Army, under Cooperative Agreement No. DAMD17-99-2-9001.The material presented does not necessarily reflect the position or the policy of the Government, and no official endorsement should be inferred.A.K.D. acknowledges funding from NIH K25 NS41291-01.
Figure 1.(A) A comparison of the flux of photons remitted from a semi-infinite medium as calculated by diffusion theory and Monte Carlo.(B) A comparison of the photon fluence within a semi-infinite medium.The index-matched surface is at a depth of 0 mm.

Figure 2 .
Figure 2. (A) A comparison of diffusion theory and Monte Carlo for the temporal response to a pulse of light as measured on the surface and within the medium.The black points is for the remitted flux of light 15 mm from the source.The red points indicate the fluence within the medium at a depth of 10 mm and displaced laterally 15 mm.(B) The comparison of the isocontours for diffusion theory and Monte Carlo within the medium at 0.1, 0.5, 1.0, 1.5, and 2.0 ns after the pulse of light.

Figure 3 .Figure 4 .
Figure 3. (A) The geometry for the Monte Carlo simulation with an absorbing inclusion.(B) The relative decrease in the detected photon flux is shown as the absorption coefficient of the inclusion is increased.A comparison of 4 methods is made.See text for details.

Figure 7 .
Figure 7.An estimate of the SNR by running 10 8 photons on the 3D human head model.