An Integrated Framework Based on Full Monte Carlo Simulations for Double-Scattering Proton Therapy

Purpose: We developed an integrated framework that employs a full Monte Carlo (MC) model for treatment-plan simulations of a passive double-scattering proton system. Materials and Methods: We have previously validated a virtual machine source model for full MC proton-dose calculations by comparing the percentage of depth-dose curves, spread-out Bragg peaks, and lateral profiles against measured commissioning data. This study further expanded our previous work by developing an integrate framework that facilitates its clinical use. Specifically, we have (1) constructed patient-specific applicator and compensator numerically from the plan data and incorporated them into the beamline, (2) created the patient anatomy from the computed tomography image and established the transformation between patient and machine coordinate systems, and (3) developed a graphical user interface to ease the whole process from importing the treatment plan in the Digital Imaging and Communications in Medicine format to parallelization of the MC calculations. End-to-end tests were performed to validate the functionality, and 3 clinical cases were used to demonstrate clinical utility of the framework. Results: The end-to-end tests demonstrated that the framework functioned correctly for all tested functionality. Comparisons between the treatment planning system calculations and MC results in 3 clinical cases revealed large dose difference up to 17%, especially in the beam penumbra and near the end of beam range. The discrepancy likely originates from a variety of sources, such as the dose algorithms, modeling of the beamline, and the dose metric. The agreement for other regions was acceptable. Conclusion: An integrated framework was developed for full MC simulations of double-scattering proton therapy. It can be a valuable tool for dose verification and plan evaluation.


Introduction
A full Monte Carlo (MC) proton-dose algorithm takes into account the detailed physical processes for both primary protons and generated secondary particles during the entire particle-transport path. It is widely considered more accurate than other dose-calculation methods, such as analytical and pencil-beam algorithms [1][2][3][4][5][6]. However, the MC dose algorithm is not available in current commercial treatment planning systems (TPSs) for passive double-scattering proton therapy. A MC simulation is usually more difficult for the http://theijpt.org double-scattering system than for the pencil-beam scanning system. The full MC dose algorithm for the double-scattering system requires simulation of major mechanical components, such as beam scatterers, range modulation wheels (RMWs), a range shifter assembly, snout, and patient-specific brass collimator and range compensator [7,8]. This requires the detailed design information of geometry and material composition of these mechanical components, which is usually not available from the manufacturer because of proprietary concerns.
In our previous work [8], we overcame that obstacle by developing an innovative approach. We called it virtual machine source model (VMSM), to distinguish it from the virtual analytical source model. Within the VMSM, a virtual machine or beamline is numerically built for all considered mechanical components to represent a real physical beamline so that the simulated results using the virtual machine are closely matched to the measured reference data. This approach is based on 2 observations: (1) the exact MC model is impossible or may not be necessary and, therefore, tuning the physical parameters is not uncommon in full MC modeling, even in a situation in which the design information is provided by the manufacturer; and (2) the acceptability and correctness of the model can be generally approved by comparing the calculation results with measurements. The agreement between model predictions and measurements provides validation of the model, as widely used in the commissioning of photon radiotherapy. We have applied this approach to our passive double-scattering proton therapy machine, the Mevion S250 (Mevion Medical Systems, Littleton, Massachusetts), and finally determined reasonable geometric configuration of the beam collimation components with numerous trial-and-error iterations. Good agreement was achieved between the calculated results and commissioning beam data for the percentage of depth-dose curves, spread-out Bragg peaks (SOBPs), and lateral profiles measured in a water phantom.
The purpose of the present work was to extend our established model for simulation of real patient treatment plans and to develop an integrated framework to facilitate the clinical use of the model. Specifically, we have added the following components into our previous model and system: (1) a patient-specific applicator and compensator were reconstructed numerically from the patient plan data and incorporated into the beamline for full MC dose calculations, (2) a calibration procedure was developed to convert the MC calculation dose in the unit of dose per simulated particle into the absolute dose in units of gray, and (3) an integrated framework was developed to provide a seamless connection between the TPS and the MC simulation, such as generating a patient phantom, a geometrical transformation with the conservation of the isocenter position, and preparing the input for the MC calculation. In addition, an automation process to execute the MC program in a parallel fashion was implemented to speed up the calculation.

