First test of the prompt gamma ray timing method with heterogeneous targets at a clinical proton therapy facility

Accelerated charged particles such as protons have been proposed for cancer treatment since 1946 (Wilson 1946). These particles stop at a well defined depth (the particle range), where the ionization density reaches its maximum (the Bragg-peak) and falls sharply behind it. In theory, the dose can be focused on the tumour and the damage to healthy tissue is minimized compared to photon therapy.

Conversely, ion beam therapy is a more expensive (Goitein and Jermann 2003) and less mature irradiation technique. Is there enough evidence of clinical superiority (Freeman 2012) to justify its high costs? This controversy (Orton and Hendee 2008), chapters 2.11–2.13, arisen in the medical community spins strongly around the intrinsic range uncertainties (Deasy 1994).

In contrast to photon therapy, patient positioning, organ motion, minor density and anatomical changes during or between treatment fractions, stopping power uncertainties, among other factors (Andreo 2009, Grün et al 2013), are prone to yield high dose deviations with respect to the original treatment plan.

In the scientific community, diverse methods have been proposed during the last decades for in vivo particle range assessment (Knopf and Lomax 2013), but are not yet at the stage of widespread application in clinical routine. One approach is the measurement of the spatial distribution of the prompt gamma rays produced in nuclear reactions between the projectile ions and nuclei of the irradiated tissue. The prompt gamma rays are emitted promptly, exit the patient without severe attenuation and their emission distribution is closely correlated with the dose deposition map of the incident ions, as shown experimentally by Min et al (2006). Such correlation is dependent on prompt gamma energy and tissue composition (Moteabbed et al 2011, Gueth et al 2013, Polf et al 2013, Janssen et al 2014).

In the recent years, there is a trend towards systems with less expensive hardware and, probably, faster translation into clinical practice (Testa et al 2014, Assmann et al 2015). In this context, at OncoRay and Helmholtz-Zentrum Dresden-Rossendorf (HZDR), an innovative method for range assessment complementary to the actively and passively collimated prompt gamma ray imaging systems has been proposed (Golnik et al 2014): the prompt gamma ray timing (PGT). Based on the measurable transit time of ions through matter, the detection times of prompt gamma rays encode essential information about the depth-dose profile. Hence, a single detector in a timing spectroscopy setup aims to verify the range at minimum expense and without collimation, just measuring the time between particle entrance in the target and reception of the prompt gamma. As usual, the accelerator radio frequency can be used as time reference.

The proof of principle experiment carried out at a research accelerator (Golnik et al 2014) raises expectations about the future of the PGT method. To evaluate its potentials as well as its limitations, a dedicated experiment at a clinical therapy facility with heterogeneous targets and a pencil beam is mandatory.

The aim of this paper is:

• To test the robustness of PGT at a clinical proton accelerator against background and bunch time stability, and examine the factors limiting its applicability and precision.
• To characterize the bunch time structure of the proton beam and its dependence on the energy.
• To analyse if the PGT method is able to detect range variations due to the presence of heterogeneities.
• To estimate the precision of the retrievable range differences as a function of the collected statistics.
• To assess the mandatory next steps for translation into clinical practice.

Altogether, this work aims at exploring the problems and weak points that may arise when a technique proven at a research accelerator is translated into a real medical scenario, and proposes some strategies to circumvent or solve the encountered limitations. In other words, we do not intend to present a final clinically applicable prototype, but rather investigate and mark the way towards it.

2.1. Detectors

Three monolithic scintillation detectors are used to compare the PGT spectra at different geometries, as well as to check the robustness and compliance of the method. The experimental data are acquired in parallel for all detectors. Their characteristics are summarized in table 1.

Table 1. Monolithic scintillation detectors available in the experiment.

Alias	Material	⌀ $\times$ Length	PMT model	Rise time / ns	Rationale
B1	BaF2	$\left[25:38\right]\times 30$ mm$^{2}$	R2059	$2\pm 1$	Time resolution
B3	BaF2	$48\times 31$ mm$^{2}$	R2059	$2\pm 1$	Time resolution
L0	LaBr3:Ce	${{2}^{\prime\prime}}\times {{2}^{\prime\prime}}$	R2083	$8\pm 1$	Energy resolution

Note: B3 and L0 are cylindrical, whereas B1 is a tapered cone (for optimum time resolution). All PMTs are from Hamamatsu. The rise time refers to the anode signal.

2.2. Electronics

2.3. Accelerators

For the verification of the PGT method, we conduct an experiment at the Westdeutsches Protonentherapiezentrum Essen (WPE), Germany, a clinical proton therapy facility comprising a Cyclone$^{\circledR}$ 230 (C230) isochronous cyclotron of IBA (Louvain-la-Neuve, Belgium). Complementarily, the detector time resolution is characterized at the Electron Linac for beams with high Brilliance and low Emittance (ELBE) facility at HZDR using the same detector set and electronics. The methodology developed for such analysis is published in Hueso-González et al (2014).

2.3.1. Setup at WPE proton accelerator.

