Implementation of an improved dose‐per‐MU model for double‐scattered proton beams to address interbeamline modulation width variability

Because treatment planning systems (TPSs) generally do not provide monitor units (MUs) for double‐scattered proton plans, models to predict MUs as a function of the range and the nominal modulation width requested of the beam delivery system, such as the one developed by the MGH group, have been proposed. For a given nominal modulation width, however, the measured modulation width depends on the accuracy of the vendor's calibration process and may differ from this nominal value, and also from one beamline to the next. Although such a difference can be replicated in our TPS, the output dependence on range and modulation width for each beam option or suboption has to be modeled separately for each beamline in order to achieve maximal 3% inaccuracy. As a consequence, the MGH output model may not be directly transferable. This work, therefore, serves to extend the model to more general clinic situations. In this paper, a parameterized linear‐quadratic transformation is introduced to convert the nominal modulation width to the measured modulation width for each beam option or suboption on a per‐beamline basis. Fit parameters are derived for each beamline from measurements of 60 reference beams spanning the minimum and maximum ranges, and modulation widths from 2 cm to full range per option or suboption. Using the modeled modulation width, we extract the MGH parameters for the output dependence on range and modulation width. Our method has been tested with 1784 patient‐specific fields delivered across three different beamlines at our facility. For these fields, all measured outputs fall within 3%, and 64.4% fall within 1%, of our model. Using a parameterized linear‐quadratic modulation width, MU calculation models can be established on a per‐beamline basis for each double scattering beam option or suboption. PACS number: 87.53.Qc

Because treatment planning systems (TPSs) generally do not provide monitor units (MUs) for double-scattered proton plans, models to predict MUs as a function of the range and the nominal modulation width requested of the beam delivery system, such as the one developed by the MGH group, have been proposed. For a given nominal modulation width, however, the measured modulation width depends on the accuracy of the vendor's calibration process and may differ from this nominal value, and also from one beamline to the next. Although such a difference can be replicated in our TPS, the output dependence on range and modulation width for each beam option or suboption has to be modeled separately for each beamline in order to achieve maximal 3% inaccuracy. As a consequence, the MGH output model may not be directly transferable. This work, therefore, serves to extend the model to more general clinic situations. In this paper, a parameterized linear-quadratic transformation is introduced to convert the nominal modulation width to the measured modulation width for each beam option or suboption on a per-beamline basis. Fit parameters are derived for each beamline from measurements of 60 reference beams spanning the minimum and maximum ranges, and modulation widths from 2 cm to full range per option or suboption. Using the modeled modulation width, we extract the MGH parameters for the output dependence on range and modulation width. Our method has been tested with 1784 patient-specific fields delivered across three different beamlines at our facility. For these fields, all measured outputs fall within 3%, and 64.4% fall within 1%, of our model. Using a parameterized linear-quadratic modulation width, MU calculation models can be established on a per-beamline basis for each double scattering beam option or suboption.
Output calculation methods for different proton therapy systems, using either empirical models (1)(2)(3)(4)(5)(6)(7) or Monte Carlo simulations, (8)(9)(10) have been described in several publications. These methods have been used for MU determination and/or as independent checks of measured output. Although less sophisticated, the advantage of empirical models over Monte Carlo simulations is their explicit form.
An analytical expression for the depth-dose of a spread-out Bragg peak (SOBP) at infinite source-to-axis distance (SAD) was derived by Bortfeld and Schlegel in 1996. (11) Kooy et al. (1) extended Bortfeld and Schlegel's analysis to a specific model for the IBA (Ion Beam Applications, SA, Louvain-La-Neuve, Belgium) double scattering proton system at MGH (Massachusetts General Hospital, Boston, MA, USA). In that model, a relationship between the output and a single factor r = (R-M)/M, which is a function of the distal range, R, and modulation width, M, of the SOBP, was established. Kooy et al. (2) improved their model by adding a correction factor that takes into account the shift of the effective source position as a function of proton distal range due to the change of fixed scattering materials. Engelsman et al. (12) further refined the model by redefining the modulation width to be the distance between the proximal 98% dose level and the distal 90% dose level, rather than between the proximal and distal 90% dose levels, as the position of the proximal 98% point is well defined and has less uncertainty than the position of the proximal 90% point. In the current implementation, the MGH group can predict outputs to within 1.4% (one SD) of measurements. (12)(13)(14) As the University of Pennsylvania (UPenn) uses an IBA proton therapy system that is similar to the one at MGH, there was interest in commissioning the MGH model in our clinic. The main difference between the MGH and the UPenn systems as it pertains to output stems from the MGH group's freedom to adjust the beam current modulation (BCM) of their system. While this enables the MGH group both to fine-tune the flatness of their SOBP distributions and to bring measured modulation widths into line with nominal modulation widths, (13)(14) this is not something that is permissible contractually on the UPenn system, nor on IBA systems installed elsewhere. Due to this, we found that desirable output prediction accuracy could not be achieved by applying the MGH model directly to the UPenn proton system. To apply the methodology to our center (and, by extension, to others), it is necessary to introduce a linearquadratic transformation from the nominal modulation width to the measured modulation width. In this paper, we describe a method to overcome the problem that arises when implementing the MGH-type semiempirical MU calculation procedure if these two widths differ appreciably. We first present how the model parameters are determined from limited measurements of systematic outputs, and then compare outputs predicted by this extended model with patient-specific field measurements.

