Robustness of sweeping-window arc therapy treatment sequences against intrafractional tumor motion.

PURPOSE
Due to the potentially periodic collimator dynamic in volumetric modulated arc therapy (VMAT) dose deliveries with the sweeping-window arc therapy (SWAT) technique, additional manifestations of dosimetric deviations in the presence of intrafractional motion may occur. With a fast multileaf collimator (MLC), and a flattening filter free dose delivery, treatment times close to 60 s per fraction are clinical reality. For these treatment sequences, the human breathing period can be close to the collimator sweeping period. Compared to a random arrangement of the segments, this will cause a further degradation of the dose homogeneity.


METHODS
Fifty VMAT sequences of potentially moving target volumes were delivered on a two dimensional ionization chamber array. In order to detect interplay effects along all three coordinate axes, time resolved measurements were performed twice--with the detector aligned in vertical (V) or horizontal (H) orientation. All dose matrices were then moved within a simulation software by a time-dependent motion vector. The minimum relative equivalent uniform dose EUDr,m for all breathing starting phases was determined for each amplitude and period. Furthermore, an estimation of periods with minimum EUD was performed. Additionally, LINAC logfiles were recorded during plan delivery. The MLC, jaw, gantry angle, and monitor unit settings were continuously saved and used to calculate the correlation coefficient between the target motion and the dose weighed collimator motion component for each direction (CC, LR, AP) separately.


RESULTS
The resulting EUDr,m were EUDr,m(CCV) = (98.3 ± 0.6)%, EUDr,m(CCH) = (98.6 ± 0.5)%, EUDr,m(APV) = (97.7 ± 0.9)%, and EUDr,m(LRH) = (97.8 ± 0.9)%. The overall minimum relative EUD observed for 360(∘) arc midventilation treatments was 94.6%. The treatment plan with the shortest period and a minimum relative EUD of less than 97% was found at T = 6.1 s. For a partial 120(∘) arc, an EUDr,m = 92.0% was found. In all cases, a correlation coefficient above 0.5 corresponded to a minimum in EUD.


CONCLUSIONS
With the advent of fast VMAT delivery techniques, nonrobust treatment sequences for human breathing patterns can be generated. These sequences are characterized by a large correlation coefficient between a target motion component and the corresponding collimator dynamic. By iteratively decreasing the maximum allowed dose rate, a low correlation coefficient and consequentially a robust treatment sequence are ensured.


