A filtering approach for PET and PG predictions in a proton treatment planning system

The filtering approach postulates that there is a function f(z) that, convolved with a dose distribution D(z), yields a positron emitter (PE) distribution P(z) (Parodi and Bortfeld 2006):

$$\begin{equation} P(z) = (f{\ast}D)(z) \end{equation} \tag{ 1 }$$

It was shown that such a function exists, with f(z) being denominated as a filter. f(z) can be described in terms of one or more convolutions of a Gaussian and a powerlaw function, known as ${\tilde{Q}_\nu}$ function (Parodi and Bortfeld 2006):

$$\begin{equation} {\tilde{Q}_\nu}(z) = G(z){\ast}P_\nu(z) = \frac{1}{\sqrt{2\pi}}\exp\left(-\frac{z^2}{4}\right)\mathcal{D}_{-\nu}(-z) \end{equation} \tag{ 2 }$$

where G(z), $P_\nu(z)$, $\mathcal{D}_{-\nu}(-z)$ are a Gaussian function, a powerlaw function with shape parameter ν, and a parabolic cylinder function, respectively. The convolution of two or more ${\tilde{Q}_\nu}$ functions with scaling and shifting is still a ${\tilde{Q}_\nu}$ function and its parameters can be calculated based on the parameters of the original ${\tilde{Q}_\nu}$ functions (Parodi and Bortfeld 2006):

$$\begin{equation} \tilde{Q}_{\nu_1}\left(\frac{z-a_1}{\sigma_1}\right)\ast\tilde{Q}_{\nu_2}\left(\frac{z-a_2}{\sigma_2}\right) = \frac{\sqrt{\sigma_1^2 + \sigma_2^2}^{\nu_1+\nu_2-1}}{\sigma_1^{\nu_1-1}\sigma_2^{\nu_2-1}}\tilde{Q}_{\nu_1+\nu_2}\left(\frac{z-a_1-a_2}{\sqrt{\sigma_1^2 + \sigma_2^2}}\right) \end{equation} \tag{ 3 }$$

where a_i and σ_i are the offset and width, respectively, of the ${\tilde{Q}_\nu}$ function i.

The filter could be estimated via deconvolution methods, but these methods are usually ill-posed, which could render the filtering approach useless. However, if f(z) and D(z) are assumed to be ${\tilde{Q}_\nu}$ functions, then P(z) will be a ${\tilde{Q}_\nu}$ function with known parameters as defined in (3). By using the same mathematical formulations, it is possible to use the parameters of P(z) and D(z) to obtain the ${\tilde{Q}_\nu}$ parameters of f(z), relying only on the mathematical properties of ${\tilde{Q}_\nu}$ functions without dealing with a deconvolution problem.

The mathematical formalism of the filtering approach has been extensively used in several studies (Parodi et al 2007a, Attanasi et al 2011, Frey et al 2014a, Schumann et al 2016, Hofmann et al 2019a, Hofmann et al 2019b). Nevertheless, none of the proposals so far has implemented this method in a research version of a commercial TPS and, for the case of PG monitoring, to apply it to arbitrary tissues (only homogeneous cases of water targets were shown, not applicable to patients (Schumann et al 2016)) and energy windows (to support also spectroscopy investigations).

The filter development requires the estimation of the dose and the quantity of interest (PE or PG source distributions), both as laterally-integrated distributions. The filter is developed in water-equivalent space since the dose calculation in the TPS is based on those conditions. However, using water to develop the filters would lead to PE/PG production from oxygen only. Therefore, an arbitrary reference material constituted by the six most abundant elements in tissues (hydrogen, carbon, nitrogen, oxygen, phosphorus and calcium (Attanasi et al 2011)) was created. This reference material is based on one of the materials of the RayStation database and its properties can be seen in table 1. The dose distribution and the corresponding PE/PG source distributions obtained for the reference material will inevitably be shorter with respect to the ones in water due to its higher stopping power. Therefore, the distributions are stretched to their water equivalent thickness using the stopping power ratio. The goal is to have PE/PG distributions for all relevant nuclear interactions in a water-equivalent system (in terms of beam range). The simulations are performed using the Geant4 toolkit (version 10.02.p03) (Agostinelli et al 2003).