The experiment is carried out in the gantry treatment room TR4 with a single incidence angle (horizontal) pencil beam (no scanning) and manual control mode. The patient table is aligned to the gantry room isocenter. A target holder mounted on a linear stage with remote control is placed on the patient table, see figure 1 (bottom).

Zoom In Zoom Out Reset image size

Protons bunched with $\left(105.98\pm 0.09\right)$ MHz are shot onto a target with an energy adjustable between 70 and 230 MeV. Two hollow polymethyl methacrylate (PMMA) half cylinders with an external diameter of 15 cm, an internal of 5 cm, a length of 20 cm and a density of ${{\rho}_{\text{PMMA}}}=\left(1.18\pm 0.01\right)$ g cm$^{-3}$ are available. Cylindrical slices of 5 cm diameter and different thicknesses can be inserted into the hollow cavity left inside the two joined half cylinders. Available slice materials are PMMA, bone-equivalent tissue or air (hollow slice). Custom target configurations can thus be constructed by joining slices heterogeneously. The enclosing half cylinders are used to reproduce the therapeutic irradiation conditions concerning neutron production and moderation in surrounding tissue. The spatial dimensions are measured with 0.2 mm precision.

The bone-equivalent material corresponds to SB3 cortical bone, model 450 from Gammex-RMI (Middleton, USA). The measured density is ${{\rho}_{\text{bone}}}=\left(1.809\pm 0.009\right)$ g cm$^{-3}$. The electron density, effective atomic number and elemental composition are listed in (Helmbrecht 2015, appendix A).

Various detectors are arranged at the desired angle and distance as described in figure 1 (top) to measure the prompt gamma rays exiting from the target. A photograph of the experimental setup is shown in figure 1 (bottom). For a given proton energy, we label a homogeneous target as full target if the protons are stopped completely inside it.

The positioning inside the gantry room is based on the built-in laser room system and the pencil beam spot size is measured with an EBT3 gafchromic film from International Specialty Products (Wayne, New Jersey) placed on the front face of the target (from beam point of view). The position of the target with respect to the gantry room isocenter is described in section 3.4.

3.1. Detector and module settings

3.2. Proton bunch phase stability

3.3. Proton bunch time structure

3.4. Systematic measurement program

3.5. Data acquisition rate

3.6. Data analysis

3.7. Modelling of PGT spectra

4.1. Intrinsic detector time resolution

Table 2 summarizes the energy resolution and coincidence resolving time (CRT) of all detectors as measured with a $^{60}$Co source.

Table 2. Detector energy resolution and CRT for photons of a $^{60}$Co source.

Alias	Energy resolution	CRT versus B1
B1	$\left(9.8\pm 0.2\right)$%	Reference
B3	/	$\left(290\pm 20\right)$ ps
L0	$\left(3.1\pm 0.1\right)$%	$\left(420\pm 20\right)$ ps

Note: The energy resolution (FWHM) is measured at the 1.173 MeV photopeak. The time resolution (FWHM) is measured always with respect to the B1 detector and an individual CFD threshold of 1.0 MeV.

At the ELBE accelerator, the result of the intrinsic detector time resolution as a function of the photon energy is shown in figure 2. The electronic resolution (precision of the time stamp) with which the ELBE RF is measured, $\left(90\pm 10\right)$ ps FWHM, is already subtracted.

Zoom In Zoom Out Reset image size

4.2. Illustrative energy over time spectra

Figure 3 shows three illustrative cases of energy over time spectra for 230 MeV protons. In (a), corresponding to an upstream detector, we can distinguish clear prompt gamma lines corresponding to 6.1 MeV and 4.4 MeV, as well as single and double escape peaks. In addition, we can identify the time uncorrelated gamma line of neutron captures on hydrogen at 2.2 MeV (horizontal line). Conversely, (b) evinces a more compressed PGT distribution, no energy-resolved gamma lines (B3 has worse energy resolution) as well as a considerable increase of background. The bump located at (4 to 6 ns; 1 to 4 MeV) may be related to fast neutrons (secondary gammas), whereas the vertical structure at ($\sim$7 ns; 3 to 8 MeV) could be associated with scattered protons leaving the target. (c) corresponds to a thin PMMA target, and vertical structures due to scattered particles are also visible.

Zoom In Zoom Out Reset image size

The reason for the compressed PGT profile is the angle $\alpha$: for downstream detectors, the protons and gammas fly both downwards, and the gamma time of flight to the detector for the entrance point is larger than for the stopping point. This reduces the spread of arrival times at the detector. On the contrary, for upstream detectors, the gammas travel backwards, in opposite direction to the protons. The gamma time of flight is larger at the stopping point than at the entrance point, and the emission time spectrum is stretched. For $\alpha ={{90}^{\circ}}$ and large d, all gamma rays need the same time of flight to the detector, and the PGT spectrum reflects just the emission time distribution.

4.3. Proton bunch phase stability

The PGT spectra shown in figure 4 (left) correspond to independent measurements with the B3 detector at $\alpha ={{40}^{\circ}}$ for 230 MeV protons impinging a full homogeneous PMMA target. The mean value of the depicted distributions is shown in figure 4 (right).

Zoom In Zoom Out Reset image size

