In-vivo range verification analysis with in-beam PET data for patients treated with proton therapy at CNAO

Morphological changes that may arise through a treatment course are probably one of the most significant sources of range uncertainty in proton therapy. Non-invasive in-vivo treatment monitoring is useful to increase treatment quality. The INSIDE in-beam Positron Emission Tomography (PET) scanner performs in-vivo range monitoring in proton and carbon therapy treatments at the National Center of Oncological Hadrontherapy (CNAO). It is currently in a clinical trial (ID: NCT03662373) and has acquired in-beam PET data during the treatment of various patients. In this work we analyze the in-beam PET (IB-PET) data of eight patients treated with proton therapy at CNAO. The goal of the analysis is twofold. First, we assess the level of experimental fluctuations in inter-fractional range differences (sensitivity) of the INSIDE PET system by studying patients without morphological changes. Second, we use the obtained results to see whether we can observe anomalously large range variations in patients where morphological changes have occurred. The sensitivity of the INSIDE IB-PET scanner was quantified as the standard deviation of the range difference distributions observed for six patients that did not show morphological changes. Inter-fractional range variations with respect to a reference distribution were estimated using the Most-Likely-Shift (MLS) method. To establish the efficacy of this method, we made a comparison with the Beam’s Eye View (BEV) method. For patients showing no morphological changes in the control CT the average range variation standard deviation was found to be 2.5 mm with the MLS method and 2.3 mm with the BEV method. On the other hand, for patients where some small anatomical changes occurred, we found larger standard deviation values. In these patients we evaluated where anomalous range differences were found and compared them with the CT. We found that the identified regions were mostly in agreement with the morphological changes seen in the CT scan.

1 Introduction

Proton therapy is a type of radiation therapy that uses protons to treat cancer. The advantage of proton therapy with respect to conventional radiotherapy (X-rays and electrons) is related to the characteristic depth dose profile of charged particles (Bragg peak) (1). The accuracy of proton therapy strongly depends on the determination of the Bragg peak position. Uncertainties in the knowledge of the proton range can affect the dose distribution. These uncertainties include anatomical changes (physiological or morphological, organ motion, tumour regression, weight loss arising during the course of treatment), patient setup uncertainties, and range errors from uncertainties in CT Hounsfield units (HU), conversion of HU into particle stopping power, and reconstruction artifacts (2–6). In patients where anatomical changes are expected to occur, generally a control CT is acquired at some point during the treatment course (4–9). The scheduling of the control CT is variable and based on the clinical experience of the radiation oncologist. Based on the CT outcome, the radiation oncologist may decide for treatment replanning.

In-vivo range monitoring is desirable in order to support the radiation oncologist in the decision on when to perform a control CT (7, 10). One of the most consolidated monitoring techniques is Positron Emission Tomography (PET) (10–13). Nuclear interactions of the therapeutical beam with the tissue result in the production of β⁺ -isotopes, which decay emitting a positron, that annihilates into a 511 keV photon pair. The detection of these photon pairs by means of a PET system yields an activity image, that is indirectly correlated with the dose. Of all PET data acquisition modalities, in-beam (IB) data acquisition is generally considered the most attractive, providing real-time information about the treatment (10, 14, 15).

INSIDE (16) (INnovative Solution for In-beam Dosimetry in hadronthErapy) is a bimodal imaging system installed at the National Center of Oncological Hadrontherapy (CNAO), in Pavia, Italy (17). It consists of a particle tracker called Dose Profiler (18) and an IB-PET scanner (16). This bimodal architecture allows for the detection of annihilation photons with a PET detector, as well as for the detection of charged particles, also produced as a result of nuclear interactions. Since 2019 INSIDE is under clinical trial. During this trial (ClinicalTrials.gov ID: NCT03662373), we acquired IB-PET data for eight patients that received proton therapy treatments, fractionated in 6 weeks (about 30 sessions) along the entire course of their treatment. The first phase of the trial was completed in March 2020.

The goal of this study is twofold. First, we wish to investigate the level of experimental fluctuations in inter-fractional range differences observed in the INSIDE system, which drives the sensitivity to detect range differences. Such fluctuations can be due to the limited statistics of the acquired PET images, differences in irradiation and data acquisition conditions, PET cart setup errors, fluctuations due to random patient setup errors, and so on. This will be done by studying inter-fractional range shifts for patients that had no morphological changes, using the Most-Likely-Shift (MLS) from Frey et al. (19), where it was applied to offline PET data. To establish the efficacy of the MLS method for IB-PET images, we made a comparison with the Beam’s Eye View method (20). Second, we will compare the results with patients that did show morphological changes, and we will investigate whether it is possible to identify the regions affected by the changes with the MLS method. It must be noted that this is the first study that includes all the available data of the patients treated with proton therapy and monitored with INSIDE. It is also the first time that the MLS method will be applied to IB-PET data.

2 Materials and methods

2.1 The INSIDE in-beam PET scanner

The design, construction and initial clinical tests of the INSIDE IB-PET system are described in detail in previous works (16, 21–24), and only the most relevant information is given here. The system consists of two planar heads placed at 30 cm from the isocenter, for a total distance between them of 60 cm. Each head has an active area of 10 × 25 cm² and consists of 2 × 5 modules produced by Hamamatsu Photonics (25). Each module is a square matrix of 16 × 16 scintillating crystals of lutetium fine silicate (LFS) (26), each with dimension 3 × 3 × 20 mm³. Thus the total number of channels in one PET head is 2560. LFS is a commercial name of the set of Ce-doped silicate scintillation crystals (26, 27). The density (7.35 g/cm³) and light yield (80% with respect to NaI:Tl) are comparable to lutetium oxyorthosilicate and LYSO but with improved time performance: it has a ~36-ns decay constant. These crystals are optically coupled one-to-one to silicon photomultipliers (SiPMs). Their temperature was maintained at 18°C with the help of an integrated cooling system.

