Localization of anatomical changes in patients during proton therapy with in-beam PET monitoring: A voxel-based morphometry approach exploiting Monte Carlo simulations.

PURPOSE: In-beam positron emission tomography (PET) is one of the modalities that can be used for in vivo noninvasive treatment monitoring in proton therapy. Although PET monitoring has been frequently applied for this purpose, there is still no straightforward method to translate the information obtained from the PET images into easy-to-interpret information for clinical personnel. The purpose of this work is to propose a statistical method for analyzing in-beam PET monitoring images that can be used to locate, quantify, and visualize regions with possible morphological changes occurring over the course of treatment.
METHODS: We selected a patient treated for squamous cell carcinoma (SCC) with proton therapy, to perform multiple Monte Carlo (MC) simulations of the expected PET signal at the start of treatment, and to study how the PET signal may change along the treatment course due to morphological changes. We performed voxel-wise two-tailed statistical tests of the simulated PET images, resembling the voxel-based morphometry (VBM) method commonly used in neuroimaging data analysis, to locate regions with significant morphological changes and to quantify the change.
RESULTS: The VBM resembling method has been successfully applied to the simulated in-beam PET images, despite the fact that such images suffer from image artifacts and limited statistics. Three dimensional probability maps were obtained, that allowed to identify interfractional morphological changes and to visualize them superimposed on the computed tomography (CT) scan. In particular, the characteristic color patterns resulting from the two-tailed statistical tests lend themselves to trigger alarms in case of morphological changes along the course of treatment.
CONCLUSIONS: The statistical method presented in this work is a promising method to apply to PET monitoring data to reveal interfractional morphological changes in patients, occurring over the course of treatment. Based on simulated in-beam PET treatment monitoring images, we showed that with our method it was possible to correctly identify the regions that changed. Moreover we could quantify the changes, and visualize them superimposed on the CT scan. The proposed method can possibly help clinical personnel in the replanning procedure in adaptive proton therapy treatments.

INTRODUCTION

Particle therapy is the treatment of solid tumors by charged particles, mostly protons or carbon ions. Thanks to their characteristic dose deposition profile (Bragg peak), charged particles deposit their dose much more locally than the conventionally used photons, allowing for a better sparing of healthy tissue. However, the steep dose gradients also make it more sensitive to errors in dose delivery. , Among the main error sources are interfractional patient variations, occurring because of setup errors and anatomical modifications in patients during the treatment course. , Such errors can possibly compromise the accuracy and affect therapeutic outcomes of particle therapy treatments in particular.

For some patient groups, like head‐and‐neck patients, typically one or more control computed tomography (CT) scans are made along the course of treatment in order to ensure that the treatment plan is still appropriate. , , , , , , This patient group is particularly prone to interfractional morphological changes mostly from presence of cavities, tumor regression, inflammation, organ motion, and weight loss. , , If the plan is judged no longer appropriate because of dose deteriorations, a new plan is created (“replanning” or “adaptive planning”).

The control CTs are typically made at predefined intervals over the course of treatment, driven by clinical experience for a given patient group and personal assessment of the radiation oncologist. However, given the wide variability in morphological changes between patients throughout radiation therapy, the need for dose optimization for each patient, and the time‐ and resource‐consuming nature of adaptive interventions, it would be preferred to decide in an automated manner when to perform the control CT and eventually adapt the plan. Undoubtedly, cone beam computed tomography (CBCT) can provide valuable information about the anatomy and setup of the patient before treatment, , however it does not give indications about possible dose modifications. Online imaging techniques, that is, techniques that can be used during dose delivery, that allow to raise red flags in case of suspected morphological changes associated with dose modifications are thus highly desired to pursue patient‐tailored treatment optimization.

In this context, various noninvasive in vivo treatment monitoring techniques have been proposed in the past decades that are based on detecting reaction products of nuclear interactions of the beam in the patient, including PET, prompt gammas, and charged fragment detection (see Refs. 12, 13, 14 and references therein). Treatment monitoring with positron emission tomography (PET) imaging is based on the measurement of the activity, induced by nuclear interactions between the beam and the patient tissue. Typically, measured PET images are compared with predicted images, either from Monte Carlo (MC) simulations or from previous measurements, to estimate whether the dose was delivered correctly or whether a new repeat CT is needed.

Among all PET data acquisition modalities, in‐beam PET monitoring is the only option allowing to assess the dose delivery in real time. , , , , In‐beam PET imaging is based on detecting the activity from isotopes like O (half life: 122.24 s), N (half‐life: 2.84 min), C (half life: 19.29 s), and C (half life: 20.334 min). At the National Center for Hadrontherapy (CNAO) in Pavia, currently a dual imaging modality system (an in‐beam dual head PET scanner and a charged fragment tracker) named INSIDE (INovative Solutions for In‐beam DosimEtry in hadrontherapy) is being used. INSIDE is currently in a clinical trial phase to test the clinical feasibility and effectiveness in detecting morphological changes in head‐and‐neck patients.

