A complete workflow for utilizing Monte Carlo toolkits in clinical cases for a double‐scattering proton therapy system

Abstract The methods described in this paper allow end users to utilize Monte Carlo (MC) toolkits for patient‐specific dose simulation and perform analysis and plan comparisons for double‐scattering proton therapy systems. The authors aim to fill two aspects of this process previously not explicitly published. The first one addresses the modeling of field‐specific components in simulation space. Patient‐specific compensator and aperture models are exported from treatment planning system and converted to STL format using a combination of software tools including Matlab and Autodesk's Netfabb. They are then loaded into the MC geometry for simulation purpose. The second details a method for easily visualizing and comparing simulated doses with the dose calculated from the treatment planning system. This system is established by utilizing the open source software 3D Slicer. The methodology was demonstrated with a two‐field proton treatment plan on the IROC lung phantom. Profiles and two‐dimensional (2D) dose planes through the target isocenter were analyzed using our in‐house software tools. This present workflow and set of codes can be easily adapted by other groups for their clinical practice.


| INTRODUCTION
Monte Carlo techniques are important tools in radiotherapy due to their ability to accurately calculate dose in heterogeneous mediums. 1,2 This is particularly true in proton therapy where there are few commercially available algorithms capable of accurately handling the particle transport through heterogeneities and calculating secondary neutron dose. 3 Monte Carlo methods therefore have widely been used for dose calculation in the field of medical physics. One of the Monte Carlo packages commonly used, particularly in proton therapy, is GEANT4, which has been packaged for proton therapy in the TOPAS platform. 2 The Geant4 simulation toolkit, originally created for the study of high-energy particles, currently has applications in many fields that include nuclear physics, astrophysics, high-energy physics, aerospace, and of course medical physics. 4 TOPAS wraps and extends the toolkit so as to provide an easily accessible tool for end users in the radiation oncology and medical physics fields to simulate proton doses. It operates by simulating particle transportation through a user-defined geometry and calculates the radiation dose delivered. In Our institution operates the MEVION S250 double-scattering proton system. 5,6 With this treatment unit, each field has two unique components: a compensator and an aperture. The compensator is an acrylic object designed to conform the distal edge of each beam to the target while the aperture collimates the edges of the beam to the target. These components are designed in the treatment planning software. To calculate dose to a given patient or phantom, our treatment planning software, Eclipse -Varian Medical Systems, Palo Alto, utilizes a pencil beam convolution algorithm. 7 This algorithm has limitations in its accuracy in highly heterogeneous tissue such as lung. 8 It is therefore often appropriate, for research and clinical purposes, to calculate the dose using an alternative method. 4 In this paper, we present a methodology of using 3D Slicer 8 in conjunction with MATLAB (The MathWorks, Inc., Natick, MA, USA) to bring field-specific components generated by the treatment planning software into the simulation space. We then present in-house codes developed to visualize and compare the calculated doses between the simulations and the actual treatment plan. Previous works have provided the guiding information of using the MC toolkits for creating general geometries for simulation purpose. 9 However, the authors feel there are still gaps to be addressed in order to facilitate clinical practice.
We believe that the present procedure will offer end users a seamless workflow and sets of tools for both clinical and research purposes.

| MATERIALS AND METHODS
Two procedures were developed in the present study. The first provides a detailed workflow for the creation of 3D models for fieldspecific components, which can be imported into the simulation geometry for double-scattering systems. The second describes the steps of visualization and comparison of doses generated by MC toolkits and treatment planning system.
To start the simulation, we have first modeled our proton machine in TOPAS and benchmarked against measured commissioning beam data. 5 Then the apertures and compensators for each of the treatment fields were designed with the treatment planning software, Varian Eclipse, and exported as part of the DICOM RN plan file. TOPAS and Geant4 have the built-in capability to import these structures into simulation geometry once they are in stereolithographic (STL) format. 10 There is therefore a need for a solution that converts each aperture and compensator file to this format. The method to convert from DICOM to STL is shown schematically in Figs. 1 and 2 for compensators and apertures, respectively. To address the gantry rotation in the simulation, we found it more straightforward to rotate the CT data set to match the angle of incidence of the proton beam. This was due to the restriction of TOPAS allowing phase space sources to be translated or rotated in the simulation space. Each of the fields was simulated separately with the corresponding patient-specific components in place. The composite dose from the fields was F I G . 1. Process for loading compensators into TOPAS. F I G . 2. Process for loading apertures into TOPAS.
3. An approximation of geometric projection of the field from the source (right side) to the nominal isocenter. The size of components is given by their projection at isocenter. The approximate position and sizes of the compensator and aperture are given at the green and brown lines respectively. When modeling these components, it is necessary to scale them to their position in the beam line. calculated and displayed via summation of the outputs from two DICOM dose data sets.

2.A | Compensator
The information necessary to define each compensator was expressed though several tags in the RN Plan file. Specifically, the size of the compensator was given by the corner position, number of rows, number of columns, and pixel dimensions, and the thickness was given as a dimensionless 2D array. The method, detailed in the Appendix A, uses this information to define a 3D array which is then converted into "solid" and then saved in STL format.
When inserting the device models into the beam line geometry, the size of the compensator needed to be scaled to their actual size determined through its position in the beam line as shown in Fig. 3.
The dimensions of the compensator in DICOM format were given as projection at the isocenter of the treatment delivery system. To account for this, the size was geometrically scaled to the location of the component in the beam line based on the virtual source axis distance (VSAD) (300A, 030A), the "applicator" or snout position (300A, 030D), and the aperture thickness (300A, 0100), as shown in Eq. (1).

