Dosimetric impact of intrafraction motion under abdominal compression during MR-guided SBRT for (Peri-) pancreatic tumors

Twenty patients were included that underwent online adaptive MR-guided SBRT for upper abdominal malignancies between January 2021 and September 2021. Consecutive informed consent was provided for registry in the prospective Multi-OutcoMe EvaluatioN of radiation Therapy Using the MR-linac (MOMENTUM) study (NCT04075305). A summary of the patient characteristics is given in table 1.

All patients were treated with a hypofractionated stereotactic regimen (5 × 8 Gy) on a Elekta Unity (Elekta AB, Stockholm, Sweden) magnetic resonance linear accelerator (MR-Linac), a 7 MV linear accelerator combined with a 1.5 T wide bore MRI scanner. The treatment protocol was as follows: for each fraction, patients were positioned on a custom fitted vacuum cushion, and secured with a custom fitted Neofrakt abdominal corset (Spronken Orthopedie NV, Genk, Belgium) (Heerkens et al 2017). This corset is fastened with Velcro straps. For the first fraction, the tightness of the straps was marked with a pen to ensure reproducibility in abdominal compression for the remaining fractions. The treatment was conducted in an online adaptive fashion (Daamen et al 2021), where every treatment session began with a 3D T₂-weighted (T2w) pre MRI (FOV: 451 × 451 × 220 mm, voxel size: 0.64 × 0.64 × 2.00 mm, TE/TR: 124/1300 ms, FA: 90°) to visualize the anatomy. The gross tumor volume (GTV) and OAR contours were adapted to this scan by the attending radiation oncologist. An isotropic GTV to planning target volume (PTV) margin of 3 mm was used. Our institutional planning objectives and constraints are shown in supplementary data A (available online at stacks.iop.org/PMB/67/185016/mmedia). After contouring, a 9–14 beam step-and-shoot intensity modulated radiotherapy plan was generated, and the treatment delivery was started after approval by the physician. During treatment delivery, interleaved coronal and sagittal T₂/T₁-weighted bFFE cine MRIs were acquired, centered on the GTV. Both slices were acquired at a frequency of 2.8 Hz for the entire beam-on time. A detailed overview of the cine MRI sequence parameters is given in supplementary data B.

Deformable image registration (DIR) with EVolution (de Senneville et al 2016) (α = 0.35, patch size = 5 × 5) was used to warp every dynamic of the cine MRI set to a reference dynamic. The resulting 2D deformation vector fields were used to warp the GTV annotation to the dynamics, resulting in a contour that moved along with the anatomy in 2D + time in both the sagittal and coronal plane. The displacement of its center of mass in the CC, LR (respectively the y and x components of the contour centroid in the coronal plane), and AP (the x component in the sagittal plane) direction was tracked over time, resulting in motion profiles of the tumor centroid with a high temporal resolution of 2.8 Hz (supplementary video (1)). These motion profiles were used for the dose accumulation described in the following section. In order to test the validity of the estimated tumor motion, an additional study was conducted to compare the motion profiles from EVolution to the motion profiles extracted with optical flow (Zachiu et al 2015) and Elastix (Klein et al 2009), two commonly used methods for DIR. A high degree of agreement between EVolution, optical flow and Elastix was found, with the mean differences in center of mass position under one millimeter. Details and results of this study can be found in supplementary data C.

The intrafraction motion of a fraction was characterized by the respiratory amplitude and maximum baseline drift. These were quantified by applying a high-pass and low-pass filter to isolate the motion caused by respiration and baseline drift, respectively. The high-pass filter was an elliptic infinite impulse response filter with a cutoff frequency of 0.05 Hz, and the low-pass filter was a moving average filter with a sliding window of 25 s. The respiratory amplitude was calculated as the difference between the 5th and 95th percentile of the respiratory motion profile, and the maximum baseline drift as the largest signed value in the drift motion profile. Finally, the magnitude of the respiratory amplitude and drift was calculated as the respective vector length (2-norm) over the individual values for CC, AP, and LR, that is: ${\mathrm{resp}}_{2}=\sqrt{{\mathrm{resp}}_{CC}^{2}\,+\,{\mathrm{resp}}_{AP}^{2}\,+\,{\mathrm{resp}}_{LR}^{2}}$ and ${\mathrm{drift}}_{2}=\sqrt{{\mathrm{drift}}_{CC}^{2}\,+\,{\mathrm{drift}}_{AP}^{2}\,+\,{\mathrm{drift}}_{LR}^{2}}.$

