First full-beam PET acquisitions in proton therapy with a modular dual-head dedicated system

During particle therapy irradiation, positron emitters with half-lives ranging from 2 to 20 min are generated from nuclear processes. The half-lives are such that it is possible either to detect the positron signal in the treatment room using an in-beam positron emission tomography (PET) system, right after the irradiation, or to quickly transfer the patient to a close PET/CT scanner. Since the activity distribution is spatially correlated with the dose, it is possible to use PET imaging as an indirect method to assure the quality of the dose delivery. In this work, we present a new dedicated PET system able to operate in-beam. The PET apparatus consists in two 10 cm × 10 cm detector heads. Each detector is composed of four scintillating matrices of 23 × 23 LYSO crystals. The crystal size is 1.9 mm × 1.9 mm × 16 mm. Each scintillation matrix is read out independently with a modularized acquisition system. The distance between the two opposing detector heads was set to 20 cm. The system has very low dead time per detector area and a 3 ns coincidence window, which is capable to sustain high single count rates and to keep the random counts relatively low. This allows a new full-beam monitoring modality that includes data acquisition also while the beam is on. The PET system was tested during the irradiation at the CATANA (INFN, Catania, Italy) cyclotron-based proton therapy facility. Four acquisitions with different doses and dose rates were analysed. In all cases the random to total coincidences ratio was equal or less than 25%. For each measurement we estimated the accuracy and precision of the activity range on a set of voxel lines within an irradiated PMMA phantom. Results show that the inclusion of data acquired during the irradiation, referred to as beam-on data, improves both the precision and accuracy of the range measurement with respect to data acquired only after irradiation. Beam-on data alone are enough to give precisions better than 1 mm when at least 5 Gy are delivered.


