Benchmarking and validation of a Geant4–SHADOW Monte Carlo simulation for dose calculations in microbeam radiation therapy

Microbeam radiation therapy (MRT) is a synchrotron-based radiotherapy modality that uses high-intensity beams of spatially fractionated radiation to treat tumours. The rapid evolution of MRT towards clinical trials demands accurate treatment planning systems (TPS), as well as independent tools for the verification of TPS calculated dose distributions in order to ensure patient safety and treatment efficacy. Monte Carlo computer simulation represents the most accurate method of dose calculation in patient geometries and is best suited for the purpose of TPS verification. A Monte Carlo model of the ID17 biomedical beamline at the European Synchrotron Radiation Facility has been developed, including recent modifications, using the Geant4 Monte Carlo toolkit interfaced with the SHADOW X-ray optics and ray-tracing libraries. The code was benchmarked by simulating dose profiles in water-equivalent phantoms subject to irradiation by broad-beam (without spatial fractionation) and microbeam (with spatial fractionation) fields, and comparing against those calculated with a previous model of the beamline developed using the PENELOPE code. Validation against additional experimental dose profiles in water-equivalent phantoms subject to broad-beam irradiation was also performed. Good agreement between codes was observed, with the exception of out-of-field doses and toward the field edge for larger field sizes. Microbeam results showed good agreement between both codes and experimental results within uncertainties. Results of the experimental validation showed agreement for different beamline configurations. The asymmetry in the out-of-field dose profiles due to polarization effects was also investigated, yielding important information for the treatment planning process in MRT. This work represents an important step in the development of a Monte Carlo-based independent verification tool for treatment planning in MRT.

1. Introduction

Microbeam radiation therapy (MRT) is an exciting development in the treatment of cancer. It is based on the apparent ability of healthy tissue, and inability of tumour tissue, to withstand high doses of ionizing radiation when delivered in a spatially fractionated manner (Slatkin et al., 1992; Bouchet et al., 2010; Crosbie et al., 2010). MRT is currently limited to synchrotron radiation sources and spatial fractionation is achieved using a high-precision multi-slit collimator (MSC) with a nominal slit width of 50 µm and a distance between slits of 400 µm. The low divergence of the synchrotron radiation ensures that spatial fractionation is maintained with depth in the patient; moreover, due to the high fluence provided by insertion devices, the irradiation takes place in a sufficiently short time so as to avoid blurring of the spatial fractionation caused by movement of the vasculature. The ID17 biomedical beamline at the European Synchrotron Radiation Facility (ESRF) has pioneered synchrotron-based MRT, with over a decade of pre-clinical research demonstrating the feasibility of the technique (Bräuer-Krisch et al., 2010; Bouchet et al., 2010, 2013; Schültke et al., 2008). The treatment of spontaneous tumours in veterinary oncology patients (cats and dogs) is imminent.

As for any radiotherapy modality, the ability to accurately predict the distribution of dose in the patient is critical. The most accurate method of doing so is the Monte Carlo (MC) technique; however, due to prohibitively long computation times, analytical techniques and simpified MC methods are often used by treatment planning systems (TPS) to iteratively optimize the plan. Nonetheless, MC codes are commonly used for independent patient-specific pre-treatment verification of TPS dose calculations (Bush et al., 2008). Such independent calculations are required for mitigating the risk of calculation errors, as highlighted in IAEA Technical Reports Series No. 430 (Sharpe, 2006). In addition to treatment plan verification, MC codes are useful for optimizing the design of radiation dosimeters used for quality assurance in radiotherapy (Lo Meo et al., 2011; Lian et al., 2011; Rosenfeld et al., 2005).

MC calculations in MRT require careful consideration of a number of factors which do not apply to conventional megavoltage radiotherapy, namely the photon beam in MRT is highly polarized in the electron orbital plane (Hugtenburg et al., 2010), the photon beam energy is significantly lower (in the range 30–500 keV) (Siegbahn et al., 2006), and there is a need for dose calculations on the micrometre scale. The latter is required for calculation of the quantity related to the degree of spatial fraction, the peak-to-valley dose ratio (PVDR), which must be accurately known in order to ensure that normal tissue tolerances are not exceeded. Additionally, it is necessary to scan the patient vertically through the MRT beam owing to a limited beam height, representing a time-dependent geometry that must be accounted for in MC calculations. To investigate the importance of these factors, and develop reliable MC models for MRT, a considerable amount of work has been carried out to date (Stepanek et al., 2000; Siegbahn et al., 2006; Nettelbeck et al., 2009; Spiga et al., 2007; De Felici et al., 2005; Hugtenburg et al., 2010). This work has demonstrated the feasibility of MRT and enabled the design of various beamline components. The most recent (and detailed) model of the ESRF ID17 MRT beamline was developed by Martinez-Rovira et al. (2012) using the PENELOPE MC code and validated against a comprehensive dataset of radiochromic film measurements in solid water phantoms.

