Comparability of three output prediction models for a compact passively double‐scattered proton therapy system

Abstract The purpose of this study was to investigate comparability of three output prediction models for a compact double‐scattered proton therapy system. Two published output prediction models are commissioned for our Mevion S250 proton therapy system. Model A is a correction‐based model (Sahoo et al., Med Phys, 2008;35(11):5088–5097) and model B is an analytical model which employs a function of r = (R’‐M’)/M’ (Kooy et al., Phys Med Biol, 2005;50:5487–5456) where R’ is defined as depth of distal 100% dose with straggling and M’ is the width between distal 100% dose and proximal 100% dose with straggling instead of the theoretical definition due to more accurate output prediction. The r is converted to ((R‐0.31)‐0.81 × M)/(0.81 × M) with the vendor definition of R (distal 90% dose) and M (distal 90% dose‐to‐proximal 95% dose), where R’ = R‐0.31 (g cm−2) and M’ = 0.81 × M (g cm−2). In addition, a quartic polynomial fit model (model C) mathematically converted from model B is studied. The outputs of 272 sets of R and M covering the 24 double scattering options are measured. Each model's predicted output is compared to the measured output. For the total dataset, the percent difference between predicted (P) and measured (M) outputs ((P‐M)/M × 100%) were within ±3% using the three different models. The average differences (±standard deviation) were −0.13 ± 0.94%, −0.13 ± 1.20%, and −0.22 ± 1.11% for models A, B, and C, respectively. The p‐values of the t‐test were 0.912 (model A vs. B), 0.061 (model A vs. C), and 0.136 (model B vs. C). For all the options, all three models have clinically acceptable predictions. The differences between models A, B, and C are statistically insignificant; however, model A generally has the potential to more accurately predict the output if a larger dataset for commissioning is used. It is concluded that the models can be comparably used for the compact proton therapy system.


| INTRODUCTION
Current treatment planning systems (TPS) do not provide accurate monitor units (MUs) used to deliver the prescribed dose for passively scattered proton therapy systems. 1 Therefore, output prediction models are needed to convert the prescribed dose to MU. This is mainly due to the variety of models and vendors of proton therapy systems currently available. Additionally, different systems have varying proton acceleration and extraction designs, beamline designs used to shape the beam, and nozzle delivery techniques. Combinations of these variations add uncertainty in output prediction models across different systems.
Various approaches of output (cGy/MU) models have been employed for proton therapy systems. A correction-based model such as the model proposed by Sahoo et al. 2 utilizes the multiplication of correction factors based on beam parameters (option, range, and modulation width) that can independently change the output in the model. In contrast, analytical models 3-6 use mathematically derived analytical formulas preferably combined with empirical fitting parameters that are fitted to the measured data. In addition, Hotta et al. 7 used a simplified Monte Carlo simulation to calculate MU for a passively scattered proton therapy system. In the model, the output is predicted by a product of the output measured in a standard condition and a clinical beam delivery condition factor (F calc,clinical ).
The F calc,clinical is a multiplication of three factors: a beam-spreading device factor (F BSD ), a patient-specific device factor (F PSD ), and a field size correction factor (F FS ). The F BSD is measured in the phantom, while F PSD and F FS are calculated using the simplified Monte Carlo simulation.
Recently, we commissioned the correction-based model 2 and one of the analytical models 4 for our compact double-scattered S250 proton therapy system (Mevion Medical Systems, Littleton, MA, USA). In addition, a new model 8 mathematically derived from the analytical model was commissioned. The purpose of this research is to investigate the comparability of the three output prediction models for the compact proton therapy system.