The dose and PE/PG distributions were fitted with ${\tilde{Q}_\nu}$ functions to parameterize P(z) and D(z) to then obtain the filter parameters. The fit was initially performed on distributions from a single energy for a first filter approximation. The filter exhibits a small energy dependency and a previous study has tackled this issue by calculating PE distributions for every proton energy of the clinical facility considered (Frey et al 2014a). In the present study, the first filter approximation was used to set the initial conditions for an optimization where the mean squared error from distributions for different energies with the corresponding prediction were considered. The optimization ensures that the filter function will have the best overall accuracy for the entire energy range used in clinical conditions. This two-step process is needed due to the complexity of setting the initial conditions for the optimization algorithm when considering multiple energies. Such an approach allows for the use of the same filter function for any beam model implemented in the TPS, hence not requiring the calculation of the corresponding PE/PG distribution for each energy.

A correction is required when considering arbitrary tissues. Without this correction, any P(z) distribution calculated would yield the distribution in the reference material. This correction scales the yield by a factor corresponding to the amount of target nuclei in an arbitrary compound tissue with respect to the reference material, and is known as local factor g_X(z) or g-factor (Parodi and Bortfeld 2006). X refers to the target nucleus for a given PE production channel or PG emission:

$$\begin{equation} g_X(z) = \left[\frac{w_X(z)\cdot\rho(z)}{w_{X,ref}(z)\cdot\rho_{ref}}\right] \end{equation} \tag{ 4 }$$

w_X(z) is the weight fraction of X in the material at depth z with mass density ρ(z), while w_X,ref is the weight fraction of X in the reference material that has a mass density ρ_ref. This scaling is thoroughly explained in Parodi and Bortfeld (2006).

The filtering approach is particularly suited to be implemented in a TPS since the quantity governing both dose and nuclear interactions is the particle fluence spectrum Φ_E(E). The particle fluence spectrum is an extension of the definition of particle fluence Φ to realistic polyenergetic beams, where Φ is the quotient of the number of particles $\mathop{dN}$ incident on a sphere of cross-sectional area $\mathop{dA}$ (International Atomic Energy Agency 2005). The particle fluence spectrum is then defined as:

$$\begin{equation} \Phi_E(\vec{r},E) = \frac{d\Phi}{dE}(\vec{r},E) \end{equation} \tag{ 5 }$$

The terms proton fluence and proton fluence spectrum can be used instead for the present case. One of the principles of the pencil beam algorithm is the description of the lateral shape of the beam using a function or a set of functions, usually Gaussian functions. This allows to work with laterally-integrated distributions of dose and fluence, which are then multiplied by the lateral shape function:

$$\begin{equation} \Phi_E(\vec{r},E) = F(\vec{r},E)\frac{d\bar{\Phi}}{dE}(z,E) \end{equation} \tag{ 6 }$$

Where $F(\vec{r},E)$ is the lateral shape function and $\frac{d\bar{\Phi}}{dE}(z,E)$ is the laterally-integrated proton fluence spectrum at depth z. $F(\vec{r},E)$ usually considers at least a set of Gaussian functions to describe the beam spread based on the multiple Coulomb and nuclear scattering. $\frac{d\bar{\Phi}}{dE}(z,E)$ can be regarded as a proton energy spectrum, which if considered over the energy interval of interest, yields number of protons. In turn, the dose $D(\vec{r})$ at the position $\vec{r}$ can be calculated as:

$$\begin{array}{*{20}{c}} {D(\vec r)}&{ = \mathop \smallint \nolimits^ {{\rm{\Phi }}_E}(\vec r,E)\frac{1}{{\rho (\vec r)}}\left( {\frac{{dE}}{{dx}}} \right)(\vec r,E)dE}\\ {}&{ = \mathop \smallint \nolimits^ F(\vec r,E)\frac{{d{\rm{\bar \Phi }}}}{{dE}}(z,E)\frac{1}{{\rho (\vec r)}}\left( {\frac{{dE}}{{dx}}} \right)(\vec r,E)dE} \end{array} \tag{ 7 }$$