INTRODUCTION
Delivering dose to a moving target can potentially degrade the dose homogeneity in the target volume. 1,2 Various strategies to minimize the relative motion between the treatment beam and the target volume have been described for volumetric modulated arc therapy (VMAT): 3 The first approach is to treat the target volume in a static situation by using either respiratory gated radiation therapy 4 or active breath-hold techniques. 5,6 Since these techniques make use of only certain breathing phases, the treatment time is typically increased. The second approach is to perform a real-time adaptation of the multileaf collimator (MLC) aperture depending on the position of the target volume in beams-eye-view, known as dynamic multileaf collimator (dMLC) tracking. [7][8][9] While this approach has great potential to assure a fast and accurate treatment delivery, real-time information of the target position may not be available with sufficient accuracy. If neither gating nor tracking can be applied, treatments have to be performed in free breathing. For these cases, the resulting dosimetric effects have been theoretically predicted [10][11][12][13] and validated in experiments [14][15][16][17][18] and simulations. [19][20][21][22][23] Two manifestations of dosimetric errors have to be considered: dose blurring and interplay effects. A dose delivery with a static MLC while the target is moving results in a blurred (smeared out) dose distribution. A degradation of the resulting dose distribution in the target volume occurs adjacent to the gradients of the static dose distribution by the extent of the motion amplitude. In dose blurring, the expectation value of the dose in every voxel can be estimated by a convolution of the dose distribution without motion and the probability density function (pdf) of the target motion. Nonrigid organ motion 24 and the changes in radiological depth due to target motion are not considered in this model (static dose cloud approximation 12 ). According to the central limit theorem, every pdf converges to a Gaussian distribution if a sufficient number of fractions are delivered. Treatment planning in the presence of a blurred dose distribution is presented in AAPM report 91. 25 The planning target volume (PTV) is obtained by enlarging the clinical target volume (CTV) with an appropriate CTV-PTV margin. 26 For a dynamic dose delivery, interplay or interference effects between the linear accelerator's (LINAC's) MLC and the target motion may be present. In contrast to a pure amplitude determined dose blurring, dosimetric errors can occur at all segment dose gradients within the target volume. Therefore, the pdf of a voxel no longer equals the pdf of the entire motion pattern. It has to be modified to only include the time intervals where a beam is actually contributing to the dose in the voxel. The structure of the temporo-spatial delivery pattern will have the highest impact when only a small number of fractions are delivered in combination with a fast, highly modulated treatment sequence. For these delivery types, a convergence of the voxel pdf to a Gaussian distribution may not be warranted.
Sweeping-window arc therapy (SWAT) relies on a change of the motion direction of the individual MLC leaves after a constant degree of gantry angle Θ (increment factor ∆Θ). VMAT dose deliveries with the SWAT technique 27 are therefore characterized by a potentially periodic motion of the collimator. 28 This feature allows, on the one hand, for a good starting condition during the treatment plan optimization and, on the other hand, it minimizes the amount of MLC leaf travel. This reduces the treatment time even further, compared to conventional VMAT. In SWAT, different collimator angles may be used. While Cameron 27 used 90 • collimator angle and therefore an MLC leaf motion parallel to the axis of gantry rotation, Ulrich et al. 29 used 0 • and Otto 3 45 • collimator angle. Intensity modulation is predominantly performed through dose rate and leaf speed variation, 3 which allows for a close to constant angular gantry rotation velocity.
Since in VMAT typically one arc is used, the phase relation cannot change on a beam-to-beam basis, as it is the case in step-and-shoot or dMLC IMRT treatments. With a fast MLC, and a flattening filter free (FFF) dose delivery, treatment times close to 60 s per fraction are clinical reality, 30 even for hypofractionated treatments with more than 5 Gy fraction dose. Since these treatment times contain only few breathing cycles, an averaging of the breathing pattern gets less likely and the human breathing period may get close to the collimator sweeping period. In analogy to the phenomenon of interference of temporally coherent waves, this may cause additional dosimetric deviations.
These facts require a further analysis of dosimetric robustness beyond the observations presented for IMRT. 11,12,17 To evaluate motion induced dosimetric artifacts specific to SWAT, three aspects shall be addressed: First, a method to quantify dosimetric artifacts in measured dose distributions in the presence of motion is presented; second, a theoretical prediction of breathing periods with minimum dose coverage is made; and third, a sensitive parameter that can be used to mitigate dosimetric deviations due to intrafractional target motion is identified and applied.

2.A. Patient and treatment sequence collective
An Elekta (Stockholm, Sweden) Synergy LINAC with MLCi2 collimator and an Elekta VersaHD LINAC with Agility MLC were used for treatment sequence delivery. Treatment plans were generated with the Monaco 3.2/3.3 (Elekta) treatment planning system. Fifty different treatment plans for 13 different patients with lesions in lung (8), liver (2), stomach (1) or esophagus (2) were evaluated. Eight plans were generated using an MLCi2 collimator (fraction dose D fx = 2 Gy), 28 plans were generated with an Agility collimator (D fx = 1.8-3 Gy), and 14 plans were planned for a VersaHD LINAC for hypofractionated flattening filter free treatments (D fx = 5-12 Gy). All treatment deliveries were performed with one 360 • VMAT arc and increment factors ∆Θ between 10 • and 40 • . The collimator angle was set to 0 • , what resulted in a MLC leaf motion perpendicular to the axis of gantry rotation. Treatment volumes ranged from 12 to 2273 cm 3 .

