Commissioning of GPU–Accelerated Monte Carlo Code Fred for Clinical Applications in Proton Therapy

We present commissioning and validation of Fred, a graphical processing unit (GPU)–accelerated Monte Carlo code, for two proton beam therapy facilities of different beam line design: CCB (Krakow, IBA) and EMORY (Atlanta, Varian). We followed clinical acceptance tests required to approve the certified treatment planning system for clinical use. We implemented an automated and efficient procedure to build a parameter library characterizing the clinical proton pencil beam. Beam energy, energy spread, lateral propagation model, and a dosimetric calibration factor were parametrized based on measurements performed during the facility start-up. The Fred beam model was validated against commissioning and supplementary measurements performed with and without range shifter. We obtained 1) submillimeter agreement of Bragg peak shapes in water and lateral beam profiles in air and slab phantoms, 2) $<2\text{\%}$ dose agreement for spread out Bragg peaks of different ranges, 3) average gamma index (2%/2 mm) passing rate of $>95\text{\%}$ for $>1000$ patient verification measurements using a two-dimensional array of ionization chambers, and 4) gamma index passing rate of $>99\text{\%}$ for three-dimensional dose distributions computed with Fred and measured with an array of ionization chambers behind an anthropomorphic phantom. The results of example treatment planning study on $>100$ patients demonstrated that Fred simulations in computed tomography enable an accurate prediction of dose distribution in patient and application of Fred as second patient quality assurance tool. Computation of a patient treatment in a CT using ${10}^4$ protons per pencil beam took on average 2′30 min with a tracking rate of 2.9$\times$${10}^5$${\text{p}}^{+}/\text{s}$. Fred was successfully commissioned and validated against the clinical beam model, showing that it could potentially be used in clinical routine. Thanks to high computational performance due to GPU acceleration and an automated beam model implementation method, the application of Fred is now possible for research or quality assurance purposes in most of the proton facilities.

1 Introduction

In proton radiation therapy, Monte Carlo (MC) methods offer more accurate modeling of proton interactions with heterogeneous media and improve dose calculation accuracy in complex geometries with respect to analytical pencil beam algorithms [1–4]. The application of MC algorithms in treatment planning can eventually lead to a reduction in the target volume safety margins by about 2% and more accurate prediction of the treatment outcomes [5]. The state-of-the-art commercial proton beam therapy (PBT) treatment planning systems (TPS) employ MC methods for treatment plan optimization and dose calculation [6, 7], but they are still not the standard treatment planning tools in all clinically operating PBT facilities. Many proton facilities still use analytical pencil beam algorithms of limited accuracy in heterogeneous media. Also, the time performance of the MC-based TPS remains to be an issue, especially when applying robust optimization algorithms that require computing several dose distributions for one computed tomography (CT) image or in treatments of moving targets where 4D-CT consisting of a series of CT images of several motion phases of one patient are employed in treatment plan optimization [8]. In addition, proton radiation therapy quality assurance (QA) procedures are time consuming and require manpower for experimental measurements of dose distributions in phantoms, typically performed at a few depths in water for each treatment field. In fact, time needed for patient QA could be dedicated for the actual patient treatment. Therefore, reduction in the number of measurements is widely discussed among medical physicists [9–14]. Supplementing or replacing patient QA measurements with dose distribution recalculation using a second, independent, dose-calculation engines can be beneficial for PBT facilities.

In several PBT facilities, general purpose MC simulation toolkits, such as: FLUKA [15], Geant4 [16, 17], or Shield-HIT [18] as well as more user-friendly environments built on Geant4 like GATE/GATE-RTion [19–21] and TOPAS [22, 23], are used to support research activities and/or simulations for patient QA. The clinical application of general purpose MC tools is limited, mainly due to the time required to recalculate a complete plan ranging from tens of minutes to even a few hours. For this reason, the parallelization of the particle tracking on several central processing units (CPU) or general purpose graphical processing units (GPU) is of interest for radiotherapy. The PBT-dedicated GPU-based MC code gPMC implemented by Jia et al. [24] was further developed [25] and validated using clinical patient data [26]. Following the gPMC development, Wan Chen Tseung and colleagues presented a high-performance GPU-accelerated MC code, which is used for routine clinical QA and as the dose calculation engine in a clinical MC-based Intensity Modulated Proton Therapy (IMPT) treatment planning system [27]. Recently, an analytical pencil beam algorithm, the FRoG platform, was implemented on GPU for clinical investigations with different ion types [28, 29].

The commissioning and validation of the independent, MC-based dose calculation engine for research or patient QA purposes is a time-consuming process that requires knowledgeable and experienced manpower. Only recently, standards for beam modeling and beam model commissioning for MC dose calculation–based radiation therapy treatment planning were proposed [30]. The experimental characterization of the proton beam properties (longitudinal and lateral profiles as well as dosimetric calibration) as a function of primary beam energy is facility dependent because different PBT centers use different accelerators, measurement methods, and TPS. The complete implementation of passive and active beam delivery nozzle geometry was described by Paganetti et al. [1] for cyclotron-based facilities and by Parodi et al. [31] for synchrotron-based facilities. However, it was suggested later that for MC dose calculation purposes, defining the beam model following the clinical commissioning procedure and avoiding detailed simulations of the beam nozzle geometry is possible with a precision that is sufficient for clinical application [10, 32, 33].

This article reports on commissioning of the GPU-accelerated MC code Fred [34] and its validation at two cyclotron-based proton beam therapy facilities of different beam line design: Varian ProBeam in Atlanta, GA (United States), and IBA Protheus C-235 in Krakow (Poland). The software toolkit Fred (Fast paRticle thErapy Dose evaluator) [34] was developed at the University of Rome for parallelized proton beam transport simulations in heterogeneous geometry defined by the patient CT. We describe in detail Fred commissioning steps, that is, automated characterization of the beam model that describes the proton beam used for patient treatment and follows the clinical QA procedures. Finally, we validated our commissioning procedure using the optimized beam models. We simulated dose distributions in Fred and compared the results with verification measurements performed in homogeneous and heterogeneous phantoms with and without range shifters as suggested by Winterhalter et al. [35]. Such extensive experimental validation of Fred accuracy and time performance has been never reported before. To increase the confidence of the reader about the accuracy of Fred simulations, selected results were also compared with clinical TPS simulations. Eventually, we evaluated clinical cases of patient treatment plans to demonstrate the clinical applicability of Fred.

2 Materials and Methods

2.1 GPU–Accelerated Monte Carlo Code Fred

The great benefit of Fred with respect to general purpose MC codes is its computation performance achievable on a variety of different hardware without compromising the dose computation accuracy. The typical tracking rates range from 10–100 thousand protons per second using a single CPU to about million particles per second using GPU cards. Fred is equipped with an interface to convert phantom/patient geometries stored in DICOM CT images to a voxelized geometry of the patient containing the atomic tissue composition using a conversion table based on stoichiometric calibration [36]. In addition to patient geometry, user-defined geometries of specific material composition can be included enabling simulations of proton transport in passive elements like range shifter.

The physical interaction models implemented in Fred are trimmed down with respect to general purpose MC codes, such as Geant4/FLUKA within the regime that is relevant for particle therapy, in order to speed up the execution time without compromising the accuracy of dose-deposition calculations. In particular, the physics processes contributing to the dose deposited by protons in patient tissue, that is, mean energy loss, energy fluctuations, nuclear elastic and inelastic interactions with target nuclei as well as the trajectory deflection via a multiple Coulomb scattering, are implemented in Fred [34]. Moreover, Fred offers linear energy transfer (LET) and relative biological effectiveness (RBE)–weighted dose calculations by means of different RBE models, providing further information, which is not available in the state-of-the-art commercial TPS. The LET and RBE computations in Fred are out of the scope of this study.

2.2 Commissioning Measurements and Fred Simulations

Fred commissioning was performed for one gantry room of two PBT facilities of different beam line design equipped with scanned proton beams that are in clinical operation since 2016 and 2018, respectively. Krakow facility is an IBA design based on Proteus C-235 cyclotron equipped with two rotational gantries, an eye treatment room and an experimental hall. The TPS Eclipse from Varian, version 13.6, is used for treatment planning in CCB. It uses analytical proton convolution superposition (PCS) algorithm for the dose calculation and optimization [37]. EMORY PBT center in Atlanta is a ProBeam system designed by Varian and equipped with three rotational gantries and two horizontal beam lines. The TPS RayStation from RaySearch laboratories, version 8A, equipped with MC dose algorithm is used for treatment planning in EMORY [7]. The properties of proton beams and the measurement methods used for the acquisition of clinical beam model commissioning data at both facilities are listed in Table 1.

TABLE 1

