Automation of Monte Carlo‐based treatment plan verification for proton therapy

Abstract Purpose Independent calculations of proton therapy plans are an important quality control procedure in treatment planning. When using custom Monte Carlo (MC) models of the beamline, deploying the calculations can be laborious, time consuming, and require in‐depth knowledge of the computational environment. We developed an automated framework to remove these barriers and integrate our MC model into the clinical workflow. Materials and Methods The Eclipse Scripting Application Programming Interface was used to initiate the automation process. A series of MATLAB scripts were then used for preprocessing of input data and postprocessing of results. Additional scripts were used to monitor the calculation process and appropriately deploy calculations to an institutional high‐performance computing facility. The automated framework and beamline models were validated against 160 patient specific QA measurements from an ionization chamber array and using a ±3%/3 mm gamma criteria. Results The automation reduced the human‐hours required to initiate and run a calculation to 1–2 min without leaving the treatment planning system environment. Validation comparisons had an average passing rate of 99.4% and were performed at depths ranging from 1 to 15 cm. Conclusion An automated framework for running MC calculations was developed which enables the calculation of dose and linear energy transfer within a clinically relevant workflow and timeline. The models and framework were validated against patient specific QA measurements and exhibited excellent agreement. Before this implementation, execution was prohibitively complex for an untrained individual and its use restricted to a research environment.


| INTRODUCTION
Proton therapy is becoming an increasingly common treatment modality in radiation oncology. 1 As this technology matures, consensus guidelines continue to be developed for many proton therapy centers. Although currently there are clear guidelines on the use of independent dose calculations to verify the dose in 3D and intensity-modulated treatment photon therapy, 2 at the time of writing this manuscript, such guidelines did not exist for proton therapy. Due to lack of independent commercial solutions, proton centers often develop their own in-house Monte Carlo (MC) [3][4][5] or analytical 6  can fail in regions of heterogeneities. 7 When heterogeneities exist in the path of the proton beam, the Bragg peak is degraded and it affects the proton fluence and width of the distal fall-off region of the Bragg peak. 8 Goitein and Sisterson, 9 Urie et al., 10 and Sawakuchi et al. 11 reported that the main cause for degradation of the Bragg peak is multiple Coulomb scattering, which is not taken into account in PBA dose estimations. A recent publication from the Imaging and Radiation Oncology Core Quality Assurance Center, Houston (IROC-H), 7 demonstrated the inaccuracy of PBA in predicting the dose in anthropomorphic lung phantoms used in its remote auditing program. Differences observed between PBA and measured dose distributions were due to density heterogeneities close to the bone-air and tissue-air interfaces.
Unlike PBAs, MC calculations include explicit modeling of particle transport, which increases the accuracy of dose calculations. 4,12 Explicit transport of particles allows MC simulations to perform more accurate dose calculations in regions of heterogeneities. The IROC-H study 7 compared dose calculations in the lung phantom among MC, PBA, and film measurements. Dose differences between PBA and film measurements were as high as 30%, whereas the maximum dose difference between MC and film measurements was 12%.
Accuracy of the MC technique, coupled with the observed differences between clinical dose calculations, suggest that an MC-based second check system is well-suited to identify cases in which limitations of the TPS algorithms are clinically relevant.
Due to concerns of linear energy transfer (LET) correlating with treatment-related complications, 13,14 a prospective evaluation of LET distributions before treatment would be useful in establishing best practices with regard to beam arrangement, spot placement, and spot weighting. 15 LET can be easily scored in MC calculations and integrated into the independent secondary check. 5 MC-based patient-specific dose calculations require the execution of several steps before and after MC simulations. When the process is executed manually and in a stepwise manner, it can be time consuming and prone to human error. Increased computational power and improved access, and even direct integration to TPS data with scripting, have enabled workflow automation in the treatment planning process, 16 including the use of secondary dose checks. [17][18][19] In this study, we develop a framework for an automated MC dose/LET verification system that reduces the person-hours required to perform an MC calculation. We present a workflow that can be implemented at any proton therapy center that uses spot scanning.
Furthermore, we focus on studying dose and LET distributions in pediatric patients. If properly developed, an MC-based system can provide an excellent independent secondary check of the calculated dose and also calculate LET distributions before treatment with minimal workflow disruption. The following sections provide details of each step in the automation process. Figure 1 summarizes the automated steps and the environment in which the steps occur.

