Validation and clinical implementation of an accurate Monte Carlo code for pencil beam scanning proton therapy

Abstract Monte Carlo (MC)‐based dose calculations are generally superior to analytical dose calculations (ADC) in modeling the dose distribution for proton pencil beam scanning (PBS) treatments. The purpose of this paper is to present a methodology for commissioning and validating an accurate MC code for PBS utilizing a parameterized source model, including an implementation of a range shifter, that can independently check the ADC in commercial treatment planning system (TPS) and fast Monte Carlo dose calculation in opensource platform (MCsquare). The source model parameters (including beam size, angular divergence and energy spread) and protons per MU were extracted and tuned at the nozzle exit by comparing Tool for Particle Simulation (TOPAS) simulations with a series of commissioning measurements using scintillation screen/CCD camera detector and ionization chambers. The range shifter was simulated as an independent object with geometric and material information. The MC calculation platform was validated through comprehensive measurements of single spots, field size factors (FSF) and three‐dimensional dose distributions of spread‐out Bragg peaks (SOBPs), both without and with the range shifter. Differences in field size factors and absolute output at various depths of SOBPs between measurement and simulation were within 2.2%, with and without a range shifter, indicating an accurate source model. TOPAS was also validated against anthropomorphic lung phantom measurements. Comparison of dose distributions and DVHs for representative liver and lung cases between independent MC and analytical dose calculations from a commercial TPS further highlights the limitations of the ADC in situations of highly heterogeneous geometries. The fast MC platform has been implemented within our clinical practice to provide additional independent dose validation/QA of the commercial ADC for patient plans. Using the independent MC, we can more efficiently commission ADC by reducing the amount of measured data required for low dose “halo” modeling, especially when a range shifter is employed.


| INTRODUCTION
The use of Pencil Beam Scanning (PBS) is expanding rapidly in proton therapy, in large part because the approach produces highly conformal dose distributions and facilitates optimized delivery, without the requirement of field-specific hardware such as compensators or apertures, in contrast to conventional double scattering and uniform scanning delivery. At the University of Pennsylvania Roberts Proton Therapy Center, PBS delivery has been implemented for clinical treatment on two universal nozzles and one dedicated nozzle.
MC-based dose calculation is generally superior to analytical algorithms commonly used in treatment planning system (TPS) in modeling the dose distribution for PBS treatments. [1][2][3] This is particularly true when protons propagate through bone-soft tissue, soft tissue-air, and bone-air interfaces in treatment sites such as head and neck and lung, as multiple Coulomb scattering (MCS) can lead to a distortion of the field and inadequate target coverage. 2,3 While PBS eliminates most patient specific hardware, a beam modifier is still required in some situations. Various technical constraints in the current delivery systems result in a minimum proton energy limitation of between 70 and 100 MeV, 4 thus a range shifter is needed to degrade the proton range in order to treat tumors located shallower than the minimum range. 5 It is well known that the energy spread (due to energy straggling) and spot size (due to MCS) increase at the exit of a range shifter. 6 To minimize the spot broadening, the air gap between the range shifter and patient should be as small as possible, though the potential for collision with the patient often requires a gap that is larger than physically optimal. Due to the generation of secondary products as well as the particle transport within the air gap, it is difficult to model the dose calculation with a range shifter analytically given a limited measured data set, 5 thus an approach such as MC is desirable; the broader MC generated dataset is valuable for analytically approximating the low-dose halo 5 using multi-Gaussian lookup tables in water or in air after a range shifter given the magnification of potential MCS and halo calculation inaccuracies by various range shifter thicknesses and air gaps. 7,8 For any MC dose calculation, the first step is always to construct an accurate source model to parameterize the proton's distribution information in phase space (beam size, angular divergence and energy spread) at the position where it enters the simulated area. Several papers [9][10][11] have reported how to develop such beam source model by deriving source parameters through a set of simple measurements for individual beam lines. The major advantage is that this does not require knowledge of beam line or nozzle components and material compositions, and hence significantly reduces computing time without the need to model the nozzle. As the halo caused by the range shifter is intrinsically different from a halo in vacuum, a proper characterization of the halo component of the beam, below a factor of 10 −4 of the central axis dose is necessary, 12 as is a double-Gaussian fluence model, 8 to avoid dose inaccuracies. Hence, the range shifter in this paper is included as part of the simulated area as described by Grevillot 10 rather than creating an additional source model. The methodology is subsequently validated using a comprehensive set of measurements in water, both without and with range shifter to emphasize role of the low-dose halo, and also an anthropomorphic lung phantom for dose accuracy in heterogeneous medium.
The aim of this work is to develop and validate an accurate dose calculation platform based on TOPAS 13 that can be used to configure and validate both commercial 14,15 and in-house 16

