ProTheRaMon—a GATE simulation framework for proton therapy range monitoring using PET imaging

The main advantage of proton radiotherapy, compared to conventional radiotherapy using x-rays or electrons, is the ability to precisely deliver the therapeutic dose to the tumor, while sparing surrounding healthy tissues in direct vicinity of the target volume. This characteristic property is due to the depth dose distribution of proton beams, which deposit the maximum dose at the end of their range (in the Bragg peak; BP), followed by a steep dose fall-off (Wilson 1946). Although the steep dose gradient could be beneficial for sparing critical structures located distally to the tumor, it is not fully exploited in clinical practice because of the uncertainty of the BP position in depth, commonly named range uncertainty (Knopf and Lomax 2013). The range uncertainty may occur due to uncertainties in CT calibration, imaging artifacts and variations in patient anatomy during the course of the fractionated treatment, and may lead to differences between the planned dose distribution and the dose actually administered. Therefore, a major subject in proton radiotherapy research is the development of detector techniques for proton range monitoring, which would allow treatment planning and delivery to fully benefit from the physical properties of therapeutic proton beams (Parodi and Enghardt 2000, Shakirin et al 2011, Knopf and Lomax 2013, Parodi 2020).

The most widely explored approaches for range monitoring are based on measurements of secondary particles produced in nuclear reactions of the proton beam with nuclei in the patient's body. One approach concentrates on measurements of prompt gammas emitted from the excited nuclei (Kurosawa et al 2012, Krimmer et al 2018). Another focuses on the measurements of 511 keV electron-positron annihilation gammas, where the positrons are the decay product of β⁺-emitting isotopes generated by the proton beam during irradiation. These 511 keV gammas can be measured with positron emission tomography (PET) scanners (Hishikawa et al 2002, Enghardt et al 2004, Parodi et al 2005, 2007b, Nishio et al 2010). Different approaches to PET range monitoring have been investigated to overcome the limitation of incorporating the scanner in the treatment room (Kraan 2015). Commercially available diagnostic PET systems installed outside the treatment room were used for imaging patients offline, however the obtained PET images suffered low signal and biological washout (Bauer et al 2013), which limited their clinical use for range monitoring. Early developments of on-line PET range monitoring systems date to 1990s at Gesellschaft für Schwerionenforschung GSI, Darmstadt, Germany (Jäkel et al 2022) and in 2000s at the proton therapy facility of the National Cancer Center, Kashiwa, Japan (Nishio et al 2010). The application of the commercial PET scanners tested at Heidelberg Ion Beam Therapy Center (HIT), Germany, and Massachusets General Hospital (MGH), Boston, USA (Parodi et al 2008) showed insufficient efficiency of the off-line range monitoring. The mid 2010s saw the development of new era in situ PET systems capable of acquiring treatment verification images during and/or just after the irradiation (Ferrero et al 2018).

Monte Carlo (MC) simulations are a useful clinical and research tool in medical physics, enabling study of the feasibility of new treatment and imaging techniques, often supplementing complex and demanding measurements (Muraro et al 2020). In proton therapy, MC techniques are extensively used to optimize the design of new detectors and for preliminary characterization of their performance (Kraan 2015). The design and optimization of new PET scanner systems for range monitoring is particularly challenging. It requires the simulation of proton beam delivery and production of β⁺-emitting isotopes in the patient geometry, taking into account the isotope decay during patient handling in the treatment room, the simulation of 511 keV gamma transport to the PET detectors of particular design, the coincidence events scoring and PET image reconstruction. The full MC simulation in a one-step approach is not optimal for this purpose because it requires time-consuming tracking of the history of several particles and their interactions and a complex data analysis (Muraro et al 2020). Further challenges relate to the implementation of the software environment and the time performance of simulations. The simulation of the production of β⁺-emitters requires relatively large statistics of primary particles compared to simulations of dose or linear energy transfer distributions. Moreover, the effectiveness of range monitoring depends not only on the detector system performance but is also specific to treatment indication, tumor volume and location. Therefore, new PET designs need to be optimized and verified by extensive simulation studies with multiple treatment scenarios and patient treatment plans. Given the high complexity of the proton treatment and PET imaging simulations, there is need for a simulation framework which is flexible, automated, and has good time performance.

