Feasibility of quasi-prompt PET-based range verification in proton therapy

Compared to photon therapy, proton therapy allows a better conformation of the dose to the tumor volume with reduced radiation dose to co-irradiated tissues. In vivo verification techniques including positron emission tomography (PET) have been proposed as quality assurance tools to mitigate proton range uncertainties. Detection of differences between planned and actual dose delivery on a short timescale provides a fast trigger for corrective actions. Conventional PET-based imaging of 15O (T1/2 = 2 min) and 11C (T1/2 = 20 min) distributions precludes such immediate feedback. We here present a demonstration of near real-time range verification by means of PET imaging of 12N (T1/2 = 11 ms). PMMA and graphite targets were irradiated with a 150 MeV proton pencil beam consisting of a series of pulses of 10 ms beam-on and 90 ms beam-off. Two modules of a modified Siemens Biograph mCT PET scanner (21 × 21 cm2 each), installed 25 cm apart, were used to image the beam-induced PET activity during the beam-off periods. The modifications enable the detectors to be switched off during the beam-on periods. 12N images were reconstructed using planar tomography. Using a 1D projection of the 2D reconstructed 12N image, the activity range was obtained from a fit of the activity profile with a sigmoid function. Range shifts due to modified target configurations were assessed for multiples of the clinically relevant 108 protons per pulse (approximately equal to the highest intensity spots in the pencil beam scanning delivery of a dose of 1 Gy over a cubic 1 l volume). The standard deviation of the activity range, determined from 30 datasets obtained from three irradiations on PMMA and graphite targets, was found to be 2.5 and 2.6 mm (1σ) with 108 protons per pulse and 0.9 and 0.8 mm (1σ) with 109 protons per pulse. Analytical extrapolation of the results from this study shows that using a scanner with a solid angle coverage of 57%, with optimized detector switching and spot delivery times much smaller than the 12N half-life, an activity range measurement precision of 2.0 mm (1σ) and 1.3 mm (1σ) within 50 ms into an irradiation with 4 × 107 and 108 protons per pencil beam spot can be potentially realized. Aggregated imaging of neighboring spots or, if possible, increasing the number of protons for a few probe beam spots will enable the realization of higher precision range measurement.


