On the de-noising of Monte Carlo calculated dose distributions

This paper presents an algorithm for de-noising Monte Carlo calculated dose distributions for use in radiation treatment planning. The algorithm is a three-dimensional generalization of a Savitzky-Golay digital filter and uses an adaptive smoothing window size to reduce the probability for systematic bias. The paper also introduces five accuracy criteria that are relevant for the expected clinical use of Monte Carlo techniques, which can be used to evaluate the performance of smoothing algorithms. Using these accuracy criteria it is demonstrated that the smoothing algorithm presented here decreases the uncertainty of Monte Carlo calculated dose distributions. The corresponding decrease in necessary particle tracks ranges from a factor of 2 to a factor of 20, depending on the accuracy criterion used. It is shown that very short Monte Carlo simulations combined with smoothing deliver satisfactory dose distributions and may therefore be extremely valuable for the initial trial and error phase of the radiation treatment planning process.     

Monte Carlo modeling of radiation dose distributions in intravascular radiation therapy

Radiation dose distributions are developed for balloon and wire sources of radioactivity within coronary arteries. The Monte Carlo codes MCNP 4B and EGS4 were used to calculate dose distributions for photons and electrons at discrete energies around such sources, with and without the presence of a high-density atherosclerotic plaque. An interactive computer program was developed which then calculates dose distributions for many radionuclides by applying the emission spectra to the discrete energy grids calculated by the Monte Carlo codes, weighting appropriately for electron energy and abundance. Results for Re-186 and Re-188 balloon sources are shown in comparison to an Ir-192 wire source. The program provides dose distributions as well as estimates of activity levels needed to deliver prescribed doses to the vessel wall at selected distances from the lumen in a selected time interval. In addition, dose calculations are presented in this paper for other organs in the body, from photon radiation as well as from possible loss of liquid activity into the bloodstream in the case of a balloon rupture. These results, especially the interactive computer program permitting easy comparison of various radionuclides and their physical characteristics, will greatly facilitate the comparison process and aid in the selection of the best candidate(s) for clinical use.     

Evaluation of dose-volume metrics for microbeam radiation therapy dose distributions in head phantoms of various sizes using Monte Carlo simulations

This work evaluates four dose-volume metrics applied to microbeam radiation therapy (MRT) using simulated dosimetric data as input. We seek to improve upon the most frequently used MRT metric, the peak-to-valley dose ratio (PVDR), by analyzing MRT dose distributions from a more volumetric perspective. Monte Carlo simulations were used to calculate dose distributions in three cubic head phantoms: a 2 cm mouse head, an 8 cm cat head and a 16 cm dog head. The dose distribution was calculated for a 4 × 4 mm2 microbeam array in each phantom, as well as a 16 × 16 mm2 array in the 8 cm cat head, and a 32 × 32 mm2 array in the 16 cm dog head. Microbeam widths of 25, 50 and 75 μm and center-to-center spacings of 100, 200 and 400 μm were considered. The metrics calculated for each simulation were the conventional PVDR, the peak-to-mean valley dose ratio (PMVDR), the mean dose and the percentage volume below a threshold dose. The PVDR ranged between 3 and 230 for the 2 cm mouse phantom, and between 2 and 186 for the 16 cm dog phantom depending on geometry. The corresponding ranges for the PMVDR were much smaller, being 2-49 (mouse) and 2-46 (dog), and showed a slightly weaker dependence on phantom size and array size. The ratio of the PMVDR to the PVDR varied from 0.21 to 0.79 for the different collimation configurations, indicating a difference between the geometric dependence on outcome that would be predicted by these two metrics. For unidirectional irradiation, the mean lesion dose was 102%, 79% and 42% of the mean skin dose for the 2 cm mouse, 8 cm cat and 16 cm dog head phantoms, respectively. However, the mean lesion dose recovered to 83% of the mean skin dose in the 16 cm dog phantom in intersecting cross-firing regions. The percentage volume below a 10% dose threshold was highly dependent on geometry, with ranges for the different collimation configurations of 2-87% and 33-96% for the 2 cm mouse and 16 cm dog heads, respectively. The results of this study illustrate that different dose-volume metrics exhibit different functional dependences on MRT geometry parameters, and suggest that reliance on the PVDR as a predictor of therapeutic outcome may be insufficient.     