In this paper, we present a simulation framework, PROton THErapy RAnge MONitoring—ProTheRaMon, based on the Geant4/GATE Monte Carlo toolkit, the Customizable and Advanced Software for Tomographic Reconstruction (CASToR) PET image reconstruction toolkit and in-house implemented scripts. The ProTheRaMon framework was designed with the aim of characterizing and optimizing the performance of a PET scanner for proton therapy range monitoring. Our goal was to implement an automated and time-efficient simulation platform to compute the expected activity of β⁺ emitters in phantoms and patients undergoing proton therapy and to predict PET images that may be acquired after the irradiation of separate fields or after the whole treatment.

The requirement of time-efficient simulations has been addressed here by a multi-stage simulation protocol and its deployment on a computational cluster. In the methods section of this paper, we present the implementation and functionality of the framework. In the results section, we show an example application of the framework to the simulation of proton treatment and computation of the β⁺ production in a patient treated at the Cyclotron Centre Bronowice (CCB) in Kraków, Poland, as well as PET acquisition and image reconstruction using the J-PET, a novel multi-photon portable J-PET scanner based on plastic scintillator (Rucinski et al 2018, Jakub et al 2019, Moskal et al 2021a, 2021b).

So far, every institute who wants to perform PET simulations has to set up their own simulation workflow. ProTheRaMon is available to the community, and as it is modular, it can easily be adapted and expanded for the specific institution.

The ProTheRaMon framework addresses the challenge of combining software of different types and data in various formats. The clinical data include a CT scan of the phantom or patient, description of clinical structures and a treatment plan consisting of one or more irradiation fields. Accurate simulations of dose and activity produced in a phantom or a patient require a dosimetrically validated proton beam model from a specific proton facility that needs to be provided. Due to the time gap between the end of the irradiation and beginning of PET image acquisition, the activity decay is modeled for each β⁺-emitting radioisotope individually using in-house developed scripts. The simulation of the emission of the 511 keV annihilation photons towards the PET detector is followed by scoring of coincidence events in a predefined PET scanner geometry. The PET scanner geometry definition is also used during image reconstruction.

For the validity of simulations, it is important to implement a geometrical transformation between the coordinate systems specific for the irradiation and for PET scan acquisition. See figure 1 for details on relationship between those coordinate systems.

Zoom In Zoom Out Reset image size

For the simulation of the proton treatment and PET imaging, the following software packages have been used. The ProTheRaMon framework combines the GATE software package (Jan et al 2004, 2011), which is based on the general purpose Geant4 Monte Carlo simulation platform (Agostinelli et al 2003, Allison et al 2016), CASToR image reconstruction toolkit (Merlin et al 2018), GPU-accelerated Monte Carlo code FRED (Schiavi et al 2017, Gajewski et al 2021), and several in-house implemented Python and bash scripts enabling parallelization of the tasks on a computational cluster, conversion between data formats, data processing, analysis and visualization. We have created, tested and evaluated our framework using GATE version 9.0 (with Geant4 10.7 patch 1), CASToR version 3.1, FRED version 3.60, FREDtools version 0.6.54. The software framework design allows its easy adaptation to a specific proton therapy facility and PET scanner geometry as well as the needs of any research group conducting simulation-based investigations into using PET for proton therapy range monitoring approaches and development of PET detector technologies.

For MC simulations, the GATE software was selected due to its unique capability of combining proton therapy simulations in the patient, the scoring of the production of β⁺ emitters, β⁺ decay and annihilation, the propagation of the 511 keV photons, as well as detection by the PET scanner and processing into coincidence data stored in root files. For the PET image reconstruction, among several available image reconstruction software tools, e.g. STIR (Thielemans et al 2012), PRESTO (Scheins et al 2011), OMEGA (Wettenhovi et al 2021) and others, we have selected CASToR (Merlin et al 2018). The choice of CASToR is motivated by its design intrinsically compatible with the GATE output data format and the ability to import the PET geometry directly from GATE macros. The advanced PET image reconstruction options offered by CASToR include time-of-flight (TOF) processing as well as handling of multi-layer and non-cylindrical PET scanner designs often used for PET-based range monitoring, which are additional requirements for the PET detector used in our investigations. CASToR offers tools for automatic conversion of the GATE geometry macro file to the CASToR format when using the cylindricalPET system. CASToR is very flexible, users could potentially create any PET scanner geometry and describe it in a *.lut file which does not create any limitations for GATE simulations. In that case, the user would have to create a proper converter of the coincidences saved in ROOT format to the CASToR compatible format.