A TPS for MRT has been developed (Bartzsch, 2011; Debus, 2012) in preparation for forthcoming veterinary oncology trials, and eventual clinical trials in humans. The goal of the current work is to develop an MC model of the MRT beamline to be used for independent verification of TPS dose calculations. Recent modification of the ID17 MRT beamline resulted in a need to develop an updated MC model. We describe herein the development of a MC model of the beamline based on a novel interface between the Geant4 MC toolkit and the SHADOW code, porting it to high-performance computing facilities (HPC), benchmarking against the previous MC model for the ID17 beamline, and validation against dosimetry data acquired in homogeneous water-equivalent phantoms.

2. Methods

The simulation utilizes the Geant4 MC toolkit version 4.10.00 (Agostinelli et al., 2003) interfaced with the SHADOW code (Sanchez del Rio et al., 2011). Geant4 is a C++ MC toolkit originally developed for high-energy physics applications, but has since found extensive use in radiotherapy, medical imaging and radiation protection applications. Geant4 is attractive for use in MRT for a number of reasons: it contains photon–electron physics models that have been validated for energies and materials relevant to radiotherapy applications (Thiam et al., 2008; Poon & Verhaegen, 2005); it includes an interface for importing patient computed tomography (CT) data (Paganetti et al., 2004); it is able to model the time-dependent geometries of patient scanning in MRT; and it has the ability to model complex geometries via a computer-aided design (CAD) interface (Poole et al., 2012). Rather than continuing development of the previous PENELOPE model, the beamline was remodelled using Geant4 in order to take advantage of the extensibility offered by object-oriented software design. This allows for the use of existing classes for geometry definitions (e.g radiation detectors), CT interfaces, customized physics and scoring, and minimizes the time required to extend the model by future developers. SHADOW (Sanchez del Rio et al., 2011) is a ray-tracing and X-ray optics code developed in the early 1980s that has since been used extensively for simulations in beamline optics at most synchrotron facilities. It is capable of modelling synchrotron X-ray production in the ID17 insertion device, a task not possible with Geant4 at the time of writing. We have interfaced the Geant4 model with the SHADOW code via the C++ application programming interface (API) of the latter. In this way, calls to the Fortran subroutines of SHADOW are made possible from within the C++ based Geant4 application, thereby removing the need for the production of phase space files as done in previous simulations (Martinez-Rovira et al., 2012).

The simulation is carried out in three stages: the first stage models synchrotron radiation production in the insertion device, the second transports resulting photons through the ID17 MRT beamline, and the third performs high-accuracy dose calculations in the vertically scanned phantom.

2.1. Stage I: synchrotron radiation production

The SHADOW code models synchrotron radiation production in the wiggler (Lai et al., 1988) by first calculating the trajectory followed by a single electron in the magnetic field of the wiggler. An electron in the wiggler emits radiation along its trajectory which is more intense in the region of higher curvature. The code computes the number of photons generated at each point along the trajectory and the cumulated distribution function from which a MC sampling of the photon emission along the trajectory is obtained. Once the position of emission along the trajectory is chosen, the electron beam position and velocity are computed also taking into account the electron beam emittance. A photon (ray) is created for this particular electron at this position, with the energy, divergence and electric field (polarization) calculated using the well known equation for synchrotron radiation emission in a bending magnet (Green, 1976).

The parameters used to model the ID17 insertion device are summarized in Table 1 and include: the electron beam spot-size, horizontal and vertical divergence, and energy; and the period, number thereof, and deflection parameter of the wiggler. Also shown are virtual slit settings used by SHADOW for variance reduction, defining a solid angle outside which photons are rejected. During stage I of the simulation, the Geant4 application makes a call to a SHADOW subroutine via the C++ API, which subsequently returns an array of photon phase space information containing the position, momentum vector, polarization states and energy of each photon. This phase space store (PSS) resides in memory for the subsequent stage of simulation. The use of phase space files is avoided, thereby reducing the relatively slow read–write operations to the file system and reducing computation times.

Table 1
Machine, ID17 wiggler and variance reduction used in the SHADOW calculations

