A quantitative assessment of Geant4 for predicting the yield and distribution of positron-emitting fragments in ion beam therapy

Purpose: To compare the accuracy with which different hadronic inelastic physics models across ten Geant4 Monte Carlo simulation toolkit versions can predict positron-emitting fragments produced along the beam path during carbon and oxygen ion therapy. Materials and Methods: Phantoms of polyethylene, gelatin or poly(methyl methacrylate) were irradiated with monoenergetic carbon and oxygen ion beams. Post-irradiation, 4D PET images were acquired and parent $^{11}$C, $^{10}$C and $^{15}$O radionuclides contributions in each voxel were determined from the extracted time activity curves. Experiments were simulated in Geant4 Monte Carlo versions 10-11.1, with three different fragmentation models: binary ion cascade (BIC), quantum molecular dynamics (QMD) and the Liege intranuclear cascade (INCL++) - 30 combinations. Total/parent isotope positron annihilation yields were compared between simulations and experiments using normalised mean squared error and Pearson cross-correlation coefficient. Depth of maximum/distal 50\% peak position yield were also compared. Results: Performance varied considerably across versions and models, with no one best predicting all positron-emitting fragments. BIC in Geant4 10.2 provided the best overall agreement with experimental results in the largest number of test cases. QMD consistently provided the best estimates of both the depth of peak positron yield (10.4 and 10.6) and the distal 50\%-of-peak point (10.2), while BIC also performed well and INCL generally performed the worst across most Geant4 versions. Conclusions: Best spatial prediction of annihilation yield and positron-emitting fragment production during carbon and oxygen ion therapy was found to be 10.2.p03 with BIC or QMD. These version/model combinations are recommended for future heavy ion therapy research.


Introduction
One of the chief advantages of particle therapy as a treatment for cancer is the high dose gradient between the treatment area and surrounding regions [1].However, to enable the smallest possible margins, particle range must be accurately verified relative to the patient anatomy.Many authors have investigated the post-irradiation imaging of positron-emitting fragments via positron emission tomography (PET) for range verification of the primary beam during particle therapy [2,3,4,5,6,7,8,9], with several prototype in-beam PET scanners currently under active development [10,11,12,13].Monte Carlo simulation methods, implemented using platforms such as Geant4, are now widely used in particle therapy -applied clinically as a secondary quality assurance (QA) technique to verify dose during particle therapy [14,15,16] and there is active research in using it for developing and testing range models [17,18,19], PET image-to-dose models [20,21,22,23], and designing novel PET scanners optimised for range and dose verification [24,25].
Good models for nuclear fragmentation processes are essential for faithfully simulating imaging applications in particle therapy, such as PET-based dose estimation methods for quality assurance, since the production and distribution of positron-emitting radionuclide fragments directly affects the quality of the resulting image [20,21,22,26].Currently, the Geant4 Monte Carlo simulation framework offers a choice of three hadronic inelastic fragmentation models that are appropriate for particle therapybinary ion cascade (BIC), quantum molecular dynamics (QMD), and Liège intranuclear cascade (INCL++) [14,27,28] ‡.In our previous study, we evaluated these models by comparing the spatial distributions of positron-emitting radionuclides predicted following irradiation of PMMA, gelatin and polyethylene targets by monoenergetic carbon and oxygen ion beams (simulated using Geant4 10.2.p03) to equivalent results estimated from experimentally-obtained PET data [30].The BIC model provided the best estimates overall; however, none of the models provided a perfect fit in all evaluated cases, and some significant discrepancies were observed.Other authors have noted similar variations in simulation fidelity in different contexts; most recently, Wrońska et al. [31] evaluated two models, BIC and QBBC, with several versions of Geant4 for prediction of prompt gamma emission during proton therapy with a PMMA target.They concluded that Geant4 10.4p02's BIC implementation with high-precision neutron data libraries and binary ion fragmentation models provided the best agreement with the observed experimental data for this use case, with the simulations performed using other version/model combinations producing features which were absent in the experimental data.Pei et al. [7] evaluated the ability of Geant4 to predict the range of particles using Geant4 version 10.6 and the Bertini fragmentation model.They concluded that Geant4 was able to accurately predict within 1 mm the positron activity peak and distal falling edge using two carbon monoenergetic beams and a single PMMA phantom.
Since the publication of our previous study, there have been several updates to Geant4; specifically, six minor releases (versions 10.x) and one major release (version 11, which has since been updated to version 11.1).Each of these releases includes modifications to the physics models implemented in Geant4, which can affect the simulation of positronemitting fragment production in particle therapy.
In addition to the normalised mean squared error (NMSE) metric used in the previous study, three additional metrics -the Pearson cross-correlation coefficient (CCC), the depth of the positron annihilation peak, and the depth at which positron annihilation ‡ INCL++ is considered the most appropriate option for neutron spallation simulations, but is included here for completeness [29].
intensity has decreased to 50% of peak are also used to compare the shape of the predicted positron-emitting fragment distribution to the experimentally-measured distribution.The experimental absolute yields of the dominant positron-emitting fragment species ( 10 C, 11 C and 15 O) created during irradiation of gelatin, PMMA and polyethylene block phantoms with beams of 12 C ions with energies of 148.5, 290.5 and 350 MeV/u and 16 O ions with energies of 148 and 290 MeV/u were obtained across the full width at tenth maximum (FWTM) of the beam from a previous study conducted using the OpenPET scanner at HIMAC, Japan [12,30].Akamatsu2019 results in each of the three target materials and 5 ion/energy combinations were then compared to those predicted in equivalent simulations performed in Geant4 using each of hadronic fragmentation models (BIC, QMD and INCL++) across the ten evaluated Geant4 versions for a total of 150 unique target/ion/energy/version/model test conditions.The total positron yields and yields of the individual positron-emitting fragment species from each model and Geant4 version were then compared to experimental measurements using each of the chosen metrics at the entrance, build-up and Bragg peak, and tail regions.