2.B. Estimation of the dose distribution in the presence of motion
All treatment sequences were delivered on a tablemounted, two dimensional ionization chamber array (Matri-XX Evolution, IBA, Schwarzenbruck, Germany) inside a water equivalent phantom (Multicube Lite, IBA). The dose was recorded with a sampling time of ∆t = 100 ms and stored in separate matrices. To be able to detect interplay effects along the three coordinate axes, all treatment sequences were delivered on the detector, aligned once in vertical and once in horizontal direction. A purpose written  ® (Natick, MA) based script package was developed to evaluate the dose distribution in the presence of target motion. First, every matrix of the two dimensional measured static dose grid is interpolated to 2 mm spatial resolution. As shown in Fig. 1, the individual dose planes are then moved in either craniocaudal (CC), left-right (LR) or anterior-posterior (AP) direction. The nonzero component x t of the target motion vector depends on the time after the start of the delivery t. It is determined according to Lujan et al., 10 where A is the peak-to-peak amplitude, T the breathing period, and φ the starting phase with respect to the treatment delivery start. The asymmetry factor k is in the following set to The EUD is calculated with a = −25 for all N dose pixels D l with more than 70% of the maximum dose of the static dose distribution. The threshold dose value is motivated by the following consideration: For measurements in a phantom, the PTV structure is not available. Therefore, the maximum dose serves as a reference value. The dose within the target should be within 95%-107% of the prescription dose. 33 Additionally, the dose distribution in the phantom does not ideally resemble the dose distribution in the patient which can result in a further target dose inhomogeneity. Since in the patient the 95% isodose line should cover the PTV outline, a threshold value lower than 88% of the maximum dose, but still within the high dose gradient area, should be chosen to assure that all relevant voxels are being considered. The factor a in Eq. (2) allows for an explicit selection of the degree of the minimum dose component (for a = −∞) or the mean dose component (for a = 1). With the choice of a = −25, a result close to the minimum ROI dose but still sufficient dependency on the size of the underdosed volume is provided.
The minimum relative equivalent uniform dose EUD r, m is obtained by the ratio of the EUD of the dose slice if motion was present EUD dyn and the EUD in the static case EUD stat The EUD r, m is obtained over all breathing starting phases. Phase averaged, the relative EUD equates the EUD of a pure dose blurring case, i.e., without consideration of the temporal substructure. The static case typically shows a nonuniform distribution in the target volume. Therefore, to be sensitive to motion-induced effects, the ratio of dynamic over static dose distribution was considered.

2.C. Estimation of periods with a minimum EUD r,m
For motion in CC-direction, the jaws, which do not follow a sweeping window pattern, are the beam limiting devices. Motion in LR-direction was reported [34][35][36] to be generally less pronounced than motion in AP-direction. Therefore, for a collimator angle of 0 • , target motion in AP-direction is most likely prone to SWAT-specific dosimetric errors and will be considered in this section.
In the presented model, dosimetric deviations resulting from target motion parallel to the central beam axis are neglected, and solely dose gradients along the beam penumbra are considered. This assumption is motivated by the fact that for a 6 MV beam with flattening filter and a (5×5) cm 2 square field, the dose gradient along the central axis at 10 cm depth is dD/dx pdd = −3%/cm of the maximum dose while the mean gradient between the 20% and the 80% penumbra dose at the same depth is dD/dx CL = ±61%/cm of the maximum dose.
To model a SWAT beam delivery, a treatment sequence with a sweeping collimator gap that performs s sweeps with an amplitude A c and period t tot /s, while the gantry performs one 360 • considered. The dynamic of the resulting dose gradient areas is depending on two quantities: the current gantry angle and the position of the collimator gap. The impact of the gantry angle on the orientation of the penumbra with respect to the AP-direction of target motion can be described by a sine function of the current gantry angle, whereas the periodic collimator motion is modeled by another sine function. This yields an AP-component of the penumbra dynamic x p of The number of sweeps s is a function of the increment factor ∆Θ: s = 360 • /(2 · ∆Θ) defined during treatment planning. As a first order approximation for periodic breathing patterns, the target AP-motion can be assumed to perform a periodic motion. For k = 1, a constant amplitude A, and ϕ = 2φ + π/2, Eq. (1) can be rearranged to In this model, the angular frequency ω t = 2π/T corresponds to the dominant amplitude square component of any Fourier transformed irregular motion pattern. It is assumed that in SWAT, the gantry angle changes with a constant angular velocityΘ = 2π/t tot . Therefore, the gantry angle Θ 1 is reached at t = 1/4 t tot and Θ 2 at t = 3/4 t tot .
In order to achieve a maximum dosimetric deviation in the target volume, the two functions x p and x t have to be both extremal at Θ 1 and Θ 2 and both either in-phase (+) or of opposite-phase (−) around points with extremal AP-components Θ 1 and Θ 2 sgn In the opposite-phase case, the target voxels avoid the beam to a maximum degree. For a target volume that is larger than the motion amplitude, this results in an underdosed central target region and an overdosed peripheral target region along the motion direction as well as an increased dose deposit outside the target. In the in-phase scenario, the target voxels follow the MLC sweep. Therefore, they accumulate more dose in the central target region while a sufficient target coverage at peripheral target regions is not guaranteed. Between Θ 1 and Θ 2 an integer number of breathing cycles can take place [T = t tot /(2b), b ∈ N] without changing the resulting behavior. From the above conditions, it follows that s has to be a positive, odd integer number. Solving the resulting system of equations yields for the two free target motion parameters T and ϕ with Z ∈ Z and p = 0 for the opposite-phase case and p = 1 for the in-phase case. Equation (8) implies that the periods of maximum in-phase behavior are identical to the periods of maximum opposite-phase behavior.