2.1. Design and implementation of the ProTheRaMon framework

In this section we present the concept, design and implementation of the ProTheRaMon framework. The framework is executed in the following five stages that are also illustrated in figure 2:

• Stage 1 : Data pre-processing.
• Stage 2 : GATE Monte Carlo simulation of β⁺-emitting isotope production.
• Stage 3 : Activity map modeling.
• Stage 4 : GATE Monte Carlo simulation of β⁺ emission and PET acquisition.
• Stage 5 : CASToR image reconstruction.

Zoom In Zoom Out Reset image size

The motivation for the design of the framework, consisting of separated data processing, simulation and reconstruction stages is the optimization of the simulation time (Stage 2 and 4), and flexibility to study simulation and reconstruction parameters (Stage 2, 4 and 5). Dividing the framework execution into stages allows also for time efficient, parallelised use of available computing resources. Moreover, the full simulations of primary and secondary particles are overall more time consuming with respect to simulations of a limited number of physics interactions, transport of given particle type, and scoring of selected quantities. Therefore, the proposed framework has a modular construction where each stage has a separate settings, for example a separate physics list in MC simulations to optimize the processing time.

For instance, the modification of the PET scanner geometry requires repeating of the Stage 4 and 5 only, while modification of the reconstruction parameters only Stage 5. Also, Stage 2 can be replaced by other tools or user defined isotope distributions, when the β⁺-emitters production map is obtained with a Monte Carlo simulation code different than GATE, and, once data format compatibility is warranted, the new map can be used as an input of Stage 4.

2.1.1. Stage 1: data pre-processing

The goal of Stage 1 is to prepare all the necessary input data for the subsequent stages. The data required for the simulations in Stage 2 are patient/phantom geometry, treatment plan as well as facility-specific beam model data, exported from the treatments planning system (TPS) and stored in DICOM format (Pianykh 2010). The phantom/patient files include: CT scan of the phantom/patient (DICOMCT), treatment plan (DICOM RT) and structures used during treatment planning (DICOM RS). Furthermore, the proton beam model and CT calibration must be converted from the clinical TPS format to the GATE-specific format and validated as, for example, was performed in previous work (Grevillot et al 2011, Winterhalter et al 2018, 2020a, Gajewski et al 2021).

The CT image, converted from the DICOM CT images, and provided to the MC is required by GATE in one of the acceptable formats, from which we have chosen the MetaImage format (extension *.mhd or *.mha) (Johnson et al 2015). Also, it is good practice to reduce the CT image boundaries to the object of interest only (phantom or patient), by cropping the unnecessary space covered by air in the CT image. This allows greater flexibility in detector positioning while avoiding overlapping volumes of the phantom (or patient) and the scanner in GATE. For the image cropping, we used the external structure saved in a DICOM RS file, which contains the patient body and immobilizing elements. Moreover, to perform the simulations more efficiently, it is recommended to resample the CT image to the reconstucted PET resolution. In our study we used a 2.5 × 2.5 × 2.5 mm³ voxel size. This voxel size is conditioned by the simulation statistics and PET scanner properties.

The DICOM RT file contains all the information concerning the treatment plan, including definition of one or more irradiation fields, irradiation direction, gantry and couch rotations, isocenter position, etc. A dedicated Python script getPlan reads the required data from the DICOM file and converts it into an ASCII macro file in a format that can be read by GATE. Additionally, another ASCII file (CSV format) is prepared containing the information about each field of the treatment plan and consisting of the following columns: DeliveryNo—(integer) the order of delivery of the fields; RS—(boolean) determine the usage of the range shifter in the field; GantryAngle_deg—(floating-point number) angle of the gantry rotation for the field; CouchAngle_deg—(floating-point number) angle of the couch rotation for the field; PencilBeams—(integer) number of individual pencil beams in the field; ProtonNo—(integer) number of protons to be delivered in the field.

