Monte Carlo investigation of the characteristics of radioactive beams for heavy ion therapy

This work presents a simulation study evaluating relative biological effectiveness at 10% survival fraction (RBE10) of several different positron-emitting radionuclides in heavy ion treatment systems, and comparing these to the RBE10s of their non-radioactive counterparts. RBE10 is evaluated as a function of depth for three positron-emitting radioactive ion beams (10C, 11C and 15O) and two stable ion beams (12C and 16O) using the modified microdosimetric kinetic model (MKM) in a heterogeneous skull phantom subject to a rectangular 50 mm × 50 mm × 60 mm spread out Bragg peak. We demonstrate that the RBE10 of the positron-emitting radioactive beams is almost identical to the corresponding stable isotopes. The potential improvement in PET quality assurance image quality which is obtained when using radioactive beams is evaluated by comparing the signal to background ratios of positron annihilations at different intra- and post-irradiation time points. Finally, the incidental dose to the patient resulting from the use of radioactive beams is also quantified and shown to be negligible.

A large number of annihilation photons must be detected in order to obtain a PET image of sufficient quality for useful treatment QA. The cross-sections for inelastic ion collisions depend on several parameters, including incident ion species and energy, and the density and composition of the target 17 . These factors determine the mix of fragments produced, which, in turn, determines the number and distribution of positron-emitting radionuclides resulting from each beam spill. To improve image quality, several authors have proposed the use of positron-emitting radioactive nuclei (such as 11 C, 15 O or 10 C) as the primary particle in the heavy ion beam. Most primary particles will survive intact to decay via positron emission at their stopping point, corresponding to the location of the Bragg peak. Therefore, for radioactive beams, the spatial distribution of the stopping points of primary particles is the dominant component of the PET image, while positron-emitting target and beam fragments making up a secondary component.
Beamlines capable of producing beams of radioactive ion species such as 11 C, 10 C and 15 O with sufficient dose rates and beam purity for therapeutic use are currently under development at the National Institutes for Quantum and Radiological Science and Technology (NIRS, QST) in Japan and other facilities around the world [18][19][20][21][22][23] .
In order to perform proper treatment planning with positron-emitting radioactive beams, and to understand how their use will impact image-based QA, it is necessary to address three key research questions: 1. How does the relative biological effectiveness (RBE) of polyenergetic radioactive beams vary as a function of depth within a spread out Bragg peak, and how does this compare to the corresponding stable ion species? 2. What quantitative differences are expected between the maps of positron annihilation resulting from treatment with stable and positron-emitting radioactive ion beams, and how will these impact the use of PET images as an intra-treatment or post-treatment QA mechanism? and finally, 3. What additional dose will be received by the patient if a positron-emitting radioactive beam is used instead of a stable beam?
In this work, simulations of a simple treatment plan (consisting of a flat biological dose in a rectangular-prismatic primary treatment volume inside a human skull phantom) are performed for five primary nuclei (three positron-emitting and two stable) using the Geant4 Monte Carlo toolkit. The values of RBE 10 (RBE at 10% survival fraction) are estimated across a range of depths along the beam path (in the entrance, SOBP and tail regions) using Kase's modified microdosimetric kinetic model (MKM) [24][25][26] . The validity of using Monte Carlo simulations to evaluate RBE using the MKM has previously been established by Bolst et al. 27,28 ; however, to our knowledge, this is the first time that this approach has been applied to estimate the RBE 10 of a polyenergetic radioactive beam. The method can easily be extended to other homogeneous or heterogeneous targets and heavy ion species, and is a convenient and cost-effective alternative to in vitro experiments.
Monte Carlo simulation-based 2D maps of positron yield obtained in a skull phantom using a spread out Bragg peak (with the same flat biological dose (in Gy(RBE)) delivered throughout the planned treatment volume) are compared across all beam types. The distribution of positron production in the target volumes, as measured during the beam-off periods during irradiation of the phantom with the radioactive and corresponding stable heavy ion beams were measured, and the resulting signal to background ratios (SBRs) estimated. The chosen physics models in the simulation are validated via experimental work conducted at NIRS's HIMAC facility.
Finally, the additional dose to the patient resulting from the use of radioactive beams is estimated to determine whether it poses any significant risk to the patient compared to the use of a stable ion beam.
The remainder of this paper is organised as follows. A summary of key related work, including a description of the modified MKM which is adopted in this paper, is presented in Section 2. Details of the Monte Carlo simulations, including the phantom, physical and biophysical models used and the experimental validation of the selected physics models, the implementation of a pseudo-clinical beamline and treatment plan for stable and their corresponding radiactive ion beams are discussed in Section 3. Simulation results and analysis of the RBE 10 values of stable and radioactive beams, the resulting positron yield maps and the incidental dose resulting from the use of the radiactive beams are presented and discussed in Section 4. Conclusions and proposed future work are presented in Section 5.