There is an evident drift of the PGT spectrum at a scale of hours: the mean value shifts around 400 ps in a time span of four hours. Albeit a constant target position and proton energy, the PGT spectrum is not stable. The origin of this deviation is attributed to the phase shift of the proton bunch with respect to the RF signal.

The results of the complementary experiment to analyse the stability of the bunch phase with respect to the RF when slightly varying the cyclotron main coil current are shown in figure 5 (left). The dependence of the distribution's mean value on the magnet current is plotted in figure 5 (right) and is more or less linear in this range.

Zoom In Zoom Out Reset image size

A shift in the time distribution of about 14 ps mA$^{-1}$ at an average current of $\sim$740 A is measured. This reflects the extremely high sensitivity of the PGT spectrum mean value to cyclotron instabilities: relative current variations of one per million result in measurable center of gravity shifts and thus affect or limit the long-term precision of the PGT method if no bunch phase monitor is deployed.

4.4. Proton bunch time structure

In figure 6 (left), the dependence of the bunch width on the proton energy is demonstrated indirectly when measuring the prompt gamma rays produced in a thin PMMA target.

Zoom In Zoom Out Reset image size

The presence of a tail right from the main peak in the spectra is related to secondary charged particles, mainly scattered protons, as the thin target does not stop them completely and can thus generate delayed but correlated background (see figure 3(c)). This is not suppressed even though we increase the gamma energy threshold to 3 MeV for figure 6 (statistics are not critical in this case). The bunch spread is measured with different detectors and detailed in table 3 by determining the FWHM of the PGT distributions and subtracting the proton transit time as well as the intrinsic detector and RF time resolution.

Table 3. Proton bunch time spread (FWHM) at the gantry TR4 of the WPE accelerator for different detectors (${{\alpha}_{\text{B}1}}=-{{140}^{\circ}}$, ${{d}_{\text{B}1}}=20.0$ cm, ${{d}_{\text{B}1,100\,\text{MeV}}}=30.0$ cm, ${{\alpha}_{\text{B}3}}={{40}^{\circ}}$, ${{d}_{\text{B}3}}=60.0$ cm, ${{d}_{\text{B}3,160\,\text{MeV}}}=45.0$ cm, ${{\alpha}_{\text{L}0}}={{140}^{\circ}}$, ${{d}_{\text{L}0}}=30.0$ cm), three proton energies and the usual slit closing (25 mm).

Alias	100 MeV	160 MeV	230 MeV
B1	$\left(2.1\pm 0.2\right)$ ns	$\left(1.3\pm 0.1\right)$ ns	$\left(310\pm 30\right)$ ps
B3	$\left(2.1\pm 0.1\right)$ ns	$\left(1.0\pm 0.2\right)$ ns	$\left(380\pm 50\right)$ ps
L0	$\left(2.1\pm 0.1\right)$ ns	$\left(1.2\pm 0.2\right)$ ns	$\left(300\pm 40\right)$ ps

Note: The fit is done with a prompt gamma energy threshold of 3 MeV. The proton transit time, intrinsic detector time and RF resolution are already subtracted.

Figure 6 (right) shows the resulting bunch time spread at 100 MeV for different values of the momentum slit, see section 3.3. All measurements of section 4.5 are done with a momentum slit opening of 25 mm, i.e. the usual clinical setting.

4.5. Systematic measurement program

4.5.1. 230 MeV proton energy

Stacked target. Figure 7 (left) evinces the potential of the PGT method at an energy of 230 MeV compared to 160 or 100 MeV due to the very low bunch time spread, see figure 6 (left). Individual characteristics can be identified along the spectrum shape, instead of focusing on the distribution momenta, as done in the experiment at the AGOR accelerator (Golnik et al 2014).

Zoom In Zoom Out Reset image size

In the case of the full PMMA target, the proton stopping point is at 28.0 cm depth. A smooth leading edge is identified in the PGT spectrum, which is associated to the front face of the target (beam entrance point). A trailing edge related to the distal fall-off is found, but with lower relative intensity. Despite the expected increased gamma production rate within the Bragg peak, the measured profile does not reproduce it directly due to the following factors:

• Solid angle effects for upstream detectors, i.e. Bragg peak is much farther away than the entrance point.
• Larger path through PMMA (attenuation) for the prompt gammas emitted at the distal edge.
• Transformation between measured time, gamma emission time and nuclear interaction depth. As the particles' velocity is not constant, the conversion from space to time is not linear and the histogram is stretched over more time bins close to the Bragg peak (slowest protons) than at the beam entrance point (fastest).

The removal of a slice from the target's rear face links to a decrease in the area and shift to the left in the mean value of the PGT distribution over the whole thickness range studied (except for the 400 mm case). For the thinner targets (thickness $\leqslant$50 mm), the reduction in the area correlates with a decrease in the maximum of the curve, as the proton transit time effect has the same order of magnitude as the system time resolution. Conversely, for the thicker targets, the effect of removing a slice can be clearly identified as a decrease of the prompt gamma production at a given depth, which changes the right leaf of the distribution but not the maximum. Comparing the 300 mm and 400 mm cases, we observe no change at all. This confirms the intuitive idea that adding material beyond the particle range (280 mm) does not have any effect on the PGT profiles, as the protons are completely stopped at equal depth in both cases and the region of prompt gamma production is equivalent.