Despite the broad interest in PET monitoring, there is still no straightforward method to translate the information obtained from the PET images acquired during a treatment session into useful and easy‐to‐interpret information for clinical personnel. In fact, the obtained PET data, arising from nuclear physics processes, cannot be directly related to the delivered dose, determined mostly by electromagnetic processes. Moreover PET data cannot provide direct information about changes in the anatomy of the patient. Various methods to solve this issue have been proposed. The first attempts focused on translating measured activity into dose. , , , Later, several empirical methods were developed to detect variations in particle range by studying the activity profile, including the middle‐point analysis, the shift approach, the most‐likely‐shift approach, and range evaluation based on iso‐activity surface extraction. , , The latter has resulted in the clinical validation of a data‐driven method to identify the pencil beams that pass through regions where morphological changes were identified using the INSIDE PET data.

The goal of the current work is to propose an innovative and fundamentally different approach to apply to in‐beam PET monitoring images, that is inspired by the voxel‐based morphometry (VBM) method widely used in neuroimaging. In contrast to all previously published works, the proposed method does not consider particle range, but gives clear voxel‐by‐voxel statistical information about the probability that a morphological changes occurred. We test the method on a series of simulated in‐beam PET monitoring images based on the INSIDE system, with artificially introduced morphological changes.

MATERIALS AND METHODS

Patient data

This work included data from a 70‐year‐old patient treated in 2018 for squamous cell carcinoma (SCC) with proton therapy at CNAO. A planning CT scan was available with corresponding structure set. The resolution of the CT was 0.98 0.98 2 mm. Figure 1(a) displays a slice of the planning CT scan, with the clinical target volume (CTV) highlighted. The CTV was divided in a low dose region (CTV‐low) receiving 60 Gy in 30 fractions and a high dose region (CTV high) receiving 66 Gy in 33 fractions. Figure 1(b) shows the same CT slice, but now with the region‐of‐interest, including the sinonasal cavity, highlighted. An intensity modulated proton therapy (IMPT) treatment plan based on the planning CT scan was available (Syngo Treatment Planning system version VB10, Siemens), which had three orthogonal irradiation fields. The field at 270° IEC (International Electrotechnical Commission) was selected for this work, delivering 2.1 protons of the total number of 7.0. The relevant characteristics of the treatment plan are summarized in Table 1.

FIGURE 1

(a) Slice of the planning CT of the SCC patient with the CTV‐low in yellow and the CTV‐high in red. (b) Same, but now without target structures and with the region of interest, indicated as orange rectangle, containing the sinonasal cavity. (c) Corresponding slice of the first control CT. The sinonasal cavity is somewhat emptied, as can be seen by inspecting the region inside the orange rectangle. (d) Corresponding slice of the second control CT, where further emptying of the sinonasal cavity is visible

TABLE 1

Relevant treatment plan characteristics of the SCC patient. Np represents the number of particles

Treatment modality	IMPT
Fractions	33
Fields	270° IEC, 0° IEC, 180° IEC
Tumor volume	CTV‐low: 139 mL
	CTV‐high: 40 mL
Dose prescription	CTV‐low: 60 Gy
	CTV‐high: 66 Gy
Dose per fraction	CTV‐low: 2 Gy
	CTV‐high: 2 Gy
Field 1 nr of protons	2.1×10^10 Np
Field 1 minimum energy	66 MeV
Field 1 maximum energy	168 MeV
Average duration of field 1	232 s

Furthermore a first control CT, taken after two fractions, and a second control CT, taken after 22 fractions, were available. A slice of the first and second control CT is displayed in Figure 1(c) and (d), respectively. While the morphological changes seen in the first control CT were minimal and did not require adaptive interventions, the second control CT demonstrated large morphological changes caused by sinonasal cavity emptying, that resulted in serious dose modifications requiring replanning. A new treatment plan was thus created and applied starting from fraction 23.

INSIDE system

The INSIDE system is a bimodal in‐beam treatment monitoring system combining an in‐beam PET system and a secondary charged particle tracker. This work focuses on the PET system, installed since 2019 at CNAO. It consists of two PET heads with a planar geometry and an active area of 10 25 cm each. Each head has 2 5 detection modules, each made of a 16 16 matrix of lutetium fine silicate (LFS) scintillating crystals, with 3.1 3.1 20 mm size and 3.2 mm pitch. The readout electronics consists of a 16 16 array of Hamamatsu Silicon Photomultiplier (SiPM) coupled one‐to‐one to each crystal. The resulting field of view (FOV) is 11.2 22.4 26.4 cm with a voxel size of 1.6 1.6 1.6 mm, that is, 70 140 165 voxels. Both the PET heads and the charged particle tracker are mounted on a submillimetric precise mechanical cart. With the current design, the INSIDE in‐beam PET is compatible with most head and neck treatments, as well as brain treatments. At present the system is undergoing a clinical trial.

Analysis methodology

To determine whether and to which extent the effects of morphological changes occurring over the treatment course are detectable with the INSIDE PET system, a statistical analysis pipeline has been set up. While the current section is limited to the methodology, the technical details of the implementation follow in Section 2.4.