2.D. Analysis of correlated motion patterns in SWAT
During plan delivery, the gantry angle, the cumulative monitor unit, and the individual MLC and jaw positions were recorded in LINAC logfiles with a frequency of 4 Hz. The Pearson correlation coefficient r(x c ,x t ) for the collimator position component x c along a specific axis (AP, LR, CC) and the target motion x t in the same direction for a discrete sequence can be determined with where i from 1 to n are the individual discrete times. Parallel to the direction of MLC leaf motion, this calculation is performed for all MLC leaves j that are not behind the craniocaudal jaws. Perpendicular to the direction of MLC leaf motion, the jaw positions are taken. The maximum correlation coefficient r max (A,T) of all relevant leaves or jaws j rel and of all starting phases φ rel is obtained for each (A, T) combination Besides this purely geometrical consideration of the two dynamic patterns, a weight factor according to the dosimetric contribution of the segment has to be added. The relative weight is a product of the number of monitor units n MU delivered in a time interval and the attenuation factor r att (FS equ ,d). This factor is the relative dose of an equivalent square field of size FS equ of the current segment at a depth d. d is the water equivalent distance traversed between patient or phantom surface and the beam isocenter. Therefore, the dose weighed mean can be expressed as and the weighed standard deviation is 3. RESULTS Figure 2 shows an analysis of two representative treatment sequences. In the upper row, the minimum relative EUD as a function of the applied breathing period for different peak-topeak breathing amplitudes is displayed, while the lower row shows the corresponding correlation coefficients between a target motion with A = 2 cm and the collimator dynamic APcomponent. In both cases, motion in AP-direction was applied with the detector array being positioned in vertical alignment. The figures in the left column were generated with a treatment sequence of a target volume in the central right lung and were delivered with an MLCi2 collimator (t tot = 197 s, ∆Θ = 20 • ), while the figures in the right column represent a treatment sequence for a target volume in the distal left lung, delivered with an Agility MLC (t tot = 61 s, ∆Θ = 20 • ). The black circles indicate the predicted EUD r, m period times according to Sec. 2.C for Z ∈ {−6,−5, ..., 3}. The minima with the lowest EUD r, m were found for Z = −5 for both cases. In the presented cases, the deviations between predicted and measured periods were (0.01 ± 0.08) s and (−0.04 ± 0.25) s.
In order to avoid high correlation coefficients within the range of human breathing periods, the maximum allowed doserateḊ max of a treatment sequence can be reduced. The result of limitingḊ max on the periods with maximum cross correlation is presented in Fig. 3. The correlation coefficient between a target motion with different (A,T) settings and a collimator dynamic was obtained for three differentḊ max : 560, 400, and 270 monitor units per minute. The collimator dynamic of the left MLC leaf bank of the treatment sequence presented in Fig. 2 (right) was used for the determination of the correlation coefficients. For the presented treatment sequences, the resulting delivery times were 61, 76, and 109 s, respectively.
A dose delivery with a partial arc can lead to a further degradation of the dose homogeneity due to less dosimetric contributions from segments that move parallel to the motion direction. Therefore, the treatment sequence used in Fig. 2 on the right was limited to gantry angles between 40 • and 160 • with an increment of ∆Θ = 20 • . The resulting analysis is shown in Fig. 4. A correlation coefficient of 0.94 resulted in a minimum EUD r, m of 92.0%.
F. 3. Correlation coefficients for different maximum dose ratesḊ max . The collimator dynamic of the left MLC leaf bank of the treatment sequence presented in Fig. 2 (right) was used.
In Fig. 5, a measured dose distribution of the Agility MLC treatment sequence used in Fig. 2 without target motion is shown on the left. A target motion with T = 6.1 s, A = 2 cm, and φ = 90 • (Z = −4) was then applied. The corresponding dose differences to the static case are shown in Fig. 5 on the right. While the middle figure was generated with motion in CC-direction, the right figure was generated with motion in AP-direction.
The effect described in Sec. 2.C that in the opposite-phase case, dose from the central PTV regions gets deposited in peripheral regions due to an opposing motion of MLC leaves and patient breathing is demonstrated in the right figure.
The dose distributions represented the minimum EUD r,m of this treatment sequence.
The EUD r, m was determined for all treatment sequences and the four different motion scenarios. A histogram of the resulting EUD r, m distribution of all simulated cases was F. 4. Treatment sequence with a partial arc from 40 • to 160 • gantry angle. Upper row: Relative maximum equivalent uniform dose EUD r, m (A, T). Lower row: Maximum correlation coefficient r (x c , x t ). created with an EUD r, m bin width of 1%. This distribution is shown in Fig. 6.
On average, the EUD r, m for the target motion directions parallel to the MLC leaf motion direction was 0.7% lower when compared to a target motion perpendicular to the MLC leaf motion direction.  Table I shows the dependency of the applied number of sweeps on the EUD r, m (AP V ). Solely, treatment sequences of target volumes larger than 50 cm 3 were considered to ensure a pronounced sweep behavior. In agreement with the model presented in Sec. 2.C for 9 sweeps, the EUD r, m (AP V ) had the lowest mean value.

