The production and delivery of a radiotherapy dose can be simulated as a series of particle interactions. Since a Monte Carlo (MC) approach is model-driven (without many of the corrections and assumptions of the calculation approaches mentioned on the treatment planning page), the precision of the dose predictions is nominally only limited by the number of particles simulated (known as histories), the accuracy of the cross-section data and the accuracy of the simulation parameters.
The simulation of particle interactions relies on accurate modelling of the geometries through which the particles are transported: for a radiotherapy treatment this means the linear accelerator and the patient anatomy. The simulation of particles in the linear accelerator provides an estimate of the fluence below the treatment head, while the simulation of particles in the patient anatomy results in a dose prediction. A number of MC-based planning systems use analytically characterised fluence maps in order to simulate dose results, though these can fail to describe effects such as multi-leaf collimator (MLC) leakage (Oh Kim et al 2001).
Accelerator Geometry
The treatment head of the linear accelerator is modelled so that a description of the radiation field incident on the patient can be found. Analytical models of the linear accelerator output have been used (Verhaegen and Seuntjens 2003) but this page will focus on complete simulations of the radiation source.
The field of radiation generated by a linear accelerator largely consists of bremsstrahlung emitted by a target material placed in the path of an electron beam. This beam is well defined, having narrow energy, angular and spatial distributions (Verhaegen and Seuntjens 2003): it is not necessary to model the electron gun or accelerating wave guide to describe it. These properties are patient-independent: they do not vary with different collimator positions or gantry angles.
The particle transport simulation begins with these electrons: once the path of one (and the paths of all particles resulting from it) has been fully simulated a new electron is generated and the process repeated, until the desired statistical variance is achieved. Thus, the first step in setting up a simulation of a radiotherapy treatment is the definition of the energy, angular and spatial distribution of this electron ‘source’. The electron beam can be considered monoenergetic (Sheikh-Bagheri and Rogers 2002), though the energy usually differs from the nominal operation energy of the linear accelerator. For example, Keall et al (2003) determined simulation energies of 6.2 MeV and 17 MeV for nominal operation energies of 6 MeV and 18 MeV respectively. The spatial distribution is commonly assumed to be a Gaussian function (Sheikh-Bagheri and Rogers 2002), where the full width at half maximum (FWHM) may be taken from the literature or determined through experimentation. Typical FWHM values include 1.0 mm (Cho et al 2005), 1.5 mm (Loewenthal et al 1992) and 0.5-3.4 mm (Jaffray et al 1993). The angular distribution of the electron source, known as the electron beam divergence, is not considered in most models. An investigation by Sheikh-Bagheri and Rogers (2002) indicated that credible divergences of up to 0.5° resulted in no observable effect on beam properties at the isocentre, when compared to simulations where the divergence was not considered. These variables are optimised, in a process known as model commissioning, until they fit experimental measurements.
The remaining components of the treatment head are defined according to the manufacturers’ specifications, with minor changes potentially made to accommodate the limitations of the MC system being used. A curved flattening filter may need to be approximated by conical volumes, and focused jaws may be treated as sliding pairs. Some components, such as the ion chamber and mirror have no significant effect on the simulated accelerator output, and may be omitted for improved performance.
Phase space descriptions contain the details of the particles resulting from a given simulation, including the particle type, energy, coordinates (x,y,z), direction cosines in the x,y,z directions (u,v,w) and the statistical weight (van Dyk 2005). These descriptions are stored in phase space files on disk, which can contain billions of entries, and can take up gigabytes of hard disk space. The contents of these files can be used as a source in further simulations. The radiation field produced by an accelerator model is recorded below the treatment head at a scoring plane for all particles incident on it.
The most computationally efficient approach, ignoring the effect of hard disk drive access times, would be to split the linear accelerator simulations into two parts: patient independent and patient dependent. This would be achieved by defining a scoring plane above the secondary collimators, where all components are identical for every treatment. The phase space above the collimator jaws, which define the field used in the treatment, could be recycled for every treatment, decreasing the required number of histories, or number of initial particles. However, as processor speeds and available memory increases, hard drive access times may limit performance.
Model Commissioning
To simulate a radiotherapy treatment accurately the accelerator model must be commissioned: the electron source parameters need to be optimised to provide the most accurate dose predictions in water. The commissioning process requires reference data, generally obtained through measurements in water phantoms as part of clinical quality assurance, for multiple field sizes (40 cm × 40 cm and 10 cm × 10 cm for example). Since the commissioning simulations are performed for very simple geometries (square fields on homogeneous water phantoms) it follows that the level of agreement here is a limit on the accuracy of simulations in more complex treatments.
The electron energy and/or FWHM of the Gaussian spatial distribution are modified iteratively for improved agreement between measurements and the dose deposition in a defined water geometry. Agreement is evaluated in three ways: depth dose profiles, off-axis dose profiles and output factors.
Depth dose profiles describe the absorbed dose in the medium along the central axis of the beam for increasing depth in water. These values are normalised using the maximum dose (which can be found within 1-2 cm depth) or the dose at a chosen depth (such as 5 or 10 cm).
Off-axis dose profiles, which can be in-plane or cross-plane (in-plane being parallel to the gantry rotation axis), allow visualisation of the penumbra and the beam broadening effect that occurs with increasing depth. The values are normalised against the central axis dose, which is not necessarily the maximum dose.
Comparison of the measured output factors (which, again, are the ratios of absorbed dose on the central axis at the maximum dose depth in a reference field to the corresponding dose value in smaller or larger fields) with the simulated values can indicate the accuracy of the model used.
The simplest method of commissioning linear accelerator models involves simulating dose distributions for a number of reasonable electron energy values (within ±1 MeV of the nominal operating energy), and then evaluating the level of agreement with reference measurements. A good metric for this is the chi-squared statistic:
where O (observed) is the simulated dose and E (expected) is the reference dose. The difference between these values is known as the residual, and the smaller these residuals are, the better the fit. To obtain a more intuitive quantification of the quality of fit a reduced chi-squared statistic can be calculated:
where N is the number of observations and σ is the statistical uncertainty associated with the observed value. Here a value of χ²red≤1 indicates excellent agreement, that is, the residuals do not exceed the uncertainties. There is no guarantee this will find the best fit: if the MC simulation results have been normalised at a point, all the results are dependent on the accuracy of the simulation at that point: the best fit for depth dose data normalised at 10 cm depth is not necessarily the best fit for the same data normalised at 5 cm depth. The technique will, regardless of this, allow the identification of a model with a good fit.
Patient Geometry
Once an accurate description of the field below the treatment head has been obtained it can be used to predict the distribution of dose in the patient. These simulations are often performed separately because the accurate mapping of dose deposition is usually not important when simulating the treatment head and this can simplify the way the geometry is defined. An accurate dose prediction requires an accurate representation of the patient and, as in treatment planning systems, is developed using CT data.
Hounsfield unit values are translated to appropriate electron- and mass- density values, in order to construct an array of material definitions with identical dimensions to the CT data. This is done using a CT conversion ramp, which is specific to the equipment and energy used to acquire the image (Walters et al 2009). The dose distribution is mapped to a dose grid which is overlayed on this data, existing as an array of dose bins: voxels in which the energy deposited is recorded. The resolution of the dose grid does not necessarily need to match the resolution of the patient images, nor do the bins need to be uniformly sized.
The dose bins which constitute the final distribution score a cumulative total of the energy that has been deposited within their boundaries during simulation. Deposition in this case occurs when a particle (generally an electron) reaches a low-energy threshold (after losing energy through interactions during transport) and is terminated. This threshold energy prevents superfluous interactions being simulated, but care needs to be taken to ensure that particles are deposited correctly, and so this energy should correspond to a path length much less than the dimension of a dose bin – Rogers et al (2009) suggest 1/3 of the smallest dimension. This value is different for electrons and photons: given equal energies, the penetration depth differs by several orders of magnitude. It has been suggested that 0.7 MeV is an appropriate threshold for electron transport and 0.01 MeV for photon transport (Rogers et al 2009).
References
- Cho SH, Vassiliev ON, Lee S, Liu HH, Ibbott GS and Mohan R, 2005. Reference photon dosimetry data and reference phase space data for the 6MV photon beam from Varian clinac 2100 series linear accelerators, Med Phys, 32(1): 137-148.
- Jaffray DA, Battista JJ, Fenster A and Munro P, 1993. X-ray sources of medical linear accelerators: Focal and extra-focal radiation, Med Phys, 20(5): 1417-1427.
- Keall PJ, Siebers JV, Libby B and Mohan R, 2003. Determining the incident electron fluence for Monte Carlo-based photon treatment planning using a standard measured data set, Med Phys, 30(4): 574-582.
- Loewenthal E, Loewinger E, Bar-Avraham E and Barnea G, 1992. Measurement of the source size of a 6- and 18-MV radiotherapy linac, Med Phys, 19(3): 687-690.
- Oh Kim J, Seibers JV, Keall PJ, Arnfield MR and Mohan R, 2001. A Monte Carlo study of radiation transport through multileaf collimators, Med Phys, 28(12): 2497-2506.
- Rogers DWO, Walters B and Kawrakow I, 2009. BEAMnrc Users Manual.
- Sheikh-Bagheri D and Rogers DWO, 2002. Sensitivity of megavoltage photon beam Monte Carlo simulations to electron beam and other parameters, Med Phys, 29(3): 379-390.
- van Dyk J, 2005. The Modern Technology of Radiation Oncology.
- Verhaegen F and Seuntjens J, 2003. Monte Carlo modelling of external radiotherapy photon beams, Phys Med Biol, 48: R107-R164.
- Walters B, Kawrakow I and Rogers DWO, 2009. DOSXYZnrc Users Manual.