Up to now, carbon ions have shown the most favorable physical and radiobiological properties for radiation therapy of, for example, deep-seated radioresistant tumors. However, when carbon ions penetrate matter, they undergo inelastic nuclear reactions that give rise to secondary fragments contributing to the dose in the healthy tissue. This can cause damage to radiosensitive organs at risk when they are located in the vicinity of the tumor. Therefore, predictions of the yields and angular distributions of the secondary fragments are needed to be able to estimate the resulting biological effects in both the tumor region and the healthy tissues. This study presents the accuracy of simulations of therapeutic carbon ion beams with water, with the 3D MC (Monte Carlo) general purpose particle and ion transport code PHITS.
Simulations with PHITS of depth-dose distributions, beam attenuation, fragment yields, and fragment angular distributions from interactions of therapeutic carbon ion beams with water are compared to published measurements performed at Gesellschaft für Schwerionen Forschung (GSI).
The results presented in this study demonstrate that PHITS simulations of therapeutic carbon ion beams in water show overall a good agreement with measurements performed at GSI; for example, for light ions like H and He, simulations agree within about 10%. However, there is still a need to further improve the calculations of fragment yields, especially for underproduction of Li of up to 50%, by improving the nucleus-nucleus cross-section models.
The simulated data are clinically acceptable but there is still a need to further improve the models in the transport code PHITS. More reliable experimental data are therefore needed.
Carbon ion therapy has shown to be very promising for treating a wide variety of tumors, especially deep-seated and radioresistant ones, and by the end of 2016 nearly 20 000 patients had been treated with carbon ions worldwide . This success is due to a number of advantages when using heavy ions for radiation therapy compared to photon or proton radiation therapy. The main benefit of using accelerated heavy charged particles for radiation therapy stems from their favorable physical characteristics. When energetic heavy charged particles penetrate an absorbing medium, they are slowed down by losing their kinetic energy mainly through ionization of the medium. The energy loss of the particle during its travel through a material is called Bragg curve after William Henry Bragg, who discovered this phenomenon in 1903. The specific ionization increases with decreasing particle velocity, giving rise to a sharp maximum (Bragg peak) in ionization near the end of range. This property makes it beneficial to use charged particles, for example, protons or heavier ions, to treat deep-seated tumors. To treat an extended target, the Bragg peak is spread out to cover the required volume by superimposing many pristine Bragg curves through modulating the energy of the particles. One additional advantage of heavier ions compared to protons is that multiple scattering and range straggling vary approximately inversely to the square-root of the mass of the particle . The sharpness of the lateral dose fall-off, often called apparent penumbra, is of clinical importance because the radiation exposure to normal tissues adjacent to the target volume often limits the therapy dose, especially for a tumor close to critical organs.
Heavy ions have higher linear energy transfer (LET) than photons. They therefore cause larger ionization density when travelling through human organs and thus a considerably increased relative biological effectiveness (RBE), particularly in the Bragg peak region where the tumor is located, in comparison to therapeutic photons and protons. Following exposure to ionizing radiation, if molecular oxygen is present, organic peroxide and other reactive oxygen species are produced, thus making permanent the damage incurred by the DNA. Under hypoxic conditions, however, DNA damage induced by low LET radiation can be more readily repaired. The densely located ionization events from heavy ions cause more complex clustered DNA damage than the sparsely distributed events from photons and are therefore less influenced by the oxygen content in the cells. The so-called oxygen enhancement ratio is hence low. High LET radiation is therefore more efficient to kill hypoxic and radioresistant tumor cells than other types of radiation therapy. Ionizing radiation with low oxygen enhancement ratio is also less dependent on the cell cycle, suggesting that it may be effective for treating photon-resistant late–S-phase cells. S-phase is the part of the cell cycle in which the DNA is replicated, occurring between the G1 and G2 phases. The modes of carbon ion–induced cell death and inactivation include apoptosis, necrosis, autophagy, premature senescence, accelerated differentiation, delayed reproductive death of progeny cells, and bystander cell death . Since high LET radiation is less dependent on reoxygenation of the tumor environment and reassessment of cells into radiosensitive phases of the cell cycle between the dose fractions, high LET radiation is also less dependent on fractionation.
When high energetic ions with Z > 1 penetrate tissue they undergo inelastic nuclear reactions, which give rise to secondary projectile and target fragments. The heavier the projectile, the more fragments are created. These fragments travel about in the same direction as the primary beam. Owing to the dependency on A/Z2 of the range of these particles at the same velocity, many will penetrate further than the tumor volume, therefore causing damage to the healthy tissue. When treating tumors very close to radiosensitive organs at risk (OARs), such as the optic nerve and brain stem in brain cancer or the rectum in prostate cancer, it is especially important to achieve high tumor-dose “conformality” with low dose deposited to the OARs for reduction of treatment toxicities. The first clinical trials with energetic heavy ion beams for cancer treatment were performed at the Lawrence Berkeley Laboratory between 1977 and 1992 [4, 5]. These trials were mainly performed with helium and neon ions, but also with carbon ions. Some trials were also done with argon and silicon ions, but owing to the observed severe side effects, treatment trials with these beams were stopped. In the beginning of the 1990s, the Heavy Ion Medical Accelerator in Chiba (HIMAC) was constructed at the National Institute of Radiological Sciences in Chiba, Japan. The first clinical trials at HIMAC started in 1994. During the construction of HIMAC, the RBE of different ions, including carbon, neon, silicon, argon, and iron at different LET were studied. It was clearly found that tumor cells vary in their LET-dependent responses. Carbon ions were chosen for cancer therapy at HIMAC because they appeared to have the optimal biological and physical properties. For carbon ions, the strongly elevated DNA damage is restricted to the end of the particle range, where there is a lot of irreparable clustered DNA damage. In the entrance channel, where healthy tissue is located, most of the DNA damage can be repaired. For heavier ions like neon, silicon, and argon, the irreparable damage becomes more common in the entrance channel and can therefore create normal tissue complications. More fragments, which can damage OAR beyond the tumor region, are also created when heavier ion beams are used. The maximum in RBE for carbon ions nearly coincides with the Bragg peak, which is also an advantage for tumor therapy. However, RBE depends in a very complex manner on physical and biological parameters, including particle type, ion beam energy, LET, dose, dose rate, as well as the cell and tissue type. In particular, such complex dependence requires a reliable description of the mixed radiation field, which is produced by fragmentation reactions when heavy ion beams penetrate a medium.
To model in detail the physical interactions of ion beams in tissue, and couple this information with proper models for estimation of biological response, Monte Carlo (MC) codes offer the most accurate computational engines. However, owing to the scarcity of experimental fragmentation data for the materials and energy range of therapeutic interest, uncertainties remain in the description of the nuclear processes underlying the production of the fragments, which have varying biological effectiveness and are likely to deposit dose in the OARs and healthy tissue outside the tumor. Simulations of fragment yields in water of a therapeutic carbon ion beam, using the MC particle transport codes FLUKA and Geant4, and comparisons with the measurements performed at Gesellschaft für Schwerionen Forschung (GSI) by Haettner et al  and Haettner , have been previously published . In this article, simulations of the same experiment using the 3D MC Particle and Heavy Ion Transport code System (PHITS) version 2.88  is presented. We show comparisons of depth-dose distributions, fragment yields, and angular distributions from the interactions of 400-MeV/u carbon beam in water. We also compare simulated and measured beam attenuation of 200- and 400-MeV/u carbon ions in water. These simulations are fundamental to estimate accurately the RBE in complex mixed radiation fields for reliable patient treatment planning.
Materials and Methods
The simulations of beam attenuation, depth-dose, angular distributions, and fragment yields presented in this study are compared to measurements performed by Haettner et al  and Haettner  at the GSI. During the experiments, 200 and 400 MeV/u 12C were delivered in slow extraction mode from the heavy-ion synchrotron SIS-18 to the biophysics area and passed through a thin aluminum vacuum exit window to the target area. The pencil-like beams with approximately Gaussian-shaped horizontal and vertical profiles were focused to a spot size of about 5 mm full width at half maximum at the target position. A water target, enclosed in a polymethyl-methacrylate box of variable thickness, was used at room temperature for all measurements. The depth-dose distributions of 12C ions in water were measured with 2 ionization chambers: a small one, placed behind the beam exit window for beam current normalization, and a second one located downstream of the water target, as schematically shown in Figure 1A. The energy spectra, yields, and angular distributions of secondary fragments emerging from the water target were investigated by simultaneous measurements of the energy loss (ΔE) and time-of-flight. A thin scintillation paddle (NE102A, 10 × 10 cm2, 1.5 mm thick) was used as beam counter and start detector. For the (ΔE) measurement and as a stop detector, a plastic scintillation detector (NE102A, 40 × 40 mm2) coupled to a Hamamatsu (Hammamatsu, Japan) H6779 photomultiplier tube was used. A detailed description of the measurements is described elsewhere [6, 7].
The simulations of the experiments [6, 7] presented in this study were performed with the 3-dimensional multipurpose particle and heavy ion MC transport code PHITS, version 2.88 . The geometry and conditions of the experimental setup used by Haettner et al  and Haettner  at GSI were reproduced in PHITS. Corrections for beam spread between the target and the detectors were made in the PHITS simulations, based on the measurements presented in Haettner et al . In PHITS, the Kurotama model [10, 11] was used for the total nucleus-nucleus reaction cross section, and the JQMD-2.0 (JAERI Quantum Molecular Dynamics v. 2.0 ) [12–14] for the first step of the nucleus-nucleus interactions. The nucleus de-excitation algorithm used in the simulations was the general evaporation method, or GEM .
Comparisons of Depth-Dose Distributions and Beam Fluence Attenuation
In Figure 2, the Bragg curves of 400-MeV/u 12C ions in water are displayed as predicted by PHITS together with experimental data [6, 7]. The simulated contributions from the primary carbon ion beam and the secondary particles are shown, together with the total depth-dose distributions. In the lower part, the depth-dose distributions for the fragments are separately shown for each charge (Z). PHITS is found to reproduce the whole Bragg curves very well.
In Figure 3, the measured attenuations of 200- and 400-MeV/u 12C ion beams in a thick water absorber are shown, together with the simulated attenuations. As can be seen, PHITS reproduces the measured attenuations well with an average difference between the measurements of less than 5%.
Comparisons of Fragment Yields at Different Depths in Water
In Figure 4, the measured fragment yields from the reactions of 400 MeV/u 12C in water as predicted with PHITS are shown together with measured values. As can be seen, the PHITS simulations underestimated the yields of Li and B, and slightly overestimated the yield of He. The yields of H and Be were reproduced quite well. The estimated differences between the measured and calculated values are listed in the Table.
Comparisons of Angular Distributions of Charged Fragments
In Figure 5, the fragment angular distributions as predicted with PHITS are shown, together with measured values [6, 7] from the reactions of 400-MeV/u 12C ions after (A) 15.9 cm and (B) 31.2 cm of water. The estimated differences between the measured and calculated values are listed in the Table.
In this article, PHITS simulations of depth-dose distributions, fragment yields, and angular distributions from the interactions of therapeutic (400 MeV/u) carbon ion beams with water are compared to measurements performed at GSI by Haettner et al  and Haettner . Comparisons of simulated and measured attenuation of 200- and 400-MeV/u carbon ions in water are also presented. Simulations of the measured [6, 7] fragment yields from the interactions of 400-MeV/u carbon in water, using FLUKA and Geant4, have previously been published . That study showed that the used nuclear models of FLUKA and GEANT4 could reproduce the fragment yields with reasonable accuracy. However, discrepancies in the order of some tens of percent were found and the conclusion was that the models should be improved for their application to carbon ion therapy.
From the study presented in this article, it can be seen that PHITS reproduces the measured Bragg curves quite well, but there are some deviations from the measured values for the fragment yields (up to 40%) and the fragment angular distributions (up to 50% for Li). However, the uncertainty of the experimental data for fragment production, dominated by systematic errors such as ambiguous identification and angular resolution, is within 25%. Hence, the level of agreement between the measured and simulated data is acceptable and is of the same order as the previously published results with FLUKA and Geant4 . Since the dose contribution from the fragments is low, the simulated data are clinically acceptable but there is still a need to further improve the calculations of fragment yields by improving the nucleus-nucleus cross-section models. To improve and benchmark the transport codes and their models, more experimental data are however needed.
ADDITIONAL INFORMATION AND DECLARATIONS
Conflicts of Interest: The authors have no conflicts of interest to disclose.
Acknowledgments: This work was partly financially supported by funding from European Community's Seventh Framework Programme (EURATOM) contract FP7-295970 (ANDANTE). T.T. acknowledges funding from DFG (KFO214).