Abstract
Background/Aim: The Purpose of this study was to develop a Monte Carlo (MC) model for the Agility multileaf collimator (MLC) mounted and to validate its accuracy. Materials and Methods: To describe the Agility MLC in the BEAMnrc MC code, an existing component module code was modified to include its characteristics. The leaf characterization of the MC model was validated by comparing the calculated interleaf transmission and tongue-and-groove effect with EBT2 film and diode measurements and IMRT and VMAT calculations with film measurements. Results: Agreement between mean calculated and measured leaf transmissions was within 0.1%. The discrepancy between MC calculation and measurement in a static irregular field was less than 2%/2 mm. Gamma analysis of the comparison of MC and EBT2 film measurements in IMRT and VMAT fields yielded pass rates of 99.1% and 99.5% with 3%/3 mm criteria, respectively. Conclusion: Our findings demonstrate the accuracy of the MC model using an adapted BEAMnrc component module for the Agility MLC.
Intensity modulated radiation therapy (IMRT) and volumetric modulated arc therapy (VMAT) have been introduced in clinical radiation oncology to improve the control of the tumor by delivering a high dose to the tumor and/or to reduce the risk of normal tissue injury (1). These techniques are commonly implemented by means of multileaf collimators (MLCs). To achieve optimal dose distribution delivery, it is important that the characteristics of the MLCs are well-designed, and new MLCs have been developed by various manufacturers with specific design characteristics for each MLC (2, 3). The Agility MLC is a recent MLC produced by Elekta AB (Stockholm, Sweden) and consists of 160 interdigitating tungsten alloy leaves with a width projected to the isocenter of 5 mm and a maximum field size at the isocenter of 40×40 cm2 (4). The smaller leaf width and faster leaf speed of the Agility MLC as compared to more conventional MLCs are characteristics designed to deliver a better modulated dose distribution with advanced techniques such as IMRT and VMAT. Since using these advanced techniques is complicated, accurate commissioning for both treatment planning and quality assurance (QA) is essential (5, 6). QA for VMAT is a recent area of active research because VMAT dose delivery uses the most complex technique in modern radiation therapy including continuous variations in gantry speed and dose rate during delivery (7, 8). The temporal description of these correlated variations must therefore be an integral part of the modeling.
Consequently, it is necessary to make it possible for Monte Carlo (MC) calculations to be used with this most recent MLC as a verification tool for IMRT and VMAT radiation treatments. MC simulation is recognized as the most accurate dose calculation algorithm currently available as it can take into account all relevant aspects of photon and electron transport even in heterogeneous situations. In addition to measurement-based validation for IMRT and VMAT dose distributions, several investigators have reported that MC simulations could be essential for patient specific QA or treatment planning system (TPS) commissioning (9, 10). Complete and detailed models of linear accelerator head components such as target, flattening filter, MLC and jaws are required for these simulations. Accurate modeling of the MLC is especially necessary to fully account for specific effects such as tongue-and-groove and leaf transmission and their impact on dose distributions in patients. In the case of IMRT and VMAT, these effects are not limited to the outer field edges but also affect the in-field dose.
BEAMnrc (11) is the most widely used MC code to simulate radiation transport through megavoltage linear accelerators and makes it possible for the treatment head of the linear accelerator to be modeled by independent component modules. Several of these modules are already available for modeling various designs of MLCs [MLCE (12), MLCQ (13), VARMLC (14) and DYNVMLC (15)]. Several investigators have tried to model other MLC designs by using existing component modules with modifications or approximations (15, 16) made with the BEAMnrc code. Due to its unique design, no component module is available yet that can directly simulate the Agility MLC.
The purpose of this study was to develop an MC model of the Agility MLC as part of an Elekta Synergy linear accelerator with BEAMnrc based on an existing component module, to include its temporal behavior and to validate its accuracy in static and intensity modulated (IMRT and VMAT) 6 MV fields by comparison with measurements.
Materials and Methods
Accelerator head. Figure 1 shows a schematic overview of the linear accelerator equipped with the Agility MLC (both from Elekta AB). The following independent treatment components were modeled: target, primary collimator, flattening filter, ionization chamber, back scatter plate and mirror. The MLC and jaws are parts of the linear accelerator that are included depending on the treatment. In contrast to the Elekta SLi series, no upper jaws (back-up jaws) are placed under the MLC with the Elekta Synergy S. The MLC leaves travel along the X-axis and the jaws travel along the Y-axis in accordance with the IEC 61217 coordinate system.
Agility MLC. The Agility MLC consists of 160 interdigitating tungsten alloy leaves with each leaf projecting to 5 mm width at the isocenter. Maximum field size is 40×40 cm2 at the isocenter. Each leaf bank is connected to a dynamic leaf guide and the combined speed of these two components can reach 6.5 cm per second–approximately twice as fast as other commonly used MLCs. For a satisfactory approximation, the sides of the leaves can be assumed to be flat and gaps between the leaves are effectively dealt with through tilting the leaves by shifting their focal point to reduce inter-leaf leakage radiation (Figure 2).
Monte Carlo simulations. The BEAMnrc/EGSnrc (17) (version 4-2-3-1) code was used to simulate radiation transport through the linear accelerator head from the top edge of the target to the bottom edge of the jaws. Dimensions, materials and densities for each component module were in accordance with the specifications provided by the manufacturer. At the bottom of the linac head, i.e. below the jaws, phase space data (PSD) was scored which contains information about the position, direction and energy for each particle. Subsequently, the PSD was used as input for DOSXYZnrc/EGSnrc (18) to calculate the dose for a patient or phantom. These calculations were performed on a workstation at the University Medical Center Groningen (UMCG), which consists of 12 cores, and also on the Millipede cluster at the University of Groningen, which consists of 252 multi-core nodes. For both BEAMnrc and DOSXYZnrc, the cutoff energies of ECUT and PCUT used in this study were set at 0.7 MeV and 0.01 MeV, respectively. As a rule, more than 10 million particles were stored in the PSD files for each simulation.
Schematic overview with IEC 61217 Y-Z view of the Elekta Synegy linear accelerator head with Agility MLC.
Modeling of the MLC. To model the characteristics of the Agility MLC, we found that the existing original VARMLC component module code was sufficiently universal to be adapted for our purpose. Because the original VARMLC component module models the leaf sides divergent with respect to the radiation source, we changed the focal point of the leaf sides so that they could be focused on a virtual source instead of the real one. The density and composition of the tungsten of the leaves were set to 18.0 g cm−3 and 95% tungsten, 3.75% nickel and 1.25% iron, respectively, to agree with the vendor's specifications. The leaf tip was round with a radius of 17.0 cm and the leaf thickness was 9.0 cm. The interleaf air gaps were set to 0.009 cm. Dynamic leaf guides were not considered for this model.
After the Agility MLC was modeled, its geometry was visualized by means of a ray tracing method which traces particles and records their position each time a leaf boundary is crossed. Subsequently, these positions were plotted in a 2D projection and checked visually.
MC simulation for IMRT and VMAT. To simulate the motion of MLC leaves and jaws for IMRT and VMAT fields, the BEAMnrc code was modified to sample the leaf positions for each particle history as described by Liu et al. (19). With their method, a random number between 0 and 1 is assigned to a particle when it is initiated in a simulation. The leaf positions between two consecutive segments are defined by linear interpolation using a random number on the assumption that each leaf moves at a constant speed within each segment. For example, the random numbers of 0 and 1 represent the start and end positions, respectively, of the leaves and jaws in the entire field. The sequence files which include the MU index and leaf and jaws positions for each segment from a treatment planning system were converted to input files of BEAMnrc by using in-house software (MATLAB R2012a; MathWorks, Natick, MA, USA).
Schematic drawing of the leaf tilt design. Leaf sides are focused on a virtual source which is shifted from the radiation source perpendicular to the leafs' travel direction to reduce inter-leaf leakage.
BEAMnrc and DOSXYZnrc are independent tools for a linear accelerator head and a phantom simulation, respectively. To synchronize the motion of leaves and jaws in BEAMnrc with gantry motion in DOSXYZnrc, which is necessary for VMAT calculations, the random numbers used for the BEAMnrc simulation were recorded in an extension of the PSD file in addition to information of position, direction and energy for each particle (20). Subsequently, these random numbers were read in during the DOSXYZnrc simulations. Since the gantry angle was also controlled by the MU index, these parameters were defined for each particle with the same random numbers as for BEAMnrc on the assumption that the gantry and collimator move at a constant speed within each segment. MC simulations were performed for each arc in VMAT followed by the accumulation of segment dose distributions to acquire the total dose distribution.
Validation of the implemented MC model. The validation of the implemented linear accelerator model was performed in two steps in this study. First, we validated the head model by comparing MC calculations with measurements for open square fields. The electron beam energy spectrum at the vacuum exit window and the full width at half maximum (FWHM) of the radial intensity profile were validated by comparing simulations with measurements in relative depth dose profiles and lateral dose profiles for various open square fields, i.e., 5×5, 10×10, 20×20 and 40×40 cm2. The voxel size was set at 0.5×0.5×0.5 cm3 for these MC simulations.
Second, we validated the accuracy of the MLC model including MLC leakage and intensity modulation by comparison of MC calculations with measurements in closed, irregular and intensity modulated fields. Leaf transmission was studied with a field with all leaves ‘closed’ at a 15 cm offset from the central axis and jaws fully retracted. In line with the minimal requirements for our clinical configuration, in this ‘closed’ state the leaf tips were still 6 mm open at the isocenter plane. To minimize interaction of the phantom scatter from radiation through the leaf gap with the phantom or treatment couch, the phantom edge was offset from the isocenter which resulted in a minimal distance of 2 cm from the projected gap to the phantom edge and treatment couch. In an MC simulation, a completely closed leaf gap was also simulated to evaluate pure leaf transmission without scatter dose from the leaf gap. Dose profiles were calculated in a phantom and normalized to the central point dose of a 10×10 cm2 field. To validate modulation of the leaf edges including the tongue-and-groove effect, dose distributions for an irregularly blocked field were determined. The voxel sizes were set at 0.1×0.5×0.5 cm3 for these simulations.
Step-and-shoot IMRT and VMAT fields were simulated in order to validate the overall accuracy of the implemented model. A seven-field step-and-shoot IMRT and a dual-arc VMAT treatment plan for head and neck cancer were newly created with the aid of our clinical treatment planning system (Pinnacle version 9.0; Philips Medical Systems, Fitchburg, WI, USA). The voxel size was set at 0.2×0.2×0.2 cm3 for these MC simulations.
We ensured that statistical uncertainties were within 2.0% for all simulations.
Measurements. Measurements were performed using an ionization chamber and a diode detector in a water phantom and films in a solid water phantom.
A thimble-type ionization chamber (CC04; IBA Dosimetry, Schwarzenbruck, Germany) with a sensitive volume of 0.04 cm3 was used to measure the relative depth and lateral dose profiles in the water phantom. These profiles were measured for 5×5, 10×10, 20×20 and 40×40 cm2 open square fields at 90 cm source surface distance (SSD) and at 10 cm depth. Although measuring penumbra regions with a CC04 introduces effects associated with detector size, we considered this effect acceptable since a comparable effect occurred in the simulation due to the chosen voxel size. This chamber can be considered as a reference chamber for a field as small as 2×2 cm2 (21).
Since irregularly blocked fields may produce very steep dose distributions, a diode detector (EFD-3G: IBA Dosimetry AB, Uppsala, Sweden) was used for measurement. The profile was measured at 90 cm SSD and at 10 cm depth.
Radiocromic films (Gafchromic EBT2; Ashland Specialty Ingredients, Wilmington, DE, USA) were used for film measurements placed in a coronal orientation in a 30×30×10 cm3 solid water phantom (Gammex-RMI, Middleton WI, USA). Film measurements for closed field and IMRT-VMAT fields were performed at 98.5 cm SSD and 1.5 cm depth and at 95 cm SSD and 5 cm depth, respectively.
The day after irradiation, the films were scanned on a commercial flatbed scanner (V700; Epson America Inc., Long Beach, CA, USA) with a resolution of 75 dots/inch. The method used for converting optical density to dose has been described in recent publications (22, 23) and implemented in an in-house modification of Doselab (24). The dose calculated by means of MC and the dose as measured by films was interpolated from a 0.2 cm and 0.034 cm resolution, respectively, to a 1 mm grid size. A gamma analysis with 3%/3 mm criteria and a 50% hot area dose threshold was used for both step-and-shoot IMRT and VMAT fields. Dose profiles were also analyzed for comparison of MC calculation with film measurements.
Comparison of the MC calculations with ionization chamber measurements for 5, 10, 20 and 40 cm2 square open fields for 6 MV. (a): depth dose profiles, lateral dose profiles at 10 cm depth along (b) the X-axis (parallel to MLC travel direction) and (c) the Y-axis (perpendicular to MLC travel direction). The profiles were normalized to the central dose of a 10×10 cm2 field at 10 cm depth. The lines correspond to the measurements and the squares to the MC calculations.
Ray tracing of the central 4 leaves modeled in BEAMnrc. The particles were traced and their positions recorded when they crossed the leaf boundaries. The leaves in Y-Z view were visualized by plotting their outline.
Results
Square open fields. The parameters of the mean energy of the electron spectrum hitting the Bremsstrahlung target and the FWHM of the radial intensity distribution in MC simulations were tuned to get best fitted depth dose and lateral behavior. The values of 6.0 MeV for the first parameter and 0.15 cm, for the second parameter showed the best agreement with measurements. Figure 3a shows the results of a comparison between the relative dose profiles of MC calculations and ionization chamber measurements for various square fields at an SSD of 90 cm. Each profile was normalized to the central dose of the 10×10 cm2 field at 10 cm depth. The MC calculations of the relative dose profiles for each of the fields corresponded well with the measured curves. Discrepancies beyond the build-up regions were within 1.5%.
Lateral dose profiles along the X-axis and the Y-axis for each field are also shown in Figure 3. Each of these profiles was determined at an SSD of 90 cm at 10 cm depth and normalized to the central dose of the 10×10 cm2 field. The discrepancies were within 2.5% for each field in both the X and Y direction. The correspondence between calculations and measurements in the penumbra regions was considered to be an indicator of the accuracy with which the rounded leaves and jaws tip could be modeled because open field edges in the X and Y direction were defined by leaf tips and jaws, respectively, as shown in Figure 1.
Comparison of relative film measurements with MC calculations at 1.5 cm depth in a closed leaf field with a 6 mm gap remaining (upper curves). MC calculation of a fully closed leaf field (lower curve) is also included for comparison. All profiles were normalized to the central axis point dose in a 10×10 cm2 open field at the same depth.
Ray tracing. To validate the MLC geometry modeled in BEAMnrc, the ray tracing method was employed. This method allows for a direct visualization of the geometry and it provides a first proof of the correct modeling of the MLC. A macro in the code allows particles to be traced and their positions to be recorded each time they cross a leaf boundary. These positions were then plotted in a 2D projection in the Y-Z view as shown together with the central four leaves in Figure 4. The leaf tilt design is visible from the leaf sides that were focused not on the radiation source located at Y=0 but on the virtual source. The distance between radiation source and virtual source was set at 3.25 mm. It can be seen in this figure that all traced positions were within the 2D projected boundaries of traversing leafs.
Leaf transmission. A comparison between composite leaf transmission, leaf reflection and interleaf leakage of MC calculations and film measurements is shown in Figure 5. Agreement between mean calculated and measured leaf transmissions with a 6 mm leaf gap remaining and fully opened jaws, normalized to the center dose of a 10×10 cm2 field at the same depth, was within 0.1%. It was found that the positions of the calculated and measured profile peaks showed good correspondence. The repetition distance of 5 mm corresponds to the leaf width. When the gap was completely closed, leakage was lower than 0.5%.
Comparison of the MC calculated dose profiles with diode measurements at SSD of 90 cm and at 10 cm depth for an irregular field. Both profiles were normalized to the central point dose of a 10×10 cm2 open square field at the same depth. The blue lines correspond to the measurements and the red squares to the MC calculations.
Irregular field. An irregular MLC field with alternating sections of closed and opened leaves was used to assess the tongue-and-groove effect. This irregular field was also found to be sensitive to the leaf tilt angle because the leaf width projected onto the isocenter depends on this angle. Figure 6 shows comparisons of MC calculated and measured profiles through the central axis, demonstrating a good agreement with discrepancies below 2%/2 mm.
IMRT and VMAT fields. To validate the overall accuracy of the MLC model, a seven-field step-and-shoot IMRT and a dual full arc VMAT plan for head-and-neck cancer were calculated and measured. Both calculations and measurements were performed in the center of a 30×30×10 cm3 solid water phantom. The Monte Carlo calculations were then normalized to obtain optimal overall agreement to measurements.
Figure 7 shows comparisons between the dose profiles and the gamma index plots of MC calculation and film measurement for the IMRT treatment plan. The discrepancies were within 3% for the profiles at the isocenter along the Y-axis. The gamma analysis showed a 99.1% pass rate with 3%/3 mm criteria. The deviations above gamma index1 are shown in color.
Figure 8 shows a comparison between the dose profiles and the gamma index plots of MC calculations and measurements for the VMAT treatment plan. For this complex plan, too, the calculation showed good agreement with film measurement. The discrepancies were within 3% for the profiles at the isocenter along the Y-axis. The gamma analysis shows a 99.5% pass rate with 3%/3 mm criteria.
Discussion
MLCs have become essential components of linear accelerators for modern radiation therapy to deliver treatment plans for high-dose concentration on tumors while avoiding surrounding normal tissues as much as possible. It can be said that the quality of the dose distribution that can be achieved depends on MLCs specifications such as leaf width, leaf speed, maximum field size and leakage. New MLCs with sophisticated designs have therefore been developed recently to better satisfy these conditions.
The MC model of the Agility MLC mounted on the Synergy linear accelerator was developed based on an existing component module and the study reported here has validated its accuracy since our results showed good agreement with measurements in not only static fields but also intensity modulated fields.
In general, good agreements between the MC calculations and measurements using an ionization chamber were obtained for the square fields for various field sizes. The discrepancies between the relative dose profiles beyond the build-up regions and the lateral dose profiles were within 1.5% and 2.5%, respectively for each field size. The International Commission on Radiation Units and Measurements (25) recommended that the overall dose delivery accuracy should be within 5% including many considerations of uncertainty, for example, the accuracy of computed dose distributions should be less than 2%. Over the years, much effort has been devoted to MC simulation of different types of linear accelerators (26, 27). These also showed good agreements within 2% and it was concluded that the key parameters for matching the MC simulations to the measurements were the mean energy and the radial intensity distribution of the incident electron beam. We were able to show that both parameters were well tuned for our accelerator evaluated in the current study.
Even without the back-up jaws that characterized the previous MLC design, the leaf transmission of our accelerator was very low, less than 0.7%, and thus about two times lower than that of the previous type. A previous study (28) reported a leaf transmission of the Agility MLC of less than 0.5%. Our findings showed a slightly higher leaf transmission, most probably because of the inclusion of additional radiation through the leaf tip gap and reflection at the leaf tips that represents the minimum gap size in our clinical configuration, as was illustrated by our calculation of MC with artificially completely closed leafs.
Comparison of MC calculations and measurements of (a) the dose profiles at isocenter (position indicated by the dashed line in the insert) and (b) the gamma index plots for a step-and-shoot 7-field head-and-neck IMRT treatment plan. The line corresponds to the film measurement and the red squares to the MC calculations. In the gamma index plots, the deviations above gamma index=1 (3%/3 mm) are shown in color.
Comparison of MC calculations and measurements of (a) the dose profiles at isocenter (position indicated by the dashed line in the insert) and (b) the gamma index plots for a two-arc head-and-neck VMAT treatment plan. The line corresponds to the film measurement and the red squares to the MC calculations. In the gamma index plots, the deviations above gamma index=1 (3%/3 mm) are shown in color.
The DYNVMLC component module designed to model the Varian Millennium MLC can simulate the step-and-shoot IMRT field using a linear interpolation method (29). We simulated the IMRT technique with the modified VARMLC component module using the same method. The gamma analysis showed a 99.1% pass rate with 3%/3 mm criteria. VMAT as a treatment technique was developed to improve therapeutic efficiency with continuous motion of MLC and dose rate modulation during gantry rotation and it is increasingly being used for head and neck cancer (30, 31). VMAT planning for head and neck is one of the most complex techniques in modern radiation therapy that allows for large field size and simultaneously high dose gradients. The importance of QA for VMAT dose delivery has been discussed and several investigators reported that MC simulations could be used as independent verification tools (32). However, these methods approximated the gantry and MLC motions by discretization of the effective gantry angle so that the gantry arc rotation is simulated as a series of static gantry positions. Subsequently, a number of these MC simulations are accumulated as a VMAT field which represents an approximation for the simulation of dynamic motion of a linac. Taking into account the increased maximum speed of the Agility, we considered this approximation would be less appropriate for this new type of MLC. Our MC simulation therefore also used a linear interpolation method for the DOSXYZnrc code, so that the dynamic movement of the gantry during VMAT could be simulated without discretization. Compared with the results of film measurements, the gamma analysis of the results of the VMAT calculation for head and neck cancer showed a 99.5% pass rate with 3%/3 mm criteria.
Belec et al. reported that the modified BEAMnrc and DOSXYZnrc code could calculate VMAT 4D dose distributions for lung stereotactic treatment using an Elekta linear accelerator with a previous MLC model (33). An extended PSD file of recorded random numbers used for VMAT simulation in this study could also be applied to 4D dose calculations by synchronizing linac and patient motion.
The Agility MLC with its fast, complex movement and small leaf width could be used to provide superior dose distributions also to moving tumors. Further study is, however, needed to expand the current study to include evaluation of the risk of such complex movements possibly introducing large dose deviations due to interplay effects. To allow such studies by MC simulations, accurate modeling of the applied MLC is essential. Our study may prove useful for this purpose since the BEAMnrc and DOSXYZnrc codes are freely available open source and flexible, and we could arrange its original cord as we want.
Finally, we note that the MLC parameters in this study were specific to the MC model of our Agility MLC and Synergy accelerator. We were, however, able to tune five similar accelerators at our department to nearly identical parameters, typically within a 2% overall dose difference, and therefore all these accelerators were judged to be fully interchangeable for treatment plans. Accelerators commissioned elsewhere might need slightly different beam model parameters (29, 34), but with the adaptations our conclusions should then be also valid.
Conclusion
The Agility MLC mounted on the Synergy linear accelerator for 6 MV could be successfully MC modeled based on an existing component module and the high accuracy of its calculations could be validated against ionization chamber and radiochromic film measurements. A BEAMnrc and DOSXYZnrc based Monte Carlo model could thus be developed as an accurate tool for patient-specific QA and commissioning of treatment planning systems, even for the most complex current treatment techniques such as IMRT and VMAT.
Acknowledgements
The Authors wish to thank Korevaar EW, Wauben DJL and van't Veld AA (University Medical Center Groningen) for assistance with this project. They also thank Elekta AB for providing the specifications needed for the linac simulations. This study was partially supported by JSPS KAKENHI Grant (Grant-in-Aid for Young Scientists 19K17285).
Footnotes
Authors' Contributions
Research design: Ohira S, Takegawa H; experimental work: Ohira S, Takegawa H; statistical analysis: Ohira S, Takegawa H; manuscript writing: Ohira S, Takegawa H, Miyazaki M, Koizumi M, Teshima T; review & revision: Ohira S, Takegawa H, Miyazaki M, Koizumi M, Teshima T.
This article is freely accessible online.
Conflicts of Interest
The Authors have no conflicts of interest to declare regarding this study.
- Received May 27, 2020.
- Revision received June 18, 2020.
- Accepted June 22, 2020.
- Copyright© 2020, International Institute of Anticancer Research (Dr. George J. Delinasios), All rights reserved