The procedure was inspired by the statistical parametric mapping (SPM) techniques, which are widely used in neuroimaging studies, mainly involving PET and magnetic resonance imaging (MRI) techniques in the study of a large variety of disease conditions. , , SPM maps aiming to highlight possible morphometric differences between two groups of subjects (e.g., subjects affected by a particular pathological condition and control subjects) can be generated according to the VBM approach. In that case, three‐dimensional brain scans of two groups of subjects, once properly spatially aligned, are compared at the single‐voxel level according to a statistical model. Typically, the general linear model is used as a parametric model to describe the voxel signals in terms of a number of predictors. Statistical tests (e.g., two‐sample t‐test) are conducted at each voxel, and a map of the ‐value for each voxel to contain significant distinguishing information between the two cohorts is derived. This map, which has the same dimension of the original images, can be overlaid on the brain scan of a group‐representative subject to highlight the brain areas where significant differences between the two populations have been detected. The SPM map thus contains an easily viewable information, whose interpretation depends on the input data and on the applied preprocessing steps.

In our analysis, rather than comparing two groups of subjects, we examined data from a single subject, that changes over time. In particular, we wanted to examine whether a given observed PET treatment‐monitoring image, affected by morphological changes, acquired during a certain fraction along the treatment course, was compatible within statistical fluctuations with the PET monitoring image acquired at the beginning of the treatment, that is, in absence of morphological changes. To this purpose, we set the problem as a null hypothesis statistical test, where empirical ‐values were obtained with the help of MC simulations. We generated the null reference distribution of the voxel intensities (our test statistics) for the PET image in absence of morphological changes, by running the MC simulation of the treatment based on the planning CT with N random seeds. Moreover, we generated a series of observed1 PET monitoring images, affected by morphological changes. Then, for each voxel the null hypothesis being tested was that the intensity value in that voxel in the observed monitoring PET image was not different from the null reference one, that is, without morphological changes. Rejection of null hypothesis implied that the voxel intensity had significantly changed.

In general, for each voxel in an observed image we can calculate the empirical probability value that its intensity is compatible with the null hypothesis. , , In our problem, the observed voxel intensity could depart from the null distribution in both directions, that is, it could either be more than or less than the average simulated voxel intensity. In fact, a range overshoot, caused for instance by a region that is emptied such as a cavity, could result in less activity than expected in the emptied region, and more activity behind it. A two‐tailed implementation of the null hypothesis statistical test was thus necessary. In this case, the null hypothesis was rejected when the ‐value was sufficiently extreme and thus unlikely to be the result of chance. When providing a significance level , the critical regions would exist on the two tail ends of the 1‐D distribution of voxel intensity values. Let us define as the intensity value in voxel in a PET image, as the average intensity value over all replicates, and as the intensity value in the observed PET image . The empirical ‐values at each voxel could then be computed as  , , : where was the number of replicates and was the number of replicates that produced a test statistic of either if , or if .

The resulting ‐values evaluated for each voxel of the th PET image can be reported as 3‐D voxelized maps. In particular, by setting the conventional significant threshold , two separate masks were obtained, highlighting only the voxels for which the null hypothesis was rejected in favor of the alternative hypothesis, that is, that was either significantly larger or significantly smaller than the reference intensity. These could be overlaid onto the CT of the subject. The two masks indicated the areas, where the activity was less‐then‐expected or more‐than‐expected with respect to the reference treatment situation.

The proposed statistical analysis can thus be summarized as follows:

Obtain the null reference distribution of intensities for each voxel in the PET image by repeating the simulation of the treatment on the planning CT times with different random seeds.

Generate a series of observed PET images with MC, simulating the PET signal from intermediate treatment fractions based on artificially modified CT scans. As mentioned before, we stress that these observed images are generated by MC simulations. Although the method should naturally be extended to observed PET images from experimental data, for the scope of validating the statistical method and assessing its sensitivity in detecting morphological changes, the simulated observed PET images are sufficient.

For each observed PET image , do a voxel‐wise two‐tailed null hypothesis test, yielding for each voxel a ‐value. Significant ‐values, that is, , can be visualized with appropriate colors for the case where (significantly more activity than expected) and (significantly less activity than expected), respectively. These colored 3D maps can be overlaid onto the original CT image, allowing directly to localize the areas where changes in patient's anatomy occurred, that may have caused an over/under dose exposure.

Apply the above method to the real first and second control CT of the patient that were available, and identify the regions that changed.

Technical implementation

MC simulations of null reference distributions of PET voxel intensities

The null reference distributions of PET intensities in each voxel were obtained from MC simulations by running the same patient simulation many times independently. The MC simulation tool developed for INSIDE , was used for this purpose. The tool is based on the FLUKA MC code, suited for particle therapy treatment plan simulations. The following issues were considered in the simulation of a treatment:

Beam line. The last part of the CNAO beam line, including exit window and nozzle were simulated. The pencil beam transverse dimension and shape were reproduced as in measurements for all pencil beam energies. This was done with a previously developed simulation FLUKA tool, that was also used for the dosimetric commissioning and quality assurance of proton and carbon beams at CNAO. , ,