The denoising of Monte Carlo dose distributions using convolution superposition calculations

Monte Carlo (MC) dose calculations can be accurate but are also computationally intensive. In contrast, convolution superposition (CS) offers faster and smoother results but by making approximations. We investigated MC denoising techniques, which use available convolution superposition results and new noise filtering methods to guide and accelerate MC calculations. Two main approaches were developed to combine CS information with MC denoising. In the first approach, the denoising result is iteratively updated by adding the denoised residual difference between the result and the MC image. Multi-scale methods were used (wavelets or contourlets) for denoising the residual. The iterations are initialized by the CS data. In the second approach, we used a frequency splitting technique by quadrature filtering to combine low frequency components derived from MC simulations with high frequency components derived from CS components. The rationale is to take the scattering tails as well as dose levels in the high-dose region from the MC calculations, which presumably more accurately incorporates scatter; high-frequency details are taken from CS calculations. 3D Butterworth filters were used to design the quadrature filters. The methods were demonstrated using anonymized clinical lung and head and neck cases. The MC dose distributions were calculated by the open-source dose planning method MC code with varying noise levels. Our results indicate that the frequency-splitting technique for incorporating CS-guided MC denoising is promising in terms of computational efficiency and noise reduction.     

Characterization of a novel micro-irradiator using Monte Carlo radiation transport simulations

Small animals are highly valuable resources for radiobiology research. While rodents have been widely used for decades, zebrafish embryos have recently become a very popular research model. However, unlike rodents, zebrafish embryos lack appropriate irradiation tools and methodologies. Therefore, the main purpose of this work is to use Monte Carlo radiation transport simulations to characterize dosimetric parameters, determine dosimetric sensitivity and help with the design of a new micro-irradiator capable of delivering irradiation fields as small as 1.0 mm in diameter. The system is based on a miniature x-ray source enclosed in a brass collimator with 3 cm diameter and 3 cm length. A pinhole of 1.0 mm diameter along the central axis of the collimator is used to produce a narrow photon beam. The MCNP5, Monte Carlo code, is used to study the beam energy spectrum, percentage depth dose curves, penumbra and effective field size, dose rate and radiation levels at 50 cm from the source. The results obtained from Monte Carlo simulations show that a beam produced by the miniature x-ray and the collimator system is adequate to totally or partially irradiate zebrafish embryos, cell cultures and other small specimens used in radiobiology research.     

Denoising of electron beam Monte Carlo dose distributions using digital filtering techniques

The Monte Carlo (MC) method has long been viewed as the ultimate dose distribution computational technique. The inherent stochastic dose fluctuations (i.e. noise), however, have several important disadvantages: noise will affect estimates of all the relevant dosimetric and radiobiological indices, and noise will degrade the resulting dose contour visualizations. We suggest the use of a post-processing denoising step to reduce statistical fluctuations and also improve dose contour visualization. We report the results of applying four different two-dimensional digital smoothing filters to two-dimensional dose images. The Integrated Tiger Series MC code was used to generate 10 MeV electron beam dose distributions at various depths in two different phantoms. The observed qualitative effects of filtering include: (a) the suppression of voxel-to voxel (high-frequency) noise and (b) the resulting contour plots are visually more comprehensible. Drawbacks include, in some cases, slight blurring of penumbra near the surface and slight blurring of other very sharp real dosimetric features. Of the four digital filters considered here, one, a filter based on a local least-squares principle, appears to suppress noise with negligible degradation of real dosimetric features. We conclude that denoising of electron beam MC dose distributions is feasible and will yield improved dosimetric reliability and improved visualization of dose distributions.     

Characterization of scattered radiation in kV CBCT images using Monte Carlo simulations