Overview of the Framework and Clinical Implementation for Proton Therapy
We have developed a home grown system called VeriPlan, which focuses on convenient tools for dose verification of patient treatment plans for various treatment modalities, including TomoTherapy (Accuray, Sunnyvale, California) [9] and Gamma Knife (Elekta, Stockholm, Sweden) stereotactic radiation surgery [10,11]. The system was written in the Cþþ programming language. A couple of open-source packages are used in the system; for example, the Visualization Toolkit (Kitware, Clifton Park, New York) is used for 3-dimensional (3D) image processing and visualization, and the Grassroots Digital Imaging and Communications in Medicine (DICOM; National Electrical Manufacturers Association, Arlington, Virginal) is used for handling the DICOM image input and output. The system consists of a core framework and an additional plug-in architect that allows for extending the capabilities for external applications. Compared with other comprehensive toolkits, such as 3D slicer [12] and its extension SliceRT [13], the system is more compact, and the emphasis is on providing accurate dose calculations using MC simulations. It is an integration of our previous research work with a friendly graphical user interface with several tools, such as image processing, visualization, and result analysis.
More specifically, the framework has 5 main components: (1) an imaging module that provides utilities for importing and exporting DICOM radiation therapy (RT) data from, and to, a TPS and for image visualization; (2) a source-modeling module that is implemented as a plug-in to provide parameters for various radiation-treatment modalities for MC simulations; (3) an MC dose engine that can be either the internal MC code [22] for photon and electron calculations or for external packages, such as the Geant4 (CERN, Geneva, Switzerland) code; (4) an analysis module that provides tools such as dose comparison and dose-volume histogram (DVH) calculations; and (5) an in-memory database module that provides storage and retrieval of patient treatment-related data.
For current work, the MC simulation for proton therapy was implemented as a plug-in to the framework. For a treatment plan, the computed tomography (CT) image, structures, plan, and dose in the format of DICOM RT are imported into the patient database. The patient image is used to generate the geometric volume, with each pixel represented by a material with the help of a CT number-to-density and material table. The structures and dose distributions are used for DVH calculations and for dose comparisons. The RT plan contains information, such as the gantry angle and the shape of beam aperture and compensator. For each beam in the RT plan, the beamline parameters are determined according to the field size, range (distal 90% dose level of the SOBP), modulation (distance between the distal and proximal 90% dose level of the SOBP), and the position of the snout. The field size and beam range together determine the beam option to be used. The chosen option and the desired modulation together determine the stop angle of the RMW using the information from the dosimetry table. The nozzle position determines the positions of the applicator and compensator. The beam applicator and compensator are represented by voxelized volumes, and their shapes are determined by the milling machine data. Because the Geant4 MC calculation is performed in the machine-coordinate system, the patient volume needs to be transformed according to the gantry and couch angle. With all the necessary input data being prepared, multiple MC runs are executed in a simple, parallel fashion (no dependency between the parallel tasks) by specifying independent seeds for random number generation. All the dose distributions from the subtasks are accumulated and then transformed to the patient coordinate system for the field. The total dose from all treatment fields is finally calculated with the beam weight as the multiplier.
The transport of proton particles was performed with Geant4 code [14] and our virtual beamline for the Mevion S250 compact proton therapy system, which was described in our previous work [8]. The included physics processes are energy loss and straggling in electromagnetic interaction with atomic electrons, multiple scattering because of coulombic interactions with atomic nuclei, elastic scattering, inelastic scattering, and nuclear reaction with the nuclei. The modular physics list emstandard_opt4 is used for electromagnetic interactions, and the physics lists G4DecayPhysics and G4RadioactiveDe-cayPhysics are included for simulation of the decay of subatomic particles and radioactive decay of the nuclei. For hadronic physics, we include G4HadronPhysicsQGSP_BIC_HP and G4HadronElasticPhysicsHP for the physics of inelastic and elastic scattering of hadrons, respectively, with the high precision (HP) model. Ion hadronic interactions are included using the physics list of G4IonBinaryCascadePhysics. For the beamline geometric design, we digitally constructed a virtual doublescattering proton machine [8], which includes all major beam-shaping components that interact with the proton beam along the beam path. Guided by the machine configuration data, the RMWs, first and second scatterers, and the characteristics of the proton source were generated and tuned to reproduce the measurement as closely as possible. The thickness of the wheel steps was derived from the pristine Bragg curves and the full-modulation SOBP. The validation of the RMW design was performed by comparing the measured SOBPs with the MC calculations. In our model, the RMW numerically rotates, and the desired modulation is generated by setting a stop angle and by modulating the beam fluence. The model has been validated by comparing the percentage of depth-dose curves, the SOBPs, and the lateral profiles against measured beam data.

