Detection of range shifts in proton beam therapy using the J-PET scanner: a patient simulation study

Objective. The Jagiellonian positron emission tomography (J-PET) technology, based on plastic scintillators, has been proposed as a cost effective tool for detecting range deviations during proton therapy. This study investigates the feasibility of using J-PET for range monitoring by means of a detailed Monte Carlo simulation study of 95 patients who underwent proton therapy at the Cyclotron Centre Bronowice (CCB) in Krakow, Poland. Approach. Discrepancies between prescribed and delivered treatments were artificially introduced in the simulations by means of shifts in patient positioning and in the Hounsfield unit to the relative proton stopping power calibration curve. A dual-layer, cylindrical J-PET geometry was simulated in an in-room monitoring scenario and a triple-layer, dual-head geometry in an in-beam protocol. The distribution of range shifts in reconstructed PET activity was visualized in the beam’s eye view. Linear prediction models were constructed from all patients in the cohort, using the mean shift in reconstructed PET activity as a predictor of the mean proton range deviation. Main results. Maps of deviations in the range of reconstructed PET distributions showed agreement with those of deviations in dose range in most patients. The linear prediction model showed a good fit, with coefficient of determination r 2 = 0.84 (in-room) and 0.75 (in-beam). Residual standard error was below 1 mm: 0.33 mm (in-room) and 0.23 mm (in-beam). Significance. The precision of the proposed prediction models shows the sensitivity of the proposed J-PET scanners to shifts in proton range for a wide range of clinical treatment plans. Furthermore, it motivates the use of such models as a tool for predicting proton range deviations and opens up new prospects for investigations into the use of intra-treatment PET images for predicting clinical metrics that aid in the assessment of the quality of delivered treatment.