The reconstruction of the delivered dose was based on Kontaxis et al (2020). For each fraction, the corresponding log file was saved from the machine. These log files include the current state of the linac, e.g. gantry angle, multi-leaf collimator positions, and beam energy, logged at an interval of 40 ms. Each log file was split into partial log files of 360 ms, containing a single linac state. The splits were made on the absolute acquisition times of the cine MRI dynamics. Then, each partial log file was combined with a pseudo-CT volume, generated from a bulk density assignment of the pre MRI, to calculate a partial dose plan with GPUMCD (Hissoiny et al 2011). All dose calculations were run on a dedicated workstation equipped with an NVIDIA GeForce Titan X GPU. An isocenter shift was applied to each partial dose plan based on the displacement measured in the cine MRI dynamic that corresponds to the same time point as the partial dose plan. In other words, a displacement was only applied to the part of the treatment that was delivered at the same time as when the displacement occurred in the body. Finally, all translated partial dose plans were summed to obtain the delivered dose plan.

The method of dose reconstruction is illustrated in figure 1.

Zoom In Zoom Out Reset image size

Aside from the motion as measured in the imaging data, we simulated intrafraction motion profiles with a wide range of motion amplitudes and drift magnitudes. The rationale for this exercise is that by doing so, we could break down the dosimetric outcome for each clinical motion profile to the individual origins. And reversely, we could assess the pattern of the impact of the periodic respiratory motion separately from the motion drift in anterior, posterior, cranial, caudal, left and right directions. To mathematically model intrafraction motion, we adopted a commonly used equation, which assumes respiratory motion follows an asymmetric cosine pattern in 1D and was first coined by Lujan et al (1999):

$$\begin{eqnarray}&&z\left(t\right)={z}_{0}-r{\cos }^{2n}\left(\displaystyle \frac{\pi t}{\tau }\right).\end{eqnarray} \tag{ 1 }$$

With z₀ the position at exhale, r the respiratory peak-to-peak amplitude, τ the respiratory period and n a parameter to set the asymmetry of the motion profile: a larger value for n leads to a longer time spend in exhale position and steeper motion profiles during a respiratory cycle. The starting phase φ is excluded in equation (1). Visually, we determined that τ = 5 s and n = 3 best fit the general shape of the motion profiles from our measured data. Furthermore, we set z₀ to the average position of the first 30 s, similar to how the mid-ventilation position was defined in the measured motion profiles. Equation (1) was extended from 1D to the three principal directions of motion: cranio-caudal (CC), anterior–posterior (AP) and left–right (LR). An additional component was added to equation (1) to simulate a linear drift over the course of the motion profile, with total magnitude d. The value of r was set as the respiratory amplitude in the CC direction, and the respiratory amplitudes in the AP and LR directions were set to r/2. Because of the relation between r and z₀, that essentially means that the LR and AP motion profiles are half that of the CC motion, outside of the drift term.

Overall, the model for the 3D motion vector ${\bf{z}}\left(t\right)=\left[\begin{array}{ccc}{z}_{CC}\left(t\right) & {z}_{AP}\left(t\right) & {z}_{LR}\left(t\right)\end{array}\right]$ was described with the following equations:

$$\begin{eqnarray}&&{z}_{CC}\left(t\right)={z}_{0}-r{\cos }^{6}\left(\displaystyle \frac{\pi t}{5}\right)+d\displaystyle \frac{t}{T}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{z}_{AP}\left(t\right)=\displaystyle \frac{1}{2}\left({z}_{0}-r{\cos }^{6}\left(\displaystyle \frac{\pi t}{5}\right)\right)+d\displaystyle \frac{t}{T}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{z}_{LR}\left(t\right)=\displaystyle \frac{1}{2}\left({z}_{0}-r{\cos }^{6}\left(\displaystyle \frac{\pi t}{5}\right)\right)+d\displaystyle \frac{t}{T},\end{eqnarray} \tag{ 4 }$$