Kilovoltage (kV) cone beam computed tomography (CBCT) images suffer from a substantial scatter contribution. In this study, Monte Carlo (MC) simulations are used to evaluate the scattered radiation present in projection images. These predicted scatter distributions are also used as a scatter correction technique. Images were acquired using a kV CBCT bench top system. The EGSnrc MC code was used to model the flat panel imager, the phantoms, and the x-ray source. The x-ray source model was validated using first and second half-value layers (HVL) and profile measurements. The HVLs and the profile were found to agree within 3% and 6%, respectively. MC simulated and measured projection images for a cylindrical water phantom and for an anthropomorphic head phantom agreed within 8% and 10%. A modified version of the DOSXYZnrc MC code was used to score phase space files with identified scattered and primary particles behind the phantoms. The cone angle, the source-to-detector distance, the phantom geometry, and the energy were varied to determine their effect on the scattered radiation distribution. A scatter correction technique was developed in which the MC predicted scatter distribution is subtracted from the projections prior to reconstruction. Preliminary testing of the procedure was done with an anthropomorphic head phantom and a contrast phantom. Contrast and profile measurements were obtained for the scatter corrected and noncorrected images. An improvement of 3% for contrast between solid water and a liver insert and 11% between solid water and a Teflon insert were obtained and a significant reduction in cupping and streaking artifacts was observed.     

Adaptive anisotropic diffusion filtering of Monte Carlo dose distributions

The Monte Carlo method is the most accurate method for radiotherapy dose calculations, if used correctly. However, any Monte Carlo dose calculation is burdened with statistical noise. In this paper, denoising of Monte Carlo dose distributions with a three-dimensional adaptive anisotropic diffusion method was investigated. The standard anisotropic diffusion method was extended by changing the filtering parameters adaptively according to the local statistical noise. Smoothing of dose distributions with different noise levels in an inhomogeneous phantom, a conventional and an IMRT treatment case is shown. The resultant dose distributions were analysed using several evaluating criteria. It is shown that the adaptive anisotropic diffusion method can reduce statistical noise significantly (two to five times, corresponding to the reduction of simulation time by a factor of up to 20), while preserving important gradients of the dose distribution well. The choice of free parameters of the method was found to be fairly robust.     

Correlation between CT numbers and tissue parameters needed for Monte Carlo simulations of clinical dose distributions

We describe a new method to convert CT numbers into mass density and elemental weights of tissues required as input for dose calculations with Monte Carlo codes such as EGS4. As a first step, we calculate the CT numbers for 71 human tissues. To reduce the effort for the necessary fits of the CT numbers to mass density and elemental weights, we establish four sections on the CT number scale, each confined by selected tissues. Within each section, the mass density and elemental weights of the selected tissues are interpolated. For this purpose, functional relationships between the CT number and each of the tissue parameters, valid for media which are composed of only two components in varying proportions, are derived. Compared with conventional data fits, no loss of accuracy is accepted when using the interpolation functions. Assuming plausible values for the deviations of calculated and measured CT numbers, the mass density can be determined with an accuracy better than 0.04 g cm(-3). The weights of phosphorus and calcium can be determined with maximum uncertainties of 1 or 2.3 percentage points (pp) respectively. Similar values can be achieved for hydrogen (0.8 pp) and nitrogen (3 pp). For carbon and oxygen weights, errors up to 14 pp can occur. The influence of the elemental weights on the results of Monte Carlo dose calculations is investigated and discussed.     

Clinical CT-based calculations of dose and positron emitter distributions in proton therapy using the FLUKA Monte Carlo code

Clinical investigations on post-irradiation PET/CT (positron emission tomography/computed tomography) imaging for in vivo verification of treatment delivery and, in particular, beam range in proton therapy are underway at Massachusetts General Hospital (MGH). Within this project, we have developed a Monte Carlo framework for CT-based calculation of dose and irradiation-induced positron emitter distributions. Initial proton beam information is provided by a separate Geant4 Monte Carlo simulation modelling the treatment head. Particle transport in the patient is performed in the CT voxel geometry using the FLUKA Monte Carlo code. The implementation uses a discrete number of different tissue types with composition and mean density deduced from the CT scan. Scaling factors are introduced to account for the continuous Hounsfield unit dependence of the mass density and of the relative stopping power ratio to water used by the treatment planning system (XiO (Computerized Medical Systems Inc.)). Resulting Monte Carlo dose distributions are generally found in good correspondence with calculations of the treatment planning program, except a few cases (e.g. in the presence of air/tissue interfaces). Whereas dose is computed using standard FLUKA utilities, positron emitter distributions are calculated by internally combining proton fluence with experimental and evaluated cross-sections yielding 11C, 15O, 14O, 13N, 38K and 30P. Simulated positron emitter distributions yield PET images in good agreement with measurements. In this paper, we describe in detail the specific implementation of the FLUKA calculation framework, which may be easily adapted to handle arbitrary phase spaces of proton beams delivered by other facilities or include more reaction channels based on additional cross-section data. Further, we demonstrate the effects of different acquisition time regimes (e.g., PET imaging during or after irradiation) on the intensity and spatial distribution of the irradiation-induced beta+-activity signal for the cases of head and neck and para-spinal tumour sites.     