All the necessary conversions and manipulations of the DICOM information are performed using the FREDtools package described in section 2.4. Some exemplary files are available at the project github webpage ProTheRaMon (Borys et al 2022).

Additional files required for the simulation are the beam model and the CT calibration. The beam model describes energy-dependent parameters of pencil beams in terms of polynomials. The parameters required are: mean energy and energy spread, parameters describing lateral propagation, separately for the X and Y directions, and a scaling factor converting the monitor units (MU) retrieved from the DICOM RT plan to the number of protons. In our study we defined the polynominals for the entire available nominal energy range from 70 to 220 MeV. Additionally, we have implemented a range shifter (RS) used clinically, which is used in the simulations if requested in the DICOM RT plan. In our study, the beam model was prepared and validated for two gantry rooms available at CCB by converting the experimentally validated beam model implemented in our previous work for the FRED MC code (Gajewski et al 2021) to the ProTheRaMon framework. The CT calibration files contain the information about composition and density for each Hounsfield unit (HU) present in the CT image. This allows GATE to convert each voxel of the CT image to specific material. The facility-specific clinical CT calibration curve was prepared by means of stoichiometric calibration proposed by Schneider et al (1996). The CT calibration curve contains information on the composition and density of 421 materials.

The description of the PET scanner geometry must be provided in the form of a macro file in a format specific to GATE. This file will be used in Stage 4, for PET acquisition simulation, and in the Stage 5 for definition of the reconstruction geometry specific to CASToR.

2.1.2. Stage 2: GATE simulation of β⁺-emitting isotope production

The main objective of Stage 2 is to perform the Monte Carlo simulation from which we obtain production maps of β⁺-emitting radioisotopes. In our framework, we score production maps for 7 isotopes with the longest half-life and the largest contribution to positron emission using the ProductionAndStopping actor, which are listed in the table 1. The generated isotope production maps are used in further steps for PET imaging simulation.

Table 1. Isotopes scored in the form of production maps.

isotope	half-time [s]
15O	122.2
14O	70.6
13N	597.9
11C	1223.4
10C	19.2
30P	149.9
38K	458.2

The set of required input files for this step includes: CT image in MetaImage format and field statistics in CSV format, as well as the treatment plan, CT calibration, beam model, all in GATE macro format. The simulation results of β⁺-emitter production maps are stored in ASCII files, with a simplified header containing basic information about the matrix size, resolution, voxel size and the number of voxels. The production maps have the same resolution as the input CT and output reconstructed PET image, which was 2.5 × 2.5 × 2.5 mm³ in our study, corresponding to the reconstructed resolution of the J-PET scanner (Moskal et al 2021b).

Simulations in Stage 2 exploit the QGSP_BIC_HP_EMY physics list, that includes the Binary Cascade (QGSP_BIC) to model hadronic interactions, high precision library of neutron interactions (HP) below 20 MeV and electromagnetic physics component in option 3 (EMY ) (Zacharatou Jarlskog and Paganetti 2008, Winterhalter et al 2020b). Use of the QGSP_BIC model does not directly calculate the cross sections, but instead samples the possible products of each nuclear interaction from the binary cascade. It is possible to instead use user defined cross sections, e.g. experimentally measured cross sections, through use of look-up tables as applied by Parodi et al (2007a), and McNamara et al (2022). Additional limiters and cuts are used to speed up the calculations while not affecting β⁺ production accuracy in a patient or phantom. Because the world region, in which all the simulation is 'immersed', is filled with vacuum, the step limiter and production cuts for electrons, positron and gammas are set to 1 cm. The particles in the patient or phantom CT volumes are tracked with the step limiter (10 mm), larger than the voxel size, hence the maximum step size was limited by the voxel size and physics interactions. Because the simulations are devoted to production of β⁺-emitting isotopes, the production of secondary electrons, positrons and gammas is irrelevant and their production cuts are set to large value (5 m in our study). On the other hand, the production cut of secondary protons is set to 10 μm. Additionally, in order to simulate properly a pencil beam scattering on the RS, the step limiter and the production cuts were set to 1 mm in the RS region.

