Motion effects in proton treatments of hepatocellular carcinoma—4D robustly optimised pencil beam scanning plans versus double scattering plans

Tina Pfeiler; Dalia Ahmad Khalil; Myriam Ayadi; Christian Bäumer; Oliver Blanck; Mark Chan; Erik Engwall; Dirk Geismar; Sarah Peters; Sandija Plaude; Bernhard Spaan; Beate Timmermann; Jörg Wulff

doi:10.1088/1361-6560/aaecfc

1. Introduction

Proton therapy (PT) takes advantage of the finite range of proton beams to create steep dose gradients in beam direction. Thus, the amount of dose deposited in healthy tissue is reduced but range uncertainties, setup errors and specifically organ motion can cause changes of the intended dose distribution leading to over- and underdosage. As a consequence, moving targets in the thorax or abdomen pose a challenge for proton therapy.

The actual dose distribution in the presence of organ motion depends on the delivery technique and the dynamics of the machine. In active pencil beam scanning (PBS), a narrow proton beam scans the target three-dimensionally. Intensity modulations enable multi-field optimisation (MFO), which allows for a very conformal dose delivery. However, the interplay of organ motion and beam scanning might result in a displacement of the individual beamlets ('spots') with respect to the planned positions, which causes additional dose distortions (Phillips et al 1992, Lambert et al 2005, Bert et al 2008, Seco et al 2009). Passively scattered beams, such as double scattering (DS) beams, are more robust against motion (Slopsema 2012). In DS, the narrow proton beam is spread laterally by a pair of scattering devices. Patient-specific beam shaping devices, in the form of apertures and range compensators, adjust the field laterally and distally to the shape of the target. The modulation of the beam energy is achieved by a fast spinning range modulator wheel, producing a spread-out Bragg peak (SOBP), which is stable over the time scale of organ motion. On the downside, DS deposits unnecessary dose proximal to the target, does not allow for MFO strategies and cannot be applied for all target geometries and sites due to organ at risk (OAR) dose limits.

In the past, moving targets were treated with passive scattering techniques, firstly because the PBS technology was developed later and secondly because DS is more motion-resistant than PBS. In recent years, the first centres started to irradiate tumours in liver and lung with pencil beams and today a majority of new particle therapy centres is equipped with active scanning systems only (Knopf et al 2016). The trend towards PBS requires, even more than DS, adequate motion mitigation techniques. This could be realised by 4D robust optimisation (Graeff 2014, Liu et al 2016, Yu et al 2016) alone or in combination with other motion mitigation techniques. Current clinical 4D robust optimisation tools do not explicitly consider interplay effects but take into account anatomical changes in different respiratory phases. The advantage over commonly used methods, such as rescanning or gating (Bert and Durante 2011, Kubiak 2016), is that the method is less time-consuming and that no additional hardware must be purchased, respectively. Indeed, 4D robust optimisation involves longer computation times, but the time needed for patient preparation and treatment is shorter. The application of 4D robust optimisation is recommended by recently published guidelines of the particle therapy co-operative group (PTCOG) for thoracic malignancies (Chang et al 2017). The recommendation, however, is based on two exploratory studies only (Liu et al 2016, Yu et al 2016).

Technical requirements for the safe implementation of 4D robust optimisation into clinical routine are 4D dynamically accumulated dose (4DDD) analyses for different tumour sites and optimisation algorithms. 4DDD computations depend on the knowledge of the time structure of the beam delivery. Therefore, they require access to the scanning control system calculations or an appropriate, facility-specific beam time model to allow for prospective evaluations. An over-simplified beam time model, assuming for example a constant scanning speed, spot delivery time and energy layer switching time, would cause significant differences between 4DDD computations and experiments (Pfeiler et al 2018).

This article deals with the question whether treatment planning based on 4D robust optimisation can help to reduce motion effects and the dose to normal tissue in PBS PT for hepatocellular carcinoma (HCC). HCC is the most common type of primary liver cancer and the second leading cause of cancer-related death worldwide (Ferlay et al 2015). It was chosen as a model system for the study of moving targets because liver tumours are less affected by range errors than tumours in the thorax where high density gradients and the microstructure of lung tissue lead to additional uncertainties (España and Paganetti 2011, Titt et al 2015, Baumann et al 2017). An experimentally validated interplay effect routine including a facility-specific beam time model served to compute 4DDD distributions in the treatment planning system RayStation (RaySearch Laboratories AB, Stockholm, Sweden) (Pfeiler et al 2018). The clinical target volume (CTV) coverage of 4D robustly optimised plans was compared to conventional 3D non-robustly optimised PBS plans and DS plans. In 3D PBS and DS plans, the planning target volume (PTV) concept was applied to achieve treatment plan robustness. In addition, a conservative single field uniform dose (SFUD) approach was used for 3D non-robustly optimised PBS plans as they are less resistant to disturbances. In contrast, MFO was utilised for 4D robustly optimised PBS plans to exploit the full potential of intensity-modulated proton therapy (IMPT). Further, the influence of 4D robust optimisation on the dose to the normal liver and ribs was investigated.

2. Material and methods

2.1. Patient collective

The current study was based on 4D CT data sets from nine HCC patients including eleven CTVs, which were defined as the gross tumour volume (GTV) for the plan comparison. Seven patients, collected for a previous study (Chan et al 2016), exhibited good tumour contrast due to the selective uptake of lipiodol by tumour cells after trans-arterial chemo-embolization (TACE). The tumour motion was limited by abdominal compression for these patients (patient ID 1–7). For patient 8 and 9, 4D imaging was carried out under free breathing conditions. Patient 8 had two 4D CT scans, one with and one without contrast agent injection, whereas patient 9 had implanted fiducial markers. All 4D CTs were acquired with a Philips Brilliance Big Bore CT scanner (Philips Medical Systems, Cleveland, OH, USA) using retrospective spiral correlated imaging. The CT images were reconstructed from 10% time intervals equally spaced over the respiratory cycle (phase binning). 0% corresponds to the end-inspiration and 50% to the end-expiration phase. The slice thickness of a CT image was 3 mm (ID 1–7), 1 mm (ID 8) and 2 mm (ID 9), respectively.

