TOPAS Simulation of the Mevion S250 compact proton therapy unit

Abstract As proton therapy becomes increasingly popular, so does the need for Monte Carlo simulation studies involving accurate beam line modeling of proton treatment units. In this study, the 24 beam configurations of the Mevion S250 proton therapy system installed recently at our institution were modeled using the TOolkit for PArticle Simulation (TOPAS) code. Pristine Bragg peak, spread out Bragg peak (SOBP), and lateral beam profile dose distributions were simulated and matched to the measurements taken during commissioning of the unit. Differences in the range for all Percent Depth Dose (PDD) curves between measured and simulated data agreed to within 0.1 cm. For SOBP scans, the SOBP widths all agreed to within 0.3 cm. With regards to lateral beam profile comparisons between the measured and simulated data, the penumbras differed by less than 1 mm and the flatness differed by less than 1% in nearly all cases. This study shows that Monte Carlo simulation studies involving the Mevion S250 proton therapy unit can be a viable tool in commissioning and verification of the proton treatment planning system.


| INTRODUCTION
Success in radiation therapy is dependent on maximization of the tumor control probability while minimizing the normal tissue complication probability. The characteristics of the proton Bragg peak help to deliver large doses to the target, lower entrance doses proximal to the target and almost no doses distal to the target. 1 Advancements in Monte Carlo simulation and radiation transport calculations have enabled more accurate characterizations of the radiation fields during proton treatments that can benefit patients. 2 The cases where the greatest advantages of proton therapy could be realized are in targets near critical organs and treatments involving pediatric patients. [3][4][5] In order for the full potential of proton therapy to be utilized, there is a need for an accurate dose calculation in proton treatment plans. Monte Carlo simulation has traditionally been shown as a prominent method for conducting various research topics that include dosimetric, linear energy transfer (LET), and commissioning studies. [6][7][8] In a study done by Paganetti et al., it was shown that the modeling of the IBA proton treatment head at the Northeast Proton Therapy Center resulted in simulated data matching with measured beam data within millimeter accuracy for beam range and 3 mm for SOBP width. 9 This acceptable tolerance helps generation of detailed simulated beam data for use in commissioning of a treatment planning system. 10 Monte Carlo simulation has also been used to calcu-from nuclear reactions in the treatment system on patients that have undergone proton craniospinal irradiation treatments. 11 Paganetti has also reported that Monte Carlo simulations can even improve proton beam range uncertainty by up to 2.2%. 12 A passive double scattering compact proton therapy unit has recently been installed at our facility. The Mevion S250 is the first proton therapy system of its kind, delivering a pulsed beam with a nominal energy of 250 MeV and utilizing an in-room superconducting synchrocyclotron mounted on a gantry that allows for 180°rotation. The compact nature of the unit has many cost and therapy treatment advantages, and eliminates the need for a complex beam transport system. 13 The main beam shaping components in the nozzle of the unit include a lead first scattering foil present at the cyclotron exit for initial beam spread. The beam then passes through a bimaterial staircase type range modulator wheel (RMW) consisting of one track made of graphite to modulate the range of each Bragg peak, and a second track made of lead to ensure uniform scattering power over all steps of the wheel. The last step of each wheel has a brass wedge to completely stop the beam. A bilayer contoured second scattering foil made from lead and Lexan is incorporated downstream of the RMW to further spread and flatten the beam. The last components the beam passes through are a graphite absorber for energy degradation and a post absorber for fine tuning the beam range, followed by two ion chambers to monitor beam output. For more information on the initial clinical experience with the system, the reader is referred to Zhao et al. 13 The Mevion S250 has 24 different beam configurations, divided into large, deep, and small groups. The large group utilizes a large beam nozzle with an uncollimated field size of 25 cm in diameter where the deep and small group share a small beam nozzle with the field size of 14 cm in diameter. The distinction between deep and small groups occurs in their range capabilities and modulation widths. The deep group has a depth range from 20.1-32 cm with a maximum modulation width of 10 cm, whereas the small group has a shallower depth range of 5-20 cm but has a maximum modulation width of 20 cm. For any given group, a special beam configuration results from the unique order and combination of different beam line components as summarized in Fig. 1. In total, the Mevion system offers 12 beam configurations in the large group, 5 for the deep group and 7 in the small group which are composed from 18 first scattering foils, 14 RMWs, and 3 secondary scattering foils. Table 1 shows the beam characteristics for all configurations in each group.
The objective of the this study has been to model the beam line components for all 24 beam configurations of the Mevion S250 system using Monte Carlo simulation in order to accurately characterize the radiation field exiting the treatment head, and benchmark the simulation results with the measured beam data obtained during commissioning of the machine. The diode detector provided the necessary spatial resolution needed for capturing the sharp penumbra seen in lateral beam profile scans.