Because of the relatively long time required to simulate the production maps, it is good practice to parallelize the calculations and use a computing cluster. In our study, all GATE simulations were performed on the Ziemowit cluster located at the Laboratory of Bioinformatics and Computational Biology at the Silesian University of Technology in Gliwice, Poland.

The GATE package includes tools, namely the GATE job splitter (gjs), to split a simulation into separate processes. Using the previously prepared macro files, files for individual processes are generated. It is important to notice that the gjs script does not change the number of protons in each sub-simulation but just divides the simulation time by the number of threads or machines. Therefore, it is necessary to know the number of threads or machines used for the simulation when generating the master script by the getPlan, before the simulation is split into many sub-simulations, in order to achieve the desired total number of primary protons.

Given the relatively long simulation time and limited available computational resources, it is possible to simulate only a fraction of the total number of protons defined for a field and stored in the field statistics CSV file. In our study, we chose to simulate 10% of protons after determining that this level of statistics resulted in acceptable variance of the isotope production in simulated production maps, while maintaining an acceptable simulation time. This approach requires multiplying the activity values during the simulation of β⁺ emission by an appropriate value—in our case 10.

Within this work, Monte Carlo simulation of the production of β⁺ emitters has been performed using GATE. However, the modular workflow allows other Monte Carlo codes and scoring techniques to be implemented, such as FLUKA or TOPAS (Augusto et al 2018, Onecha et al 2022). The GPU accelerated Monte Carlo code FRED (Gajewski et al 2021) included the implementation of the scoring of the production of isotopes within a patient during proton therapy (McNamara et al 2022). The work validated the predicted isotope production scored in FRED against the isotope productions predicted in this work by GATE. No statistically significant deviations were found, providing confidence that the isotope prediction from FRED may be reliably used. The GPU accelerated implementation allows for simulation of the isotope production maps for the isotopes in table 1 within a fraction of the time necessary for the GATE simulations, with an average time for simulation of a field being 2.4 min. The integration of FRED into the ProTheRaMon workflow reduces the computational requirements of future simulation studies. However, the development of this was concurrent to the development of the ProTheRaMon framework.

Additional to the production maps of the β⁺-emitting radioisotopes, it is possible to score dose and LET maps at this stage, depending on the settings for the getPlan script, with the appropriate physics settings adjusted to these needs.

2.1.3. Stage 3: activity map modeling

The aim of this stage is to aggregate and post-process the production maps generated in the previous step into a single activity map, taking into account the imaging protocol.

The imaging protocol takes into consideration the times necessary for the irradiation and patient setup in the PET detector. Depending on the specific needs of the user, different time structures may be implemented for the treatment as well as in-beam, in-room or off-line imaging protocols. For the purpose of our studies, we chose 90 s delay before the start of irradiation of each field, except the first one. This estimates the duration of exposure (plan delivery) as taking 60 and 30 s is estimated as the time necessary to set up the gantry and beam for the next field. After the last field, the imaging was done immediately, without any additional delay. The time structure was chosen to reflect a clinically relevant protocol and is necessary to take into account the half-lives of individual isotopes when final production maps are generated. More realistic beam times, using log files from the real treatment planning, can easily be integrated if needed.