Machine parameters
σy (mm)	0.057
σz (mm)	0.0103
∊y (rad mm)	3.9 × 10^−6
∊z (rad mm)	3.9 × 10^−8
E (GeV)	6.04
Wiggler parameters
λ (m)	0.15
Nperiods	11
Keff	19.5
Variance reduction parameters
Distance (m)	40.5
ΔY (mm)	60
ΔZ (mm)	4

2.2. Stage II: beamline transport

Following synchrotron radiation production in the wiggler, the beam is transported through the ID17 MRT beamline. There are two beamline configurations used at ID17 which are referred to in this study. The first, referred to herein as Condition 1, has been used for radiobiological studies and commissioning work and is described in detail in the literature (Martinez-Rovira et al., 2012). The second, referred to as Condition 2, will be used for the forthcoming veterinary oncology trials. The following description of the beamline is made with reference to Fig. 1; distances of each component from the source, d, are given in parentheses. The white beam originating from the wiggler source passes through a copper diaphragm (d = 21650 mm) with a rectangular aperture of dimensions 24 mm × 15 mm, the purpose of which is to reduce heat load on downstream elements. The beam then traverses a krypton gas filter (d = 26520 mm) with a dimension in the beam direction of 2191.5 mm. This redundant filter is present only in Condition 2; its role is to complement and protect a downstream energy filter train (Requardt et al., 2013). Following the gas filter are the copper primary slits (d = 29000 mm), the horizontal component of which can be adjusted to define the horizontal limits of the irradiation field at the patient position. The following series of carbon (1.42 mm), aluminium (0.28 mm and 1.24 mm) and copper (0.35 mm and 0.69 mm) filters (d = 29920 mm) provide redundant filtering of low energies along with the aforementioned krypton filter. Either an ionization chamber (IC) or Compton chamber (CC) beam monitor (d = 32600 mm) is then traversed prior to exiting the in vacuo part of the beamline; the purpose of these chambers is to monitor beam stability. The IC (Condition 1 only) is an air-filled device of length 140.5 mm with two Be windows of thickness 0.5 mm, whereas the CC (Condition 2 only) comprises two devices, each containing Al plates of 0.5 mm coated on both sides with 2 × 10⁻⁴ mm of Au.

Figure 1 Geometry of the beamline and phantom considered in the Geant4–SHADOW simulation, showing diaphragm (DIA); krypton filter (KF) (Condition 2 only); horizontal and vertical primary slits (PS); low-energy filters (FLT); either Compton chamber (CC) (Condition 2 only) or ionization chamber (IC) (Condition 1 only); exit windows (EXT); vertical slits (VS); multi-slit collimator (MSC); clinical ionization chamber (IC0bis); and vertically scanned water-equivalent phantom (PHAN). Dimensions not to scale.

The in-air component of the MRT facility comprises a set of tungsten slits (d = 38800 mm) which define the vertical height of the MRT beam, with a nominal aperture height of 0.5 mm. This is followed by the high-precision tungsten-alloy MSC (d = 39300 mm) (Brauer-Krisch et al., 2009) which provides the spatial fractionation required for MRT. The nominal slit width of 50 µm and distance between slits of 400 µm was considered for this work.

Following the MSC, and prior to the phantom/patient, is a pair of parallel-plate Bragg peak chambers (PTW, Freiburg, Germany) that are used to monitor beam fluence (d = 4 × 10⁴ mm). These are present for Condition 2 only and are modelled in the simulation as two discs of polymethyl methacrylate (PMMA), each of thickness 95.0 mm in the beam direction (PTW, personal communication).

Stage II of the simulation concludes with photons traversing a phase-space scoring plane situated prior to the patient or phantom. Here their position, momentum, energy and polarization information are stored in a second PSS for the final stage of simulation.

2.3. Stage III: dose calculation

The beam then impinges the patient/phantom at a distance of 40.5 m from the source. Vertical scanning of the patient is required in order to achieve tumour coverage due to the fact that the vertical beam dimension at the patient position (nominal height of 520 µm) is smaller than most tumour sizes. By exploiting the ability of Geant4 to model time-dependent geometries (Paganetti et al., 2004), this vertical scanning is considered in the simulation by modifying the phantom vertical offset as a function of particle history, which is used as a surrogate for time. Because of the phantom scanning, it is possible to recycle the PSS and reduce the time required for the dose calculation; however, care must be taken in order to avoid artefacts in the dose distribution caused by the interference of events which are no longer statistically independent (Sheikh-Bagheri & Rogers, 2002). This PSS recycling number (N_r) was optimized based on a previously published method (Kawrakow & Fippel, 2000), with a maximum simulation efficiency obtained at N_r = 500. This corresponds to the PSS being recycled every 40 µm, corresponding to 3.8% of the beam height.