The CTV volume ranged from $\mathrm{3~cm}^{3}$ to $\mathrm{640~cm}^{3}$ (figure 1(a)). The tumour motion amplitude $d_{\rm iCTV}$ was defined as the average vector length of the deformation vector field (DVF) in-between end-expiration and end-inspiration (Inoue et al 2016):

$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle d_{\rm iCTV}=\frac{\sum\limits_{i}\sqrt{x_{i}^{2}+y_{i}^{2}+z_{i}^{2}}}{N},\label{eq:tumour motion} \nonumber \end{align} \tag{ 1 }$

where the internal CTV (iCTV) represents the union of all CTVs on the different 4D CT phases, x_i, y_i and z_i are the components of the deformation vector of voxel i and N is the total number of voxels within the iCTV. The deformable image registration algorithm used for the DVF computation is presented in section 2.2. $d_{\rm iCTV}$ varied between 0.4 cm and 1.5 cm (figure 1(b)).

**Figure 1.** CTV volume (a) and tumour motion amplitude (b) as defined in equation (1) for the selected patient collective. Patient IDs labelled with the letters 'a' or 'b' refer to different CTVs of one patient. The homogeneity category of a CTV is stated above each bar in plot (a): (I) iCTV encompassed by liver tissue, (II) iCTV extending into other parts of the abdomen and (III) iCTV adjacent to the ribs or lung.
Download figure:
Standard image High-resolution image

In addition, all CTVs were divided into three homogeneity categories depending on the density of the iCTV and its surroundings. Category I comprised all iCTVs entirely encompassed by liver tissue (ID 4 and 6). Category II included iCTVs extending partly into other parts of the abdomen such as the stomach or the kidney (ID 2, 7, 8b and 9a). Category III contained iCTVs adjacent to the ribs or the lung (ID 1, 3, 5, 8a and 9b).

2.2. Treatment planning

Treatment plans were created in a research version of RayStation 6 operated within the clinical installation using a stoichiometric CT calibration curve, the pencil beam dose algorithm and a uniform dose grid of 2 mm. For complex 4D computations, the resolution of the dose grid had to be set to 3 mm or 4 mm for the optimisation process due to the high data volume which could otherwise not be processed. Although the Monte Carlo dose engine is known to be more accurate, the pencil beam dose algorithm was explicitly chosen for the treatment planning study since at time of writing RayStation was not supporting MC calculations for DS in a clinical release. The median dose prescribed to the CTV was 63 $\mathrm{Gy_{RBE}}$ (D_presc) in 15 fractions according to Bush et al (2011). A generic relative biological effectiveness (RBE) value of 1.1 was used to calculate the equivalent biological dose for protons relative to high-energy photons.

Contours for the CTV, the body and the external were created manually on every 4D CT phase for each patient. OARs were only delineated on the planning CT. The PTV was constructed by adding margins for setup and range uncertainties to the iCTV (see sections 2.2.1 and 2.2.2). In this study, the 50% phase served as planning CT since it is better reproducible and less affected by organ motion than other phases (Berbeco et al 2006, Keall et al 2006).

The anatomically constrained deformation algorithm Anaconda was used to perform deformable image registration (DIR) in RayStation (García-Mollá et al 2015, Weistrand and Svensson 2015, Kadoya et al 2016, Zhang et al 2018). Anaconda pursues a hybrid approach by combining intensity and geometric based algorithms. Here, the CTV and body contour were included as controlling structures and the body as focusing region to improve the accuracy of the DIR.

All treatment plans had to fulfil the clinical goal $V_{95}(CTV)=100\%$ , where $V_{95}$ is the percentage volume receiving at least 95% of the prescribed dose. In addition, a robustness evaluation was performed for each plan including a density perturbation, a shift of the isocentre and several respiratory phases. Since it was not feasible to check all possible perturbation scenarios for every respiratory phase, a set of representative scenarios was chosen and evaluated using the scripting interface of RayStation. The permutation of the following cases resulted in eight perturbation evaluations: (1) an isocentre shift of $dx=dy=dz={\pm}2~\mathrm{mm}$ , (2) a density perturbation of ${\pm}5\%$ and (3) a CT at end-expiration and end-inspiration, respectively. For large tumours, $V_{95}(CTV)>$ 99% was demanded in each perturbation scenario, whereas for very small tumours with a size smaller than 20 ${{\rm cm}^{3}}$ $V_{95}(CTV)$ > 98% was required. Exceptions to this rule were made when it compromised the sparing of close-by OARs. Once a plan passed these criteria and was approved from physicists and physicians, it underwent a detailed interplay effect analysis (PBS plans) or a less complex 4D dose analysis (DS plans), respectively (section 2.3).

2.2.1. Pencil beam scanning plans

Two PBS treatment plans were generated for each patient, a 4D robustly optimised plan and a non-robustly optimised, margin-based plan for comparison purposes. The same gantry angles were chosen for both techniques using two treatment fields per plan. The initial proton energy ranged from 100 MeV to 195 MeV. Range shifters were applied to reach shallow-seated tumours, when necessary. The air gap between a range shifter and the patient varied between 2 cm and 5.4 cm. Distances between energy layers were automatically scaled by RayStation depending on the width of the Bragg peak resulting in energy spacings of about 2 MeV–4 MeV. The spot spacing was typically 6 mm–8 mm and was scaled automatically depending on the radial spread in the Bragg peak for an energy layer. The spot size at the isocentre ranges from approximately 6.6 mm for 227 MeV to 13.4 mm for 100 MeV in air (full-width at half-maximum).

In non-robustly optimised PBS plans, uncertainties were taken into account by beam-specific PTVs. The beam specific PTVs comprised the iCTV, a uniform setup margin of 2 mm and a geometric margin accounting for range uncertainties in beam direction (2 mm plus 3.5%). An SFUD approach was applied and dose computation was based on the 50% phase. In the following, non-robustly optimised SFUD plans are referred to as 3D PBS plans.