Modeling of Patient-Specific Aperture and Compensator
The milling machine data for the aperture and compensator were extracted from the treatment DICOM RT plan. For the aperture, a 2-dimensional data array describes the shape of the aperture. To create the virtual hardware of the aperture incorporated in the beamline, 2 steps were taken because no geometrical primitives exist in the Geant4 library that can be directly used for modeling the arbitrary curve shape of the aperture. First, a rectangle encompassing the aperture was found using the maximum and minimum value of the 2-dimensional array. This rectangle was then discretized by uniform voxels of 2 mm. The voxels outside the aperture (eg, represented by the white portion in the pink rectangle shown in Figure 1, left) were assigned the material of copper, and the inside of the aperture was assigned air. Second, a composite structure was created Figure 1. Illustration of the modeling of a patient-specific beam aperture (left) and compensator (right). The yellow color in the aperture and the portion of pink color in the inner rectangle is copper. The green color in compensator and the portion of pink color in the inner rectangle is Lucite.
using the Geant4 element G4SubtractionSolid by subtracting a rectangle G4Box from an object representing a copper cylinder (G4Cons). The rectangle had the same size as the rectangle from the first step. By correctly setting the center and position of the rectangle, the combination of the created objects in these 2 steps can produce a virtual aperture in the beamline. The same strategy was used to model the compensator, except that a 3D data array was extracted to describe the shape of the compensator. For example, the compensator shown in Figure 1, right has a dimension of 36, 44, and 13, with a spatial resolution of 2.271 mm, 2.271 mm, and 2 mm. The voxels outside of the cavity of the compensator were assigned the material of Lucite, and the cavity of the compensator was assigned air.

Representation of Patient Geometry
The patient phantom was created from the DICOM CT image exported from the TPS. The CT image contains the Hounsfield unit (HU) information for each voxel. For patient CT images, the dimensions of a slice are usually 512 by 512, with the pixel spacing being about 0.78 mm. The slice thickness is usually from 1 mm to 3 mm. If the original CT resolution is used, it will take a long time for the MC calculation to reach a clinically acceptable level of statistical uncertainty (, 2%). Therefore, in our program, we provide a user option for choosing the calculation resolution by resampling from the original CT. For this work, we used the grid resolution of 2 mm 3 2 mm 3 2 mm for all the cases, the same as in TPS plans.
The patient phantom was transformed to the machine coordinate system to use Geant4 for MC simulations. The resampled CT image is in the patient coordinate system in which the direction of positive x-axis is from the patient's right to left, the direction of positive y-axis is from the patient's anterior to posterior, and the direction of positive z-axis is from patient's inferior to superior. In the design of our virtual Mevion proton machine, however, we define the positive z-direction as from room ceiling to floor, the positive x-direction as toward the gantry, and the positive y-direction as determined by the right-hand rule. The gantry angle and couch angle from a patient treatment plan must be aligned to the patient geometry in the machine. A transform matrix was created to perform the alignment by rotating 908 on the x-axis, the couch angle on the z-axis, and the gantry angle on the y-axis. The calculated dose distribution was then transformed back to the patient coordinate system using the inverse of the transform matrix. Note that careful manipulations are needed to ensure the beam isocenter stays the same in these transforms.
For MC calculations, the material composition, defined by mass density and elemental weights for the patient phantom, is required. We adopted the method suggested by Schneider et al [15] to assign the HU numbers in the CT to air, lung, soft tissue, and bone by dividing the HU number from À1000 to þ1600 into 24 groups. Within each group, the elemental composition and weight were the same; however, a new material was created using a density bin index of (D À D c )/dD, where D is the voxel density, D c is the defined group density in the table, and the bin interval dD is 0.005 g/ cm 3 . Therefore, thousands of materials may be created in an MC calculation. A method to dynamically assign the mass density to the material during particle transport was developed by Jiang et al [16]. In this work, however, we used a simple approach, that is, the material was assigned for each voxel statically before the particle transport started in the calculation.