2.A | DICOM import and export
The automation process starts when the user exports DICOM files associated with a patient plan from the clinical TPS to the server. By passing DICOM data through a server, the identity of the user (sender) is recorded and used to notify the user of results of the data integrity check via e-mail. During the data integrity check, patient data are verified and any missing data required for the MC simulation are identified. Scripts on the server automatically detect the result of the integrity check and launch the next stage of the process. After the simulation is complete, output data and their relevant metadata are sent back to the server and stored in a database that can later be accessed through a web interface or imported to the dose visualization software of choice (e.g., TPS).

2.C | Monte-Carlo simulation
The TOol for PAarticle Simulations (TOPAS) 20 is an MC tool specially built for proton therapy simulations based on the general purpose MC code GEANT4. 21 TOPAS has been used for all simulations in the automated MC system. The previously developed two-stage MC model for simulating treatment with the PROBEAT-V spot scanning proton therapy system was used as the computational engine for the automated system. 22 In the first step, the Hitachi PROBEAT-V beamline and nozzle were simulated, and phase-space files were generated with TOPAS v2. 22 DICOM RT plan data were used to create a set of input files for TOPAS simulations, with the MATLAB script written in-house. To simulate the beamline the TOPAS input contains the scanning magnet parameters, beamline energy parameters, and spot fluence parameters upstream of the scanning magnets.
Each proton field was simulated by using 2 × 10 8 proton histories. For the given number of histories, the statistical dose uncertainty at the center of a 10 cm spread out Bragg peak was estimated to be 0.6% at one standard deviation. The total number of proton histories was split among 50 parallel jobs to save time. To change the event sequence the initial seed was altered for each of the jobs. At the end of the beam nozzle, a phase-space scorer was used to obtain the position, particle type, energy, and momentum of particles crossing that surface and was saved in binary format.
In the second step, the generated phase-space particles were used in the calculation on the CT-based patient geometry located downstream of the nozzle. TOPAS v3.1.3 was used to simulate particle transport by using patient geometry. Information such as the location of CT images, isocenter, and CT machine-dependent density correction factors were also included in the TOPAS input files. To score physical dose and dose-averaged LET (LET d ) in the CT-based patient geometry, usually split into 512 × 512 × (number of CT slices) voxels, "DoseToMedium" and "ProtonLET" tallies were used respectively. At the end of the simulation, each parallel job generated corresponding dose and LET d files in binary format.

2.D | Postprocess
The system detected the end of the simulation stage by tracking the job ID issued by the HPC facility when the job was originally submitted. After calculations were completed, the postprocess stage was automatically initiated. First, the dose and LET d files created from each parallel simulation or batch, typically 50 files per field, were read into dose and LET d grids. LET d and average dose of a voxel F I G . 1. Automation Workflow -User interaction is limited to the scripting application programming interface in the treatment planning systems. All actions within the staging server and highperformance computing environment are automated. DICOM CT and RS are the imaging and structure files, respectively, while DICOM RN is the plan file were calculated by Eqs. (1) and (2) respectively. Dose i and LET i values represent dose and LET d values of the ith batch.
The conversion of Monte Carlo scored dose, usually tallied as either energy per voxel or energy per voxel per particle, requires a correction factor that scales the scored dose by the number of particles simulated and MU in the plan file. 23

2.E | Integration of software
Because in a generalized scenario the independent calculations will require at least two different applications or software packages (i.e., the TPS and the independent calculation), there needs to be in place a method to get each piece of the automation process to communicate with the previous and subsequent steps. We developed a set of

2.F | Validation
The accuracy of the system and integrity of the data exchange in the automation was validated against patient-specific QA measurements. Validation of the beamline model (i.e., comparison of ranges, spot width, etc.) was previously performed and reported on Ref. [24,25]. Briefly, the system validation and patient-specific QA comparison performed in this work entails measurements conducted for at least two depths per patient, per treatment field. The setup con- of greater than 95% of tested pixels passing.