Time structure of the delivered protons, including time, intensity, and duration of each single spot.

Patient geometry. This was obtained by importing the CT scan with its original resolution in a voxel FLUKA geometry. Tissue composition and densities were determined according to the CT Hounsfield Units method by Schneider.

Radiation transport and interactions of protons in the patient. Although it was time consuming, we used the maximum precision defaults of the physics models options of FLUKA, and the energy threshold for charged particle transport was set to 100 keV.

Generation and decay of the 3‐D spatial and time distribution of emitting isotopes. The slowest part of the simulation was the proton transport from the beam line to the production of the emitter inside the patient. Therefore a two‐step strategy was implemented. In the first step, a fraction of the entire statistics of the treatment field, that is, protons, was simulated up to the production of the emitter. A user scoring routine in FLUKA was adopted to extract the spatial coordinates and production time of the emitting isotope (C,O, and so on). This collection of isotopes is referred to as the isotope production map (IPM). Then, in the second step, the IPM was used as positron generator. Each produced emitting isotope was used =100 times to decay into a positron according to the specific decay spectrum of the given isotope and with random direction. The acollinearity angle of the photons was also considered. Accurate photon transport in the patient was simulated down to 100 keV.

Geometry of the INSIDE PET detector, including the individual pixelated LSF crystals. The patient geometry with the PET planes is given in Figure 2.

FIGURE 2

Simulation geometry with beam line (purple), range shifter (green), patient and PET planes, obtained by using the FLAIR graphical interface to FLUKA

Energy depositions of the pairs of 511 keV photons in the detector. Coincidences were selected, that occurred between the start of treatment up to about 6 min. Previously developed software was used to model energy depositions in the detector, and energy and time resolutions were included by Gaussian smearing.

Creation of line‐of‐responses (LORs) and subsequent image reconstruction. This was done with an iterative maximum likelihood estimation maximization (MLEM) procedure.

A 3‐D PET image was obtained, with a FOV of cm, with mm voxels, just like in data (see Section 2.2). Due to the partial angular coverage of the detector, the obtained images suffered from reconstruction artefacts in the direction perpendicular to the two PET detector planes. We applied a median filter of 1.6 mm (1 voxel) in all directions in order to decrease statistical fluctuations.

Since the MC simulations of an entire treatment are computationally heavy, we relied on parallelization. A number of independent replicates was produced, each time changing the random number in the FLUKA simulation. This resulted in 120 IPM and 120 corresponding PET images, making up the null distribution in all voxels. These will be referred to as IPM and PET images, respectively. An illustration of the procedure to create the null distribution for one of the voxels in the PET image is given in Figure 3.

FIGURE 3

Illustration of procedure to create the null distribution in a voxel, starting from 120 PET images, each of which is created with a unique random seed. The color scale used to represent the voxel intensities is also displayed

MC simulation of PET images including morphological changes

We had to estimate the response of the in‐beam PET monitoring images to the day‐by‐day evolution of the emptying of the sinonasal cavity between the planning CT and the control CTs. For this purpose, we created a set of observed PET images representing various situations with morphological changes. Thus, a set of artificial CTs had to be created, that mimicked this evolution. The control CTs were used as qualitative references for mimicking volumes and directions of this process. An in‐house developed image processing tool in Python was used to modify the CT scans as follows.

First, we used rigid coregistration to align the planning CT with the first and the second control CT. The region where emptying of the sinonasal cavity occurred was identified, referred to hereafter as cavity emptying region. Then, in a first layer of voxels in the sinonasal region, the tissue was replaced with air (Hounsfield Unit) in the control CT, and an intermediate CT scan was generated. This was repeated gradually, emptying step‐wise the entire left sinonasal cavity, each time generating a new intermediate CT scan, up to the moment that the cavity was empty. This latter situation was very similar to the situation of the second control CT. The planning CT together with three examples of the intermediate CTs obtained with the above procedure are given in Figure 4.

FIGURE 4

Examples of artificially created CTs with different sizes of cavity tissue volume reductions: 0 mL (the planning CT, (a)), 3.8 mL (b), 7.3 mL (c), and 13.1 mL (d). The region of interest containing the left sinonasal cavity is highlighted with an orange rectangle

In total, six intermediate treatment sessions (= 1 to 6) were simulated. A list of the reduction in tissue volume in the cavity corresponding to each generated intermediate CT is reported in Table 2. Thus, the cavity in the modified CTs evolves from being totally filled (0 mL change in tissue, i.e., no morphological modifications, first table entry) to being totally emptied (13 mL change in tissue, i.e., large morphological changes, last table entry). The sample with will be referred to hereafter as the control sample. Regarding the size of the volume that was changed, we remind that the CT image voxels had a physical volume of 0.002 mL (0.98 0.98 2 mm), so the number of CT voxels that were modified ranged from zero to a several thousands (13.1/0.002) voxels when the cavity was emptied. In terms of PET voxels (1.6 1.6 1.6 mm = 0.004 mL), a complete emptying of the cavity would correspond to about 3 voxels (13.1/0.004).