TABLE 1. Selected properties of CCB and EMORY PBT centers and measurement methods used for the proton beam model commissioning.

The commissioning measurements that include depth dose distribution measurements in water phantom, measurements of the lateral profiles (without range shifter) in air, and absolute dose measurements in a water phantom were used to build parameter libraries characterizing the Fred beam model for Krakow and Atlanta facilities. The water phantom and in-air setup used for commissioning measurements are schematically illustrated in Figures 1 A and B respectively. The figure indicates how the proton beam is transported from the nozzle toward the detector/phantom. During irradiation, the beam is deflected vertically and horizontally by scanning magnets and crosses a position sensitive ionization chamber (IC23), which is used for beam lateral position and size measurement. The procedure of the commissioning data acquisition is not described here in detail as it is out of the scope of this article.

FIGURE 1

FIGURE 1. Experimental and simulation setups for water phantom (A), in-air scintillating screen measurements (B), and measurements in a solid water slab phantom (C). On the left, beam nozzle elements (scanning magnets and position sensitive beam monitor (BM) chambers), not taken into account in MC simulations, are shown (gray scale). In MC simulations, the primaries are generated in Monte Carlo virtual source and transported through range shifter (RS) to phantoms/detectors (blue). The figure is not to scale.

The Fred simulation setup mimics the commissioning measurements setup shown in Figures 1 A and B. The virtual beam source is located at the position of the scanning magnet located closer to the isocenter because at this position, the deflection of the beam in both X and Y directions is defined. The different position of the X and Y scanning magnets is taken into account, while calculating the direction of a single pencil beam. The beam propagation in the IC23 is omitted in the simulations and is taken into account by adjusting beam source parameters, in such a way that the beam size fits the results of beam size measurements in air performed with scintillating screen (Lynx). The proton beam was propagated without and with range shifter. Fred simulations in water were performed in 40 $\times$ 40 $\times$ 40 ${\text{cm}}^3$ virtual phantoms of 1 $\times$ 1 $\times$ 1 ${\text{mm}}^3$ voxel size (Figure 1A). The ionization potential of water was set to 80 eV [38]. Fred simulations of the in-air setup used for beam model validation were performed in a virtual air phantom. The total time of Fred MC simulations includes tracking time, time needed for memory allocation, and the file writing. The tracking rate of simulation is given as the number of protons tracked per second (${\text{p}}^{+}/\text{s}$).

2.3 Beam Model Parameters

The beam model parameters characterize longitudinal and lateral pencil beam profiles as well as dosimetric calibration. Two parameters, energy (E) and energy spread (E_σ), characterize proton pencil beam depth dose distribution (longitudinal) profile. One further parameter, monitor units (MU) to the number of particles conversion factor (${\text{SF}}_{\text{MU}}$), characterizes integral dose distribution (IDD) dosimetrically, by means of dose measurement at 2 cm depth, following TPS commissioning protocol and other references [37, 39, 40]. The lateral propagation of the proton pencil beam can be characterized by a quadratic model by means of modeling beam emittance or bilinear model by defining virtual point source. In fact, the bilinear model is an approximation of a quadratic model in a limited range. The virtual point source approach can be applied when the waist of the quadratic function of emittance model is far enough from the isocenter to approximate lateral beam propagation behind the nozzle exit by a bilinear function. Fred is capable of handling lateral beam propagation using both virtual point source or emittance approaches.

For characterizing the lateral propagation, the lateral beam profiles measured during facility commissioning in air at different Z positions (cf. Figure 1B) were fitted using the Gaussian fit, and its σ(z) was calculated. Additionally, the σ(z) measured with the beam monitor chambers in the nozzle can be used [41]. This improves the quality of the lateral beam propagation modeling, especially in cases where the waist of the beam is located between the nozzle and the first measured point in air. Fitting σ(z) to commissioning data from both facilities at different distances from the isocenter using bilinear and quadratic functions indicated that the emittance model is appropriate for Krakow facility, whereas the virtual point source model can be used for EMORY.

For characterizing the beam lateral propagation in Krakow, six emittance model parameters ($\mathit{ϵ},\alpha ,\beta$), three in X direction and three in Y direction, were used. The Twiss parameters $\mathit{ϵ},\alpha ,\text{and}\beta$ were obtained according to the following formula [42]:

${\sigma }^2\left(z\right)=\mathit{ϵ}\cdot \left(\beta -2\cdot \alpha \cdot z+\frac{1+{\alpha }^2}{\beta }\cdot z^2\right),(1)$

where the emittance ϵ corresponds to the area in the X/Y position–velocity phase space and is assumed to be constant over the beam propagation in air. The Twiss parameter α is related to the focusing/defocusing of the beam, whereas β characterizes the length over which the beam changes its transverse shape.

For characterizing the beam lateral propagation in Atlanta, four parameters, two in X direction and two in Y direction, specific for a bilinear approximation were used. The parameters were obtained according to the following formula:

$\sigma \left(z\right)=S\cdot z-VSD,(2)$

where the S is the function slope and corresponds to the rate of the spot size variation and $VSD$ stands for virtual source distance and corresponds to the distance from the virtual source to the isocenter. Note that for both approaches, virtual point source and emittance model of lateral beam propagation, particles are transported starting from the position of the scanning magnets regardless of the position the emittance waist and $VSD$.

For TPS exploiting analytical pencil beam algorithm, the emittance model is defined for configurations with and without range shifter, whereas in MC-based TPS and in Fred, only the configuration without range shifter is defined, and proton transport in range shifter is simulated according to its model parameters (material composition, density, physical thickness).

2.4 Generation of Beam Model Parameter Library

We implemented a set of software tools that calculate beam model parameters in three automated steps (see Figure 2). The beam model parameter libraries were generated in the entire proton beam energy range in 10 MeV steps (Table 1) for both facilities. Figure 2 schematically illustrates how the Fred MC commissioning procedure uses the facility commissioning measurements as the input to obtain beam model parameters per nominal energy, that is, beam energy E, energy spread E_σ, MU scaling factor ${\text{SF}}_{\text{MU}}$, and six emittance or four virtual point source parameters. The procedure is automated and does not require any interaction with the user, except preparation of the measurement data. Fred simulations of single pencil beams were performed using ${10}^8$ primary protons.

FIGURE 2

FIGURE 2. A flow chart illustrating Fred commissioning and validation steps. Simulations steps benefiting from GPU-accelerated MC simulations are indicated (GPU).

Step 1. In the first step (Figure 2: Step 1), the emittance or virtual point source model (Eqs 1 and 2) was fitted to the measured beam spot size (${\sigma }_{x/y}$) as a function of the position along the beam (see Section 2.3). For Krakow beam model, in addition to the beam size measurements performed with Lynx (pixel size 0.5 $\times$ 0.5 ${\text{mm}}^2$), the beam size measurements performed during irradiation with IC23 (resolution 5 mm in $X/Y$ directions) installed close to the nozzle exit were used to fit the emittance model (see Section 2.3). In this way, emittance model parameters ($\mathit{ϵ},\alpha ,\beta$) or virtual point source parameters (S, $VSD$) were obtained for X and Y directions and each energy.

Step 2. In the second step (Figure 2: Step 2), beam energy (E) and energy spread (E_σ) were obtained. The measured and simulated IDD profiles were fitted using a formalism proposed by Bortfeld [43, 44]. Using the fit and semiempirical relations proposed by Bortfeld [43], the initial energy and energy spread of protons producing an IDD distribution were computed. The Bragg peak range (${\text{R}}_{80\text{\%}}$) defined as 80% of the maximal value at the distal falloff and the Bragg peak full width at half maximum (FWHM) were numerically calculated from the fitted curve. The E, E_σ, ${\text{R}}_{80\text{\%}}$, and FWHM parameters were calculated for experimental data and each Fred simulation. An automated iterative optimization procedure was developed to find such E and E_σ values in Fred, which minimize the absolute difference of Bragg peak range ($\left|{\text{ΔR}}_{80\text{\%}}\right|$) and FWHM ($\left|\text{ΔFWHM}\right|$) between simulation and measurement. The dependence of $\left|{\text{ΔR}}_{80\text{\%}}\right|$ and $\left|\text{ΔFWHM}\right|$ on E and E_σ is a continuous function with a single global minimum. The optimization procedure was implemented in Python exploiting the Nelder–Mead simplex algorithm [45]. The initial guess of energy and energy spread was estimated from the Bortfeld curve fitted to measured data. Each consecutive step of the optimization algorithm included the following: 1) new simulation of a depth dose distribution in water with energy and energy spread computed by the optimization algorithm, 2) Bortfeld curve fit and estimation of ${\text{R}}_{80\text{\%}}$ and FWHM for the simulated curve, and 3) estimation of $\left|{\text{ΔR}}_{80\text{\%}}\right|$ and $\left|\text{ΔFWHM}\right|$ comparing measurement and new simulation. The Fred beam energy (E) and energy spread (E_σ) are considered optimal when $\left|{\text{ΔR}}_{80\text{\%}}\right|$ and $\left|\text{ΔFWHM}\right|$ are less than or equal to 0.05 mm.

