Cancer is the second leading cause of death for children and one of the treatment options for this disease is radiotherapy. Children treated with radiotherapy using photon beams are more likely to develop secondary neoplasms. Proton therapy can reduce the probability of neoplasm formation by up to 50%. Recent studies propose the use of ultra high dose rates as a treatment option. From the threshold of 40 Gy/s it is possible to reach the FLASH effect. This technique protects healthy tissue while maintaining tumor control. The effect was validated in vivo using a proton beam and, therefore, it will be available as a new treatment option. On the other hand, the proposal for FLASH treatment with a proton beam would not use the Bragg peak located in the target volume, which is the differential of proton radiotherapy. In addition, the increase in the intensity of the beam and the energy of the particles, lead to the generation of a greater amount of neutrons. The objective of this work is to evaluate the dose due to the neutrons generated in the interaction with the accelerator components in FLASH proton therapy in relation to conventional proton therapy. The dose evaluation was performed through Monte Carlo simulations, using a water phantom, with the code TOPAS MC. The results found show that the dose of neutrons in the FLASH technique would be about 100 times greater than the dose in the conventional technique. Still, it would be below 1% of the prescribed dose.

Cancer is the second leading cause of death among children. 25% of them will need radiotherapy treatment [1]. Furthermore, young patients, from 1 to 19 years old, who receive this treatment modality have a 2.7 times greater risk of developing secondary neoplasms. On the other hand, proton beam radiotherapy can reduce the incidence of these radio-induced neoplasms by 50% [2]. Therefore, countries that have proton therapy centers tend to treat some pediatric tumors with this method [3].

The use of protons in medicine is nothing new. One of the first articles dates from 1946, when Wilson proposed using high-energy beams. Previously it was barely used due to techitsleogical limitations in particle acceleration, not allowing greater beam penetration. However, several improvements in the accelerators, part of them using the radar technology developed during the Second World War, allowed an immense advance in the final energy of the accelerated particles [4].

Currently, proton beam treatments are performed using passive scattering, uniform scattering or pencil-beam scanning (PBS) techniques. The existence of three techniques suggests, on one hand, an evolving field and, on the other hand, that none of these methods is clearly superior to the others for all clinical indications.

Proton therapy has enormous potential, which is easily shown under quality control reference conditions in radiotherapy, such as a water tank. However, in many important aspects like image orientation, workflow, and others, it is not a mature technique yet. This conclusion is most evident when dealing with treatment-related uncertainties in patients, especially those that affect the proton range [5].

The concern with the secondary dose from neutrons and photons is a recurring issue in the literature. These radiations are generated by the interaction of the proton beam with the structure of the beam line or with the patient itself.

Matsumoto et al. used the Particle and Heavy Ion Transport code System (PHITS) to model and simulate the treatment line at the National Cancer Center Hospital East. In addition to the MC code, they also used a hybrid computer simulator. As a result, the dose deposited by neutrons originating in the collimator is about 30 times greater than the dose delivered by those originating in the beam transport line. Therefore, the study suggests that the doses produced by secondary neutrons can be avoided with a more careful design of the beam transport line structures [6].

Schneider et al. performed simulations in Monte Carlo using the FLUKA code and measurements in a water simulator. The study aimed to analyze the dose delivered by secondary neutrons during a proton therapy treatment. The results showed that the equivalent dose from neutrons in the treatment volume corresponds to about 1% of the treatment dose. Furthermore, at medium and large volumes (outside the treatment region), the neutron doses were approximately 0.004 Sv and 0.002 Sv, respectively, per treatment. Thus, the study showed that the doses in the treatment region can be neglected and that the highest dose in the treatment region, when compared with photon radiotherapy, is still insufficient, considering that the proton beam treatment has a lower dose in healthy tissues [7].

A new treatment technique has recently emerged: FLASH radiotherapy. This technique uses ultra-high dose rates (D > 40 Gy/s) to preserve healthy tissue adjacent to the target volume while maintaining tumor control [8-10]. The scientific community was still skeptical about whether it would be possible to achieve the FLASH effect utilizing a proton beam, even though it had already been proven with electrons and with photons.