Dose Calibration
As we know, the dose distribution calculated by the MC method is in units of dose per particle. To obtain an absolute dose, a calibration factor is needed to convert it to a dose in grey units. In our institution, Pinnacle TPS (Philips Healthcare, Amsterdam, the Netherlands) is used for treatment planning. Note that the dose distribution calculated by the TPS is relative. It needs to be rescaled to the prescribed dose at the prescribed target location, which is usually the plan isocenter. To obtain the machine monitor units (MUs) that deliver the desired dose, an output factor is required. At our institution, 103.8 centigray (cGy) was measured with an ion chamber if 100 MUs were delivered under the beam configuration of a range of 15 cm, a modulation of 10 cm, and a snout position of 20 cm. This beam configuration is referred to as the dose-calibration beam configuration, labeled c. The output factor is defined as the ratio of the dose at the center of the SOBP for the measured beam and the calibration beam. To obtain the dose distribution in a patient, 2 MC calculations were performed. One calculation was performed with the same beam settings as that for the patient treatment plan, except no compensator was inserted in the beamline, and the irradiated volume was the water phantom instead of a patient. This beam configuration is referred to as the reference configuration, labeled r. It is the same configuration as in the output factor measurement. Another calculation was performed according to the actual treatment beam configuration, labeled s. With these definitions, the MC-simulated absolute dose and MU can be calculated as follows: where w is the beam weight, D p is the prescribed fraction relative biological effectiveness (RBE) dose, N i is the number of protons simulated, and M i is the MC-calculated dose result for beam configuration i (i ¼ c, s, r). To calculate the MU, the fraction RBE dose D p needs to be divided by the proton RBE factor (1.1) to convert to the fraction physical dose.

Validation of the Framework with End-to-End Tests
The MC simulation for a patient proton-therapy treatment plan is a complicated process, which requires data exchange between the TPS and the MC program, accurate geometrical transformation of the patient phantom in different coordinate systems, and correct descriptions of components in the beam-delivery system. To demonstrate the correct implementation of the framework, end-to-end tests of simple cases in a pure-water phantom are desirable to reduce the uncertainties coming from the conversion of the CT HU to materials.
To make the whole procedure the same as using real-patient CT images, we created a series of pseudo-CTs for the water phantom. The phantom had a volume of 40 cm 3 and a dimension of 200 3 200 3 200 with a resolution of 2 mm 3 2 mm 3 2 mm. This image was imported into the TPS as handled for clinical patients. The treatment plans were generated with specified beam parameters and then exported for MC simulations. Two sample plans were created. In our model, the beam configuration (setting) is fully determined by the field size, range, modulation, and position of the snout from our MC program. To verify the correctness of the algorithm for choosing the right beam option (there are a total of 24 beam options) for a given beam range and modulation, we created a simple plan with only 1 field using the calibration beam configuration, as discussed earlier. No compensator and aperture were used. The MC calculation for this beam setting also served for the absolute dose calibration. The calculated depth dose was compared with the TPS result, as shown in Figure 2. We can see that they are in good agreement. The slight fluctuation of the plateau near the distal edge is due to the composition of the first couple of individual, pristine Bragg peaks of different ranges.
Another end-to-end test was generated in the TPS by creating an artificial planning tumor volume (PTV) in the water phantom, and 1 beam was used to irradiate the PTV as handled in a clinical plan. The shape of the aperture and the compensator were determined by the TPS. This test was to verify the correctness of the geometric transformation and numerical representation of the aperture and compensator. The plan had a couch angle of 1458 and a gantry angle of 308. The beam parameters were range (19.25 cm), modulation (9.41 cm), field diameter (14.89 cm), and snout position (42 cm). For the quantitative comparison of 3D dose distributions, the c-index method was used. For all cases in this work, we considered only the dose points with values . 20% of the prescription dose. We used our previously developed method [17] for the c-index calculation, that is, a k-d tree was first built with the reference dose distribution (TPS calculated), and then, the c-index for the  Figure 3. We can see that the agreement is good, except for the region near the end of range, where the MC result showed statistical noise, with even 10 8 proton particles simulated. The difference was likely caused by differences in algorithms and head modeling used by the TPS and MC method, considering that the same water phantom was used.