2.A | Proton therapy system
The S250 proton therapy system is a single-room proton therapy system consisting of a synchrocyclotron which accelerates proton beams at a nominal energy of 250 MeV on a rotating gantry, a double scattering field shaping system (FSS), and a six-degree-of-freedom (6-DOF) robotic patient couch system. The beam delivery system has 24 total options with clinical ranges (R) from 5.0 to 32.0 g cm À2 . The R is defined as depth of the 90% dose point on the distal falloff in a normalized percent depth dose (PDD) scan and the modulation width (M) is defined as the width between the distal 90% and proximal 95% dose points of the PDD curve.
These options are separated into three general field types. The first option group is denoted "large" (options 1 to 12) with R from 5.0 to 25.0 g cm À2 mainly delivered by a large applicator with a maximum field size of 25 cm diameter. This option group uses six modulation wheels where every two options share the same modulation wheel, five high-Z first scatterers, and one second scatterer. The ranges are shifted using two sets of final absorber plates after the second scatterer (coarse in steps of 0.5 g cm À2 and fine in steps of 0.1 g cm À2 ) whose step resolution is 0.1 g cm À2 with a maximum change in range within an option 2.4 g cm À2 . The M for the large option group ranges from 2.0 to 20.0 g cm À2 in steps of 0.1 g cm À2 .
The second option group is denoted "deep" (options 13-17) with R from 20.1 to 32.0 g cm À2 delivered by a small applicator of the maximum field size of 14 cm diameter. This group shares one modulation wheel and a unique combination of one first scatterer and one second scatterer. The ranges are shifted using a static absorber wheel, with varying thicknesses by means of a circular wedge downstream of the modulation wheel set in the beam line for each option, and the final absorber plates. The M for the deep group ranges from 2.0 to 10.0 g cm À2 in steps of 0.1 g cm À2 .
The third group is denoted "small" (options 18 to 24) with R from 5.0 to 20.0 g cm À2 delivered by the small applicator. The M for the small group ranges from 2.0 to 20.0 g cm À2 in steps of 0.1 g cm À2 .
Each small option has a unique combination of a modulation wheel (total 7 wheels), one first scatterer (0.0 g cm À2 thickness for our system), and one second scatterer. The ranges are shifted using a static absorber wheel upstream of the modulation wheels rotated to a particular thickness corresponding to a range. Note that the combination of modulation wheels and scatterers for options is machine specific.

2.B | Correction-based model
A correction-based model was proposed by Sahoo et al. 2 and modified to predict the output for our Mevion system as follows: where W A (originally defined as (d/MU) wnc ; W A is used instead for comparison with the other models) is the output (cGy/MU) in water without range compensator, W o is the absolute output for a reference field (this is added to the original model since the output for our reference field for option 20 with R = 15.0 g cm À2 and M = 10.0 g cm À2 is 1.06 cGy/MU), ROF is the relative output factor normalized to the reference field, SOBPF is the spread-out Bragg peak (SOBP) factor, RSF is the range shifter factor, SOBPOCF is the SOBP off-center factor, OCR is the off-center ratio, FSF is the field size factor, ISF is the inverse square factor, and GACF is the gantry angle correction factor that is unique to the Mevion S250 system.
The correction-based model employs a large dataset obtained by fixing the beam parameters at reference conditions then incrementally changing one or more parameters while measuring the output at each increment. The output is then converted into a factor by taking the ratio of the measured output to the output of the reference conditions. The factors are then tabulated and either 1-D or 2-D linear interpolation is performed for data points that are not directly measured.

2.B.1 | Relative output factor (ROF)
The ROF is an option-specific parameter that accounts for relative changes in output when any option other than the calibration option (option 20 with R = 15.0 g cm À2 , M = 10.0 g cm À2 ) is used. It is the ratio of the measured output of a selected set of R and M for an option to that of the calibration option. Changes in the output for options are largely due to different beamline configurations. In this study, the deepest R for each option is used with M of 10.0 g cm À2 for ranges larger than 10.0 g cm À2 or a modulation width slightly less than the range for ranges less than 10.0 g cm À2 (a total of 24 measurements; see Table 1).
2.B.2 | Spread-out Bragg peak factor (SOBPF) The SOBPF accounts for changes in output for modulation widths that differ from the ROF modulation width within an option. It is measured for each option by fixing the R (the largest R within the option) and varying the M from 2 g cm À2 to the maximum M for the option in steps of 2 g cm À2 . The SOBPF is then found by taking the ratio of the output for the desired M to that of the reference ROF modulation width within the option. Due to relatively large output change between 2 and 4 g cm À2 modulation widths a SOBPF for 3 g cm À2 is also obtained (a total of 182 measurements).