where T represents the total duration of the simulation motion, which was set to the duration of the measured motion profile of the fraction in question. Equations (2)–(4) were sampled at the frequency of the measured motion (2.8 Hz). In each simulation, drift was implemented in one direction only: if d was set to a nonzero value in the one direction, d was set to zero in the other directions. The model parameters r and d were set to r ∈ {0, 5, ..., 20} mm and d ∈ {−5, 0, 5} mm. Negative values indicate drift in the caudal/anterior/right direction, and positive values indicate drift in the cranial/posterior/left direction. The instance (r, d) = (0, 0) was excluded, since this is equal to the planned dose. In sum, this yielded 34 different motion simulations. For every fraction, the delivered dose was reconstructed for all 34 simulated motion patterns, in the same manner as the dose was reconstructed for the measured motion of that fraction (as described in section 2.3).

The impact of the measured and simulated intrafraction motion was dosimetrically evaluated by comparing the planned to the delivered dose based on key dose-volume histogram (DVH) parameters in the GTV and multiple OARs. The OARs of interest were the duodenum, stomach and small bowel. For the GTV, we evaluated ratios of the D_99%, i.e. the minimum dose to 99% of the volume. For the OARs, we evaluated the ratios of the D_0.5cc , i.e. the minimum dose in 0.5 cm³ of the volume. The ratios were calculated as the delivered dose divided by the planned dose. For example: a fraction with a relative D_99% < 1 means the delivered D_99% of the GTV was lower than the planned D_99%. Only OARs were included for dose evaluation if they were within a 2.5 cm isotropic margin of the GTV, because the dose level outside this margin was too low to be dosimetrically relevant.

The average (SD) relative D_99% of the GTV was 0.98 (0.03) (figure 2). In more than 85% of the fractions (77/89), the delivered D_99% of the GTV was more than 95% of the planned D_99%. For only two fractions, the delivered D_99% of the GTV was lower than 90% of the planned D_99%. Subanalysis revealed that the dose distributions for some patients were more prone to deterioration with the respect to planned target coverage than other, see supplementary data D.

Zoom In Zoom Out Reset image size

The average (SD) relative D_0.5cc of the duodenum, small bowel, and stomach was 0.99 (0.03), 1.00 (0.03), and 0.97 (0.05), respectively. All hot spots of the delivered dose to the OARs remained below 110% of the planned D0.5cc.

Overall, it was observed that fractions with larger respiratory amplitudes generally also inhibit larger drifts. For these fractions, the variability of dosimetric impact on both the GTV and OARs was larger than in fractions where smaller motion had been measured. This can be seen in figure 2, as the relative DVH parameter values fan out with larger respiratory amplitude and drift, while still centered around 1.00.

The GTV D_99% was very robust to respiratory amplitudes up to 10 mm (figure 3). For drifts, asymmetric patterns were observed for the different directions. Drifts in the anterior and right direction induced a decreased tumor coverage, while drifts in the cranial, caudal, posterior and left direction hardly impacted delivered dose.

For the OARs, the dose decline as a function of respiratory amplitude increased was less steep than for the GTV. On average, each OAR was only influenced by a drift in a single direction. This was LR drift for the duodenum, AP drift for the small bowel and CC drift for the stomach. The stomach had a much more pronounced relation to drift, albeit with larger uncertainties. With a respiratory amplitude of under 10 mm, a caudal drift of 5 mm led to the delivered D_0.5cc of the stomach to be over 110% of the planned D_0.5cc on average. Vice versa, a cranial drift of 5 mm caused the delivered D_0.5cc to drop to 90% of the planned D_0.5cc . The 5 mm LR drifts in the duodenum and AP drifts in the small bowel only caused deviations of no more than 5% above or below the planned D_0.5cc .