Step 3. In the third step (Figure 2: Step 3), the dosimetric calibration from TPS MU to the number of particles (SF${\text{SF}}_{\text{MU}}$) was obtained for each nominal energy, mimicking the measurement setup. For this purpose, a monoenergetic 10 $\times$ 10 ${\text{cm}}^2$ field in water was simulated with spot spacing 2.5 mm, 1 MU per spot and unitary MU scaling factor. The dose in the uniform field center at 2 cm depth in water, ${\text{D}}_{2\text{cm}}$, was derived from the simulation. The MU scaling factor (${\text{SF}}_{\text{MU}}$) was obtained as the ratio between ${\text{D}}_{2\text{cm}}$ obtained from commissioning measurement and Fred MC simulation.

The output of the characterization procedure is a list of beam model parameters per nominal energy and is stored in a text file. We developed a software tool that converts clinical TPS treatment plan into Fred input files using the beam model library (cf. Figure 2: Conversion and calculation of treatment plans). The parameters in between nominal energies are linearly interpolated, mimicking the procedures applied by TPS and beam line control system.

2.5 Validation in Homogeneous Media

This section describes how the beam model library was validated by comparing Fred simulations with measurements performed at each facility. We compared 1) lateral propagation of proton pencil beams, 2) treatment plans of dose cubes, and 3) patient QA treatment plans. The beam model validation steps are schematically illustrated in Figure 2 (lower row). The treatment plans were exported from TPS and converted from DICOM to Fred input file format. The QA treatment plans were simulated in Fred using ${10}^5$ protons per pencil beam. After simulation, the dose from each spot was scaled to the actual number of particles optimized in the treatment plan using dosimetric calibration (${\text{SF}}_{\text{MU}}$). This approach warranties the same statistical precision of calculation of dose delivered by each pencil beam, regardless of its weight in the treatment plan.

Lateral propagation of proton pencil beams. The measurements of lateral profiles of proton pencil beams at 100, 150, and 200 MeV were performed using the Lynx scintillating screen (IBA Dosimetry) in air for CCB and EMORY [46] at five positions behind the range shifter. The beam lateral profiles in solid phantoms were measured with Lynx in RW3 slab phantom for beam energies 100, 150, and 200 MeV at CCB and in PMMA slab phantom for beam energies 130, 180, and 240 MeV at EMORY.

Fred simulations for pencil beams were performed at the corresponding positions behind the range shifter in air and in solid phantoms. The transverse shape of the beam in X and Y directions was fitted with a single Gaussian fit, and the σ obtained from measurements and simulations were compared.

Spread Out Bragg Peak (SOBP). The longitudinal profiles of dose cubes (SOBPs) were measured 1) at CCB using a dosimetrically calibrated plane-parallel Markus chamber placed in a water phantom (sensitive volume 0.055 ${\text{cm}}^3$) with variable 0.1‐1 cm step length and 2) at EMORY using the Zebra detector (IBA Dosimetry) without dosimetric calibration. The QA treatment plans of dose cubes were optimized in clinical TPS aiming at achieving homogeneous biological dose of 1 Gy (RBE) and 4 Gy (RBE) at CCB and EMORY, respectively. All cubes had a lateral size of 10 $\times$ 10 ${\text{cm}}^2$. At CCB, dose cubes of 5 cm length (modulation) and variable range of 10, 15, 20, 25, and 30 cm without range shifter were optimized and evaluated. At EMORY, dose cubes of 10 cm length (modulation) and constant range of 15 cm without and with three range shifters of different thickness were investigated. For each measurement, the isocenter position in water was in the middle of SOBP, causing that the measurements were performed with air gaps ranging from 5 to 32 cm for EMORY and from 11.2 to 29.4 cm for CCB. Simulations of the SOBP plans were performed in a virtual water phantom. The measured SOBP dose profiles were compared with the profile extracted from three-dimensional (3D) dose calculation obtained from Fred MC simulations. Absolute dose comparison was performed for Markus chamber measurements conducted at CCB, whereas relative dose comparison was performed for Zebra measurements conducted at EMORY.

Patient QA. To evaluate the accuracy of Fred simulations, patient QA treatment plans were simulated in a virtual water phantom and compared with patient QA measurements routinely performed in the clinic. The comparison of TPS vs. measurement is also shown.

In CCB and EMORY, the MatriXX PT (IBA Dosimetry) is currently in use for patient QA [47]. MatriXX is a two-dimensional (2D) array of 1020 plane-parallel ionization chambers of 4 mm diameter arranged in a 32 $\times$ 32 grid with the distance between chambers of 7.62 mm. In both facilities, the MatriXX detector was calibrated to dose in water according to protocol proposed by the manufacturer. Patient QA measurements are typically performed at 3-5 depths at CCB and at 1-3 depths at EMORY. The measurement depths are selected by a medical physicist during the QA preparation process for each patient individually, to cover the entire treatment field. For EMORY, the air gap ranges from 5 to 22 cm, whereas for CCB, it ranges from 21.7 to 27.7 cm. The patient QA treatment plans of 74 patients (1077 measured layers, 967 without and 110 with range shifter) treated in Krakow and 13 patients (56 measured layers) treated in EMORY were evaluated. The dose distributions obtained from TPS and Fred calculations were compared to measured data by means of dose profile and gamma index (GI) analysis [48]. GI calculation tools implemented in PyMedPhys Python package [49] were used for evaluation. The 3D GI test (2 mm distance-to-agreement and 2% of local dose difference criteria, with the dose cutoff at 2% of the maximum dose) was used to compare 2D slice of dose field measurement (reference) with 3D Fred dose distribution calculation (evaluation).

2.6 Validation in Heterogeneous Media

The end-to-end experimental validation of Fred physics models, beam model, and CT calibration using a heterogeneous CIRS head-and-neck phantom (model 731-HN) [50] was performed in Krakow. The experimental setup is shown in Figure 3. The CIRS phantom consists of five materials equivalent to the following tissues/organs: brain, bone, larynx, trachea, sinus, teeth, and nasal cavities. One half of the phantom consists of single piece, and the other is sliced into three segments as shown in Figure 3A. The CIRS phantom was positioned in the treatment room using orthogonal X-ray imaging system and the phantom CT scan, following the clinical patient positioning procedure applied in Krakow. The irradiation plans of 10 $\times$ 10 ${\text{cm}}^2$ monoenergetic fields at nominal energies 100, 150, and 200 MeV were prepared in clinical TPS with and without range shifter. The dose distribution downstream from the CIRS phantom was measured using the MatriXX detector placed in the DigiPhant water phantom (IBA Dosimetry, see Section 2.5). Data were acquired in 5 mm water-equivalent steps yielding 3D dose distribution with lateral resolution of 7.62 mm and longitudinal resolution of 5 mm. Dose distributions were measured behind half CIRS head in water for nominal energies 150 and 200 MeV (cf. Figure 3B). The dose distribution was measured behind 1/6 slice of CIRS head in water-equivalent RW3 slab phantom using 100 MeV proton beam (IBA Dosimetry; cf. Figure 3C) because 100 MeV protons have insufficient range to traverse the half-head phantom to acquire dose distribution in water using MatriXX (with and without range shifter).

FIGURE 3

FIGURE 3. Schematic illustration of CIRS phantom (A) and setup used for experiment and Fred MC simulations (B, C). (A) CIRS head phantom sliced into one-piece half-head and the other half sliced further into three segments; (B) setup with half-head CIRS and MatriXX detector placed in water phantom; (C) setup with one slice of CIRS and MatriXX detector placed between water-equivalent RW3 solid phantom. Setup (B) was irradiated with monoenergetic field at nominal proton beam energy 150 and 200 MeV, whereas setup (C) at 100 MeV, all with and without range shifter.

The measurements were compared to Fred simulations of the experimental setup performed in the CT image of the CIRS and water phantoms. The CT image of CIRS phantom was acquired using the CT scanner (Siemens SOMATOM) calibrated for treatment planning in Krakow. The comparison of measured and simulated 3D dose distributions was performed using a 3D GI method.

2.7 Patient Data

A retrospective patient study was performed to investigate time performance of Fred as an independent, MC-based, proton dose computation tool and demonstrate its applicability for patient QA in the clinic. For this purpose, we referred our results to the TPS computations.