Clinical Case Studies
We used 3 cases (spine, head and neck, and lung) to demonstrate the application of the framework for clinical usage. The 3 patients have been treated with our double-scattering Mevion proton therapy system. The resolution of the CTs was 0.78 3 0.78 3 1 mm 3 . In the MC dose calculation, we resampled the CT image slices into a phantom volume with a resolution of 2 3 2 3 2 mm 3 to decrease the computing time and reduce the statistic noise; 10 8 proton particles were generated and simulated. The calculated dose distribution was converted to the RT dose format and imported into SureCalc software (MIM Software, Cleveland, Ohio) for data analysis, such as computing the dose difference and DVHs. The first case was a treatment of left paraspinal mass at the region of spine T12. The prescribed dose was 30 Gy (RBE) with 12 fractions. Two posterior beams were used. The beam parameters were a range (15.48 cm and 17.38 cm), modulation (8.28 cm and 8.72 cm), field diameter (9.67 cm and 9.67 cm), and snout position (26.12 cm and 28.73 cm) for the posterior-anterior beam (gantry angle 1808) and the posterior beam (gantry angle 1508), respectively. The fraction dose was 2.5 Gy (RBE) and prescribed to the same isocenter. The field relative weight was 40% and 60%, respectively. Figure 4 shows the MC-calculated dose distributions for the 2 beams in 1 axial and coronal slice, compared with the dose distributions calculated by the TPS. The dose difference is also given in the figure. We can see that large differences (6 2 Gy, up to 17% of the prescription) exist in the lateral penumbra regions. The effect of the bony anatomy is clearly shown to be the cause of the large differences in the distal fall-off region.
The second example was a case of head and neck treatment. The patient had biphenotypic sinonasal sarcoma. Three fields were delivered: 1 anterior beam and 2 posterior beams. The prescribed dose was 50 Gy (RBE) in 25 fractions. The beam parameters for the 3 fields were range (11.95 cm, 16.09 cm, and 18.56 cm), modulation (10.19 cm, 7.85 cm, and 10.01 cm), field diameter (10.3 cm, 12.97 cm, and 11.18 cm), snout position (19.7 cm, 27.55 cm, and 29.04 cm), couch angle (2708, 1808, and 08), and gantry angle (3558, 1358, and 1358), respectively. The fraction dose was 2 Gy (RBE). The field relative weight was 30%, 37%, and 33%, respectively. Figure 5 shows a comparison of the total dose distribution calculated by the TPS and the MC. Note that the dose outside the patient is shown in the MC calculation. Differences can be clearly seen at the region of nasal cavity. Figure 5, third row shows the DVHs for the PTV and other organs at risk. The agreement is reasonable. The large difference for optical nerve, lens, and optical chiasm may be due to the small volume of these structures because small changes in the dose distribution may have a large influence for these structures in the DVH calculations. The lesser dose coverage of the PTV mainly came from the nasal cavity, which was also included in the PTV.
The third case was a patient with stage 3 non-small cell lung cancer. Two anterior beams with the same gantry angle (458) but different couch angles (2008 and 1608) were used. The prescribed dose was 64 Gy (RBE) in 32 fractions. The beam parameters for the 2 fields were range (10.03 cm and 10.77 cm), modulation (8.41 cm and 9.49 cm), field diameter (11.15 cm and 9.33 cm), and snout position (19.28 cm and 21.12 cm). The fraction dose was 2 Gy (RBE). The field relative weight was 50% and 50%. Figure 6 shows a comparison of the total dose distribution calculated by the TPS and the MC. Again, the dose outside the patient is shown in the MC calculation, which can be trimmed out in future work. Differences can be seen in the beam penumbra and near the end of range. Figure 6, bottom shows the DVHs for PTV, carina, right lung, heart, and spinal cord. The agreement is quite good, except for the cord. The difference apparently came from the lung-dose calculation and shows at the region of the distal range within the vertebrae, which is the area that the pencil-beam algorithm cannot calculate accurately because of range degradation.