Diffenderfer et al., designed, implemented, and validated, in vivo, the FLASH technique using a proton beam. Some changes were implemented in the therapeutic beam line, from the double scattering technique, such as the use of scattering foils, the increase in proton energy, and the use of the portion prior to the Bragg peak for the treatment. These modifications justify the resumption of the discussion of dose by neutrons generated by the materials of the beamline [11].

Therefore, the objective of this work is to perform neutron dose simulation in organs and tissues adjacent to a hypothetical target volume submitted to FLASH proton therapy treatment for head and neck tumors using the Monte Carlo method. Also, compare the simulation values ​​with a conventional proton therapy passive scattering technique.

All simulations were performed using the Monte Carlo method, with the TOPAS MC code. The software is free for academic use and was created with the purpose of wrapping Geant4 and expanding its use among the medical physicists community, without the need of C++ knowledge [12].

The code was primarily designed to simulate clinical proton beams and then it was expanded to other radiation particles. Thus the default TOPAS physics list was used for the simulations. The default list (QGSP_BERT_HP) has been shown to work well for proton therapy simulations, both for the primary particle and secondary particle. The Bertini model produces more secondary neutrons and protons with a better agreement to experimental data. Furthermore it was used the high precision neutron library (neutronHP), which includes precise cross section data for neutrons below 20 MeV to thermal energies. No modifications or energy cuts were made to the physical processes included in the library. The validation of physics list settings for low-energy proton therapy (up to 250 MeV) has been studied by Jarlskog and Paganetti [13].

As the purpose of this work is to compare two types of proton therapy, two configurations and geometries were used:

A water box, measuring 100x100x400 mm³ (X, Y, Z), was choosen as phantom. It was divided in 51x51x100 voxels, respectively. Both X and Y axis were divided in an even number in order to evaluate a central column of voxels in the beam propagation direction (Z). The box composition was set as water accordingly to the Geant4 material library: G4_Water.

For each simulation, one may store the desired physical quantities in each voxel of the geometry or in the entire volume, if it was not previously divided. In this work, the physical quantities were scored in the water box voxels. The output data are presented as a 4 dimension array, where each element has a value associated with that coordinate.

The conventional passive double scattering proton therapy line was base on the reference study about the production of neutrons from a proton line in Leuven, Belgium [14]. The study provided detailed geometry data using the same double scattering technique used in proton therapy treatments and therefore could be used as a reference. The geometry and configuration of the beam line was reproduced in order to obtain the neutrons and protons energy fluence as well as the absorbed dose by neutrons and protons in voxels representing the planning target volume and adjacents organs or tissues at risk. It is worth mentioning that the paper used as reference also contains simulated data about the neutron spectrum that were used to evaluate the simulation results of this work. The schematic diagram of the simulation can be seen in the figure 1. More details, as dimensions and materials, have been described elsewhere [14].

The simulation starts with a 200 MeV proton beam placed into a vacuum tube. Immediately after the beam it was placed the first scatterer, or range modulation wheel (RMW), and then the second scatterer, which has also the function of flattening the beam profile. The environment was configuredas atmospheric air, in the same way it is in the experimental beam. Finally, the collimators set and the snout were inserted in the simulated line.

Additionally to the voxel scorer, in this simulation, it was also used two dosimeters rings in order to measure the spectrum of neutrons generated by the interaction of the proton beam with the beam line. The first one with an inner radius of 20 mm and an outer radius of 40 mm and the second one with an inner radius of 120 mm and an outer radius of 140 mm. Both rings were positioned after the first collimator. The dosimeter material was set as vacuum so that the neutrons energy fluence could be measured without the influence of their interaction with the dosimeter itself, even though it would not be possible in a real measure. The neutron energy fluence was obtained, in each dosimeter ring, normalized by the number of simulated protons. The neutron spectrum data obtained were compared with literature data in order to evaluate the simulation.