The 122 treatment plans (including boost plans) of 90 head and neck as well as brain patients treated at CCB from 2016 to 2018 and an example treatment plan of a patient treated in EMORY in 2019 [7] were simulated in Fred on CT geometries. The clinical CT images were sampled down to 1.5 $\times$ 1.5 $\times$ 1.5 ${\text{mm}}^3$ voxel size. The facility-specific clinical CT calibration curve obtained from stoichiometric calibration [36] was implemented in Fred. The CT calibration curve used in Fred contains information on the composition, relative stopping power (RSP) of protons, radiation length, and density of 93 materials. The density and RSP of CT numbers between 93 predefined points are linearly interpolated. The CT images of the patient anatomy and delineated contours were used for the optimization of plans in clinical TPS using an analytical intensity modulated proton therapy (IMPT) optimization algorithm. Depending on the target size and the number of fields, the number of pencil beams in a treatment plan varied from 1,378 to 32,290 with the median value of 10,989. ${10}^4$ protons per pencil beam were simulated for each patient treatment plan recalculated in Fred, and the obtained dose was scaled to the actual number of particles optimized in the treatment plan. In order to investigate the impact of PTV volume on the Fred dose calculation accuracy, we divided treatment plans of patients treated in CCB into three subgroups distinguishing 12 plans with small PTV volume (${\text{V}}_{\text{PTV}}<$ 50 ml), 60 plans with medium PTV volume (50 ml $\le {\text{V}}_{\text{PTV}}<$ 200 ml), and 50 plans with large PTV volume (${\text{V}}_{\text{PTV}}\ge$ 200  l). PTV volumes ranged from 28.5 ml to 1,010 ml

An example treatment planning study on 122 plans included a comparison of dose distributions obtained from Fred and from clinical TPS. We evaluated four parameters based on dose volume histogram (DVH) that characterize the quality of dose distribution. 1) The mean dose (${\text{D}}_{\text{mean}}$) is related to the prescribed dose (${\text{D}}_{\text{p}}$). 2) The homogeneity index (HI) characterizes the slope of the DVH; hence, the uniformity of the dose distribution in the PTV. The HI is defined as $\text{HI}=\left({\text{D}}_{2\text{\%}}-{\text{D}}_{98\text{\%}}\right)/{\text{D}}_{\text{p}}$, where ${\text{D}}_{2\text{\%}}$ and ${\text{D}}_{98\text{\%}}$ are the doses received by 2% and 98% of the PTV, respectively [51]. 3) The conformity index (CI) describes how much dose prescribed to the planning target volume (PTV) is delivered outside the PTV, possibly to organs at risk. The CI is defined as $\text{CI}={\text{V}}_{95\text{\%}}^{\text{body}}/{\text{V}}_{95\text{\%}}^{\text{PTV}}$, where ${\text{V}}_{95\text{\%}}^{\text{body}}$ and ${\text{V}}_{95\text{\%}}^{\text{PTV}}$ are the volumes of the body and PTV, which receive at least 95% of the prescribed dose ${\text{D}}_{\text{p}}$ [52]. 4) The relative mean square error (RMSE) characterizes the deviation of a DVH from the prescribed dose ${\text{D}}_{\text{p}}$. It was calculated at the slope of a DVH, in a range between $D_{5\text{\%}}$ and $D_{95\text{\%}}$, and it is defined as $\text{RMSE}=\sqrt{{{{\displaystyle \sum }}^{\text{}}}_5^{95}{\left({\text{D}}_{\text{x\%}}-{\text{D}}_{\text{p}}\right)}^2/90}$.

3 Results

3.1 Generation of the Beam Model Parameter Library

The beam model parameter libraries characterizing the proton beam model for CCB and EMORY facilities were generated using an automated procedure (cf. Section 2.3) and are illustrated as a function of nominal proton beam energy in Figure 4. Using the beam model library, the nominal primary proton beam energy for each pencil beam from the treatment plan is used to define the initial parameters of the pencil beams used by Fred simulations. Figure 4 (top-left panel) shows a linear relation between the nominal proton beam energy used by TPS and Fred. The energy spread values fluctuate within 1 MeV and are slightly smaller for Krakow than for Atlanta proton center. Figure 4 (top-right panel) shows the dosimetric scaling factors used to convert MU to the number of primary particles per pencil beam spot. The bottom panels of Figure 4 show the six parameters of emittance model used for Krakow (bottom-left panel) and the four parameters of VPS model used for Atlanta facility, characterizing the lateral beam propagation (bottom-right panel). The lateral asymmetry of the pencil beams in X (filled circles) and Y (empty circles) directions is taken into account in the beam model characterization.

FIGURE 4

FIGURE 4. The parameters characterizing proton beam model used in CCB and EMORY facilities at the entire primary proton beam energy range. Nominal energy corresponds to energy used by clinical TPS. Top-left panel: Beam energy and energy spread; Top-right panel: dosimetric calibration; bottom-left panel: emittance model parameters used in CCB; bottom-right panel: VPS model parameters used in EMORY.

The IDD profiles of single proton beams in water for three nominal energies: 100, 150, and 200 MeV are given in Figures 5 A and B for the Krakow and Atlanta facilities. The profiles are in agreement with the commissioning measurements: the range (${\text{R}}_{80\text{\%}}$) of the pencil beams agrees within 0.02 mm, the relative dose difference along the pencil beam profile is below 4%, the FWHM of the Bragg peak agrees within 0.05 mm, the distal falloff width between 80% and 20% Bragg peak dose agrees within 0.04 mm, and the peak-to-plateau ratio agrees within 0.11.

FIGURE 5

FIGURE 5. Examples of longitudinal proton beam propagation in water (top panels) and lateral proton beam propagation (σ) in X and Y directions in air (bottom panels) for CCB (left) and EMORY (right) facilities at three proton beam energies: 100, 150, and 200 MeV. Depth dose distribution profiles of proton pencil beams simulated with beam model parameters in Fred (FRED Bragg peak) and obtained experimentally during the facility commissioning (measured Bragg peak) for CCB (panel A) and EMORY (panel B). The transverse shape and velocity evolution of the proton beam represented by means of the emittance model for CCB (panel C) and VPS model for EMORY (panel D).

The fitted single beam sizes in air obtained in commissioning measurements, described by ${\sigma }_{x/y}$ of lateral pencil beam profiles is shown in the Figures 5 C and D for three nominal energies: 100, 150, and 200 MeV for the Krakow and Atlanta facilities, respectively. The maximum absolute difference between fitted and measured beam sizes ranging from −20 to 20 cm (CCB) and −30 to 5 cm (EMORY) in Z direction with respect to the isocenter is smaller than 0.05 mm. We deem this sufficiently accurate to model lateral beam propagation in clinical applications. The quadratic and linear shape of the fit justifies the use of the emittance (Figure 5C) and VPS (Figure 5D) model for the Krakow and Atlanta facilities, respectively.

Dose computation time for a single pencil beam at 100, 150, and 200 MeV simulated with ${10}^8$ primary protons was 36, 44, and 53 s, respectively. The corresponding tracking rate is 10.1 $\times$${10}^6$, 5.7 $\times$${10}^6$, and 3.6 $\times$${10}^6$${\text{p}}^{+}/\text{s}$. The tracking rate decreases with the beam range as more interactions must be processed.

The total computation time needed to determine the beam model parameters for all reference energies following the automated procedure described in Section 2.3 was approximately 12 h. Within this time, 1) the parameters characterizing beam lateral propagation were fitted (Figure 2 step 1; total time: few seconds), 2) simulations required for E and ${\text{E}}_{\sigma }$ optimization were performed and the optimization procedure itself was executed (Figure 2 step 2; total time: approximately 10 h), and 3) simulations of monoenergetic 10 $\times$ 10 ${\text{cm}}^2$ fields required for ${\text{SF}}_{\text{MU}}$ calculation were performed (Figure 2 step 3; total time: approx. 2 h). For CCB, full-beam model characterization required a total of 303 Fred MC simulations, including 286 simulations for E and ${\text{E}}_{\sigma }$ optimization and 17 simulations for SF${\text{SF}}_{\text{MU}}$ calculation (average time of single simulation was approximately 2 and 7 min, respectively).

3.2 Validation in Homogeneous Media

Lateral propagation of proton pencil beams. The lateral propagation of pencil beams in air behind range shifter of different thickness (Figure 6) and in slab phantoms (Figure 7) was simulated in Fred and compared with the beam size ${\sigma }_{x/y}$ of lateral pencil beam profiles obtained experimentally. Note that the comparison was performed at different positions/depths and for different primary proton beam energies at CCB and EMORY facilities.

FIGURE 6