| RESULTS AND APPLICATIONS
The developed automated MC system can be used for verifying the treatment plan dose as well as for dose-weighted LET calculations.
Exact reductions in human-hours are difficult to quantify due to the dependence on user familiarity with the preautomation workflow.
Before implementation, a familiar user spent 1-2 hr on necessary data preparation and processing. For the unfamiliar user, the process can be prohibitively complex and as a result unusable. After implementation, user interaction was reduced to 2-3 min, all of which were spent within the TPS.
In the following sections, we present validation results and illustrate two patient cases to demonstrate the two applications-dose verification and LET calculations-of the automated system.

3.B | Dose comparison
A pediatric patient with Ewing sarcoma of the infratemporal fossa treated with two lateral beams was chosen to calculate the MC dose. This tumor site was selected because of presence of heterogeneities near the target. The MC dose was calculated using the automated MC system and compared with the TPS dose. Figure 3 shows an axial slice of the patient's CT with an overlay of MC and TPS dose calculations.

3.C | LET d calculations
Although LET distributions in proton radiation therapy have been extensively studied, 13,14 no study has quantitatively correlated LET values, volume effects, or spatial distribution with the likelihood of patients developing radiation necrosis or other radiation-induced toxicities. In the absence of quantitative guidelines on LET distributions, we propose using a qualitative approach to minimize LET-related biological uncertainty. Figure 4 shows sagittal images of a pediatric patient with an atypical teratoid rhabdoid tumor (ATRT) treated with proton therapy.
Images are color wash maps of dose (a and b) and LET d (c and d) distributions for two unique treatment plans (1 and 2) designed to deliver the same physical dose and meet the same clinical constraints.
The average 3D gamma passing rates for a ± 3%/3 mm criteria are reported for each measurement scenario. The gantry, fixed, and fixed mini-beamlines are distinct beamline models. The nozzle mounted range shifter was explicitly modelled in the TOPAS simulations. The only change in plan was the beam angles used to cover the target. In the images 1c and 1d, areas of increased LET d are present in the posterior part of the brainstem. After the beam orientation is changed, images 2c and 2d show a reduction in LET d in the brainstem. This reduction in LET d was achieved without compromising the physical dose distribution traditionally used for the plan evaluation.

Model
We envision a workflow in which treatment plans can be prospectively assessed for biological robustness (i.e., how uncertainties related to LET d maps may affect plan quality) during the treatment planning process. Individual plans can be evaluated, or suitable plans can be compared to evaluate differences in LET distributions.
In our early use of the tool for these types of comparisons we have focused on identifying areas of overlapping high dose and high LET d , defined by us as dose greater than 80% of the prescription dose and LET d greater than 6 keV/micron. When the overlapping these two overlapping areas is in an area of high functional importance (e.g., the brainstem), we will choose beam angles or manipulate spot placement to mitigate end of range effects while preserving dosimetric quality of the plan. As we continue to collect data through the use of the automated Monte Carlo tool, we plan to build a repository of LET d distributions for use in retrospective studies that can improve our understanding of biological risks associated with LET d .

| DISCUSSION
We developed an MC-based The dosimetric comparison between TPS and MC calculations did not address actionable levels at which intervention is recommended. We feel this is best handled on a case-by-case basis through a discussion between the dosimetrist, physicist, and physician involved in the plan. When the discrepancy is due to heterogeneous tissues, the ability to correct for TPS calculation accuracy is limited as that is driving the optimization of spot placement and weighting during the planning process. In these scenarios we typically consider alternate beam paths that might limit the extent to which heterogeneous tissue is traversed.
Alternatively, the renormalization of the dose can be used to ensure minimum coverage or maximum hot spot constraints are met.
The second case (Fig. 4) was that of a pediatric patient with ATRT. The choice of beam delivery angle was evaluated based on resultant LET d maps, without compromising dose distribution. Even though this result was not clinically used, it is a good example to demonstrate that critical organs or sites could be spared by taking LET d distributions into account during treatment planning. The concept of biologically robust planning by using different optimization strategies or treatment angles has been previously explored. 15 Our example demonstrates how biological robustness can be evaluated within the clinical workflow.

| CONCLUSION
The main purpose of developing the automated MC system was to