Related Work
The use of positron-emitting radioisotopes for heavy ion therapy has been investigated by a number of authors. In 2001, Urakabe et al. demonstrated that a positron-emitting 11 C scanned spot beam could be directly used as the therapeutic agent 29 . However, the estimate of RBE 10 used to obtain a flat biological dose was based on an extrapolation of previously-reported results for 12 C in water, which was assumed to extend to human tissue 30 . Iseki et al. at NIRS used low-intensity monoenergetic 10 C probe beams with between 10 4 and 10 5 particles per spill to estimate the depth of the therapeutic 12 C beam's Bragg peak, while keeping the dose received during the range measurement under 100 mGyE (a few percent of therapeutic dose) 31 . RBE of the radioactive beam was estimated via simulation using the one-dimensional HIBRAC beam transportation code from Sihver et al. combined with Kanai's RBE model 30,32,33 . However, this work only considered monoenergetic 11 C ion beams, and ignored the effects of low-LET fragmentation products, which resulted in an overestimation of the RBE for 11 C. Augusto et al. used the FLUKA Monte Carlo toolkit to investigate the use of 11 C beams either alone or in conjunction with 12 C 34 . It was found that for beams with equivalent energy per nucleon incident on the same water phantom, 11 C and 12 C beams produce very similar fragmentation products, with the main differences being the relative yield of helium ions and several boron isotopes. While this study demonstrated the potential of using 11 C in heavy ion therapy, it only considered monoenergetic beams of 11 C at a fixed depth (100 mm) in a homogeneous water phantom. The www.nature.com/scientificreports www.nature.com/scientificreports/ composition of the phantom, the isotope and the specific beam energy are important factors affecting the fragmentation processs and the spatial distribution of positron-emitting nuclei which results 35,36 .
These works demonstrate the potential for using positron-emitting beams both for radiotherapy and for range verification. However, in order to conclusively establish their clinical utility, it is necessary to quantify their RBE and evaluate the quality of the resulting PET image in a clinically relevant configuration, through the use of heterogenous tissue-equivalent phantoms and polyenergetic ion beams.
Relative biological effectiveness (RBE) is an empirically-derived ratio which can be used to predict the physical dose of a specific type of radiation which will result in the same cellular survival fraction as a reference dose (typically a 200 keV X-ray beam) 37,38 . The complex dependencies of RBE on the energy and type of radiation, as well the location of the target and the specific tissue types present, require the use of biophysical methods for accurate theoretical estimation of RBE [39][40][41] . The Microdosimetric Kinetic Model (MKM), proposed by Hawkins et al., is a widely-used method for estimating RBE in which the microdosimetric spectrum (f(y)) is measured through the use of a tissue-equivalent proportional counter (TEPC) 24 . It was subsequently extended by Kase et al. to relate the saturation-corrected dose-mean lineal energy (y ⁎ ) to the radiation sensitivity coefficient α of the linear quadratic model (LQM, measured in units of Gy −1 and Gy −2 ), such that the method can be applied to therapeutic heavy ion beams 25,26,42 . This modified MKM has been extensively validated for carbon ion therapy, and also extended to proton and helium ion therapy 25,26,[42][43][44] .
The RBE 10 for an ion beam, defined as the ratio of the physical dose from a 200 kVp X-ray beam required to achieve a cellular survival fraction of 10% (D (10,R) ) to the ion beam dose resulting in the same cell survival fraction, can be derived using the microdosimetric spectra f(y), using (1), (2) and (3): For human salivary gland (HSG) tumour cells, the dose resulting in a survival fraction of 10%, D (10,R) is 5 Gy for 200 kVp X-rays; the LQM radiation sensitivity coefficient values are α 0 = 0.13 Gy −1 and β 0 = 0.05 Gy −2 . ρ and r d are the density and the radius of the sub-cellular domain, and assumed to be 0.42 μm and 1 g/cm 3 , respectively 25 .
In this work, RBE 10 is estimated using an extension to the modified MKM proposed by Bolst et al., whereby the mean path length l path < > of the charged particles that cross the sensitive volume is introduced to account for the directionality of the radiation field when deriving the microdosimetric spectra f(y) in a non-spherical sensitive volume, as opposed to the average chord length used in isotropic fields 27,28 .
Although estimates of the RBE 10 for radioactive beams have been reported previously, these have been calculated using simplified analytic models with parameters interpolated/extrapolated from limited experimental data from beams of stable isotopes in homogeneous targets 45,46 . The assumption that the RBE of radioactive ion species can be estimated from its stable analog has not been previously demonstrated in the literature.

