If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
Dipartimento di Elettronica, Informazione e Bioingegneria - Politecnico di Milano, Via Ponzio, 34/5, 20133, Milan, ItalyINFN - sezione di Milano, Via Giovanni Celoria, 16, 20133, Milan, Italy
INFN - sezione di Milano, Via Giovanni Celoria, 16, 20133, Milan, ItalyDipartimento di Energia - Politecnico di Milano, Via Lambruschini, 4, 20156, Milan, Italy
Dipartimento di Elettronica, Informazione e Bioingegneria - Politecnico di Milano, Via Ponzio, 34/5, 20133, Milan, ItalyINFN - sezione di Milano, Via Giovanni Celoria, 16, 20133, Milan, Italy
INFN - sezione di Milano, Via Giovanni Celoria, 16, 20133, Milan, ItalyDipartimento di Energia - Politecnico di Milano, Via Lambruschini, 4, 20156, Milan, Italy
Dipartimento di Elettronica, Informazione e Bioingegneria - Politecnico di Milano, Via Ponzio, 34/5, 20133, Milan, ItalyINFN - sezione di Milano, Via Giovanni Celoria, 16, 20133, Milan, Italy
Monte Carlo feasibility study on Prompt Gamma Imaging in carbon ion therapy.
•
The prompt-gamma signal is correlated with the depth-dose curve produced by C-ions.
•
The technique is unfeasible for one irradiation spot due to low counting statistics.
•
A range sensitivity of 4 mm is achievable with our set-up for an irradiation spill.
Abstract
Purpose: In-vivo range verification has been a hot topic in particle therapy since two decades. Many efforts have been done for proton therapy, while fewer studies were conducted considering a beam of carbon ions. In the present work, a simulation study was performed to show whether it is possible to measure the prompt-gamma fall-off inside the high neutron background typical of carbon-ion irradiation, using a knife-edge slit camera. In addition to this, we wanted to estimate the uncertainty in retrieving the particle range in the case of a pencil beam of C-ions at clinically relevant energy of 150 .
Methods: For these purposes, the Monte Carlo code FLUKA was adopted for simulations and three different analytical methods were implemented to get the accuracy in the range retrieval of the simulated set-up.
Results: The analysis of simulation data has brought to the promising and desired precision of about 4 mm in the determination of the dose profile fall-off in case of a spill irradiation, for which all the three cited methods were coherent in their predictions.
Conclusions: The Prompt Gamma Imaging technique should be further studied as a tool to reduce range uncertainties affecting carbon ion radiation therapy.
]. Today, particle therapy (PT) has become a well established option for tumour treatments, next to surgery, chemotherapy, conventional radiotherapy and emerging immunotherapy. More than one hundred facilities all over the world can provide beams of protons, carbon, or other ions for clinical treatments [
It is known that ions produce a characteristic dose distribution with a maximum, called the Bragg peak (BP), followed by a steep dose fall-off, when they interact with matter. The opportunity to have such a step in the dose distribution is the main advantage of charged particles over photons, so that it allows to precisely hit the tumour site, sparing deeper normal tissues. Moreover, tissues located upstream of the BP also get smaller dose than in conventional therapy. However, the intrinsic precision of charged hadrons can be a double-edged sword, as particle range is not exempted from uncertainties, which are closely related to the imperfect knowledge on the patient morphology (see [
] for an extensive discussion on range uncertainties). For this reason, if the BP position is mislocated with respect to the treatment plan, it can result in a non-effective treatment of the pathology, putting in danger nearby normal tissues, therefore making the advantageous feature of hadrontherapy not fully exploitable.
Due to range uncertainties, safety margins are applied in the treatment planning and this limits the full benefits of PT. To improve precision and restrain normal tissue irradiation, range uncertainties must be reduced. Thus, efforts have been made to develop instruments aiming at the precision of a few mm or better in verifying the particle range in-vivo [
Different approaches have been investigated to reduce range uncertainty through in-vivo range verification. They can be classified in direct methods, based on direct measurement of the depth-dose distribution, and indirect methods, which take advantage of the secondary emission from the patient during treatment. A general comparison and overview of the techniques used for in-vivo proton range verification are summarized in various review articles [
Online control of the particle range can be performed by means of secondary radiations emitted from nuclear reactions, following the interaction of the primary ion beam with human tissues [
]. Indeed, for some types of secondary radiation, there exists a correlation between the secondary emission and the depth-dose distribution, so that different techniques have been developed in order to monitor the particle range. In-beam Positron Emission Tomography (PET) has been already applied clinically at GSI [
]. The detection of secondary charged particles, as protons in carbon beam irradiation, has been also adopted in inter–fractional monitoring in clinical trials [
Stichelbaut F, Jongen Y. Verification of the proton beam position in the patient by the detection of prompt gamma-rays emission. In: 39th meeting of the particle therapy co-operative group, Vol. 16. 2003.
], with measured inter-fractional global range variations in the order of . A second prototype of prompt gamma camera utilizing a knife-edge collimator design was tested, as reported in [
In recent years, there has been a substantial progress in the application of PGI mostly for what concerns proton therapy. On the other hand, fewer studies have been made in the use of PGI verification within Carbon Ion Radiation Therapy (CIRT). Indeed, even if there are undoubted advantages of C-ion over proton and photon therapy [
], it is also true that the very same factors that make C-ions such interesting may hinder the application of PGI monitoring. For example, although C-ions have a higher prompt gamma yield with respect to protons [
], they also have a higher neutron yield. In addition to this, for ions heavier than protons, the number of incident projectiles needed to provide a given physical dose is littler than for protons, due to the dependence of energy loss (Z being the atomic number) and the smaller mean multiple scattering angle of heavier ions. Moreover, the Relative Biological Effectiveness (RBE) leads to a further reduction of the number of ions necessary to deliver the desired biological dose with higher-Z ions. Therefore, the total number of generated PGs for a given pencil beam spot is lower when considering heavier ions, strongly affecting their counting statistics [
]. The aforementioned factors make the optimization of detection efficiency and gamma-neutron discrimination even more critical for the application of PGI to CIRT. So far, Time-of-flight (TOF) discrimination technique has been successfully applied to the observation of PGs in carbon ion irradiation [
Monte-Carlo simulations play not only a crucial role in the design and optimization of detection systems, including those dedicated for PGI, but also in techniques for range shift detection. GEANT4 is the most used framework dedicated to this field, even if modifications to PG emission modelling in newer software versions have shown to increase the discrepancy between the simulation results and the experimental data [
]. On the other hand, the FLUKA simulation tool has been also used to reproduce the PG emission and good agreement was found both in the spectrum shape, and in the energy-integrated PG depth profile [
]. The results of the Monte Carlo (MC) simulations showed a weak correlation between PG and BP profiles. Anyway, the author himself questioned whether the MC code he used (MCNPX version 2.7.E) was mature enough for the simulation of prompt gamma emission during carbon ion irradiation at clinical energies. Hence, we decided to take the precious legacy of the work conducted by Smeets and, in the present work, we investigated whether with updated Physics models inside the general purpose FLUKA MC code results come out differently. We thus conducted a numerical study to explore the feasibility of using a knife-edge slit prompt-gamma camera with a beam of C-ions at clinical energies. Firstly, we simulated the response of the detector to the delivery of a single pencil beam spot to a uniform soft-tissue phantom. In the second part of the study, we examined the cumulative signal obtained from the delivery of multiple spots within the same transversal layer, seeking the chance to still reconstruct the position of the BP. The final goal of this work is to explore the PG fall-off retrieval capability with a well-established slit camera configuration, already tested clinically with proton irradiation [
]. Specifically, we scored the response of a pixelated knife-edge slit camera to the secondary particles emitted by a phantom irradiated with a C-ion beam.
Although the irradiation geometry was slightly varied for the last simulation set, the generic setup is showed in Fig. 1. A mono-energetic pencil beam of C-ions with a 150 energy was directed along the -axis of the frame of reference towards one of the faces of the phantom, a square cylinder of 30 cm diameter made of ICRP soft tissue. The energy of choice fits within the interval of therapeutic energies used in CIRT (at CNAO the beam energy range is 120–400 [
]. Thus, the irradiated face of the phantom was positioned at 5 cm on the -axis, so that the dose peak was almost centred on the origin of the reference frame (see Fig. 1(a)). We will label this as the “no-shift” position. The detection system of reference for the simulations was a gamma camera similar to the one already developed and characterized in previous studies for proton therapy [
]. Such camera consists of a knife-edge collimator and a pixelated detector. The former is made of a tungsten alloy, 4 cm thick on the radial direction and 30 cm wide along the beam axis. The knife-edge slit is 6 mm wide and has a 63°aperture angle. The detector is an array of 16 pixels (slabs) made of Lu1.8Y0.2SiO5 (LYSO) scintillator, distributed along the beam axis with a pixel pitch of 6 mm. Each pixel is a 10 cm wide slab with a 3 cm thickness. As shown in Fig. 1, the slit and the face of the detector were placed at a distance of 25 cm and 50 cm from the -axis, respectively, giving a field-of-view of 10 cm on the phantom. The detector and the slit were centred on the origin of the frame of reference to maximize the PG signal. The whole setup was defined inside a spherical room of 4 m radius filled with air. See Table 1 for details about the materials.
Fig. 1Schematic representation of the simulation geometry employed in this work: (a) 2D representation on the xz-plane; the depicted position of the phantom is labelled as the “no-shift” position; (b) 3D representation.
We used the FLUKA “DETECT” estimator to score the spectrum of deposited energy in the detector on an event-by-event basis, considering each pixel as an independent detector. Only the energy deposition events in the interval between 3 and 7 MeV were scored. The most probable PGs (4.4 MeV of 12C, 5.2 MeV of 15O and 6.1 MeV of 16O) are emitted within this energy window, as reported by Kelleter et al. [
]. The spectrum was scored in 1024 energy bins, having a width of about 4 keV. For each pixel, the total number of events (or counts) in the scored energy interval was considered as the detector signal. Since DETECT does not discriminate the type of particle interacting with the detector, and since no particle was discarded in the simulation, the signal was due to all secondary particles able to reach the detector, i.e., mainly photons and neutrons. Anyway, the FLUKA subroutine called “MGDRAW.f” allows to save in a phase space file all information about a particle (type, energy, position, direction) crossing a certain geometry region. For instance, recording all transport details about photons entering the detector region in a file and the same for neutrons in a second file, it is possible to separately transport them evaluating which are the two contributions to the correlation between dose distribution and secondary emission particles following carbon ion irradiation. Moreover, all FLUKA outputs were given per unit incident C-ion.
Specific transport and production thresholds were introduced to cut off useless particle histories without affecting the physics of the problem. In particular, transport thresholds inside air and collimator were set at 3 MeV for photons and electrons. Moreover, production thresholds of photons and electrons inside the target were set at 3 MeV. Furthermore, fluorescence was turned off in the whole simulation volume, except for the detector. In the latter, the electron and photon transport thresholds were set at 40 keV. This energy is related to the K-shell of the element with the highest Z inside the detector, which results in a mean free path of the particle (electron or photon) considerably smaller than the pixel width of 6 mm. Neutron generation by Fermi break-up of nuclei in the phantom was also introduced.
For each case of the simulation sets described in the following three paragraphs, primary particle histories were launched. Such number was chosen as it allows to have a relative standard deviation lower than 1% in the quantity of interest, i.e., the number of energy deposition events between 3 and 7 MeV in each pixel of the detector. We also performed one simulation with the target in the “no-shift” position, using the “MGDRAW.f” subroutine to separate the photon and neutron contribution. For this specific case, we used primary particle histories.
Table 1Composition and density of the materials assigned in FLUKA to the various regions of the simulation volume: the composition of ICRP soft-tissue is given in mass fraction, while the remaining elements are expressed in atomic percentage.
2.2 Particle range sensitivity determination by target longitudinal shifts
For this first set of simulations, a 2D representation of the set-up is reproduced in Fig. 1(b). The axis of the phantom and the beam lied along the -axis. The goal of the first set of simulations was twofold. The first one was to investigate whether the spatial correlation between the BP position and the PG profile along the beam axis can be identified even in the high neutron background of CIRT. The second one was to explore the theoretical limits in the determination of PG fall-off in C-ion irradiation with the delivery of a single irradiation spot to the phantom. To this aim, we performed several simulations using the same described configuration, but shifting the position of phantom along the -axis, effectively shifting the position of the BP. This is a typical procedure that is used in experiments to characterize the sensitivity of a given PGI system in determining the particle range (see [
]). From the “no-shift” position, the phantom was translated with a 2 mm step both in the positive and negative direction of the -axis, and the same simulation was repeated. Such a step was selected to investigate the feasibility of discriminating a 2 mm shift of the BP, the minimum shift we hoped to find. In addition to this, a step of 2 mm allowed to reduce the number of simulations in case of 1 mm-step and to get a denser study with respect to 3 mm-step.
For each shift of the phantom, we found that primary particle histories were enough to obtain a low Monte Carlo statistical error.
2.3 Delivery of multiple irradiation spots at constant energy
In the second set of simulations, we evaluated the response of the gamma camera to the delivery of multiple irradiation spots within a single layer of the phantom, shifting the beam axis along the xy-plane, but keeping its energy constant. This emulates what happens in pencil beam scanning (PBS) systems of synchrotron based treatment facilities, in which the spots of a given layer are delivered in a single “spill” of the accelerator (more information about spills in [
Peters A, Cee R, Haberer T, Hoffmann T, Reiter A, et al. Spill structure measurements at the Heidelberg ion therapy centre. In: European particle accelerator conference. 2008, p. 1824–6.
]). We will thus refer to this case as a “spill irradiation”. As shown in Fig. 2, the nine selected spots have the same z-coordinate but are distributed on the xy-plane with a pitch of 1 cm, covering a 2 × 2 cm2 area centred around the axis of the cylindrical phantom. Such an area was selected as it represents a realistic transversal section of a tumour volume (e.g., at CNAO the 2 × 2 cm2 is the minimum field size [
]). We chose to simulate only the central spot and those delivered to the edges of the tumour, as the latter are the most critical for range retrieval when the signal is integrated in the whole spill.
Fig. 2Illustration of the irradiation spots (red) selected for the second set of simulations. The spots have the same z-coordinate, and cover a 2 × 2 cm2 area on the xy-plane inside the soft-tissue phantom (yellow). The pitch between the spots is 1 cm. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
2.4 Tilted target simulations to account for non-planar shape of patient surface
For the last set, we still wanted to examine the detector response to a spill irradiation. However, in a real treatment, the covered spots do not have all necessarily the same distal coordinate. Indeed, even though the beam energy is constant, the non-planar shape of the patient surface and the inhomogeneities introduce a variation in the beam range from spot to spot. Instead of adopting an anthropomorphic phantom, the cylindrical phantom was inclined with respect to the -axis by a 22°angle (Fig. 3). In this way, the spots have a different z-coordinate, and the geometry of the setup can be easily replicated in an experimental session. Moreover, the inclination angle of the phantom was chosen to get, when moving the beam of 1 cm along the xy-plane, a particle range variation of 4 mm, which is the maximum shift we selected in the first set of simulations. Hence, good results for this kind of geometry would justify studies with a more realistic one. The detector geometry and the spot spacing in the x-y plane were still those described in the previous subsection.
Fig. 3Simulation geometry of the third set of simulations: the target was inclined by a 22°angle with respect to the -axis, so that the spots of Fig. 2 have a different z-coordinate. This takes into account the variations of the beam range in a patient from spot to spot, even with the beam at constant energy.
We performed an analysis to determine the accuracy in range retrieval of the PG camera, and its dependence on the number of C-ions delivered to the phantom. The algorithm was based on the approach used by Roellinghoff et al. [
Firstly, we obtained the signal curve from the FLUKA data. This curve represents the total number of counts per unit primary particle in each pixel as a function of the pixel position on the -axis. Then, we chose a given number of incident C-ions Ni and multiplied the counts in the signal curve by this number. The result of this operation was labelled as the “reference curve”. It is worth noting that the statistical uncertainty of the reference curve is only affected by the number of primary histories simulated in FLUKA () and not on the chosen number Ni. To account for the effect of a given number of C-ions on the counting statistics, we performed a MC propagation of the error on the reference curve. Namely, for each pixel we sampled a given number Nc of values from a Poisson distribution having a mean value equal to the number of counts of the reference curve. In this way, we obtained a given set of Nc noisy curves, each one estimating the number of counts and the statistical fluctuation of the signal corresponding to the selected Ni. Lastly, every noisy curve was interpolated with cubic smoothing spline method, to get a continuous curve (see Fig. 4 for some examples of the interpolated noisy curves).
The final step of data analysis was to determine the precision of our detection system in retrieving the position of the dose fall-off, or its shift from the nominal position. To do this, we compared the obtained noisy curves with the reference one by employing three different methods to quantify the closeness of the curves. An example of Method 1, which is similar to the one used by Roellinghoff et al. [
], is shown in Fig. 5. For each noisy curve, the associated reference curve was shifted from its original position with a step of 0.02 cm for a total path of 1 cm, both in the positive and negative direction of the -axis. To retrieve the shift for which the reference and the selected noisy curves were superposed the most, a root mean square error (RMSE) was defined as figure of merit to quantify the degree of superposition:
(1)
where represents the generic point of the reference curve with shift and is the point of the generic noisy curve . is the number of points which constitute the reference and each noisy curve.
Fig. 4Graphical representation of the generated noisy curves for different phantom shifts and numbers of C-ions: (A) , (B) , (C) , (D) .
We then found the optimal shift of the reference curve that resembled the most the noisy curve under analysis by minimizing the RMSE (Fig. 6). The explained process was repeated for all generated noisy curves, obtaining a histogram of the optimal shifts, which could be approximated to a normal distribution, according to the central limit theorem. We defined the precision in retrieving the BP fall-off as twice the standard deviation of the resulting normal frequency distribution.
Fig. 5Graphical explanation of Method 1: the reference curve (blue) obtained by MC simulations is shifted up to 1 cm, in 0.02 cm steps, on the forward (green) and backward (purple) direction of the -axis. The shifted curves are then compared with the investigated noisy curve, and the best match is found by a minimization of the RMSE curve (Fig. 6). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 6Typical RMSE as a function of shift for the four different statistics and the phantom shift of +2 mm. The first 500 curves out of are plotted for each C-ion number.
For Method 2, refer to Fig. 7. The range was associated with the inflection point of reference curve (the same holds for Method 3). Each noisy curve (yellow) underwent a linear interpolation (purple) in the fall-off region, where the code selected the two y-coordinates of the noisy curve cm far from the x-coordinate of the inflection point. The latter was found as the local minimum of the reference curve derivative, obtained via the MATLAB “gradient” function. The choice of the linear interpolation interval was the outcome of several attempts. Choosing a range between cm and cm the results were coherent, but the thinner the interval the higher the dependence on the local shape of the noisy curve. Instead, the selection of an interval wider than cm brought to a line unrepresentative of the fall-off region because the maximum of the curve ended in the interpolation region. Then, the algorithm was written to find the abscissa (green) of the intersection between the linear interpolation and the ordinate (blue) of the inflection point of the interpolated reference curve (red). The distance between the green and black line (x-coordinate of the inflection point) was chosen as the error in the similarity between the specific noisy curve and the reference one.
Repeating this procedure for each noisy curve, the code produced a frequency distribution of positions in which the interception has occurred. As before, twice the standard deviation of the resulting distribution gave the uncertainty on finding the particle range.
Fig. 7Graphical explanation of the Method 2: a generic noisy curve (yellow) is interpolated in the expected fall-off region (pink), which is then crossed by the y-coordinate of the reference curve inflection point (blue). The difference between the x-coordinates of the crossing point (green) and of the inflection point is chosen as the error in the fall-off estimation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Method 3 is explained graphically in Fig. 8. The idea was to find the intersection between the noisy curve (yellow) and the ordinate of the inflection point (blue) of the reference curve (red). The difference between the abscissas of the inflection point (black) and of the intersection point (pink) was chosen as the reconstruction error. Repeating this procedure for all the noisy curves, a normal distribution was again obtained.
We found that noisy curves were enough to enable the determination of the range retrieval precision with negligible statistical uncertainty for all methods.
Finally, regarding the number of incident C-ions, we performed the analysis of the first set results using Ni equal to and . As reference to typical values used in therapeutic irradiation (see [
]), the first Ni is close to the typical number of C-ions used in an irradiation spot and little lower than the minimum number of C-ions typically used in an irradiation spill (about ), while the second and the third are, respectively, in the middle and in the order of the maximum number of C-ions used in a spill. The last number is not found in a usual treatment, but it is considered as an upper bound inside our evaluation, which may be of interest for future developments in FLASH therapy [
]. For the second and third set of simulations (respectively, Sections 2.3, 2.4), we performed the analysis using the average number of C-ions delivered in a spill, i.e., .
Fig. 8Graphical explanation of Method 3: the algorithm finds the x-coordinate (green) corresponding to the interception between a generic noisy curve (yellow) and the y-coordinate (blue) of the inflection point of the reference curve (red). The error in the fall-off estimation was defined as the distance between the x-coordinates of the interception and of the reference curve inflection point. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
3.1 Particle range sensitivity determination by target longitudinal shifts
Fig. 9(a) shows the profile of the total counts per unit incident C-ion for each pixel, with the phantom in the “no-shift” position (see Section 2.1). As for all the reference curves presented in the Results section, the error bars refer to the MC statistical uncertainty given by FLUKA. The reference curve (blue) is also compared with the mirrored BP profile (orange) obtained with FLUKA for the same simulation geometry and beam energy. The signal has a global maximum followed by a fall-off, that is clearly correlated to the BP. It is worth noticing that no discrimination was applied on the detector, so the signal is mainly due to the superposition of the PGs and the neutron background. Indeed, Fig. 9(b) separately shows the photon and the neutron contribution to the total detection signal. It is clear that there exists a correlation between dose distribution and photon profile, while neutron energy deposition does not follow the step decrease of dose distribution. Without considering the possible losses due to the coupling between scintillator and readout in a real detector, the expected detection efficiency of the camera is of the order of .
The signal profiles for the “no-shift” position and the symmetric translations of the phantom by 2 mm and 4 mm along the -axis are displayed in Fig. 10. The curves are clearly distinguishable, and both their maximum and fall-off move together with the BP. Considering that these are ideal curves, as they do not include the error due to the counting statistics, a 2 mm shift is clearly distinguishable at least on a theoretical basis.
Fig. 9Results of the correlation between the PG detection profile and the Bragg peak profile (orange), obtained with FLUKA simulations for the phantom in the “no-shift” position: (a) Total signal detection profile (blue); (b) Photon contribution to the detection profile (blue) and neutron contribution to the detection profile (black). Notice that the photon and neutron contributions were simulated separately, sampling the target secondary emission with the FLUKA subroutine “MGDRAW.f”. In this specific simulation, primary histories were simulated. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Results of the data analysis algorithm applied to the reference curves are presented in the following figures and are summarized in Table 2. The application of Method 1 to these noisy curves produced the distributions illustrated in Fig. 11, while Fig. 12 refers to Method 2. We did not show the results of the Method 3, as they were qualitatively identical to those of Method 2 (see Section 2.5 for a description of the methods). An advantage of the latter two methods is the opportunity to see how the frequency distributions associated to the sets of noisy curves are located along the detector axis.
Fig. 10Range shift sensitivity: reference curves obtained by MC simulations for various shifts of the phantom along the -axis.
As expected, the precision in retrieving the fall-off of dose distribution is strictly connected to the number of incident C-ions. As reported in Table 2, the three methods gave similar results and Method 1 involved the lowest except for case (A). Fig. 13 shows the plot of the precision retrieved by Method 1 versus the number of C-ions on a double logarithmic scale, considering the averaged precision over the shifts. The trend-line of the double-logarithmic data (i.e., log(Ni) vs. log(2)) is a line with slope of about -, thus the accuracy is found to be inversely proportional to . Considering the average results given by the three methods, a number of ions of about would be necessary to achieve an accuracy of 2 mm in a single spot irradiation. Finally, all Gaussian curves are perfectly distinguishable at , showing that with a very high number of particles even a sub-millimetre precision is achievable.
Fig. 11Results of Method 1. The frequency distributions pertaining to the various shift of the phantom are presented for different numbers of C-ions: (A) ; (B) ; (C) ; (D) . Note the different scale in the -axis of the plots.
Fig. 12Results of Method 2. The frequency distributions pertaining to the various shift of the phantom are presented for different numbers of C-ions: (A) ; (B) ; (C) ; (D) . Unlike for Method 1, the centroids of the Gaussian curves shift together with the phantom position. Note the different scale in the -axis of the plots.
Table 22 with associated uncertainties of the Gaussian frequency distributions resulting from the data analysis of the first set of simulations. Results are reported in mm, for the various numbers of C-ions: (A) , (B) , (C) , (D) .
Table 22 with associated uncertainties of the Gaussian frequency distributions resulting from the data analysis of the first set of simulations. Results are reported in mm, for the various numbers of C-ions: (A) , (B) , (C) , (D) .
Fig. 13Range retrieval precision (blue) obtained by Method 1 plotted against the number of incident C-ions on a log–log scale. The precision is averaged over the shifts of the phantom. The linear interpolation of the double logarithmic data (i.e., log(Ni) vs. log(2)) has a slope of about . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
3.2 Delivery of multiple irradiation spots at constant energy
Fig. 14 illustrates the results of the second set of simulations. The profiles of the signal per unit incident C-ion are represented for the nine selected spots. The curves appear grouped in subsets of curves pertaining to a specific shift of the beam along the -axis. Indeed, since the detector was only pixelated on the -axis, the signal was almost independent on the shift of the beam along the -axis. In contrast, moving the beam along the -axis, the distance between the beam and detector was varied, thus changing the geometrical acceptance and signal attenuation (negligible for PGs with energy higher than 3 MeV in soft-tissue) without affecting too much the shape of the curve. This is further demonstrated in Fig. 15, in which the curves were compared after a peak normalization, showing that the normalized profiles are very similar for all the curves.
As stated in the Methods section, we were interested in the signal coming from the delivery of an accelerator spill, to study the feasibility of applying PGI on a layer-by-layer basis. Thus, we summed the curves corresponding to the nine spots, obtaining a reference curve for the whole spill. We then performed the analysis described in Section 2.5 to the reference curve, considering a total incident number of ions equal to . Consequently, we assumed that these ions were delivered homogeneously in the nine selected spots. The results of the analysis are shown in Fig. 16. The accuracy is given again by twice the standard deviation of the frequency distribution. In this case, we can state that the accuracy is referred to the capability of retrieving the range of a given energy layer covered by the beam, in the simplifying hypothesis of a perfectly planar transverse layer. The estimated accuracy given by the average of the three methods is .
Fig. 14Graphical results of simulations for the nine selected spots (see Fig. 2). The signal is varied only by a shift in the -direction.
Fig. 16Results from range retrieval precision analysis. The Gaussian distributions are found applying the aforementioned methods, whose precision is summarized in the legend.
3.3 Tilted target simulations to account for non-planar shape of patient surface
In this subsection, we present the results of the simulations of C-ion irradiation on the tilted phantom. The detected profiles are plotted in Fig. 17. As for the previous set, we calculated the integral signal of all the considered spots, for a total number of incident ions. We then performed the analysis, whose results are shown in Fig. 18. The frequency distributions of Method 2 and 3 are shifted with respect to the previous case, as a consequence of the target tilt. However, the accuracy was not worsened in the case of a tilted layer, as the average precision derived from the three methods is .
Fig. 17Reference curves obtained with FLUKA for the nine selected spots, in the case of the tilted target.
Fig. 18Results from range retrieval precision analysis in the case of the tilted phantom. The Gaussian distributions are found applying the aforementioned methods, whose precision is summarized in the legend.
Results about PGI fall-off determination for different numbers of C-ions incident on a soft-tissue phantom have been obtained using MC simulations with the FLUKA code, which is now considered a reliable tool for simulations in hadrontherapy, both for the prediction of dose delivery and PG emission [
]. The most recent FLUKA models for PG production showed a good reproduction of the relative shapes of experimental photon profiles as a function of depth as well as the absolute photon yield per primary C-ion, with an accuracy of about 15%–20%.
MC simulations showed that, for the beam energy and the detection system of interest, the signal profile due to the secondary particles emitted from the delivery of a single spot to the phantom was correlated with the BP. Moreover, we demonstrated that the shape of the profile was clearly determined by the prompt-gamma signal over an almost uniform neutron background. To probe the theoretical limitations of our gamma camera, we numerically reproduced an experiment for the measurement of the sensitivity in particle range retrieval. The reference curves for 2 mm shifts of the phantom were clearly distinguishable, opening the possibility to detect 2 mm shifts of the BP with a sufficiently high number of incident ions. Indeed, in accordance with Krimmer et al. [
], we obtained that the accuracy of the system in range retrieval scales with the number of incident ions Ni as 1/. However, the number of ions delivered in a single spot has proved to be too small for reaching that level of precision. From the trend-line of data fitting, we found that at least around particles were needed to go below a 2 mm accuracy with the present setup, without considering sources of experimental error or factors that may affect the detection efficiency in a real detection system. For instance, the main factor we expect is the dead time of the readout electronics and data acquisition system which, in case of high rates, may limit the counting efficiency and linearity of the detector response. Thus, we do not expect this technique to be effective on a spot-by-spot basis due to the low counting statistics of prompt-gammas, at least in the considered knife-edge configuration. However, different improvements could be introduced to increase the detection efficiency of the camera. For example, a first step could be to enhance the geometrical acceptance by using multiple cameras. Then, the intrinsic efficiency of the detector could be improved by removing the mechanical collimator, as under study in the development of Compton Cameras. Nonetheless, particles is actually the maximum amount used in an irradiation spill at CNAO (see [
]), delivered by their synchrotron to cover a given layer of the tumour volume. Thus, we proceeded to verify whether the distal position of the beam in an energy layer could be retrieved with a good accuracy. For this, we used ions, i.e., the average number of particles in a spill according to the interval found in [
]. Consequently, we combined the PGs signals obtained from 9 spots covering a 2 × 2 cm2 area, assuming that the ions were divided equally among the selected spots. We found that a 4 mm accuracy is achievable with the present camera. Even introducing a beam range variation within the same layer by tilting the phantom (see Section 2.4), the accuracy was not worsened.
Although our results suggested that range verification on a spot-by-spot basis seems not to be feasible, it may be worth to further investigate PGI technique for CIRT, at least on the single spill basis. Even though a layer-by-layer approach in range verification is not strictly a real-time technique, it could still be useful to reduce the uncertainties during the treatment and also for post-treatment verification. Differently, having in mind that our results were dominated by the counting statistics, for the number of C-ions used in a single spot of FLASH therapy, there could be a chance to use real-time PGI during a FLASH therapy treatment. Finally, the conclusions drawn in the present work will be tested in a future experimental campaign we are currently organizing.
We would like here also to highlight some limitations of our study. First of all, our simulations did not consider additional sources of neutron background, like the contribution of scattered neutrons from walls in the room. However, simulations performed including a 2 m thick spherical concrete wall with a 4 m inner radius did not show a significant loss of accuracy in range retrieval. Regarding the selected beam energy, a higher neutron background is expected for irradiation at higher ion energies, which could affect the capability to distinguish a PG fall-off. Therefore, future simulation studies will be necessary to explore the feasibility in the whole interval of therapeutic energies and, for the upper part of such interval, an additional effort to adopt gamma-neutron discrimination techniques might be useful and decisive.
Concerning the detector signal, we did not include the coupling between the scintillator and the readout system, which may reduce the detection efficiency, for instance due to count rate limits. Moreover, a 4 keV energy resolution was assumed for the detector, which is much better than the one of the real system, but having considered the events in an energy range of , we do not expect this to play a critical role. Finally, in this work we have considered only a homogeneous phantom. Therefore, further studies with inhomogeneities in the composition of the phantom should be carried out, in order to compare the results with the accuracy found in the homogeneous case. Regarding this topic, some ideas were elaborated for proton therapy, which could be inspiring for the same aim in CIRT (see [
In this work, we presented the first simulation study using the Monte Carlo code FLUKA to explore the feasibility of applying Prompt Gamma Imaging with a knife-slit camera to Carbon-ion Radiation Therapy. We took inspiration from previous studies, selecting the PG camera that was already successfully employed in proton therapy. We have shown that, for a typical treatment energy of C-ions, 150, a spot-by-spot range verification is unfeasible with the current set-up. However, combining the signal of various spots covering a single tumour layer, a fall-off precision of about 4 mm can be achieved. In the near future, measurements will be planned to check the simulation results found in this work. An experimental validation of the theoretical findings would open the possibility of using PGI monitoring in CIRT on a layer-by-layer basis.
CRediT authorship contribution statement
Andrea Missaglia: Performed Monte Carlo simulations and data analysis, Prepared the original draft. Aicha Bourkadi-Idrissi: Performed data analysis, Revised the original manuscript. Francesco Casamichiela: Contributed to setting up Monte Carlo Simulations, Prepared the original draft. Davide Mazzucconi: Contributed to setting up Monte Carlo Simulations, Revised the original manuscript. Marco Carminati: Conceptualization, revised the original manuscript. Stefano Agosteo: Contributed to setting up Monte Carlo Simulations, Revised the original manuscript. Carlo Fiorini: Conceptualization, Revised the original manuscript.
Stichelbaut F, Jongen Y. Verification of the proton beam position in the patient by the detection of prompt gamma-rays emission. In: 39th meeting of the particle therapy co-operative group, Vol. 16. 2003.
Peters A, Cee R, Haberer T, Hoffmann T, Reiter A, et al. Spill structure measurements at the Heidelberg ion therapy centre. In: European particle accelerator conference. 2008, p. 1824–6.
00031-5&id=gr1.jpg){kind=link}
00031-5&id=gr2.jpg){kind=link}
00031-5&id=gr3.jpg){kind=link}
00031-5&id=gr4.jpg){kind=link}
00031-5&id=gr5.jpg){kind=link}
00031-5&id=gr6.jpg){kind=link}
00031-5&id=gr7.jpg){kind=link}
00031-5&id=gr8.jpg){kind=link}
00031-5&id=gr9.jpg){kind=link}
00031-5&id=gr10.jpg){kind=link}
00031-5&id=gr11.jpg){kind=link}
00031-5&id=gr12.jpg){kind=link}
00031-5&id=gr13.jpg){kind=link}
00031-5&id=gr14.jpg){kind=link}
00031-5&id=gr15.jpg){kind=link}
00031-5&id=gr16.jpg){kind=link}
00031-5&id=gr17.jpg){kind=link}
00031-5&id=gr18.jpg){kind=link}