4D dose calculation for pencil beam scanning proton therapy of pancreatic cancer using repeated 4DMRI datasets

4D magnetic resonance imaging (4DMRI) has a high potential for pancreatic cancer treatments using proton therapy, by providing time-resolved volumetric images with a high soft-tissue contrast without exposing the patient to any additional imaging dose. In this study, we aim to show the feasibility of 4D treatment planning for pencil beam scanning (PBS) proton therapy of pancreatic cancer, based on five repeated 4DMRI datasets and 4D dose calculations (4DDC) for one pancreatic cancer patient. To investigate the dosimetric impacts of organ motion, deformation vector fields were extracted from 4DMRI, which were then used to warp a static CT of the patient, so as to generate synthetic 4DCT (4DCT-MRI). CTV motion amplitudes  <15 mm were observed for this patient. The results from 4DDC show pronounced interplay effects in the CTV with dose homogeneity d5/d95 and dose coverage v95 being 1.14 and 91%, respectively, after a single fraction of the treatment. An averaging effect was further observed when increasing the number of fractions. Motion effects can become less dominant and dose homogeneity d5/d95  =  1.03 and dose coverage v95  =   within the CTV can be achieved after 28 fractions. The observed inter-fractional organ and tumor motion variations underline the importance of 4D imaging before and during PBS proton therapy.