Materials and Methods
This section describes the methods used for obtaining and quantitatively comparing the experimental and simulated positron annihilation profiles.The general approach is similar to that used in our previous study (see [30]); however, it has been extended to include a much wider range of Geant4 versions, and additional comparison metrics are introduced.
The experimental methods used to estimate the total positron annihilation profile and activity of the dominant positron-emitting fragment isotopes ( 11 C, 10 C and 15 O) are briefly summarised in Section 2.1.Equivalent simulation configurations were constructed for each Geant4 version under test, and the total positron annihilation profile and activity of 11 C, 10 C and 15 O were predicted for each beam ion/energy, target material, hadronic inelastic fragmentation model and Geant4 version; the design and parameters of these simulations are described in detail in Section 2.2.
The simulated annihilation profiles are compared with the experimental annihilation profiles using the following metrics in each of the entrance, build-up and Bragg peak and tail regions: • Normalised mean squared error (NMSE); and Additionally, the depth of positron annihilation peak and the depth of the distal point at which the magnitude of the positron annihilation profile decreases to 50% of the peak value are evaluated.All metrics are described in detail in Section 2.3.

Experimental configuration
The experimental data obtained in our 2019 paper were used as the ground truth for this simulation study; a detailed description of the experimental procedures is presented in that paper [30].In summary, phantoms constructed from either pure PMMA, polyethylene or gelatin (encased in a thin-walled PMMA container), each with dimensions of 100 mm×100 mm×300 mm, were irradiated with monoenergetic carbon or oxygen ion beams of various energies -three for carbon ions and two for oxygen (see Table 1).Positron annihilation profiles (with respect to depth in the target) were obtained using the whole-body DOI-PET scanner prototype developed at QST [12].These profiles were decomposed into the individual population of each of the dominant parent positron-emitting fragments ( 11 C, 10 C and 15 O) at t = 0 (end of irradiation period) by fitting the observed time-decay curves in each voxel to a multiexponential decay model.
For each version of Geant4, three alternative hadronic ion fragmentation models were evaluated -Binary Ion Cascade (BIC), Quantum Molecular Dynamics (QMD) and Liège Intranuclear Cascade (INCL) models § [27,28].All simulations modelled electromagnetic § The INCL model was developed specifically for spallation reactions but is included in this study as it can also model fragmentation.
interactions using the standard option 3 list (G4EmStandardPhysics_option3).The remaining physics processes, including hadronic physics models, are listed in Table 10.
The location of each positron annihilation, as well as the identity of the parent isotope which decayed to emit each positron (principally 11 C, 10 C and 15 O), were scored with a resolution of 1.5 mm 3 to match the voxel dimensions of the experimental OpenPET image reconstruction output.The pristine positron annihilation profiles were convolved with a 2.3 mm FWHM Gaussian filter to simulate the point spread function of the PET system.