In 4D robustly optimised treatment plans, uncertainties were incorporated in the optimisation process to increase the robustness of the treatment plan. This is done in RayStation by the minimax optimisation which minimises the objective functions considering the worst case scenario (Fredriksson et al 2011, Fredriksson 2013, RaySearch 2017). Each scenario represents a respiratory phase, setup shift and density perturbation. Here, the robust optimisation settings were chosen to achieve plan robustness against 2 mm setup error, 5% range uncertainty and morphological changes in ten respiratory phases. The applied perturbation parameters originate from the setup uncertainties and range uncertainties clinically used in our centre. For the given energy range, a range uncertainty of 2 mm and 3.5% corresponds to a percentage uncertainty of about 5%. The spots of all fields were simultaneously optimised using MFO. Objectives referring to OARs were only evaluated for the nominal scenario. The 4D robust MFO-IMPT plans are termed 4D PBS plans hereinafter.

2.2.2. Double scattering plans

A three-field concept was necessary to achieve adequate dose conformality for DS plans. The target was defined as the iCTV plus a 2 mm setup margin. In addition, 2 mm and 3.5% range uncertainty were incorporated into the dose computation and the smearing radius of the compensator was set to 3 mm. Like for the 3D PBS plans, the dose was computed on the 50% phase. The SOBP range varied between 8.0 cm and 19.2 cm and the modulation width of the SOBP ranged from 3.5 cm to 13.6 cm.

2.3. 4D dynamic dose accumulation

2.3.1. DS

In DS mode, the IBA range modulator wheel spins with a frequency of about 600 revolutions per minute (Slopsema 2012). Compared to the typical length of a respiratory cycle of 3–5 s, the exact time structure of the beam delivery could be neglected for 4DDD computations assuming a uniform distribution of the dose over all motion phases. Therefore, the 4DDD equals the 4D accumulated dose (4DD) in DS. It is determined by computing the dose for the original plan settings on all respiratory phases with a weighting factor of one over the number of phases, performing deformable image registrations between each respiratory phase and the planning CT and summing up all doses on the planning CT.

2.3.2. PBS interplay effect evaluation

In PBS, the beam delivery is more complex in terms of timing. The time structure has to be taken either from electronic irradiation protocols (log files) or from a dedicated beam time model. The advantage of a model-based computation is that possible fluctuations of the beam delivery parameters can be simulated and considered for treatment planning studies.

A routine tailored to 4DDD computations of irradiations with an IBA Proteus^®Plus proton therapy system (IBA, Louvain-La-Neuve, Belgium) has been developed using the scripting interface of RayStation. In the following, the steps of the interplay routine are briefly summarized.

Based on an empirical beam time model (Pfeiler et al 2018), the delivery time of each spot can be computed and assigned to a respiratory phase of the 4D CT. For this purpose, the respiratory cycle length and the starting phase at the beginning of the irradiation must be known. The resulting fraction dose d_fx,CTi, delivered on respiratory phase i, is calculated and subsequently mapped to the reference phase using deformable image registrations. This procedure is repeated for all respiratory phases. In the end, the deformed doses $d_{fx, CTi}^{Ref}$ and the dose delivered within the reference phase are added up to the fractional 'interplay dose distribution':

$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle d_{fx}=d_{fx,CT1}^{Ref}+d_{fx,CT2}^{Ref}+...+d_{fx,CTp}^{Ref}, \nonumber \end{align} \tag{ 2 }$

where fx denotes the fraction number and p the number of respiratory phases. The absorbed dose delivered over the whole treatment course is given by: $D=\sum\limits_{fx=1}^{N}d_{fx}$ .

In order to distinguish between motion effects in general and interplay effects in particular, 4DD distributions were computed additionally to the 4DDD for PBS plans (analogous to DS computations).

2.3.3. Evaluation settings and criteria

Contrary to DS, the 4DDD depends, in PBS, on the the respiratory cycle length, the start phase at the beginning of the irradiation and the energy layer switching time of the system. Since the timing was assumed to be unknown for the treatment planning study, the outcome of a single fraction was analysed based on 60 simulations per PBS plan with varying starting phases and energy layer switching times according to Pfeiler et al (2018). Due to the high CPU and memory load, the outcome of a whole treatment course was modelled with a reduced set of 5 simulations per fraction, yielding 75 simulations per plan. A constant respiratory period of 4 s was assumed for all evaluations.

Because of their sensitivity to interplay effects the homogeneity index HI

$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle HI=\frac{D_{5}-D_{95}}{D_{presc}} \nonumber \end{align} \tag{ 3 }$

and the percentage over- and underdosage $V_{107/95}$

$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle V_{107/95}=V_{107}+(100-V_{95}) \nonumber \end{align} \tag{ 4 }$

were chosen as evaluation criteria for the plan comparison. Representative for close-by OARs, dose-volume histogram (DVH) metrics of the normal liver tissue (liver volume minus CTV) and ribs were analysed. The risk of radiation-induced liver disease (RILD) was calculated based on the Lyman–Kutcher–Burman (LKB) normal tissue complication probability (NTCP) model (Lyman 1985, Kutcher and Burman 1989, Allen Li et al 2012) as a function of the generalized equivalent uniform dose (gEUD) (Luxton et al 2008) using the scripting interface of RayStation. The used LKB model parameters for RILD were TD₅₀ = 39.8 $\mathrm{Gy_{RBE}}$ (tolerance dose leading to a complication probability of 50% for a uniform irradiation of the entire organ), m = 0.12 (steepness of the NTCP curve) and n = 0.97 (tissue-specific value indicating the dependence on the irradiated partial volume) (Dawson et al 2002). Differences in the fractionation scheme were taken into account by a correction of the DVH according to the linear-quadratic model. In addition to the risk of RILD, the mean dose, $V_{42Gy_{RBE}}$ , $V_{33Gy_{RBE}}$ and the volume receiving less than 15 $\mathrm{Gy_{RBE}}$ were calculated for the normal liver tissue. $V_{60Gy_{RBE}}(rib)$ was chosen to determine the risk of rib fracture after hypofractionated proton beam therapy (Kanemoto et al 2013). The three most exposed ribs per patient were considered for the evaluation. The dose to ribs ranked fourth or higher with respect to $V_{60Gy_{RBE}}(rib)$ was clinically not relevant. The alimentary tract, the heart and the kidney, which were located adjacent to or even within the iCTV for some of the patients, were excluded from the analysis since sufficient CTV coverage in the static plan was of primary importance to examine motion-induced dose distortions. Compromises in the CTV coverage to spare these OARs would have significantly impacted the magnitude of motion effects and thus the outcome of the study. The following limits were used as evaluation criteria for the OARs: $V_{42Gy_{RBE}}(normal\, liver)$ = 33%, $V_{33Gy_{RBE}}(normal\, liver)= 67\%$ , average normal liver dose = 24 $\mathrm{Gy_{RBE}}$ and $V_{60Gy_{RBE}}(rib)= 4.48~{{\rm cm}^{3}}$ (Kanemoto et al 2013, Toramatsu et al 2013, Crane and Koay 2016).