2.B | Aperture
Whereas the compensator is defined as an array of thicknesses, the aperture shape was given as a series of coordinates defining perimeter of the cutout region. These coordinates were reshaped and then scaled with Eq. (2) to its respective position in the beam line. RTDose file. In both cases additional software is required to display the dose and perform quantitative analysis. There are a number of tools which can be helpful for the purpose of evaluating and comparing doses. 11 Among these are isodose lines, gamma analysis, DVH analysis, and plan normalization. In evaluating the dose to a given organ, DVH analysis is critical in assessing the likelihood of given biological endpoints. Additionally, comparing doses calculated using different programs and algorithms can be beneficial. For this, gamma analysis is a commonly used tool. 12,13 The open source software 3D Slicer with the radiotherapy extension provides these tools through the extension "Slicer RT." Doses generated in TOPAS and saved as DICOM RT Dose files need only to be normalized and have the position corrected before they were loaded directly into 3D Slicer.  The "Number of Frames" is the number of slices in the z direction and the "Grid Frame Offset Vector" is a vector containing the distance from the "Image Position" in the z direction for each slice.
The "series description" tag was updated to allow for easy identification of the correct file when loading it. These tags were defined for a given geometry and saved using the "dicomwrite" function with Matlab image processing toolbox version 2011b or later.

| RESULTS
The procedure described above was performed on an example case using the IROC lung phantom. 15 The phantom is composed of tissue, lung, bone, and heart equivalent material, and has a target in the center which moves during treatment delivery. To account for the motion a 4D CT scan was done with our CT Discovery using the Varian RPM system. A treatment was plan with two fields was made using the maximum intensity projection from the 4D scan to delineate the target. The compensators and apertures were converted to STL format and were shown with Netfabb in Figs. 4 and 5, respectively. All simulations were done on our Dell PowerEdge T420 F I G . 8. Profiles through the target in each orthogonal direction. The film, TOPAS, and Eclipse doses are shown with blue triangles, green diamonds, and red squares, respectively. In each profile there is a better agreement between the film and TOPAS than between the film and Eclipse calculations. The film was processed by IROC.
computing server equipped with two Xeon E5-2400 CPUs and 64 GB memory. For each field, 2 × 10 8 initial proton particles were used and the computing time lasts from 2 to 3 h.
Normalized-and filtered-simulated doses and Eclipse-calculated doses were then loaded into 3D Slicer. The dose from each was shown side by side in Fig. 6(a). A three-dimensional (3D) gamma analysis was performed using criteria of 5 mm/7% with a 10% threshold and the gamma passing rate was found to be 90.68%. The passing points (points with a gamma index value between zero and 1) were shown in Fig. 6(b). Two-dimensional gamma analyses were also performed in planes (axial, coronal, and sagittal) through isocenter using an in house Matlab GUI. The gamma passing rates were tabulated in Table 1 and the overlaid isodose lines and dose difference histograms were shown in Fig. 7.
The TOPAS-simulated doses and Eclipse-calculated doses were compared with the film placed through each plane of the target in the IROC phantom. Gamma analysis of the films was performed by IROC with criteria of 7% and 5 mm. IROC uses this criteria to account for the effects of film quenching. 15,16 Results from the gamma analysis were shown in Figs. 6 and 7. Profiles from the film, Eclipse, and TOPAS were shown in Fig. 8. In this particular case, we observe acceptable agreement between Eclipse doses and film measurements as per IROC standards of ≥80% gamma passing rate. TOPAS agreement with the film was better than that for Eclipse as is expected in a highly heterogeneous phantom. In both cases, the accuracy of the dose calculation is affected by the fact that target motion during dose delivery is not accounted for. This limitation of the modeling gives a better agreement between the TPS and TOPAS than between either and the film.

| CONCLUSION
In this paper we detail procedures for modeling treatment field-specific components and performing MC simulation for double-scattering proton therapy system previously unpublished. This procedure also allows for users to easily evaluate and compare the doses between commercial treatment planning software and MC simulations. All Matlab code used in this publication is available freely so other in the community can use these methods. It can be found at https://www.mathworks.com/matlabcentral/fileexchange/68660toolkits-for-monte-carlo-dose-simulation-and-visualization

ACKNOWLEDGMENT
We would like to acknowledge the Imaging and Radiation Oncology Core (IROC) and the MD Anderson Cancer Center in providing the lung phantom and performing analysis of the film and TLD results.

CONF LICT OF I NTEREST
The authors certify that they have no affiliations with or involvement in any organization or entity with any financial interest or nonfinancial interest in the subject matter or materials discussed in this manuscript.

R E F E R E N C E S COMPENSATOR
For a given field (i) in a DICOM plan "info", the DICOM tags for compensator are as follows: Position ¼ (300A,00EA) = info.IonBeamSequence.Item (i).

APPENDIX B
The following appendix provides the main details and explanations about data manipulation performed in Matlab. The full code will be made available on the Mathworks website.

APERTURE
The aperture is defined by its thickness, diameter, and area to be cut out. The DICOM tags for this information are given as follows for a plan named "info". The surface defined by the x, y, and AP variables is used to define the solid shape using triangle faces and vertices. This solid shape is in turn saved in STL format The produced STL file will have points defined across the entire square of the meshgrid. This points need to be cleaned up before loading the aperture into the MC simulation. We found that the easiest way to do this is with the Autodesk program Netfabb, which is free for academic purposes. This model of the aperture can now be saved and loaded into the MC simulation. The aperture file size can be reduced by manipulating the mesh to reduce the number of triangles used to define the object shape.

APPENDIX C
The conversion from CopyID to a DICOM format dose grid is given below. NX ¼ numberofcolumns; NY À numberofrows; NZ ¼ numberofslicess