Where $\frac{1}{\rho(\vec{r})}\left(\frac{dE}{dx}\right)(\vec{r},E)$ is the mass stopping power in the material at $\vec{r}$ with density $\rho(\vec{r})$. Considering the simplest case (i.e. large field, no range-modifying nor blocking devices) for an inhomogeneous volume, the dose in a TPS can be calculated with (Schaffner et al 1999, Soukup et al 2005, RaySearch Laboratories AB 2017):

$$\begin{array}{*{20}{c}} {D(\vec r)}&{ = \mathop \smallint \nolimits^ F(\vec r,E) \cdot \frac{{d{\rm{\bar \Phi }}}}{{dE}}(z,E) \cdot IDD(z,{E_0})dE}\\ {}&{ = \mathop \smallint \nolimits^ F(\vec r,E) \cdot \frac{{d{\rm{\bar \Phi }}}}{{dE}}(z,E) \cdot ID{D_w}({z^ * },{E_0})dE} \end{array} \tag{ 8 }$$

Where IDD(z, E₀) is the laterally-integrated depth-dose profile with mean initial proton energy E₀, while ${IDD_w}(z^*,E_0)$ is the laterally-integrated depth-dose profile in water in the radiological depth z^*. The radiological depth is calculated using the stopping power ratios with respect to water for the voxels in the beam path.

Regarding PE/PG, the production yields $P(\vec{r})$ in a volume depend on the number of target nuclei $n(\vec{r})$ at the position $\vec{r}$ and the cross section σ_X(E) for the nuclear reaction X and energy E. Using the same approach as in (5), (6) and (7), this can be written as:

$$\begin{array}{*{20}{c}} {P(\vec r)}&{ = \mathop \smallint \nolimits^ n(\vec r){\sigma _X}(E){{\rm{\Phi }}_E}(\vec r,E)dE}\\ {}&{ = \mathop \smallint \nolimits^ n(\vec r){\sigma _X}(E)F(\vec r,E)\frac{{d{\rm{\bar \Phi }}}}{{dE}}(z,E)dE} \end{array} \tag{ 9 }$$

The same way that the IDD_w is a representation of the stopping power of protons in a thick water target, a laterally-integrated depth-PE/PG profile (IDP_w) is a representation of σ_X(E) in a thick water target. Therefore, (9) can be rewritten as in (8):

$$\begin{equation} P(\vec{r}) = \int{n(\vec{r}) \cdot F(\vec{r},E)\cdot\frac{d\bar{\Phi}}{dE}(z,E)\cdot{IDP_w}(z^*,E_0)\mathop{dE}} \end{equation} \tag{ 10 }$$

Equation (10) shows that it is possible to use the pencil beam algorithm implementation in a TPS to calculate PE/PG distributions by exploiting the built-in modelling of $F(\vec{r},E)\cdot\frac{d\bar{\Phi}}{dE}(z,E)$. The beam model and the entire dose propagation workflow can be used, provided the number of target nuclei is accounted for and that IDP_w is an input to the TPS. By convolving the IDD_w required for the treatment plan with the filter function, we can obtain the corresponding IDP_w, and then perform the correction of number of target nuclei with the g-factor. Corrections such as fluence reduction due to blocking devices or radiological depths calculation are naturally taken into consideration. Hence, PE and PG sensitivity to heterogeneities and other effects are identical to those of the TPS dose calculation.

For the PET implementation in RayStation, the six most important nuclear reaction channels in proton therapy were included (Parodi et al 2007b, Frey et al 2014a) with corresponding filter functions: 12 C(p,pn)11 C, 16 O(p,3p3n)11 C, 16 O(p,pn)15 O, 16 O(p,2p2n)13 N, 31P(p,pn)30P, 40Ca(p,2pn)38 K. To provide accurate results, experimental cross sections (Bauer et al 2013b) were used in the simulations coupled with a track-length scorer (Parodi et al 2007b), which combines the proton energy fluence with said cross sections. These cross sections were empirically determined by matching corresponding MC predictions to PE yields deduced from dynamic PET scans of phantoms irradiated with protons