For the benchmarking simulations, a homogeneous phantom of dimensions 200 mm × 300 mm × 300 mm was modelled, composed of either water-equivalent plastic or water. Dose deposition is then recorded in a Cartesian coordinate system with a variable spatial resolution. Although photons and electrons are transported throughout the entire phantom, the scoring of dose only takes place in a reduced volume of the phantom, the region of interest (ROI). Doing so has the effect of reducing simulation times, because the number of boundary crossings within the phantom is drastically reduced; moreover, scoring in a limited volume of the phantom avoids the production of very large output files. The extent of this ROI varies depending on the scenario being simulated.

High-performance computing facilities were used to improve simulation efficiency, with each simulation being executed on a separate processor and seeded with a unique seed for its respective random number generator. The mean of all dose distributions is then calculated along with an uncertainty given by the standard error within 95% confidence limits.

2.4. Physics

The Lawrence Livermore polarized physics list included with the Geant4 MC toolkit was used for this study. This list models photoelectric, Rayleigh and Compton interactions for photons. Low-energy corrections are applied to the Klein–Nishima formula to account for electron binding. Polarization effects for photoelectric and Compton interactions are also considered. Multiple Couloumb scattering, ionization interactions and bremsstrahlung production are modelled for electrons. The reader is directed to the Geant4 Physics Reference Manual for full details (GEANT4 Collaboration, 2007).