4.A. Suitability of the evaluation method
Various methods to evaluate the effect of target motion on the resulting dose distribution in IMRT and VMAT, such as Monte Carlo simulations 19 or measurements in suitable phantoms 14,16,18 with representative parameter settings have been performed. However, in SWAT, an appropriate selection of the parameters (A,T,Φ) is crucial for an accurate estimation of the worst possible treatment sequence. As shown in Fig. 2, the treatment plan robustness varied for different breathing periods. Increasing the target motion amplitude on the other hand showed a similar behavior over the entire period range.
Mean and standard deviations, as well as gamma analysis, 8 as measures for dosimetric errors are insensitive if the deviations affect only a small volume. If one considers exclusively minimum doses on the other hand the size of the underdosed volume is not taken into account. DVH based methods are either multidimensional 37 or use volumetric and dosimetric scalar quantities, such as V 50 (%) or D 95 (Gy), to evaluate the data. This kind of evaluation is typically being performed 38 for several dose or volume levels. In order to be able to observe the explicit dependence of the breathing period and amplitude on the treatment sequence robustness, the scalar EUD r, m serves as an appropriate measure.
The Pearson correlation coefficient as presented in Eq. (11) is a measure of a linear correlation between the target and collimator motion. Nonlinear behavior will therefore lead to a reduced maximum correlation coefficient. However, if a random-like segment order was present, the correlation coefficient was close to 0 and a relatively low correlation coefficient threshold of 0.5 was sufficient to detect minima in EUD r, m for all observed cases.
The use of a solid water phantom as a substitute for lung tissue will affect the dosimetric outcome. Therefore, the relative minimum EUD serves in this context more as a relative measure to indicate nonrobust frequencies than an absolute measure for the actual target EUD. The same applies for the choice of the parameter a in Eq. (2). Furthermore, in this study, solely the dosimetric consequences of a moving target on the target volume were considered. In a realistic patient scenario, a nonrigid organs at risk (OAR) deformation will additionally affect the OAR doses, as well.