TABLE 2

List of the generated intermediate CTs, with the corresponding absolute and relative reduction in tissue volume in the cavity

mod. CT j	ΔV(mL)	ΔV(%)
0	0	0
1	1.9	15
2	3.8	29
3	5.7	44
4	7.3	56
5	10.2	78
6	13.1	100

Using the same framework as discussed in Section 2.4.1, for each volume change from Table 2 the IPM and PET images were generated and stored. We will refer to each observed image as IPM and PET images, respectively.

To estimate the statistical error in the evaluation of the volume changes due to the MC simulations, the latter were repeated 10 times on each artificially modified CT.

Besides the PET distributions based on the artificially modified CTs, we also simulated the PET signal for the first (no replanning) and second (replanning) control CT. Rigid registration was applied before we simulated the PET signal. These images will be referred to as “” and “.”

Voxel‐wise null hypothesis test

Using the null reference distributions of PET intensities, we calculated for each voxel in each observed monitoring PET image, , from Table 2 the ‐value as in Equation (1). The significant threshold was , so for either tail it was . For each PET image the two‐tailed statistical test lead to two voxel‐wise ‐values maps identifying voxels whose intensity values were unlikely to occur by chance: The same was done for the and  images.

a blue map indicating voxels whose intensity value was classified as significantly low with respect to the null reference distribution;

a red map indicating voxels whose intensity value was classified as significantly high with respect to the null reference distribution.

Assessment of the capability to identify regions with morphological changes

Qualitative assessment

To validate that the statistical method can correctly identify regions where morphological changes occurred, the 3‐D distribution of obtained ‐values related to each intermediate CT was overlaid onto the planning CT. The red and blue areas of the maps were then compared to the artificially modified CT to confirm that the regions identified as significantly different were indeed related to the regions where morphological changes occurred. The same was done for the and images.

Quantitative assessment

We also checked whether the method can give quantitative information about the morphological changes. Knowing that the emptying of the cavity should result in voxels with an intensity value that is significantly reduced with respect to expectations, we investigated whether the size of the regions identified in the blue map were related to the size of the modified cavity region in the planning CT. Thus, for each of the observed PET monitoring images ( = 1 to 6) plus the control sample ( = 0) we compared the volume that was actually changed (emptied) in the CT scan, , with the size of the volume that was classified as significant in the blue (lack of activity) color map of the PET image analysis, . The latter was determined by summing the volume of the PET image voxels in the blue ‐value map, that fall within the CTV‐low. The value for was determined by summing the volume of the CT voxels, whose densities were changed to air, and that were contained inside the CTV‐low.

Influence of PET image artifacts

We also assessed to what extend the method was affected by the typical image artifacts, that were present in our PET images due to the partial angular coverage of the detector and limited statistics. This was done by repeating the analysis from Section 2.5, but using the 3‐D IPM maps rather than the PET images, that is, an ideal theoretical situation without detector and image reconstruction effects. Regarding the quantitative analysis in Section 2.5.2, we now compared the size of the volume that was classified as significant in the blue color map of the IPM image analysis, , with , calculated as for the PET image analysis. The obtained results were compared with results that were based on the analysis of the PET images.

Uncertainty analysis

In order to verify the robustness of the analysis results, the assessments in Sections 2.5 and 2.6 were repeated for each of the 10 repetitions of the simulations performed for each tissue volume change in Table 2. Each time the mean and standard deviation was evaluated. Furthermore, the impact of various additional error sources was evaluated for the case with the largest volume change of Table 2, 13.1 mL, including the following:

Experimental uncertainties in the cross sections for emission in the MC simulation. These were changed by +10% and −10% and the analysis of Section 2.5 was repeated. This represents an extremely pessimistic situation, where all cross sections are changed in one direction.

Uncertainties in the number of histories in the MC simulations. We varied the fraction of the number of protons in the MC simulation (see Section 2.4.1), which was . The simulation was repeated with and the analysis of Section 2.5 was repeated.

The position of the INSIDE dual head PET chart with respect to the patient. In all relevant directions the chart was moved with 1.6 mm and the analysis of Section 2.5 was repeated.

RESULTS

Visualization of PET and IPM images

We show in Figure 5 an example of a simulated PET image in case of no sinonasal emptying (one of the 120 replicates) overlaid onto the planning CT in the coronal (a), axial (b), and sagittal (c) planes. The corresponding IPM images are also displayed in (d), (e), and (f), respectively. A few observations can be made. First, strong deformation in the axial and coronal planes can be seen in the PET image, which is a result of the partial angular coverage of the detector. The coronal plane, parallel to the PET planes, is the least affected and therefore in the following selected for visualization of the analysis results. Second, the maximum activity generated was not located in the target region, but much in front with respect to the beam entrance, as expected (see for instance Refs. 25, 30).

FIGURE 5