FIGURE 6. Spot sizes in air in X (blue) and Y (red) directions for CCB and EMORY without range shifter and behind the range shifters used at facility (single range shifter (RS) of thickness 36.7 mm for CCB and RS2, RS3, and RS5 of thickness 20, 30, and 50 mm, respectively, for EMORY). The measured spot sizes are shown as points with error bars ($\pm 0.1$ mm), and the solid lines show the simulation results.

FIGURE 7

FIGURE 7. The transverse shape evolution (σ) of proton pencil beam measured and simulated in water equivalent slab phantom.

The lateral propagation of the beam in range shifter and in slab phantom is accurately modeled in Fred. The values of ${\sigma }_{x/y}$ obtained from measurements agree with simulated values mostly within 100 $\mu \text{m}$, as indicated by error bars in Figures 6 and 7. The results in air and in slab phantoms are within the spot size QA acceptance criterion of $\pm$0.6 mm used by CCB therapy center.

Spread Out Bragg Peak (SOBP). Depth dose distribution profiles of cubic volumes obtained from measurements and Fred simulations are shown in Figure 8 for CCB in the top panels and for EMORY in the bottom panels. The results obtained for CCB are absolute dose, whereas they are relatively normalized to the dose value in the middle of the SOBP for EMORY. Because the treatment plans were optimized in clinical TPS, the obtained physical dose differs from the prescribed biological dose by the RBE factor of 10%.

FIGURE 8

FIGURE 8. Dose profiles of cubic volumes of SOBP obtained from Fred MC calculations (solid line) and measurements (dots) for CCB (top panel) and EMORY (bottom panel) facilities. The relative dose difference between the measurement and simulation is illustrated by crosses.

Good agreement between Fred MC simulations and dose measurements along the SOBP profiles was obtained. The maximum relative dose difference is 2% for most of the measurement points. The largest relative dose differences are observed at the distal falloff, that is, a high-dose gradient region, and result from the detector positioning uncertainties, estimated to be about $\pm$0.3 mm. Small variations between the measurements and simulations are present at the beginning of the plateau and in the SOBP of cubes between the range of 25 and 30 cm. They are potentially related to the implementation of the nuclear interaction model in Fred for the highest beam energies. This accuracy is acceptable for the scope of the presented clinical application.

The tracking rate of the dose cube simulation ranged from 4.5 $\times$${10}^6$ to 2.0 $\times$${10}^6$${\text{p}}^{+}/\text{s}$ and the complete dose computation time for a single dose cube was up to 10 min, with the statistics ${10}^5$ primaries per pencil beam.

Patient QA. 2D transversal dose maps obtained from measurements performed with the MatriXX detector in water phantom were compared with Fred and TPS simulations of patient treatment plans using the GI analysis. Data from 1077 measurements performed at CCB and 52 measurements performed at EMORY were investigated, and the results of the comparison are summarized in Figure 9. The average GI passing rate obtained comparing all simulated and measured layers was 97.83% (4.94) (1σ) for CCB and 95.51% (3.88) (1σ) for EMORY. Of 1,077 layers evaluated for CCB, 1,022 fulfilled the requirement for the GI passing rate (%GP) to be greater than 90%. For EMORY, 47 of 52 investigated layers fulfilled this requirement. Figure 10 shows an example of a transversal dose field layer extracted from Fred MC simulation and the corresponding dose distribution measured with MatriXX at the same depth in water, as well as in the GI map.

FIGURE 9

FIGURE 9. A transversal 2D dose distribution layer measured with an array of ionization chambers in water phantom (left panel), obtained from Fred MC simulations (middle panel) and a GI map computed comparing Fred simulation and measurement using GI (2%/2 mm) method (right panel). GI passing rate is $99.53\text{\%}$ for the CCB case shown in the top panels and $95.95\text{\%}$ for EMORY case shown in bottom panel.

FIGURE 10

FIGURE 10. Evaluation of gamma index passing rate (%GP) for 2D dose maps obtained from patient QA of 1,077 layers measured in CCB (left and right panels) and of 52 layers measured in EMORY (middle panel). Red and blue box plots correspond to the distribution of %GP obtained from the comparison of measurements to TPS and Fred calculations, respectively. In the left and middle panels, we compared the layers planned with range shifter (RS) and without range shifter (NRS), whereas in the right panel, small (${\text{V}}_{\text{PTV}}<$ 50 ml), medium (50 ml $\le {\text{V}}_{\text{PTV}}<$ 200 ml), and large (${\text{V}}_{\text{PTV}}\ge$ 200 ml) PTV volumes. Green numbers labeled as “pass” stand for the number of cases passing %GP$>$90% criterion, whereas “total” is the population of a given group.

For a patient verification treatment plan, the average tracking rate and complete dose computation time were 3.4(0.4)$\times$${10}^6$${\text{p}}^{+}/\text{s}$ (1σ) and 2’34 (1’38) min (1σ), respectively.

3.3 Validation in Heterogeneous Media

The experimental validation of Fred accuracy was performed by comparing 3D dose distributions behind the heterogeneous phantom obtained experimentally and from Fred simulations (cf. Section 2.6). An example of the comparison of Fred simulation against the experimentally acquired data is shown in Figure 11. Two 3D dose measurements, one with and other without range shifter, were performed for each of the investigated energies (100, 150, 200 MeV). An excellent agreement between Fred simulations and measurements was achieved. For all the investigated cases, the 3D GI (2%/2 mm) is greater than 99%. Comparing the clinical (analytical) TPS simulation and the measurements, the GI passing rate is ${93.2}_{76.3}^{98.0}\text{\%}$ (σ = 8.4%). See the Supplementary Material of the article for detailed results of other measurements performed at 100 and 200 MeV, with and without range shifter.

FIGURE 11

FIGURE 11. The experimental validation of Fred simulations in heterogeneous CIRS phantom. Panel (A): measurement of 3D dose distribution in water phantom performed using MatriXX. Panel (B): Fred simulation of 3D dose distribution. Panel (C): 2D GI map (2%/2 mm) obtained comparing experiment to Fred simulations. The color maps on panels A–C are overlaid on CT scan of CIRS and water phantom. Panels (D) and (E) show longitudinal and lateral profiles, respectively, obtained from measurements (dots) and simulations (solid line). See Supplementary Material for the complete report of the validation.

3.4 Example Clinical Application of Fred

As an example, dose distributions, dose profiles, and DVHs recalculated with Fred and clinical TPS, for one patient case from CCB and one from EMORY, are shown in Figure 12. For CCB patient case (Figure 12 top panels), dose distributions computed with Fred are less uniform compared to the analytical TPS calculations. This is also observed analyzing the dose profiles and the DVH for PTV and results in the reduction of the mean dose in PTV and organ at risk. For EMORY patient case (Figure 12 bottom panels), the differences in dose distributions are less visible as MC-based TPS was used for the dose optimization and calculation. The observed differences between Fred and RayStation MC-based TPS are similar to the results obtained comparing RayStation with ECLIPSE MC algorithm reported by Chang et al. [7].

FIGURE 12

FIGURE 12. The evaluation of the treatment plan of patient treated at CCB (top panels) and at EMORY (bottom panels). On the left panels, dose distributions computed with clinical TPS and Fred are shown. PTV (black solid line) and 95% isodose (blue dashed line) are delineated. The corresponding dose profiles and DVHs are shown in top-right and bottom-right panels, respectively.

Analysis of 122 treatment plans of patients treated at CCB was performed to quantify the time performance and demonstrate the clinical applicability of Fred dose computations for patient QA. Comparing dose distributions in PTV, we observed that the ratio ${\text{D}}_{\text{mean}}/{\text{D}}_{\text{p}}$ obtained with Fred is more dispersed than the one obtained with analytical TPS, while the effect is more pronounced for small targets. The average relative difference in median value ranges from 3% for small targets, through 1.5% for medium size target volumes, to 1% for large target volumes, as shown in Figure 13 (left panel). The analysis of HI in PTV is shown in Figure 13 (middle-left panel). On average, the median HI is 0.11 and 0.16 for clinical TPS and Fred, respectively. Independently on the target volume, the HI in PTV calculated with Fred is higher, that is, dose distribution is less homogeneous than the HI calculated with analytical TPS. Figure 13 (middle-right panel) shows the CI distributions, which present no substantial difference between both, Fred and TPS calculations (median CI is 1.26 and 1.23 for TPS and Fred, respectively). In general, for both, Fred and TPS calculations, dose distributions of small PTV are less conformal with respect to dose distributions for large PTV. The comparison of DVH in PTV by means of RMSE analysis confirms the conclusions from ${\text{D}}_{\text{mean}}/{\text{D}}_{\text{p}}$ ratio and HI analysis. The histogram of RMSE for TPS distribution is narrower with smaller mean value, whereas for Fred, the RMSE distribution is wider with slightly greater mean value. This is because the dose distributions calculated with Fred are less uniform in PTV, as indicated by HI analysis, and the mean dose in PTV differs from the dose in PTV calculated with TPS, as indicated by ${\text{D}}_{\text{mean}}/{\text{D}}_{\text{p}}$ ratio analysis.