4.B. Dosimetric deviations in SWAT
A minimum in the EUD r, m results from a high correlation coefficient of the delivery pattern parallel to the target motion. If an extremal collimator elongation is opposed by an elongation in the same direction, an integer number of collimator sweeps between the two opposing beams will result in another minimum in EUD r, m . It was shown in Fig. 5 that the manifestation of interplay effects changes in SWAT from spotlike artifacts, as presented in Fig. 5 (middle) or by Bortfeld et al. 11,12 for IMRT, to widespread underdosed areas as presented in Fig. 5 (right). Since all treatment plans were generated with 0 • collimator angle, the sweeping windows moved parallel to the direction of gantry motion and caused dose deviations for target motion perpendicular to the gantry rotation axis to be on average 0.7% higher than the deviations resulting from motion parallel to the gantry rotation axis. The presented approach focuses on dose homogeneity in the PTV, whereas the CTV is the relevant structure that will affect the treatment efficiency. For standard IMRT, the location of the underdosed volume within the PTV remains unknown with this approach. In the case of SWAT, a nonrobust breathing period will, provided an unfortunate starting phase is chosen, generate an underdosed area in the central region of the PTV [ Fig. 5 (right)]. As CTV-PTV margins are usually chosen symmetrically around the CTV, (Ref. 25) the underdosage will occur in the CTV.

4.C. Robust treatment planning in SWAT
While the motion amplitude is used to determine the PTV margin during treatment planning, 33 the breathing period is typically not considered. One option to minimize interplay effects in SWAT is to increase the increment parameter to reduce the number of sweeps. The human breathing is typically between 9 and 23 breathing cycles per minute. 36 Provided that the rotational velocity of the gantry does not exceedΘ ≤ 2π/60 s in combination with the results from Table I, care should be taken if choosing an increment factor of ∆Θ = 20 • .
Using multiple arcs with different increment factors or collimator angles during treatment delivery will further desynchronize the delivery pattern. For treatment plans, in which the treatment volume is located laterally, most of the dose will originate from a subarc around 90 • or 270 • gantry angle, where the motion component in AP-direction is maximal. These cases have to be considered with special care, since an averaging of opposing beamlets and dosimetric contributions from beamlets, where the target moves parallel to the central axis of the beam, will not be present.
For standard IMRT, Jiang et al. 14 recommend to use a low dose rate when treating a moving tumor without any motion mitigation techniques. The presented approach can be considered as a quantitative method to obtain the fastest robust treatment sequence for VMAT with a sweeping window collimator dynamic. An iterative decrease of the maximum dose rate can be used to alter the temporal behavior of the delivery pattern x p . For x t , a realistic target trajectory can be obtained from 4D-CT data. Due to interfractional changes in the breathing pattern, it is recommended to test against additional realistic breathing patterns. This procedure provides the possibility to produce robust sequences without having to reoptimize the treatment sequence. As the presented minima in EUD r, m are typically temporally local, a small decrease in treatment time will suffice to move out of a minimum. Chin et al. 9 adjusted the increment factors to match the patient breathing, obtained from 4D-CT data. While this is the optimum method if the patient breathing trajectory is known, the method presented in this work does not rely on an exact knowledge of the target dynamic during treatment delivery. Although it does not warrant a reduction of safety margins, it assures a fast convergence to a pure dose blurring situation.

CONCLUSION
With the advent of fast delivery techniques less robust treatment plans within the human breathing range can be generated. In SWAT, the correlation coefficient between the collimator pattern and realistic motion trajectories serves as an estimate for the plan robustness. By iteratively decreasing the maximum dose rate, a robust treatment sequence can be warranted. No degradation of the relative minimum EUD above 8.0% compared to the static case was found.

ACKNOWLEDGMENT
This work was supported by an Elekta research grant held by F.W. and F.L.