Figure 7 (right) shows the PGT distributions based on the simBox model. The general shape along the depth axis is reproduced, however, near to the right fall-off, the model fails to incorporate the rectangular edge of the 300 and 400 mm targets, in which the protons stop completely. This is related to the fact that the simBox model underestimates the gamma production rate close to the Bragg peak. Furthermore, the background is not included in the model.

Air cavities. An extensive study about air cavities of different thicknesses and locations inside the full PMMA target is accomplished for 230 MeV protons, as shown in figure 8. By comparing the PGT spectra, one can visually identify a deficit in the gamma production of a magnitude correlated to the cavity thickness at the time position corresponding to its location, as well as a shift to the right in the position of the trailing edge (overshoot). Air cavities and thus range differences down to 2 mm can be retrieved with the gathered statistics and without any sophisticated numerical analysis.

Zoom In Zoom Out Reset image size

These experimental spectra are in significant agreement with the modelled ones, if we exclude the region near the Bragg peak. The centroid of the cavity, namely the time point with the largest deficit in counts with respect to the homogeneous case, is depicted in figure 9 for different locations z of the 10 mm cavity. The values expected by the simBox model are in reasonable agreement with the experimental ones.

Zoom In Zoom Out Reset image size

Bone inserts. A bone-equivalent slice of 20 mm thickness is inserted at different locations inside a PMMA target for 230 MeV protons. The PGT distributions measured with the L0 detector are shown in figure 10. The bone insert correlates evidently with an increase of the prompt gamma production with respect to the homogeneous case at the bone location. Moreover, a shift to the left in the falling edge (undershoot) is visible. These effects are basically attributed to the higher density of bone with respect to PMMA. Hence, about 8 mm range shifts due to a bone heterogeneity (according to stopping power data) can be at least retrieved with the collected statistics. Note that the bone at 269 mm depth is a special case, as the proton stopping point is inside the bone insert.

Zoom In Zoom Out Reset image size

Figure 10 (right) shows a tentative conversion of PGT spectra to depth profiles, after corrections for solid angle effects, attenuation, gamma time of flight and transformation of time to space with the stopping power differential equation as in (Golnik et al 2014, equation (5)). In other words, a simple backprojection is applied to recover the depth profile information from the experimental data. Note that we apply an iterative background subtraction algorithm (Ryan et al 1988). However, at this stage, the space origin is elected arbitrarily and we do not apply a deconvolution with the system time response, which explains the absence of a sharp leading edge at $z=0$.

For all bone insert measurements (the special case at 269 mm aside), the trailing edge is $\left(7\pm 2\right)$ mm left from the homogeneous phantom, while the expected range difference is $\sim$8 mm. It is worth pointing out that the a priori assumption of the target geometry, stopping power and time to space conversion curve predetermines (biases) the respective position of the falling edge in figure 10 (right).

The centroid of the bone, namely the space point with the largest surplus in counts with respect to the homogeneous case, is depicted in figure 11 for all locations z of the 20 mm insert. The values expected by the simBox model are in considerable agreement with the experimental ones.

Zoom In Zoom Out Reset image size

4.5.2. 100 MeV proton energy

Stacked target. For the stacked target experiment at the WPE accelerator (figure 12), analogous to that at the AGOR facility (Golnik et al 2014), one can identify a shift to the right in the peak centroid correlated to the target thickness, as well as an increase in the prompt gamma production (area under the curve). Compared to figure 7, the PGT profile is highly blurred due to the much larger bunch time spread at 100 MeV (see table 3). Nonetheless, we can rely on the analysis of distribution momenta, which retain essential information about the proton range.

Zoom In Zoom Out Reset image size

Figure 12 (right) evinces the linear relationship between centroid and target thickness until a saturation or plateau is reached, i.e. when the protons are already completely stopped (proton range is 6.5 cm). As in Golnik et al (2014), a reasonable agreement between modelling and experiments is found: when fitting the slope of the left part of the graph, we obtain $\left(5.1\pm 0.1\right)$ ps mm$^{-1}$ for the experiment versus $\left(5.03\pm 0.03\right)$ ps mm$^{-1}$ for the simBox model. Note that the election of the global time offset in the experiment relative to the model is arbitrary. We select the value maximizing the agreement with all the experimental points.

Air cavities. The influence of tissue heterogeneities is illustrated in figure 13 (left): air cavities of different thicknesses are inserted at a fixed position inside the full PMMA target. On the one hand, the spectrum maximum decreases with increasing thickness of the cavity. On the other hand, the falling edge of the spectrum shifts slightly to the right as a consequence of the overshoot. Thus, at least 5 mm range differences can be distinguished at 100 MeV. The shift to the right in the centroid of the distribution for the 5 mm cavity is $\left(28\pm 1\right)$ ps for the experiment and $\left(27.1\pm 0.2\right)$ ps for the simBox model, while for the 10 mm cavity, the experimental shift with respect to the homogeneous case is $\left(65\pm 2\right)$ ps versus $\left(53.6\pm 0.2\right)$ ps for the model.