2.B.3 | Range shifter factor (RSF)
The RSF accounts for output changes within an option due to range  Table 2). The RSF is determined for an option by 2-D linear interpolation of R and M in Matlab (Mathworks, Natick, MA; version R2013b). For some options, the maximum M is greater than for the minimum R. For instance, the minimum R for option 12 is 5.0 g cm À2 ; however, the maximum M is 6.0 g cm À2 . If the modulation exceeded the range, it would extend outside of the patient (e.g., for R = 10.0 g cm À2 and M = 11.0 g cm À2 , a high dose region of 1.0 g cm À2 will be outside of the patient). For this reason an R < M is not an allowable parameter entry into the Mevion system. Because of this the output cannot be measured for this set, but T A B L E 1 Relative output factors (ROFs) for all options. Additionally, the range (largest range in the option), modulation, and applicator field sizes are shown that were used for measurements.

Option
Range Field size (cm 2 ) ROF

2.B.4 | SOBPOCF and ISF
The SOBPOCF takes into account the relative changes in output when the measurement point is located longitudinally (along the beam line) away from the center of the SOBP. The SOBFOCF for a fixed source to detector distance (SDD) is given: where PDD is the percent depth dose normalized at the center of the SOBP, SSD is the source to surface distance used for the PDD measurement, d p is the depth of MU calibration, and d c is the depth of center of the SOBP. This factor is the dose at the calibration point relative to center of the SOBP with a fixed SDD (detector at the isocenter for both measurements in the published model). The ISF takes into account an MU calibration point that is located away from the isocenter: where SAD is the source to axis (isocenter) distance and SPD is the source to point of MU calibration distance. We propose combining SOBPOCF and ISF in this study for several reasons. First, the PDD is cumbersome to measure for all different sets of R and M with a fixed SSD. Second, the SAD for our system is option-specific and thus different effective source to axis distance (ESAD) should be used for each option. Third, to minimize the setup error in our institution, the MU calibration is not performed with a fixed SDD with off-center distance, but always at the isocenter and center of the SOBP (i.e., SDD = ESAD and chamber depth = R-M/2) and corrected for off-center and off-isocenter changes. In theory, the outputs at p c (point at center of the SOBP) and p MU (off-center point) are the same in the isocentric setup (isocenter at the depth of R-M/2) when the p MU is on the flat SOBP [Ψ(at p c )=Ψ(at p MU )] as shown in Fig. 1.
The output at p' MU (MU calibration point in patient plan from the TPS) is inversely proportional to the distance squared from the source between p MU and p' MU . The output at p' MU is given by:  | 111 and center of the SOBP for each option with the same R and M for ROF (a total of 24 measurements; Table 1). The profile data are then tabulated to be used for 2-D linear interpolation once the distances (x, y) from the central axis are found in a beam's eye view from the TPS.

2.B.6 | Field size factor (FSF)
The field size factor accounts for relative scatter changes due to changing the collimated open field size. Our proton system does not have rectangular variable collimators, and thus the FSF entirely relies on the brass block cutout. It is found by taking the ratio of the output of a field size to the reference field size (10 9 10 cm 2 for small and deep option groups and 20 9 20 cm 2 for large option group).
The outputs for three options with the deepest R and intermediate