II. MAtErIALS And MEtHodS
The IBA double scattering technique utilizes eight treatment "options", designated B1-B8, each of which is defined by a unique combination of second scatterer, range modulator wheel track, and BCM, and is applicable to a limited span of beam ranges. In order to be able to use a range modulator wheel track over a wide proton energy range, IBA applies different BCM for three subspans of range within an option. Such a subspan within an option is called a "suboption" and is designated by the suffix _1, _2 or _3 (e.g., B1_1, etc.).
The MGH MU model is derived fundamentally from SOBPs of ideal modulation width, where nominal and measured values are identical, and is not, therefore, directly applicable to practical situations where these values differ. In each IBA beamline, a polynomial fit is performed by the vendor prior to customer acceptance that relates the window of time during which the beam is on as the modulator wheel rotates to the resulting modulation width. Where these fits perform less well, differences in the nominal (i.e., fitted) and measured widths can exceed 1 cm in length or 3% in dose at the nominal proximal 90% dose point in some extreme cases, although IBA can generate flat SOBPs to make the dose at this point fall within 88% to 92% dose for typical modulation widths. For instance, Table 1 shows, for an example beam range of 17.5 cm, that measured and nominal modulation widths agree to within 2 mm for 5 and 10 cm nominal modulation widths across the three proton double scattering beamlines at our facility (named P1, P4, and P5), but that there are marked differences for 2 cm and full modulation widths. For the shortest modulation width, this translates into ~ 18% interbeamline variation in output (Table 1). By comparison, the consistency of measured range and output within the same day is better than 0.5 mm and 0.5%, respectively. Moreover, just as the Eclipse treatment planning system (Varian Medical Systems, Inc., Palo Alto, CA) can be configured to account for the difference between nominal and measured modulation widths, (15) we seek to relate these two widths for the purpose of output prediction. We propose to do so through a linear-quadratic model: where r nominal is related to a beam's nominal range, R, and nominal modulation width, M, via The constant 0.91 in Eq. (2) is a theoretical value used for converting our definition of modulation width (proximal 90% to distal 90%) to the original definition of proximal 100% to distal 100% by Bortfeld and Schlegel, (11) and was derived according to Eq. 8 in their paper. By propagating the MGH model, the output at the center of the SOBP is then given as (3) where CF is a constant to correct for the output change per option, s is a fit parameter to account for the variation of effective SAD within a beam option, and R m is the minimal range of the option. Equations (2) and (3) follow from the work performed at MGH; (1)(2) Eq. (1) is newly formulated here. Coefficients b 0 , b 1 , and b 2 are to be determined for each beam option using reference beams of the midrange suboption. Maximal R 2 or minimal residual sum of squares of the difference between the fitted and measurement data are used to determine the optimal parameters. After b 0 , b 1 , and b 2 are determined, a 1 and a 2 parameters are derived using the midrange suboption reference beams. However, when the midrange suboption's output dependence on modulation width does not represent that of the other two suboptions within an option, model parameters must be derived per suboption in order to fit the measured output data to within 2%. After b 0 , b 1 , b 2 , a 1 , and a 2 are determined, parameters s and CF are subsequently determined for the source position change with beam range and overall output constant using all the reference beams in all the three suboptions of each beam option.  Table 2 lists 60 reference beams that were used to derive the model coefficients in Eqs. (1) and (3). For each option, reference ranges were chosen at the two extremes and approximately midway between. For the minimal and maximal ranges, reference modulation widths (10 cm for B5-B8, 5 cm for B2-B4, and 3 cm for B1) were selected, while for the midrange suboption several modulation widths from 2 cm to full modulation width were utilized. Of these 60 beams, 47 were used initially to fit for the model parameters. We observed that dedicated fitting of the B2 and B6 suboptions was necessary to achieve 2% output accuracy. Therefore, 13 additional beams from the B2 and B6 low-and high-range suboptions are included for the fitting. The performance of the output model was initially validated with 28 reference beams from low-and high-range suboptions of B1, B3, B4, B5, B7, and B8 options ( Table 3) and subsequently tested with 1784 patient-specific fields.
Output measurements were made in a water phantom with SAD geometry using a PPC05 ionization chamber (IBA Dosimetry, Schwarzenbruck, Germany) aligned to isocenter and the center of the SOBP. Previous reports have investigated the output dependence on field size and snout position (5,7,16) and therefore the dependence of output on field size or snout position is not reported in this paper. Instead, a 10 × 10 cm 2 field and an air gap of 15 cm were used in this work to represent the average scatter condition, which minimizes the discrepancies introduced by the difference of field size and snout position in patient fields from the reference conditions. We restrict our application of the model to patient fields above 5 × 5 cm 2 , as the measured outputs of smaller field size would often fall below 2% of the modeled outputs and need patient-specific measurement. As a nominal rate of 2 Gy per minute is always used for our double scattering delivery, dose rate dependence was not investigated. Table 2. 60 reference beams used to derive the beamline-specific and option-specific MGH model parameters. SOBP RxMy has range x cm and modulation y cm. A span of SOBP is called beam option B#. A subspan within an option is called a "suboption" and is designated by the suffix _1, _2 or _3 (e.g., B1_1, etc.