Introduction
Proton beam radiotherapy is increasingly being considered as a treatment option for oncological patients worldwide. The dosimetric advantage of irradiating a tumor with protons, rather than the more widely used high-energy photons, lies in the characteristic depth-dose profile of protons interacting with matter, the socalled Bragg peak. This depth-dose profile is characterized by a pronounced peak in a localized region at the end of the proton range, which enables excellent dose conformation to the tumor while sparing healthy tissues in the vicinity of the target volume (Durante and Löffler 2010). However, the full translation of the dosimetric properties of protons into clinical benefit is potentially limited by proton range uncertainty. The main sources of the uncertainty are inaccuracies in the estimation of the relative proton stopping power (RPSP) of different human tissues, which is calculated from a planning x-ray CT, as well as even small inter-and intra-fractional variations in patient position and anatomical changes. A reliable method of in vivo verification of the delivered dose, during or shortly after the irradiation, would be greatly beneficial for greater clinical exploitation of the potential of protons for achieving high dose conformation.
Since the primary protons stop at the end of their range, in vivo verification techniques rely on secondary signals correlated with the dose distribution. Most such methods involve the imaging of secondary radiation, resulting from nuclear reactions between the incident protons and atomic nuclei of patient tissues. For an overview of these nuclear techniques, the reader is referred to several reviews (Knopf and Lomax 2013, Krimmer et al 2018, Parodi and Polf 2018, Parodi 2020, Tashima and Yamaya 2022. The most clinically mature imaging method used for ion beam range monitoring is positron emission tomography (PET). For protons, this technique involves the detection of 511 keV annihilation photons following the decay of β + -emitting target fragments from the nuclear interactions between the irradiated tissue and the proton beam. The PET signal resulting from the beam delivery may be acquired during or immediately following the irradiation, using a PET scanner mounted in situ in the treatment position (Enghardt et al 2004, Nishio et al 2010, or following some delay, using a scanner located in or near the treatment room (Hishikawa et al 2002, Parodi et al 2007a, Zhu et al 2011. For the latter approach it is possible to make use of a conventional PET scanner, while having the scanner in situ would require a dual-head (Enghardt et al 2004, Ferrero et al 2018b or otherwise modified PET geometry (Crespo et al 2006, Tashima et al 2016. These unconventional geometries present a challenge during image reconstruction, while the radiation background produced when the beam is on places a high demand on the beam-on data acquisition (Crespo et al 2005, Sportelli et al 2014 as well as the radiation hardness of in situ detectors . However, having the PET scanner mounted in the beam position offers several advantages. Scanning immediately after the beam delivery allows for the acquisition of a larger signal originating from 15 O (the most abundantly produced isotope, with half-life of 2.04 min) and minimizes the effect of biological washout, whereby the activity is displaced by physiological processes away from where it was originally produced (Bennett et al 1978). Furthermore, acquiring data when the beam is on raises the possibility of almost real-time monitoring of the proton range (Ozoemelam et al 2020).
The distribution of PET activity created by the proton beam is correlated to the delivered dose, but it is not identical to it. Dose deposition depends on the inelastic interactions of the proton beam with atomic electrons in the target material, while β + emitters are produced by nuclear interactions between the proton beam and atomic nuclei. Therefore, the abundance and type of β + emitters produced are dependent on the elemental composition of the tissue in a way that dose deposition is not. The varied elemental composition of human tissue and the wide range of possible nuclear reactions of the incoming protons on each element, each with very different crosssections as a function of proton energy, weakens the correlation between β + activity and dose. Furthermore, the finite threshold energy for reactions that produce β + emitters results in a PET activity whose distal edge does not coincide with that of the dose distribution. The implication of these factors is that the reconstructed PET images cannot be compared to the dose distribution from the treatment plan. Instead, the typical approach taken when using PET for proton range monitoring is to use shifts in PET activity distribution along the direction of the beam path as a surrogate for dose deviations, by comparing the measured PET distribution with a reference image.
One possibility is to use a PET image acquired during a previous fraction of the treatment as a reference (Nishio et al 2010). This approach allows the identification of day-to-day variations in dose delivery due to deviations in patient positioning and anatomy, as well as the sudden filling or emptying of organs and cavities. However, deviations from the treatment plan that do not change from day to day, such as errors in the calibration of Hounsfield unit (HU) to RPSP, are not detected using this approach. A PET reference image that more accurately represents the treatment plan may be calculated analytically (Miyatake et al 2011, Frey et al 2014 or using Monte Carlo simulations (Enghardt et al 2004, Parodi et al 2007a, Kraan et al 2014 by including models for proton transport, nuclear production cross-sections, treatment and imaging time structures and the PET detector response. Methods of measuring range shift in PET distribution include visual inspection as well as several proposed quantitative techniques, consisting of comparing the position of the distal fall-off of PET activity along the beam direction (Knopf et al 2009, Helmbrecht et al 2012, Kuess et al 2013, Frey et al 2014, Fiorina et al 2021. Because the PET activity distribution precedes that of the dose along the beam, the distal edge of the activity from one radiation field becomes blurred inside the activity from an opposing radiation field. One way to overcome this challenge is to acquire the PET image after each or selected field by employing an in situ setup (Fiorina et al 2021). If imaging after the delivery of several fields, one proposed solution is only to consider patients with a maximum angle of 90 o between two beam directions . Another approach has successfully extracted the activity due to one field from the residual activity produced by a previous irradiation (Ferrero et al 2018a).
Given the fundamental differences between PET activity and dose distributions, combined with the limitations to using range analysis methods following irradiation by opposing fields, insight into the correlation between range shifts measured in PET activity and those in dose, in a variety of treatment scenarios, may be of interest. A measure of this correlation may be used for evaluating the precision of a particular PET scanner in detecting deviations in dose delivery during proton therapy in a broad range of patients. Furthermore, for a given scanner, knowledge of such a correlation may be used for the development of a prediction model which may serve as an assistance tool for clinical application.
The jagiellonian positron emission tomography (J-PET) scanner (Moskal et al 2021a(Moskal et al , 2021b, a novel PET scanner technology based on plastic scintillators, is being considered for proton range monitoring applications (Rucinski et al 2018, Baran et al 2019. The J-PET was originally proposed as an inexpensive technology for the construction of total-body diagnostic PET scanners (Moskal et al 2014, Niedźwiecki et al 2017, Moskal et al 2021c capable of multi-photon and positronium imaging (Moskal et al , 2021a(Moskal et al , 2021b. The same characteristic makes it an attractive option for implementation in proton therapy centres as a dedicated tool for range monitoring. In addition, the modular construction of the J-PET (Moskal et al 2016(Moskal et al , 2021c and the solely digital programmable and triggerless electronics (Korcyl et al 2018), facilitate its adaptation to either a full-ring, dual-head or other unconventional geometries required for in situ PET. Furthermore, the timing performance of the J-PET system, with 140ps CRT (FWHM) (Moskal et al 2021c) expected with BC-408 scintillator strips, in a strip geometry as used for this work, should allow the possibility of beam-on imaging in the long-term.
In this work, the application of the J-PET scanner for proton therapy range monitoring is investigated by means of a simulation study involving treatment plans for 95 head-and-neck as well as brain patients treated at the Cyclotron Centre Bronowice (CCB) in Krakow, Poland. An in-beam imaging protocol using a dual-head J-PET geometry and an in-room protocol using a full-ring geometry are considered. The sensitivity of the J-PET scanners to artificially introduced range shifts was investigated by means of the correlation of the measured shifts in PET activity to those in simulated dose. A method for using the correlation data from a wide range of treatment plans to form a prediction model, which may be utilized in the development of a tool for assisting clinical decision making, is presented and applied to the data obtained using the J-PET scanners. This study is an important step toward experiments with the J-PET prototypes and real patients undergoing proton therapy.

Materials and Methods
2.1. J-PET geometries and imaging protocols J-PET technology utilizes plastic scintillator strips for detecting 511 keV coincidence photons via Compton scattering (Moskal et al 2014). The strips may be arranged axially and read out at the ends, allowing for a multilayer geometry (Moskal et al 2016(Moskal et al , 2021c) that makes it possible to compensate for the lower efficiency of the plastic for detecting 511keV photons, compared to crystals such as LSO or BGO (Vandenberghe et al 2020). The current modular design consists of 13 rectangular strips of BC-404 plastic scintillator, each with a cross-section of 6 × 24 mm 2 and 500 mm in length, arranged side-by-side with 7 mm pitch, in detector modules read out on either end using silicon photomultipliers (SiPMs).
The J-PET scanner design incorporates a layer of wavelength shifter (WLS) strips, with a thickness of 3 mm, applied perpendicularly to the modules for the purpose of reconstructing the interaction position of the annihilation photon along the length of a scintallator strip (Smyrski et al 2014(Smyrski et al , 2017. Such an arrangement is expected to result in a spatial resolution of FWHM = 5 mm along the length of the strips (Moskal et al 2021c). For the simulated system, an energy resolution of E E E 0.044 MeV ( ) ( ) s = was assumed, as measured for a single strip J-PET prototype (Moskal et al 2014), as well as a conservative estimate of 500 ps coincidence resolving time (CRT) (Niedźwiecki et al 2017). Coincidence events were selected in a 3 ns time window and a 200-389 keV energy window, chosen to include the Compton edge for 511 keV photons. The spatial resolution along the length of the scintillator strips was modeled by blurring the axial position of each photon interaction by a Gaussian filter with FWHM = 5 mm, as per the expected resolution.
For the purpose of range monitoring in proton therapy, two scanner geometries were considered, as shown in figure 1: (a) a double-layer full-ring geometry for an in-room imaging protocol and (b) a triple-layer dualhead geometry mounted in situ for in-beam imaging.
The cylindrical J-PET scanner consisted of 48 modules, arranged in two layers, forming a regular 24-sided polygon with an inner diameter of D = 74.0 cm. The double-layer cylindrical geometry was chosen with costeffectiveness in mind. It has been shown through Monte Carlo simulations (Baran et al 2023) that the addition of a second layer to the current single-layer full-ring experimental prototype would increase system efficiency by a factor of three, while a third layer would result in a diminished increase of a factor of 1.7 over the two-layer geometry. The in-room imaging protocol, which typically uses a stand-alone PET scanner positioned within the treatment room (Zhu and Fakhri 2013), included a 1 minute delay following the delivery of the treatment plan, allowing time for the patient to be moved on a treatment coach from the treatment position into the scanner, followed by a two minutes PET scan.
The cylindrical J-PET scanner could also be utilized in an off-line protocol, for which it would be located outside of the treatment room for PET imaging following a delay of up to 30 min after the treatment (Zhu and Fakhri 2013). Such a protocol is interesting from the perspective of economics and clinical workflow aiming to reduce the treatment room occupancy. However, this work focuses on in-room PET due to several advantages it offers over the off-line approach. The in-room protocol allows for a shorter time gap between the delivery of the proton treatment and the PET acquisition, taking advantage of a larger signal from 15 O and reducing biological washout. Furthermore, the method of image analysis used in this work relies on having limited patient repositioning errors and anatomical changes between the treatment and PET positions, which are greatly reduced using in-room PET.
The dual-head geometry consisted of three layers of four side-by-side J-PET detector modules in the gantry position, with the scintillator strips running perpendicular to the beam direction so as to achieve maximum spatial resolution for measuring range. The detector heads were placed 40 cm apart, one in front of the patient's face and the other underneath the treatment couch. The triple-layer dual-head geometry corresponds to the current experimental prototype, which is under investigation at CCB. Monte Carlo simulation studies have shown that additional detector layers provide a diminishing benefit to the system efficiency (Baran et al 2023). For all plans, the first delivered field was directed from either the right or the left side of the head, within 30 degrees of horizontal, thus avoiding a collision of the gantry with the proposed dual-head geometry. Such a configuration was also capable of fully containing the largest cropped patient CT from the cohort, while avoiding overlaps between the detectors and the patient location. The in-beam protocol, where the PET acquisition is started immediately following the irradiation (Zhu and Fakhri 2013), consisted of a two-minute scan immediately following the delivery of the first field of the treatment plan using the dual-head scanner in the treatment position.
For both investigated protocols, the in-room and in-beam, the patients were positioned in such a way that the treatment isocentre was axially centered in the scanner.

Patient data from CCB proton therapy facility
The CCB proton therapy facility in Krakow (Poland) utilizes an IBA Proteus C-235 cyclotron equipped with two rotational gantries with dedicated scanning nozzles, an eye treatment room and an experimental hall. Patient data for this study was collected from the institution's anonymised clinical database from patients who underwent intensity-modulated proton therapy (IMPT) between 2016 and 2018 on both gantries (Garbacz et al 2021, McNamara et al 2022. All patients signed informed consent, allowing the use of treatment data for research purposes and statistical analysis. To evaluate the effectiveness and robustness of the J-PET systems for detecting proton range deviations, a set of treatment plans from 95 consecutive head and neck as well as brain patients was assembled. These included planning target volumes (PTV) ranging from 28.5 to 1010 ml (median volume of 185.5 ml), irradiated using 1378-32 290 pencil beams (median of 2843 pencil beams). There were 9 patients treated with two fields, 29 treated with 3 fields, 56 with 4 fields and one with 6 fields.There were 50 head and neck and 45 brain patients. The treatment plans contained 1 × 10 10 to 1.91 × 10 11 total protons.

Monte Carlo simulation and image reconstruction workflow
This study made use of the ProTheRaMon framework (Borys et al 2022), developed using Gate v9.1 (Jan et al 2011, Sarrut et al 2021. For simplicity, the framework may be considered as consisting of three main stages: (a) simulation of the delivery of the proton treatment and production of β + -emitting isotopes, (b) simulation of β + emission and PET acquisition and (c) image reconstruction.
The proton beam delivery was simulated using a beam model developed to match measurements taken at CCB (Gajewski et al 2021). The 95 plans from the clinical treatment planning system (TPS) were converted into Gate input files using a dedicated software tool, incorporating the beam model (Borys et al 2022). Each treatment plan was simulated on a phantom created from the corresponding planning CT image using a facility-specific clinical CT calibration curve obtained from stoichiometric calibration (Schneider et al 1996) and implemented in Gate (Gajewski et al 2021). Production maps for seven β + -emitting isotopes, produced on C, N, O, P and Ca, were generated: 10 C, 11 C, 13 N, 14 O, 15 O, 30 P and 38 K. In addition to simulating the prescribed treatment, deviations were artificially introduced into the phantom setup. To model day-to-day variations in patient positioning, offsets of 2, 4, 7, and 10 mm in the positive and negative x, y and z directions of the CT coordinate system were applied to each treatment plan. Additionally, to model uncertainties in the procedure of using the planning CT for estimating the RPSP of patient tissues, a positive and negative offset of 2% was applied to the RPSP calculated from each patient CT value. The total number of simulations per patient was 27, including the delivery of the prescribed treatment, 24 patient mispositioning errors and two calibration errors.
Simulating individual activity production and emission for each of 27 configurations was very computationally demanding. Even when running the simulation on a computational cluster, with a reduced fraction of primaries, using a large step size for particle tracking and large production cuts for particles and components of the geometry that were not of interest for the scoring of β + -emitters, the simulation could not be performed in an acceptable time-frame. In order to simulate isotope production in a large number of patients within a reasonable time-frame, the GPU-accelerated Monte Carlo code FRED (Gajewski et al 2021) was used for the majority of the first stage simulations. The FRED simulations used 10 5 primaries per delivered pencil beam, so that on average 1.55% of the total number of planned primaries were simulated for each field of the patient plans (McNamara et al 2022). In this way, a field consisting of 4 × 10 8 primaries could be simulated in approximately 5 min on two NVIDIA TITAN X GPUs, compared to 4 h required by Gate using 400 CPUs and 10% of the total number of planned primaries, as required in order to achieve acceptable statistics. FRED was previously validated against Gate by performing simulations of all of the treatment plans used in this study (McNamara et al 2022).
The positron emission and PET acquisition was simulated using Gate v9.1. The production maps of each of the seven β + -emitting isotopes obtained in the previous simulation step were used to generate β + activity sources, scaled according to the half-life of each isotope and the time structures of the treatment and imaging protocols. The β + energy spectrum of each isotope was simulated. To mimic the treatment protocol timing, a delay of 90 s was introduced between fields, assuming 60 s are required for irradiation and 30 s for each gantry rotation, values based on gantry rotation speed and irradiation duration times from CCB. The imaging protocol time structure is described in section 2.1. The medium for positron annihilation and propagation of the 511 keV annihilation photons was generated using the planning CT of each patient. The full-ring and dual-head J-PET geometries were simulated using activity distributions corresponding to the in-room and in-beam imaging protocols, respectively.
Image reconstruction was performed using the CASToR toolkit (Merlin et al 2018). Five iterations of the ML-EM algorithm were used to reconstruct all images, including time-of-flight (ToF) information with a 500 ps (FWHM) CRT. Coincidence events involving 511 keV photons that were scattered in the patient typically reach approximately 40% of the total and were removed from the data. A Gaussian filter with σ = 5 mm was used for post reconstruction smoothing of all images.

Biological washout
Correction for biological washout is an essential consideration for in vivo PET range verification of proton therapy. However, modeling the process is a major challenge. Well-established compartmental models used in conventional tracer imaging cannot be applied here due to the unknown molecular form of the positron emitters created during the irradiation (Parodi et al 2007b). Instead, models based on animal studies are generally used, typically in the form of a three-component model (Mizuno et al 2003), describing a fast, medium and slow washout component. Following the formalism of Ammar et al (2014), the three-component model describes the effective decay in a living subject as: where A phys is the physical activity of a radionuclide, decaying exponentially according to the half-life of the particular isotope, and C bio is the biological component, decaying according to The biological decay model for the remaining isotopes was biased towards the faster of the two models, namely that of oxygen. The washout model was applied in simulation by scaling activity maps to the value of average activity present during the PET acquisition, taking into account the time between the delivery of each field and the beginning of the PET scan, as well as the duration of the scan.
The effect of biological washout on the sensitivity of the J-PET scanners to dose range shifts was examined by including a washout model in the simulations of the post-treatment in-room and in-beam PET scans of 10 patients and comparing the results to those obtained when washout was ignored.

Range shift analysis
In order to investigate the sensitivity of the J-PET to deviations in proton range, PET images acquired after each artificially introduced range shift were analyzed in terms of how well they reflected the compliance of the delivered dose with the prescribed treatment. For this purpose, the methodology developed by Fiorina et al (2021) was adapted and expanded to include range analysis of the dose maps that were obtained from simulation. The analysis was implemented in Python using FREDtools (Gajewski 2022).
As a visual tool, beam's-eye-view (BEV) maps of the range shifts in PET activity were obtained for each altered plan for each patient using the in-room and in-beam protocols. Such maps may be overlaid with the patient CT to give a visual representation of regions where the activity range differs from that expected from the prescribed treatment. To demonstrate that such a representation reflects changes in dose range distributions, the corresponding BEV maps of the range shifts in dose were obtained for comparison. For the in-room protocol, some degradation of the correlation between the distal edges of the PET and dose distributions is expected, due to multiple and overlapping radiation fields. However, it can also be expected that a shift in the range of protons of a given field (directly reflected as a shift in the dose range) will also result in a shift in the profile of the distal falloff of the PET activity along the BEV of that field. The shift is expected to be most strongly correlated to that of the dose if there is a small angle between the delivered fields. The correlation should also be strongest along the last delivered field, because less of the PET activity due to that field will have decayed.
The formation of the BEV maps of range shifts, in both reconstructed PET activity and dose, differs somewhat in implementation from that proposed by Fiorina et al for PET images. First, for both dose and PET distribution analysis, a region of interest (ROI) was selected corresponding to the area containing a dose above 50% of the prescribed dose. A map of single-fraction dose was transformed to the BEV, maintaining the 2.5 × 2.5 × 2.5 mm 3 voxel size. For the in-room protocol, the cumulative dose from all fields was summed and the BEV chosen along the last delivered field. For the in-beam protocol, only the dose from the first field was considered, with the BEV being along that field. In BEV coordinates, where the beam direction is defined to be along the z¢ direction, the resulting dose distribution is D x y z , , BEV ( ) ¢ ¢ ¢ . Line profiles were then taken through D x y z , , BEV ( ) ¢ ¢ ¢ along z¢ at every voxel x y , ( ) ¢ ¢ . For the range measurement, ten thresholds t were defined in 1% intervals, centred at the standard 50% of the distal falloff of each profile. The depth z t ¢ , at which each profile crossed these thresholds was used to form a series of images R x y , t D ( ) ¢ ¢ , each representing the range distribution of an iso-dose surface at a given threshold t.
The same method was applied to the reconstructed PET images, but with the ten thresholds centred at 35% of the distal falloff of each profile. These thresholds were set significantly higher than those used by Fiorina et al (2021) in order to explore an activity level above the activity background created by preceding radiation fields. At the same time, the thresholds were kept low because positions near the tail of the PET distribution are more closely associated with the end of the range of the distal edge protons. The same thresholds were chosen for both in-room and in-beam protocols in order to remove threshold level as an additional variable in the comparison of the two methods. In this way, a series of images R x y , t PET ( ) ¢ ¢ of the range distribution of iso-activity surfaces were formed. The process of identifying R D and R PET is illustrated in figure 2 for a single pair of dose and PET profiles.
The procedure was repeated for all 27 irradiation scenarios for each patient. For both dose distributions and reconstructed PET images, maps of range differences between the prescribed treatment and the 26 altered treatment scenarios were obtained by averaging the difference images obtained at each threshold t. If we let R x y , ( ) D ¢ ¢ be the map of averaged range differences, for either the dose distribution or the reconstructed PET image, the calculation is as follows: ) ¢ ¢ represents the range from the prescribed treatment, at the respective thresholds t defined above, while R x y , t ( ) ¢ ¢ is the corresponding range resulting from an alteration introduced into the irradiation scenario.
In order to quantitatively relate J-PET measurements to deviations in proton range, the correlation between range shifts in dose and those in reconstructed PET activity was investigated. As a quantitative metric, the mean of range shift values ( R x y , ( ) D ¢ ¢ ) was considered appropriate, as it expresses the magnitude of the overall deviation in the position of the distal edge of both dose and PET activity. A prediction model was constructed where the mean of R x y , PET ( ) D ¢ ¢ , calculated from measured data, was used as the predictor of the mean of The linear regression model was tested by computing the coefficient of determination (r 2 ) for the individual patient cases, as well as for data from the entire cohort. In order to study the effect of tumor size on the predictive power of the model, the patient cohort was divided according to the target volume, and the correlation study was repeated for two individual groups of patients: those with medium (50-200 ml) and large (>200 ml) target volumes. In order to exclude the statistical effect of cohort size on the correlation results, each of the groups contained the same number of patients (39).
In order to assess the influence of statistical variations during PET imaging on the prediction of the mean proton range deviation, the PET simulations, image reconstruction and analysis were repeated for a group of 10 patients with the in-beam protocol.

Range analysis of a selected treatment plan
In order to illustrate the proposed protocols for range monitoring using the J-PET scanner geometries, and in particular the method for analyzing the reconstructed PET images, results are first presented for a patient arbitrarily chosen from the cohort using the in-room protocol. In this case, a brain patient treated with three fields with the gantry at 0, 90 and 320 degrees, and corresponding couch angles of 0, 0 and 90 degrees, is considered.
Reconstructed PET images, which serve as the input for the analysis, are shown in figure 3(a) as transverse, coronal and sagittal slices through the treatment isocentre, overlaid on the planning CT. The PTV is shown for reference. Also shown are the directions of the three treatment fields, the last to be delivered is highlighted in red to indicate the direction of the BEV analysis. The corresponding simulated dose distribution is shown in figure 3(b). Figure 3. Section of (a) the distribution of the PET activity generated by the simulated treatment plan and reconstructed from data acquired using the J-PET full-ring geometry, and (b) the dose deposited during the simulated treatment. The distributions are overlaid on the patient's planning CT. The arrows indicate the three radiation field directions. The red arrow indicates the last delivered field of the treatment plan, chosen as the BEV for further analysis of the images from the in-room protocol. The PTV is shown in green. Figure 4 illustrates the initial steps of the analysis method for the same patient. Figures 4 (a) and (  0.89 (ranging from 0.08 to 0.99) for the in-room and in-beam protocols, respectively. The distribution of r 2 values for individual patients is shown in figure 9 for both geometries.
For the individual patients for whom the study was repeated, the coefficient of determination for the linear regression showed a standard deviation 0.04%, or 4% of the mean, between independent realizations. This resulted in a change of 0.02 mm in the residual standard error of the prediction model for those patients and a change of 10% of the coefficient of determination for the cohort.
For the in-beam protocol, variations in the coefficient of determination and residual standard error between patients with large and medium sized tumor volumes were within the limits set by expected statistical variation. A larger than expected difference in the residuals was found for the in-room protocol, with medium sized tumors giving a value of 0.41 mm compared to 0.22 mm for large tumors. In this case, the difference between the coefficient of determination was within statistical limits. The relatively large difference in residuals indicates that the precision of the model may be influenced by PET activity due to opposing beams in the case of medium sized tumors using the in-room protocol.
The inclusion of a model for biological washout in the in-room protocol had minimal effect on the fit of the prediction model: r 2 decreased from 0.85 for the simulations that did not take washout into account, to 0.84 with washout. On the other hand, residual standard error increased from 0.22 to 0.32 mm with the inclusion of washout, a difference beyond what was expected from statistical variation. This suggests that washout has a noticeable effect on the prediction model generated with in-room PET data. However, the prediction power of the model degraded by washout is still good.

Discussion
The J-PET technology, configured in cylindrical and dual-head geometries, was investigated in terms of its effectiveness in detecting proton range shifts in a large selection of simulated patient treatment plans that included artificially introduced deviations from the prescribed treatment. To quantify the sensitivity of the proposed geometries to changes in proton range in a clinical setting, reconstructed PET images were compared to dose distributions by means of analyzing BEV range shift maps.
The range shift maps obtained from the reconstructed PET images may be used as a visual tool to aid in the identification of range shifts in dose that arise from deviations from the treatment plan. From the distributions in figures 5 and 6, it can be seen how the magnitude of range shifts in R x y , D ( ) D ¢ ¢ is reflected in that seen in R x y , PET ( ) D ¢ ¢ . Information regarding the approximate location of where the dose range shifts occur is also provided by R x y , PET ( ) D ¢ ¢ . The patient chosen for the analysis shows a representative PET-Dose correlation for the in-room protocol (r 2 = 0.95, same as cohort median) and better than average PET-Dose correlation for the in-beam protocol (r 2 = 0.95, cohort median r 2 = 0.89). The range difference map obtained from a PET scan acquired during or following a treatment delivery may be used to identify the specific regions where a deviation in proton range becomes critical for the treatment outcome. This may be done by overlaying the range difference map with the patient CT and identifying regions where an overshoot in proton range may lead to overdosing an organ at risk. Similarly, regions of undershooting may be assessed for the possibility of underdosing the tumour. For this purpose thresholds may be established to assess if a range deviation is significant. For example, in the work of Fiorina et al (2021), a compliance interval of −4 to 4 mm was reported, within which the detected range difference may be considered as a statistical fluctuation, with 95% confidence level. Thresholds may be set according to such values to select regions where the range deviations are clinically significant Figures 8(a) and (b) indicate that average shifts in proton range of up to ±5 mm have been measured in individual patients in this study. Using the residual standard error as a measure of the precision of the measured shift using PET imaging, values of 0.33 mm and 0.23 mm may be given here for the precision of the in-room/ cylindrical and in-beam/dual-head protocols respectively. Another study that correlated shifts in PET distribution to those in dose (Handrack et al 2017) measured shifts of up to ±3 mm well correlated to dose data to roughly within 1.8 mm.
The cohort-based prediction model described in this work may serve in the development of an assistance tool for clinical application. For a given PET measurement, the model may aid in the identification and quantification of proton range variations by providing an average value of the overall range difference, either after the delivery of an entire fraction (in-room protocol) or the first delivered field of that fraction (in-beam protocol). The measure of mean range difference with respect to the planned treatment may be used in the way of a red flag, as an initial indication that additional CT imaging and dose recalculation may be required. In the case the flag is raised, the BEV range maps may be consulted. If the model's predictions are confirmed by the BEV images and critical range differences are observed, further clinical decisions, including, e.g. therapy adaptation, can be made. The difference in r 2 between the in-room and in-beam models is partly due to statistical variations. However, several differences may be noticed in the shape of the residuals plots as well as the slope of the linear model. Such differences are likely the result of a convolution of a number of effects that stem from the different characteristics of the two protocols. It is possible that in the in-room protocol, the superior angular coverage and sensitivity of the full-ring J-PET geometry play a role. It should be noted that a 500 ps CRT results in spatial resolution of ∼75 mm, several times lower than the reconstructed voxel size, and does not completely compensate for missing projection angles. The differing delay times following the treatment in the two protocols result in different contributions to the signal from the various isotopes. Finally, overlapping radiation fields in the case of the inroom protocol will influence the PET-dose correlation. Overall, the correlation between dose and PET ranges in the BEV is strong for both of the PET scanner geometries and imaging protocols investigated. Future studies of clinical imaging protocols may include deconvolution of the effects of these individual characteristics. However, the clinical constraints for PET imaging in a busy schedule of daily proton therapy treatment should first be established by clinical staff.
Regarding uncertainties in the results presented in this work, it should be noted that these affect the reconstructed PET images and, to a lesser extent, the dose distributions. The main sources of uncertainty in the PET images come from the Monte Carlo simulations and include: uncertainties in the beam model (Gajewski et al 2021), uncertainties in the nuclear interaction models used in Geant4 (Kraan 2015), and uncertainties in modeling the J-PET detector response. Uncertainties in dose distributions have been discussed in Gajewski et al (2021). The impact of all the uncertainties mentioned above is somewhat mitigated by the fact that the results are based on the comparisons of two measurements performed using identical assumptions regarding beam model, nuclear cross-section, etc. The error originating from independent realizations of the PET scan simulation resulted in an uncertainty of 4% (standard deviation over mean) in the coefficient of determination of individual patients. This gave rise to a small change in the residual standard error of the prediction model of 0.02 mm. These uncertainties must be kept in mind when using the linear model for predicting mean proton range deviations.

Conclusions
The application of the J-PET technology as a tool for range monitoring in proton therapy has been investigated by means of a detailed Monte-Carlo study, involving treatment plans from a large cohort of patients. The ability of the scanner to detect deviations artificially introduced into the simulated treatment delivery was studied by considering full-ring and dual-head J-PET geometries.
A method for analyzing the PET images has been proposed, based on using the distribution of reconstructed PET activity to estimate the degradation of dose conformation cause by discrepancies between the prescribed and delivered treatment. Beam's-eye-view maps of range shifts in PET activity, which can be overlaid on the patient CT, indicate the degree and location of deviations in dose range. The mean of range deviation values in the PET maps was used to construct prediction models for obtaining a quantitative estimate of the degree of the deterioration of dose conformation at the distal edge. The ability to construct such models, using simulated dose maps and images obtained using the J-PET scanner, is evidence of the capability of the J-PET technology to detect small deviations in proton range.
For the prediction models, the mean value of dose range shifts was proposed as a convenient measure of the degradation of dose conformation. However, the prediction of more clinically relevant metrics based on the dose volume histogram (DVH), such as Homogeneity Index (HI), will be the subject of future studies that will include other tumour regions beyond the brain and head and neck. Furthermore, substantial motivation for bringing PET-based proton range monitoring into routine clinical practice would be provided by a correlation between parameters extracted from PET images and the clinical outcome of patients (i.e. pattern of failure of tumors). The possibility of extracting the failure pattern from a clinical data set, such as used in this study, should be explored.