The 2 × 2560 PET detector channels are read by front-end electronics (FE), based on 64 TOFPET ASICs channels (28), which gives the time stamp of the event encoded through a time-to-digital converter with a resolution of 50 ps. The energy information was obtained with the time over threshold (TOT) method. The signals from the FE ASICs are processed locally by Field Programmable Gate Array (FPGA) boards. The channel dead time is imposed by the FE and is below 300 ns for a 511 keV event (28). The energy resolution for 511 keV photons was ~13%, as measured with a β⁺ source of ²²Na. The measured coincidence time resolution (sigma) was 450 ps. The acquisition is performed all along the treatment during both irradiation (spill) and beam pauses (inter-spill). However, in this analysis we have only used the data acquired during the inter-spill beam pauses, as well as data acquired a short time after the delivery of the field (see Section Pre-processing).

The 3D image reconstruction is performed on-the-fly with an iterative multicore Maximum Likelihood Expectation Maximization (MLEM) (29) algorithm, and its results are provided online during the treatment. Advanced analysis of the PET images is done offline (see Section Analysis). The reconstructed image FOV is 224 × 112 × 264 mm³, with a voxel size of 1.6 × 1.6 × 1.6 mm³. The IB-PET reference system is reported in Figure 1A, where the z axis is parallel to the beam direction. Figure 1B shows a PET image in the three main projections. Due to the partial angular coverage of the PET system, the images suffer from artefacts in the direction perpendicular to the planes (30), as well as from limited statistics [see also (31)].

FIGURE 1

Figure 1 Reference frames. (A) Sketch of the INSIDE setup and reference frame. (B) Example of the three orthogonal projections of a reconstructed PET image for patient 002P (see 2.3). (C) The same reconstructed PET image in the CT reference frame. In (B, C) the beam axis is indicated with the red arrow. (D) The color-scale of the PET images.

2.2 Dose delivery at CNAO

At CNAO protons are delivered with a modulated pencil beam scanning technique in a fixed beam-line, while the patient couch can rotate. The Dose Delivery System (DDS) (32) guides and monitors the proton beams accelerated by the synchrotron, and distributes the dose with a full 3D scanning technique. The dose results from the superposition of a large number of beams conveyed to the tumour volume. The irradiated volume is divided into several energy layers (slices), each of which is irradiated by various overlapping iso-energetic beams (spots) arranged in a grid of 3 mm pitch. A typical proton therapy treatment has 20 to 30 fractions. Each treatment fraction is delivered through two or more fields, each consisting of many thousands of spots to cover the target. For more details about the beam specifications, scanning technique and the CNAO particle accelerator we refer to previously published works (17, 32).

The DDS provides a log-file in which the information about the actually delivered particles is stored, including the number of delivered primaries, energy and the lateral position of the spot. The information from this log-file was converted into a binary mask in the same reference frame and with the same voxel size as the PET images (see Figures 1A, B), in which the non-zero values define a volume-of-interest, where PET activity can be present. This volume-of-interest will be referred to as the expected activity mask. The dose delivery time for the patients was typically a few minutes. The temporal structure of the beam extraction was 1 s of spill, during which the protons are delivered, and 2 s of pause between spills (inter-spill periods). The INSIDE PET system is mounted on a cart that is positioned before data acquisition. The position accuracy is about 1 mm, thanks to a dedicated cart positioning system. The presence of the cart has limited the number of patients that could be enrolled in the trial: it was only possible to monitor patients with beam angles without mechanical incompatibilities of the INSIDE setup with the patient couch movements.

2.3 Patient dataset

We analyzed the IB-PET data of the patients treated with protons. A patient treated with protons at the CNAO typically receives a total dose of 60 Gy divided into fractions of 2 Gy each.

Table 1 reports the patients analyzed, with the patient ID, the pathology, the number of fields N_fields, the number of treatment fractions delivered N_frac, the number of treatment fractions with PET data that could be included in the analysis, N_PET, the control CT outcome, as evaluated by the radiation oncologist, and the monitored treatment angle. The control CT is acquired approximately during the third week from the start of treatment. Two of the monitored patients (006P and 007P) showed morphological changes, nevertheless, they were not replanned because the treatment plan quality was still deemed acceptable. In this work we used the PET data for each patient in Table 1, that were acquired during the first field delivered. In this way contamination of the activity from the other fields was avoided. The choice of the field was made by clinical personnel.

TABLE 1

Table 1 Patients treated with protons analyzed in this work, with patient ID, pathology, number of fields, number of fractions delivered, number of fractions monitored, the control CT outcome, as evaluated by the radiation oncologist, and the monitored treatment angle.

Figure 2 reports the planning CT scan (upper row) and the control CT (lower row) of patient 006P. Figure 3 reports the same for patient 007P. In this case there are two anatomical modifications: a mild emptying in the GTV for a possible tumor response and also some inflammation reaction (mucositis) in the CTV within the nasal cavity.

FIGURE 2

Figure 2 Planning CT (upper row) and control CT (lower row) for patient 006P. (A) Axial, (B) sagittal and (C) coronal views for the planning CT, and (D) Axial, (E) sagittal and (F) coronal views for the control CT. The morphological change is highlighted with a red arrow. Also, the CTV (Clinical Target Volume) area is highlighted with a green line.