III. rESULtS
Tables 4, 5, and 6 displays the model coefficients extracted from Eqs. (1) and (3) for each option or suboption of the three beamlines. It was found that the coefficients of Eq. (1) had to be derived separately among the suboptions of B2 and B6 in each case in order to achieve output accuracy within 2%, but that suboption-specific parameterizations were not required for the other options. Depending on the sign of b 2 , r model will depart upward or downward from the linear relationship with r nominal (Fig. 1) and this upward or downward departure could be different for large and small modulation widths (small and large r, respectively). If the 2% accuracy of the fit of Eq. (3) could not be achieved for all modulation widths within an option or a suboption, r nominal was further broken into large and small modulation width components (options B5-B8 (Tables 4, 5, 6)). As an example, Fig. 1 shows the linear-quadratic relationship between r nominal and r model for each suboption of B6 in beamlines P1, P4, or P5. r model is similar to r nominal when r nominal is smaller than 2 (i.e., large modulation width cases (M / R > ~ 37%)). However, these terms diverge from one another for the different beamlines and different suboptions when r nominal is larger than 2 (i.e., in the case of small modulation widths (M / R < ~ 37%)). Further, it can be observed that measured and nominal modulation widths were well matched by the vendor for all three suboptions in beamline P1, but that there was imperfect optimization for beamlines P4 and P5, and also variation by suboption. Without the linear-quadratic transformation of nominal to measured modulation width (Eq. (1)), R 2 between the linear-fit r model and the nominal r nominal was above 99% for P1, but between 89% and 91% for P4 and P5. Using the linear-quadratic fit, all the R 2 are above 99%. Figure 2(a) shows the fit of output to nominal modulation width for the reference beams of each of options B1, B3, B4, B5, B7, and B8 in all beamlines; Fig. 2(b) shows the fit for the three suboptions of both B2 and B6. From these, the necessity for per-beamline modeling of the output can be seen. For instance, for option B5, P1 has significantly smaller outputs than P4 and P5 for small modulation widths, while for all suboptions of B6 the converse is true. This is because P1 measured modulation widths are longer than nominal values for option B5 beams when the modulation width is smaller than ~ 3 cm (hence higher MUs are required for  the same mid-SOBP dose), and because P4 and P5 measured modulation widths are longer than nominal values for option B6 beams when the modulation width is smaller than ~ 4 cm. Other suboptions also show interbeamline variations in output over some or all of the span of modulation widths. The maximal 7% (B5) and 18% (B6) interbeamline output variation at small modulation can be modeled individually to within 2% of measurement for each beam option (B5) or each suboption (B6) by using linear-quadratic correction. The measured modulations are close to the nominal values for P4 and P5 of B5 option and P1 of B6 option. Without a linear-quadratic correction of the nominal modulation width to the measured modulation width, the difference of modeled and measured B5 outputs in P1 and B6 outputs in P4 and P5 at small modulation can be modeled below the 7% and 18% interbeamline output variation, but at the expense of disagreement at medium and large modulation.  Figure 3 shows the difference between the model output prediction and measurement for 1784 patient-specific fields for beamlines P1, P4, and P5. The modeled output is within 2% of the measurements for more than 95% of these fields, and for only two of these fields does it exceed 3% (one with 3.05% in P5 and the other with -3.07% in P1). The distribution of the fields amongst options is shown in Table 7. Because P1 treats primarily brain tumors and pediatric patients, B4 and B5 are dominant options, whereas in P5 where most treatments are for prostate cancer, B8 is dominant.