3. Results

3.1. Treatment planning

3.1.1. CTV coverage and robustness evaluation

All plans, except one, passed the robustness evaluation with $V_{95}(CTV)$ values of at least 99% or 98% (very small tumours). Only for CTV 9b, which is located adjacent to the lung, it was not possible to achieve sufficient plan robustness using a 3D PBS approach. Four out of eight perturbation scenarios exhibited a $V_{95}(CTV)$ value below 98%, the lowest 83.2%. Figure 2 shows exemplarily one of the examined perturbation scenarios for DS, 4D PBS and 3D PBS.

**Figure 2.** Sagittal view of the perturbed dose (overlayed as colour wash) for CTV 9b on CT50% for a density perturbation of 5% and an isocentre shift of dx = dy = dz = −2 mm: $V_{95}(CTV)$ is 100.0% for DS, 99.3% for 4D PBS and 86.6% for 3D PBS. (a) DS plan. (b) 4D PBS plan. (c) 3D PBS plan.
Download figure:
Standard image High-resolution image

3.2. 4D dynamically accumulated dose analysis

Three examples for 4DDD distributions are shown in figure 3, each representative of one treatment technique. DS does not exhibit any over- or underdosage of the CTV within the presented CT slice. In contrast, 4D and 3D PBS plans suffer from cold and hot spots, displayed in yellow (90%–95% of D_presc), dark red (107%–110% of D_presc) and pink (>110% of D_presc). Another detail which can be seen on the images is that 4D PBS is superior in sparing OARs as ribs or the normal liver tissue.

Figure 4 visualises motion interplay effects in DVHs using the example of 4D PBS. The first DVH belongs to a large CTV with a small motion amplitude (ID 7) and the second DVH to a small CTV with a large motion amplitude (ID 9a). Both DVHs include the static case, the 4DDD evaluations of a whole treatment course (15 fractions) and 60 4DDD single fraction evaluations forming a DVH band. The small tumour with the large motion amplitude exhibits a much broader CTV DVH band, an indication for a low robustness level. In addition, the difference in the normal liver dose between the static and the dynamic evaluation is slightly larger for this tumour.

**Figure 4.** 4D dynamically accumulated DVHs of the CTV and the normal liver tissue (liver minus CTV) representative for a large CTV volume with a small motion amplitude (a) and a small CTV volume with a large motion amplitude (b) based on a 4D PBS plan: DVH (a) corresponds to patient ID 7 with a CTV volume of about 640 ${{\rm cm}^{3}}$ and a motion amplitude of 0.37 cm. DVH (b) corresponds to patient ID 9a with a CTV volume of 6 ${{\rm cm}^{3}}$ and a motion amplitude of 1.45 cm. The 4DDD outcome of a single fraction ('1fx') was simulated 60 times with varying starting phases and energy layer switching times. In addition, the static case and a 4DDD evaluation of the whole treatment course ('15fx') were included in the plot. (a) DVH for patient ID 7. (b) DVH for patient ID 9a.
Download figure:
Standard image High-resolution image

The following sections give quantitative measures of CTV coverage and OAR sparing.

3.2.1. CTV coverage

Looking at a single fraction, the HI values of the CTV were smallest for DS plans with an average value of 2.0% and an SD of 0.9% (figure 5(a)). For 4D PBS, a significantly higher average HI value of 8.1% ± 1.3% was observed, closely followed by 3D PBS plans with 8.4% ± 1.3%. The maximum HI, corresponding to the worst case out of 60 interplay simulations, reached 7.9% to 18.4% for PBS plans depending on the patient. Figure 5(a) shows that 4D plans exhibited slightly lower HI values than 3D non-robustly optimised plans for category II and III in 8 out of 9 cases, whereas they yielded higher values for category I (patient ID 4 and 6). Simulating the whole treatment course over 15 fractions, the HI decreased to below 3% for 4D and 3D PBS (figure 5(b)), close to the static plan results computed on the reference CT, which were in-between 0.6% and 3.1% for all treatment techniques.

No over- or underdosage could be detected for DS plans in 4DDD analyses. Unlike DS plans, 4D PBS plans exhibited on average a percentage over- and underdosage of 5.9% ± 3.4% in a single fraction. 3D PBS plans were slightly behind with 6.9% ± 5.7%. Again, 4D robust optimised plans yielded smaller average values for category II and III than 3D non-robustly optimised plans in 8 out of 9 cases (figure 6). In the worst case scenario, $V_{107/95}$ was in-between 5.2% and 26.8% for patients 1–8. For patient 9, CTV a and b even exhibited a $V_{107/95}$ of 51.0% to 84.3% in the worst case. However, $V_{107/95}$ approached zero for both, 3D and 4D PBS plans, when averaging over 15 fractions. There was no over- or underdosage in any static plan.

**Figure 6.** Percentage over- and underdosage ( $V_{107/95}$ ) of the CTV for a single fraction based on 4D dynamically accumulated dose computations: the average $V_{107/95}$ is plotted for 4D PBS and 3D PBS plans. Error bars represent the standard deviation of 60 interplay simulations. The patient ID is marked by a 'I' if the iCTV belongs to the homogeneity category I, i.e. it is completely encompassed by liver tissue. This figure does not include DS or the overall outcome for 15 PBS fractions since they did not exhibit any over- or underdosage.
Download figure:
Standard image High-resolution image