FIGURE 13

FIGURE 13. The parameters characterizing the quality of 122 dose distribution obtained from patient treatment plans computed with clinical TPS (blue) and Fred (red) for small (${\text{V}}_{\text{PTV}}<$ 50 ml), medium (50 ml $\le {\text{V}}_{\text{PTV}}<$ 200 ml), and large (${\text{V}}_{\text{PTV}}\ge$ 200 ml) PTV. The left panel shows the ratio ${\text{D}}_{\text{mean}}/{\text{D}}_{\text{p}}$, the middle-left panel shows the homogeneity index (HI), the middle-right panel the conformity index (CI), and the right panel shows relative mean square error between a prescribed dose and DVH in PTV computed with Fred and clinical TPS (the dark blue area depicts overlapping of two histograms).

For a treatment plan, the total simulation time varied depending on the complexity of the plan, that is, the total number of pencil beams and the presence of range shifter in the plan. For the simulations in CT geometry rescaled to 1.5 $\times$ 1.5 $\times$ 1.5 ${\text{mm}}^3$ voxels, the computation time ranged from 21 s to 6′26 min (average value 2′28 (1′25) min (1σ)) with the average tracking rate of 2.9 (1.1)$\times$${10}^5$${\text{p}}^{+}/\text{s}$ (1σ).

4 Discussion

We have built a proton beam model libraries for Fred MC code according to the QA protocols, and we accomplished acceptance tests required for beam model validation in a commercial TPS at proton therapy facilities. We performed MC commissioning avoiding the nozzle geometry modeling, similar to the work presented by other groups [10, 32, 33]. The beam model library parameters containing the information on initial proton energy and energy spread, lateral beam propagation, and dosimetric calibration were identified in 10 MeV energy steps in the therapeutic energy range to best fit the commissioning measurements of proton pencil beams (cf. Section 3.1). A submillimeter agreement between simulated and measured Bragg peaks shape and range in water and lateral beam sizes in air and in solid phantoms was obtained with and without range shifter for beam model of two facilities of different beam line design.

In the study, we assumed the uncertainty of single pencil beam and SOBP depth dose profile measurements to be $\pm$3%. The uncertainty of positioning of the ionization chamber in the water phantom is about 0.3 mm. The uncertainty of the lateral pencil beam size measurement performed with scintillating screen (Lynx detector) in air and in the RW3/PMMA slab phantom is $\pm$0.1 mm, whereas the measurement with IC23 has 0.5 mm uncertainty [32]. We estimate the uncertainty of the slab phantom positioning at 1 mm, but it has negligible impact on the beam lateral profile measurements.

We performed beam model commissioning and validation using the proton per pencil beam statistics that it is required to assure no impact of the statistical uncertainty on these results. For single pencil beams, ${10}^8$ protons per beam offer statistical uncertainty below 1% in 3σ distance from the beam core, when simulations are performed in 1 $\times$ 1 $\times$ 1 ${\text{mm}}^3$ grid. Lower statistics can be used for recalculation of treatment plans in water in the same 1 $\times$ 1 $\times$ 1 ${\text{mm}}^3$ grid because dose distribution is obtained from superposition of hundreds of pencil beams. We found that for treatment plan recalculation in water, the statistical uncertainty below 1% can be achieved using ${10}^5$ protons per beam for small fields. For clinical application of Fred, the limiting factor is the time of simulations. We found that for resampling the patient geometry in CT to 1.5 $\times$ 1.5 $\times$ 1.5 ${\text{mm}}^3$, ${10}^4$ primaries per pencil beams can be used, achieving statistical uncertainty of about 2%. We consider this setting as a good compromise between simulation time and simulation accuracy, allowing treatment plan recalculation in CT scan within a few minutes. No statistical uncertainty of the dose calculated with analytical TPS used as CCB was considered, whereas the dose was calculated with the statistical uncertainty of 0.5% in MC-based RayStation TPS on 2 $\times$ 2 $\times$ 2 ${\text{mm}}^3$ grid in water phantom and $3\times$3$\times$3 ${\text{mm}}^3$ grid in patient CT.

The comparison of Fred simulations to QA measurements in water presented in Section 3.2 indicates that, on average, Fred dose distributions agree better with measurements than the prediction made by TPS pencil beam algorithm used in CCB (Figure 9, left panel); however, Fred dose distributions are comparable to predictions of commercial MC-based TPS used in EMORY (Figure 9, middle panel). Analysis of CCB patient QA data shows that for small, medium, and large PTV volumes, on average, the dose distributions computed by Fred agree better with measurements when compared with dose distributions computed with pencil beam algorithm (Figure 9, right panel). We have not observed substantial differences in Fred dose calculation accuracy between different PTV volume categories. Note that small PTV volumes ranging from 28.5 to 50 ml were investigated for CCB. In Section 3.3, we presented the results of end-to-end Fred validation of Fred simulations. For various beam energies, large air gaps, and setups with and without range shifter, we compared Fred simulations with measurements of 3D dose distributions behind anthropomorphic CIRS head phantom containing high-density gradients on the boundary between head bones and nasal cavities. The high accuracy of the Fred dose calculations was confirmed in the results of GI tests better than 99% for all of the investigated cases.

Comparison of experimental results in homogeneous media and anthropomorphic phantom with Fred simulations (cf. Sections 3.2 and 3.3 and Supplementary Materials) indicates that fast dose recalculations in patient CT performed with Fred (cf. Section 3.4) is a very accurate simulation of proton treatment. A retrospective treatment planning study and the statistical evaluation of DVH parameters are example of routine clinical application of Fred for patient QA. The dose nonuniformities in PTV shown in an example CCB patient case recalculated with Fred (Figure 12) are also observed in the analysis of ${\text{D}}_{\text{mean}}/{\text{D}}_{\text{p}}$ and HI for 122 patient cases summarized in Figure 13. The differences of the mean dose delivered to PTV structures, calculated by Fred and predicted by TPS are more pronounced for small PTV volumes (Figure 13, left panel). Fred calculations predict dose nonuniformity for small, medium, and large PTV volumes, which cannot be calculated with analytical TPS used in Krakow. In general, dose distributions are less conformal in small targets than in large targets because it is predicted both by Fred simulations and by TPS pencil beam algorithm calculations (Figure 13, right panel). Note that these clinical results, both from TPS and Fred, include uncertainties related to acquisition of commissioning data, beam model implementation, CT calibration, and the like. On the other hand, the distribution of ${\text{D}}_{\text{mean}}/{\text{D}}_{\text{p}}$, HI, and CI indicate that overall, the dose distribution calculations performed with both clinical TPS and Fred are within the clinically relevant acceptance.

In clinical practice, additional information about dose, LET and RBE-weighted dose distributions calculated with Fred can be an indication for medical physicists to revise the treatment plan optimization or to perform additional experimental validation, when the results deviate from the predictions of TPS exceeding acceptance criteria. The time performance of Fred enables to obtain this information within about 2.5 min. Fred is currently adapted to be executed as a stand-alone library, which will enable its easy integration with commercial TPS (eg, Eclipse or RayStation) and dedicated software tools for patient QA (eg, MyQAion).

Schiavi et al. [34] reported that simulation of dose deposition in a water phantom induced by ${10}^6$ primary protons can be reduced from 22 min required by FLUKA MC code to 0.5 s when employing Fred running on two GPU modules [34]. Regarding dose distribution simulation in patients, Grassberger, Anthony Lomax, and Paganetti [33] reported that the patient simulation for the head and neck took 371 min (${10}^6$ primaries simulated) on single CPU using TOPAS (Geant4), which corresponds to a tracking rate of 45 ${\text{p}}^{+}/\text{s}$, whereas the average tracking rate obtained with Fred is 2.9 $\times$${10}^5$${\text{p}}^{+}/\text{s}$ in patient CT rescaled to 1.5 $\times$ 1.5 $\times$ 1.5 ${\text{mm}}^3$ using two GPUs. The time performance results presented in this article can be linearly scaled as a function of the number of GPU cards applied [34]. Note that the simulation time depends on the number of primaries simulated per pencil beam, tumor depth (i.e. the beam energy), and scoring resolution used for the simulation. The most accurate dose calculations in tissue heterogeneities can be obtained performing the simulation in original CT grid. In order to achieve the statistical uncertainty below 1% on CT grid used at CCB 0.7 $\times$ 0.7 $\times$ 1.2 ${\text{mm}}^3$, ${10}^5$ primaries per pencil beam should be simulated. The average simulation time for the patient group investigated in Section 3.4 in original CT resolution is ${31.8}_{3.5}^{161.8}$ (σ = 23.8) min.