Zoom In Zoom Out Reset image size

Bone inserts. Figure 13 (right) compares a full target PGT spectrum to the case with a 20 mm bone insert at 29 mm depth. The distributions are normalized to their respective maximum, as the higher amount of gammas produced in the heterogeneous case shadows the comparison of the trailing edge of both spectra. We can conclude that a bone insert increases the gamma production rate as well as decreases the width of the time distribution due to the undershoot effect ($\sim$8.5 mm range difference according to stopping power data), which results in an earlier average detection time of the prompt gammas. The shift to the left in the distribution's centroid with respect to the homogeneous case is $\left(28.6\pm 0.6\right)$ ps for the experiment and $\left(17.9\pm 0.4\right)$ ps for the simBox model.

Motivation. A main purpose of the present experiment is the estimation of the minimum range difference that generates a measurable change in the shape of the PGT spectrum, rather than the quantitative retrieval of the range. The setup is optimized for timing but is not intended for high rate measurements or long-term detector gain stability. This is considered as a further step once the suitability of PGT in clinical environments is shown at all.

Systematic program. Proton range variations due to a stacked PMMA target as well as cavity and bone heterogeneities are detected. In most cases, an energy threshold of 1 MeV is applied to the data for maximizing the statistics, despite the fact that the correlation of the gamma emission profile with the dose is higher for an energy threshold of 3 MeV (Verburg et al 2013). For 100 MeV proton energy, range differences of 5 mm can be successfully retrieved (without any sophisticated numerical analysis). These are identified qualitatively based on the distinctness of the spectra compared to the homogeneous case. The quantitative assessment relies on the comparison of distribution momenta between experiment and simBox model, with overall reasonable agreement. There are some systematic deviations in certain cases, either due to bunch-to-RF phase drift or because of the oversimplified model.

For an energy of 230 MeV, the much lower bunch time spread allows to recover information across the shape of the PGT spectra and not just through distribution momenta. This could be deployed for depth-emission profile reconstruction with a single detector, a milestone in the field of in vivo proton dose verification. Figures 9 and 11 demonstrate the potential for differentially correlating time to depth, and imply a qualitative leap forward of the PGT method with respect to Golnik et al (2014), notwithstanding the need of more advanced conversion tools.

First tested at the AGOR research accelerator, the robustness of the PGT method is examined in terms of background, bunch time spread and bunch stability during the experiments in clinical facilities. In this scenario, we encounter several factors limiting the accuracy of PGT in its simplest form. First attempts to understand and quantify these effects are described, and strategies to circumvent them are proposed.

Bunch structure. The bunch spread is measured indirectly with a thin PMMA target. Together with the prompt gammas, there is a high background of scattered protons. This motivates the use of dedicated charged particle detectors in upcoming experiments to determine the bunch width (and phase) directly from protons, e.g. those scattered in the beam exit window.

The increase of the bunch time spread at lower energies is an intrinsic hurdle for the PGT method at the widespread C230 (and, in general, at other clinical fixed-energy cyclotrons relying on a degrader to brake the ions), compared to the AGOR research accelerator. As a consequence, PGT spectra revealing valuable profile information at 230 MeV due to a low bunch time spread may blur completely at 160 or 100 MeV, which are proton energies much more frequent in an actual therapeutic irradiation. Strategies like the reduction of the momentum slit opening suggest that there is still room for improving the system time resolution, or even for optimizing the time spread in future designs of clinical isochronous cyclotrons.

Some research cyclotrons (like AGOR) do not require an energy selection system, as they vary magnetic field and radio frequency to provide just the desired beam energy, but are not usual in the clinics. Rather, fixed-energy isochronous cyclotrons are a widespread choice in hospitals. Synchrotrons or synchrocyclotrons are also available (or planned) in a few treatment centers. A priori, the application of the PGT method at these synchronous accelerators, where the micro bunch time spread may be large and the instantaneous current higher, is not straightforward and has to be explored.

Phase drift. The main burden that affects the experimental data presented in this paper is the bunch phase drift with respect to the RF signal, which prevents us from making further quantitative analysis on centroid shifts. To overcome this instability experimentally, we are currently developing a bunch phase monitor based on phoswich detectors (Wilkinson 1952), which could identify and select by energy the protons scattered at the beam exit window, and measure their crossing time. Then, the actual shift between bunch and RF can be subtracted from the PGT spectra. Additionally, the measurement of the temporal shape of the beam bunch (non-Gaussian shape) is essential for the modelling, in order to apply a convolution with the proper system response.

Within a clinical pencil beam scanning session (few seconds), the effect of long-term bunch drift is expected to be negligible when comparing consecutive spots. Rather, one has to correct the mean drift (and variation of bunch time spread) when comparing the experimental results of two different treatment fractions (one day separation). Moreover, for the comparison of PGT spectra with modelled ones, one has to discriminate the different energy layers within the fraction, as the crossing time of the beam exit window with respect to the RF signal depends on the beam energy (and beam line length). This has to be calibrated accurately beforehand, or monitored with the aforementioned phoswich detector for each energy layer in real-time.