Top figures: PET signal in the coronal (a), axial (b), and sagittal (c) planes. Deformation in the axial and coronal plane can be seen. Bottom figures: IPM signal in the coronal (d), axial (e), and sagittal (f) planes. The beam direction is indicated with an arrow in the coronal and sagittal plane, whereas it is orthogonal to the displayed plane in the axial view

Identification of regions with morphological changes: Qualitative assessment

In Figure 6 we display examples of the simulated PET signal (one of the 10) for volume changes of 0 mL (a), 3.8 mL (b), 7.3 mL (c), and 13.1 mL (d), overlaid onto the planning CT. In fact, this would represent the kind of information that clinical personnel would observe if looking directly at the PET signal. However, the insets at the left bottom of Figure 6(b), (c), and (d) represent how the situation in the CT actually changed. On one hand, by looking carefully at the PET image, we see that the distribution somewhat changes in the region‐of‐interest: for instance, the activity in that region is somewhat enlarged in Figure 6(c) and (d) with respect to Figure 6(a) and (b). Thus, we confirm that the activity signal is sensitive to morphological changes. On the other hand, it is very difficult to interpret the information: It is not feasible to estimate the location and cause of the change in the PET signal. It must be noted that simply subtracting the expected PET images from an observed PET image is not appealing, as shown already in various works. , Since the significant differences in PET signal are often located in regions where the activity is only moderate, statistical fluctuations cause the highlighted regions to appear all over the activated zone, mostly uncorrelated with the region where morphological changes occurred.

FIGURE 6

Evolution of the PET signal with the gradual emptying of the sinonasal cavity, overlaid on the planning CT, for cavity tissue volume reductions of 0 mL (a), 3.8 mL (b), 7.3 mL (c), and 13.1 mL (d). The small insets on the left bottom in (b), (c), and (d) represent the real underlying situation in the artificially modified CT

In Figure 7 we show three different slices of the ‐value maps overlaid on the planning CT, resulting from the statistical analysis of the control sample (Figure 7(a), (e), and (i)), and three intermediate PET images, for 2, 4, and 6, corresponding to emptied volumes of 3.8 mL (Figure 7(b), (f), and (j)), 7.3 mL (Figure 7(c), (g), (k)), and 13.1 mL (Figure 7(d), (h), and (l)), respectively. First of all, we observe that the region of morphological changes was well identified through the blue zones (lack of activity), and corresponded well to the regions that were actually emptied in the CT scan. The latter can be seen from the small insets in the left bottom of the slices, showing the corresponding region in the artificially modified CT scan. Second, red zones are visible, located behind the cavity region with respect to the beam entrance, representing voxels with an excess of activity with respect to expectations. This characteristic pattern, first blue and then red when following the beam direction, is a strong indicator for a significant morphological changes resulting in a beam overshoot, here caused by a replacement of tissue by air, resulting in a lower than expected PET signal in the tissue emptying region, and in a large than expected activity behind it due to the increased beam range. Third, we note in Figure 7 that some noise is visible across the slice, that is not associated to anatomical changes. We come back to this in Section 3.3.

FIGURE 7

‐value maps for emptied volumes of 0.0 mL (a, e, i), 3.8 mL (b, f, j), 7.3 mL (c, g, k), and 13.1 mL (d, h, l). The upper, medium, and lower four plots correspond to a different CT slice. The small insets in the slices show the region that corresponds to the observed ‐value map

Figure 8(a) and (b) show the ‐value map resulting from the statistical analysis performed on the and images corresponding to the first and second control CT, respectively. In Figure 8(a), corresponding to minimal changes, we see that the ‐value map highlights some small red zones behind the emptied region, while no blue voxels can be seen. In Figure 8(b), corresponding to the situation requiring adaptive intervention (see Section 2.1), the significant regions are remarkably well visible through the blue‐red pattern. Thus, the method clearly identifies the regions with morphological changes also for the PET signal simulated on the real CTs.

FIGURE 8

‐value maps overlaid on the planning CT, resulting from the statistical analysis of the PET images corresponding to the first control CT (a) and the second control CT (b). The small insets in the slices show actual situation in the control CTs that gave rise to the maps

Identification of regions with morphological changes: Quantitative assessment

Figure 9 shows the relationship between the volume that was classified as significantly different from the null reference distribution, , and the actually changed volume, . In this figure, the bullet points represent the average values of the observed volume differences obtained over the 10 simulations made for tissue volume change (see Section 2.4.2), and the error bars are the standard deviation. A linear fit is carried out, obtaining the following parameters: and , with an = 0.99. The fitted line is given in the figure, where the line is also given.

FIGURE 9

Observed volume reduction in PET images, , versus actually introduced volume reduction in the intermediate CT scans,

From Figure 9 we see that even small morphological changes in the simulated images can be detected through the VBM‐like analysis we implemented. At the same time we note that most of the points are located above the line . This is an expected effect of our statistical analysis, due to the multiple testing problem. A large number of statistical tests are actually carried out simultaneously, as large as the voxel extent of the CTV‐low (approximately = 35 000 voxels in our case). Once we fix the significance threshold for each test to , for the lower‐than‐expected part of the signals we will obtain up to voxels that are by chance erroneously classified as positive (False Positives or Type I errors). This effect, causing also the noise in Figure 7, can be mitigated by implementing suitable strategies to correct for multiple comparisons which control the proportion of Type I errors.