Geant4 implements a variance reduction technique known as `range cuts' in order to minimize calculation times by reducing the number of secondaries that are tracked in electromagnetic cascades. If the residual range of a secondary particle is less than the range cut value, the particle is not tracked and is assumed to deposit its energy at the point of generation, otherwise it is tracked to zero energy. A global range cut of 10 µm was used throughout the geometry for secondary photons and electrons. This corresponds to a cut-off in electron energy of 61.7 keV, 71.2 keV and 14.1 keV in copper, MSC tungsten alloy and water, respectively. In addition to range cuts, electron tracking cuts of 100 eV were implemented in the phantom ROI.

2.5. Benchmarking: Geant4 versus PENELOPE

Owing to the difficulties of performing accurate dosimetry in MRT with micrometre spatial resolution, reference dosimetry measurements are performed with a broad-beam irradiation, i.e. with no MSC present in the path of the beam. In this situation, the horizontal extent of the field size is defined by the primary horizontal slits, whereas the vertical extent is defined by the scan limits of the goniometer upon which the phantom is mounted. The ability of the code to accurately calculate relative dose profiles in homogeneous phantoms subject to irradiation by broad-beams was assessed via benchmarking against the PENELOPE calculations of Martinez-Rovira et al. (2012). In turn, these PENELOPE calculations were validated against radiochromic film measurements, full details of which can be found in the literature (Martinez-Rovira et al., 2012). As the SHADOW API does not presently provide the number of electrons incident on the insertion device, a comparison of absolute dose calculated by both codes was not possible. Instead, the comparison was made in relative dose distributions normalized to the dose at a reference depth in the phantom of 20 mm. Three-dimensional dose distributions were calculated for irradiation of a 200 mm × 300 mm × 300 mm RW3 (PTW Freiburg, Freiburg, Germany) solid water phantom by broad-beams of field sizes 10 mm × 10 mm, 20 mm × 20 mm and 30 mm × 30 mm. A scoring resolution of 2 mm × 0.2 mm × 0.2 mm was used because the dose gradients are greater in the horizontal (y) and vertical (z) directions compared with along the beam direction (x). The comparison between Geant4–SHADOW and PENELOPE calculations was made using the γ index, a hybrid distance-to-agreement/relative dose difference metric (Low & Dempsey, 2003). The γ index is the de facto method of comparing dose distributions in radiotherapy (Solberg et al., 2008; Stasi et al., 2012; Pappas et al., 2008). This test considers two dose distributions: a reference distribution D_r(r) and an evaluation distribution D_e(r), where r is a point in space. Acceptance criteria for distance-to-agreement, Δd, and dose values, ΔD, are considered. This results in an ellipsoid described by the following equation,

where δD = [D_r(r) − D_e(r)]/D_norm and the normalization value D_norm can be taken either to be the global maximum in the reference dose distribution or the local D_r(r) reference dose value. Local normalization was chosen for the current study as it is more sensitive to discrepancies in low-dose regions (Bresciani et al., 2013). These are of particular importance for MRT because the valley dose and out-of-field dose must remain below normal tissue tolerances and therefore must be accurately modelled.

For the evaluated distribution to agree with the reference distribution, there must be at least one point lying within this ellipsoid, i.e. one point for which

The γ test is particularly useful for distributions with high gradients, as is the case in MRT, for which a very small spatial shift can result in a very high relative difference in dose.

Benchmarking against the PENELOPE code was also carried out for microbeam irradiations of a 200 mm × 300 mm × 300 mm RW3 solid water phantom in order to determine the ability of the MC code to accurately calculate PVDRs in water. This benchmarking was done for field sizes of 10 mm × 10 mm, 20 mm × 20 mm and 30 mm × 30 mm. A comparison was made using calculated central microbeam depth-dose distributions and central microbeam PVDRs as a function of depth, along with radiochromic film measurements of both as sourced from the literature (Martinez-Rovira et al., 2012). A reduced ROI was used for scoring in order to minimize calculation times, the dimensions of which were 200 mm × 2 mm × 2 mm, with a scoring resolution of 1 mm × 0.05 mm × 0.1 mm.

2.6. Validation

Additional experimental validation was carried out by simulating depth-dose distributions in a 300 mm × 300 mm × 300 mm water phantom, irradiated by a 20 mm × 20 mm broad-beam field, for both Condition 1 and Condition 2. These were calculated using a scoring resolution of 2 mm × 2 mm × 2 mm and compared with those measured using a PinPoint (PTW Freiburg, Freiburg, Germany) ionization chamber.

2.7. Investigation of polarization effects

Due to the highly polarized nature of the MRT beam, predominantly in the electron orbital plane of the synchrotron, photons undergoing Compton interactions in the phantom/patient will be preferentially scattered in the vertical direction. As a consequence, the out-of-field dose is greater in the vertical direction than in the horizontal direction for the same offset (Hugtenburg et al., 2010). The ability to predict this out-of-field dose asymmetry is of great importance to MRT because it will ensure that organs-at-risk are not exposed to doses above the normal tissue tolerance. In order to assess the accuracy of the MC model in predicting this out-of-field dose asymmetry, the dose distribution in an RW3 phantom irradiated by a 20 mm × 20 mm reference field was simulated (Condition 1 only). Since a high spatial resolution is not required for out-of-field dose calculations, a scoring resolution of 2 mm × 2 mm × 2 mm with an extended ROI of dimensions 150 mm × 150 mm × 200 mm was considered. The calculated vertical to horizontal out-of-field dose ratios were compared with those measured by EBT2 (Internationl Specialty Products, Wayne, NJ, USA) radiochromic film under the same conditions.

Preparation and analysis of EBT2 film were carried out based on a previously published protocol (Hartmann et al., 2010). Calibration was performed using 30 mm × 30 mm pieces from a single sheet placed at a depth of 20 mm in a RW3 solid water phantom and irradiated with a 20 mm × 20 mm broad-beam to doses of 30, 50, 75, 100, 200 and 300 Gy. Measurement pieces, taken from the same sheet and with dimensions of approximately 100 mm × 50 mm, were placed at depths of 20 mm and 40 mm in a solid water phantom and irradiated with a 20 mm × 20 mm field. The total dose delivered at the reference depth was approximately 40 kGy; this was required to ensure that the out-of-field dose fell within the calibration range. Readout was performed six days post-irradiation using an EPSON Perfection V700 scanner, 72 d.p.i. resolution (equivalent spatial resolution of 0.35 mm) and 48-bit colour. ImageJ was then used to isolate the red channel of the images, and export these pixel value images as text files. An automated Python script was subsequently used to apply the dose calibration and sample the out-of-field dose at distances in the vertical and horizontal direction. This dose was calculated using the average pixel value over a 1 mm² area and applying the calibration curve. Uncertainties were calculated using the standard error of these readings, within 95% confidence intervals. The experiment was repeated and results averaged over the two sets of measurements.

3. Results

Fig. 2 shows a plot of Geant4 and PENELOPE simulated depth-dose profiles in an RW3 water-equivalent phantom irradiated by broad-beams of given field sizes. The total simulation time for each field was approximately 5 h running in parallel on one-hundred 2.67 GHz cores [Intel(R) Xeon(R) X5650 CPU]. The percentage of total time spent for stages I, II and III were 2%, 2% and 96%, respectively. Doses were normalized to a depth of 20 mm. Also shown are γ indices for distance-to-agreement/relative dose difference criteria of 3 mm/3%; points below a value of unity indicate agreement between the two codes. There was excellent agreement, with 97%, 99% and 94% of points passing the γ test within a depth of 150 mm for 10 mm × 10 mm, 20 mm × 20 mm and 30 mm × 30 mm fields, respectively.

Figure 2 Comparison between Geant4- and PENELOPE-calculated depth-dose profiles for broad-beam irradiation of an RW3 water-equivalent phantom for given field sizes. The γ indices for criteria of 3 mm/3% are shown as crosses.

Fig. 3 presents the corresponding horizontal dose profiles. For all field sizes the out-of field dose differed by more than 0.3 mm/3% between the two codes, with Geant4 results lying closer to the experimentally measured values. For the in-field dose component, the results for the 20 mm × 20 mm and 30 mm × 30 mm broad-beam showed significant difference between Geant4 and PENELOPE toward the edge of the field; once more, Geant4 results were closer to experimental values. Although both Geant4 and PENELOPE simulations account for the polarization of the source, the latter used an approximate source model for which position, direction and energy of photons were sampled from parametric probability distributions based on phase space files scored prior to the phantom. These distributions included approximations such as sampling the horizontal position of photons with a uniform probability over the field size (Martinez-Rovira et al., 2012). In order to investigate whether such assumptions were the cause of the discrepancy observed in the current study, identical assumptions were made in the Geant4–SHADOW code. Results are shown in Fig. 4, with agreement in the horizontal dose profiles for the 30 mm × 30 mm field size once these approximations were implemented.

Figure 3 Comparison between Geant4-simulated, PENELOPE-simulated and experimentally measured horizontal dose profiles at 20 mm depth in an RW3 water-equivalent phantom irradiated by given broad-beam field sizes. The γ indices for criteria of 0.3 mm/3% are shown as crosses.

Figure 4 Comparison between Geant4- and PENELOPE-simulated horizontal dose profiles at 20 mm depth in an RW3 water-equivalent phantom irradiated by a 30 mm × 30 mm broad-beam. Geant4 results were obtained with the same source parameterization used in the PENELOPE beam model.

Agreement was obtained in-field for the vertical profiles given in Fig. 5, with disagreement in out-of-field doses. Contrary to the horizontal profiles, PENELOPE results were closer to the film measurements than Geant4 results for the out-of-field doses.

Figure 5 Comparison between Geant4-simulated, PENELOPE-simulated and experimentally measured vertical dose profiles at 20 mm depth in an RW3 water-equivalent phantom irradiated by given broad-beam field sizes. The γ indices for criteria of 0.3 mm/3% are shown as crosses.

Fig. 6 presents a sample of Geant4-calculated horizontal dose profiles for MRT microbeam irradiation of a water-equivalent phantom by a field size of 10 mm × 10 mm, with results normalized to the peak dose at 3 mm. The total simulation time for all fields was approximately 10 h running in parallel on one-hundred 2.67 GHz cores [Intel(R) Xeon(R) CPU X5650]. The effect of the ROI on computation time was marked for this case, with a 34-fold decrease compared with the case of scoring in the entire phantom. The decrease in peak dose with depth in the phantom can be seen, as can the build-up of valley dose between 3 mm and 20 mm depth, followed by the reduction with depth in accordance with attenuation of the primary beam. Fig. 7 presents central microbeam peak depth-dose (PDD) distributions for various field sizes as calculated by Geant4 and PENELOPE, along with experimental measurements obtained with radiochromic film. Simulation and experimental results all agreed within uncertainties of the experimental results. However, relative experimental uncertainties were up to 14.7%, highlighting the need for improvement of the precision of dosimetry in MRT. The valley depth-dose (VDD) for the same simulation is shown in Fig. 8. The valley dose initially increased with depth, followed by reduction with depth in accordance with attenuation of the primary beam. Moreover, the valley dose increased with an increase in field size due to the concomitant increase in scattered photons in the phantom. Both codes agreed with each other, and with the experimental data, within uncertainties.

Figure 6 Horizontal dose profiles of a microbeam array at given depths in a water-equivalent phantom irradiated by a 30 mm × 30 mm field; calculated by the Geant4–SHADOW code.

Figure 7 Geant4- and PENELOPE-calculated and experimental depth-dose profiles of the central microbeam for microbeam irradiation of an RW3 water-equivalent phantom for given field sizes.

Figure 8 Geant4- and PENELOPE-calculated and experimental valley depth-dose profiles of the central microbeam for irradiation of a water-equivalent phantom at given field sizes.

Results of the corresponding central microbeam PVDRs are given in Fig. 9. The PVDR decreased with depth owing to the increase in valley dose cause by photons scattered out of the primary beam. At the same depth, the PVDR decreased with an increase in field size due to the increased contribution from scattered photons. Both codes agreed with experimental measurements within the experimental uncertainty (up to 28% for the 10 mm × 10 mm field size). Once more, higher-precision dosimetry is required to conclusively determine which of the two simulations was more accurate in calculating PVDRs.

Figure 9 Geant4- and PENELOPE-calculated and experimental PVDRs of the central microbeam for microbeam irradiation of an RW3 water-equivalent phantom for given field sizes.

Fig. 10 shows depth-dose profiles in a water phantom subject to irradiation by the 20 mm × 20 mm broad-beam, normalized to dose at 20 mm depth. Slight beam-hardening was observed for Condition 1 owing to greater attenuation of the low-energy component by the additional beam-monitoring elements, namely the CC and Bragg peak chambers. Geant4-simulated and experimental results agreed within 2.3% for both Condition 1 and Condition 2.

Figure 10 Percentage depth-dose as calculated by Geant4–SHADOW code and measured by a PinPoint IC in a water phantom irradiated by a 20 mm × 20 mm reference field for both Condition 1 and Condition 2 beamline configurations.

The asymmetry of out-of-field doses in the vertical and horizontal direction is presented in Fig. 11, as calculated by Geant4 and measured using EBT2 radiochromic film. At a given depth, the dose asymmetry increased with increasing distance from the field edge. With increasing depth in the phantom, the asymmetry reduced, presumably due to the increasing proportion of scattered photons having randomly oriented polarization vectors. The relative differences between experiment and simulation (relative to simulation) were within 3.2% at 20 mm depth and 11.5% at 40 mm depth.

Figure 11 The ratio of vertical-to-horizontal dose as a function of distance from the beam axis for depths of 20 mm and 40 mm in an RW3 water-equivalent phantom; calculated by the Geant4–SHADOW code and measured using EBT2 radiochromic film.

4. Discussion

Good agreement was obtained between the Geant4–SHADOW and PENELOPE simulations of broad-beam irradiation of an RW3 phantom when comparing depth-dose profiles within a depth of 150 mm, as well as horizontal and vertical dose profiles. For larger field sizes, there was discrepancy between codes toward the edge of the field for horizontal profiles, with Geant4 results lying closer to film measurements. This discrepancy was shown to be due to the parameterized source model used in the PENELOPE simulation, demonstrating the importance of detailed MC source modelling in MRT when high-accuracy dose calculations are required.

Additional experimental validation of the Geant4–SHADOW code was carried out by simulating depth-dose profiles in a water phantom and comparing with those measured with a calibrated IC under different beamline configurations, for which agreement within 2.3% was observed. This agreement confirms the ability of the code to accurately predict relative depth-dose distributions in homogeneous phantoms subject to irradiation with the MRT broad-beam.

Notwithstanding the high uncertainties of the film measurements of central microbeam PDDs and PVDRs, agreement was seen between both codes and experimental results. The lack of precision in experimental data prevents a definitive assessment of which code is more accurate in the calculation of microbeam parameters, once again highlighting the need for higher precision dosimetry in MRT. Silicon microstrip detectors offer the ability to provide high-precision measurements of valley doses (Petasecca et al., 2012); however, energy dependence at low energies is pronounced (Cheung et al., 2009). This energy dependence could lead to an over-response in valley dose due to the dominance of lower-energy photons scattered out of the primary beam (Siegbahn et al., 2006). Such an artefact could be mitigated by incorporating the detector geometry and composition directly into the MC model (Rosenfeld et al., 2006). A direct comparison between PVDR values obtained in the present study and those published in the literature is difficult owing to significantly different beam conditions (different beam energies, field sizes and phantom materials) used for each.

The accuracy of the Geant4–SHADOW code in predicting out-of-field dose asymmetry in MRT was also determined. The importance of polarization has been the topic of discussion in the MRT research community; its effect was found to be negligible at the centre of a 30 mm × 30 mm microbeam field, yet increasing to 16% at the field edge for a 200 keV monoenergetic beam (De Felici et al., 2005). In a separate study (Hugtenburg et al., 2010), clinically relevant dose volume histograms were calculated for an organ-at-risk in a hypothetical MRT treatment, in which case the secondary scatter was found to have a significiant effect. The current findings support the notion that the effect of beam polarization on out-of-field doses must be taken into consideration in any TPS system for MRT. Based on these results, ignoring asymmetry in the out-of-field dose profiles could lead to errors of up to 150%. Since an intended application is for paediatric patients suffering glioblastomas, the importance in the context of mitigating the risk of secondary malignancies is amplified (Newhauser & Durante, 2011). The present study has demonstrated the ability of Geant4 to predict the dose asymmetry resulting from beam polarization, albeit with limited precision and accuracy. The challenge remains to develop a more precise dosimetry system in conjunction with refining the theoretical models of Compton scattering present in Geant4 to obtain very high accuracy in out-of-field dose calculations.

The deflection parameter used in the current work (K = 19.5) is markedly lower than that used in the work of Martinez-Rovira et al. (2012) (K = 22.3). The deflection parameter of ID17 is based on a parameterized fit of peak magnetic (B_o) field versus wiggler gap, where the magnetic field is assumed to be sinusoidally varying along the length of the wiggler. In reality, the wiggler comprises a total of 21 poles including 17 internal poles with a peak field B_o (described by the aforementioned parameterized formula), two intermediate poles at 0.86B_o and two extremity poles at 0.53B_o (Chavanne, 1998). This uncertainty in the deflection parameter supports its use as a quasi-free parameter that may be tuned in order to obtain agreement between experimental results and the Geant4–SHADOW code.

The clinical implications of this research lie in the need for independent pre-treatment plan verification tools in MRT (as for any radiotherapy modality) as recommended by the IAEA Technical Report Series No. 430 (Sharpe, 2006). MC represents the most accurate method of calculating dose in patient geometries, and we have developed a tool for dose calculations in homogeneous phantoms. This tool has been benchmarked against an MC model of the previous beamline configuration and validated against experimental data acquired with both beamline Condition 1 and 2. It enables event-by-event simulation from photon production in the wiggler through to dose deposition in the phantom, within time constraints deemed suitable for a radiotherapy modality in the pre-clinical stage of development. In its current form, the simulation will be able to provide accurate dose estimates for radiobiological experiments performed on the beamline. It may also be used in the optimization of radiation detector design by incorporating relevant composition and geometry into the simulation (Othman et al., 2010; Lian et al., 2011).

Prior to use as a TPS verification tool, this application will need to be validated by simulating experimental data acquired in heterogenous phantoms using the CT interface of Geant4. The increase in simulation time is not expected to be excessive. Dose calculations in the CT dataset will involve a similar number of scoring voxels as the current simulations. It is expected that simulation times of 10 h running on 100 cores of a high-performance computing facility will provide acceptable statistics. This will be acceptable for the purpose of pre-treatment plan verification in the early clinical stages of MRT, for which patient turnover will be low. Experimental validation will be carried out by irradiation of film, silicon devices and ionization chambers embedded in an anthropomorphic phantom. The current study used normalized dose distributions to benchmark and validate the code; future studies will be conducted to determine the ability of the code to determine absolute dose values in MRT. Moreover, future studies will compare full two-dimensional and three-dimensional dose distributions of microbeam irradiations, rather than the one-dimensional profiles used in the present study.

Comparison of dose distributions was performed using the γ index method. Use of the γ index for comparing dose distributions in MRT is unprecedented, with previous authors implementing a variety of metrics. Martinez-Rovira et al. (2012) implemented a piece-wise relative dose difference metric with tolerance varying depending on whether the dose points are determined to be in-field, in the beam penumbra, or out of the treatment field (Venselaar et al., 2001). Although all such metrics are useful in validating dosimetry calculations in phantoms, their clinical relevance has recently been brought into question (Nelms et al., 2011; Stasi et al., 2012); i.e. the ability of a dose calculation engine to predict one-dimensional and two-dimensional dose distributions in phantoms has little correlation with its ability to accurately predict dose to anatomical structures. A more clinically relevant metric, such as a dose volume histogram for target volumes and organs at risk, should be used to perform routine quality assurance of patient treatment plans.

5. Conclusions

We have developed an up-to-date MC model of the ID17 MRT beamline based on a novel interface between the Geant4 toolkit and SHADOW code, allowing for event-by-event modelling from wiggler to phantom. The MC code was ported to HPC facilities and benchmarked by comparison with PENELOPE simulations. Additional experimental validation of the code was carried out via comparison against IC measurements in a water phantom, as well as a study of the out-of-field dose asymmetry in a water-equivalent phantom subject to broad-beam irradiation. This interface removes the need for parameterized source models based on phase space files, leading to more accurate dose calculations. The time-dependent geometry of the application will allow for the simulation of novel time-resolved dosimetry systems and patient scanning in MRT. The accuracy of Geant4 in the calculation of dose in MRT was demonstrated, paving the way for the creation of an independent treatment plan verification tool.