Normalization. The uncertainty of the normalization factor to number of protons is also a major hurdle to the presented data, but is a circumstantial limitation of the dead time model. We expect this provisional problem to be solved with upgraded electronics.

Background. Significant background radiation due to scattered charged particles, neutrons and material activation is observed in the experiments, specially for high proton energies. The PGT method has proven robust for retrieving range shifts in most cases, but the pedestal present in the time histograms is still an important concern. A low signal to noise ratio is observed for the thin target measurements. For thicker targets, the background has a particular (not flat) time structure. Efforts have to be undertaken to either reduce physically the background, or to model and subtract it for minimizing the bias on the retrieved range shift.

Detector geometry. For detectors placed downstream, suppression of charged particles, e.g. by pulse shape discrimination, is advisable to reduce the high relative background generated by scattered protons (as observed in the energy over time spectra). Nevertheless, an upstream positioning provides in principle a better performance, as the gamma time of flight to the detector extends the PGT spectrum and thus allows to resolve more details. In addition, the inherent fast neutron background is lower in backward direction (Gunzert-Marx et al 2004). The drawback is that a high percentage of the measured gammas stem from the target's front face rather than from the Bragg peak position due to a different solid angle coverage. In cases where the accumulated statistics is low, this may be a limitation if a reconstruction of the particle range based on the trailing edge of the PGT spectra is attempted.

Throughput. The measurements described in this paper are made with a pencil beam current between 10 and 100 pA and a duration of 5 to 7 min. Let us give some numbers for a concrete case: the irradiation of a full PMMA target at 160 MeV proton energy. The acquisition lasts 6.5 min and the delivered charge is 23 nC, corresponding to $1.4\times {{10}^{11}}$ protons and a beam current of 59 pA. Assuming an isotropic prompt gamma yield of 0.16 as in Golnik et al (2014), the emission rate of prompt gammas with energies higher than 1 MeV is $59\times {{10}^{6}}\,\gamma$ s$^{-1}$. The solid angle covered by the L0 detector at 30 cm distance from the ring center is $\text{d} \Omega /\Omega =1.8\times {{10}^{-3}}$, the minimum detection efficiency is ${{\epsilon}_{\text{D}}}\sim 0.58$ (Berger et al 2010) and the gamma attenuation inside the target is of the order of 45% (Hubbell and Seltzer 2004). The resulting overall efficiency ${{\epsilon}_{\text{p}\gamma}}$ is $\sim {{10}^{-4}}$ and the measured detector trigger rate is 33 kcps. In our particular acquisition system, as the dynamic dead time is 34%, the acquisition rate is 22 kcps and the total number of registered events is $8.5\times {{10}^{6}}$. Thus, the dynamic efficiency is $\epsilon _{\text{p}\gamma}^{*}\sim 7\times {{10}^{-5}}$. Note that we dismiss background contributions for this calculation.

Besides the bunch time spread, detector throughput is the critical parameter for PGT and has to be addressed for obtaining a reliable range measurement in real-time and in vivo. The electronics required to cope with high count rates up to 1 Mcps without detector overload are under development.

Statistics. What would be the differences in a clinical treatment plan compared to this experiment? Concerning a realistic 2 Gy irradiation fraction with $\sim {{10}^{10}}$ protons (in total) delivered to a small area of the brain with pencil beam scanning at an IBA C230 machine, the instantaneous current at nozzle exit during spot delivery is around 2 nA, i.e. $\sim {{10}^{10}}$ p s$^{-1}$ (Perali et al 2014). A typical pencil beam spot is thus delivered within a few milliseconds. Note that the overall treatment fraction duration includes the times necessary for magnet sweeping between the spots and for energy switching between the layers, that are dependent on the actual treatment plan. Therefore, we base our rate estimates on proton beams of 2 nA peak current measured at nozzle exit.

If we consider the L0 detector setup described above (${{\epsilon}_{\text{p}\gamma}}\sim {{10}^{-4}}$) without dead time effects, the resulting acquisition rate is $\lesssim$1 Mcps (no overload expected with upgraded electronics), i.e 500 kcps nA$^{-1}$. The planned number of protons per spot is $\leqslant 2\times {{10}^{7}}$ for $\sim$75% of the spots. In other words, $\lesssim$2000 gammas are measured for each raster point, so that the analysis of the PGT spectrum seems not realistic on a single spot basis due to insufficient statistics. Rather, one should collect data from several detectors or accumulate gamma rays from neighbour spots in a single spectrum.

How do we quantify this qualitative statement? The statistics necessary for detecting a range shift with 2 mm precision in a PGT spectrum depending on beam current and proton range has been analysed in Golnik et al (2014) based on the theoretical error of distribution momenta. Equivalent conclusions could apply here for 100 and 160 MeV proton energies. In the case of low bunch spread (230 MeV protons), one may want to retrieve a range difference based on the shift of the falling edge (see figure 8), which corresponds to the Bragg peak region, instead of on the mean value or distribution's width. What is the minimum number of gammas in the PGT spectrum for detecting a range shift in this case? There is no general answer for this question, as the detector orientation, target attenuation and amount of background may play a significant role on the shape of the falling edge.