IV. dIScUSSIon
We described a procedure to implement the MGH model for calculating the output of proton double-scattered beams as a function of range and measured modulation width. Since the nominal and measured modulation widths on our systems are different due to limitations in the vendor's ability to establish this correspondence over the full span of modulation widths for all beam ranges, we introduced quadratic parameters b 0 , b 1 , and b 2 to take into account this difference. All the model parameters can be derived from the 60 reference beams listed in Table 2. Subsequent measurements of 1784 patient-specific field outputs demonstrate agreement to within 3% of the model prediction with < 1.1% standard deviation. Such a good agreement for a large cohort of patient fields is at least comparable to the best results reported by the MGH group (1.4% standard deviation), and highly suggestive that the MGH output model can be implemented at institutions that may not have full control of how the SOBP is achieved. The additional parameters b 0 , b 1 , and b 2 model imperfections in the SOBP, which were unexpected in the MGH model, as parameters a 1 , a 2 , s, and CF are fundamentally related to the original formulation with empirical corrections. (11)

V. concLUSIonS
A linear-quadratic transformation of the nominal to the measured modulation width is essential to the clinical implementation of the MGH MU calculation model in order to account for imperfectly matched SOBP widths and achieve 3% output prediction accuracy. A method to derive the linear-quadratic coefficients b 0 , b 1 , and b 2 is established.