Introduction
Dose range in proton therapy is subject to several sources of uncertainty, including the incompleteness of physical models (Pia et al 2010) and the effects of registration and patient movement (Hui et al 2008). These uncertainties worsen the already limited capabilities in predicting relative biological effectiveness for all tissues, especially those close to the Bragg peak. As a result, the use of beam angles that would place the distal edge close to critical structures is discouraged (Grassberger et al 2011). This partially neutralizes the advantage of charged particles with respect to photons.
One of the main actions that are being taken in order to overcome these problems is to research effective methods for in vivo verification of the beam delivery and, in particular, of the proton range in the patient. Since the sum of range uncertainties typically exceeds 2 mm, except in some cases, such as lung or deep-seated treatments, where it might be bigger (Paganetti 2012), an in vivo range verification tool has to provide both millimetric accuracy and precision in order to effectively contribute to the clinical practice. If effective, an in vivo monitoring would allow better treatment planning by assuring the detection of discrepancies between the planned and the actual delivery close to critical organs. Also, it could support experimental activities for better understanding uncertainties and their sources, thus contributing to increase the reliability of safety margins and eventually to enhance the number and diversity of beam directions.
The most used in vivo range verification method is positron emission tomography (PET) (Enghardt et al 2004, Knopf et al 2009, 2011. PET monitoring is based on the fact that during proton therapy several positron emitters are generated in the irradiated region (Bennett et al 1978, Litzenberg et al 1992. These emitters, e.g., 11 C, 13 N and 15 O, have half-lives that allow to examine the subject during and after the irradiation with a PET scan (Parodi et al 2002). The resulting images can then be used to measure indirectly the delivered dose (Del Guerra et al 1997, Parodi andEnghardt 2000). Another similar method exists, where the detected radiation comes from prompt gamma emission after nuclear excitation by the protons in tissue (Polf et al 2009a(Polf et al , 2009b. Prompt gamma detection has the main advantage of a more direct correlation between the photon signal and the dose range with respect to PET. However, it still suffers from low detector efficiencies, thus leaving PET as the only practical method at this time (Moteabbed et al 2011).
The PET signal results from inelastic nuclear reactions with energy thresholds ranging from 3 to 20 MeV. Therefore the distal activity fall-off is correlated, but not directly matched to the dose distribution. Thus, there is much interest in determining how the β + activation can be used to estimate the dose distribution (Oelfke et al 1996, Parodi and Bortfeld 2006, Attanasi et al 2008a, Knopf and Lomax 2013.
The established method is to use a Monte Carlo-simulated distribution of the positron emitters and compare the predicted image with the measured image (Parodi et al 2007a. Other methods use analytical forward or back-filters that transform the axial dose profile into the activity profile (Parodi and Bortfeld 2006) and vice versa (Attanasi et al 2008a(Attanasi et al , 2011. A successful estimation of the dose range from the activity depends on the accuracy of underlying cross section data  and requires in any case PET images of excellent quality (Remmele et al 2011, Aiello et al 2013. However, achieving high image quality is not always possible. PET imaging for proton therapy is subject to geometrical and temporal constraints that are absent in clinical PET applications (Crespo et al 2006, Nishio et al 2010. For systems placed around the immobilized patient in the beam line, i.e., in-beam, the angular coverage is limited by the free passage needed for the beam port. Conversely, fullring systems need to be placed out of the beam line, thus requiring a mechanical movement between the irradiation and the monitoring. This reintroduces registration uncertainties, partial waste of the emitted radiation and worsens the effects of biological washout. A systematic comparison of these modalities shows that off-beam monitoring is not recommended, except for some specific tumours or unless the PET system is placed in the same treatment room (i.e., in-room) (Shakirin et al 2011).
In-beam monitoring is the first clinical implementation of PET for monitoring of radiotherapy. It was launched in 1997 at the experimental carbon ion therapy facility at GSI Helmholtzzentrum für Schwerionenforschung (centre for heavy ions research), Darmstadt, Germany (Enghardt et al 2004). The PET installation is completely integrated into the therapy facility and has the potential to provide the best quality feedback for treatment monitoring, since it requires no additional time to move the patient to the PET scanner. The ability to detect alignment errors is lower than in other modalities. However, the artefacts due to incomplete angular coverage can be mitigated, at least with respect to the beam direction, by increasing the detectors size and reducing the gaps between the two heads (Crespo et al 2006). Additionally, counting statistics can be improved if the acquisition is performed also during the dose delivery.
The main technological problem consists in filtering out data acquired during particles extraction, because the strong beam-induced background noise floods the PET detectors (Enghardt et al 2004. Such noise might originate from the decay of β + -emitters with half-lives in the millisecond range (e.g., 8 B, 9 C and 12 N) with correspondingly high β + -endpoint energies, positron generating processes accompanying nuclear reactions and random coincidences with the γ -decay of excited nuclear levels (Pawelke et al 1997). Therefore, the effect is either a broadening of the β + activity spatial distribution or a paralysis of the acquisition system due to excessive random coincidence data that do not carry any information regarding the β + distribution (Enghardt et al 1999, Parodi et al 2002. Thus, state of the art in-beam PET solutions are only applicable to synchrotron-based facilities and their effectiveness is limited by the accelerator duty cycle (Nishikido et al 2010, Tashima et al 2012, Nishio et al 2010, Shakirin et al 2011. Based on previous experiences with the first dedicated in-beam DoPET system (Vecchio et al 2009), we developed a new in-beam PET, with wider detectors and modularized acquisition for better counting performances (Sportelli et al 2011a). In the following we present the first reconstructed PET images of the activity produced during the continuous irradiation of a polymethyl methacrylate (PMMA) phantom at the CATANA cyclotron-based proton therapy facility. The new architecture and the improved dead time at the front-end (Belcari et al 2007) allowed to mitigate the effects of random coincidence rates and to acquire valid PET data in the high noisy environment produced by the irradiation (Sportelli et al 2012).
The paper is organized as follows: in section 2, the distribution of PET activity in time is described to give an idea of the trade-off between acquisition duration and available statistics. The PET system, the experimental set-up and the software algorithms used for image reconstruction and range determination are described in section 3. Section 4 reports the obtained results: counting statistics, energy spectra, in-beam images and activity ranges. In section 5 we discuss and conclude the work.

PET activity and range verification quality
Positron emitters production rate during irradiation is proportional to the current of the accelerator, i.e., to the treatment dose rate. The production of radioisotopes competes with the decay process. For most of the observed radioisotopes the irradiation condition is far from the secular equilibrium and the result is a continuous increase on the β + activity. After the irradiation, the β + decay is the only process occurring and the expected behaviour of the activity over time is a sum of exponential decays which is the convolution of the contribution of the various isotopes present in different quantities and decaying at different speeds.
We are interested in quantifying the amount of statistics available from most abundant isotopes produced by a continuous irradiation (beam-on) as compared to those available after the end of it (beam-off). To do so we assume that the PET detectors have sufficiently low dead time, in order not to paralyze during the irradiation because of the high background noise. In this kind of comparison, the proportion between beam-on and beam-off data is independent from the dose rate, i.e., it only depends on the half-life of each produced isotope and the involved time intervals. For simplicity let us consider two main isotopes produced by protons on PMMA, i.e., 15 O and 11 C, in a likely scenario in which the irradiation lasts 3 min and the acquisition 10 min. Assuming that the production is a constant process, we can describe these behaviours in arbitrary units for each separated isotope. To give an idea of the proportions between beam-on and beam-off available data, we report in figure 1 (left) the activity rates calculated as the result of the production-decay processes, with the corresponding integral contribution of data acquired during and after the irradiation (right).
Since one of the main limiting factors in state of the art PET for in vivo monitoring is the low statistics, we aim at demonstrating that in-beam data, taken while the beam is on, can effectively contribute to improve the statistical accuracy of the reconstructed image.
In the following discussion we use three dataset types: we refer to the data acquired during the irradiation as 'beam-on' data; all the remaining data is referred to as 'beam-off' data. We call 'full-beam' data the union of both datasets for a given acquisition. In order to characterize the trade-off between acquisition duration and reconstruction noise, beam-off data are reconstructed at 60, 180, 300 or 600 s after the end of the irradiation. To evaluate the impact of beam-on data in range verification, we will compare the activity ranges calculated on the reconstructed images and their statistical precision for a given monoenergetic irradiation and for each dataset.

The PET system
A planar positron imaging system was arranged at the beam line at CATANA proton therapy facility (INFN, Catania, Italy) (Cirrone et al 2004) as in previous experiments (Attanasi et al 2008b). In comparison to the previous system (Vecchio et al 2009), the two detector heads were increased four-fold, i.e., from 5 cm × 5 cm to 10 cm × 10 cm . The number of detector modules per detector head was increased from 1 to 4. The gap between each module is 6 mm. Each module is composed of 23 × 23 LYSO crystals of 1.9 mm × 1.9 mm × 16 mm. The energy response of the used LYSO crystals is characterized by a high-energy tail due to optical effects (Bonifacio et al 2010). The PET system was installed on the patient chair and the centre of its detection area was aligned with the isocentre of the beam line in the treatment room (Sportelli et al 2012. The distance between the two opposing detector heads is adjustable from 7 to 30 cm. The field of view (FOV) is 10 cm × 10 cm × D cm, where D is the distance between the two detector heads. The point source resolution at the centre is 1.2 mm × 1.2 mm × 6.0 mm full width at half maximum.
The acquisition system features a custom modular board able to handle up to 18 modules (Sportelli et al 2011a). Coincidence detection is performed by a synchronous processor operating at 240 MHz and providing a coincidence resolution of 3 ns (Sportelli et al 2011b). The maximum data collection rate for detected coincidences is about 1 million counts per second (cps), though this limit has been never reached with eight modules, because of the counting limits of single photons at the front-end. The information of the on-off time points of beam irradiation is currently inferred from the count rates and double checked with the beam monitors with a resolution of 1 s. For future acquisitions it is planned to record the beam status from an external signal synchronized with the accelerator.
The detection efficiency was calibrated by using a 11 cm × 11 cm × 0.3 cm planar plastic container filled with 18 F-fluorodeoxyglucose (FDG) placed at the middle of the FOV and parallel to the detectors surface. The calibration is used for a correction of the imaging uniformity and the detection sensitivity. The maximum coincidence counting rate observed with the planar source filled with 26 MBq of FDG is 300 kcps.

Experimental set-up
Eye melanoma treatments at CATANA use typically 62 MeV proton beams with dose rates that vary from 10 to 20 Gy min −1 (Cuttone et al 2011). Our measurements were performed with monoenergetic passively collimated proton beams of 62 MeV. Only pristine Bragg peaks were used, which are characterized by a peak-to-plateau ratio of 5. We did not apply any range modulation in order to simplify the analysis at the activity distal endpoint. This differs from the clinical practice, where a spread out Bragg peak (SOBP) obtained with a PMMA modulator wheel is generally used to cover the whole extension of the tumour in depth with the prescribed dose. In our experiment variable accelerator currents were delivered on the PMMA phantom with dose rates at the low end of the typical values, i.e., ranging from 0.7 to 9.7 Gy min −1 . A homogeneous PMMA block phantom was placed at the centre of the FOV and irradiated with a beam of 34 mm diameter. The distance between the two heads was set at 20 cm. A list of the four performed irradiations is reported in table 1. The value of the absolute dose corresponds to the dose at the entrance of the phantom.
Each measurement was performed sequentially on initially inactive phantoms of PMMA. The dose and the beam current were preset at the beginning of the irradiation. The phantom was irradiated continuously until reaching the preset dose in the beam monitor. The dose rate was calculated as the total dose divided by the irradiation time. Activity position and intensity were measured during the irradiation plus 600 s immediately after the end of the irradiation.
Coincidence data were stored in a list mode format. Single, coincident and random counts rates were recorded in a separate file with binning intervals of one second. Random coincidence rates were measured, but their spatial distribution was not acquired.

Image reconstruction
Images were reconstructed using an iterative Maximum Likelihood Expectation Maximization (MLEM) reconstruction algorithm. The implemented algorithm makes use of a system model computed using a multi-ray variant of the Siddon algorithm (Moehrs et al 2008). Images and profiles used in this work were all obtained by performing five MLEM iterations. Each iteration took 5 s to be computed on a 64 bit Intel i5-2400 based machine with four cores running at 3.1 GHz. There are also plans for performing the same reconstruction on a GPUbased implementation (Sportelli et al 2013). Data were reconstructed using a LOR based data structure and applying a 350-850 keV energy window. The size of the reconstructed image is 100 × 100 × 100 voxels, 1 mm × 1 mm × 1 mm each, thus including the complete FOV up to 5 cm from the detectors surface. At present, no random or attenuation correction was applied. A problem in estimating the activity distal range on voxels lines is that such lines have much lower statistics with respect to macroscopic profiles over larger regions. The low statistics can cause spikes that make it difficult to properly determine the threshold level. Our approach has been to first detect the proximal and distal edges of each line as the two longest monotonic paths in the profile. In doing so, any relatively small spike with respect to the detected path was removed (figure 2). The position of each edge was then determined as the threshold level over the local path maximum. The proper choice of the range verification position in the fall-off region in the clinical practice can be controversial, because ranges obtained with high thresholds are affected by the width of the SOBP while those obtained with low thresholds are generally more sensitive to image noise and to the beam energy spread (Parodi et al 2007a). Thus, thresholds ranging from 20% to 50% are generally used (Knopf et al 2008). For our purposes a 50% threshold value is a reliable choice because in this study we have do not use SOBP. Also, lower thresholds would be more problematic due to the low statistics of the datasets involved in this study.

Estimation of the activity distal fall-off
For irradiations performed with monoenergetic beams, like in the present set-up, one can suppose that the proton range is the same for each profile and then a constant depth for the distal activity fall-off can be expected. Thus, the precision of the range measurement was estimated by calculating the standard deviation of the measurements performed in each profile (σ depth ). The accuracy of the measurement of the distal fall-off of the activity profile can be estimated as the deviation of the average depth (μ depth,t ), as calculated from all the 1D profile measurements with data acquired up to a time t, from the estimation of the average depth, calculated in the same way, but using the whole statistics of 10 min of PET measurement with beam-off (μ depth,600s ). The latter is considered here as the best estimation, for each irradiation, of the depth of the activity distal fall-off, being the measurement with the lesser statistical error.

Comparison of counting statistics
Total and random coincidence rates for the four acquisitions are shown in figure 3. During the irradiation the count rates for both true and random coincidences increase abruptly, indicating the presence of fast nuclear processes. The random to total coincidence ratio decreases during the irradiation. The maximum observed random to total coincidence ratio is 24.5%, at the beginning of the 9.7 Gy min −1 irradiation. The acquisition rates at the beginning, at the end and right after the irradiation are reported in table 2.

Energy distribution of the acquired data
Energy histograms of the acquired coincidences are reported in figure 4. The black line corresponds to the coincidences acquired during the irradiation (beam-on). The grey line corresponds to the coincidences acquired after it (beam-off). Dashed lines mark the 350-850 keV energy window used for image reconstruction. The spectrum of beam-off data  remains almost unchanged with the different dose rates, except for the 307 keV intrinsic emission peak of 176 Lu, which is more evident at low dose rates, as expected. A high-energy tail is present in both beam-off and beam-on data due to optical attenuation in the used scintillating crystals (Bonifacio et al 2010). This tail goes from the 511 keV peak to roughly 1 MeV. With higher dose rates, beam-on data show an increasing degradation of the annihilation photo-peak and a generally higher base level at all energies. This behaviour was also expected due to the radiation noise during the delivery. The fraction of events within the selected energy window is reported in table 3. While the beam is on, the fraction of events with energies >850 keV is much higher and is increasing with the dose rate.

Distal fall-off during and after the irradiation
The reconstructed images for acquisition 1 (2 Gy) and 4 (20 Gy) are shown in figures 5 and 6, respectively, including beam-on only data (a), 60 s of beam-off data (b), 60 s of full-beam data (c) and 600 s of full-beam data (d). The ROIs used for range determination are marked with a dashed line. The ROIs are cylindrical with 2 cm diameter and 9 cm height. Profiles calculated on such images are shown in figure 7(a)-(d), i.e., 2 Gy, and compared with images of the fourth acquisition (e)-(h), i.e., 20 Gy. The dose range is reported as a black dashed line in each plot. Dose range has been calculated given the depth of the Bragg peak in homogeneous PMMA, i.e., 29 mm (Oelfke et al 1996), starting from the activity rising edge. In figure 7 it can be seen how the image background noise is higher when the reconstruction is performed with low statistics, i.e., low dose (a)-(d), rather than with high dose rates (e)-(h). The noise level in the image is reflected in the distribution of the calculated distal ranges. figure 8 shows the variation of the average depth and the relative standard deviation of the measurement of the distal fall-off related to the activity profile for each of the four acquisitions here considered. For each acquisition μ depth and σ depth are reported as estimated from beam-on data only (t = 0, dotted line), beam-off data only (at t = 60 s, t = 180 s, t = 300 s and t = 600 s), and full-beam data, including both beam-on and beam-off data. Values for μ depth,t and σ depth,t are also reported in table 4. Plots show that the value of μ depth,t derived with beam-off data only (grey line) is slightly underestimated with respect to the reference value of μ depth,600s (highest statistics) when the acquisition time is short, while it converges to the latter as the time increases. The μ depth,0s value (beam-on only) is overestimated with respect to μ depth,600s in all cases but the 2 Gy irradiation. The inclusion of beam-on data (full-beam) always improves the convergence speed. As for image noise, a higher standard deviation is obtained for low doses and thus a less precise range is determined.
In each of the four acquisitions, the largest mean range deviation from the reference value is seen in the images reconstructed with only 60 s of beam-off data. The same effect has been observed in all the acquisitions with less than one million events (figure 9), thus suggesting a systematic range determination error for very low statistics.    ((a), (e)), first 60 s of beam-off data ((b), (f)), beam-on data plus 60 s of beam-off data ((c), (g)), beam-on data plus 600 s of beam-off data ((d), (h)). The plot is obtained by superimposing the profiles of each voxel line included in the ROI (see figures 5 and 6). The white dashed line is the mean of all voxel lines. The z coordinate of the dose range is indicated with a black dashed line. Dose range has been calculated given the depth of the Bragg peak in homogeneous PMMA, i.e., 29 mm, starting from the activity rising edge.

Conclusions
We have built and optimized a modular dual-head PET system for treatment monitoring in hadron therapy. The system can operate continuously in-beam and is able to successfully acquire the produced β + activity in a field-of-view of 10 cm × 10 cm × 10 cm also during Range vs statistics 20 Gy -9.7 Gy/min, Beam off only 10 Gy -3.0 Gy/min, Beam off only 5 Gy -1.6 Gy/min, Beam off only 2 Gy -0.7 Gy/min, Beam off only the irradiation. A distinction has been introduced between beam-on, beam-off and full-beam data, referring to the acquisition output during the irradiation, after or both, respectively. The question that has been addressed is if beam-on data can be used, either alone or jointly with beam-off data, to improve the quality of the system response for monitoring purposes. A set of PET measurements was then performed at the cyclotron-based proton therapy facility at CATANA, during the irradiation of a PMMA phantom at different dose rates and stopping the acquisition 10 min after. As already observed in other studies, a characteristic high random count rate appears during the irradiation. This count rate is roughly constant with a continuous beam and increases quickly with the dose rate. Since the phantom activation increases during the irradiation, the randoms to total coincidences fraction decreases with time and becomes abruptly negligible right after the beam turns off. The maximum observed fraction at the beginning of an irradiation with 9.7 Gy min −1 dose rate is 24.5%. This value is of the same order of an acceptable randoms fraction in clinical PET imaging. However, the noise structure is different from random coincidence noise in conventional PET, thus future studies are expected to better characterize its nature with respect to useful β + activity and to improve the techniques used for reconstructing such data.
The coincidence energy spectra of beam-on and beam-off data were also compared. A small flat contamination from radiation with energies up to 2 MeV was observed in beam-on data (2 MeV is the dynamic range of the system). With dose rates up to about 10 Gy min −1 no significant degradation of the energy spectra was observed, being the percentage of acquired events within the 350-850 keV energy window up to 62.2% of the total.
For each measurement and modality we estimated accuracy and precision by studying the distribution of the activity ranges of a set of voxel lines within the activated material. Results show that the inclusion of beam-on data improves both the precision and accuracy of the proton range measurement with respect to data acquired after the irradiation only. The use of beam-on data is enough to give a precision better that 1 mm (σ ) when at least 5 Gy are delivered. Alternatively, beam-on data can be effectively used to increase the available statistics for image reconstruction when the dose is lower. For the lowest available dose rate (2 Gy min −1 ), beam-on data did not produce sub-millimetric range precision, but the precision of the first 60 s of full-beam data is half with respect to that obtained with beam-off data only.
The PET reconstructed signal is not quantitative. No developments were carried out to study its capability as a quantitative instrument, but with proper calibration we expect its response to be linear within a wide range of input rates. There are no reasons a priori that would hamper its ability to give quantitative results more than in conventional PET systems.
The results of this investigation are promising for improving the quality and promptness of PET monitoring, though there is still a big margin for improvement. The proposed in-beam PET system has in principle the same geometrical limitations as previous in-beam solutions, but uses new electronics that allow to reduce the acquisition duration. It is compatible with Monte Carlo based simulations as well as with all the range verification techniques that have been used so far. However, it requires the detectors response to be included in the simulation model so as to take into account the artefacts due to the limited angular coverage. Washout effects are expected to be less prominent than in any other PET imaging technique, since the acquisition lasts only a few minutes from the beginning of the irradiation.
Future work is planned to increase the angular coverage of the system, in order to improve random coincidences management, to refine system synchronization with the beam and to characterize the detectors response. Extending the detectors area from 10 cm × 10 cm to 15 cm × 15 cm would require no major changes. The number of detector modules would increase from 8 to 18 and the software analysis would be more complex, but the overall acquisition architecture and implementation would remain the same.