For implementation of biological and physical decays, the so-called synchrotron approach (Parodi et al 2008) was followed, relying on time parameters of the beam delivery from available machine beam records. It includes the entire beam time structure in the calculation of activity, thus providing a coherent solution for any kind of monitoring modality (in-beam, in-room or offline) in contrast to frequent simplifications for offline PET monitoring, which only assume continuous irradiation without beam pauses (Bauer et al 2013a). The formalism of Mizuno et al was employed for modelling biological washout (Mizuno et al 2003), decomposed into a slow, medium and fast biological decay based on exponential functions analogue to the physical decay, and depending on different tissues (Mizuno et al 2003, Parodi et al 2008). The response model of the PET scanner was implemented as a Gaussian kernel (Parodi et al 2007a).

PG production is linked to the target nuclei, hence a filter was created for each nuclei in the reference material. Several studies have shown advantages in using different PG energy selection windows (Min et al 2008, Verburg et al 2013, Polf et al 2015, Richter et al 2016, Xie et al 2017), hence the need to consider the energy selection window as a free parameter. Therefore, a single filter for 1 to 10 MeV was developed for each target nuclei, covering the typical PG energy selection. Without further corrections, the convolution of dose with the filters would estimate the number of PG corresponding to an energy selection of 1 to 10 MeV at a given position $\vec{r}$, $P_{1-10\ {\rm MeV}}(\vec{r})$. An LUT approach was then developed to allow for arbitrary windows. The LUT were built based on simulations using homogeneous phantoms made of a single target and an incident proton energy of 300 MeV, i.e. above the clinical range. Number and energy of PG were scored as a function of proton energy while slowing down, resulting in a LUT (see figure 1 for the case of the LUT to estimate the carbon contribution).

Zoom In Zoom Out Reset image size

The LUT-based correction per voxel to estimate the PG yields for the requested energy selection window, $P_{\textrm{min}-\textrm{max}}(\vec{r})$ was then applied, consisting of the ratio of the PG energy integrals of the data in the requested PG energy range [min, max] and the range 1–10 MeV. The PG energy spectrum, $N_{E_{\gamma}}$, to be used for the ratio depends on the proton energy, $E_{proton}(\vec{r})$ at a given position $\vec{r}$ in the patient geometry:

$$\begin{equation} P_{\textrm{min}-\textrm{max}}(\vec{r},pb) = P_{1-10\ {\rm MeV}}(\vec{r},pb)\times\frac{\int_{\textrm{min}}^{\textrm{max}}N_{E_{\gamma}}(E_{proton}(\vec{r}))dE}{\int_{1}^{10}N_{E_{\gamma}}(E_{proton}(\vec{r}))dE} \end{equation} \tag{ 11 }$$

This method is applied for every pencil beam pb in the treatment plan. The PG energy spectrum weighted by $P_{\textrm{min}-\textrm{max}}(\vec{r},pb)$ and in the range [min, max] is scored and associated to position $\vec{r}$. Number of PG and their energy spectrum are obtained for each voxel and pencil beam for each target element.

The suggested energy selection window of 3–6 MeV for the knife-edge slit camera clinical prototype (Richter et al 2016, Xie et al 2017) was used where not stated otherwise. For the case of the knife-edge slit camera, previous studies have shown that neutrons contribute with a significant flat signal (Sterpin et al 2015). However, the current study addresses PG source distributions and neutrons were not considered.