Evaluation methods
The irradiated target was divided into three separate regions for analysis since different physics processes dominate in each: the entrance region, the build-up and Bragg peak region, and the tail region.This segmentation is defined in the same way as in our previous paper [30]; in summary, the central build-up and Bragg peak region is defined as follows: • The proximal edge in the z dimension (along the path of the beam) is defined as the first point at which the dose deposited along the central axis exceeds the entrance plateau dose by more than 5% of the difference between peak dose and the entrance plateau dose; and • The distal edge in z is defined as the last point at which the deposited dose is greater than 5% of the absolute peak dose value.
The entrance region is then defined as the region proximal to the build-up and Bragg peak region, while the tail region is defined as the region distal to the build-up and Bragg peak region.
The yields of positron-emitting nuclei are defined via (1): where N (Isotope) is the yield of the isotope under study in that region and N (Primary) is the total number of primary particles.Yields were calculated in each voxel along the path of the beam.
Three different metrics were chosen to quantify the accuracy of each model in Geant4: the normalised mean squared error (NMSE), the Pearson cross-correlation coefficient (CCC), and the range (depth along the path of the beam) of both the position annihilation peak and the point beyond the peak at which positron annihilation intensity decreases by 50%.The NMSE was used to measure the absolute difference between the experimental measurements and simulation-predicted positron yields in each region: where S i and E i are the simulation and experimental yields in the ith voxel of the N reg voxels in region reg (with a lower value indicating a better match).
The Pearson CCC was used to compare the shapes of the simulation-predicted positronemitting fragment distributions to the experimental measurements: where S norm,i and E norm,i are the normalised simulation and experimental yields in the ith voxel of the N reg voxels in region reg.Normalisation is performed by dividing each S i and E i by the maximum value in its respective region.S norm and E norm are the mean values in each region.
When comparing the models, the closer the CCC between the model prediction and the experimental estimate of positron-emitting fragment distribution is to +1, the more accurate the prediction.A CCC greater than 0.95 is considered to indicate that the simulation has accurately predicted the shape of the yield.
For each version of Geant4, phantom, beam type and energy, the NMSE and CCC were calculated for both total annihilation photon yield profiles and also for the profiles of the three main positron-emitting fragment species ( 10 C, 11 C and 15 O).The calculation was repeated for each of the N reg regions (entrance, build-up and Bragg peak, and tail regions).The NMSE and the CCC were then compared across all evaluated Geant4 versions for each region, phantom material and beam type.
A total of 5 energy/ion combinations are evaluated (carbon ions at three energies and oxygen ions at two energies).For oxygen ions, three target materials (gelatine, PMMA and polyethylene) are evaluated for total positron annihilation yield and 11 C/ 10 C/ 15 O yield.For carbon ions, the same three target materials are evaluated for total positron annihilation yield and 11 C/ 10 C yield and two for 15 O yield (polyethylene is omitted since it is not possible to produce 15 O fragments with a 12 C ion beam and a PE target which only contains carbon and hydrogen).Thus, there are a total of 15 cases evaluated for total positron annihilation, 11 C yield and 10 C yield, while there are 12 cases evaluated for 15 O yield.
For range calculations, the difference between the depths at which the positron annihilation yield reached its maximum value in the experiment and simulation was calculated (see ( 4)).Additionally, the point distal to this maximum at which positron annihilation yield decreases to 50% of the maximum value was also compared between experiment and simulation.For each version and model, the mean differences between the experimental and simulation-based values, as well as the standard deviations and maximum differences were calculated across all test cases (ion species, energies and target materials).
where R x is the range (depth) of the voxel with the maximum value (or, for distal 50%of-peak, the first distal voxel to fall below 50% of the maximum value) in either the simulation or experiment.

Results and Discussion
The number of cases in which each version/model combination performed the best or equal-best in terms of each of the evaluated metrics are counted across all simulations in the entrance, build-up and Bragg peak and tail regions, and summarised in this section.Detailed results for each experiment are included in the Supplementary Materials.
For the NMSE metric, equal best is defined by firstly identifying the version/model with the best mean NMSE, and then accepting any additional version/model combination whose mean NMSE is within the confidence interval of the best-performing version/model as being equal.All confidence intervals are ±2σ (i.e.95% assuming a Gaussian distribution).