All activity maps, produced by each sub-simulation in Stage 2, are aggregated for each isotope separately and saved as single files in an Interfile format (extension *.h33/*.i33), ready to be used in the next stage of the workflow.

2.1.4. Stage 4: GATE Monte Carlo simulation of β⁺ emission and PET acquisition

The aim of Stage 4 is to perform simulations of the β⁺ emission in a patient or phantom, propagation of the 511 keV annihilation photons and signal acquisition in a predefined PET geometry.

The input data needed to prepare the macro files for GATE are: a CT file in MetaImage format, a treatment plan in DICOM RT format, with the location of the isocenter, and the activity maps for the seven isotopes listed in section 2.1.2 in the Interfile format. Also, the description of the PET scanner geometry, in the GATE macro file format, has to be prepared. Based on this data, a Python script called petPlan generates GATE macro files needed for simulations. As in the (isotope production simulation) Stage 2, this is the basis for using the GATE gjs script to split the simulation into multiple, shorter, simulations that can be executed in parallel on a multi-core computer or computer cluster. To allow different scanner geometries to be tested, the geometry name is one of the script parameters.

Simulations in Stage 4 are performed using the emlivermore_polar list of physics processes (Kowalski et al 2018) with all the default production cuts, but the user can substitute any other physics model. The list includes the interactions of electrons and photons with matter down to 10 eV, interpolated based on the Livermore library, which is sufficient for the PET response investigation. Physical processes include the photoelectric effect, Compton scattering, gamma conversion, Rayleigh scattering, ionisation and bremsstrahlung.

The simulation results are saved at this stage in ROOT (extension *.root) files (Antcheva et al 2009), where information about the individual events occurring in the detector's sensitive volume is stored. The coincidence data are saved for image reconstruction and, optionally, singles and hits may be stored for further analysis.

It must be kept in mind that the voxelised activity source has a different origin, with the image corner placed starting at [0, 0, 0] mm (x, y, z), than the voxelised phantom image (CT or phantom), which are placed in the center of their geometrical center. In order to align the activity image with a phantom volume, a proper translation is needed. Additionally, to keep the planning target volume (PTV) center (isocenter) in the center of the field of view (FOV) of the scanner, we perform an extra translation of the voxelised source in the z axis.

2.1.5. Stage 5: image reconstruction

The last stage of the workflow collects the data from the second simulation step and aims to produce the final PET images. For the image reconstruction step, the CASToR framework (version 3.1) (Merlin et al 2018) was adapted. However, any software can be used here, depending on the needs. We have used a modified version of CASToR that allows to produce sensitivity maps that include the effect of patient (or phantom) attenuation. This non-standard, for CASToR, combination of sensitivity map with attenuation map in one set is due to the non-standard (multi-layer and non-cylindrical) geometry of the scanner used in our research. In CASToR by default this information is separated and delivered as separate files to the reconstructor.

The geometrical sensitivity is specific for a given PET scanner, while the effect of attenuation is computed for an individual phantom/patient based on the information from a CT scan. CT images are converted to the attenuation map with the bi-linear function described in Carney et al (2006). The sensitivity map was calculated by performing a back-projection of counts from simulated data, weighted by the attenuation of each line of response (LOR) due to the patient. The simulations for the sensitivity map are based on the acquisition of a cylindrical air phantom uniformly distributed with the β⁺ activity (back-to-back sources) covering the whole field-of-view (FOV), resulting in high-statistics list-mode data.

The input files for this stage are: the ROOT files with the coincidence data, the attenuation map stored in the Interfile standard, and the geometry description file (GATE macro). One of the reasons to choose the CASToR image reconstruction software was the fact that it can easily use, with a simple conversion, the GATE geometry description macro file. The first two preprocessing steps at this stage are: the GATE geometry to CASToR format conversion and the simulation to generate the coincidence data from a homogeneous source for each PET scanner geometry used.

The CASToR package uses the maximum likelihood expectation maximization (ML-EM) iterative reconstruction algorithm and can include information about time-of-flight (TOF), which can be set up at the ROOT to CASToR conversion phase.

The reconstructed images are saved separately for each reconstruction iteration in Interfile format, without any information about the coordinate system, but the resulting PET images are in the frame of reference, voxel size and image size as the patient/phantom CT volume.

2.2. Implementation of the framework on a computational cluster

2.3. Example application of the framework for CCB Krakow proton facility and J-PET scanner

2.4. FREDtools for data analysis and visualization

In this section we will present the most important results from all stages of the ProTheRaMon framework. Stage 1 focuses mainly on the preparation of input data for the other stages. It results in appropriate CT image format, appropriate treatment plan format for use with GATE, beam model, CT calibration and scanner geometry description. These are all key data for the next Stages. However, they are not suitable for presentation in graphical form or any kind of summary or plot drawing, therefore the presentation of results starts from Stage 2.

Stage 2 output. The main output of interest from Stage 2 are the production maps, but other output, depending on the settings, are available, such as linear energy transfer (LET) and dose. An example of the β⁺-emitters production accompanied with the dose (left) and with modeled activity (right) maps is shown in figure 3. The dose map can be considered here as a reference distribution. The two next figures give more detailed information about the β⁺-emitter production. Each field that is planned for the therapy is simulated separately and the cumulative map (for all isotopes) for each field is shown in figure 4 (top row). Even more detailed information, presented for the last field, is shown in figure 5 for each isotope separately. Under each of the subplots the colorbar presents the level of the intensity in terms of radionuclei per voxel, representing the amount of radionuclei for each field or for each isotope produced and scored during the simulation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On all images presented here and further, the PTV and Brain Stem regions are displayed, the second represents an organ at risk.

Simulations of the presented in the paper patient, that has four fields in his treatment plan, took an average of 160 min (209; 130; 118 and 181 [min] each field) using newer nodes with Haswell 400 CPU cores.

Stage 3 output. The activity modeling stage takes into account the imaging protocol, its type (in-beam, in-room etc) and the delay times between the irradiated fields.

In our simulations we have tested a scenario in which imaging was performed after the last irradiated field. This last field had the greatest impact on the production of β⁺-emitters and indirectly on the activity distribution, as it decays the least compared to the previously deposited fields.

This effect can be seen in figure 3—where the β⁺-emitters production map (image in the middle) is compared with the modeled activity from all fields (image on the right). The activity distribution in figure 3 largely corresponds to that of Field 4 shown in figure 4.

All Fields separately are presented in figure 4 as β⁺-emitter production maps (upper row) and as a modeled activity distribution after Stage 3 (bottom row). Here we can compare and clearly see what is the contribution of the β⁺-emitter production from each Field to the final activity. In the bottom row, the level of activity increases for the consecutive fields, with the greatest values for the Field 4.

Stage 4 output. The second MC simulation of β⁺ emission and PET acquisition gives us the hits in the detector volume and coincidences lists which are stored in a ROOT file format that is not feasible to present in a graphical way. Both types of computing nodes were used for this simulations and it took about 90 [min] using 200 CPU cores with older Westmere nodes or about 30 [min] with Haswell nodes.

Stage 5 output. The last figure in this section, figure 6, presents the output of the CASToR image reconstruction software as a final reconstructed image in three planes (axial, coronal and sagittal). The images presented were postprocessed by using a Gaussian smoothing filter with sigma = 2 voxels, in order to remove statistical noise.

Zoom In Zoom Out Reset image size

The ProTheRaMon framework provides the user with the ability to perform MC simulations of proton therapy treatment and secondary radiation detection for range monitoring using PET. However, instead of creating a single-step full simulation that provides a unique output for a given set of input parameters, we have decided to design our framework as a set of separate and independent stages. This has allowed us to prepare a more flexible and versatile tool. This design also allows the user to decide on a range of relevant parameters at each step. In the second stage, we have decided to score the seven radioisotopes listed in table 1 and also shown as an example of production maps in figure 5. However, if any of these are not of interest in the specific study, they can be switched off or other isotopes can be added. Also, a variety of different PET imaging protocols that differ depending on the capabilities of specific proton centers and specific scanner designs can be proposed. Off-line, in-room or in-beam (if time information is incorporated) imaging protocols may be implemented to model the PET signal decay in Stage 3. The protocol we have used is presented in figure 4, where the weakening (by decay) of activity from earlier fields is apparent. In our study we have not implemented a biological washout model, only radioactive decay with the assumption that all β⁺-emitter radionuclides are produced at the beginning of the field irradiation. The effect of biological washout can be added by the user in Stage 3 according their needs. An example of the method that can be additionally and easily implemented by the user into the proposed framework can be found in (Mizuno et al 2003).

Two stages of the ProTheRaMon workflow, i.e. Stage 2 and Stage 4, exploit MC simulations that are used for different purposes. In Stage 2 (see 2.1.2), simulations of a treatment plan are performed in order to obtain distributions of β⁺-emitting isotopes produced during the irradiation. The transport of photons, resulting from positron annihilations, out of the patient body and the photon energy deposition in the PET detectors is simulated in Stage 4. Therefore, the two stages require different setting of physics options which can be set individually for each stage of MC simulation. The production of β⁺-emitting isotopes during the proton irradiation requires physics models incorporating hadronic and neutron interactions, while the detection of γ particles in the PET scanner detector can be done with electromagnetic models only. The exclusion of unnecessary physics processes has the effect of decreasing the simulation time.

Stages 4 and 5 require the definition of a PET scanner geometry, which can be optimized and changed independently from the β⁺ activity simulation performed in Stage 2. By changing one macro file with the geometry description the user can easily implement many different detector designs.

The treatment and PET imaging parameters may be varied and investigated in a more time-efficient manner by using a staged architecture implemented in ProTheRaMon. For example, at the activity production simulation stage (Stage 2), we simulated 10% out of the total number of particles used in the treatment. This was motivated by the need to optimize the simulation time with limited computational resources. For our purposes, this level of statistics represented a satisfactory compromise between the variance of the isotope production and the simulation time. The other settings that will significantly affect calculations times are the resolution of the voxelized phantom and source as well as the reconstructed image resolution. Such parameters should be adjusted to the available computational resources and the performance of the PET detector.

The uncertainty of β⁺ activity modeling is essential when validating the performance of a new PET detector for proton beam range monitoring. The uncertainty of Monte Carlo simulations is intrinsically related to the uncertainty of the nuclear model used (Sihver et al 2012). The experimental validation of Geant4 and other MC codes in terms of the yields of β⁺-emitters during proton therapy has been the subject of extensive research (for a review of the subject see (Kraan 2015)) and has resulted in the tuning of the relevant cross section data. In simulation of PET activity maps in patients, further inaccuracies may result from modeling of the clinical proton beam, which for our example has been recently addressed in Gajewski et al (2021), while the uncertainty related to the modeling of patient tissue composition based on CT scans has been addressed, e.g. by Paganetti (2012). Overall, these uncertainties have a minor impact on the application of the ProTheRaMon framework to simulation studies aiming at PET detector and protocol designs because these studies are based on a relative comparison of the effectiveness of different approaches.

ProTheRaMon was developed to simplify simulation studies needed to design, test and validate PET scanner geometries and protocols dedicated to proton therapy range monitoring in a clinical environment. ProTheRaMon is an automated multi-stage protocol designed for the following tasks: to simulate the production of β⁺-emitting radioisotopes originating from nuclear reactions of the primary proton beam with the patient during proton therapy, to simulate the emission and transport of the gamma rays originating from electron-positron annihilation to the PET detector, and to reconstruct the PET image. ProTheRaMon offers complete processing of proton therapy treatment plans and patient CT geometries in the coordinate system of the treatment room and PET scanner, taking into account the activity decay related to the PET imaging protocol specific to the proton therapy center.

So far, any institute wishing to perform Monte Carlo simulations of range monitoring in proton therapy using PET must set up their own simulation workflow. The ProTheRaMon framework can simplify this task as it is available to the community via the github repository ProTheRaMon (Borys et al 2022). Moreover, the modular design of the ProTheRaMon framework allows for easy adjustment and expansion of its functionality for other proton therapy range monitoring applications, where two basic MC simulation steps, i.e. β⁺-emitters production and PET imaging, have to be performed.

ProTheRaMon was designed to execute the most computationally intensive stages on a computational cluster. Due to the time requirements for performing the simulations yielding β⁺-emitter distributions and β⁺ emission simulations, the use of a computational cluster is essential to perform the study on a large group of patient treatment plans.

This work was supported by the National Centre for Research and Development (NCBiR), Grant No. LIDER/26/0157/L-8/16/NCBR/2017. Calculations were performed on the Ziemowit computer cluster in the Laboratory of Bioinformatics and Computational Biology at the Biotechnology Centre, Silesian University of Technology, created in the EU Innovative Economy Programme POIG.02.01.00-00-166/08 and expanded in the POIG.02.03.01-00-040/13 project. Development of FREDtools was supported by the National Centre for Research and Development (NCBiR), Grant No. LIDER/43/0222/L-12/20/NCBR/2021. We acknowledge also support by the Foundation for Polish Science through the TEAM POIR.04.04.00-00-4204/17 program, the National Science Centre of Poland through grant 2021/42/A/ST2/00423, and the SciMat and qLIFE Priority Research Areas budget under the program Excellence Initative—Research University at the Jagiellonian University, and Jagiellonian University Project No. CRP/0641.221.2020. This work was also supported by the Swiss National Science Foundation, Grant No. CRSII5189969.