| ME TH ODS
This specific diode was chosen in part due to its resilience against radiation damage. McAuley et al. showed that a loss of sensitivity as a function of accumulated dose was only 1% per 100 Gy, which is well under the doses delivered during our measurements. 14 TOPAS (TOolkit for PArticle Simulation) version 2.0 was utilized to model the beam delivery system in this study. 15 TOPAS 2.0 is an extension of the GEANT4 10.1.p02. toolkits that was developed specifically as a user friendly proton therapy tool. It has been validated experimentally as a viable choice when tasked with reproducing beam data from a passive scattering proton system. 16 The dimensions   and materials for each beam component were modeled in TOPAS based on the data given by the manufacturer. For all simulations, the TOPAS default modular physics list was used, which is composed of "TsEmStandardPhysics_option3_WVI", "HadronPhysicsQGSP_BIC_HP", "G4DecayPhysics", "G4IonBinaryCascadePhysics", "G4HadronElas-ticPhysicsHP", "G4StoppingPhysics", and "G4RadioactiveDecay Physics". The range cut for secondary particles was set to 0.005 cm, with an energy cut in water of 990 eV, 57.3 keV, 5 keV, and 56.6 keV for gamma rays, electrons, protons, and positrons respectively. A water phantom with dimensions of 40 9 40 9 40 cm 3 , placed downstream of the delivery nozzle, was used. The detector mesh used to score the dose was 2 9 2 9 0.1 cm. The face of the detector mesh was large enough to allow for fewer particles to be run achieving a smooth PDD, while the fine resolution in the z direction prevented distortion of high gradient areas.

2.B | Spread out bragg peak simulations
The The pristine Bragg peaks from each step of the wheel were scaled by applying weighting factors to create a SOBP that matched to the measured data (see Eq. 1).
where D d ð Þ is the dose as a function of depth, w i is the weighting factor applied to the ith peak, and p d ð Þ i is the dose distribution associated with the ith pristine Bragg peak.
The sum of square errors between measured and simulated data for each point at 0.1 cm depth increments were minimized to give the best agreement between the calculated SOBP and commissioning data (see Eq. 2).
where SSE is the sum of square error, D d ð Þ m is the dose distribution for the measured data and D d ð Þ s is the dose distribution for the simulated data.
An example of an SOBP and corresponding pristine Bragg peaks for configuration 13 are shown in Fig. 2. Three SOBPs (one configuration from each of the group categories) were constructed using this method (detailed information is shown in Table 2).

2.C | Lateral beam profile matching
In a passive scattered proton system, scattering foils are put in place to spread the beam to a clinically relevant treatment size. For Individual peaks created from the method described in the text to be summed to create the SOBP. Each peak was assigned a specific weighting factor that, when summed together with the other peaks, created a flat SOBP shown by the dotted green curve.  Table 3 gives a summary of each of the configurations and depths used to calculate the lateral beam profile.