In addition to the statistical errors related to MC simulations and to the VBM‐like analysis, the impact on the maps and the estimated volume change of the following error sources mentioned in Section 2.7 was evaluated. All associated tissue volumes changes, which were estimated for the point with the maximum volume nominal volume change ( mL), are reported in Table 3. The corresponding maps are given in Figure S1 of the Supplementary Material. Although it is clear that the size of the identified volume is sensitive to errors, the identification of the region was seen to be robust. Moreover, most errors from Table 3 can be mitigated by appropriate MC validation and alignment procedures.

TABLE 3

Possible error sources in quantification of the tissue emptying volume, for the case where mL. The default for this case, , was 12.6 mL

Error source	Size of identified volume (mL)	ΔV(mL)=V_obs−V_default
Cross sections in MC +10%	8.7	−3.9
Cross sections in MC −10%	19.4	7.1
Fraction f=1/60	12.7	0.1
Positioning	13.9	1.4

Figure 10 shows three slices of the ‐value map overlaid onto the planning CT, based on the statistical analysis of the images. The cavity emptying region is very clearly identified and the image hardly contains noise. This is because the images are free of artifacts, and because statistical fluctuations are much smaller than in the PET images displayed in Figure 7.

FIGURE 10

‐value maps resulting from the IPM analysis, for cavity tissue volume reductions of 0.0 mL (a, e, i), 3.8 mL (b, f, j), 7.3 mL (c, g, k), and 13.1 mL (d, h, l). The upper, medium, and lower four plots correspond to a different CT slice. The small insets in the slices show the region in the modified CT that corresponds to the observed ‐value map

In Figure 11 we demonstrate the observed volume change versus the actually introduced volume change, . As in Figure 9, the bullet points represent the average values of the observed volume differences obtained over the 10 simulations made for tissue volume change (see Section 2.4.2), and the error bars are the standard deviation. Again a a linear fit with was carried out, yielding = 0.89 and = 2.15, with an = 0.99. The line is also given. When comparing Figure 11 with Figure 9, the points in the former are seen to be distributed much more regularly than those in the latter figure. Thus, the partial angular coverage of the detector and the resulting small geometrical efficiency result in much more fluctuations in the ‐value maps. Moreover a deviation can be seen from the line : the detected volume is biased toward somewhat too high values. This is caused by False Positives, discussed before in Section 3.3, but since the IPM distributions are less fluctuated than the PET distributions, the effect is small.

FIGURE 11

Observed volume reduction in IPM maps versus actually introduced volume change in the intermediate CT scans

DISCUSSION

Using MC simulations of a patient treatment with proton therapy, we demonstrated that it was possible to produce 3‐D probability maps of morphological changes based on a VBM‐inspired method. There are various issues that deserve discussion.

Application of VBM method to in‐beam PET

We saw from Figure 7 that the quality of the ‐value maps was affected by the limited statistics and the image artifacts in the PET images (Figure 5) when compared to the ideal situation in Figure 10. To investigate the fluctuations, we tested the effect of different sizes of median filters, but we observed that one voxel in all directions was sufficient to mitigate large fluctuations. To decrease fluctuations from poor statistics, caused by the detector geometry and the limited data acquisition time, it would help to combine the different irradiation fields, so more statistics could be collected. This is challenging, since the PET signal acquired during the second and third field delivery would be a mixed contribution of the fields, but strategies are being developed to overcome this problem. The time passing between the delivery of one field and the next (typically of the order of minutes) and combination of signals would have to be carefully taken into account. This is part of future work. Concerning the image artifacts, these cannot be mitigated and are an intrinsic limitation of dual‐head in‐beam PET systems. However, they did not hamper the correct identification of the zones with morphological changes.

Interpretation of color maps

Concerning the interpretation of the color maps, we remind that in our case, the tissue volume change was caused by an emptying of the cavity as a response to the radiation, so it concerned a replacement of water‐like tissue by air. Water and air having a very different activity signal, and given the beam direction, the method resulted in the blue‐red pattern indicating in this case an air cavity followed by a beam overshoot. This was just one example, we believe that the method would also be sensitive to other kinds of anatomical changes, provided that the anatomical change results in a statistically significantly different PET signal. For instance, let us consider the case of tumor shrinkage, which is a common phenomenon in clinical practice. If the tumor tissue shrinks and is replaced by air, the same blue‐red pattern can be expected and the quantification of the volume change would likely be possible. In case the tumor tissue is replaced by water‐like material, some differences may be expected too, since (tumor) cell material contains more carbon than water and thus the PET signal would change. However, maybe there would not be a blue‐red pattern, but only a local blue, or only a local red pattern. Future research and patient data would be needed to investigate the expected sensitivity of the method to reveal different types of anatomical changes. For the moment, the PET images can serve as a warning (“red flag”), that should trigger the attention of the radiation oncologist if changes are found before a control CT is planned.