Method
All Monte Carlo simulations were performed using the Geant4 toolkit (version 10.2.p03) 47,48 . The hadronic physics models used in the simulations are listed in Table 1, while electromagnetic interactions were modelled using the standard Geant4 option 3 physics constructor (G4EmStandardPhysics_option3). The hadronic physics processes and models are listed in Table 1. www.nature.com/scientificreports www.nature.com/scientificreports/ Section 3.1 details the methods used to experimentally validate the Geant4 simulation. The phantoms used in the simulations are described in detail in Section 3.2. In Sections 3.3, 3.4 and 3.5, the implementation of the modified MKM for the evaluation of the RBE 10 of pseudo-clinical, polyenergetic carbon and oxygen beams and their corresponding radioactive beams is described. Lastly, Section 3.6 describes a simulation study which examines the yield of different positron-emitting radionuclides during and after the irradiation of a skull phantom with radioactive and corresponding stable beams and introduces the metric used for the evaluation of the quality of the resulting annihilation maps.

Experimental validation of the physics models.
To validate the Monte Carlo physics models, several simulations evaluating depth-dose profiles and positron-emitting radionuclide yield were performed and compared with measurements obtained from equivalent physical experiments.
All experiments were performed at the Heavy Ion Medical Accelerator in Chiba (HIMAC), Japan, with the stable ion beams produced at the primary beam course, and the radioactive ion beams at the secondary beam course 19,23 . The peak energies of the non-radioactive 12 C and 16 O ion beams, as measured at the beamline nozzle, were 290 MeV/u and 400 MeV/u, respectively with an energy spread of σ = 0.2%. The peak energies of the radioactive 11 C and 15 O ion beams were 330 MeV/u and 290 MeV/u, respectively, each with an energy spread of σ = 5%. The nominal transverse diameter of all beams was 2 mm full width at half maximum (FWHM).
Depth-dose relationship. The experimental configuration shown in Fig. 1 was used to compare the experimental and simulation-based depth-dose curves. The deposited energy was measured using a pre-calibrated cross ionisation chamber (IC) with a sensitive volume of 36 mm 3 , inside a 300 × 300 × 300 mm 3 water phantom 23 . The IC was encased within a 0.5 mm PMMA casing and moved along the path of the beam using a motorised stage, with an accuracy of 10 μm. The energy deposited within the ionisation chamber at each point along the beam was normalised to the energy deposited at the entrance (i.e. at the front of the phantom). All depth measurements were converted to water equivalent depth.
For the simulation study, each beam was modelled using a monoenergetic incident beam with a Gaussian energy distribution, with the same peak energies and spreads as for the HIMAC beamlines. The simulated beams entered the water phantom perpendicular to its front surface (see Table 2), with an air gap of 2.5 m between the beamline nozzle and the phantom surface as per the corresponding experimental configuration. The energy deposited was scored in the water phantom using 1 mm 3 voxels and summed over a 36 mm 3 volume equivalent to the sensitive volume of the ionisation chamber used throughout the experimental measurements. Energy deposited in the sensitive volume (as a function of depth) was normalised to value observed at the entrance plateau.
Positron-emitting fragmentation product yield. The hadronic physics models of Geant4, including the Quantum Molecular Dynamics (QMD) ion hadronic inelastic scattering and Radioactive Decay physics models, were validated by comparing the simulated and experimentally estimated yields of 11 C, 10 C and 15 O, the three dominant   www.nature.com/scientificreports www.nature.com/scientificreports/ positron-emitting radionuclides generated during irradiation of a 100 × 100 × 300 mm 3 PMMA phantom by monoenergetic 12 C and 16 O beams with energies of 290 MeV/u and 400 MeV/u, respectively.
The experimental configuration is shown in Fig. 2. The phantom was positioned such that the expected location of the Bragg peak was aligned with the centre of the field of view in the OpenPET scanner 14 . 20 spills were used, with a beam intensity of 1.0 × 10 9 particles per second (pps). In each spill, the beam was on for 1.9 seconds and off for 1.4 seconds. List-mode PET data were collected intra-spill, and for 36 minutes after the final spill. Dynamic (4D) images were reconstructed using the 3D ordinary Poisson ordered-subset-expectation-maximisation algorithm (3D-OP-OSEM) with 1.5 × 1.5 × 1.5 mm 3 voxels. Temporal frame lengths were chosen so as to be able to observe decay over several half-lives of 11 C, 10 C and 15 O. Yields of each positron-emitting radionuclide were estimated by fitting the parameters of a simple analytical model to the observed time-activity curves (TACs). Total activity as a function of time t in a volume with initial activities of 11 C, 10 C and 15 O of A 0,C11 , A 0,C10 and A 0,O15 , respectively, is given by 10 15 where T C11 , T C10 and T O15 are the respective half-lives of 11 C, 10 C and 15 O. Total activity is measured as a function of time across the build-up and Bragg peak region, defined as the region from the point at which the dose profile has risen 5% above the entrance plateau to the point after which the profile is below 5% of the peak value. The individual initial activities for each radionuclide are then estimated for both the simulation results and the experimental data by fitting the model to the observed curve.
For the simulation studies, monoenergetic 12 C and 16 O beams were directed perpendicularly to the surface of a simulated 10 × 10 × 30 cm 3 PMMA with an air gap of 1.75 m between the beamline nozzle and the phantom surface, matching the experimental configuration. Density, mean excitation, ionisation potential and dimensions of simulated phantoms were chosen to match those used in the experiment. The spatio-temporal distributions of positron-emitting nuclei, positron production and annihilation were recorded with a scoring volume resolution of 1.5 mm, corresponding to the voxel dimensions in the experimental PET images. Simulated yield profiles were convolved with a Gaussian filter, with its FWHM equal to the estimated OpenPET spatial resolution (3.5 mm) 14 .
Phantom geometry. The phantoms used in the simulation were rectangular prisms with compositions as listed in Table 2. All material compositions were based on data from the National Institute of Standards and Technology (NIST) database 49 . estimation of RBe 10 for a pseudo-clinical SOBP. To evaluate and compare the RBE 10 of polyenergetic stable and positron-emitting radioactive beams, sensitive 1 mm × 1 mm × 10 μm volumes were defined every 100 μm along the path of the beam. The lineal energy deposition spectrum in each volume for all interactions (f(y)) was stored and used to calculate the RBE 10 at that point, using (3) (equivalent results for monoenergetic carbon and oxygen ion beams with the energies listed in Supplementary Table S1 are presented in Supplementary  Table S2). A correction factor 1.05 were used to account for the difference in stopping power and density of water relative to brain tissue.
A simple variance analysis method was used to estimate a sufficient number of primary particles to use in the simulations. M test simulations were conducted, each with N primary particles, with RBE estimated for each simulation and the mean and standard deviation (SD) calculated across the M simulations. The standard deviation should approach zero as N tends to infinity; therefore, in this experiment, N was progressively doubled with a fixed value of M = 50 until the ratio of standard deviation to mean was less than an arbitrary threshold of 1%. This analysis suggested that N = 10 7 would be sufficient to get a good estimate of RBE (95% probability of the estimated RBE being within ±2% of the true RBE).
Carbon. The spectrum for the simulated carbon beams was generated using an experimentally-validated model of the passively-scattered 12 C beamline at HIMAC, which is known to produce a flat biological dose across a 60 mm depth range 50 . The spectra of the positron-emitting radioactive beams ( 10 C and 11 C) were based on the www.nature.com/scientificreports www.nature.com/scientificreports/ 12 C spectrum from this beamline, by determining the energies for which the Bragg peaks of monoenergetic radioactive ion beams where located at the proximal and distal edges of the desired SOBP, and linearly mapping the weights of the energies of the 12 C SOBP spectrum to this range of energies. Finally, the SOBPs were compared and confirmed to both correspond to the planned depths.
Oxygen. Currently a validated model of the 16 O beamline does not exist. Therefore, generation of the 60 mm flat biological dose SOBP in the target depth range was achieved by performing monoenergetic Monte Carlo simulations of an 16 O beam at a range of energies (177, 237, 297, 345 and 418 MeV/u), and evaluating the RBE 10 as a function of depth for each energy using the modified MKM (see Section 3.3 and Supplementary Table S2). This RBE was used to convert the physical dose deposited in the simulations to an estimated biological dose for the 5 evaluated energies. Profiles were then generated for other intermediate energies by interpolating between the simulated values in increments of 1 MeV/u. Finally, the target flat biological dose was achieved by adjusting the weights of each of these profiles such that a flat biological dose rate of 5 Gy(RBE)/min was achieved within the target depth range. The spectra of the positron-emitting radioactive beam ( 15 O) was based on the 16 O spectrum, with energies scaled such that the SOBP was positioned in the desired depth range (as per carbon).
Positron-emitting radionuclide yield study. The impact of using positron-emitting primary beams on interspill and post-irradiation image quality was evaluated by comparing the spatial distributions of positron decays observed in the simulation over several different intervals during treatment of the skull phantom. A simple treatment plan was designed for each primary particle type, aimed at producing a constant biological dose rate of ≈5 Gy(RBE)/min in a depth range of 78-138 mm within a skull phantom. A total of 1 × 10 9 primary particles were used in each simulation. As for the experimental validation study, twenty spills were simulated, with the beam on for 1.9 seconds and off for 1.4 seconds.
The distributions of positron decays were acquired for each beam type between the first and second spill, during the first five inter-spill intervals, and in the five minutes following the final spill.
The contrast-to-noise ratios (CNRs) between the inside and outside of the proximal, distal and upper lateral edges of the SOBP are computed for each image. The CNR provides a metric for objectively comparing the specificity with which the irradiated region is delineated, and is defined as: where μ a and μ b are the mean signal amplitudes and σ a and σ b are the standard deviations of the image intensity in two regions a and b of the image 51,52 .