However, these approaches are all based on regular CT imaging, which implies additional dose to the patient and in the case of cone-beam CT, usually poor soft tissue contrast. To alleviate this issue, time-resolved volumetric magnetic resonance imaging (4DMRI) has many benefits for repeated imaging during fractionated treatments, being able to provide precise motion information (von Siebenthal et al 2007) without applying any imaging dose to the patient. Additionally, 4DMRI enables longer imaging sessions for covering and reconstructing organ and tumor motion within several breathing phases. Furthermore, MR imaging for abdominal imaging has a much higher soft tissue contrast than CT, which is beneficial for abdominal imaging.
To investigate the impact of the resulting motion information on pencil beam scanning (PBS) particle therapy, where the beam is scanned through the target volume (Lomax et al 2004), 4D dose calculation (4DDC) is needed. By means of 4DDC, the interplay effect (Phillips et al 1992, Bert et al 2008, Seco et al 2009, Kraus et al 2011, Knopf et al 2011 can be estimated by taking into account both deformations of the target volume and OARs, and timing information of the individual pencil beams.
In this study, we aim to explore the feasibility of such a 4DDC approach for PBS proton therapy of pancreatic cancer. We propose an approach for generating synthetic 4DCT based on repeated 4DMRI, which is then used in combination with 4DDC to quantify the dosimetric impacts of inter-and intra-fractional motion effects (Boye et al 2013, Bernatowicz et al 2013, Zhang et al 2016.
To the best knowledge of the authors, such a 4DDC study for proton therapy of a clinical pancreatic cancer case has not been reported so far.

4DMRI acquisition and reconstruction
Five 4DMR images were acquired for one pancreatic cancer patient, who was undergoing fractionated proton therapy at the Heidelberg Ion-Beam Therapy Center (HIT). Within five weeks, the MR images were taken using a 1.5 T Magnetom Aera MR scanner (Siemens Healthcare, Erlangen, Germany) under free breathing, using a T1weighted radial stack-of-stars sequence with golden angle spacing (Block et al 2014). The patient was positioned in a supine position on the MR table with the arms above the head, as positioned during proton therapy. One 4DMRI measurement lasted for approximately eight minutes. The sequence settings are listed in table 1.
The MR data were offline reconstructed using a modified version of an iterative 4D joint motion-compensated high-dimensional total variation (4D joint MoCo-HDTV) algorithm (Rank et al 2017). It uses the k-spacecenter signal as a self-gating signal to identify 20 overlapping breathing phases for each 4DMRI measurement. In principle, the number of reconstructed breathing phases is a free parameter in the reconstruction algorithm. In this study, it was chosen high enough to reduce intra-bin motion with distinction of inhale and exhale while keeping the reconstruction results stable and the reconstruction time reasonable (Grimm et al 2015, Rank et al 2017. The reconstruction alternates between image reconstruction and motion estimation. Therefore, the reconstruction process is divided into five resolution levels, starting with the lowest resolution. In each resolution level, the under-sampled radial k-space is reconstructed iteratively for each motion phase, followed by the estimation of the vector fields (VFs). VFs estimated in a previous resolution level are incorporated into the iterative reconstruction with increasing weight. The image registration is based on the demons algorithm (Thirion et al 1998, Vercauteren et al 2009. This algorithm uses underlying diffusion models, where object contours of the fixed image are interpreted as a semi-permeable membrane, whereas the moving image's deformation grid diffuses through these interfaces in order to match both images. Figure 1 shows slices of the 3D volume in the end-exhale and end-inhale phase from the reconstructed 4DMRI image of the patient.

Synthetic 4DCT generation
To quantify intra-fractional motion, the end-exhalation phase of the first reconstructed 4DMRI was taken as the reference phase and the VFs to all other breathing phases were calculated by means of the demons based deformable image registration using the open-source software Plastimatch (Sharp et al 2010) for the whole image region. A static 3DCT was chosen as the reference geometry, which was acquired for the purpose of treatment planning with all target volumes (GTV, CTV, ITV, PTV) and OARs contoured. The reference MR image was registered to the reference CT using rigid transformation. The rigid matching was performed manually by shifting the reference MR image to the reference CT until a visually good agreement between kidneys, spleen and liver was achieved. All VFs from 4DMRI were resampled to the CT image resolution (1 × 1 × 3 mm 3 ). A body-mask was applied on the resampled VFs to remove non-physical VF contributions outside the patient arising from the 4DMRI reconstruction artifacts. The quality assurance and validation of the VFs is described in section 2.3. In total, five synthetic 4DCT images (called 4DCT-MRI in the following) were generated by warping the reference 3DCT by means of the resulting motion VFs from each 4DMRI respectively. The workflow is described in figure 2.

QA of registration
Both motion analysis and 4DDC strongly depend on the quality and accuracy of the performed deformable image registration and the resulting vector fields. The quality assurance was performed by the use of three different measures. The tolerances are taken from Brock et al (2017).
• Inverse consistency error (ICE): the end-exhale MR image was registered to the end-inhale MR image and back to end-exhale for each measurement by means of the demons algorithm. The resulting image was subtracted from the initial end-exhale image and the mean intensity difference was calculated (called image ICE in the following). Moreover, the VFs from end-exhale to end-inhale were added to the VFs from end-inhale to end-exhale for each 4DMRI measurement to determine an upper estimate of the deviation of the respective VF pairs (called VF ICE in the following). The recommended tolerance of the VF ICE is 2-3 mm. • Jacobian determinant D: the results should not have negative values, since D 0 indicates non-physical motion, 0 < D < 1 shows volume reduction and D > 1 volume expansion. The Jacobian determinants from VFs between end-exhale and end-inhale MR images were calculated within the patient's abdominal volume and are expected to be ≈1. • Dice similarity coefficients: the liver, spleen and left kidney were manually segmented on both end-exhale and end-inhale phases for each 4DCT-MRI. The end-exhale segmentations were warped by the calculated deformation vector fields between these two phases. The overlap between the end-inhale segmentations and the warped end-exhale segmentations are calculated using the Dice coefficient. Tolerance 0.8-0.9.

CTV motion analysis
Based on the 4DCT-MRIs, VFs and the clinical target volume (CTV) contour of the static 3DCT, a tumor motion analysis was performed by warping the contoured CTV for all five 4DCT-MRIs using the VFs between the reference phase and all other breathing phases. This results in a spatial motion distribution for all voxels inside the CTV, see figure 3.

4D treatment planning
In order to evaluate the dosimetric influence of the organ motion, a 3D single-field uniform dose plan was generated on the reference 3DCT using two dorsal fields with beam angles of +20 • and −30 • , the angles used for PBS therapy of this patient; see figure 4. The 3D dose distribution was optimized to deliver 2 Gy (RBE) in the PTV, which was defined as a 5 mm extended geometric ITV that encapsulated the CTVs of all 4DCT phases. 4D dose distributions in the presence of daily respiratory motion were calculated using each of the five 4DCT-MRI datasets, respectively, and simulating 20 different respiration starting phases, each for a single fraction of dose delivery. From the motion side, breathing periods were extracted from single-slice coronal cine-MRI, acquired after the 4DMRI measurements, from which the diaphragm position was manually determined on 50 consecutive images with a time resolution of 0.22 s between images. The 4DDC considered each of the five synthetic 4DCT-MRI datasets with their respective breathing periods T (9.7 s, 8.2 s, 8.7 s, 8.4 s, 6.7 s).
For the 4D dose calculation, the PSI 4D pencil beam proton dose calculation algorithm was used (Boye et al 2013, Zhang et al 2014. It takes into account both the vector field information of the underlying CT (in this case from the 4DCT-MRI) and the density information from the synthetic CTs in the different breathing phases. The vector field information is utilized for linear interpolation between consecutive respiratory phases. The time stamps t for each pencil beam were extracted from the machine dynamics, providing information on the number of proto ns that will be delivered at each respective time point during irradiation. The 4D dose calculation is based on a deforming calculation grid with a linear interpolation of motion between two breathing phases to account for the lower time resolution between images of consecutive breathing phases, compared to the time between delivered pencil beams. The beam properties of PSI Gantry2 were assumed with energy-dependent beam widths in air of 2.5-5 mm and spacing of 4 mm between two pencil beams (Zhang et al 2015). The doses were finally warped back into the reference geometry of the 3DCT by means of the inverse vector fields and standard linear interpolation.
In detail, the 4D dose evaluation was performed in three subsequent steps, in which the plan quality in terms of homogeneity index d5/d95 and dose coverage index v95 in the CTV were calculated and compared among single fraction 4D plans, ×28-fractionated 4D plans and repeated ×28-fractionated 4D plans.
• Single fraction 4DDC For each 4DCT-MRI with its respective breathing period T, 20 breathing phases were available, which may be the breathing starting phase to the time point t when the first pencil beam is applied to the patient. These 20 phases may be the breathing starting phases for each of the two irradiation fields. Therefore, the dose distribution for all 20 · 20 = 400 possible combinations of breathing starting phases were calculated for the two dorsal fields and one fraction of PBS proton therapy. • Repeated fractionated 4DDC The full 28-fractions treatment scenario was calculated. For each individual fraction, the breathing starting phases for both dorsal fields were randomly selected. Moreover, the underlying 4DCT-MRI were randomly selected. The starting phase dependent single fraction 4D dose distribution of the two fields were accumulated into the reference CT. The robustness of the 4D dose distributions was evaluated by repeating the full 28-fraction treatment scenario 100 times, for which both 4DCT-MRI and breathing starting phases were selected randomly.

QA of registration
The image registration fulfilled the required criteria with respect to ICE, and Jacobian and Dice coefficients, see table 2. Mean absolute intensity differences (image ICE) from 0.31-0.90 (a.u.) were observed between the endexhale MR images and the matched images, registered from end-exhale to end-inhale and back to end-exhale. The calculated mean VF ICE values of 1.28-2.97 mm are within the tolerance limit of 2-3 mm. Mean Jacobian values of 0.95-0.97 were obtained. Also the Dice coefficients of liver, spleen and kidney with values from 0.84-0.94 fulfill the tolerance criteria of 0.8-0.9. Both the criteria and the physical meaning of the quantities are listed in section 2.3.  Figure 3 shows the extracted CTV motion distributions between different breathing phases for all five 4DCT-MRI (m1-m5). The results confirm the superior-inferior (SI) direction to be the main direction of motion. Mean motion amplitudes around 5-8 mm were seen for this patient between inhalation and exhalation, with maximum amplitudes of up to 15 mm and pronounced day-to-day variations of the motion amplitudes. Moreover, CTV motion amplitudes up to 4 mm and 7 mm are present in the left-right (LR) and anterior-posterior (AP) directions, respectively. While the motion patterns for m1-m4 show similar behaviors, the deviations for m5 indicate the necessity of regular 4D imaging due to day-to-day motion variations. The variations in the width of the boxes, quantified by means of the respective standard deviations of the mean motion amplitudes, indicate the amount of CTV deformation during respiration and the day-to-day variations. Standard deviations of the motion distributions of 1.9-2.5 mm were observed for this patient. Figure 4 shows the interplay effects with hot and cold spots within the target volume in the comparison of 3D and 4D dose calculation in one example single fraction PBS proton therapy, including dose-volume histograms (DVHs). CTV, ITV and PTV are contoured in green, yellow and blue, respectively.

Single fraction 4DDC
For the dose homegeneity (d5/d95) within the CTV, median values around 1.13-1.17 can be obtained for single fraction 4DDC for the individual 4DCT-MRI. However, outliers up to 1.28 were also present, as shown in figure 5(a). These values are significantly higher than in the static case (d5/d95 = 1.02), indicating the decreased dose inhomogeneity due to the presence of intra-fraction motion. The dose coverage shows the same trend with mean v95 = 87%-93% and outliers down to below 80% for 4DDC, in contrast to CTV coverage of v95 = 100% for the static case, as listed in table 3.

Repeated 28 fractions 4DDC
For each fraction, both the breathing starting phases and underlying 4DCT-MRI were selected randomly, which resulted in a different dose accumulated on the reference CT. In order to evaluate the robustness of this 4D dose calculation, the 28-fractionation schema were simulated 100 times by randomly selecting 4DCT-  MRI and breathing starting phases. An averaging effect was observed with pronounced improvements in dose homogeneity and tumor coverage, comparable to the impact of the rescanning technique (Knopf et al 2011). Figure 6 shows the repeated 4DDC for a 28-fractionation treatment scheme, where the mean dose homogeneity improved from 1.14 to 1.03 and mean CTV coverage increased from 91% to 100% after fractionation. Some outliers still occur even after 28 fractions, however close to the mean values of d5/d95 and v95. The DVH in figure 7 displays slight differences for both OARs and target volumes comparing 3DDC (dashed lines) to the presented 4DDCs (solid lines) for full treatments with 28 fractions. A zoom into the high-dose region shows the different DVHs for the CTV in the 100 times calculated 28-fraction proton treatments for the CTV, compared to the sharper DVH for the static case. Important dose quantities are also given in table 3 for CTV, PTV and OARs. The maximum dose generally increases due to the interplay effect for both target volumes and OAR. Especially for the PTV, differences in d5/d95 and v95 are present.

Discussion
In this study we showed that 4DMRI based 4DDC is feasible for PBS proton therapy of pancreatic cancer. We observed an averaging out of the interplay effect by fractionation, most prominent in the first six fractions, maintaining a sufficient dose homogeneity and coverage within the CTV after 28 fractions. However, for each individual fraction hot and cold spots occur due to intra-and inter-fractional motion, which may have an impact on the dose response in the tissues.
As this is a feasibility study, we cannot yet conclude on general treatment benefits or constraints, resulting from this study. However, for the analyzed patient, we see that for treatments with a small number of fractions the interplay effect seems to be more pronounced and may have a larger impact on the resulting dose distribution. In particular, for single fraction 4D dose scenarios we observe clear under dosage of the CTV, which would violate CTV coverage criteria. Further mitigation of the interplay effect could be obtained by means of rescanning (Zhang et al 2016(Zhang et al , 2018 or motion reduction by means of gating.  Table 3. Comparison of dose quantities for 3D/4D dose calculations for target volume (CTV, PTV) and OARs (stomach/bowel, spinal chord, kidneys). All quantities are given in % of the prescribed PTV dose after 28 fractions. For the 4D dose quantities, the mean values of 100 simulated 28-fraction treatments are given, except for the maximum dose, where the largest maximum dose of all 100 simulated 28-fractions treatments is given to present an estimation of a worst-case scenario.

Maximum dose [%]
Mean dose (%) We would like to emphasize the significant imaging dose reduction for the patient with repeated 4DMRI measurements in comparison to repeated 4DCT measurements. While 3DCT scans of the abdomen may approach 20 mGy, in 4DCT scans, the radiation exposure may be up to one magnitude higher (Li et al 2005), depending on the CT scanner and settings. For five imaging sessions, as was the case for the reported patient, this could add up to an additional imaging dose of approximately 1 Gy to the whole upper abdomen that is avoided by the 4DMRI approach. In principle, in MR-guided proton therapy, 4DMR images could be acquired before every single fraction. However, due to patient-related limitations and time constraints, this was not performed in this study.
Our approach with five repeated 4DMRI data sets of one patient holds intrinsic uncertainties, as the five measurements may not be representative samples for the breathing variations of the patient. Organ motion could be different than the acquired five different motion patterns during the 28 fractions of proton therapy. This is suggested to be included in the simulation study by additionally randomly varying both the breathing periods and the vector field amplitudes of the underlying 4DCT-MRI for each simulated treatment fraction. However, even with daily 4DMR imaging, motion during treatment may still be different from that observed during the MR measurements. Nevertheless, a basic estimation of the averaging impact of fractionation can be obtained with this sample-based method.
The proposed method of synthetic 4DCT generation comes along with potential uncertainties that arise from image registration. First of all, the registration of the reference MR images and the reference CT held intrinsic uncertainties, since the images were not acquired on the same day with possible resulting differences in patient geometry and positioning. Moreover, also uncertainties in the mono-modal deformable MR image registrations may have a considerable impact on the resulting deformation vector field and the resulting 4DDC, which empha-  sizes the importance of registration quality assurance. Nevertheless, we have seen that abdominal MR image registration fulfills relevant pre-defined QA criteria in our case. It is noted that pancreatic tumor motion is subject to pronounced intra-patient variability. We have observed mean motion amplitudes of 5-8 mm with maximum motion of 15 mm in this study. This is in the range of observations based on 4DCT measurements (mean pancreatic tumor motion amplitudes of 1.34 mm (Huguet et al 2015)/5 mm (Hallmann et al 2012)), cine-MRI (mean amplitudes 6-34 mm (Heerkens et al 2014)/20 mm (Feng et al 2009)) and from tracking systems with external and internal marks (mean amplitudes 7.3-27.3 mm (Knybel et al 2014)). As the magnitude of organ motion has a pronounced influence on the resulting interplay effect, the dosimetric impact will also be patient specific. To fully evaluate the clinical benefit and to conclude on possible treatment strategies for pancreatic cancer patients with PBS proton therapy, more 4DMRI patient data are required and currently being acquired for subsequent 4DDC studies.
In principle, for each new patient, such a 4DDC study could be performed by acquiring multiple 4DMRI data sets with high soft-tissue contrast and without applying any imaging dose to the patient on different days before the treatment. These data could be used to investigate to what extent the interplay effect would wash out for the patient-specific underlying motion pattern and may help to decide whether additional motion mitigation techniques, such as gating, may be required for the patient.

Conclusion
4D dose calculations for PBS proton therapy of pancreatic cancer, based on repeated 4DMRI datasets, are feasible. Interplay effects can become less dominant and a sufficient dose homogeneity and coverage within the CTV can be obtained by fractionation. However, for small numbers of fractions, the interplay impact is more pronounced.