Finally, we would like to stress that our method should be seen as complementary to already established imaging techniques like CBCT imaging. In fact, CBCT can give valuable information about the positioning of the patient and the tumor anatomy before treatment, but it does not give information about associated dose modifications. For instance, there can be anatomical changes that hardly alter the delivered dose, because the plan is robust enough. In contrast, in‐beam PET can give valuable information about the dose delivery during treatment. Although PET activity and dose are only indirectly correlated, the presence of significant changes in the in‐beam PET signal is undoubtedly a strong indication that the patient dose has changed. We believe that the voxel‐by‐voxel method to identify the location of changes as presented in this work can bring us a step closer to a real dose reconstruction and support further research in the interpretability of PET monitoring images.

Free parameters

Another point of discussion is the presence of free parameters in the ‐value analysis. One free parameter is the significance level for defining significant ‐values (Section 2.3), which was fixed to 0.05, a commonly used value. We repeated the analysis with and saw that the volume and especially the noise decreased somewhat (see Supplementary Material), This is as expected; only the most significant voxels remained.

Another free parameter is the number of replicates used to build the null reference distributions (Section 2.3). We generated replicates to build the null distributions of the IPM and PET image voxels, which turned out to be sufficient. We expect a larger to be somewhat more robust to statistical fluctuations. However, this brings us to a drawback of our method: the large amount of computing resources needed. It took us about 18 h to simulate one replicate, however this aspect can be greatly improved by further parallelization of the simulations or by cloud‐based computing resources. A possible implementation would be a web interface that allows medical personnel to upload all the relevant pieces of information, which would then be processed by a cloud computing service, giving back the results for subsequent analysis and interpretation. Another possibility is the usage of fast MC simulation codes, which are being more and more frequently used in clinical practice.

Other considerations

Regarding the applicability to real data, the viability and effectiveness of the method in real treatments is still to be tested. However, a few promising results were already obtained in this study, that evidenced some clear advantages of our method with respect to previously used methods:

The number of free parameters is small. Apart from and , there are no free parameters in the analysis method. It relies purely on statistical considerations, with a priori fixed significance thresholds.

The method gives voxel‐by‐voxel information. This approach allows not only to detect problems, but also to identify the regions where morphological changes have likely occurred (Figures 7, 8, and 10) and to quantify their volume (Figures 9 and 11), something that is neither immediate nor feasible with other methods. The voxel‐by‐voxel information would thus bring us a step closer to eventually translate the results into dose values that are of real clinical interest, using for instance machine learning techniques and artificially created series of customized CT scans including best guesses of morphological changes.

The method gives specific indications of what actually occurred regarding the morphological change in terms of characteristic blue‐red patterns, as we saw from (Figures 7, 8, and 10). Looking from the beam direction, for an air cavity, we expect a pattern of first a blue zone followed by a red zone. And vice versa, in case the beam passes first a zone of increased density with respect to the unmodified CT, one would observe first a red zone, followed by a blue zone. The presence of such patterns are strong indications of morphological changes, giving qualitative and quantitative indications about the changes ongoing in the patient along the beam path. In case of a gradually occurring morphological change, this signal would become more notable over subsequent treatment fractions. The information can be overlaid onto the CT to investigate the cause. Not only such information can be presented easily to medical personnel, but it could also be used to automate adaptive interventions, for instance, by raising red flags in case the pattern reaches a certain size.

Finally, the method was tested for in‐beam PET simulations, but can be applied also to other imaging modalities, like in‐room or after‐beam PET treatment monitoring, or eventually even charged particle or prompt gamma detection.

CONCLUSIONS

We presented a method to detect interfractional morphological changes in a patient during a proton therapy treatment. To the best of our knowledge, this is the first time that a VBM resembling method, used in neuroimaging, has been applied to in‐beam PET proton therapy monitoring images. Using a set of simulated PET images based on CT scans with and without morphological variations, we tested the method for a patient with a head‐and‐neck tumor where morphological changes occurred, that was treated at CNAO. Although the viability and effectiveness of the method in real treatments are still to be confirmed, a few important results were already obtained. First of all, we confirmed that a voxel‐wise statistical analysis resembling VBM is appropriate to locate morphological changes through in‐beam PET images. This was possible despite the strong image artifacts and limited statistics typically present in in‐beam PET images. Moreover, the analysis results could be clearly visualized through characteristic colored patterns that can be effortless overlaid onto the CT scans. We believe that the presence of such patterns makes the method particularly appropriate for raising red flags, that is, warning signs to trigger the attention of the radiation oncologist. Whether the method eventually allows to automate adaptive planning can only be confirmed in future research with large amounts of patient data. Finally, we demonstrated a satisfactory sensitivity to the detection of cavity volume changes and an approximately linear proportionality with the cavity tissue volume reduction, even down to small changes, which is encouraging in view of applying the method to real data.

CONFLICT OF INTEREST

The authors have no conflicts to disclose.

Supporting Information

Click here for additional data file.