Entrance region
In the entrance region, positron-emitting fragments are created by target fragmentation rather than projectile fragmentation.The projectile ions lose energy via Coulomb interactions, slowing down at an approximately constant rate as they traverse this region, with only gradual changes to projectile/target cross sections.As a result, the positronemitting fragment distributions are expected to exhibit an approximately flat depthwise profile in this region.
NMSE and Pearson CCC results between simulation and experimental total positron annihilation profiles in the entrance region are summarised in Tables 2 and 3, respectively, with corresponding figures shown in Supplementary Material Section 1.
For the entrance region, the BIC model implemented in Geant4 versions 10, 10.1, 10.

Build-up and Bragg peak region
In the build-up and Bragg peak region, positron-emitting fragments are produced via a combination of target fragmentation and projectile fragmentation.There is a rapid change in positron-emitting fragment yield with respect to depth, especially since different positron-emitting fragments stop at different distances from their point of production.
NMSE and Pearson CCC results between simulation and experimental total positron annihilation profiles in the build-up and Bragg peak region are summarised in Tables 4  and 5, respectively, with corresponding figures shown in Supplementary Material Section 2.
In the build-up and Bragg peak region, according to the NMSE metric, total positron yield is most accurately predicted by the BIC model in Geant4 version 10.Table 6 lists difference between the experimental and simulation positron peak, while Table 7 lists the difference between the 50% fall off point for the experimental and simulated positron peak.

Tail region
In the tail region, positron-emitting radionuclides are primarily produced through fragmentation of the target material caused by light fragments created upstream from the primary beam.As such, the production of positron-emitting fragments in the tail region is highly dependent on fragmentation and scattering cross sections upstream.Therefore, the yield of positron annihilation is not expected to rapidly change across this region compared to the build-up and Bragg peak region.
NMSE and Pearson CCC results between simulation and experimental total positron annihilation profiles in the tail region are summarised in Tables 8 and 9, respectively, with corresponding figures shown in Supplementary Material Section 3.
Using the NMSE metric, 10.2/BIC was the best-performing version/model combination for overall positron yield (12/15 cases), with 10.2/QMD being the second-best (10/15).Using the depth-of-maximum-yield metric, the smallest mean differences were obtained with 10.4/QMD and 10.6/QMD.These versions/models also achieved the equal-smallest maximum difference (-3 mm and +3 mm, respectively).Across all versions of Geant4, QMD demonstrated the best overall accuracy (lowest average mean difference in peak depth) and highest precision (lowest average standard deviation).INCL was the worstperforming model across all versions, with much larger maximum differences, and a consistent underestimation of depth of maximum yield across Geant4 versions, with the exception of version 11.1 (which, despite a mean difference of 0, exhibited a large standard deviation and maximum value).Standard deviations obtained using INCL were generally around double those of QMD and BIC.BIC also showed a consistent underestimation in depth of maximum yield, although the maximum differences were much smaller than for INCL.For context, the difference between the depth of the positron annihilation peak and the Bragg peak with monoenergetic ion beams is of the order of -5.6±0.8 mm for Not all models met or exceeded the set threshold of 0.95 for the Pearson CCC metric.This means that in these cases, the shape of the predicted positron distribution differs significantly from the experimental measurements.This is of particular concern if these models are to be used for dose estimation using a deconvolution approach [20,21] or for the training of machine learning models for feature extraction [23].
Finally, it is worth noting that current evaluations of fragmentation cross sections exhibit uncertainties exceeding 10%, which must be addressed in order to accurately model positron fragmentation, particularly in the case of complex fragmentation reactions such as the production of 10 C [34,35].This also impacts other Monte Carlo simulation platforms (such as FLUKA, MCNP and PHITS) which rely on accurate cross section data.