2.B.7 | Gantry angle correction factor (GACF)
The Mevion S250 system is unique in its design by the fact that the entire cyclotron rotates with the external gantry. This causes a variation in output depending on the gantry angle (total rotation of 190°f rom 355°to 185°with 0°facing down). 9 The gantry angle correction factor (GACF) is thus introduced to account for this variation by calculating the ratio of the output to the calibration beam angle (gantry angle = 0 o , facing downwards) to that of the clinical gantry angle: This is obtained using the Farmer-type chamber placed at the isocenter inserted in an acrylic compensator phantom provided by the vendor. Due to the limit of the phantom size (buildup to the chamber is 10 g cm À2 ), the GACF data points are taken for one or two options within each option group for every 15 degree (large:

2.B.8 | MU calculation
The output in this study (model A) is thus given: In order to eliminate uncertainties due to fluence perturbation by range compensator and patient scatters, the output is measured without range compensator in the water tank. As demonstrated by Sahoo et al., 2 in clinical practice a verification plan using a virtual water phantom without range compensator is generated from a patient plan with the same geometry (a calibration point located at the same water equivalent depth, that includes the water equivalent thickness of the range compensator along the ray between the source and the point). The MU to deliver the prescribed dose at the calibration point in patient plan is given: where D pp is dose at the calibration point in patient plan, Ψ p is the output at the point in patient, Ψ wnc = Ψ A , and D vpnc is dose in the verification plan without compensator at the calibration point.

2.C | Analytical models
The second model (model B) is the analytical model 2 as a function of r: where p 0 , p 1 , p 2 , p 3 , p 4 , s 2 , and s 3 are option-specific fitting parameters. The correction factors (OCR, FSF, and GACF) are added to the original analytical models to have the equivalent formula.
To fit the basic model (the first term on the right side of models B and C), the same dataset for SOBPF for model A (absolute outputs) is used. The output change by effective source shift due to range change within an option (the second term) is fit by three measurements for each option (maximum, medium, and minimum R with corresponding M for the same value of r). All the option-specific parameters for models B and C are fit by the Matlab curve fitting toolbox. Implementation of the analytical models in detail can be found in our previous publication. 8    For these same option groups, larger modulations (M > 10 g cm À2 )

3.A | Correction factors
show little changes in SOBPF. The deep option group shows very little deviation between each option. This is mainly due to sharing the same FSS devices (a first scatterer, a modulation wheel, and a second scatterer) for all options.  The field size factors (FSFs) of various field sizes for selected options are shown in Table 3. It can be seen that field sizes greater  Relatively higher readings may come from the slit scattering for the detector at shallow depths. The reference field sizes are marked in bold.
(gantry angle = 0°) while the output for small options varies up to 3.6% (maximum deviation at angle 150°) from calibrations conditions.

3.B | Validation measurements and model comparison
For the total 272 dataset, the percent differences fell within AE 3% for the three different models. Figure 6 shows histograms of the percent differences. The average differences (AE SD) were À0.13 AE 0.94%, À0.13 AE 1.20%, and À0.22 AE 1.11% for models A, B, and C, respec-

| DISCUSSION
The correction-based models, such as model A, tend to have better output predictions (smaller SD) than analytical models. This is due to the fact that the correction-based models use the measured data directly through linear interpolation. In contrast, the fitting parameters in the analytical models do not fit the measured data perfectly at all points. Therefore, the analytical models may have varying accu-

| 115
When the modulation is small, the SOBP is not entirely flat and thus usually the output measurement has higher statistical deviation. In this situation, more data collection will not significantly improve the accuracy of model fitting and extra care (phantom setup and detector position) should be taken when the output is measured.
Uncertainties can arise from a variety of reasons when commissioning the models. For instance, the placement of the chamber can cause slight output changes in the case of small modulation widths (usually M < 4 g cm À2 ) where the modulation may not be flat. For these cases, careful placement of the chamber should be taken in order to minimize the gradient effect 10 on the data that are used in the output models, which in turn will affect the output predictions.
As noted in Fig. 3 For small field sizes (≤5 9 5 cm 2 ), it is required to directly measure the output due to inaccurate calculation of scatter contributions by pencil beam algorithm used in the TPS and irregular shape of aperture for actual clinical beams. Additionally, it is highly discouraged to use an aperture with radius of <1 cm for deep seated targets since the depth dose curve of a proton beam significantly deviates from the distinct Bragg peak due to the multiple Coulomb-scattering effect, and thus the surface dose is higher than the target dose. 11 For large field sizes (> 5 9 5 cm 2 ), squarefield apertures (10 9 10 cm 2 for the small applicator and 20 9 20 cm 2 for the large applicator) can be used for clinical beams with uncertainty of up to 1% for most of the beams as shown in Table 3.
The Mevion S250 cyclotron is mounted on an external gantry located behind a wall adjacent to the treatment room. During beam preparation the external gantry rotates to match the positioning angle of the inner gantry (applicator). The gravitational sag due to the various cyclotron positions results in slight changes in the energy. These energy shifts are accompanied by SAD shifts and beam shape changes.
An energy modulation disk (EMD) inside the cyclotron system is used to correct for these changes as well as to steer the beam. It is believed that the subtle correction of the beam by the EMD, depending on the gantry angle, results in the variations of output. It is for these reasons that changing gantry angles can cause a maximum output change of 3.6%. Because one GACF (an average of several options within a group) was used for each option group, an error up to 0.5% was found.
It should be noted that these changes in output are system specific.
The outputs are measured in the water phantom without patient-specific range compensators. There is intrinsic uncertainty associated with conversion of absorbed dose in a water phantom without compensator for calibration measurement to absorbed dose in a patient. It can be accounted for using compensator and patient scatter factors (CSF and PSF; CPSF = CSF 9 PSF) as shown in Sahoo et al.'s publication. 2 For MU calculation, these are cancelled out by calculating dose in a virtual water phantom without range compensator in TPS (i.e., D pp = D vpnc ÁCSFÁPSF and Ψ p = Ψ wnc ÁCSFÁPSF in Eq. (7)). However, it has been reported that the pencil beam dose calculation algorithm had limitations in accurately modeling scatter in the range compensator and patient (2-3% differences between the measurements and the estimations). 12 Newhauser et al. has also reported that uncertainty of CPSF is approximately 1% for proton treatments of prostate cancer. 13 By considering all the sources of uncertainty above, the estimated uncertainty by the correction-based model for our system is about 2% for large field and 3.5% for small field (error propagation of all uncertainties; ROF = 0.5%, SOBPF = 0.5%, RSF = 1.0%, ISF OCF = 0.5%, OCR = 1.0% for central area, FSF = 1.0% for large field or 3.0% for small field (estimated for 2 9 2 cm 2 < field size < 5 9 5 cm 2 from a publication using photon beams; 14 a further study is required using proton beams), and GACF = 0.5%. The uncertainty of 2.0% for large field is equivalent to about 2 SDs (0.94% 9 2) of our measurements. For the analytical models, these are about 2.5% for large field (also equivalent to 2 SDs of our measurements; 1.20% 9 2) and 4.0% for small field (about 1.0% uncertainty of the basic model and source shift added to the uncertainty of ROFÁSOBPFÁRSF) as summarized in Table 4. In addition, uncertainty for CSF and PSF is considered clinically.

| CONCLUSION
For all options, all three models had acceptable predictions (<3% difference between prediction and measurement). The differences between model A, model B, and model C are statistically insignificant.
In general, the model A has a potential to more accurately predict output if a larger dataset for commissioning is used. Field sizes less than 5 9 5 cm 2 should have their output measured directly to ensure accuracy. System-specific changes in output, such as the OCR and GACF for our Mevion system, should be investigated for each individual system. It is concluded that the models can be comparably used for the compact passively scattered proton therapy system.

CONF LICT OF I NTEREST
The authors declare no conflict of interest.