Monte Carlo simulations of dose near a nonradioactive gold seed

The relative doses and hot/cold spot positions around a non-radioactive gold seed, irradiated by a 6 or 18 MV photon beam in water, were calculated using Monte Carlo simulation. Phase space files of 6 and 18 MV photon beams with a field size of 1 x 1 cm2 were generated by a Varian 21 EX linear accelerator using the EGSnrc and BEAMnrc code. The seed (1.2 x 1.2 x 3.2 mm3) was positioned at the isocenter in a water phantom (20 x 20 x 20 cm2) with source-to-axis distance = 100 cm. For the single beam geometry, the relative doses (normalized to the dose at 5 mm distance above the isocenter) at the upstream seed surface were calculated to be 1.64 and 1.56 for the 6 and 18 MV beams respectively when the central beam axis (CAX) is parallel to the width of the seed. These doses were slightly higher than those (1.58 and 1.52 for 6 and 18 MV beams respectively) calculated when the CAX is perpendicular to the width of the seed. Compared to the relative dose profiles with the same beam geometry without the seed in the water phantom, the presence of the seed affects the dose distribution at about 3 mm distance beyond both the upstream and downstream seed surface. For a pair of opposing beams with equal and unequal beam weight, the hot and cold spots of both opposing beams were mixed. For a 360 degree photon arc around the longitudinal axis of the seed, the relative dose profile along the width of the seed was similar to that of the opposing beam pair, except the former geometry has a larger dose gradient near the seed surface. In this study, selected results from our simulation were compared to previous measurements using film dosimetry.     

Absolute dose calculations for Monte Carlo simulations of radiotherapy beams

Monte Carlo (MC) simulations have traditionally been used for single field relative comparisons with experimental data or commercial treatment planning systems (TPS). However, clinical treatment plans commonly involve more than one field. Since the contribution of each field must be accurately quantified, multiple field MC simulations are only possible by employing absolute dosimetry. Therefore, we have developed a rigorous calibration method that allows the incorporation of monitor units (MU) in MC simulations. This absolute dosimetry formalism can be easily implemented by any BEAMnrc/DOSXYZnrc user, and applies to any configuration of open and blocked fields, including intensity-modulated radiation therapy (IMRT) plans. Our approach involves the relationship between the dose scored in the monitor ionization chamber of a radiotherapy linear accelerator (linac), the number of initial particles incident on the target, and the field size. We found that for a 10 x 10 cm2 field of a 6 MV photon beam, 1 MU corresponds, in our model, to 8.129 x 10(13) +/- 1.0% electrons incident on the target and a total dose of 20.87 cGy +/- 1.0% in the monitor chambers of the virtual linac. We present an extensive experimental verification of our MC results for open and intensity-modulated fields, including a dynamic 7-field IMRT plan simulated on the CT data sets of a cylindrical phantom and of a Rando anthropomorphic phantom, which were validated by measurements using ionization chambers and thermoluminescent dosimeters (TLD). Our simulation results are in excellent agreement with experiment, with percentage differences of less than 2%, in general, demonstrating the accuracy of our Monte Carlo absolute dose calculations.     

A voxel-based mouse for internal dose calculations using Monte Carlo simulations (MCNP)