The minimum and the maximum dose per fraction, $D_{99.99}(CTV)$ and $D_{0.01}(CTV)$ , remained in-between 95% and 107% of D_presc for all DS plans. For PBS plans, the average maximum dose was below 114% of D_presc for all patients considering 60 scenarios. However, the worst case ranged from 111.6% to 121.8%. The minimum dose per fraction was on average above 90% of D_presc for a single fraction in 4D and 3D PBS. But, depending on the patient, it was 77.4% to 91.3% in the worst case.

Figure 7 visualises the impact of tumour motion, CTV volume and tissue homogeneity on the interplay effect for PBS plans based on $V_{107/95}$ . The highest $V_{107/95}$ values were reached for CTV 9a and CTV 9b which have a small volume, a large motion amplitude and belong to the homogeneity categories II and III. Large volumes with minor motion amplitudes (ID 3 and 7) led to relatively small $V_{107/95}$ values independent of the tissue homogeneity for 3D and 4D PBS plans. Contrary to 4D PBS plans, 3D PBS plans performed better for very small tumours in homogeneous regions (patient 4 and 6).

4DD evaluations of PBS plans revealed no over- or underdosage of the CTV and HI values very similar to those of the static case ranging from 0.8% to 3.9%. The HI of the 4DD distribution showed no trend with regard to the treatment planning technique. It was in 6 out of 11 cases better for 3D than for 4D PBS plans.

3.2.2. OAR sparing

4D PBS plans deposited on average the smallest amount of dose to the normal liver tissue with an average dose of 11.1 $\mathrm{Gy_{RBE}}$ ± 6.8 $\mathrm{Gy_{RBE}}$ (figure 8(a)). The average normal liver dose was about 5.7% ± 6.0% higher for 3D PBS plans and 15.9% ± 10.3% higher for DS plans. In all cases, 4D plans performed either better or at least equivalent to 3D PBS and DS plans. The dose limit of 24 $\mathrm{Gy_{RBE}}$ was satisfied by all plans except for one DS plan which exhibited an average dose of 26.5 $\mathrm{Gy_{RBE}}$ (patient ID 7). 4D PBS plans were also superior to 3D PBS and DS plans regarding $V_{42Gy_{RBE}}$ (figure 8(b)). On average, $V_{42Gy_{RBE}}$ was 15.5% ± 10.2% for DS, 12.7% ± 8.4% for 4D PBS and 13.8% ± 8.9% for 3D PBS. None of the plans exceeded the 33% limit. However, two of the DS plans almost reached this limit (patient ID 3 and 7). The maximum of $V_{33Gy_{RBE}}$ was 37.0%, which is far below the 67% limit. The calculated risk of RILD was non-zero only for two patients. Patient 3 exhibited NTCP values of 19.2% (DS), 1.2% (4D PBS) and 3.4% (3D PBS) and patient 7 of 30.8% (DS), 6.7% (4D PBS) and 7.0% (3D PBS).

**Figure 8.** Average dose (a) and $V_{42Gy_{RBE}}$ (b) for the normal liver tissue based on 4D dynamically accumulated dose computations: each bar illustrates the average value for a certain patient ID and treatment plan (DS, 4D PBS or 3D PBS). PBS interplay simulations exhibited a standard deviation below 0.2 $\mathrm{Gy_{RBE}}$ (average dose) and 0.4% ( $V_{42Gy_{RBE}}$ ), respectively, for all IDs. A dashed line indicates the dose/ DVH constraint in each plot.
Download figure:
Standard image High-resolution image

**Figure 8.** Average dose (a) and $V_{42Gy_{RBE}}$ (b) for the normal liver tissue based on 4D dynamically accumulated dose computations: each bar illustrates the average value for a certain patient ID and treatment plan (DS, 4D PBS or 3D PBS). PBS interplay simulations exhibited a standard deviation below 0.2 $\mathrm{Gy_{RBE}}$ (average dose) and 0.4% ( $V_{42Gy_{RBE}}$ ), respectively, for all IDs. A dashed line indicates the dose/ DVH constraint in each plot.
Download figure:
Standard image High-resolution image

For the examined ribs, $V_{60Gy_{RBE}}$ was on average 2.3 ${{\rm cm}^{3}} \pm 2.2~{{\rm cm}^{3}}$ (DS), 0.3 ${{\rm cm}^{3}}\pm 0.5~{{\rm cm}^{3}}$ (4D PBS) and 1.2 ${{\rm cm}^{3}} \pm 1.4~{{\rm cm}^{3}}$ (3D PBS). In two cases, DS led to irradiated volumes larger than the volume constraint of 4.48 ${{\rm cm}^{3}}$ with 7.3 ${{\rm cm}^{3}}$ and 4.8 ${{\rm cm}^{3}}$ , respectively (table 1). The maximum value for 4D PBS plans was 1.4 ${{\rm cm}^{3}}$ and for 3D PBS plans 4.0 ${{\rm cm}^{3}}$ . The differences in the exposure of close-by ribs for DS, 4D PBS and 3D PBS plans are exemplarily shown in enlarged image sections in figure 3.

Table 1. Comparison of $\mathrm{V_{60Gy_{RBE}}(rib)}$ for DS, 4D PBS and 3D PBS including all patients with ribs in the high-dose region. Only the three most exposed ribs per patient were considered for the evaluation. For 4DDD PBS computations, the standard deviation among 60 simulation with varying starting phases and energy layer switching times is stated. Since the dose was distributed evenly over all respiratory phases in 4DDD DS computations, there is no standard deviation for DS.

	$\mathrm{{V_{60Gy_{RBE}}(rib)}}$ in $\mathrm{{cm^{3}}}$
	ID 1			ID 3
	rib a	rib b	rib c	rib a	rib b	rib c
DS	3.8	1.7	0.2	7.3	4.8	1.2
4D PBS	1.3 ± 0.3	0.2 ± 0.1	0.0 ± 0.0	1.4 ± 0.2	0.4 ± 0.1	0.0 ± 0.0
3D PBS	4.0 ± 0.3	1.5 ± 0.3	0.4 ± 0.1	3.7 ± 0.4	2.2 ± 0.4	0.0 ± 0.0

	$\mathrm{{V_{60Gy_{RBE}}(rib)}}$ in $\mathrm{{cm^{3}}}$

	ID 5	ID 7			ID 8a
	rib a	rib a	rib b	rib c	rib a	rib b