2.D | Spread out bragg peak modulation width adjustments
In the clinical applications of a passive scatter system, it is necessary to adjust the modulation size of the SOBP for target coverage. The  Table 4). The average differences over the deep group for the PDD (0.5) and D90 depth was 1.3% (ranged from 0.6 to 1.8%) and 0.03 cm (ranged from 0 to 0.04 cm) respectively (shown in Table 5). The average differences over the small group for the PDD (0.5) and D90 depth was 1.2% (ranged from 0.3 to 1.9%) and 0.04 cm (ranged from 0.00 to 0.07 cm) respectively (shown in Table 6). An example from each group is shown in Fig. 3, where normalized depth dose curves from beam commissioning measurements and simulated data are plotted for comparison.
To achieve agreement between the measured data and the simulated data for the pristine Bragg peak of each configuration, small adjustments were made to the manufacturer provided geometry of the beam delivery system. As described by Bednarz et al., 18 there is uncertainty as to how accurate the provided blueprints reflect the real geometry of the components in the commissioned treatment head. Therefore, tuning some of the physical parameters of the geometry was necessary assuming that the adjustments were within manufacturer tolerance. 19 The range of proton beam was first matched by fine adjustments of the post absorber thicknesses using stopping power conversion ratios. In most cases, after range correction, the entrance dose and Bragg peak distal fall off needed to be tuned. Paganetti and Bednarz showed that the beam spot size, energy, and angular spread can influence the shape of the Bragg curve and provides one of the largest sources of uncertainty, due to the difficulty in measuring these parameters. 10,18 Paganetti et al., also showed that beam energy spread has the largest influence on entrance dose and distal fall off. 10 The same trend is also found in our study, where increasing the energy spread increased the entrance dose and the distal fall off length. Therefore, in order to tune our pristine Bragg peaks to appropriate distal falloff, the energy spread was adjusted iteratively until the best agreement between measured and simulated data was reached. We found that the energy spread that achieved the best agreement between the measured and the simulated data for most of the configurations was 0.4%. The angular spread and spot size were then optimized in order to achieve further agreement in data. A 0.005 radian standard deviation in angular spread in both the x and y directions gave the best representation of commissioned data, along with a beam spot size of 8 mm for the majority of configurations.
In the case of the large configuration group, some of the shallower range configurations were not able to be tuned to match within our set criteria with adjustments to the source alone. The Mevion system has a consistent distal falloff margin regardless of the energy or configuration. We found that as we decreased the energy of the beam and increased the field size, the distal fall off margins in simulation began to degrade. This is likely due to the presence of more material in the beam that introduced large uncertainty in geometry; and influenced the beam energy distribution in a way not accounted for in our simulations. In these specific cases, small adjustments (on the order of a tenth of millimeter) to the thickness of the first scattering foil were made in order to adjust the distal gradient to agree with measurements. These adjusted first scattering foils were used in the geometry for the subsequent Bragg peak simulations.

3.B | Spread out bragg peak matching
At our institution, SOBP width is defined as the distance between the proximal 95% dose depth and the distal 90% dose depth. To evaluate the agreement between simulated results and measured data, the SOBP width, beam range, and the depth of distal 20% dose were compared. Figure 4 shows the matching of normalized percent

3.C | Lateral beam profile matching
Three of the lateral beam profiles are shown in Fig. 5 with the simulated and measured data plotted on each graph. The penumbras were calculated as the distance between the 80% and 20% dose levels. Beam flatness and symmetry was calculated using Equations 3 and 4 respectively.
where D min and D max are the minimum and maximum doses within the middle 80% of the field size.
where DL and DR are the integral doses of the left and right side of the radiation field respectively.
The full width at half maximum values (FWHM) for both measurement and simulated data were calculated and compared with each other. The absolute differences in the penumbras between simulated and measured profiles at each depth for each configuration all agreed to well within a millimeter. Flatness and Symmetry values for all 12 profiles (six configurations at two depths each) are listed in Table 8.   Using the modulation adjustment algorithm described in the methods section, a reference SOBP with various modulation widths were simulated and shown together with a pulse modulation function in Fig. 6

Mod. Width (cm)
For a desired modulation width of 10 cm, the beam stops after 40.77 pulses. Linear interpolation is used to determine the fraction of the last stop pulse to be applied when the desired modulation width is between pulses. This curve was used to provide the stop pulse to simulate the reference beam at our institution. (b) The measured data from commissioning (circles) and simulated (solid lines) normalized spread out Bragg peak for the reference SOBP. system with acceptable accuracy and possible reduction of the intensive time required in measured data collection for the commissioning and verification of the proton treatment planning system.