Alternatively, we perform a statistical analysis on experimental data of the L0 detector for 230 MeV protons and three particular target configurations: (a) a full homogeneous PMMA phantom, (b) a 5 mm air cavity inside it, and (c) a 20 mm bone insert in place of the cavity.

The original histograms (complete measurements) have around $7\times {{10}^{6}}$ entries each. Once the region of interest (the fall-off) inside the three PGT spectra is defined, the trailing edge deviation is calculated by finding the virtual time shift on (a) that minimizes the weighted χ² test between the (a)–(b) histograms for the cavity or (a)–(c) for the bone insert. These χ² tests are performed for four different statistical sample sizes (random subsets of the whole measurement): ${{10}^{6}}$, ${{10}^{5}}$, ${{10}^{4}}$ and ${{10}^{3}}$ histogram entries (prompt gammas plus background). Assuming for simplicity the aforementioned efficiency ${{\epsilon}_{\text{p}\gamma}}\sim {{10}^{-4}}$, these virtual measurements would correspond approximately to an irradiation with ${{10}^{10}}$, ${{10}^{9}}$, ${{10}^{8}}$ and ${{10}^{7}}$ protons, respectively. For each case, 100 subsets of equal sample size are selected randomly in order to determine the average time shift and the standard deviation.

Figure 14 shows the complete statistical study. One can see how the reduction of sample size (number of protons) degrades the shape of the PGT spectra. Statistical noise hampers the identification by eye of the shift of the falling edge the lower the number of protons. On the basis of the χ² numerical analysis (bottom right plot of figure 14), ${{10}^{8}}$ protons seem to be the threshold for determining a 5 mm shift of the falling edge in a single spot. Below this boundary, namely the region enclosing most of the clinical spots, the statistical deviation is too high, as illustrated by the overlapping of the coloured uncertainty areas. Looking at the error bars, a range shift down to 2 mm can be identified with ${{10}^{10}}$ protons. It is worth to remember that these considerations apply for a particularly low bunch spread scenario (230 MeV), a single detector in upstream orientation and selected target configurations.

Zoom In Zoom Out Reset image size

How could we increase the collected statistics? One may enlarge the detector volume at the price of developing electronics that sustain a higher throughput. One could also sum up data from several detectors or from neighbour spots to recover statistics, or consider a downstream detector, that maximizes the solid angle at the Bragg peak. More sophisticated numerical analysis including the whole histogram for the shift retrieval instead of just a region of interest is also advisable. For example, in our χ² tests, we waste useful information like the difference in the gamma production at the location of the heterogeneity, which is still visible in the bottom left plot of figure 14. There is undoubtedly room for improvement.

An opposite approach, which would affect clinical workflow and planning, is to deliver a benchmark spot at the distal region, namely the raster point of the plan with the highest number of protons. The rationale is to detect with confidence (along the selected incidence direction) potential range shifts due to anatomy changes between treatment fractions. The number of protons of the richest spot depends on each particular tumour location and treatment planning software, but it is not uncommon to have at least one spot with $\geqslant {{10}^{8}}$ protons, see (Smeets et al 2012), figure 19. One could deliver this benchmark spot immediately at the beginning of the treatment and check the acquired PGT spectrum. If no deviations are observed, one would consecutively repeat this step for all spots of the distal layer to cover other irradiation areas. If the results comply with the expectations, one could then deliver all remaining spots of the plan. Hence, just the spot order delivery but not the dosimetry is affected.

Another option would be to double the detector size and halve the instantaneous clinical current (for all spots) to keep the count rate at an acceptable level, but multiply the detection efficiency ${{\epsilon}_{\text{p}\gamma}}$ by a factor of two. For the particular treatment plan considered above (${{10}^{10}}$ protons at ${{10}^{10}}$ p s$^{-1}$), the actual irradiation time of the small area of the brain lasts about one second. This irradiation time is relatively small compared to the global delivery time ($\sim$15 s), with a significant dead time due to magnet sweeping and change of energy layers. Hence, halving the beam current would not increase significantly the treatment duration in the WPE treatment facility.

A further strategy, which combines the two last, is to deliver at treatment beginning the spot with the highest number of protons of the plan, e.g. ${{10}^{8}}$ protons, with a very low beam current (one proton per bunch, around 1 s irradiation at 2 pA) and measure the passage of every proton with a beam hodoscope. The added value is that the effect of bunch time spread is suppressed for any beam energy, as each detected prompt gamma can be correlated to the primary proton. The collected statistics can also be maximized as the lower beam current allows the use of larger detectors.

Quantification. For the quantitative retrieval of the range, either directly or by comparing with the expected PGT profile, an extensive simulation and modelling based on section 3.7 is under progress. More realistic PGT modelling should consider effects like energy-resolved emission profiles, detector sensitivity for complex geometries, background subtraction and non-Gaussian temporal system response.

Stopping power. The PGT spectra at 230 MeV (low bunch spread scenario) do not provide a proportional image of the gamma emission but rather a stretched one. An inverse transformation is needed as intermediate step if one wants to obtain a depth profile, see figure 10. However, this calibration depends on the accurate knowledge of the stopping power of the traversed tissue (Golnik et al 2014). One may then argue that no method for detecting range uncertainties should depend on their most common source: the variation of the tissue composition (stopping power) between clinical fractions. The absence of a daily patient imaging prevents a recalibration of the time to space relationship for each fraction. Does this mean that the PGT approach is useless?