Murine models are useful for targeted radiotherapy pre-clinical experiments. These models can help to assess the potential interest of new radiopharmaceuticals. In this study, we developed a voxel-based mouse for dosimetric estimates. A female nude mouse (30 g) was frozen and cut into slices. High-resolution digital photographs were taken directly on the frozen block after each section. Images were segmented manually. Monoenergetic photon or electron sources were simulated using the MCNP4c2 Monte Carlo code for each source organ, in order to give tables of S-factors (in Gy Bq-1 s-1) for all target organs. Results obtained from monoenergetic particles were then used to generate S-factors for several radionuclides of potential interest in targeted radiotherapy. Thirteen source and 25 target regions were considered in this study. For each source region, 16 photon and 16 electron energies were simulated. Absorbed fractions, specific absorbed fractions and S-factors were calculated for 16 radionuclides of interest for targeted radiotherapy. The results obtained generally agree well with data published previously. For electron energies ranging from 0.1 to 2.5 MeV, the self-absorbed fraction varies from 0.98 to 0.376 for the liver, and from 0.89 to 0.04 for the thyroid. Electrons cannot be considered as 'non-penetrating' radiation for energies above 0.5 MeV for mouse organs. This observation can be generalized to radionuclides: for example, the beta self-absorbed fraction for the thyroid was 0.616 for I-131; absorbed fractions for Y-90 for left kidney-to-left kidney and for left kidney-to-spleen were 0.486 and 0.058, respectively. Our voxel-based mouse allowed us to generate a dosimetric database for use in preclinical targeted radiotherapy experiments.     

Enhancement of electron beam dose distributions by longitudinal magnetic fields: Monte Carlo simulations and magnet system optimization

A Monte Carlo electron-photon transport code was developed in order to determine the effects of static, longitudinal, magnetic fields on dose distributions produced by high-energy electron beams, and to optimize the design of a superconducting magnet system. As a result of these simulations, a 20-cm-i.d., 30-cm-o.d., 15-cm-tall, single-coil, magnet system was designed that could be incorporated into a mobile treatment table for use with a standard radiation therapy accelerator. Operating at a current density of 18 kA/cm2, the magnet would produce field strengths of 1-4 T in the phantom and 0.01 T at the accelerator exit window. Magnetically enhanced dose distributions, calculated for 20- and 30-MeV electron beams, show a pronounced Bragg peak, steeper gradients to the sides and rear, and a roughly fourfold increase in the peak dose to entrance dose ratios relative to those similarly calculated without a magnetic field. These magnetically enhanced dose distributions have the potential for sparing intervening tissue when high-energy electrons are used for the treatment of deep-seated tumors.     

Monte Carlo dose verification for intensity-modulated arc therapy

Intensity-modulated arc therapy (IMAT), a technique which combines beam rotation and dynamic multileaf collimation, has been implemented in our clinic. Dosimetric errors can be created by the inability of the planning system to accurately account for the effects of tissue inhomogeneities and physical characteristics of the multileaf collimator (MLC). The objective of this study is to explore the use of Monte Carlo (MC) simulation for IMAT dose verification. The BEAM/DOSXYZ Monte Carlo system was implemented to perform dose verification for the IMAT treatment. The implementation includes the simulation of the linac head/MLC (Elekta SL20), the conversion of patient CT images and beam arrangement for 3D dose calculation, the calculation of gantry rotation and leaf motion by a series of static beams and the development of software to automate the entire MC process. The MC calculations were verified by measurements for conventional beam settings. The agreement was within 2%. The IMAT dose distributions generated by a commercial forward planning system (RenderPlan. Elekta) were compared with those calculated by the MC package. For the cases studied, discrepancies of over 10% were found between the MC and the RenderPlan dose calculations. These discrepancies were due in part to the inaccurate dose calculation of the RenderPlan system. The computation time for the IMAT MC calculation was in the range of 20-80 min on 15 Pentium-Ill computers. The MC method was also useful in verifying the beam apertures used in the IMAT treatments.     

The effect of Monte Carlo statistical uncertainties on the evaluation of dose distributions in radiation treatment planning