The proton beam delivery system of the Heidelberg Ion Beam Therapy Centre (HIT) was modelled and included in both TPS pencil beam dose engine and Geant4 (Parodi et al 2012, Bauer et al 2014). The reference materials of the TPS database were imported into Geant4 and the same computed tomography (CT) tissue material assignments were used. The density of the materials was calibrated to obtain the same proton range of Geant4 and TPS (Jiang and Paganetti 2004). The calibration factors were obtained by matching MC and TPS depth-dose curves of homogeneous phantoms made of materials derived from the same Hounsfield units in both Geant4 and the TPS.

Four patients treated at HIT with scanned proton beams followed by offline PET (same protocol as in Bauer et al (2013a)) were considered. Patients P1 and P3 received one field from a fixed horizontal beamline for an acoustic neuroma and a sacral chordoma, respectively. Patient P2 was treated for an astrocytoma with two fields from the fixed horizontal beamline, while patient P4 was treated for a cervical paraganglioma with two fields from the gantry beamline. No patient data from PG monitoring were available for this study and the same PET patients were used to assess in-silico the accuracy of the PG filtering approach compared to MC simulations.

The filtering approach was compared to MC estimations in terms of 1) longitudinal shifts and 2) production yield of PE/PG. Shifts of the PE/PG distributions may not exactly correspond to proton range shifts, but they are indicative of possible problems in the treatment delivery. Each case was analysed using the cumulative dose and PE/PG of the entire field, with one-dimensional distributions along the beam penetration direction sampled at each spot coordinate in the plane transverse to the beam direction.

The shifts between MC and TPS distributions were assessed with a fitting method. The fitting method is one of the established methods to estimate shifts from PE/PG distributions (e.g. Gueth et al (2013), Roellinghoff et al (2014), Janssen et al (2014), Pinto et al (2014), Schmid et al (2015)). Following the same approach as in Gueth et al (2013), Roellinghoff et al (2014), Pinto et al (2014), we implemented a spline function, more specifically a non-uniform rational B-spline (NURBS) function. This type of function allows for great flexibility and it is possible to create a function that reproduces very closely a given distribution. In this work, the NURBS was built using the longitudinal distributions predicted with the filtering approach and then used to fit the respective MC one. The longitudinal displacement necessary to obtain the best fit was assumed to be the shift between predicted and MC distributions. The filtering prediction was used to construct the NURBS function due to its inherently smoother distributions, minimizing interpolation issues during the fitting procedure. Positive shifts mean that MC data have longer range. Figure 2 shows some examples of PE and PG distributions from both prediction and MC, and the corresponding NURBS function.

Zoom In Zoom Out Reset image size

The PE/PG production yield differences between MC and TPS profiles were evaluated using the normalized mean absolute error (NMAE):

$$\begin{equation} NMAE = \frac{\sum_{i = 1}^{N}\left|MC_i-TPS_i\right|}{N}\times\frac{1}{\textrm{max}(TPS)}\times100 \end{equation} \tag{ 12 }$$

where MC_i and TPS_i are the i^th bin of the longitudinal MC and TPS-predicted distributions, respectively, and N the number of bins in the distributions. Normalization is done with respect to the maximum value of the considered TPS distribution. The data are compared in the CT grid.

Percentage difference at pixel (i, j), PD(i, j), was used when comparing slices:

$$\begin{equation} PD(i,j) = \frac{\left|MC_{i,j}-TPS_{i,j}\right|}{\textrm{max}(TPS)}\times100 \end{equation} \tag{ 13 }$$

where MC_i,j and TPS_i,j are the value of (i, j) of MC and TPS-predicted distributions, respectively. Only pixels with TPS dose were considered. While the NMAE metric gives differences between production yields, the percentage difference normalized to the maximum of the TPS prediction aims at identifying regions where the effects of differences between TPS and MC can have higher impact. Furthermore, the percentage difference used herein follows a similar concept to dose reporting, in which dose values are relative to the prescribed dose to the PTV (ICRU 2007). For simplicity reasons, the maximum value was used. However, low dose regions, which are intrinsically less highlighted by the percentage difference, may play an important role on the quality assurance of the treatment when considering shifts. Therefore, no low dose threshold was used for the assessment of shifts. All spots in the treatment plan, regardless of their dose, are being considered for the shift analysis.