Let us put it the other way round and calculate a PGT spectrum according to the treatment plan expectations (and uncertainties). If, unluckily, the assumed anatomy is inaccurate, the reference PGT profile will be biased. As a consequence, an overshoot of e.g. 5 mm may occur and the measured PGT profile will change with respect to the expected one due to its dependence on the spatial emission profile, but also on the stopping power itself (simultaneous dependency). Due to the unawareness on the exact anatomy and the bias introduced on the reference spectrum, one may not be able to quantify with utmost accuracy the actual range deviation. But we believe that the differences (in shape or momenta) between the reference and the measured spectra could be blamed (with enough confidence) on e.g. a shift between 3 mm up to 7 mm. In other words, the PGT method could serve as real-time safety assurance for detecting severe deviations rather than accurate range measurements. This strategy is also applicable for lower proton energies (larger bunch spread).

In any case, even for an ideal prompt gamma camera that reproduces the spatial emission profile faithfully (without distortions) and with perfect counting statistics, the amplitude of the actual range shift is biased by the treatment plan expectations. The reason is that the correlation of prompt gamma yield to particle range depends on the assumed tissue composition. Still, experimental measurements that reveal significant differences with respect to biased reference spectra are a reliable indicator that something is not going as planned.

Outlook towards the clinics. At this stage, the experimental results do not allow a comparison of the PGT method with other more developed techniques as e.g. the knife-edge slit camera, which is able to image millimetre shifts at clinical current intensities (Perali et al 2014). Our method is far beneath in the ladder of clinical translation, and one can think of many limiting scenarios in which its applicability is questionable or yet unproven. However, some trends and qualitative aspects are already observed in this work, which encourage the hope of a quantitative retrieval of the range based on the PGT method. Whereas we can not claim a spatial precision comparable to other methods on a single spot basis, we believe that the much lower expense and footprint of this technology may boost the applicability and spread of this technique in the clinics, and thus deserves further investigations.

Thanks to its real-time range assessment capability, we consider that PGT has the potential to be translated into clinical practice as a fast handle to stop a running treatment with significant deviations during the clinical fraction. The workflow for detecting severe range shifts could look like:

• A realistic spatial emission profile is simulated according to the treatment plan.
• The expected timing profile is calculated based on the detector model and patient data.
• PGT spectra are measured during the clinical fraction and corrected for bunch drifts (with respect to the previous treatment session) by a phase monitor.
• The center of gravity of the modelled and measured PGT spectra is compared in real-time.
• If the difference exceeds a clinically critical value, the treatment can be interrupted and the source of error analysed or corrected in the upcoming treatment fractions.

In vivo range verification in real-time is a must for improving the outcome of ion beam therapy. The prompt gamma rays produced in nuclear interactions are reliable signatures for reconstructing the depth-dose profile, provided that the tissue composition is well known. The prompt gamma ray timing method aims at assessing the range based on a simple detector setup without mechanical collimation.

This paper describes the first tests of PGT at a clinical proton accelerator with different phantoms and detectors. The bunch time spread is characterized for different energies and the robustness of the method concerning background and bunch stability is investigated. For ${{10}^{10}}$ protons (230 MeV energy) and a single detector, particle range differences of 2 mm are identified in predefined heterogeneous targets based on the distinctness of timing distributions. For ${{10}^{8}}$ irradiated protons, a 5 mm range shift is retrievable by numerical analysis. The experimental results also prove that PGT spectra comprise very detailed information on gamma emission along the track when the bunch spread is low.

The PGT method appears to be feasible in a clinical scenario with a pencil beam if following intermediate steps concerning hardware, software and procedure are addressed to counteract the encountered limitations:

• Addition of a bunch phase monitor which measures the long-term drift with respect to the accelerator RF as well as the bunch temporal shape for each treatment fraction. The monitor can also correlate the crossing time of the protons through the beam exit window with respect to the RF for each energy layer.
• A much higher detector and data acquisition throughput to cope with the gamma rate expected in an actual clinical treatment fraction.
• An accurate time spectrum and background modelling, and numerical analysis or even iterative algorithms for obtaining a quantitative measurement of the particle range.

Experiments with improved detectors and electronics, realistic currents, beam scanning and anthropomorphic phantoms are the next steps towards clinical translation.

In conclusion, the overall experimental results underline the potential of the prompt gamma ray timing method for range assessment in a clinical pencil beam and reassure this novel approach as a promising alternative in the field of in vivo dosimetry.

This work is supported by the German Federal Ministry of Education and Research (BMBF-03Z1NN12) and by the European Commission within the Seventh Framework Program through ENTERVISION (grant agreement number 264552). We thank M Berthel, A Dreyer, A Hartmann, K Heidel, A Junghans, M Kempe, T Kormoll, A Laso-García, M Sobiella and D Weinberger for the excellent assistance, and the staff of the ELBE and WPE facilities for their support.