This paper discusses the effect of statistical uncertainties present in Monte Carlo (MC) calculated dose distributions on the evaluation of a 'cost function' that expresses the suitability of a treatment plan for the intended treatment. The mathematical derivations given are valid for any 'well-behaved' cost function. The validity of the general expressions is demonstrated using numerical examples. It is shown that random dose uncertainties lead to statistical and systematic uncertainties on the cost function. The balance between the two types of uncertainty and the desired accuracy on the cost function presents a clear criterion for the maximum acceptable MC dose uncertainties. It is demonstrated that it is possible to remove the systematic cost function uncertainty. Finally, it is shown that when the dose distribution is close to a true optimum, MC calculations of the cost function converge to the true result as one over the number of particles simulated, i.e. much faster than the individual dose uncertainties.     

Smoothing Monte Carlo calculated dose distributions by iterative reduction of noise

A smoothing algorithm based on an optimization procedure is presented and evaluated for single electron and photon beams and a full intensity modulated radiation therapy (IMRT) delivery. The algorithm iteratively reduces the statistical noise of Monte Carlo (MC) calculated dose distributions. It is called IRON (iterative reduction of noise). By varying the dose in each voxel, the algorithm minimizes the second partial derivatives of dose with respect to X, Y and Z. An additional restoration term ensures that too large dose changes are prevented. IRON requires a MC calculated one-dimensional or three-dimensional dose distribution with or without known statistical uncertainties as input. The algorithm is tested using three different treatment plan examples, a photon beam dose distribution in water, an IMRT plan of a real patient and an electron beam dose distribution in a water phantom with inhomogeneities. It is shown that smoothing can lead to an additional reduction of MC calculation time by factors of 2 to 10. This is especially useful if MC dose calculation is part of an inverse treatment planning system. In addition to this, it is shown that smoothing a noisy dose distribution may introduce some bias into the final dose values by converting the statistical uncertainty of the dose distribution into a systematic deviation of the dose value.     

Monte Carlo simulations of neutron spectral fluence, radiation weighting factor and ambient dose equivalent for a passively scattered proton therapy unit

Stray neutron exposures pose a potential risk for the development of secondary cancer in patients receiving proton therapy. However, the behavior of the ambient dose equivalent is not fully understood, including dependences on neutron spectral fluence, radiation weighting factor and proton treatment beam characteristics. The objective of this work, therefore, was to estimate neutron exposures resulting from the use of a passively scattered proton treatment unit. In particular, we studied the characteristics of the neutron spectral fluence, radiation weighting factor and ambient dose equivalent with Monte Carlo simulations. The neutron spectral fluence contained two pronounced peaks, one a low-energy peak with a mode around 1 MeV and one a high-energy peak that ranged from about 10 MeV up to the proton energy. The mean radiation weighting factors varied only slightly, from 8.8 to 10.3, with proton energy and location for a closed-aperture configuration. For unmodulated proton beams stopped in a closed aperture, the ambient dose equivalent from neutrons per therapeutic absorbed dose (H*(10)/D) calculated free-in-air ranged from about 0.3 mSv/Gy for a small scattered field of 100 MeV proton energy to 19 mSv/Gy for a large scattered field of 250 MeV proton energy, revealing strong dependences on proton energy and field size. Comparisons of in-air calculations with in-phantom calculations indicated that the in-air method yielded a conservative estimation of stray neutron radiation exposure for a prostate cancer patient.     

Investigation of scattered radiation in 3D whole-body positron emission tomography using Monte Carlo simulations

The correction of scattered radiation is one of the most challenging tasks in 3D positron emission tomography (PET) and knowledge about the amount of scatter and its distribution is a prerequisite for performing an accurate correction. One concern in 3D PET in contrast to 2D PET is the scatter contribution from activity outside the field-of-view (FOV) and multiple scatter. Using Monte Carlo simulations, we examined the scatter distribution for various phantoms. The simulations were performed for a whole-body PET system (ECAT EXACT HR+, Siemens/CTI) with an axial FOV of 15.5 cm and a ring diameter of 82.7 cm. With (without) interplane septa, up to one (two) out of three detected events are scattered (for a centred point source in a water-filled cylinder that nearly fills out the patient port), whereby the relative scatter fraction varies significantly with the axial position. Our results show that for an accurate scatter correction, activity as well as scattering media outside the FOV have to be taken into account. Furthermore it could be shown that there is a considerable amount of multiple scatter which has a different spatial distribution from single scatter. This means that multiple scatter cannot be corrected by simply rescaling the single scatter component.