The clinical application of proton therapy and development of new treatment protocols, for example, studies on the reduction of safety margins accounting for treatment plan robustness, require treatment planning studies that can only be performed analyzing several treatment planning approaches. The total simulation time of all 122 patient cases shown in Section 3.4 was about 5 h. An example study of 10 possible treatment planning approaches on our patient group could be performed using Fred within about two days of simulation. Another application is robust optimization of treatment plans, that is, particularly relevant for treatment planning of moving targets, when several dose distributions must be computed on 4D CT. Performing such studies without the time performance offered by Fred would not be possible with any general purpose MC code in reasonable time.

In addition to its clinical applications, the time performance of Fred enables preparation of the proton beam model faster with respect to a general purpose MC codes. This is particularly useful when a new beam model must be implemented in the clinical routine due to technical modifications or maintenance at accelerator. When the facility beam commissioning measurements are available, the GPU acceleration offered by Fred allows to parametrize the beam model within about 12 h, requiring minimal manual interventions. This potentially enables easy and quick use of Fred for research and patient QA purposes in most of the proton facilities with little experimental efforts.

5 Conclusion

In this article, we share our experience on commissioning and validation of GPU-accelerated MC code Fred based on commissioning measurements of two proton beam therapy facilities of different beam line design: CCB (Krakow) from IBA and EMORY (Atlanta) from Varian. Fred passed acceptance tests required to approve TPS for clinical use. The approach we used combines the application of a new GPU-accelerated MC code, implementation of two proton beam lateral beam propagation models, automated beam model optimization method, experimental validation of beam model parameters in an anthropomorphic phantom with and without range shifter, and comparison of patient treatment plans computed with Fred and clinical TPS in patient CT. Our commissioning and validation results demonstrate the universal and accurate implementation of the physics models in Fred, allowing its flexible applications for medical physics and research purposes. The application of Fred as a secondary MC engine for patient QA in clinical routine is foreseen in Krakow proton facility. Fred is currently used for treatment planning studies evaluating radiobiologically effective dose using variable RBE.

Data Availability Statement

The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Author Contributions

JG, MG, AS, ASc, and AR developed the beam model for CCB. JG developed automated beam model library implementation method and performed data analysis to validate the beam model. JG, MG and ASc developed the emittance and virtual point source models for Fred. JG, NM and AR designed, while JG, MG, NM, AR, MR performed validation experiments with proton beams at CCB. JG performed data analysis of experiments. KC, NK and MPN supported data analysis. RK provided access to beam model commissioning and validation data from CCB. JG participated in commissioning measurements in CCB. CC and LL provided commissioning, validation, and patient data from EMORY. JG implemented beam model for EMORY and performed analysis of validation and patient data. RK and EP provided access to patient data from CCB. KK and MR exported the patient data from clinical TPS. MG and JG performed simulations and analysis of patient data. ASc and VP developed and made substantial improvements in Fred source code required to enable presented studies. MD, IR, ES, and FT provided expertise in beam modeling and medical physics. JG prepared all figures. JG and AR designed the project and drafted the manuscript. AR acquired funding. All the authors reviewed and approved the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This project is carried out within the Reintegration program of the Foundation for Polish Science cofinanced by the EU under the European Regional Development Fund—grant no. POIR.04.04.00-00-2475/16-00. MG acknowledge the support of InterDokMed project no. POWR.03.02.00-00-I013/16. This research was supported in part by computing resources of ACC Cyfronet AGH. We acknowledge the support of NVIDIA Corporation with the donation of the GPU used for this research. We acknowledge Aleksander Krempa from CCB Krakow proton therapy center for IT support during implementation of this project.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2020.567300/full#supplementary-material.

References

1. Paganetti H, Jiang H, Parodi K, Slopsema R, Engelsman M. Clinical implementation of full Monte Carlo dose calculation in proton beam therapy. Phys Med Biol 53 (2008). 4825.

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Saini J, Maes D, Egan A, Bowen SR, James SS, Janson M, et al. Dosimetric evaluation of a commercial proton spot scanning Monte-Carlo dose algorithm: comparisons against measurements and simulations. Phys Med Biol 62 (2017). 7659.

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Widesott L, Lorentini S, Fracchiolla F, Farace P, Schwarz M. Improvements in pencil beam scanning proton therapy dose calculation accuracy in brain tumor cases with a commercial Monte Carlo algorithm. Phys Med Biol 63 (2018). 145016. doi:10.1088/1361-6560/aac279

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Tommasino F, Fellin F, Lorentini S, Farace P. Impact of dose engine algorithm in pencil beam scanning proton therapy for breast cancer. Phys Med: Eur J Med Plants 50 (2018). 7–12. doi:10.1016/j.ejmp.2018.05.018

CrossRef Full Text | Google Scholar

5. Paganetti H. Range uncertainties in proton therapy and the role of Monte Carlo simulations. Phys Med Biol 57 (2012). 99–117. doi:10.1088/0031-9155/57/11/R99

CrossRef Full Text | Google Scholar

6. Langner UW, Mundis M, Strauss D, Zhu M, Mossahebi S. A comparison of two pencil beam scanning treatment planning systems for proton therapy. J Appl Clin Med Phys 19 (2018). 156–63. doi:10.1002/acm2.12235

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Chang CW, Huang S, Harms J, Zhou J, Zhang R, Dhabaan A, et al. A standardized commissioning framework of Monte Carlo dose calculation algorithms for proton pencil beam scanning treatment planning systems. Med Phys (2020). doi:10.1002/mp.14021

CrossRef Full Text | Google Scholar

8. Trnková P, Knäusl B, Actis O, Bert C, Biegun AK, Boehlen TT, et al. Clinical implementations of 4D pencil beam scanned particle therapy: report on the 4D treatment planning workshop 2016 and 2017. Phys Med 54 (2018). 121–30. doi:10.1016/j.ejmp.2018.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Guterres Marmitt G, Pin A, Ng Wei Siang K, Janssens G, Souris K, Cohilis M, et al. Platform for automatic patient quality assurance via Monte Carlo simulations in proton therapy. Phys Med 70 (2020). 49–57. doi:10.1016/j.ejmp.2019.12.018

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Fracchiolla F, Lorentini S, Widesott L, Schwarz M. Characterization and validation of a Monte Carlo code for independent dose calculation in proton therapy treatments with pencil beam scanning. Phys Med Biol 60 (2015). 8601–19. doi:10.1088/0031-9155/60/21/8601

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Winterhalter C, Fura E, Tian Y, Aitkenhead A, Bolsi A, Dieterle M, et al. Validating a Monte Carlo approach to absolute dose quality assurance for proton pencil beam scanning. Phys Med Biol 63 (2018). doi:10.1088/1361-6560/aad3ae

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Johnson JE, Beltran C, Wan Chan Tseung H, Mundy DW, Kruse JJ, Whitaker TJ, et al. Highly efficient and sensitive patient-specific quality assurance for spot-scanned proton therapy. PloS One 14 (2019). doi:10.1371/journal.pone.0212412

CrossRef Full Text | Google Scholar

13. Zhu XR, Li Y, Mackin D, Li H, Poenisch F, Lee AK, et al. Towards effective and efficient patient-specific quality assurance for spot scanning proton therapy. Cancers 7 (2015). 631–47. doi:10.3390/cancers7020631

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Matter M, Nenoff L, Meier G, Weber DC, Lomax AJ, Albertini F. Alternatives to patient specific verification measurements in proton therapy: a comparative experimental study with intentional errors. Phys Med Biol 63 (2018). doi:10.1088/1361-6560/aae2f4

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Battistoni G, Boehlen T, Cerutti F, Chin PW, Esposito LS, Fassò A, et al. Overview of the FLUKA code. Ann Nucl Energy 82 (2015). 10–8. doi:10.1016/j.anucene.2014.11.007

CrossRef Full Text | Google Scholar

16. Allison J, Amako K, Apostolakis J, Araujo H, Dubois PA, Asai M, et al. Geant4 developments and applications. Nuclear Science, IEEE Trans Nucl Sci 53 (2006). 270–8. doi:10.1109/TNS.2006.869826

CrossRef Full Text | Google Scholar

17. Jarlskog CZ, Paganetti H. Physics settings for using the Geant4 toolkit in proton therapy. Nuclear Science, IEEE Transactions on Nuclear Science 55 (2008). 1018–25.