2.A.1 | Modeling the beam optics
The model uses the IEC61217 gantry coordinate system, where the source plane is on the positive z-axis and the origin is at isocenter.
The source plane is set at the upstream surface of the range shifter to ensure that the range shifter can be correctly calculated when used in treatment (Fig. 1). A parameterization of the source model at the source plane is therefore required, including the spatial beam spread distribution (beam spot size), σ x , and the angular spread distribution (beam divergence), σ xθ , as well as the coefficient of correlation ρ x (the same relation holds for the y-direction). According to Courant-Snyder's particle transportation theory, 24 the σ-matrix of a beam's parameters at any location Z along the beam path, neglecting dissipation and diffusion processes, can be described as from which we infer that the variance of the spot size along the beam path should satisfy Spot profiles at six locations in air along the Z-axis (455, 330, 200, 100, 0, and −100 mm) were acquired using a scintillation screen/CCD camera detector (Lynx ® -IBA Dosimetry, Schwarzenbruck, Germany) with a 0.5-mm resolution 25 which is positive for a defocusing beam and negative for a focusing beam. Although σ θ increases slightly with propagation in air due to MCS, we approximate it as a constant in air between the nozzle exit and the phantom surface. The beam optic parameters above are derived to reproduce the measured spot variance in air which has taken into account the slightly increased divergence due to the scattering effect of air. The space between source plane and the simulated object, therefore, is set to vacuum in the MC simulation.
In contrast to the parallel scanning PBS system at PSI, 4

2.A.3 | Modeling protons per MU
The reference dosimetry approach proposed by Gomà et al. 29 where D Meas is the dose measured by ionization chamber. Fracchiolla et al. 26 have reported that the difference in protons per MU between this approach and that using a Faraday cup is 0.5% on average.  number of measurements for the Bragg peak curves and spot profiles. The interaction of the protons with the range shifter generates additional secondary particles resulting in a larger halo; propagation through the air gap between range shifter and patient will create significant difficulties for both experiment and simulation. 5 To address these challenges, we simulate the range shifter as an object within the beam path as described by Grevillot, 10 Fig. 2(a). We can observe that the σ in the x direction first focuses (decreases) then subsequently defocuses (increases) from upstream to downstream, while continuously defocusing in y direction. This is due to the integrated focusing effect of the two quadruples as well as less air scattering in the dedicated nozzle compared to universal nozzle. The spot sigma generally decreases with energy, and the shape is more elliptical for lower energy [ Fig. 2(b)]. From 210 to 225 MeV, however, the spot sigma unexpectedly increases. We speculate that this phenomenon is due to better beam focusing at 210 MeV than at 225 MeV. Figure 2(c) shows the dependence of σ θ on energy, which decreases from 6 mrad at 100 MeV to~3 mrad at 210 MeV. This is comparable with the values reported by Grevillot et al. 10

3.D | Validation of dose distribution in water
The measured and simulated depth doses along the central axis for three different SOBPs, both without and with a range shifter, are shown in Fig. 9. Simulations agree well with measured data, with a maximum dose difference of less than 2.2% and a clinical range agreement within 0.6 mm. Tables 1 and 2 Table 3. It can be observed that the gamma pass rate is improved significantly, from 66% to over 93% for TOPAS over the axial plane, while the sagittal and coronal plane agreements were improved from below 85%, the passing threshold, to over 98%. The output measurement results showed an overestimation of dose to the center of the target by 4% for TPS while TOPAS had a good agreement within 1% of measurement. Although TOPAS has better general agreement with measurement than the TPS, we can find TOPAS overestimates the dose in F I G . 8. Percentage differences between calculated and measured FSF for three monoenergetic fields at two depths as a function of proton energy without (a) and with (b) a range shifter. The red markers represent the results at the surface and the blue markers represent depths close to the Bragg peak.
the plateau region in Figs. 10(c) and 10(d), though the dose difference is still within 5%. Large-dose differences can be observed in the distal fall-off region of the field along the left-right direction in Fig. 10(c), which can be ascribed to a range difference. As this region is still within the lung, with a relative stopping power ratio of~0.31, the 10 mm geometrical range difference is roughly equivalent to a 3.1-mm water equivalent thickness (WET) difference, that is, 1.9% of the nominal range (166 mm) of the lateral left-right field. Therefore, F I G . 9. Central axis depth doses calculated by TOPAS (dashed line) are compared with measurements (marker) for three proton beams of varying ranges and modulation at a field size of 96 mm without (a) and with (b) a range shifter. Depth dose of R120M40 was renormalized by multiplying 105% to avoid overlap with R200M100.
T A B L E 1 Comparison of dosimetric parameters of lateral dose profiles at mid-range of SOBPs without the range shifter, measured using the Lynx in a Solid Water ® phantom.   Fig. 11(a) with the iCTV and liver (total liver minus GTV) DVHs in Fig. 11(b). The coverage (D95, the maximum dose that covers 95% of the target volume) and the overdose (D02, the maximum dose that covers 2% of the target volume) indices are within 1%    Fig. 12   scale. Figure 13 also shows the fractional integrated dose as a function of radius, normalized to an integration radius of 200 mm. One typically considers the dose contribution beyond three standard deviations of the spot size as that attributed as indirect dose contribution from halo (due to large-angle protons and secondary particles resulting from nuclear interactions). 10

CONF LICT OF I NTEREST
There is no conflict of interest to declare.