Results and Discussion
Physics model validation. Experimental and simulation-based depth-dose curves are shown in Fig. 3; the difference between the locations of the Bragg peaks obtained from the simulated and the experimental 12 C, 16 O, 11 C and 15 O depth-dose profiles were 0.8 mm, 0.24 mm, 0.37 mm and 0.43 mm, respectively. The experimental and simulation yields of 10 C, 11 C, and 15 O produced during the irradiation of a PMMA phantom using a 290 MeV/u 12 C beam and a 400 MeV/u 16 O beam are expressed as a percentage of the total positron-emitting radionuclide yield and are listed in Table 3.
The close agreement between the experimental and simulated normalised depth-dose profiles and the relative yield estimations of the positron-emitting fragments demonstrate the validity of the simulation model. The small differences between the experimental and simulated depth dose profiles for radioactive primary particles may be due to an underestimation of the initial energy spread, heterogeneity of the beryllium target leading to contamination with other fragments and systematic errors introduced by the ionisation chamber measurements. Figure 4 presents a comparison of RBE 10 as a function of depth for the positron-emitting radioactive beams and for the corresponding stable isotope beams (for clarity, RBE 10 values are shown at depth increments of 3 mm; refer to Supplementary Spreadsheet 1 for a full list of RBE 10 values evaluated at 100 μm intervals for all ion species). In each case, the mean RBE 10 s of the stable and radioactive beams are well within each others' 95% confidence interval. Radioactive-to-stable RBE 10 ratios are also shown, with the mean values remaining very close to 1.0 in the entrance and SOBP. The larger confidence intervals in the tail region are due to very little energy being deposited beyond the end of the SOBP (as expected for heavy ion beams), resulting in significant statistical noise. Table 4 lists numerical values of the RBE 10 obtained at the entrance, the beginning, middle and end of the SOBP, and tail region for each beam type. The mean and standard deviations presented are calculated over 11 consecutive 100 μm deep sample volumes along the beam path centred about the listed depth. In all cases, the mean RBE 10 s for the radioactive and corresponding stable ion beams are within one standard deviation of each other.