FIGURE 3

Figure 3 Planning CT (upper row) and control CT (lower row) for patient 007P. (A) Axial, (B) sagittal and (C) coronal views for the planning CT, and (D) Axial, (E) sagittal and (F) coronal views for the control CT. The morphological change is highlighted with a red arrows and the CTV (Clinical Target Volume) area is highlighted with a green line.

2.4 Pre-processing

In this work we used the IB-PET data acquired over the time-interval from the start of irradiation until 10 seconds after the end of the irradiation, whereby during irradiation only the inter-spill data were used. Although at the CNAO synchrotron the exact number of spills used for a given treatment is quite stable, small variations are possible from day-to-day, depending on the accelerator conditions. Therefore the reconstructed activity can depend on the beam conditions. For the same total dose delivered, the induced activity can change somewhat if the dose delivery time is different. An example of this is shown in Figure 4, where two reconstructed PET images acquired in two different days are displayed. Here the total delivered dose was the same, but the dose delivery and acquisition times were different. These differences cannot be simply corrected by normalizing, because the PET signal is time dependent. Therefore we used only those fractions with similar temporal profiles in the analysis. The number of such fractions is given in the fifth column of Table 1.

FIGURE 4

Figure 4 Coronal views of the planning CT of patient 004P: reconstructed PET images corresponding to the monitored fractions 1 and 3. The left and right image had different dose delivery and acquisition times, resulting in different activities. Fraction 3 had to be excluded because of the larger activity compared to the other fractions.

To reduce salt-and-pepper noise in the images we applied a 5 mm radius median filter. According to the couch angle, the PET images could be rotated into the CT imaging coordinate system expressed in International Electrotechnical Commission (IEC) standard. Figure 1C gives an example of a rotated image.

2.5 Analysis

The analysis was divided into two parts. First, we wished to investigate the level of experimental fluctuations in inter-fractional range differences observed by the INSIDE system, referred to as sensitivity. This was done by studying inter-fractional range shifts for patients that had no morphological changes (see Table 1). Second, we compared these results with patients that did show morphological changes. For the latter patients, we also verified whether the region of the observed range differences is correlated to the presence of real morphological changes in the CT scan.

The strategy that we followed for the analysis is based on an interfractional comparison between a reference PET image and an image acquired during a subsequent fraction j. In order to mitigate statistical fluctuations, the reference image was calculated as the average of the first two monitored fractions.

2.5.1 Most-Likely-Shift method

To compare the PET images of the subsequently monitored fractions with the reference image, we used the Most-Likely-Shift method, originally proposed by Frey et al. for off-line PET monitoring images (19) and recently applied also in (33). This is the first time that the MLS method was applied to IB-PET images. Given that such images suffer from artifacts due to the partial angular coverage and limited (see Figure 1), it is not a-priori obvious whether the MLS method works.

In summary, given a PET image of a monitored fraction j, for each pair of coordinates (x,y) in the transverse plane (in the INSIDE reference system), we compared the 1-D activity profile along the beam direction z in their distal fall-off with the corresponding reference profile.The MLS method considers profiles belonging to two different treatment fractions for the same pair of coordinates in the transverse plane (x,y) and calculates the most probable shift necessary to align these two profiles, normalised at their maximum, along the beam direction. We included only those (x,y) pairs into the analysis that were part of the expected activity mask. Moreover, to exclude profiles with low activity from the analysis, we included here only those profiles with an integrated activity of at least 30% of the profile with highest integrated activity. The reader is referred to the original paper by Frey et al. (19) for a detailed description of the MLS method.

We applied the method exactly as described in Frey et al., with three exceptions. First, the value for z_min(x,y) (the z-value of that activity profile where to start the analysis) was set at 15% of the maximum of the normalized reference profile at coordinates x,y. Second, z_max(x,y), i.e., the z-value where to end the analysis, was set at 2% of the maximum of the normalized reference profile. These modifications were done in order to focus more on the distal fall-off part of the profiles. Third, the shift value δ was limited between -16 and 16 mm, with steps of 1 mm. An algorithm was implemented that provided for each x,y pair the optimal shift distance δ_MLS. The calculation of the range difference δR_MLS is done by minimizing the absolute profile differences in the distal part (z_MLS ≤ z ≤ z_max) of two activity depth profiles, shifted against each other, as follows (19, 34):

$\begin{array}{ll}\delta R_{MLS}=\arg \underset{{\delta }_{MLS}}{\min }\left({\displaystyle \sum _{z_{MLS}}^{z_{\max }}\left|A_{ref}\left(z\right)-A_j\left(z-{\delta }_{MLS}\right)\right|}\right)& (1)\end{array}$

with A_ref and A_j corresponding to the reference activity profile and that of the fraction to be compared, respectively.

Then, a distribution of all the δR_MLS(x,y)values was created for each fraction j. These distributions were fitted with an asymmetric skew-normal distribution and its mean, μ_MLS , and standard deviation, σ_MLS , were evaluated for all monitored fractions j of a given patient. For each patient, we then evaluated:

● the average of μ_MLS over all the fractions. This average is denoted by Δ_MLS, giving an indication of the size of typical inter-fractional range shifts between the reference image and the subsequent fractions for this patient;

● the average of σ_MLS over all the fractions, indicated by ∑_MLS, giving an indication of the fluctuations in such shifts for a given patient.