Conclusion
In this study, the accuracy with which Geant4 is able to predict the distribution of total positron annihilation yield and the distributions of individual positronemitting fragmentation products ( 11 C, 10 C and 15 O) during carbon or oxygen ion therapy was compared to experimental data.Three different hadronic inelastic physics models -Binary Ion Cascade (BIC), Quantum Molecular Dynamics (QMD) and Liege Intranuclear Cascade model (INCL) were used with ten different versions of Geant4 -10.0.p04, 10.1.p03,10.2.p03, 10.3.p03,10.4.p03, 10.5.p01, 10.6.p03,10.7.p02, 11.0 and 11.1, in three different homogeneous phantoms.The simulated and experimental data were compared using two different metrics -normalised mean squared error and the Pearson cross-correlation coefficient.Additionally, the differences between the simulated and experimental depth of maximum positron annihilation yield, as well as the distal point at which positron yield declines to 50% of the peak were evaluated.It was found that the accuracy of the hadronic inelastic physics models strongly depends on the version of Geant4 in which it was implemented, and newer versions of Geant4 were not always more accurate at predicting positron-emitting fragmentation compared to older versions.Furthermore, it was found that not all version/model combinations were able to satisfactorily predict the shape of positron annihilation or positron-emitting fragment distributions, even though they could provide a good estimation of the total positron annihilation yield and range.For future simulation studies of therapeutic irradiation using carbon or oxygen ion beams, it is recommended that Geant4 version 10.2 with the BIC model be used as it is currently the version/model combination best able to replicate the experimentally-observed total positron yield and the fragmentation product distributions, while the depth of the maximum positron yield and distal 50%-of-peak point were best predicted using the QMD model from Geant4 10.4, 10.6 (peak) and 10.2 (distal 50%-of-peak).

Acknowledgements
3 and 10.4 provided the (equal) lowest NMSE of the yields of total positron annihilation in 5 out of 15 cases.The BIC model in Geant4 10, 10.1 and 10.2 also provided the (equal) lowest NMSE for 11 C fragment production (11/15 cases), whereas for 10 C the

Table 1 .
Beam parameters for each ion species and energy.The energy spread is 0.2 % of nominal energy in each case; 95% confidence intervals are given for beam flux.

Table 2 .
Number of test cases for which each Geant4 version/model combination achieved the lowest or equal-lowest NMSE in the entrance region.Bold text denotes the version/model achieving the highest (or equal-highest) number of best results for each combination of ion/energy/target.

Table 3 .
Number of test cases for which each Geant4 version/model combination achieved a CCC greater than 0.95 in the entrance region.Bold text denotes the version/model achieving the highest number of best results for each combination of ion/energy/target.
Results for individual radionuclides were also mixed, with 10/BIC, 10.1/BIC, 10.4/BIC and 10.5/INCL achieving the threshold in 4/15 cases for 11 C, 10-10.4/BICand all versions with QMD reaching the threshold in 2/15 cases for 10 C, and all versions with BIC and 10/INCL, 10.1/INCL, 10.2/INCL, and 10.5/INCL reaching the threshold for

Table 4 .
Number of test cases for which each Geant4 version/model combination achieved the lowest or equal-lowest NMSE in the build-up and Bragg peak region.Bold text denotes the version/model achieving the highest (or equal-highest) number of best results for each combination of ion/energy/target.

Table 5 .
Number of test cases for which each Geant4 version/model combination achieved a CCC greater than 0.95 in the build-up and Bragg peak region.Bold text denotes the version/model achieving the highest number of best results for each combination of ion/energy/target.

Table 6 .
Differences between the depths of the maximum positron annihilation yield in experimental and simulation results.Each voxel has a width of 1.5 mm; the maximum error is always in multiples of 1.5 mm increments.

Table 7 .
Differences between the distal depths at which the positron annihilation yield has decreased to 50% of the peak value in experimental and simulation results.Each voxel has a width of 1.5 mm; the maximum error is always in multiples of 1.5 mm increments.

Table 8 .
Number of test cases for which each Geant4 version/model combination achieved the lowest or equal-lowest NMSE in the tail region.Bold text denotes the version/model achieving the highest (or equal-highest) number of best results for each combination of ion/energy/target.