RBE and biological dose in Gy(RBE).
The significance of this result is that it indicates that the evaluated radioactive ion beams are comparable to their non-radioactive counterparts in terms of relative biological effectiveness. Heavy ion therapy with any of the radioactive ion species examined in this study should be feasible, with only minimal changes to the current treatment planning algorithms required to account for the small differences in RBE 10 .
Positron yield. Figure 5 shows the 2D annihilation maps obtained during and after the simulated delivery of 5 Gy(RBE) for each beam type to the target volume within the skull phantom. Images in the first column correspond to data acquired during the first beam-off interval (i.e. after one spill), the centre column show images www.nature.com/scientificreports www.nature.com/scientificreports/ following 5 spills, and finally, the last column shows images acquired during the five minutes (300 seconds) immediately after the completion of the 20th (and final) spill.
The CNRs of the inside and outside of the proximal, distal and upper lateral boundaries of the SOBP images in Fig. 5 are listed in Table 5.
Positron annihilation maps acquired at different stages of the treatment process clearly demonstrate the potential improvements in range-verification QA that can be obtained with radioactive ion beams. Following a single spill, the boundaries of the SOBP are very clearly visible in the cases of 10 C and 15 O (Fig. 5(a,j)), due to their short half-lives (19.29 seconds and 122.24 seconds, respectively). The images from the 10 C simulation also exhibit the the highest CNR values for all boundaries after both 1 and 5 spills (i.e. the delivery of 5% and 25% of the total planned dose) and 5 minutes after the delivery of the full treatment for the distal boundary. 15 O also exhibits an excellent CNR following a 5 minute acquisition, demonstrating the best results for proximal and distal edge. The two stable beams produce images which are indistinct in comparison to any of the radioactive beam images. Due to its half life of 20.334 minutes, only a small number of positron annihilations resulting from decays of 11 C are observed within the first beam-off period (Fig. 5(d)). The distal edge can be clearly seen, however the proximal edge is indistinct. Finally, in the long post-irradiation image acquisition (right column in Fig. 5), most primaries from the 10 C and 15 O beams have decayed, resulting in very similar high-contrast images. A substantial number of primaries have now decayed in the case of 11 C, resulting in the emergence of a well-defined edges to the SOBP; it is expected that a 11 C beam with a post-irradiation image acquisition of 20 minute or more will result in very   Table 3. Relative yields of positron-emitting nuclei in experiment and simulation.
www.nature.com/scientificreports www.nature.com/scientificreports/ high CNRs due to its longer half-life. By contrast, after a 5 minute acquisition, the distal and proximal edges of the SOBP remain indistinct in the case of 12 C. 16 O exhibits a more well-defined distal edge to its SOBP compared to 12 C, however, the proximal edge is again poorly defined.
The images also demonstrate one of the key differences between the radioactive and stable beams. For radioactive beams, positron annihilations principally occur in the vicinity of the stopping point of the primary particle. The intensity of the decay radiation observed in a PET image is therefore proportional to the number of primary particles which have arrived at that particular depth. The energy weightings required to achieve a flat biological dose have a bias towards higher energies (since more deeply-penetrating high-energy particles also deposit an entrance dose which is added to the dose deposited by lower energy beams). Therefore, the distal edge of the SOBP can be expected to be much brighter than the proximal edge, as is clearly evident in the images from the radioactive beams. By contrast, the contribution of primary or target fragmentation, which is relatively minor for the radioactive beams, is the only source of positrons in the case of the stable beams, and positron-emitting fragmentation products are produced to a varying extent along the entire length of the beam path (see Supplementary  Tables S3 and S4). Therefore, the stable beams exhibit a flatter (although not completely flat) activity distribution in the SOBP, and weaker contrast between the SOBP and the entrance region.
Radiation dose to patients. Given the superiority of positron-emitting radioactive beams for intra-and post-treatment QA imaging, it is also important to consider whether or not the use of such beams would have any unintended side effects for the patient. From this perspective, the main difference for the patient is that an additional radiation dose will result from the use of a radioactive beam. The dose resulting from the decay of a positron-emitting radionuclide includes the kinetic energy of the positrons together with the 511 keV gamma photons resulting from their eventual annihilation; for a 11 C beam, a 70 Gy(RBE) dose delivered to a 100 mm cubic treatment volume would require approximately 2.3 × 10 11 particles, distributed throughout the treatment volume. This corresponds to an initial activity concentration of 1.3 MBq/cc, which is comparable to tissue concentrations of radiotracer which would be used in diagnostic 11 C clinical PET imaging, and would deliver a biological dose' within the treatment volume of the order of 3-10 mSv. The additional dose rapidly falls off outside the treatment volume, and would be insignificant compared to the dose due to lateral scattering of particles. . Biological dose, physical dose and RBE 10 for positron-emitting radioactive beams, together with the ratio of radioactive-beam RBE 10 to stable-beam RBE 10 , each shown as a function of depth within the phantom. The objective is a uniform dose within a 60 mm SOBP, from 78 to 138 mm depth. For carbon, 12 C is shown in red, with 11 C and the ratio of 11 C: 12 C shown in green, and 10 C and the ratio of 10 C: 12 C shown in blue. 16 Table 4. Means and standard deviations of the RBE 10 for each beam evaluated at five depths (entrance, start, middle and end of SOBP, and tail). At each depth, RBE 10 is evaluated in 11 adjacent sensitive volumes (every 100 μm along the path of the beam) and the mean and standard deviation calculated.