Introduction
Compared to radiotherapy with photons, proton therapy offers the benefit of delivering a more conformal dose to the tumor volume while reducing the dose to co-irradiated normal tissues. This benefit is derived from the finite range of protons and the high dose deposition at the end of the proton trajectory. This in principle enables accurate targeting of the tumor. However, uncertainties in the proton range, due to factors including patient immobilization and setup errors, patient and organ motion, imaging artifacts and uncertainties in the conversion of CT numbers to relative stopping power can lead to deviations of the delivered dose distribution with respect to the planned one, causing inadvertent exposure of organs at risk and/or under-dosage of the tumor volume (Goitein 1985, Lomax 2008, Paganetti 2012. Thus, in clinical practice, treatment plans are designed to be robust against these uncertainties through addition of safety margins (van Herk et al 2000, Albertini et al 2011, Paganetti 2012. Alternatively, robust optimization can be used to design plans which are robust with respect to typical uncertainty scenarios without the need for an additional safety margins (Unkelbach et al 2007, Pflugfelder et al 2008, Fredriksson et al 2011, Li et al 2015, Inoue et al 2016and review by, Unkelbach et al 2018. Depending on the magnitude of the considered range uncertainty scenarios, the adoption of these robust planning techniques may still lead to the delivery of extra dose to normal tissues and thereby compromise the reduction of radiation-induced side effects which is the main rationale for proton therapy.
The essentially 3D nature of the dose deposition by a particle beam makes in vivo treatment verification more important than for highly penetrating MV x-ray photon beams, whose dose deposition is basically 2D in nature. With the increasing complexity of treatment techniques employed in photon radiotherapy, the use of an Electronic Portal Imaging Device (EPID) for in vivo treatment verification has gained considerable traction (Meertens et al 1985and review by, van Elmpt et al 2008, Mccurdy et al 2017, Mijnheer 2019. While these devices measure the transmission of the MV x-rays, the finite range of particle beams precludes this method in particle beam radiotherapy. A number of in vivo range verification techniques for particle beam radiotherapy have been proposed, aiming at the reduction of range uncertainties and thereby the required safety margins (See reviews by Parodi 2011, Knopf and Lomax 2013, Parodi and Polf 2018, potentially resulting in a reduction of normal tissue complications. When adopting robust planning techniques, in vivo range verification technique can allow a reduction of the magnitude of the range uncertainty as well as the dose to surrounding normal tissues. In vivo verification techniques rely mostly on the detection of secondary radiation resulting from the interaction of the particle beam with the tissues: positron annihilation photons stemming from the decay of radioactive nuclides (Maccabee et al 1969, Paans and Schippers 1993, Nishio et al 2008, Enghardt et al 2004 review papers by Studenski and Xiao 2010, Fiedler et al 2012, Zhu and El Fakhri 2013, prompt gamma photons emitted by excited nuclei (Min et al 2006, Polf et al 2009, Testa et al 2008 (Amaldi et al 2010, Henriquet et al 2012, Gwosch et al 2013, Rucinski et al 2018. Beam-induced acoustic waves (Hayakawa et al 1995, Patch et al 2016, Kellnberger et al 2016, Lehrack et al 2017, Jones et al 2018 are also being investigated for this purpose.
Positron emission tomography (PET) verification is based on the imaging of positron emitters which are created through nuclear reactions between the particle beam and the nuclei in the traversed tissues. PET imaging can be performed using the following implementation strategies: during the irradiation (in-beam) (Enghardt et al 2004, Tashima et al 2016, Yoshida et al 2017, Pennazio et al 2018, inside the treatment room (in-room) (Nishio et al 2010, Zhu et al 2011 and outside the treatment room (offline) (Parodi et al 2007, Nishio et al 2008, Knopf and Lomax 2013, Nischwitz et al 2015. Preclinical and clinical studies, which mostly image the longer-lived positron emitters (with half-lives from 2 to 20 min), have demonstrated promising results for in-vivo verification with PET (Parodi et al 2005, 2007, Zhu et al 2011, Knopf and Lomax 2013. However, the imaging of these long-lived positron emitters has a limited clinical value due to factors related to their physical half-lives and biological washout. The relatively long half-lives imply that long scanning times are required to gather sufficient counts, precluding feedback during the irradiation itself. Furthermore, the precision of range verification based on the imaging of these nuclides is limited due to their higher susceptibility to biological washout which leads to a redistribution of activity and reduction in the image quality.
The in-beam PET implementation strategy provides an attractive option towards improving the verification accuracy and feedback-time. This strategy is the least influenced by biological washout. When imaging the very short-lived positron emitter 12 N (T 1/2 = 11 ms), there is no loss due to biological washout. In addition to reduced impact of biological washout, the strategy allows maximum counting statistics in a short imaging time related to the half-lives of the decaying positron emitters. During an in-beam PET measurement, a wide range of nuclides, with half-lives up to about 20 min, will contribute to the detected PET activity. PET-based fast feedback on the proton range can only be achieved through imaging of the very short-lived positron emitters. The provision of such prompt feedback will provide a trigger for implementation of actions, such as daily adaptive proton therapy (Albertini et al 2019), to correct deviations from the treatment plan. Beyond conventional fractionated proton therapy, the implementation of hypofractionation schemes (Grewal et al 2019) and high dose rate treatment techniques such as flash irradiations (Favaudon et al 2014, Vozenin et al 2019 might benefit from fast feedback at short/ultrashort timescale using a probe fraction of the full dose. The most abundantly produced very short-lived nuclide in proton therapy, 12 N, has been shown to dominate the beam-induced PET decay in carbon-rich tissues (Dendooven et al , 2019 in the early stage of dose delivery. A proof-of-principle study on the use of this nuclide for near real-time proton range verification using a small PET scanner (6.5 × 6.5 cm 2 modules) shows that a precision in the range measurement of 3 mm is obtained for the detection of 4000 12 N PET counts during irradiations with 2.5 × 10 10 protons delivered over a 120 s period (Buitenhuis et al 2017). For clinical implementation of 12 N-based range verification in pencil beam scanned delivery of proton therapy, it is highly desirable to monitor the range with ≤ 2 mm precision using about 10 8 protons delivered within a few tens of milliseconds into the irradiation. In the present study, we report on the precision of range measurement for the delivery of single spots by imaging the 12 N activity with a larger scanner at clinical beam fluences (10 8 protons per spot).

Irradiation setup
The experiment was performed at the AGOR cyclotron of the KVI-Center for Advanced Radiation Technology (KVI-CART), University of Groningen. Figure 1 shows the experimental setup. A beam of 150 MeV protons was delivered through a horizontal beam line to the irradiation cave. The intensity of the beam was monitored with an air-filled ionization chamber (beam intensity monitor (BIM)) placed after the exit window of the beam line. The number of protons corresponding to a BIM monitor unit (MU) was determined via a calibration of the BIM, at low beam intensity, with a scintillation detector which counts individual protons traversing the BIM. The width of the beam at the target position was 7 mm FWHM as measured using a harp-type (wire grid) beam profile measurement system. Two target materials, graphite (ρ = 1.72 g cm −3 ) and PMMA (ρ = 1.19 g cm −3 , elemental composition: Carbon (33 at%), Oxygen (13 at%), Hydrogen (51 at%)), were irradiated. Table 1 gives a summary of the positron emitters produced during the irradiation of these targets. The dimensions of the graphite and PMMA targets (width × height × length) were 120 × 120 × 150 mm 3 and 120 × 120 × 200 mm 3 respectively, with the long side parallel to the beam direction. The targets fully cover the range of 150 MeV protons: 104 mm in graphite and 136 mm in PMMA (PSTAR 2019). The graphite and PMMA targets were mounted such that the center of the field of view (FoV) of the scanner corresponded to a depth of 100 mm and 130 mm in the two targets, respectively. The central axes of the target and the scanner FoV coincided with the beam direction. The distance between the centers of the two scanner modules was 25.2 cm.
A recent upgrade of the AGOR irradiation facility enables spot scanning. Dipole magnets for horizontal and vertical beam deflection are located about 1.95 m upstream from the target position. This functionality was applied to simulate a broader irradiation field, as used in clinical irradiations. The spot locations were calibrated with a system consisting of a CCD camera and scintillation screen. The maximum deflection used at the target position was ± 29 mm. The scanning magnets were controlled by a dual-channel waveform generator.
The beam consisted of pulses of 10 ms beam on and 90 ms beam off. The beam pulsing was realized through a pulse generator controlling the voltage on electrostatic deflection plates installed in the injection Table 1. Properties of positron emitters produced during proton irradiation of PMMA and graphite targets.The data on the production channel and threshold energy were retrieved from NNDC (2020).

Nuclide
Half - line of the cyclotron. Switching the beam on was faster than 0.1 ms; when switching the beam off, it took 3.0 ms to reach the 3% beam intensity level.

Target configurations and irradiations
The following target configurations were investigated: PMMA and graphite targets in the reference position; shifts of the targets in the downstream beam direction by 3 mm and 6 mm and degrading of the range in the target using 4.7 mm and 9.8 mm PMMA blocks affixed to the proximal surface of the targets. In each configuration, the targets were irradiated separately with 10 8 and 10 9 protons per pulse of 10 ms for a total irradiation time of 60 s (600 pulses). Data was acquired during the 90 ms pauses between the beam pulses and for an additional 70 s after the end of the irradiation. During the post-irradiation data acquisition, the detector pulsing (described in section 2.3) was still active. The changes in the proton range due to the modifications with respect to the reference configuration were determined using the differences in the measured 1-dimensional 12 N activity profiles of each target configuration relative to that of the reference configuration.
Next to homogenous target configurations, the detection of range shifts in a mixed range target configuration was investigated. A 9.8 mm PMMA range shifter was affixed to the PMMA block such that the range in the target was reduced for the upper half of the irradiation field. The beam was scanned vertically over the target height with 12.9 mm distance between spots.
It is important to mention that each irradiation was performed with a 'fresh' target. Targets previously irradiated were reused after a waiting period of at least 200 min (10 half-lives of 11 C). This criterion ensured that residual activities in the target were reduced to no more than 0.1% of the activity at the end of the irradiation.

PET scanner description
A modified Siemens Biograph mCT clinical scanner (Karlberg et al 2016) was used to image the beam-induced activity. The scanner consists of two opposing panels of block detectors (a full-ring scanner contains 12 such panels). Each panel measures 21 × 21 cm 2 and is composed of a 4 × 4 array of block detectors. A block detector comprises a 13 × 13 array of 4 × 4 × 20 mm 3 LSO scintillation crystals read out by four photomultiplier tubes (PMTs). The energy signals are transmitted through Cat 6A twisted pair cables with RJ-45 jacks to two Detector Electronics Assemblies (DEAs) which encode position, energy and time of arrival of the photons. A coincidence unit receives the processed signals from the DEAs, determines valid coincidence events and transmits the data to a data acquisition computer. A coincidence window of 4 ns and an energy discrimination window of 435-650 keV were used, respectively. The panels were installed such that they are curved around the vertical direction.
Because of the high radiation levels near the proton beam, two modifications to the Biograph mCT hardware were implemented. To prevent radiation damage to the electronics and data acquisition computer, the standard short signal cables were replaced by 7 m long cables (Ethernet Shielded Twisted Pairs Cat 6A). This enabled the installation of the electronics and data acquisition computer in the basement below the irradiation cave. A block detector from a Siemens Biograph scanner, very similar to the unmodified version of the detectors used in the present work, was tested in-beam with 150 MeV protons at the AGOR cyclotron of KVI-CART by Hueso-González et al (2015). At a distance of 30 cm from a proton beam stopped in a graphite target, a detector count rate of 10 5 s −1 was measured for a beam intensity of ∼10 pA. Scaling this to the situation of the present experiment (10 8 protons in a 10 ms pulse and a block detector at 12.5 cm from the beam) results in a count rate of about 10 8 s −1 , which is more than the detector and electronics can handle. For this reason, the PMT voltage dividers were modified such that the detectors can be effectively switched off during the beam on periods. The detector pulsing is controlled by a TTL signal synchronized with the beam pulsing. An oscilloscope screenshot of the PMT signal of a block detector demonstrating the switching of the PMT during tests with beam and radioactive sources is shown in figure 2(left) and 2(right), respectively. The PMT signal vs time observed during both tests is similar. After the detector is switched on, a period of 300 µs is required to reach the operational state, while a shorter time of 130 µs is required to switch the PMT off. The short period of higher gain after switching on the detector, seen in both online and offline tests, is an understood feature of the way the PMT is pulsed. Due to the tail in beam intensity when switching off the beam, a 3 ms delay between the beam-off signal and the detector-on signal was introduced as a compromise between exposure of the detector to the prompt beam-related radiation flux and fast recovery of the detector after switching it on. The measurements in coincidence mode with the scanner, however, show a much slower recovery of the coincidence rate (see section 3.2).
To understand the cause of this slow recovery of the coincidence rate, we investigated the stability of the energy signal when pulsing the detectors. As the Siemens scanner does not provide energy information in the list-mode data, a different acquisition system, a mesytec MDPP-16 time and amplitude digitizer, was used for this investigation. Two block detectors from the Siemens scanners and a pulsing of 10 ms off/90 ms on were used. The energy spectrum was calibrated using a 22 Na source. Coincidence data were acquired using a 68 Ge source. A 2D plot of the energy spectrum as a function of time for coincidence acquisition between the two detectors is shown in figure 3 (top). The data show that after switching on the PMT it takes about 30 ms for the gain to recover to a stable value. The red broken lines indicate the energy window of 435-650 keV; the black broken line indicates an energy of 460 keV. The energy window for selection of events as implemented in the scanner is based on the pulse height. The decrease in PMT gain results in a lower number of counts within the window. In figure 3 (bottom) the count rate recovery as observed during the post-irradiation (10 8 and 10 9 protons per pulse) data acquisition in the actual experiment and that of the test experiment are displayed. The post-irradiation data shows the relative counts obtained from the sum of 200 detector-on periods starting 50 s after the end of the irradiation when only long-lived nuclides are observed. The PET count in this period results from the decay of 10 C, 15 O and 11 C which have half-lives that are much longer than 87 ms. Thus, a constant count rate is expected in the detector-on period of 87 ms. Similarly, the 68 Ge source has a constant activity, so a constant count rate vs time is expected (dashed line). However, the count rate is seen to increase towards a stable value in the first 30 ms, as in the actual experiment. To obtain a good match between the recovery observed during the actual experiment and in the test experiment the lower threshold of the energy window had to be increased from the nominal 435 keV to 460 keV. This 6% difference is very likely due to a small difference in the energy calibrations (gain and/or offset) of the Siemens scanner and the mesytec MDPP-16 system. The agreement of the measurements with beam-induced activity and callibration sources thus shows that the recovery of the coincidence count rate over a 30 ms period after switching on the detectors is the same regardless of the count rates, irradiation condition and can be explained by the PMT gain stabilization taking about 30 ms.
A recovery of PMT gain when used in gated mode has been observed in other studies (Hara et al 2013, Etile et al 2018. Furthermore, the same phenomenon has been observed when using PMT-based detectors for prompt gamma timing (PMT) range verification (Werner et al 2019). Given the need to acquire data as soon as possible, two solutions can be envisioned to correct for the time independent detector response. Firstly, a correction factor can be applied on an event-by-event basis to shift the measured energy as a function of time. Such corrections on the energy spectrum are implemented in Werner et al (2019)   height to account for the gain restoration time of a gated PMT. Alternatively, the measured counts can be corrected using normalization factors (relative count rates) accounting for the time variation of the count rate for an otherwise constant intensity source. Such a normalization factor can be obtained from the relative count rates of scanner measurements of beam-induced long-lived nuclides. The second approach is used for the correction applied in section 3.2 since we do not have access to the energy information of the individual events.

Scanner sensitivity measurement
The absolute system sensitivity in the center of the FoV was measured with a 1.9 MBq 68 Ge point source (∅ = 1 mm embedded in a plastic disk of ∅ = 25.4 mm and thickness of 5.3 mm). The absolute sensitivity was calculated as the ratio between the coincidence counting rate and the source activity. For measurements of the sensitivity across the FoV, identical 7.9 MBq 68 Ge line sources (active ∅ = 1.52 mm and length 182 mm enclosed in a steel rod) were used. Two such sources were connected lengthwise to make a longer line source covering the entire scanner FoV. Two separate measurements of this longer line source at the midplane between the two panels were done: one with the source along the horizontal axis of the scanner, i.e. in the beam direction, and one with the source along the vertical axis of the scanner, i.e. perpendicular to the beam direction. The sensitivity profile was obtained as the ratio between the true coincidence counts of each image pixel and the coincidence count in the center of the FoV. The 2D sensitivity map can be constructed as a product of the normalized horizontal and vertical sensitivity profile.

Data analysis
List-mode data files were recorded containing the coincidence crystal IDs and difference in detection time between the two block detectors involved. Time tags were inserted into the data stream at 1 ms intervals. The time tag values are relative to the beginning of data acquisition. There was no synchronization between the beam pulsing and the beginning of data acquisition. For each measurement, data acquisition was started a few seconds before the irradiation. For implementation of the image reconstruction algorithm, time references indicating the beginning of irradiation and the beginning of each pulse were determined from the list-mode data (see section 2.5.3).

Image reconstruction
Image reconstruction is adapted from the method described in Buitenhuis et al (2017). The main idea is to provide fast reconstruction of the profile of the 12 N activity. Fast feedback and low-count statistics preclude the use of 3D reconstruction based on sophisticated algorithms such as the 3D maximum-likelihood expectation-maximization (MLEM) (Shepp and Vardi 1982); these methods also generally require significant processing time which is not compatible with real-time feedback. A recent study by Muraro et al (2019) comparing the image reconstruction technique used here and the MLEM approach shows no significant difference between the range verification accuracy of the two techniques. Figure 4 illustrates the image reconstruction procedure. One row of crystals in each scanner head is represented by the arcs ABC and DEF, respectively. Coincidence events hitting crystals at locations P 1 and P 2 in opposing heads create a Line of Response (LoR) that intersects the image plane defined by the beam axis and the vertical symmetry axis of the scanner. The 2D histogram of the intersection points for all LoRs then gives a 2D image of the positron annihilation distribution. The 2D image has an area of 208 × 208 mm 2 segmented into pixels of 4 × 4 mm 2 . The accumulated 2D histograms were corrected for the non-uniform scanner sensitivity by dividing the histogram by the 2D sensitivity maps measured as described in section 2.4. The sensitivity correction is important because it makes the shape of the PET activity profile independent of the location of the Bragg peak within the scanner FoV. To account for the slow recovery of the coincidence rate after the PMTs are switched on, the image pixel count for each LoR is divided by the relative coincidence count rate (see section 3.2) associated with the time of the coincidence event.
Note that the coincidence resolving time of the scanner is 550 ps, which corresponds to a spatial resolution due to time-of-flight of just over 8 cm. As this is larger than the transverse width of the positron annihilation distributions being imaged and comparable to the target size, time-of-flight information cannot significantly improve the images and thus is not used. No corrections were necessary for random coincidences and dead-time due to the low activity concentration, and thus count rate, achieved with beam-induced activity in proton therapy, of 0.2-10 kBq/cm 3 /Gy (Parodi 2011). We observed a true-to-random coincidence ratio of 20 in the first 10 s of irradiation. In standard nuclear medicine, a true-to-random coincidence ratio of 5-10 is observed for brain scans, while it may become smaller than 1 when high activities are used (Cherry et al 2012). Similarly, no corrections were made for scatter coincidences and attenuation in the targets. In a clinical implementation, the measured 12 N image will be compared with the image calculated on the basis of the treatment plan. Scatter and attenuation effects will be included in the calculation and not in the image reconstruction of the measured data. It is relevant to note that the results we present would not change by implementing these corrections because the range precision is almost completely determined by counting statistics. Implementing attenuation and scatter corrections does not change the counting statistics.

12 N profiles
To retrieve the 12 N profile, two images, early and late, which are differentiated on the basis of the event time, were reconstructed. The early images were reconstructed using the events from the time period 1-59 ms and contain contributions from both 12 N and longer-lived nuclides. The late images were reconstructed using the events from the period 60-86 ms and contain essentially only counts from longer-lived nuclides. The 12 N contribution was obtained by subtraction of the late images from the early images. Prior to the subtraction, the late image was weighted with a factor 2.19, the ratio of the duration of the early to late time windows.

Determination of pulse time reference
The selection of early and late 12 N events is performed relative to a time reference corresponding to the start of the detector-on periods (see section 2.3). The list-mode data contains no information on such time reference. In the absence of this information, we resorted to an evaluation of the data trend to obtain an accurate determination of the start of the irradiation and the start of the first detector-on period. The starting times of the subsequent detector-on periods were then determined relative to the first one by adding multiples of the 100 ms beam/detector pulsing period. The time spectrum of the first five pulses during an irradiation of graphite with 10 8 protons per pulse is shown in figure 5 (top). The boundary between the low and high count rate region, which corresponds to the start of the first detector-on period after the first beam-on pulse (hereafter called the time reference), was determined by a method which involves taking the time difference between consecutive counts (figure 5, bottom). Large time differences between consecutive counts are seen before the start of the irradiation due to the low count rate. As the irradiation starts, a count rate pattern reflecting the detector pulsing is observed: there are virtually no counts for 13 ms (the sum of the beam-on period and detector-on delay, see section 2.3), followed by 87 ms of relatively high count rate. This generates a regular pattern with a period of 100 ms in the time between consecutive counts. Thus, we set the time at which this pattern appears as the time reference.

Detection of range shifts
The pixel counts of 2D images reconstructed using the method described in section 2.5.1 were integrated along the vertical axis of the image into 1D projection profiles. Using these 1D activity profiles, the activity range R A , defined as the 50% activity fall-off position, was estimated by fitting a Region of Interest (RoI) around the distal edge of the activity profile with a sigmoid function given as where s (x) is the value in bin x of the sigmoid function with amplitude s (max)and growth rate r. Figure 6 shows the procedure for defining the RoI (shaded area in figure 6). The RoI covers all data points distal to and within seven bins proximal to an automatically determined reference point. The choice of seven bins ensures that the fit covers the flat region of the profile. A linear segmentation technique (Killick et al 2012) was used to partition the cumulative sum of the 1D profile in two parts. The position of the partition which maximizes the differences in the mean value of the two data parts was used as the reference point for defining the RoI.

Scanner sensitivity
The absolute system sensitivity at the center of the FoV was measured to be 2.2%. The coincident counts as a function of the bin value (bin width of 4 mm is equal to the scintillator pixel width/height) in the 1D projection of the images of a line source placed in the center of the FoV along either the X-axis or the Y-axis of the scanner are shown in figure 7. Within the measurement accuracy, the sensitivity is the same in both orthogonal directions of the imaging plane. This indicates that the effect of the curvature of the panels is very small. The spike seen between pixels 11 and 15 in the X profile and pixels 30 and 34 in the Y profile represents the points where the line sources were connected and partially overlapped. The average of seven measurements with different positions of the joint was calculated; the spikes were excluded in the data before averaging. The 2D sensitivity map, as shown in figure 7 right, was constructed as the product of the normalized image plane X and Y sensitivity profiles.

12 N nuclide identification
The coincidence count versus time for the first five pulses during an irradiation of a PMMA target with 10 8 protons per pulse is shown in figure 8(a). It can be seen that the time spectrum is dominated by a fast decaying contribution. The periodic absence of counts corresponds to the detector-off periods. Figure 8(b) shows the time spectrum for the sum of counts from the first 50 pulses of the irradiation, representing a total irradiation time of 5 s. During the irradiations, the detectors were switched on with a delay of 3 ms after beam off such that the beam intensity had sufficiently decreased (see section 2.1). As discussed in section 2.3, the PMTs exhibit a lower gain immediately after switching on which then recovers to the nominal value over a period of about 30 ms. As a consequence, a fraction of the actual coincidences between 511 keV photons   Figure 8(d) shows the time spectrum from figure 8(b) after correction for the coincidence rate recovery. A fit of this spectrum in the period 1-86 ms shows a contribution with half-life of 11.2 ± 0.2 ms, a value consistent with the 12 N half-life. During image reconstruction, the image pixel count for each LoR is divided by the relative coincidence count rate associated with the time of the coincidence event. The total number of 12 N counts in the spectrum, for a 3 ms counting delay, following the irradiation of graphite with 10 8 protons in 10 ms is 5.5 × 10 2 counts or 5.5 × 10 -6 counts per proton and 4.3 × 10 2 counts or 4.3 × 10 -6 counts per proton with and without corrections for coincidence rate recovery, respectively; a loss by a factor 1.28 without the correction.

Imaging of 12 N
The images of the positron activity for a single pulse during the irradiation of graphite with a pencil beams are shown in figure 9 left and right columns for 10 8 and 10 9 protons per pulse, respectively. The procedure described in section 2.5.2 was used to separate the 12 N contribution (figure 9 bottom row) from that of the longer lived ones (figure 9 middle row). A broader spread, resulting from the large 12 N positron range, of the activity profile in the Y direction is observed in the 12 N image. 1D projection profiles of the 12 N, long-lived contribution (late images) and both long-and short-lived contributions (early images) for the sum of five pulses with 10 9 protons per pulse are shown in figure 10. A comparison of the 50% activity level indicates that the late profile is shifted by about 12 mm proximal to the 12 N profile. The shifts in the profile may be explained by a convoluted effect of the large positron range of 12 N and, to a lesser extent, the energy dependence of the cross sections for production of the main positron emitters-12 N (Rimmer and Fisher 1968), 11 C and 10 C (Matsushita et al 2016). In addition, the late image profile features an extended distal slope that may be attributed to the energy dependence of 10 C production, in particular to the large energy difference between the threshold (20 MeV) and maximum yield (60 MeV) (Matsushita et al 2016). . 2D reconstructed PET images for irradiation of graphite with 10 8 protons per pulse (left column) and 10 9 protons per pulse (right column). The images shown are for the 10th pulse of the irradiation. The beam direction is indicated by the leftward arrow. The images in the top and middle row were reconstructed using events occurring between detector-on time 1-59 ms and 60-86 ms respectively. The bottom row shows the 12 N image after the scaled subtraction of the middle images from the top images. All images were corrected for the scanner sensitivity (section 3.1) and the coincidence detection efficiency recovery factor (section 3.2).

Range measurement using 12 N
For the range measurement, the 1D 12 N activity profiles during the irradiation of various target configurations were investigated. Figure 11 shows, for the PMMA target configurations described in section 2.2 and various numbers of protons, the 12 N profiles within the first second of irradiation where the contribution of long-lived isotopes is very small for all pulses. Consequently, for evaluation of the performance of 12 N imaging we consider the data acquired after each beam pulse to be equivalent to a single, independent, irradiation with 10 8 or 10 9 protons. Thus, for each of the three target configurations, there are 10 equivalent profiles with 10 8 and 10 9 protons, respectively. The activity ranges and measured range shifts in graphite and PMMA targets for different configurations are presented in table 2. The following steps were taken to obtain the activity range, measured range shift and measurement precision: (a) Bootstrap sampling (Efron 1979) was used to generate 50 sampling sets of 10 samples each. Each sample contains all PET counts measured after the detectors are switched on. The samples are selected with replacement (i.e. the PET count data of an individual pulse can be present more than once in a sampling set) from the measured PET counts of the first 10 pulses of an irradiation. In such bootstrap sampling, many sampling sets are generated using random sampling with replacement from a limited number of measured data (the first 10 pulses in our case) to retrieve robust measures of statistical properties, such as the mean and standard deviation, that are derived from the data. (b) When selecting with replacement the samples in the sampling set from the measured PET counts from the first 10 pulses of irradiations, information is obtained for 10 8 and 10 9 protons per pulse. To obtain count statistics consistent with irradiations with (5, 8, 10) × 10 8 or (5, 8, 10) × 10 9 protons per pulse, the samples in the sampling sets are the sums of five, eight and ten samples selected with replacement from the 10 pulses. (c) The PET images of the data in each sampling set were reconstructed using the method described in section 2.5.1. (d) Each 1D profile was fitted with a sigmoid as described in section 2.5.4, the R A parameter being a measure of the activity range. The change in the R A parameter with respect to the reference case provides a measure of the shift of the proton range due to changes in the target configuration (downstream shifts by 3 mm and 6 mm; range shifting by attachment of 9.8 mm and 4.7 mm PMMA blocks to the proximal surface of targets). (e) The activity range and precision presented in table 2 were determined as the mean and its standard deviation for the 50 bootstrapped sampling sets.
The average precision of the activity ranges from 3 separate irradiations, with the target in the configurations given in table 1, is 2.5 and 2.6 mm with 10 8 protons per pulse and 0.9 and 0.8 mm with 10 9 protons per pulse for PMMA and graphite targets, respectively.
The change in precision of the activity range determination using 12 N profiles as an irradiation progresses, and as function of the number of protons, is shown in figure 12. At each time point considered, data from a 1 s period as specified in the figure, containing 10 proton pulses, is evaluated. The values were determined from the standard deviation of 50 bootstrapped samples as was used for the first second of irradiation. The data points were fitted with a power law function of the form: BN −0.5 + C, where N is the number of protons. The exponent 0.5 represents the effect of Poisson counting statistics and C is a parameter accounting for the increased contribution of the long-lived nuclides and systematic effects. The parameters of the fits are given in table 3. The higher precision as the number of protons increases is clearly seen and conforms quite well to the fitting function. The precision, however, reduces as the irradiation progresses. This is due to the increasing significance of the contribution of the long-lived nuclides in the image subtraction process.

One-dimensional spot scanning
The mixed range target configuration described in section 2.2 was used to investigate the shifts in range across an inhomogeneous interface. Figure 13 shows the 1D 12 N projection profiles following the irradiation of spots across a homogenous target and a mixed range target. The ordering of the spots is indicated in the insets of figure 13. For the mixed range target, the differences in the R A parameters of the sigmoid fit of the 1D profiles of spot 4, chosen as reference, relative to spots 1, 2, 3 and 5, averaged over the similar spots delivered in the first 20 s of the irradiation, are 9.6 ± 1.4, 10.5 ± 0.5, 3.0 ± 1.7, 1.0 ± 1.6 mm respectively. While the range difference of spot 1, 2, 4 and 5 are consistent with the results obtained in section 3.4 for the same range shifter, spot 3 shows a smaller shift which is attributed to range mixing at the interface of the two target thicknesses. Table 2. Precision of range verification for various configurations. The measured shifts are relative to the nominal (0 mm) configuration. The uncertainty in the shifts is propagated from the uncertainties (1σ) on the activity range which was determined from 50 bootstrapped samples of 10 profiles from the 10 measured profiles within the first second of the irradiation. The value shown for 10 × 10 8 is for a combination of 10 profiles from irradiation with 10 8 protons per pulse while that shown for 1 × 10 9 is for a single profile from irradiation with 10 9 protons per pulse. The differences in the 0 mm position in the graphite targets for irradiations with 10 8 and 10 9 protons are due to different positioning of the target.

Discussion
In this study, we investigated the feasibility of near real-time range measurements for individual pencil beam spots on the basis of beam-induced 12 N activity using two modules of a commercial PET scanner for irradiations with clinical beam fluences and intensities. The implementation of a custom switching of the PMTs of the scanner detector enables the elimination of performance degrading effects due to the strong radiation flux during beam delivery and its potentially damaging effects on the PMTs. During the measurements it was observed that the measured coincidence rates for long lived beam-induced activity increased to a stable value over a 25-30 ms period after switching on the PMTs. This phenomenon is due to  recovery of the PMT gain after the PMTs are switched on, as was demonstrated in measurements with a 68 Ge source. A time-dependent correction factor of the measured coincidence rates was applied to account for this effect.
To estimate the gain in accuracy in case of prompt recovery of the coincidence count rate immediately after the beam is switched off, we compare the extrapolated fit to 3 ms before the start of the detector-on period (see figures 8(b)-(d)) for scenarios with and without coincidence recovery corrections. This analysis shows a loss in 12 N counts by a factor of 1.57 without correction for such recovery after the beam pulse and starting acquisition with a 3 ms delay after switching the beam off. The factor derived from this comparison is consistent with the product of gains arising from two contributions. Firstly, a gain by a factor of 1.28 due to coincidence recovery correction. Such correction allows a more accurate subtraction of images from two windows to retrieve the 12 N image. Translation to a gain in precision by the same factor can only be obtained if the correction is implemented on the energy of each event prior to coincidence sorting. Secondly, the removal of the 3 ms delay before the detector is switched on (see section 3.2) in case the pulsing of the beam can be realized much faster will result in, a gain in the 12 N counts by a factor of 1.21.
A shorter beam-on period will also result in a higher precision activity range measurement based on the decay of 12 N. As the PET system is switched off during the beam-on period, the beam-induced decays during this period are not detected. Thus, a relatively short beam-on period ensures that more PET decays are measured during the beam-off period. Reducing the beam-on period from 10 ms, as used in this study, to 3 ms (which is quite typical in patient irradiations), a gain of a factor of 1.27 in the measured PET counts will result. For scenarios with no delay before the detector is switched on (a gain by a factor of 1.21), correction for PMT gain recovery (a gain by a factor of 1.28), and with a reduction of the beam-on period from 10 ms to 3 ms (a gain by a factor of 1.27 due to the reduction of spot delivery period), a gain in the precision by a factor √ (1.21 × 1.28 × 1.27 = 1.40 is expected. This study shows that the activation range can be measured with a precision of 0.9 mm (1σ) for irradiations with 10 9 protons at the 50 ms time-scale. Irradiations with 10 8 protons are about three times less precise. By extrapolating the fit of the precision as a function of number of protons to smaller values, we estimate that a precision of 3.4 mm can be realized for an individual spots of 4 × 10 7 protons. The precision mentioned does not take into account the potential gains with an optimized system, irradiation and prompt data acquisition. Buitenhuis et al (2017) investigated the feasibility of 12 N-based range verification of single spots using a scanner with a peak absolute detection efficiency at the center of the FoV of 0.27%. The authors report a range shift detection precision of 3 mm for the irradiation of graphite with 2.5 × 10 10 protons over an irradiation time of 120 s. Given the dependence of the precision on the number of counts, a 120 s irradiation with a 25 times smaller number of protons, i.e. 10 9 protons as used for a single spot in this work, the experiment of Buitenhuis et al would show a precision of 3 × √ 25 = 15 mm. The precision of range shift detection presented in this study is 1.2 mm (the average of the four entries in table 1 for 10 9 protons per pulse on a graphite target), 12.5 times better than deduced from Buitenhuis et al (2017). This better precision is attributed to the following factors. Firstly, the sensitivity of the scanner used in this study is larger by a factor of about 23 accounting for both the 8 times higher sensitivity in the center of the FoV and an extra factor due to the faster drop of sensitivity towards the edges of the FoV in the scanner used in Buitenhuis et al (2017) relative to the one used here. This gives an improvement in precision of √ 23 = 4.8. Secondly, the smaller contribution of the long-lived nuclides relative to all detected PET counts in the early image (86% for Buitenhuis et al versus 26% in our measurement) gives an improvement in precision of a factor of 3.0. The beam-on period of our measurement is 3 times shorter which, when also taking into account the 3 ms delay between the end of the beam-on period and the start of PET data acquisition, leads to a 1.37 times higher 12 N count. Due to differences in the target dimension of 12 × 12 × 15 cm 3 used here compared to 5 × 5 × 5 cm 3 used in Buitenhuis et al (2017), the number of detected coincidences from 12 N decay is a factor of 2.8 less in our experiment because of attenuation of 511 keV photons in the target. Altogether these factors result in a gain in precision by a factor 4.8 × 3 x √ 1.37/2.8 = 10, satisfactorily close to the factor 12.5 deduced from the experimental results.
A limitation to the use of PET is the small amount of secondary radiation, and thus fewer counts compared to conventional imaging applications. This limitation motivates the development of a high sensitivity detection system allowing the detection of the small beam-induced positron-emitter activity. The requirements for beam delivery nozzle access and unrestricted patient positioning in the treatment field, however, imply that a more sensitive closed-ring geometrical configuration, as in conventional PET systems, is impractical for in-beam monitoring applications without modifications. One such modification to a full ring system is the OpenPET system. The system uses either 2 full-rings with an axial gap (Yoshida et al 2017) or a slanted full ring (Tashima et al 2016) to enable unrestricted access to the patient. In addition to modified full ring systems, a number of partial ring PET systems including those developed at GSI Darmstadt for monitoring of carbon ion therapy (Enghardt et al 2004), at the National Cancer Centre in Kashiwa, Japan (Nishio et al 2010) for monitoring proton beams, and more recently by the INSIDE project at CNAO, Italy (Ferrero et al 2018) have been implemented as solution to the geometrical constraints.
The measurements presented here were performed with 2 panels, each having a dimension of 21 × 21 cm 2 , installed in a partial ring geometry. The scanner was installed with a separation distance of 25 cm between the centers of the opposing heads. The sensitivity at the center of the scanner FoV is 2.2%. In the following, we estimate the influence of scanner geometry and separation distance on the range precision for irradiations with 4 × 10 7 , 10 8, 5 × 10 8 and 10 9 protons. The coincidence counting rate R, measured by a scanner system with intrinsic efficiency ε and a panel separation distance d, exposed to a radioactive point source with activity ξ, is given as where Ω is the solid angle fraction of 4π, Λ represents the net peak fraction of the 511 keV peak, A is the scanner axial width, f is a factor accounting for absorption within the source or between the source and detector and η is the fraction of a full ring covered by the scanner geometry. The influence of scanner geometry on the precision of range detection is given in table 4. The number of protons per spot values 4 × 10 7 and 10 8 are consistent with dose per spot delivered to the distal layer of an exemplary treatment plan shown in Smeets et al (2012). The values indicated for irradiation with 5 × 10 8 and 10 9 protons correspond to aggregating a limited number of spots in an ideal situation where all aggregated spots have the same range.
We assumed an object with the same attenuation as the 12 cm thick PMMA (linear attenuation coefficient for 511 keV photons = 0.11 cm −1 ) used in this experiment, while changes in the size of the scanner are assumed not to alter the intrinsic efficiency and the net peak fraction. Given that the precision scales approximately with R −0.5 , the values of the range precision shown were obtained from a translation of the experimental values obtained with the scanner geometry from this study to potentially realizable geometries. Before the translation, the experimental values were corrected for a potential gain in precision by a factor of √ 1.28 = 1.13 resulting from correction of PMT gain recovery. Increasing the scanner size by a factor of 2 for scanner panels installed 25 cm apart (dΩ = 57% of 4π) as well as reducing the spot delivery duration and prompt switching of the detectors after the beam pulse potentially allows range verification with precision of 2.0 mm (1σ), 1.3 mm (1σ) and 0.5 mm (1σ) for irradiation with 4 × 10 7 10 8 and 10 9 protons at the beginning of an irradiation. For a wider separation distance of 50 cm, sufficient to cover most body regions, a 3 times larger scanner (dΩ = 29% of 4π) is required to give similar range verification precision as obtained with the scanner geometry used in this study.
In a clinical treatment plan, the majority of the spots contain less than 10 8 protons, with the most weighted distal layer spots having a maximum number of protons of about 2 × 10 8 (Zeng et al 2018, Kooy and Grassberger 2015, Smeets et al 2012. In range verification, no matter which method is used, the distal edge of the irradiation field is most important: confidence on the range accuracy of more proximal layers can be obtained if that of the distal energy layer is as expected. When relating our results to clinical practice, we thus focus on distal spots. To obtain a more precise measurement, the analysis of the PET counts can proceed through aggregation of data from multiple spots as used in Xie et al (2017) and Hueso-González et al (2018). One should note that in the general case of a curved distal edge of the dose distribution, the most distal spots are not all at the same proton beam energy. Thus, range determination when summing spots is a trade-off between a higher number of counts (better statistics) and averaging information over a larger irradiated area/volume. A general conclusion on the effect of summing spots on the possibility to assess the correctness of the measured proton range cannot be drawn. The benefit of summing spots depends on the specifics of each treatment plan, such as the contours of the distal edge, the anatomy of the irradiated volume and the order in which the spots are delivered. As an alternative to spot aggregation, the number of protons and the corresponding dose delivered by selected probe spots in the middle of the target (Chen et al 2018) and spots located in the most distal layer (Tian et al 2018) could be increased to enhance the precision of range measurements while ensuring no deterioration of the overall dosimetric quality of the irradiation plan. Because 12 N is produced only on carbon and not on oxygen (Dendooven et al , 2019, an additional consideration with respect to clinical implementation is the achievable precision when irradiating body tissues with a different fraction of carbon than the materials used in this study. The amount of carbon is obviously less in body tissue (an average of 10at% in soft tissue) than in graphite. As the 12 N counts scale with the carbon content, the precision of range determination will scale with the inverse of the square root of the carbon content, similar to the dependence on the number of protons shown in figure 11. So the precision will be a factor √ 33/10 = 1.8 times worse in soft tissue compared to PMMA .

Conclusion
In this work, the precision of range measurement on a spot-by-spot basis through the measurement of the 12 N decay was investigated. Graphite and PMMA targets were irradiated with pulsed pencil beams (10 ms on, 90 ms off) of 10 8 and 10 9 protons per pulse. A dual panel PET scanner (one-sixth of a clinical Siemens Biograph mCT PET scanner), with the scanner heads installed 25 cm apart, was used to image the positron annihilation profile during the beam-off periods. A peak scanner sensitivity in the center of the FoV of 2.2% was realized in this configuration. The scanner features custom-made switching of the PMTs to protect the detectors from damage due to the high photon flux during the beam-on period and to provide faster PMT gain recovery after the beam is switched off. For the irradiation of PMMA with 10 8 protons per pulse, typical for the most weighted distal layer spots in a patient irradiation, we find a range precision of 2.5 mm (1σ) at the very beginning of the irradiation. Due to the very short half-life of 12 N, a range measurement can be performed within 50 ms of spot delivery. Aggregating spots in an ideal situation, i.e. where all spots have the same range, to a total of 10 9 protons gives a range precision of 0.9 mm (1σ). We show that the range precision scales approximately with the square root of the number of protons, or equivalently, the number of 12 N-related coincident counts.
Shorter beam pulses, faster switching off of the beam and an event-by-event time-dependent PMT gain correction have the potential to increase the number of 12 N-related coincidence counts by a factor of 2. Additionally, in an optimum implementation, a more efficient scanner than used in this work is possible. The precision for range detection using a scanner with a solid angle coverage of 54% with spot delivery time much smaller than the 12 N half-life (11 ms) and optimum beam and detector switching is estimated to be 2.0 mm (1σ), 1.3 mm (1σ) and 0.5 mm (1σ) for irradiation with 4 × 10 7 , 10 8 and 10 9 protons.
The results obtained and the potential gains in precision with more optimized systems suggest that the imaging of 12 N enables the retrieval of fast and valuable feedback that may allow a further optimization of treatment plans resulting in less normal tissue being irradiated.