CrossRef Full Text | Google Scholar

18. Henkner K, Sobolevsky N, Jäkel O, Paganetti H. Test of the nuclear interaction model in SHIELD-HIT and comparison to energy distributions from GEANT4. Phys Med Biol 54 (2009). N509.

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Jan S, Santin G, Strul D, Staelens S, Assié K, Autret D, et al. GATE: a simulation toolkit for PET and SPECT. Phys Med Biol 49 (2004). 4543.

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Jan S, Benoit D, Becheva E, Carlier T, Cassol F, Descourt P, et al. Gate V6: a major enhancement of the GATE simulation platform enabling modelling of CT and radiotherapy. Phys Med Biol 56 (2011). 881.

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Sarrut D, Bardiès M, Boussion N, Freud N, Jan S, Létang J, et al. A review of the use and potential of the GATE Monte Carlo simulation code for radiation therapy and dosimetry applications. Med Phys 41 (2014). 64301. doi:10.1118/1.4871617

CrossRef Full Text | Google Scholar

22. Perl J, Shin J, Schumann J, Faddegon B, Paganetti H. TOPAS: an innovative proton Monte Carlo platform for research and clinical applications. Med Phys 39 (2012). 6818–37. doi:10.1118/1.4758060

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Testa M, Schümann J, Lu HM, Shin J, Faddegon B, Perl J, et al. Experimental validation of the TOPAS Monte Carlo system for passive scattering proton therapy. Med Phys 40 (2013). 121719. doi:10.1118/1.4828781

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Jia X, Schümann J, Paganetti H, Jiang SB. GPU-based fast Monte Carlo dose calculation for proton therapy. Phys Med Biol 57 (2012). 7783–97. doi:10.1088/0031-9155/57/23/7783

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Qin N, Botas P, Giantsoudi D, Schuemann J, Tian Z, Jiang SB, et al. Recent developments and comprehensive evaluations of a GPU-based Monte Carlo package for proton therapy. Phys Med Biol 61 (2016). 7347.

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Giantsoudi D, Schuemann J, Jia X, Dowdell S, Jiang S, Paganetti H. Validation of a GPU-based Monte Carlo code (gPMC) for proton radiation therapy: clinical cases study. Phys Med Biol 60 (2015). 2257.

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Wan Chan Tseung H, Ma J, Beltran C. A fast GPU-based Monte Carlo simulation of proton transport with detailed modeling of nonelastic interactions. Med Phys 42 (2015). 2967–78. doi:10.1118/1.4921046

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Mein S, Choi K, Kopp B, Tessonnier T, Bauer J, Ferrari A, et al. Fast robust dose calculation on GPU for high-precision 1H, 4He, 12C and 16O ion therapy: the FRoG platform. Sci Rep 8 (2018). 14829. doi:10.1038/s41598-018-33194-4

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Choi K, Mein SB, Kopp B, Magro G, Molinelli S, Ciocca M, et al. FRoG - a new calculation engine for clinical investigations with proton and carbon ion beams at cnao. Cancers 10 (2018). doi:10.3390/cancers10110395

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Ma CMC, Chetty IJ, Deng J, Faddegon B, Jiang SB, Li J, et al. Beam modeling and beam model commissioning for Monte Carlo dose calculation-based radiation therapy treatment planning: report of AAPM Task Group 157. Med Phys 47 (2020). e1–e18. doi:10.1002/mp.13898

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Parodi K, Mairani A, Brons S, Hasch BG, Sommerer F, Naumann J, et al. Monte Carlo simulations to support start-up and treatment planning of scanned proton and carbon ion therapy at a synchrotron-based facility. Phys Med Biol 57 (2012). 3759.

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Grevillot L, Bertrand D, Dessy F, Freud N, Sarrut D. A Monte Carlo pencil beam scanning model for proton treatment plan simulation using GATE/GEANT4. Phys Med Biol 56 (2011). 5203.

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Grassberger C, Lomax A, Paganetti H. Characterizing a proton beam scanning system for Monte Carlo dose calculation in patients. Phys Med Biol 60 (2015). 633. doi:10.1088/0031-9155/60/2/633

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Schiavi A, Senzacqua M, Pioli S, Mairani A, Magro G, Molinelli S, et al. Fred: a GPU-accelerated fast-Monte Carlo code for rapid treatment plan recalculation in ion beam therapy. Phys Med Biol 62 (2017). 7482–504. doi:10.1088/1361-6560/aa8134

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Winterhalter C, Aitkenhead A, Oxley D, Richardson J, Weber DC, MacKay RI, et al. Pitfalls in the beam modelling process of Monte Carlo calculations for Proton pencil beam scanning. Br J Radiol (2020). 20190919. doi:10.1259/bjr.20190919

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Schneider U, Pedroni E, Lomax A. The calibration of CT Hounsfield units for radiotherapy treatment planning. Phys Med Biol 41 (1996). 111.

PubMed Abstract | CrossRef Full Text | Google Scholar

37.Varian Medical System Inc. Proton algorithm reference guide (Eclipse). Tech. Rep. August, Varian Medical Systems, Inc., Palo Alto (2013).

CrossRef Full Text | Google Scholar

38.ICRU. Report 90: Key data for ionizing-radiation dosimetry: measurement standards and applications. J Int Comm Radiation Units Measure 14 (2016). 1–118. doi:10.1093/jicru/ndw029

CrossRef Full Text | Google Scholar

39. Zhu RX, Poenisch F, Lii M, Sawakuchi GO, Titt U, Bues M, et al. Commissioning dose computation models for spot scanning proton beams in water for a commercially available treatment planning system. Med Phys 40 (2013). 041723. doi:10.1118/1.4798229

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Langner UW, Eley JG, Dong L, Langen K. Comparison of multi-institutional Varian ProBeam pencil beam scanning proton beam commissioning data. J Appl Clin Med Phys 18 (2017). 96–107. doi:10.1002/acm2.12078

CrossRef Full Text | Google Scholar

41. Almhagen E, Boersma DJ, Nyström H, Ahnesjö A. A beam model for focused proton pencil beams. Phys Med 52 (2018). 27–32. doi:10.1016/j.ejmp.2018.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Twiss RQ, Frank NH. Orbital stability in a proton synchrotron. Rev Sci Instrum 20 (1949). 1–17. doi:10.1063/1.1741343

CrossRef Full Text | Google Scholar

43. Bortfeld T. An analytical approximation of the Bragg curve for therapeutic proton beams. Med Phys 24 (1997). 2024. doi:10.1118/1.598116

PubMed Abstract | CrossRef Full Text | Google Scholar

44.[Dataset] Gajewski J. Bragg peak analysis (2017).

CrossRef Full Text | Google Scholar

45. Nelder JA, Mead R. A simplex method for function minimization. Comput J 7 (1965). 308–13. doi:10.1093/comjnl/7.4.308

CrossRef Full Text | Google Scholar

46. Lin L, Ainsley CG, Mertens T, De Wilde O, Talla PT, McDonough JE. A novel technique for measuring the low-dose envelope of pencil-beam scanning spot profiles. Phys Med Biol 58 (2013). doi:10.1088/0031-9155/58/12/N171

CrossRef Full Text | Google Scholar

47. Lin L, Huang S, Kang M, Solberg TD, McDonough JE, Ainsley CG. Technical Note: validation of halo modeling for proton pencil beam spot scanning using a quality assurance test pattern. Med Phys 42 (2015). 5138–43. doi:10.1118/1.4928157

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Low DA, Harms WB, Mutic S, Purdy JA. A technique for the quantitative evaluation of dose distributions. Med Phys 25 (1998). 656–61. doi:10.1118/1.598248

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Biggs S, Jennings M. PyMedPhys python package (2019).

CrossRef Full Text | Google Scholar

50. Albertini F, Casiraghi M, Lorentini S, Rombi B, Lomax AJ. Experimental verification of IMPT treatment plans in an anthropomorphic phantom in the presence of delivery uncertainties. Phys Med Biol 56 (2011). 4415–31. doi:10.1088/0031-9155/56/14/012

PubMed Abstract | CrossRef Full Text | Google Scholar

51.ICRU. Report 83. Prescribing, recording, and reporting photon-beam intensity-modulated radiation therapy (IMRT). J Int Commission Radiation Units Measure 10 (2010). 112.

CrossRef Full Text | Google Scholar

52. Pathak P, Vashisht S. A quantitative analysis of intensity-modulated radiation therapy plans and comparison of homogeneity indices for the treatment of gynecological cancers. J Med Phys 38 (2013). 67–73. doi:10.4103/0971-6203.111309

PubMed Abstract | CrossRef Full Text | Google Scholar