Conclusion
This work aimed to quantitatively evaluate the therapeutic potential of positron-emitting radioactive heavy ion beams; in particular, with regard to the relative biological effectiveness of the beams compared to their non-radioactive counterparts, the spatial distribution of the positron-emitting annihilations generated during and after irradiation of the target, and the incidental dose to the patient. Monte Carlo simulations of heavy ion therapy using a pseudo-clinical spread out Bragg peak constructed with positron-emitting radioactive beams of 11 C, 10 C and 15 O as well as stable 12 C and 16 O were undertaken with the Geant4 toolkit. The simulation physics model was validated through a comparison of depth-dose curves for monoenergetic 11 C, 12 C, 15 O and 16 O beams and relative yield estimations of the positron-emitting fragments produced within the build-up and the Bragg peak region with experimental data for 12 C and 16 O obtained from the HIMAC facility in Japan. The maximum difference between the location of maximum dose in the simulation and experimental data was 0.8 mm, while the maximum difference in mean relative yields of the secondary positron-emitting fragments was 2%.
The radiobiological effectiveness (RBE 10 ) of each beam was calculated for an SOBP extending from depths of 78 to 138 mm in a skull phantom using the modified microdosimetric kinetic model (MKM). The RBE 10 of the radioactive ion beams was found to be within one standard deviation of the corresponding non-radioactive ion beams for all energies, indicating that the therapeutic efficacy of such beams should be very similar to beams of the corresponding non-radioactive ion.
Finally, the additional dose to the patient resulting from the use of radioactive beams was estimated to determine whether it poses any unreasonable risk to the patient compared to the use of a stable ion beam. The additional dose was found to be comparable to that received during diagnostic clinical PET, and therefore negligible compared to the dose delivered to the target volume or surrounding tissues during the radiotherapy procedure.
In summary, positron-emitting radioactive heavy ions are approximately equivalent to the corresponding stable isotope with respect to expected therapeutic properties in heavy ion radiotherapy, while being greatly superior to non-radioactive beams in terms of the potential for accurately imaging the treatment volume during and after treatment. The substantial increase in positron yield offered by positron-emitting radioactive beams for the same biological effective dose will allow the boundaries of the spread out Bragg peak in a PET image to be unambiguously identified, making the use of positron-emitting radioactive ions a compelling choice for heavy ion therapy.

Data Availability
All data generated or analysed during this study are included in this published article (and its Supplementary Information Files) or are available from the corresponding author on reasonable request.