DS	0.8	2.3	1.7	0.1	3.4	0.1
4D PBS	0.1 ± 0.1	0.2 ± 0.0	0.0 ± 0.0	0.0 ± 0.0	0.4 ± 0.12	0.0 ± 0.0
3D PBS	0.7 ± 0.16	0.4 ± 0.1	0.4 ± 0.0	0.1 ± 0.3	1.1 ± 0.4	0.0 ± 0.0

4. Discussion

4.1. Robustness evaluation

DS plans yielded the best CTV coverage in the perturbation scenarios for several reasons: three fields were used instead of two, smearing of compensator was applied, a uniform fluence is delivered and the amount of dose deposited proximal to the target is higher for DS than for PBS. These characteristics make DS plans more robust at the expense of OAR sparing. 4D PBS plans achieved comparatively good CTV coverage for inhomogeneous target regions due to MFO. For example, spots located directly behind a rib (beams eye view) received lower spot weights. In return, the amount of dose deposited by the second field was increased for this region. In 3D PBS plans, neither a third field nor MFO were applied. Density uncertainties were taken into account by adding uniform margins proximal and distal to the target, which cannot adequately compensate for both, range changes in liver and lung tissue, at the same time. As a consequence, it was not possible to create a robust 3D PBS plan for CTV 9b which is directly adjacent to the lung. For all other cases, robustness against the specified density perturbations, shifts of the isocentre and recomputations on the end-inhale phase could be achieved for all three planning techniques.

4.2. 4D dynamically accumulated dose analysis

4.2.1. CTV coverage

On average, 4D robust optimised PBS plans exhibited slightly less dose distortions than non-robustly optimised 3D PBS plans. However, interplay effects were still present and the observed differences of HI and $V_{107/95}$ between 4D and 3D PBS were small in relation to the standard deviation. Therefore, it is difficult to decide whether the slight improvement arises from statistical fluctuations or from an actual decrease of motion effects due to the robust optimisation. One argument against a pure coincidence is that an improvement was observed for the majority of CTVs in inhomogeneous iCTV regions whereas 3D PBS plans performed better for the two homogeneous target regions (IDs 4 and 6). One reason might be that the beam-specific PTV margins in 3D PBS plans were chosen rather conservatively resulting in a larger volume covered by the 95% isodose line for homogeneity category I. The larger safety margins are expected to reduce the magnitude of motion effects at the edge of the CTV. Besides, the SFUD approach increased the robustness if no major density gradients were involved. Contrary, 4D PBS plans had an advantage over 3D PBS plans for homogeneity category II and III since the shift of high density gradients due to motion, e.g. in the vicinity of ribs, was taken into account by the robust optimisation.

On the downside, 4D robust optimisation is associated with an increased computation time. Considering ten respiratory phases and the range and setup uncertainties specified above (resulting in 210 optimisation scenarios), the optimisation usually takes about one hour for target volumes smaller than 500 ${{\rm cm}^{3}}$ using the pencil beam dose algorithm and a resolution of 4 mm. In comparison, it took on average only four minutes to optimise 3D PBS plans with the pencil beam dose algorithm on a 2 mm dose grid. The prolonged computation time might become an issue in clinical routine since usually several iterations are required to find the objectives and constraints which lead to a reasonable treatment plan.

The additional 4DD evaluations, which neglected the time structure of the beam delivery, showed no over- or underdosage and HI values very similar to the static case for both, 4D and 3D PBS plans. This indicates, that the detected over- and underdosage and dose inhomogeneities in 4DDD distributions resulted from interplay and not from other motion effects. The small, positive effect of 4D robust optimisation on interplay, observed for homogeneity categories II and III, could be explained by the fact that 4D robust optimisation takes different types of uncertainties into account. Even though 4D robust optimisation does not consider the time structure of the beam delivery, the robustness against organ motion, setup shifts, and range uncertainties might reduce the magnitude of interplay as a side effect. For example, the consideration of motion-induced range changes may guide the optimizer toward a solution with shallower 4D accumulated dose gradients which improves the 4DDD distribution. Further, MFO was fully exploited in 4D PBS, whereas SFUD was applied in 3D PBS for reasons of precaution because it was not optimised robustly.

A proof of concept study recently published by Engwall et al (2018) investigated the inclusion of uncertainties of the time structure in 4D robust optimisation. The study only considered variations in the respiratory pattern but, in principle, uncertainties in the delivery time could be incorporated as well. For three non-small cell lung cancer patients, Engwall et al demonstrated an improvement in $V_{95}$ , $V_{107}$ and HI for the new optimisation approach compared to the standard 4D robust optimisation. The largest effect was observed for the patient with the highest motion amplitude. The study was limited to evaluations of interplay effects and did not take into account range and setup uncertainties. Thus, a higher number of optimisation scenarios must be considered for clinical treatment plans. The effect is comparable to an increase of safety margins. In order to protect normal tissue and realise acceptable computation times, supplementary motion mitigation techniques such as gating have to be applied in order to reduce the number of scenarios.

In a study on non-small lung cancer, Liu et al reported a significant improvement of the target dose homogeneity and coverage using standard 4D robust optimisation (Liu et al 2016). The greater benefit is probably attributable to the larger density variations in the thorax. The results are not directly comparable since Liu et al used 3D robust optimisation as a reference. Yu et al compared 4D robust IMPT plans with conventional IMPT plans for distal oesophageal cancer (Yu et al 2016). They also observed a notable reduction of dose deviations arising from organ motion and propose to implement robust optimisation for IMPT planning. The outcome of the studies and the reported benefit strongly depends on the patient cohort (target location) and the treatment technique used for the comparison, i.e. if a MFO-IMPT or SFUD approach is applied and if conventional margins or 3D robust optimisation are utilised for the 3D PBS plan.

In this study, PBS 4DDD distributions were additionally compared to DS, the previous standard for moving targets. DS plans are not affected by motion interplay effects due to the almost continuous, fast energy changes. The HI values of the investigated DS plans were even slightly better for 4DDD analyses than for the static case. This is due to the fact that small dose inhomogeneities of the static plan were compensated by other 4D CT phase evaluations in the 4DDD dose. Looking at a single fraction, DS plans clearly outperform 4D and 3D PBS plans. But for the whole treatment course similar results were achieved for all techniques due to the averaging effect of 15 treatment fractions. Reasons for the high homogeneity degree of DS plans under motion are the fast energy changes and the factors listed in section 4.1.

