Prediction of positron emitter distributions for range monitoring in carbon ion therapy: an analytical approach

Range verification is one of the most relevant tasks in ion beam therapy. In the case of carbon ion therapy, positron emission tomography (PET) is the most widely used method for this purpose, which images the -activation following nuclear interactions of the ions with the tissue nuclei. Since the positron emitter activity profile is not directly proportional to the dose distribution, until today only its comparison to a prediction of the PET profile allows for treatment verification. Usually, this prediction is obtained from time-consuming Monte Carlo simulations of high computational effort, which impacts the clinical workflow. To solve this issue in proton therapy, a convolution approach was suggested to predict positron emitter activity profiles from depth dose distributions analytically. In this work, we introduce an approach to predict positron emitter distributions from depth dose profiles in carbon ion therapy. While the distal fall-off position of the positron emitter profile is predicted from a convolution approach similar to the one suggested for protons, additional analytical functions are introduced to describe the characteristics of the positron emitter distribution in tissue. The feasibility of this approach is demonstrated with monoenergetic depth dose profiles and spread out Bragg peaks in homogeneous and heterogeneous phantoms. In all cases, the positron emitter profile is predicted with high precision and the distal fall-off position is reproduced with millimeter accuracy.


Introduction
The physical advantage of carbon ion therapy lies in its high ballistic precision in longitudinal and transversal direction, but it also makes it more sensitive to deviations from the treatment plan. Since carbon ions undergo nuclear reactions in tissue and create positron emitters like 11 C (T 1/2 = 20.39 min) and 15 O (T 1/2 = 2.03 min) (Tomitami et al 1993, Parodi et al 2002, positron emission tomography (PET) has been introduced as an in vivo treatment verification tool (Enghardt et al 1999, Shakirin et al 2011, Bauer et al 2013a. Due to different underlying physical processes, dose deposition and positron emitter activation are not directly comparable to each other, which prevents direct treatment verification. However, with the support of simulation platforms, a prediction of the positron emission profile can be obtained and its comparison to the measurement allows for treatment verification, when physical decay and biological washout are taken into consideration (Bauer et al 2013b). If the two distributions agree within a small range of uncertainty, the treatment is assumed to have been delivered correctly. In contrast, differences between the two distributions indicate deviations between the delivered and planned treatment and the consequential need for a revision and refinement of the treatment plan or patient setup.
Since dose is primarily delivered in electromagnetic processes and positron emitters are created in nuclear reactions, a reliable prediction of the PET image currently requires complex and time-consuming Monte Carlo (MC) simulations that consider the initial ion beam properties and all electromagnetic and nuclear reactions in tissue. To overcome this drawback and make range verification with PET imaging more convenient for clinical practice, Parodi and Bortfeld showed for proton therapy that there must exist a function whose convolution with the depth dose distribution yields the positron emitter distribution (PED) (Parodi and Bortfeld 2006). With an appropriate filter kernel obtained from fitting procedures, the one-dimensional positron emitter activity distribution can be predicted via convolution from a depth dose distribution with millimetre accuracy in very short time. The feasibility of this procedure for clinical implementation has been shown successfully in three dimensions for offline and in-room PET imaging in proton therapy (Parodi et al 2007a, Attanasi et al 2011, Frey et al 2014. Another approach for analytical prediction of PET monitoring distributions was suggested by Priegnitz et al (2012). Therein, the authors acquired experimental positron emitter yields in different phantoms of simple composition (graphite, polyethylene (PE), water and polymethyl methacrylate (PMMA)) to then create an analytical model to predict PET distributions.
In contrast to protons, carbon ions undergo projectile fragmentation which leads to different properties of the positron emitter activity profile. In particular, while positron emitters in proton therapy are created solely from target nuclei, the projectiles in carbon ion therapy fragment into positron emitters themselves as well, giving rise to a characteristic positron emitter profile maximum close to the Bragg peak of the dose distribution (Enghardt et al 1992).
In this work, we extend and modify the filtering approach to make it suitable for carbon ion therapy and predict the positron emitter distribution from the depth dose distribution analytically. The method is applied to depth dose distributions and its results are compared to positron emitter distributions which are created by MC simulations described in section 2.1. The convolution approach in proton therapy and our analytical approach to predict positron emitter distributions from dose profiles in carbon ion therapy in one dimension are presented in section 2.2. In section 3, the results of the application of this method to data from monoenergetic carbon ion beams of therapeutic energies and spread-out Bragg peaks (SOBPs) impinging on homogeneous and heterogeneous phantoms are demonstrated.
Finally, a discussion and outlook of this work is given in section 4.

Monte Carlo simulations
The simulation of positron emitter production in carbon ion therapy requires accurate models of nuclear interactions and their cross sections, as well as an accurate description of secondary fragment fluences. Since empirical data about these values are limited, MC simulations struggle to precisely predict secondary particle production in arbitrary cases (Böhlen et al 2010). For the development of PET predictions via convolution, MC simulations from Geant4 version 10.02.p02 are used in this work (Agostinelli et al 2003), under the assumption that they describe measured PET distributions within acceptable ranges of uncertainty (Böhlen et al 2010). For the modelling of physical processes, the physics list QGSP_BIC_HP (quark-gluon string precompound, binary light ion cascade, high precision neutron package (NeutronHP)) is chosen.
In the simulation, a monoenergetic carbon ion pencil beam impinges on a rectangular phantom of solid material, and positron emissions and deposited energy are scored in the phantom in slices in depth (binning = 0.5 mm). Since 11 C and 15 O are the most relevant positron emitters in carbon ion therapy (Enghardt et al 1999, Parodi et al 2002, only these two positron emitter isotopes are considered and tracked. Other positron emitters such as 10 C also play a role in the PET signal, but they were not considered for this study due to their lower production yields and shorter half-lives. Nevertheless, the methods presented herein can be easily extended to those cases. PMMA, water, PE, bone and a reference material serve as phantom materials. The so-called reference material is created for the sole purpose of this work and includes all human tissue elements leading to significant production of positron emitters. It is used as reference space in which filters are created and convolution is performed (section 2.2.2), where other materials could lack some of the considered elements. Properties of the phantom materials are displayed in table 1. Biological washout, detector responses and the time course of a potential measurement are not considered in this work. Parodi and Bortfeld (2006) suggested a framework in proton therapy to predict the positron emitter activity P(z) along depth z in tissue via a convolution of the depth dose distribution D(z) with a filter kernel f (z),

The convolution framework in proton therapy
Herein, the filter kernel f (z) is obtained analytically via deconvolution, using so-called Q ν functions that allow for a closed-form solution of the deconvolution problem in (1). Therefore, the dose distribution and positron emitter activity distribution along depth are fitted with a linear combination (one, and one to three terms, respectively) of these Q ν functions, which are built from a convolution of a Gaussian and a powerlaw function (Bortfeld 1997).
The filter parameters are refined via minimization of differences between the convolution of the filter kernel with dose distributions at different energies and MC simulated PED.
For each interaction channel of the primary ion with a target nucleus resulting in a positron emitter, one filter is created in the reference space. This framework is generalized to inhomogeneous media with a path length operator F that converts-analogously to the water equivalent path length formalism used for treatment planning-the depth z in a specific medium into depth z * in a reference space via where ρ el (z )/ρ el,ref is the electron density at depth z in the medium of interest relative to the one in the reference material, and the integration is performed over z along the depth in real space from the phantom entrance at depth z = 0 mm to the depth z of interest (Parodi and Bortfeld 2006). After transformation of the depth dose distribution in real space via the path length operator F into the reference space of 'reference' material (properties and composition see table 1), the convolution of the dose distribution with filter kernels is applied in the homogeneous reference space. Resulting distributions are transformed back to real space with inhomogeneous target materials using the inverse operator F −1 . Since the nuclear density of the target nuclei and the tissue characteristics influence the strength of the PET signal, the convolved positron emitter distributions are multiplied with weighting factors according to the tissue composition at each depth. This local factor g i (z) is defined via where w i (z)/w i,ref is the fraction by weight of the target nucleus involved in the reaction channel at depth z of the target material relative to the one in the reference material, and ρ(z)/ρ ref and ρ el (z)/ρ el,ref are the mass and electron density of the target medium at depth z relative to the one of the reference space, respectively. Finally, the total positron emitter distribution is calculated via where i sums over the contribution of N different positron emitters.

The application of the convolution framework to carbon ion therapy
In carbon ion therapy, positron emitter activity distributions show significant differences to the ones from proton therapy, since positron emitters do not only originate in the target nuclei from the tissue but also in the projectile particles themselves. 11 C for example can be produced as a fragment of the projectile 12 C as well as a fragment of a target nucleus, e.g. 12 C, 14 N, or 16 O. Consequently, the following main differences in the procedure of convolution in carbon ion therapy can be observed when comparing to the one of protons.
(i) Four Q ν functions are required to create a filter able to describe the positron emitter activity distribution reasonably well, which leads to a more challenging fitting procedure due to the larger parameter space.
The calculated filter is strongly energy-dependent and does not reproduce the positron emitter activity distribution (especially the peak height and width) when it is applied to other energies. This is illustrated in figure 1 where the filter developed for carbon ions impinging on the reference material with an energy of 200 MeV n −1 is applied to the distributions of the same carbon ion beam at different energies. Additionally to the peak, the build-up and tail are not reproduced well either, in contrast to the distal-fall off position which is predicted with deviations below 1 mm.
Due to the strong influences of the target material and ion energy on the PED, the successful application of the procedure with one filter for each positron emitter is not possible. In this work, we have developed a procedure that accounts for the different properties of positron emitter production in carbon ion therapy, i.e. the projectile and target fragmentation, the energy dependence, as well as the material dependence. Separation of the projectile contribution (positron emitters 11 C as fragments from the projectile 12 C) and the target contribution (positron emitters 11 C and 15 O as fragments from nuclei in the target medium) allows the prediction of positron emitter activity distributions in homogeneous phantoms. With various filters and the convolution method, the distal fall-off position of the positron emitter activity distribution is predicted, while the rest of the positron emitter activity distribution is fully characterized by the properties of the target medium, the energy of the incoming carbon ion, and empirically obtained analytical functions.

Separation of the target and projectile contributions
After a nuclear reaction resulting in a positron emitter, it is not possible to extract information in Geant4 whether the created positron emitter originated in the projectile or target nucleus. However, energy and momentum conservation considerations allow for a reasonable prediction. In particular, carbon ions of therapeutical energies only transfer a small amount of kinetic energy to the target nucleus at rest in the tissue in peripheral nucleus-nucleus interactions (Pshenichnov et al 2006). This process leads to the typical positron emitter activity distribution in carbon ion therapy: a broad build-up and tail from positron emitters created mainly from target nuclei that stay at their site of production, and an activity peak close to the Bragg peak position where the positron emitters formed from the projectile particle come to rest.
A plot of the kinetic energies of fragment particles 11 C after nuclear interactions of projectile particles 12 C with target nuclei in MC simulations are plotted against the kinetic energy of the projectile at the time of interaction in figure 2 (red dots) and reveals two groups of 11 C fragments. One group has a low kinetic energy (below 0.2 GeV) which is independent of the kinetic energy of the projectile in the interaction. The kinetic energy of the other group is higher and rises with the projectile energy. Based on the aforementioned considerations, it is reasonable to assume that the first group of low kinetic energy originates in target atoms at rest when they are hit by the projectile particle, while the high energy group corresponds to fragments of the projectile itself.
To separate the two groups, a 70% threshold of the incident kinetic energy of the particle is selected (shown in the figure in blue). Positron emitters with an energy above this threshold are considered to originate in the projectile particle ( 11 C proj ), while the ones with an energy below the threshold are considered to originate in target atoms ( 11 C target ). Since the two groups are well separated, any other energy threshold value between 50% and 90% leads to a very similar result. Small deviations occur only in the low energy region but do not have a significant impact on the final result, since only relatively few particles occur in this region.

The projectile contribution
The projectile contribution 11 C proj to the positron emitter profile is an order of magnitude larger than the one in the build-up and tail region and allows to determine the range of the positron emitters formed from the projectile particles. Therefore, it provides information correlated with the range of the incident ions in the medium, the most important quantity of interest in treatment verification tasks (Maccabee et al 1969, Kraan 2015. To predict the distal fall-off position of the positron emitter distribution, the convolution procedure is used. Therefore, positron emitter distributions and depth dose distributions from MC simulations are fitted with Q ν (z) functions as it has been described in the previous section, and filters are created analytically for each positron emitter and each target nucleus. The filtering approach is solely used to predict the distal fall-off position and complemented by other methods used to estimate the shape of the PED distribution.
With the distal fall-off position fixed, the next step is to estimate the PED peak width. This is accomplished calculating the range R f of 11 C isotopes in the target material, fixing the beginning of the positron emitter peak. To calculate the range R f , the concept of residual range is used. The residual range r f (E z * ) of the 11 C isotope depends on the residual energy E z * at the depth of interaction z * and can be estimated from the residual range r p (E z * ) of the projectile particle at this depth and the mass and charge numbers of projectile and fragment A p,f and Z p,f , respectively, via (Schardt et al 2010): For 11 C isotopes in 12 C therapy, this formula reduces to r11 C (E z * ) = 11 12 r12 C (E z * ). The depth R f at which a positron emitting fragment comes to rest can be estimated via a sum of the depth of interaction and the remaining range of the fragment particle, i.e.
where r p (E 0 ) denotes the residual range of the projectile at the beginning of its path with initial energy E 0 . To verify this linear relationship, the MC simulated number of 11 C particles created at different depths is plotted against their depth of positron emission in figure 3 and displays a good match between the calculated range and MC simulated ones which are smeared by a Landau distribution corresponding to the Landau distribution of the energy loss of the particles (Wilkinson 1996). As a consequence, it is possible to determine the number of positron emissions at each depth in the peak region of the 11 C proj profile based on a linear mapping of the material properties of their production sites along the path of the beam. The proximal margin of the positron emitter peak is determined by the 12 = 92% of the incident particle range, which in turn is determined by the 80% distal fall-off position of the MC simulated depth dose distribution. The proximal and distal limits of the 11 C proj peak define a mapping window in which the Kinetic energy E kin,frag of 11 C fragments after the interaction of the projectile 12 C of energy E kin,proj with a nucleus in the target material. Two different groups can be observed, separated by a 70% threshold line (blue). All fragments situated in the plot below the threshold line are considered to be target fragments. In contrast, all fragments above the line are considered to originate in the projectile particle. Data are obtained in Geant4. materials crossed during the ion path will be mapped into. For example, the distal limit of the mapping window for 300 MeV n −1 12 C ions in the reference material with a range of 12.3 cm is determined via the convolution method to be 12.3 cm. The shortest range of 11 C proj , and therefore the proximal limit of the mapping window, is for 11 C proj created when the nuclear interaction leading to their production occurs at the very beginning of the phantom (i.e. z * = 0 in (6)) and is therefore at a depth of 11 12 · 12.3 cm = 11.3 cm. The linear mapping of the total range to this window decreases the binning to a factor of 12.3 cm−11.3 cm 12.3 cm = 0.08. Since the binning in depth for the phantoms used herein was 0.5 mm, the contribution of each one of those bins to the 11 C proj peak will be squeezed into 0.04 mm.
With the limits of the mapping window determined, the aforementioned contribution must be calculated. This takes the form of three separate factors based on the materials along the path of the ions and their energy. The first factor is the g-factor, g i (z) (see (3)), the second is the fluence loss of 12 C, and the third is an energydependent scaling factor to account for the height of the positron emitter peak at different depths based on the fluence loss of the 11 C proj .
Regarding the fluence loss correction for 12 C, it shows an exponentially decreasing behaviour in depth: where Φ 0 is the fluence at the entrance of the phantom, N is the density of atomic nuclei in the tissue and σ is the reaction cross section. This behaviour can be approximated by a linear function of negative slope along the depth in tissue, where R is the particle range and m is a material dependent factor (Lee et al 1993).
On the other hand, the fluence loss of the 11 C along the depth in tissue was assessed with MC simulations and it was found to have a cubic dependency on the energy and is characterised by four coefficients h i,1−4 for each target nucleus i. This correction makes the positron emitter activity peak decrease in depth. Finally, to account for energy loss fluctuations, a Gaussian G(z) of empirically determined width σ G = 1.2 mm in the reference material is used.
Combining all corrections, the distribution of the predicted positron emitter profile of the projectile part 11 C proj that is to be mapped into the peak region defined above is given by An example of the projectile distributions predicted with this approach in comparison to the ones from MC simulations is shown in figure 4 for the interaction of the projectiles with carbon (figure 4(a)) and oxygen (figure 4(b)) targets.  (6) is indicated with a red dashed line, and the MC simulated 11 C projectile distribution is shown in solid green.

The target contribution
The contribution 11 C target to the positron emitter distribution originating in the target nuclei has a characteristic behaviour, as it is shown in figures 4(c)-(e), and is described best by analytical functions. In the part of the positron emitter distribution before the Bragg peak of the depth dose distribution, positron emitters are created in interactions from the projectile with the target nuclei. With increasing depth interactions of secondary particles with the target nuclei in the tissue also contribute to the number of positron emitters, leading to the build-up-like behaviour which is independent of the energy of the projectile. Therefore, this build-up part of the distribution is described for each positron emitter species i separately by the analytical function (11) that approaches the theoretical maximum number P i,b,max of produced positron emitters in tissue at large depths, where α i is determined by the number of positron emitters at the entrance of the phantom, and c i is a factor considering the build-up rate. For each target nucleus i and positron emitter species, this function is fitted to simulated data and characterized by the three different coefficients. It is evaluated from the beginning of the phantom until the maximum position of the convolved positron emitter distribution from the previous section, which is not equal to the maximum of the predicted projectile PED (compare figures 1 and 4(a) and (b)).
In the distal part of the distribution beyond the peak, the model does not describe the simulation results anymore, since the production of secondary particles declines. In this case, a hyperbola function was found to yield a good agreement to simulated results and is correlated with the energy E via where β i is a scaling parameter. The gap between the end of the build-up function (maximum position of convolved PED) and beginning of tail (10% distal fall-off position of the convolved PED) is interpolated linearly. An example of these target positron emitter profiles for projectile energies between 100 MeV n −1 and 400 MeV n −1 is shown in figures 4(c)-(e).

Combination of projectile and target contributions
To obtain the prediction of the total positron emitter distribution, the previously presented approaches are combined. For each element present in the target material that leads to significant production of positron emitters (here: hydrogen, carbon, nitrogen, oxygen, phosphorus, and calcium), a filter is created to determine the distal fall-off position of the predicted positron emitter distribution, the local factors g i (z) of the material are calculated and mapped to the depth of possible ranges of the 11 C and the aforementioned corrections are applied to find the projectile contribution of 11 C of each element. Subsequently, the target contribution is calculated for each element and positron emitter (here: 11 C and 15 O) using the analytical functions introduced in the previous paragraphs.

Validation
To show the feasibility of this approach, it is first applied to distributions in the case of monoenergetic irradiation of the homogeneous reference material (introduced in section 2.1 and table 1). Since SOBPs are clinically more relevant than monoenergetic irradiation, a SOBP is imitated by overlaying 25 monoenergetic carbon ion beams of small energetic difference (∆E = 2 MeV n −1 ) in the homogeneous reference material and in a slab phantom consisting of four layers of different materials: water (10 cm), bone (2 cm), polymethyl metacrylate (PMMA, 10 cm), and reference material (10 cm). Slab phantoms give the possibility to directly obtain laterally integrated distributions and the plots depicted in the figures shown are normalized to the volume of each slab. All analytically predicted positron emitter distributions PED conv (z) are compared to MC simulated ones PED MC (z) and evaluated via the mean relative error (MRE) to the maximum of the MC simulated positron emitter distribution for depths between z = 0 mm and z = 300 mm.

Monoenergetic irradiation of homogeneous phantoms
The 11 C (projectile and target contribution) and 15 O (target contribution) distributions of monoenergetic carbon ion beam irradiation at 100 MeV n −1 , 200 MeV n −1 , 300 MeV n −1 and 400 MeV n −1 of a homogeneous phantom are predicted analytically within less than a second and shown in figure 5 (red solid). For comparison, MC simulated, laterally integrated distributions of dose and positron emitters are shown in dashed green and solid blue, respectively. The predicted distributions match the simulation results for all energies. The MRE from the analytical description and MC simulation to the maximum of each MC simulated distribution is below 1% for the 11 C projectile contribution and below 3% for the 11 C and 15 O target contributions for all energies. The distal fall-off positions are predicted with a maximum deviation of 0.6 mm. A quantitative evaluation for each energy is shown in table 2.
Discrepancies for the target component are observed, especially for lower energies. This was expected since the parameterization of the build-up part was based on the parameter P i,b,max , which is the theoretical maximum number of positron emitters at large penetration depths. This approach ensures that the parameterization is applicable to any energy, with the consequence of an impact on accuracy for lower energies.

Spread out Bragg peak in homogeneous and heterogeneous phantoms
The predicted total 11 C distributions from the analytic approach in comparison to the ones from MC simulations are displayed in figure 6 for a spread out Bragg peak (242 MeV n −1 -290 MeV n −1 and 230 MeV n −1 -280 MeV n −1 ) in the reference material and a slab phantom, respectively. In the homogeneous reference material (plot (a)), the analytically predicted distribution match the MC simulated results well. Small deviations occur in the height of the positron emitter peak, giving rise to a MRE of 1.2% for the total 11 C distribution. When applying this procedure to the SOBP in the inhomogeneous slab phantom (plot (b)), the build-up and tail of analytically predicted and MC simulated distributions match again very well, while the height of the peak is underestimated with the analytical approach. The MRE amounts to 1.8% in this case. Such an underestimation is attributed to the small deviations that exist for each single energy distribution that pile-up when combined to form a SOBP. Figure 6 depicts also the 15 O distributions obtained from the target oxygen, which is the origin of most of the 15 O contribution. A systematic good agreement is obtained for the homogeneous reference material phantom (plot (c)) and slab phantom (plot (d)). Larger deviations are observed in plot (d) around the transition from water (medium I) to bone (medium II), which can be explained by the end of the range of some of the beams at or shortly after the transition. The prediction in water (medium I) is consistently lower than the MC simulations in the case of slab phantom (plot (d)), while such discrepancy is not observed in the case of the reference homogeneous phantom (plot (c)). This is attributed to the lower percentage of oxygen in the reference material with respect to water (38% versus 89%, see table 1), which makes any inaccuracy in the parameterization of the target contribution more noticeable.
The last row of figure 6 shows the sum of the first two rows, with similar trends to the 11 C distribution of the first row due to its major contribution to the PED.
It should be noted that the shape of the SOBP is not flat to account for realistic slope when delivering a biological flat SOBP.

Discussion and conclusion
In this work, the feasibility of an analytical approach to predict positron emitter distributions from depth dose distributions in carbon ion therapy was shown. Due to the strong dependence of the positron emitter distribution on the target properties and the incident carbon ion beam energy, it was not feasible to apply the filtering approach in the same way as in the case of proton therapy. Instead, the convolution approach was used to predict only the position of the positron emitter peak and its distal fall-off, while the positron emitter distribution itself is described with analytical functions. Therefore, a distinction is made between the contributions coming from different target nuclei involved in the nuclear reaction resulting in the positron emitters, the different positron emitters themselves (here 11 C and 15 O), as well as their origin, i.e. if the positron emitter is a projectile fragment or a target one. The positron emitter peak is generated by the projectile contribution of 11 C created in nuclear reactions with different target nuclei and is limited in depth. While its beginning is fixed to 11 12 of the projectile 12 C range, its distal fall-off position is determined via the 10% distal fall-off position of the convolved positron emitter distribution from the filtering approach. Local g-factors quantify the positron emitter distribution within the peak, considering the electron and mass density of the so far traversed medium, the fluence loss of the carbon ions and isotopes in the medium, and a Landau distribution blurring the edges. In contrast, the contributions from positron emitters originating in target nuclei shape mainly the build-up and tail region of the positron emitter distribution. While the build-up from the beginning of the phantom to the positron emitter peak is characterised by a sum of asymptotically growing exponential functions for each positron emitter and each target nucleus, the tail beyond the peak is estimated-again for each positron emitter and each target nucleus-with a hyperbolic function that includes a cubic dependency on the energy of the primary particles.
With this approach, positron emitter distributions in homogeneous and heterogeneous phantoms have been predicted with satisfactory agreement to MC simulated results. Mean relative errors to the maximum of the distributions were below 4% in all tested cases, and distal fall-off positions matched well with maximum deviations of 0.6 mm for homogeneous and 1.3 mm for inhomogeneous phantoms. While total distributions of monoenergetic irradiation of homogeneous phantoms are predicted with very high precision and only small deviations in the build-up region, also SOBPs are predicted with good agreement to MC simulated results in homogeneous as well as inhomogeneous phantoms. Acceptable deviations occur only in the absolute height of the positron emitter peak, while the shape is described well, also in proximity of slab interfaces.
Although the results shown in this work are satisfying, they still need to be taken under careful consideration of uncertainties in the systematic principle. Since all analytically predicted positron emitter distributions have been based on MC simulations, their benefit comes only into effect, if MC simulated distributions accurately and reliably describe measured distributions. However, different physics models in Geant4 already predict deviations in the positron emitter production and can deviate from measurements. Therefore, the analytical method is limited by the simulation accuracy of the MC simulation platform. Other groups work on the improvement of the reliability of MC predicted PET distributions which will make this method more feasible to real applications. Additionally, a method similar to the one from Bauer et al (2013b) is conceivable to refine the PET predictions to clinic-specific conditions.
The work of Priegnitz et al (2012) proposed an interesting approach. Although a thorough comparison between the aforementioned study and the one herein is very difficult due to the different goals and metrics involved, some discussion can be done. Herein, mean relative errors were below 4%, while in Priegnitz et al the normalized root mean square error (NRMSE) showed errors below 20%, when comparing predicted and measured distributions. The two metrics cannot be directly compared since the squaring of the residuals in the case of NRMSE will give more weight to larger deviations. Nevertheless, comparing the different plots herein and in Priegnitz et al (2012), the level of agreement in signal amplitude is similar. However, the striking differences between the two approaches are the simplicity and the work needed for its development. Priegnitz et al explicitly state that they want to avoid the filtering approach of Parodi and Bortfeld (2006) since it relies on cross section data with the use of Monte Carlo tools. However, their method relies on experimental yield data that, although easier to obtain than experimental cross section data, requires an extensive experimental campaign, and is impacted by the detector response and unavoidable uncertainties in separating the different isotope contributions with decay fitting of measured activity. If the yield data are obtained with Monte Carlo tools, inaccurate cross sections may have to be used, which (Priegnitz et al 2012) meant to avoid. However, their approach needs a considerable amount of tabulated data, which also resembles the work of Sterpin et al for prompt-gamma prediction in proton therapy (Sterpin et al 2015). In the present study, some parameterizations are done and several approximations are used, namely for the estimation of the signal amplitude of the distributions (build-up, peak height and width, and tail), hence introducing some complexity. Despite this disadvantage, the present work shows maximum range deviations up to 1.3 mm and, for this estimation alone, only six filter functions (i.e. to predict the 11 C distributions for the six target nuclei considered) and the dose distribution are required to predict the distal fall-off position. Finally, the method presented herein uses MC-based data but, if better cross sections, or other type of experimental or MC-based data are available, it can be further improved using such new inputs without much effort or time. Once the baseline filters and parameterizations are obtained, further improvements can be done with fine tuning of the filter functions parameters or the parameters of the parameterizations, thus making future developments considerably more practical than the approach of Priegnitz et al (2012), which is impacted by the amount of experimental data required and granularity of the methods employed.
More recently, the work of Pennazio et al (2018) shows the use of an innovative Monte Carlo tool to predict PET images in carbon ion therapy monitoring, and the comparison of its outcomes with measured PET data. A direct comparison with the results herein is difficult since different metrics were used and the study in Pennazio et al is with patients, which is more challenging. Nevertheless, based on the results shown, it is safe to consider that the tool proposed by Pennazio et al is more accurate in terms of both yields and prediction of distal fall-off positions than the one we propose herein. However, our solution offers some advantages. First, the computing time. In Pennazio et al (2018) the computing time required is not mentioned and, even if it performs faster than fullblown Monte Carlo packages, one can assume that the analytical method proposed here is considerably faster, providing results within a few seconds. Additionally, the solution herein has the potential to be further improved in terms of accuracy and extended, for example, for inclusion in a treatment planning system (TPS), as it was already demonstrated for the proton case using the filtering approach (e.g. Frey et al (2014)). Such an implementation in a TPS also allows for a relatively straightforward translation of the prediction tool to any therapy centre since it makes it depending solely on the built-in beam model. In the case of Pennazio et al, the approach requires an accurate description of the beamline in the simulation setup, as well as several calibrations and validations (Pennazio et al 2018).
The demonstrated results show the possibility of predicting PET profiles within a few seconds in different materials, which arises promising prospects for potential future applications in more complex cases. In the next steps, this method will need to be fully extended to three dimensions to be able to consider also the spatial distribution of positron emitters in three dimensions of an incoming ion beam of finite size, as well as extensive tests with patient CT data. Additionally, the results of the current filtering method and its parameters will need to be evaluated and potentially refined with measured data. In potential applications to living tissues, biological washout will blur the signal and will need to be taken into consideration (Toramatsu et al 2018), as well as the timecourse of irradiation and the decay of the positron emitters over time (Parodi et al 2007a). Finally, this method has the potential to replace computationally-intensive MC simulations to predict the PET images in a PET monitoring workflow. This can be achieved by applying the irradiation delivery timing parameters to obtain activity, and then couple the activity maps with either analytical methods to model the PET scanner response or simulate the transport of the 511 keV photons to the scanner. Regarding the aforementioned scanner response models, a possibility is to approximate the scanner response by convolving the production data with a three-dimensional Gaussian kernel, which is especially suitable to full-ring PET scanners (e.g. Parodi et al (2007b), Attanasi et al (2011) and Bauer et al (2013a)). In the case of explicit simulations of the photons, the time-consuming modelling of ion transport is circumvented, while the photon transport can be done very quickly. A realistic number of photons to be transported in the context of this work require seconds or a few minutes, while ions require hours or even days, depending on the uncertainty level. The method presented herein can also contribute to clinical decisions regarding treatment plan parameters that may allow for better PET monitoring, or be used as input to other methods, such as reconstructing dose distributions based on PET monitoring (e.g. Hofmann et al (2019)).