A total of 1 billion histories were used, divided into 8 simulations. Each simulation corresponds to one step on a total of 8 steps of the range modulation wheel (RMW), in order to obtain a spread out Bragg peak (SOBP). The steps correspond to a SOBP of 50 mm. The dimensions of the RMW materials for each of the steps can be seen in Table 1. The PMMA material information was not available in Geant4 material library, thus its composition was set as : 0.5998 carbon, 0.3196 oxygen and 0.0806 hydrogen [15]. The density of PMMA was set as 1.180 g/cm³ [16]. The simulation outputs were particles energy fluence and absorbed dose, both for protons and neutrons, for each RMW step. Additionally the data from each step were multiplied by a weight factor regarding the contribution of each curve on the spread out Bragg peaks formation. The method used for calculating the weight factors is described elsewhere [17] .

The proton beam was used to obtain the physical quantities in each phantom voxel, specially in three volumes of interest: a hypothetical treatment target corresponding to the SOBP of 10x20x50 mm³ - composed by several voxels - and two volumes representing an organ or tissue at risk at 2 cm and 5 cm from the treatment target, on the X axis. Both risk volumes are located on the X-Y plane at the maximum neutron dose depth, i.e. 1 cm on the Z axis. Absorbed doses from protons and neutrons were obtained separately. The simulated doses were normalized by the dose delivered in the target volume by protons, thus it represents the absorbed dose in a 1 Gy prescription.

The configuration of the beamline used in the simulation for FLASH proton therapy was based on the same one used by the group Diffenderfer et al., in which the authors validated the FLASH effect in vivo (figure 2). In this modality, the FLASH effect was obtained using the portion of the beam prior to the formation of the Bragg peak.

The proton beam has initial energy of 230 MeV. After the beam exit window, a first scatterer was inserted, composed of 2 mm of lead. Then a second 5 mm scatterer, also made of lead. Finally, the primary collimator, made out of brass, with 1 cm horizontal and 2 cm vertical openings. The primary collimator is thick enough to stop the highest energy protons CSDA range data obtained from PSTAR/Nist software [18]. For lead and brass, the following materials from the Geant4 library were used: G4_Pb and G4_BRASS, respectively.

In the case of the FLASH modality, the Bragg peak was not used, since the technique validated in vivo uses only the transmitted part of the proton beam. For the same reason, a RMW was not used to produce a SOBP. The simulations scored, in each voxel, mean values ​​and standard deviations of: absorbed dose due to neutrons and protons as well as the neutron absorbed dose binned by the incident neutron energy (from 0 to 10 MeV, in 100 bins of 0 .1 MeV).

For comparison purpose, the volumes of interest were determined in the same location as in conventional proton therapy, 2 cm and 5 cm from the center of the phantom, on the X axis, and 1 cm deep on the Z axis. The neutron doses were normalized by the proton dose at a depth of 1 cm in the central voxel of the phantom. Thus, the doses found, as in conventional therapy, also reflect the neutron dose values ​​as a function of 1 Gy of proton dose prescription.

The simulation results are presented separately for traditional proton therapy, with the passive scattering technique of protons of 200 MeV and FLASH proton therapy, with the technique of double scattering and using the transmission of the proton beam of 230 MeV.

The acquired neutron spectra for both dosimeters rings were superimposed to previous results from the paper used as reference. It was not possible to compare it quantitatively due to lack of reference numerical data. The graphics showing the qualitative comparison can be seen in figures 3 and 4.

For evaluation of the neutron absorbed dose, the percentage depth dose can be visualized in Figure 5. The map of relative absorbed dose for neutrons in the X-Z plane, centered on Y can be see in figure 6.

The results of neutron dosimetry in conventional proton therapy in the volumes of interest adjacent to the target volume, for a dose prescription of 1 (reference value), 30 and 60 Gy, can be seen in Table 2.

For beam evaluation, the relative dose map of the FLASH proton therapy simulation for the X-Z ​​plane centered on the Y axis and also the relative proton fluence map, both normalized by the maximum value, can be seen in figures 7 and 8, respectively.

The relative dose for depths less than 5 cm is shown in figure 9. This depth was chosen because it is equivalent to the depth studied in the proton beam FLASH effect validation paper.