Local overdoses in single fractions were assumed to be non-critical as long as they were located within the CTV. Extreme cases, where the minimum dose in a single fraction was as low as 77% of D_presc, occurred not only in patients with relatively high $V_{107/95}$ values. This is an indication that extreme underdosage does not necessarily compromise the coverage of the whole target and occurs rather locally in the form of small cold spots. Besides, scenarios with such small doses had a low likelihood, i.e. were single outliers among 60 simulations, and the minimum dose was on average still above 90% for a single fraction. In stereotactic body radiotherapy (SBRT), it is common to allow underdosage of the GTV and PTV up to a certain degree even in the static plan (e.g. $V_{90}({\rm GTV})\mathrm{>98\%}$ and $V_{70}({\rm PTV})\mathrm{>98\%}$ in an SBRT plan comparison for liver tumours (Moustakis et al 2018)). In contrast to SBRT, the over- and underdosage here is randomly introduced by the motion interplay which is why it averages out after some fractions. Therefore, both, PBS and DS, produce acceptable 4DDD distributions for the majority of investigated patients. It has been shown by several studies that conventional fractionation schemes can mitigate interplay effects to a large extent (Grassberger 2013, Li 2014, Inoue et al 2016). However, the current trend in radiotherapy towards hypofractionation reduces the averaging effect (Paganetti 2017). In addition, one may note that the impact on the biological effect in a single fraction is not yet sufficiently well known. This becomes an issue for small CTV volumes in combination with large motion amplitudes and inhomogeneous tissue, such as ID 9a and 9b, for which interplay effects are particularly pronounced. Since 4D robust optimisation alone cannot eliminate interplay effects in a single fraction, it should be combined with other motion mitigations techniques, such as respiratory gating or rescanning, for these indications. The difficulty, however, is to decide a priori whether a patient can be safely treated using solely 4D robust optimisation. Since the severity of interplay effects depends on various parameters, such as CTV volume, motion amplitude, tissue heterogeneity or field weighting, a very large patient cohort must be studied in order to establish a correlation between the individual parameters and their impact on interplay effects. So far, it is not possible to specify certain limits for each parameter beyond which supplementary motion mitigation techniques become necessary and individual 4D analyses are required for each patient.

The 3D scatter plot (figure 7) showed that $V_{107/95}$ was comparatively low for large CTVs with small motion amplitudes independent of the homogeneity category. It is not possible to conclude from the available data whether the small magnitude of interplay effects is solely attributable to the low motion amplitude or whether a large CTV volume is also associated with reduced interplay effects. For instance, motion effects might be more severe at the edge of the target where steep dose gradients and tissue inhomogeneities appear. ID 4 and 6 exhibited slightly less over- and underdosage for 3D than for 4D PBS plans. Both cases belong to the homogeneity category I and the corresponding PTV margins were large compared to the CTV volume. Therefore, the coverage of the 3D PBS plans was more conservative compared to the 4D PBS plans at the expense of OAR sparing.

4.2.2. OAR sparing

The question which treatment technique, i.e. DS, 4D PBS or 3D PBS, performs best in sparing OARs cannot be answered in general since it depends on the individual positions of OARs relative to the target and gantry angle. But looking at the normal liver tissue and ribs, PBS plans yielded lower OAR doses than DS, with 4D PBS being better than 3D PBS. OARs proximal to the target from beam's eye view, such as ribs, could be spared well with 4D PBS, whereas large parts of the OAR volume received more than 95% of the prescribed target dose using DS or 3D PBS (see table 1 and figure 3). A reason is that DS deposit unnecessary dose proximal to the target and that the third irradiation field causes a wider overlap of the fields in the entrance region. 3D PBS plans, on the other hand, led to higher OAR doses than 4D PBS plans due to the SFUD approach and uniform PTV margins which do not allow for any flexibility of the optimizer to spare those structures. Therefore, 4D robust optimisation could prevent side effects as rib fractures or RILD for HCC patients.

Although the sharp lateral dose fall-off in DS has the potential to better protect OARs next to the tumour, such as the kidney or the heart, this advantage is often cancelled out by the third irradiation field. On the one hand, the third field lowers the entrance dose of each field but on the other hand more healthy tissue is exposed to irradiation, especially when the third field angle lead to long path lengths through the body.

In PBS, the lateral dose fall-off, and thus the normal tissue sparing capabilities, could be further improved by mounting supplementary apertures on the nozzle as shown for example by Yasui et al (2015), Moteabbed et al (2016) and Bäumer et al (2018).

For two patients with large CTV volumes of almost 500 ${{\rm cm}^{3}}$ and 650 ${{\rm cm}^{3}}$ , DS exceeded the OAR dose constraints and reached unacceptable high probabilities to develop RILD. The risk of RILD for the other tumours having less than 270 ${{\rm cm}^{3}}$ was zero. Hence, the size of the tumour is an important criterion for the choice of the treatment technique and large CTV volumes might not be safely treated with DS. The correlation between the size of a hepatocellular carcinoma and the risk of RILD has been investigated for non-robustly optimised PBS proton plans and intensity-modulated photon therapy (IMRT) plans by Toramatsu et al (2013). The study reported much lower risks of RILD for protons than for photons, particularly when the diameter of the tumour was larger than 6.3 cm (average RILD risk: 94.5% for IMRT and 6.2% for PBS PT). It is expected that the risk would even be lower applying 4D PBS PT.

In this study, dose limits for close-by OARs, such as the stomach or bowel, were ignored in the optimisation if the OAR was located adjacent to or within the iCTV in order to maintain the target coverage. Otherwise, the interplay evaluation would be impacted. In clinic, the PTV (3D PBS and DS) and the robustness settings (4D PBS), respectively, would have to be compromised after some fractions. Another solution could be to place a surgical spacer in-between the liver and the gastrointestinal tract (Komatsu et al 2010).

4.3. Limitations of the treatment planning study