Then, to investigate the level of fluctuations in our PET data in absence of any morphological changes, we extracted for the six patients without any morphological changes the average of ∑_MLS over these six patients, yielding〈σ_MLS〉. Such inter-fractional range fluctuations are accidental and can be due to the low statistics of the images, patient setup-errors, INSIDE cart setup errors, small differences in irradiation conditions, and so on.

Then, for each patient and each analyzed fraction, a three-dimensional outliers map O_MLS(x,y,z) of anomalously large shifts was created in order to visualise the largest range difference obtained. The map O_MLS(x,y,z) was filled for each voxel with coordinate x,y,z, which had z ≤ z_min as:

$\begin{array}{ll}\{\begin{array}{l}{O}_{\text{MLS}}\left(x,y,z\right)=\delta R_{\text{MLS}}\left(x,y\right),\text{if }\delta R_{\text{MLS}}\left(x,y\right)\le -2·\langle {\sigma }_{MLS}\rangle \text{or\hspace{0.17em}if\hspace{0.17em}}\delta R_{\text{MLS}}\left(x,y\right)\ge +2·\langle {\sigma }_{MLS}\rangle \hfill \\ O_{\text{MLS}}\left(x,y,z\right)=0\text{\hspace{0.17em}otherwise}\hfill \end{array}& (2)\end{array}$

For z > z_min (close to the end of range), O_MLS(x,y,z) = 0. Thus, the outliers maps are represented as 3-dimensional distributions, but each point along the z-axis was filled by copying δR_MLS(x,y) from z=0 up to z = z_min. These maps graphically identify the position of the distal fall-off part of the activity distribution, and they highlight the entire region along a pencil beam path that is possibly affected by a range shift with respect to the reference. These maps were re-oriented on the patient’s CT refersence frame, in order to verify if the real morphological change seen in the CT was identified.

The implementation was validated as in (19), i.e., by creating artificially shifted PET distributions and verifying that the MLS analysis correctly identified the introduced shift.

2.5.2 Comparison with Beam’s Eye View method

To confirm the validity of the Most-Likely-Shift method for our PET images, we compared the result with another existing range verification method that performs interfractional comparisons based on 1-D activity profiles: the Beam Eye View method. Details can be found elsewhere (20, 22, 35). This method is based on a multi-threshold approach to extract iso-activity surfaces. Thresholds from 2% up to 8% on the maximum of the entire PET image at 0.5% steps were considered. Briefly summarizing, similar to the MLS method, for each pair of coordinates x,y in the transverse plane (in the INSIDE reference system), the profile belonging to treatment fraction j was compared with the profile of the reference image A_ref. For a given threshold t, the shift between these two profiles, δR_BEV,t(x,y), was defined to be the difference between the maximum depth value reached at threshold for the reference profile, $R_t^{ref}\left(x,y\right)$ and that for a monitored fraction, $R_t^j\left(x,y\right)$. To obtain the final value of the range difference between the profiles at a given x,y point, δR_BEV(x,y), the average was taken over the thresholds. In other words:

$\begin{array}{ll}\delta R_{\text{BEV}}\left(x,y\right)=\frac{1}{N}{\displaystyle \sum _{t=1}^N\delta R_{\text{BEV},\text{t}}\left(x,y\right)=\frac{1}{N}}{\displaystyle \sum _{t=1}^N\left(R_t^{ref}\left(x,y\right)-R_t^j\left(x,y\right)\right)}& (3)\end{array}$

where N is the total number of thresholds considered, which was 13 in our case.

The remaining procedure is exactly as done for the MLS method. For each patient a distribution, containing all the δR_BEV(x,y) value, was created for each fraction j. These were fitted with an asymmetric skew-normal distribution, yielding the mean μ_BEV and standard deviation σ_BEV , for each monitored fraction j. For each patient we then extracted:

● the average of μ_BEV over all the fractions, denoted by Δ_BEV.

● the average of σ_BEV over all the fractions, indicated by ∑_BEV.

We extracted for the six patients without any morphological changes the average of ∑_BEV, yielding <σ_BEV> .

For each patient and each analyzed fraction, a three-dimensional outliers map O_BEV (x, y, z) was created, whose interpretation is the same as for the MLS maps. For each voxel (x,y,z), the map was filled when $z\le R_t^j\left(x,y\right)$, with t=8%, as:

$\begin{array}{ll}\{\begin{array}{ll}{O}_{\text{BEV}}\left(x,y,z\right)=\delta R_{\text{BEV}}\left(x,y\right)\hfill & \text{if} \delta R_{\text{BEV}}\left(x,y\right)\le -2·{\mathrm{<\sigma }}_{BEV}>\text{\hspace{0.17em}or\hspace{0.17em}if\hspace{0.17em}}\delta R_{\text{BEV}}\left(x,y\right)\ge +2·<{\sigma }_{BEV}>\hfill \\ O_{\text{BEV}}\left(x,y,z\right)=0\hfill & \text{otherwise}\hfill \end{array}& (4)\end{array}$

For $z>R_t^j\left(x,y\right)$, with t=8%, we put O_BEV(x,y,z) = 0. Then the O_BEV(x,y,z) maps were superimposed to the CT images, so that they could be compared with the O_MLS(x,y,z) maps as well as with the CT.

The BEV implementation was validated by creating artificially shifted distributions, and verifying that the analysis correctly identified the shift.

3 Results

In Figure 5A we show an example of a distribution of range differences δR_MLS with , respect to the reference for a few of the available fractions in a patient that does not change (005P). Each histogram entry represents a value of δR(x,y). The same is given for the BEV method in Figure 5B.

FIGURE 5

Figure 5 Examples of the distributions of inter-fractional range differences δR(x,y) of various fractions with respect to the reference (M) in patient 005P for the MLS method (A) and the BEV method (B). The dashed lines represent the 95% confidence interval tailored for, the MLS (5.0 mm) and the BEV method (4.6 mm) as 2<σ>.

Figure 6A gives the interfractional standard deviation σ_MLSfor each patient with MLS method as a function of the fraction compared with the reference. Squares and triangles stand for the standard deviation obtained for the set of unchanged and changed patients, respectively. For the patients that do not change (squares), we observe that the values for σ_MLS are somewhat lower than 3 mm, the spatial resolution of the INSIDE PET scanner (16). Moreover, we see for these patients that the standard deviation is roughly constant during the course of the treatment. For patients subject to morphological changes (triangles), the standard deviation values mostly exceed 3 mm. In particular for patient 007P, we see that after the 14th fractions σ_MLS is larger than 5 mm.

FIGURE 6

Figure 6 Inter-fractional standard deviation of the range difference distribution as obtained with the (A) MLS and (B) BEV method for each patient (y-axis) as a function of the number of the fraction compared with the reference (x-axis). The triangles stand for those patients where morphological changes where identified in the control CT, while the squares represent the patients that did not show changes.

Figure 6B is the same as Figure 6A, but now for the BEV method. We see that the results are in agreement with Figure 6A: the unchanged patients (squares) had generally a lower standard deviation than the changed patients (triangles), the values are mostly below 3 mm, and the values for the patients that changed are larger.

Table 2 summarizes the results for each patient obtained by the MLS method in terms of average range activity difference observed over all the fractions with respect to the reference. In the first, second and third column, we report the patient ID, Δ_MLS and ∑_MLS, respectively. The reported error for Δ_MLSand ∑_MLS is the standard deviation calculated over the various fractions, considering all fractions as independent measurements. The average over the fractions for all the patients that did not present anatomy changes, <σ_MLS> (see Section Most-Likely-Shift method), was calculated to be 2.5 mm. For the two patients that changed, the values for ∑_MLS were found to be larger.

TABLE 2

Table 2 For each patient, the average inter-fractional range difference Δ_MLS and average standard deviation ∑_MLS obtained with the MLS method, together with the corresponding values Δ_BEV and ∑_MLS for the BEV method.

The corresponding values for the BEV method, Δ_BEV and ∑_BEV, are given in the fourth and fifth column of Table 2, respectively. The average over the fractions for all the patients that did not present anatomychanges, <σ_BEV> (see Section Comparison with Beam’s Eye View method), was calculated to be 2.3 mm. Again, for the two patients that changed, the values for ∑_BEV were found to be larger than those for patients that did not change.

To establish whether the group of changed patients was statistically different from the group of unchanged patients, we performed the Student’s t-test over the values of Δ and ∑. While the Δ values were not statistically different (for the MLS method t = 0.27, p = 0.79, and for the BEV method t = 0.67, p = 0.53), for the ∑ values we found a significant difference (for the MLS method t = -3.90, p = 0.008, and for the BEV method t = -4.19, p = 0.006).

Figure 7A displays the outliers O_MLS(x,y,z) maps obtained for patient 006P by the MLS method. This is a patient which was subject to small morphological changes, as observed by the radiation oncologist and demonstrated in Figure 2. The colored area represents the pencil beam paths that had anomalous range differences with respect to the reference, see Eq. 2. The red color indicates a positive range difference with respect to the initial situation (overshoot). We see that the red color becomes darker as the number of the treatment fractions increases, strongly indicating a range overshoot that increases along the treatment course. In the penultimate fraction (fraction 33) of the treatment course, we observe range differences as large as 16 mm, i.e., a substantial range overshoot. The cause of the range overshoot becomes clear when looking at Figure 2. Comparing Figures 2A–C (planning CT) with Figures 2D–F (control CT), a cavity is seen that is emptied in the control CT. The fact that the regions highlighted in Figure 7A cover those of the anatomical change in the control CT is a strong indication that our PET images are sensitive to such changes. We therefore believe that the MLS method can be successfully applied to IB-PET images to detect morphological changes like cavity emptying.

FIGURE 7

Figure 7 Coronal views of the control CT of patient 006P, with the (A) MLS outliers map O_MLS(x,y,z) and the (B) BEV outliers map O_BEV(x,y,z) superimposed, obtained by comparing fraction 17, fraction 25, and fraction 33 with the reference. The colored areas are the pencil beam paths that lead to anomalous range differences with respect to the reference, with red indicating a range overshoot with respect to the initial situation. The overshoot is especially pronounced in fraction 25 and 33. The colormaps were obtained as described in 2.5.1 and 2.5.2 for the MLS and the BEV method, respectively.

Figure 7B shows the outliers maps O_BEV(x,y,z) obtained for patient 006P by the BEV method. The colored areas are seen to be similar to those of Figure 7A. Thus, the MLS and BEV method are roughly in agreement.

In Figure 8A we display the outliers maps O_MLS(x,y,z) for patient 007P by the MLS method. In this case, the map shows both range overshoots with respect to the reference (identified by the red color) and undershoots (identified by the blue color). By comparing Figure 3 with Figure 8A, we see that the region where the cavity has emptied is only weakly highlighted (orange). The blue zones that are identified may correspond to inflammation effects (see discussion).

FIGURE 8

Figure 8 Coronal views of the control CT of patient 007P, with superimposed the (A) MLS outliers map O_MLS(x,y,z) and (B) BEV outliers map O_BEV(x,y,z) , obtained by comparing fraction 14, fraction 17, and fraction 24 with the reference. The colored areas are the pencil beam paths that lead to anomalous range differences with respect to the reference, with red and blue indicating a range overshoot and undershoot with respect to the initial situation, respectively. Most clear are the beam undershoots (the blue areas). The colormaps were obtained as described in 2.5.1 and 2.5.2 for the MLS and the BEV method, respectively.

Finally, in Figure 8B we demonstrate the outliers maps O_BEV(x,y,z) , obtained for patient 007P by the BEV method. Also in this case the zones that are highlighted are those affected by beam overshoots (identified by the red color) and undershoots (identified by the blue color) with respect to the reference. They are similar to the regions that were identified with the MLS method, but the colored regions are somewhat larger (see discussion).

4 Discussion

In this work we analyzed for the first time all of the available IB-PET patient data for proton therapy treatments monitored during the INSIDE clinical trial. To evaluate the experimental level of inter-fractional fluctuations in IB-PET images in absence of anatomical changes (sensitivity), we studied six patients that did not show morphological changes. Using the MLS method we found for these patients an overall standard deviation <σ_MLS> in activity range difference of 2.5 mm. This is smaller than 3 mm, which is the spatial resolution of the INSIDE PET scanner  (16).

Regarding the observed differences between changed and unchanged patients, we saw that the average standard deviations ∑_MLS were larger for changed patients than those for unchanged patients. Moreover, the outliers maps in Figures 7A, 8A include approximately the regions that were also seen to change in the control CT. In patient 007P, the situation is complicated because the beam goes through highly heterogeneous tissue. Still, the results are encouraging, given the small size of the anatomical changes (order of a few ml).

We evaluated to what extend the values in Table 2 for the MLS analysis were affected by the choice of the parameters chosen (see Section Most-Likely-Shift method). Varying the threshold value for inclusion of the profiles from 30% (default) to 20% and 35%, the results did not change significantly. Also, varying the z_min (default 15%) between 10% and 20%, no significant range differences were observed. Regarding z_max (default 2%), this value should be as small as possible in order to include as much activity in the fall-off region as possible. However the noise level was found to be 2% so using lower values resulted in inconsistent results. Thus, the results were sufficiently robust.

Comparing the MLS and BEV methods, we saw that the MLS method gave results that were similar to those obtained with the BEV method for all patients. Thus, the methods are mostly in agreement with each other. An example where the methods gave slightly different results was in patient 007P (see Figures 8A, B), where the identified zones were partly different. The advantage of the MLS method with respect to the BEV method is that it evaluates point by point the most suitable shift to align two profiles, and does not only look at the most distal values like the BEV method.

We also verified that our results for the BEV method were in agreement with those published previously for three of the eight patients (20). Small differences could be attributed to differences in the fit procedure and in the choice of the reference PET image, which was in the case of (20) defined to be the first acquired PET fraction, while we took the mean of the first two acquired images as reference.

Several aspects in the analysis deserve more discussion. First, the number of analyzed patients and the number of fractions that could be included in the analysis was small (Table 1). The differences between changed and unchanged patients could in principle have been caused by other factors, including the small statistics of the PET images, differences in patient setup, tumor type, tumor size, depth of the irradiated region, irradiation, etc. The influence of these sources of uncertainty should be investigated in more detail by monitoring more patients. At the same time, the outliers maps agreed approximately with the zones that were seen to change in the control CT. Thus, we believe that the INSIDE IB-PET images can give indications about inter-fractional range differences. The MLS and BEV methods can both be applied for this purpose. This was especially visible in Figures 7, 8, where the colors indicating beam overshoots and undershoots become more intense as the number of treatment fractions increased. However, the effect is small and more patient data are needed to confirm whether such trends through the course of treatment are observable.

Second, it must be noted that we monitored only one of the fields, i.e., the first field delivered. In some cases we had N_fields = 3 (see Table 1), resulting in a very limited overall statistics of the PET signal. This can possibly be improved by assuring that the field delivered first is the one with largest activity, or by combining fields (36).

Finally, the size of the morphological changes of patient 6 and 7 was small (order of a few ml). In fact, these changes were considered small enough to not require replanning, since the recalculation of the treatment plan on such modified anatomy did not show any clinically significant modification of the DVHs both for target and OARs. Increasing the number of patients with larger expected changes is foreseen for the second phase of the INSIDE clinical trial. This can better confirm the validity of the INSIDE in-beam PET system in detecting morphological changes.

Despite the small number of patients and fractions that could be included in this study, the above results are encouraging. For future in-beam PET data acquisitions, we suggest to only compare PET images which are acquired under approximately the same irradiation conditions, to monitor the first field delivered to avoid contamination from other fields, and to monitor as many fractions as possible.

5 Conclusion

In this work we performed an inter-fractional range difference analysis including all the available patient data for proton treatments acquired during the INSIDE clinical trial: six patients not subject to morphological change, and two patients subject to morphological change. We applied the Most-Likely-Shift (MLS) method for detecting inter-fractional range differences, which, to the best of our knowledge, was the first time it was applied for IB-PET images. It was compared with the Beam Eye View (BEV), that was previously applied to a subset of the patients included in the clinical trial.

When putting together all fractions and all patients that did not have morphological changes, the standard deviation of the range variations in activity profiles, <σ_MLS>, was found to be 2.5 mm with the MLS method. The corresponding value for the BEV method was <σ_BEV> =2.3 mm. On the other hand, for patients where anatomical changes occurred, we found larger standard deviation values.

For the patients with anatomical changes we created outliers maps to indicate anomalous range variations. These were superimposed on the control CT, and the regions affected by an anomalous range difference approximately covered those with real anatomical variations observed in the control CT. The results are encouraging and suggest that the MLS method can possibly be used as a support tool by clinical personnel to detect anomalous situations during the treatment course and to guarantee the effectiveness of the treatment plan.

Data availability statement

The data that support the findings of this study are partly available on request from the corresponding author. The data are not publicly available due to privacy or ethical restrictions.

Ethics statement

The studies involving human participants were reviewed and approved by Comitato etico del Policlinico San Matteo di Pavia. The patients/participants provided their written informed consent to participate in this study.

Author contributions

The analysis and manuscript was prepared by MMog, ACK and MGB. All authors have contributed to the publication, being variously involved in the design and the construction of the detectors, in writing software, calibrating subsystems, operating the detectors, acquiring data and finally analysing the processed data. The INSIDE Collaboration members discussed and approved the scientific results reported in the submitted document. The work was subject to an internal collaboration-wide review process.

Funding

The INSIDE project was funded by the Italian Ministry of Education (PRIN MIUR 2010P98A75, 2013–2016), the Italian Institute of Nuclear Physics (RDH and INFN-RT2 PETRA 172800 projects, since 2016), the Historical Museum of Physics and the Enrico Fermi Study and Research Center, the Tuscany Government (POR FSE 2014–2020, INFN-RT2 PETRA 172800 Project, 2018-2019), and the CNAO Foundation (INSIDE2, since 2017).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. All authors contributed to the article and approved the submitted version.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Loeffler J, Durante M. Charged particle therapy–optimization, challenges and future directions. Nat Rev Clin Oncol (2013) 10:411–24. doi: 10.1038/nrclinonc.2013.79

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Paganetti H. Proton therapy physics. Boca Raton, Florida:CRC Press (2012).

Google Scholar

3. Paganetti H, Botas P, Sharp G, Winey B. Adaptive proton therapy. Physics in Medicine & Biology 66 (2021). doi: 10.1088/1361-6560/ac344f

CrossRef Full Text | Google Scholar

4. Albertini F, Matter M, Nenoff L, Zhang Y, Lomax A. Online daily adaptive proton therapy. Br J Radiol (2020) 93:20190594. doi: 10.1259/bjr.20190594

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Placidi L, Bolsi A, Lomax AJ, Schneider RA, Malyapa R, Weber DC, et al. Effect of anatomic changes on pencil beam scanned proton dose distributions for cranial and extracranial tumors. Int J Radiat Oncol Biol Phys (2017) 97:616–23. doi: 10.1016/j.ijrobp.2016.11.013

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Cubillos-Mesías M, Troost EG, Lohaus F, Agolli L, Rehm M, Richter C, et al. Including anatomical variations in robust optimization for head and neck proton therapy can reduce the need of adaptation. Radiother Oncol (2019) 131:127–34. doi: 10.1016/j.radonc.2018.12.008

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Mackay R. Image guidance for proton therapy. Clin Oncol (2018) 30:293–298. doi: 10.1016/j.clon.2018.02.004

CrossRef Full Text | Google Scholar

8. Morgan H, Sher D. Adaptive radiotherapy for head and neck cancer. Cancers Head Neck (2020) 5. doi: 10.1186/s41199-019-0046-z

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Kraan AC, van de Water S, Teguh DN, Al-Mamgani A, Madden T, Kooy HM, et al. Dose uncertainties in impt for oropharyngeal cancer in the presence of anatomical, range, and setup errors. Int J Radiat Oncol Biol Phys (2013) 87(5):888–96. doi: 10.1016/j.ijrobp.2013.09.014

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Parodi K, Polf J. In vivo range verification in particle therapy. Med Phys (2018) 45:e1036–50. doi: 10.1002/mp.12960

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Enghardt W, Crespo P, Fiedler F, Hinz R, Parodi K, Pawelke J, et al. Charged hadron tumour therapy monitoring by means of pet. Nucl Instruments Methods Phys Res Section A: Accelerators Spectrometers Detectors Associated Equip (2004) 525:284–88. doi: 10.1016/j.nima.2004.03.128

CrossRef Full Text | Google Scholar

12. Nishio T, Miyatake A, Ogino T, Nakagawa K, Saijo N, Esumi H. The development and clinical use of a beam on-line pet system mounted on a rotating gantry port in proton therapy. Int J Radiat Oncol Biol Phys (2010) 76:277–86. doi: 10.1016/j.ijrobp.2009.05.065

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Knopf A, Lomax A. In vivo proton range verification: a review. Phys Med Biol (2013) 58:R131–160. doi: 10.1088/0031-9155/58/15/R131

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Kraan AC. Range verification methods in particle therapy: Underlying physics and monte carlo modeling. Front Oncol (2015) 5:150. doi: 10.3389/fonc.2015.00150

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Dendooven P, Buitenhuis T, Diblen F, Heeres P, Biegun A, Fiedler F, et al. Short-lived positron emitters in beam-on pet imaging during proton therapy. Phys Med Biol (2015) 60:8923–47. doi: 10.1088/0031-9155/60/23/8923

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Bisogni G, Attili A, Battistoni G, Belcari N, Camarlinghi N, Cerello P, et al. Inside in-beam positron emission tomography system for particle range monitoring in hadrontherapy. J Med Imaging (2016) 4(1)011005–12. doi: 10.1117/1.JMI.4.1.011005

CrossRef Full Text | Google Scholar

17. Rossi S. The status of cnao. Eur Phys J Plus (2011) 126:78–116. doi: 10.1140/epjp/i2011-11078-8

CrossRef Full Text | Google Scholar

18. Traini G, Mattei I, Battistoni G, Giuseppina M, De Simoni M, Dong Y, et al. Review and performance of the dose profiler, a particle therapy treatments online monitor. Physica Med (2019) 65:84–93. doi: 10.1016/j.ejmp.2019.07.010

CrossRef Full Text | Google Scholar

19. Frey K, Unholtz D, Bauer J, Debus J, Min CH, Bortfeld T, et al. Automation and uncertainty analysis of a method for in-vivo range verification in particle therapy. Phys Med Biol (2014) 59:5903–19. doi: 10.1088/0031-9155/59/19/5903

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Fiorina E, Ferrero V, Baroni G, Battistoni G, Belcari N, Camarlinghi N, et al. Detection of interfractional morphological changes in proton therapy: A simulation and in vivo study with the inside in-beam pet. Front Phys (2021) 8:578388. doi: 10.3389/fphy.2020.578388

CrossRef Full Text | Google Scholar

21. Piliero M, Belcari N, Bisogni M, Camarlinghi N, Cerello P, Coli S, et al. First results of the INSIDE in-beam PET scanner for the on-line monitoring of particle therapy treatments. J Instrumentation (2016) 11. doi: 10.1088/1748-0221/11/12/c12011

CrossRef Full Text | Google Scholar

22. Ferrero V, Fiorina E, Morrocchi M, Pennazio F, Baroni G, Battistoni G, et al. Online proton therapy monitoring: Clinical test of a silicon-photodetector-based in-beam pet. Sci Rep (2018) 8. doi: 10.1038/s41598-018-22325-6

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Kostara E, Sportelli G, Belcari N, Camarlinghi N, Cerello P, Guerra AD, et al. Particle beam microstructure reconstruction and coincidence discrimination in PET monitoring for hadron therapy. Phys Med Biol (2019) 64:035001. doi: 10.1088/1361-6560/aafa28

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Piliero MA, Pennazio F, Bisogni MG, Camarlinghi N, Cerello PG, Guerra AD, et al. Full-beam performances of a PET detector with synchrotron therapeutic proton beams. Phys Med Biol (2016) 61:N650–66. doi: 10.1088/0031-9155/61/23/n650

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Hamamatsu, Photonics. Photon is our business. Japan:Vendor's site (2016).

Google Scholar

26. Zagumennyi A, Zavartsev Y, Kutovoi S. Scintillation substances (variants). patent US 7132060. (2004).

Google Scholar

27. Doroud K, Williams M, Zichichi A, Zuyeuski R. Comparative timing measurements of lyso and lfs-3 to achieve the best time resolution for tof-pet. Nucl Instruments Methods Phys Res Section A: Accelerators Spectrometers Detectors Associated Equip (2015) 793:57–61. doi: 10.1016/j.nima.2015.04.056

CrossRef Full Text | Google Scholar

28. Rolo MD, Bugalho R, Gonçalves F, Mazza G, Rivetti A, Silva JC, et al. TOFPET ASIC for PET applications. J Instrumentation (2013) 8:C02050–0. doi: 10.1088/1748-0221/8/02/c02050

CrossRef Full Text | Google Scholar

29. Camarlinghi N, Sportelli G, Battistoni G, Belcari N, Cecchetti M, Cirrone GAP, et al. An in-beam PET system for monitoring ion-beam therapy: test on phantoms using clinical 62 MeV protons. J Instrumentation (2014) 9:C04005–5. doi: 10.1088/1748-0221/9/04/c04005

CrossRef Full Text | Google Scholar

30. Crespo P, Shakirin G, Enghardt W. On the detector arrangement for in-beam pet for hadron therapy monitoring. Phys Med Biol (2006) 51. doi: 10.1088/0031-9155/51/9/002

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Kraan A, Berti AC, Retico A, Baroni G, Battistoni G, Belcari N, et al. Localization of anatomical changes in patients during proton therapy with in-beam pet monitoring: A voxel-based morphometry approach exploiting monte carlo simulations. Med Phys (2022) 49:23–40. doi: 10.1002/mp.15336

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Giordanengo S, Garella M, Marchetto F, Bourhaleb F, Ciocca M, Mirandola A, et al. The cnao dose delivery system for modulated scanning ion beam radiotherapy. Med Phys (2015) 42:263–275. doi: 10.1118/1.4903276

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Zhang J, Lu Y, Sheng Y, Wang W, Hong Z, Sun Y, et al. A comparative study of two in vivo pet verification methods in clinical cases. Front Oncol (2021) 11:617787. doi: 10.3389/fonc.2021.617787

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Knopf A, Parodi K, Paganetti H, Cascio E, Bonab A, Bortfeld T. Quantitative assessment of the physical potential of proton beam range verification with PET/CT. Phys Med Biol (2008) 53:4137–51. doi: 10.1088/0031-9155/53/15/009

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Pennazio F, Battistoni G, Bisogni MG, Camarlinghi N, Ferrari A, Ferrero V, et al. Carbon ions beam therapy monitoring with the INSIDE in-beam PET. Phys Med Biol (2018) 63:145018–145029. doi: 10.1088/1361-6560/aacab8

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Ferrero V, Bisogni MG, Camarlinghi N, Fiorina E, Giraudo G, Morrocchi M, et al. Double-field hadrontherapy treatment monitoring with the inside in-beam pet scanner: Proof of concept. IEEE Trans Radiat Plasma Med Sci (2018) 2:588–93. doi: 10.1109/TRPMS.2018.2870478

CrossRef Full Text | Google Scholar