Discussion and Conclusions
Because proton therapy is increasingly used to treat many anatomic sites, investigation on the accuracy of the treatment planning algorithm is of importance. Currently, most commercial TPSs use a simple analytic pencil-beam algorithm for dose calculations. Several research groups [6,23] have demonstrated that the dose calculated by MC algorithms is more accurate than the dose calculated by an analytic algorithm comparing measurements, especially for lung treatments. The difference can be as large as 46% within the planning target volume in a lung phantom measurement, and the average difference between the pencil beam algorithm and the measurements was 8% from 1 study [6]. Implementation of an MC-based proton-dose algorithm is highly recommended for thoracic applications; otherwise, the efficacy of proton therapy could be adversely compromised. In the present work, we made a great effort toward this goal by developing an integrated framework to facilitate the clinical use of the MC dose algorithm for proton therapy. The tool we developed, although time consuming (with the The MC modeling for a double-scattering proton system is more complex than a pencil-beam scanning system [18,19]. In a pencil-beam scanning system, a phase-space file (PSF) can be generated by tracking particles through the treatment head by recording particle type, energy, position, and direction on a plane at the exit of the treatment head. Particles can be then sampled from that PSF to track though the patient geometry. That PSF can be reused for different patients to save computation time. Unfortunately, that approach is not possible for a double-scattering proton system because the treatment head configuration is highly dependent on the patient field. There is no common treatment head setup that can be used to generate a PSF and be used by other beam setups. In the Mevion system, there are 24 beam configurations from different combinations of the first scatterer, RMWs, and upstream and downstream absorbers. Even when the same RMW is used, the segment of the RMW in the beam path, which is determined by the modulation of the treatment field, will be different. In addition, the patient-field-specific aperture and compensator make the PSF approach impractical. This served as motivation for us to develop a full MC model based on a VMSM for the double-scattering system.
To facilitate clinical usage of the model, we have developed an integrated framework, including a friendly user interface. We presented calculation results for 3 treatment plans (7 fields) to demonstrate the feasibility of its application. Differences between the TPS calculations and the MC results were revealed, especially in the beam penumbra and near the end of range. The discrepancy likely originates from a variety of sources [20]. First, the dose algorithm used in the TPS is the pencil-beam algorithm [2], in which the inhomogeneity correction is handled by stretching the precomputed dose-deposition kernel (calculated in water) based on the changes in density and the superposition of multiple pencil beams. Because of approximations in the pencil-beam algorithm, it may not provide sufficient accuracy in heterogeneous media. Second, the reported dose in the TPS is dose-to-water because, in a pencil-beam algorithm, human tissues are represented as only water with various densities, whereas MC simulations calculate dose-to-tissue because the MC calculations are based on the exact material composition. It is not straightforward converting the dose-to-tissue in an MC to dose to-water [19,21] in the TPS for comparisons for proton beams because, in addition to ionization, nuclear interactions and multiple Coulomb scattering states also contribute to the energy deposition. A crude approximation for the conversion is to use the ratio of the mass-stopping powers. In our work, we used another approach, that is, a patient phantom was created by assigning only water materials but with various densities, according to the HU number. The purpose of this additional calculation was to investigate whether the result, based on this artificial phantom, agreed better with TPS calculation. The comparison (not included in the article) did not show a significant improvement in agreement. Therefore, we have only reported the MC dose-to-tissue results in this work. As pointed out by Paganetti et al [19], reporting dose to tissue may become more favorable with an increasing use of MC dose calculations. The other obvious source for the discrepancy comes from the modeling of the treatment head. In our full MC simulation, all major beam-shaping components were included, and the model simulates interactions of the beam with rotating RMWs. The beam configurations, such as the beam option, the stop angle of the RMW, and the positions of the applicator and compensator, are determined by field size, range, modulation, and the position of the snout. In Pinnacle TPS, however, a model and measurement-based algorithm is used, and the source modeling is based on parameters such as effective sourceto-axis distance, virtual source-to-axis distance, and effective source size, which are used to model the intensity losses of the proton beam from divergence, the divergence of the beam profile, and the geometrical shape of proton source, respectively. With all these factors mingled together and without a benchmark, it is difficult to make a definite quantitative comparison between the TPS and MC calculations. In future work, we will perform experiments in an anthropomorphic phantom to further evaluate the accuracy of our MC system.
Our framework includes the process, such as importing the treatment plan in the DICOM format, constructing the patientspecific beam compensators and apertures, and accumulating the calculation results to obtain the final dose distribution. It has the potential to be used for quality assurance and dose verification for our proton therapy treatment plans. It also has the potential to be developed into a TPS system. However, although the framework is capable of executing multiple MC runs in a parallel fashion, the calculation is still time consuming. For example, for each field, 10 8 proton particles were simulated, and it took about 3 hours using 20 computing nodes (approximately a total of 60 central processing unit hours). This time may be acceptable for quality assurance and dose verification, but it is too slow for treatment planning. Future work will include investigation of the runtime optimization and techniques of variance reduction to reduce the calculation time to a more acceptable level for clinical applications.

ADDITIONAL INFORMATION AND DECLARATIONS
Conflicts of Interest: The authors have no conflicts of interest to disclose.