For the evaluation of neutrons in the phantom, the dose map was acquired for the X-Z plane centered on the Y axis. The dose values ​​were normalized by the maximum value and are shown in figure 10.

In addition, the relationship between the neutron dose and the depth in the simulator and the lateral position, that is, the distance from the central axis of propagation of the proton beam, were evaluated. For the first evaluation, the percentage depth dose (PDD) curve was obtained along the Z axis (proton beam propagation axis). The PDD curve for neutrons can be seen in Figure 11. For the second evaluation, ie the lateral dose deposition, the average neutron absorbed dose profile curve over the entire phantom was obtained and can be seen in Figure 12.

Finally, the results of neutron dose as a function of proton dose for the simulation, considering FLASH proton therapy with prescriptions of 30 and 60 Gy, can be seen in Table 3.

The authors who validated FLASH proton therapy did so using the transmission of a proton beam [11]. Among several justifications, the Bragg peak would not be useful for superficial treatments. In addition, the time required for some active scattering techniques may lead to lower dose rates than the required to trigger the FLASH effect. However, the use of the proton transmission beam de-characterizes all the advantages that proton therapy has over the use of photon beams, for example. Therefore, a study that takes into account more factors, even economic ones, to investigate the real feasibility of FLASH proton therapy is needed.

The data from the simulations referring to the evaluation of the neutron spectrum showed differences, mainly in the high energy range (above 100 MeV). In the first place, it may be that this difference between high energies does not change the final result significantly, since the neutrons of the energy range in question have a relatively lower cross section than the lower energy neutrons. However, it would be necessary to assess it quantitatively in order to evaluate the effect of the spectrum differences on dose values. One way of evaluating this contribution of higher energy neutrons would be to obtain the dose values ​​per neutrons with separation in energy bins.

As expected, the doses evaluated in the volumes of interest, by the passive scattering method, scale with the increase in the prescribed dose. The evaluation performed is totally dependent on the therapeutic technique. Alternatevely active scattering techniques can be studied in the future.

The results achieved for the FLASH technique beam profile are not compatible with the experimental results presented by Diffenderfer et al [11]. The depth doses presented by the group remain constant, contrary to the result found in the simulation, which declines it. The investigation carried out by studying the fluence of protons in the phantom suggests that the divergence of the proton beam would be responsible for the dose drop along the depth. Which, in turn, suggests that the geometry indicated in the reference article is not fully clarified or not reproducible. The use of a concave-shaped compensator could contribute to the flatness of the beam, reducing its divergence in the material.

Because the dose at depth is not constant, as can be seen in figures 7 and 9, the reference dose of protons was chosen as the one at 2 cm depth, since it is approximately the average between the dose at depths 0 cm and 4 cm. For the calculations performed, this fact is not of great importance, since the interest resides almost exclusively in the large order of the values and the divergence of the numbers is approximately 7.5% more or less.

Simulated absorbed doses are estimates of physical quantities and do not necessarily reflect biological effects. The evaluation of other quantities would be important, mainly due to the strong interaction between neutrons in the energy range from 100 keV to 2 MeV and biological tissue, with a weight factor for equivalent dose 20 times greater than photons. In addition, the simulation was performed in a water simulator. Heterogeneity issues and radiosensitivity factors of the organs and tissues adjacent to the tumor that can have a great influence are not taken into account.

The neutron doses of the FLASH technique are 100 times greater than those of the conventional technique. Even so, the neutron doses for both techniques evaluated are within the limits described in the literature. Furthermore, in an extreme case, they would be below 1% of the prescribed dose value. The doses are higher when close to the phantom center, contrary to the conventional technique.

The neutron dose is directly proportional to the prescribed dose of protons. In addition to being directly proportional to the size of the target volume in the conventional technique, the technique uses spreader sheets to increase the beam cross-section and RMW to form the SOBP.

For future work, it will be necessary to evaluate the dose taking into account the weighting factors of the type of radiation and the weighting factors of radiosensitivity of the evaluated tissues and organs. In addition, it may be an option to use anthropomorphic virtual simulators. Finally, future works should also evaluate and compare dose deposition with other application techniques.