Table 9 .
Number of test cases for which each Geant4 version/model combination achieved a CCC greater than 0.95 in the tail region.Bold text denotes the version/model achieving the highest number of best results for each combination of ion/energy/target.For total positron annihilation yield, 10/BIC, 10.1/BIC, 10.2/BIC, 10.6/INCL, 10.7/BIC, 10.7/INCL, 11/BIC, 11/INCL, 11.1/BIC and 11.1/INCL all exceeded the target threshold for 12/15 cases.Even the worst-performing version/model combinations still exceeded the threshold in 10/15 cases.For 11 C yield, 11/INCL and 11.1/INCL reached the threshold in 13/15 cases (with the worst-performing combination scoring 9/15 wins).Fewer wins were seen with 10 C; the best results were obtained with 10/BIC, 10/INCL, 10.1/BIC, 10.1/QMD, 10.1/INCL, 10.2/BIC, 10.2/QMD and 10.2/INCL (4/15 cases).Finally, 15 O yield was best predicted by 10/INCL, 10.1/INCL, 10.6/INCL, 10.7/INCL, 11/INCL and 11.1/INCL (9/12 cases) -again, in this case, even the worst-performing version/model combinations exceeded the threshold in 7/12 cases.In the build-up and Bragg peak region and tail region, the results are more clearcut.Results for the NMSE metric conclusively show that version 10.2/BIC is clearly the best choice for total positron yield, 11 C production and 15 O yield, while 10.5-11.1/INCLperformed the best for 10 C. Pearson CCC results are more mixed, but again, 10.2/BIC gives the best results for total positron annihilation yield, with most versions of Geant4 with BIC performing well.10.6/QMD performed the best for 11 C, 10/INCL and 10.1/INCL performed the best for 10 C, and there was no clear winner for 15 O.
The Pearson CCC results in the tail region were all very similar across Geant4 versions, with only a few wins separating the best and worst-performing version/model combinations in most instances.All version/models exceeded the threshold of 0.95 for a clear majority of cases for total positron yield as well as11C and15O production.
[32,33,6]6.6±0.8 mm for16O[32,33,6].Results were generally better for the the distal depth at 50% of peak metric. In, 10.2/QMD, 11.1/BIC and 11/INCL all achieved a mean of zero, with 10.2/QMD also having the equal-lowest maximum value of 1.5 mm (a depth difference of one voxel).QMD's maximal values were slightly smaller overall compared to BIC, and INCL's were the largest at ±3 mm for all versions.INCL tended to consistently overestimate the depth of this point, with both mean and maximum differences being negative in most cases.BIC and QMD both tended towards underestimating the 50%-of-peak depth, with the exception of version 11.1 (negative maxima for both, and means of 0 and -0.1 mm, respectively).Standard deviations were quite small for all versions and models (with the maximum standard deviation being 1.68 mm, for 11.1/INCL).Finally, in the tail region, Geant4 10.2 with BIC and QMD again provided the best prediction of total positron and 11 C yield in terms of NMSE, while 10.6/INCL performed the best for 10 C and15O.All version/model combinations performed well for total positron annihilation, 11 C and15O yield according to the Pearson CCC metric, while no version/model performed especially well for 10 C. Across all regions, ion species, beam energies, and target materials evaluated, the combination of Geant4 version 10.2 and BIC is best able to reproduce experimental results as evaluated using the NMSE and Pearson CCC metrics -especially in the buildup and Bragg peak region and tail region.Since the build-up and Bragg peak region is the location where (1) the majority of the dose resulting from carbon or oxygen ion beam irradiation in heavy ion therapy is deposited, and (2) where the strongest positron annihilation signal is observed, the results in this region are the most relevant to PET image-based QA simulation work.Version 10.2 also provided the best estimate of the depth of the distal depth at which positron yield decreased to 50% of peak, although this was obtained with QMD rather than BIC; the most accurate estimate of the depth of the peak itself was also achieved with QMD, but with Geant4 versions 10.4 and 10.6.As QMD exhibited the best accuracy and precision across Geant4 versions, it is the recommended model if the depth of the yield peak is critical.The BIC model implemented in Geant4 version 10.5 suffered from a run-time stability error which resulted in it being unable to simulate all test scenarios; therefore, we recommend that this version/model combination should be avoided for future studies.In the evaluation of individual positron-emitting fragment yield profiles, predictions of 10 C distribution were generally the least accurate in terms of both the NMSE and Pearson CCC.Interestingly, the INCL model often performed the best for prediction of 10 C fragment yield, although it rarely performed the best for total positron annihilation and 11 C or15O.Therefore, INCL should be considered for studies focusing on 10 C fragmentation, with the caveat that range estimation will be less accurate with this model.