One of the limitations of the study is the small size of the patient cohort, which restricts the validity of the outcome to the underlying conditions. The availability of more suitable HCC patient data sets would enable to divide the patient cohort in subgroups of sufficient size to investigate different influencing factors systematically and get statistically meaningful results. Due to the high number of interplay relevant parameters, not all influencing factors could be studied within this work. For instance, a fixed respiratory cycle length was assumed. The magnitude of interplay effects, however, may differ for other respiratory periods because of the relative timing of dose application and organ motion. Further, another spot size might be more beneficial regarding interplay effects and OAR sparing (Grassberger 2013, Liu et al 2018).

A critical point is the accuracy of the DIR which was checked on a sample basis by calculating the target registration error (TRE) and the dice similarity coefficient (DSC) (Varadhan et al 2013, Brock et al 2017). The comparison with former DIR validation studies (Weistrand and Svensson 2015, Kadoya et al 2016, Zhang et al 2018) showed that Anaconda performed well for the selected HCC patients. The TRE and DSC values were very similar to and in some cases even better than previously reported values. Only when comparing to lung studies, the quality was slightly worse. This was expected due to the higher CT contrast in the case of lung which enables a better DIR quality. The average three-dimensional TRE for the maximum deformation (50% phase to 0% phase) was 0.20 cm ± 0.12 cm standard deviation (SD), with the superior-inferior direction contributing most to the error. The error is smaller when including all respiratory phases in the evaluation. The average DSC value for the liver, considering the deformation between the 50% and 0% phase, was 0.96 ± 0.03, where 1 would indicate a perfect match and zero no overlap at all. The geometric errors induced by DIR can significantly impact 4DDD evaluations and depend on the DIR method as recently shown for three liver cancer patients by Ribeiro et al (2018). All six DIR algorithms under investigation in that study underestimated dose inhomogeneities. The smallest dosimetric errors were observed for a DIR method which included the liver contour as controlling ROI (mean error of $V_{95}(CTV)$ : 0.1% to 2.0%). The magnitude of the error decreased with the number of irradiation fields (one or three) and increased with the motion amplitude (8 mm–21 mm). A mean error of up to 10.6% and an SD of 14.1% were observed for a single field using a DIR method which did not include any controlling or focussing ROI. The results cannot be transferred to the current work in the form of error bars due to differences in the DIR method, the 4DDD computation and the number of fields. But they can help to interpret the 4DDD outcome. In this study, the liver was used as focussing region and multiple fields were applied which should reduce the dosimetric uncertainties. Since the CTV served as controlling ROI, the registration error is expected to be much smaller within the CTV area than in the rest of the liver. As a consequence, the impact of deformation errors on the DVH metrics used for the interplay analysis, namely the HI and $V_{107/95}$ of the CTV, is rated comparatively low. Despite all measures, the 4DDD evaluations must be treated with great caution since the exact extent of residual uncertainties is not known for the given conditions.

Following the example of a previously published study on scanned proton treatments of HCCs (Toramatsu et al 2013), the risk of RILD was calculated based on the LKB-NTCP model parameters determined by Dawson et al for patients treated with 3D conformal photon beams (Dawson et al 2002). However, the use of biological models is prone to uncertainties. Even though the higher biological effectiveness of protons compared to photons was considered in the calculations and DVHs were normalised to correct for differences in the fractionation scheme using the LQ model, the results only represent estimates. Adequate follow-up data for the institution specific patients would be required to assess the validity of a model for deviating conditions (patient cohort, technique and fractionation scheme). Despite all uncertainties, the calculated NTCP values are well-suited for a relative comparison of 4D PBS, 3D PBS and DS plans.

The presented results are only valid for the used machine parameters and each institute has to develop its own beam time model or needs direct access to the scanning control system. Even then, treatment interruptions due to interlocks could prevent 'accurate' predictions. Furthermore, all 4DDD computations rely on the organ motion recorded by the 4D CT scan. Any changes in the breathing pattern during treatment can potentially invalidate them. Also the treatment conditions, as for example the application of contrast agents, could influence the results. For the clinical implementation, this means that the knowledge gained from 4DDD planning studies must be interpreted with great caution and one must bear in mind that individual 4DDD evaluations for patients might suggest a false sense of security. Nevertheless, such studies are crucial to test the suitability of motion mitigation techniques and to identify which types of tumours and target regions benefit most.

5. Conclusion

4D robustly optimised PBS plans did not significantly improve the dose homogeneity of the CTV in proton treatments of HCC patients, but spared OARs to a better extent than non-robustly optimised 3D PBS plans or DS plans. For small, strongly moving tumours, it is advisable to make use of supplementary motion mitigation techniques because both, 4D and 3D PBS plans, yielded rather high $V_{107/95}$ and HI values in a single fraction for those targets compared to DS. DS did not cause any over- or underdosage and exhibited homogeneity indices smaller than 4% in the 4DDD distribution. Therefore, it is the method of choice in terms of target coverage for moving tumours. However, DS is often not a clinical option, particularly for large CTV volumes, because of the less conformal dose distribution which causes OAR dose limits to be exceeded. In addition, DS is not available in most of the recently built particle therapy centres. Considering that 4D PBS plans yielded (1) homogeneous 4DDD distributions looking at the sum dose of all fractions, (2) on average a minimum dose above 90% within a single fraction and (3) the lowest OAR doses, 4D robust optimisation represents a good alternative to DS for most HCC patients. Contrary to commonly used motion mitigation techniques, such as gating or rescanning, 4D robust optimisation does not increase the treatment time and does not need any additional hardware unlike gating. Currently, individual 4DDD analyses are needed to decide whether a patient would benefit from 4D robust optimisation and whether supplementary motion mitigation technique is needed. The categorisation of patients according to tumour motion, size and tissue homogeneity into different treatment groups could replace these time-consuming evaluations but requires the examination of a huge patient cohort. Improvements in computing power and software could facilitate such studies in the next few years.

Acknowledgments

We would like to thank Professor Dr Jürgen Dunst from the University Clinic Schleswig-Holstein for his support. Thanks also goes to the local IT team at WPE.

Motion effects in proton treatments of hepatocellular carcinoma—4D robustly optimised pencil beam scanning plans versus double scattering plans

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction