If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
UCLouvain, Institut de Recherche Experimentale et Clinique (IREC), Center of Molecular Imaging, Radiotherapy and Oncology(MIRO), 1200 Brussels, BelgiumKULeuven Department of Oncology, Laboratory of Experimental Radiotherapy, 3000 Leuven, Belgium
Directional MidP method incorporates motion directionality in PTV definition.
Higher 4D D95% and D99% are achieved with the proposed method.
PTV volume is similar compared to classical MidP strategy.
Planning target volume (PTV) definition based on Mid-Position (Mid-P) strategy typically integrates breathing motion from tumor positions variances along the conventional axes of the DICOM coordinate system. Tumor motion directionality is thus neglected even though it is one of its stable characteristics in time. We therefore propose the directional MidP approach (MidP dir), which allows motion directionality to be incorporated into PTV margins. A second objective consists in assessing the ability of the proposed method to better take care of respiratory motion uncertainty.
11 lung tumors from 10 patients with supra-centimetric motion were included. PTV were generated according to the MidP and MidP dir strategies starting from planning 4D CT.
PTVMidP dir volume didn’t differ from the PTVMidP volume: 31351 mm3 IC95% [17242–45459] vs. 31003 mm3 IC95% [ 17347–44659], p = 0.477 respectively. PTVMidP dir morphology was different and appeared more oblong along the main motion axis. The relative difference between 3D and 4D doses was on average 1.09%, p = 0.011 and 0.74%, p = 0.032 improved with directional MidP for D99% and D95%. D2% was not significantly different between both approaches. The improvement in dosimetric coverage fluctuated substantially from one lesion to another and was all the more important as motion showed a large amplitude, some obliquity with respect to conventional axes and small hysteresis.
Directional MidP method allows tumor motion to be taken into account more tightly as a geometrical uncertainty without increasing the irradiation volume.
]. Radiation therapy of mobile tumors therefore requires dedicated methodologies. In contrast to sophisticated motion mitigation techniques based on the synchronization between beam delivery and instantaneous tumor location (e.g. gating or tracking), margin strategies can be easily implemented on conventional linear accelerators (linacs) without any addition of specialized device. This approach thus remains the most widespread in clinical routine [
Based on a four-dimensional computed tomography (4D-CT), two methods are commonly used to compute safety margins. The first involves the creation of an Internal Target Volume (ITV) encompassing all tumor positions in the respiratory cycle. To this volume is added a margin compensating for other geometrical uncertainties in order to generate the Planning Target Volume (PTVITV) [
], the second method, including both Mid-Position (MidP) and Mid-Ventilation (Mid-V) strategies, handles breathing motion as a random error that smears the dose distribution. In the Mid-V method, delineations are performed on the 4D CT phase closest to the time-weighted average of all tumor positions [
]. This frame of the 4D CT data set is however subject to motion-related artefacts degrading image quality. Moreover, a systematic error is introduced because of the hysteresis of tumor motion. A refinement of the concept has consequently been devised, the MidP method [
]: target volumes are contoured on the MidP-CT which represents the patient's anatomy at its time-weighted average position over the respiratory cycle. A single volumetric expansion is then performed to generate the PTVMidP. It has the advantage of being on average 6.5% to 13.9% smaller than PTVITV, thus lowering the probability of radiation-induced complications in surrounding healthy tissues [
]. Among simplifications operated in the van Herk’s margins model, tumor motion directionality is not considered. Indeed, only motion variances along the three conventional spatial dimensions are integrated in the calculation of geometrical uncertainties.
A concept of directional margin was firstly proposed by Alan McKenzie [
]. Briefly, his method involves successive dilations of the target volume with ellipsoids whose semi-diameters are aligned with the eigenvectors of the covariance matrices of systematic and random uncertainties. The length of their semi-diameters is equal to the square-root of the eigenvalue of its associated eigenvector. However, breathing motion is handled differently than the classical MidP approach. Actually, its contribution to the margin recipe is added linearly to other geometric errors combined in a quadratic way [
]. To our knowledge, no study has evaluated the potential benefit of a more refined description of the probability density function of tumor motion. The first goal of this work is therefore to propose a solution allowing integration of motion directionality in the MidP margin recipe. The second objective aims at quantifying the gain of our method in terms of tumor dosimetric coverage compared to the classical PTVMidP approach.
2.1 Definition of geometry
Let (x,y,z) be an orthogonal basis whose origin (0,0,0) corresponds to the mid-position of CTV centers of mass (CTVCM) among all 4D CT phases. The x, y and z axes correspond to the left-to-right, antero-posterior and cranio-caudal directions, respectively.
2.2 Consistent integration of motion directionality in margin recipe
The N systematic errors are presumed to be independent of each other and the same applies to the relationship between their spatial components. The probability density function of geometrical uncertainties is supposed to be Gaussian. The first CTV extension compensating for systematic errors can thus be carried out in the same way as van Herk et al [
] by defining a confidence ellipsoid. This volume, called systematic PTV (PTVs), ensures CTV coverage for a given proportion of the patient population.
Margins m can be calculated as follows in the orthogonal basis (x,y,z) described in section 2.1:
Coefficient α determines the percentage of patients whose CTV receives the minimal dose level that is allowed. Σ(i) represents the standard deviation of a given systematic error.
Concerning random errors, preliminary assumptions are hypothesized to be identical as systematic one. In the classical model, perfect dose conformality is assumed for a spherical CTV of radius . Along direction defined by an arbitrary vector quantifying systematic errors, the minimum dose level then happens at the point and can be calculated, as demonstrated by van Herk et al [
], by the convolution product between the probability density function of random errors and planned dose named Dnominal :
with S, a step function whose vector delimits an ellipsoid defining the volume of the 50% isodose:
refer to the Gaussian probability density functions of combined random errors.
This approach leads, in clinical practice, to the computation of anisotropic margins compensating for random errors in the classical orthogonal basis regardless of the motion directionality. However, if tumor motion shows some obliquity, its spatial components cannot be treated as statistically independent. The mathematical corollary implies that the probability density function of combined random errors cannot be approximated by the product of its marginal distributions due to non-null covariances between these dimensions. It is therefore necessary to decorrelate these components. Let Atm be the covariance matrix of tumor positions along its trajectory. It is calculated from the sample of CTVCM positions reconstructed by the 4D CT:
The covariance matrix resulting from the combination of the K independent random errors can therefore be written as follows:
Since tumor motion is the only random uncertainty considered as directional, we can simplify as follows:
It must therefore be possible to define another orthogonal basis in which becomes diagonal. Such an orthogonal basis corresponds to the one defined by the eigenvectors of :
with ,, being the eigenvalues of the covariance matrix defined in (x,y,z) associated with their corresponding eigenvector : , and respectively.
]. We thus arrive at an identical margin recipe but in a different orthogonal basis:
Coefficient defines the minimum dose level delivered to the CTV.
As the penumbra of the dose distribution is in good approximation isotropic, the following equality applies as demonstrated in Appendix 1:
11 lung lesions from 10 patients treated by stereotactic body radiotherapy (SBRT) from 2018 to 2020 for stage I non-small cell lung carcinoma were included. A 4D-CT (Aquilion LB, Canon medical systems corporation, Japan) without contrast injection and retrospectively reconstructed in 10 temporally equivalent phases was performed. The slice thickness of the 4D-CT frames was set to 2 mm. The planning CT was performed in supine position with hands above the head and with a vacuum bag as immobilization device. Abdominal compression (AC) was used in some cases depending on the instructions of the patient's physician. Patients underwent audio-coaching during CT acquisition based on their spontaneous respiratory rate and their external breathing signal was captured using the GateCT surface guided system (Vision RT Ltd, UK). The same patient-specific audio-coaching was applied throughout treatment delivery. Lesions features are summarized in Table 1.
Table 1Summary of lesions characteristics. Angle = angle between major motion direction and conventional axis associated with the largest amplitude. Linear coefficient = root mean square of linear correlation coefficient of tumor trajectory projected in the three spatial planes. * = Application of abdominal compression.
]. Delineation of organs at risk (OARs) and GTV were performed in Raystation Planning system (clinical version 9B, RaySearch Laboratories, Stockholm, Sweden) on the MidP CT by the same radiation oncologist.
As shown in Fig. 1A for a simulated motion using the equation of Lujan et al [
], PTVMidP was then generated by expanding the GTVMidP via safety margins calculated according to van Herk's formula as:
where Σ and σ refer to the standard deviations of systematic and random errors, respectively. Labels “TM”, “BL”, “SETUP” and “D” correspond to tumor motion, baseline shift, set-up and delineation, respectively. Except for tumor motion, their values were taken from the literature and adapted to our institution protocols specificities [
Regarding tumor motion, corresponds to the residual error of the non-rigid co-registration algorithm and was equal to 0.5 mm. coincides with the standard deviation of tumor motion. This parameter was individualized for each patient based on the standard deviation of the GTV center of mass coordinate set derived from all 4D CT phases.
Coefficient α was fixed at its conventional value of 2.5 ensuring that 90% of the patient population receives a given dose level. Coefficient β was determined according to Thomas et al [
]: an overdosage was applied within the PTV so that the prescribed dose was at most equal to 90% of the maximum dose level. Coefficient β equals 1.28 and guarantees that the GTV receives at least 90% of the prescribed dose in all planning scenarios.
2.4.2 PTV directional MidP
The same MidP CT and GTV contours as for the PTVMidP were used. A first dilation of GTVMidP in the classical orthonormal axes was operated in order to generate the PTVs. Margins were derived from the first square root of the van Herk’s formula multiplied by the coefficient α. This intermediate volume is identical between PTVMidP and PTVMidP dir.
As shown in Fig. 1B, an expansion of PTVs was then performed using an ellipsoidal structuring element whose axes are aligned with the eigenvectors of . The length of ellipsoid axes were equal to the eigenvalues associated with their respective eigenvectors to which the axes were parallel. This second volumetric extension leads to the creation of the PTVMidP dir. The above-described operations were executed using a Matlab script developed in the Open Reggui software environment.
Treatments were planned by the same medical physicist in the treatment planning system Raystation (from Raysearch labs, version 9B) for a conventional linac (Infinity Elekta linac). Once optimized, treatment plans were computed on the average CT thanks to a convolution collapsed cone dose calculation algorithm. Voxel intensities of the average CT correspond to the mean of the 4D-CT phases, thus faithfully representing the average lung density during respiratory cycle [
], planning objectives were the following: D99% (dose received by 99% of the PTV volume) above 90% of the prescribed dose, D95% equal to 100% of the prescription dose and D2% between 110 and 140% of the prescribed dose. Dose constraints to organs at risk were kept below literature recommendations [
], was determined by the accumulated dose to the extended GTV along its trajectory. It corresponds to the dilated GTV according to van Herk’s formula for all geometrical uncertainties except tumor motion (ΣTM = 0 and σTM = 0). For this purpose, 3D dose maps were accumulated in the Mid Position using a hybrid method combining rigid and non-rigid methods. Dose accumulation workflow is illustrated in Fig. 2.
As the GTV is an anatomical structure, dose received by GTVi associated with phase i was summed at GTVMidP location by the deformation vector field i (DVFi) obtained with the same deformable image registration (DIR) algorithm as to reconstruct corresponding MidP CT. Contribution of this distorted dose outside GTVMidP was set to 0. Since this accumulation critically depends on the DIR accuracy [
], a correction was applied if centers of mass of GTVi deformed by DVFi and GTVMidP were not identical: a vector ci matching the center of mass of the deformed GTVi to that of GTVMidP was applied to the DVFi deformed dose. It is known that lung tumors “track the dose” along their trajectory when type B dose engines are used for treatment planning [
]. Hence, 3D dose map has been particularized for each phase by recomputing the original plan on each 4D CT phase thanks to the Raystation’s convolution collapsed cone algorithm.
As for the rest of the extended GTV volume, dose distribution was rigidly shifted by a translation vector vi joining the center of mass of GTVi and GTVMidP. Actually, the Mid position is shifted towards expiratory phases and the selected tumors have a large motion amplitude. Highly mobile anatomy in the close environment of the tumor could consequently lead to an underestimation of the delivered dose to the extended GTVMidP with an exclusively non-rigid 4D dose accumulation. Indeed, substantial distension of lung tissue during respiratory cycle may result in the accumulation of dose voxels located outside the extended GTVi within the extended GTVMidP. The extended GTV being a purely geometrical volume, this kind of accumulation has no physical meaning. Contribution of this rigidly translated dose within GTVMidP was nullified. Finally, the distorted dose has here not been individualized on each phase and corresponds to the original one computed on the average CT. Density of its voxels is indeed closer to that of the GTV than on a single 4D CT phase and thus offers an acceptable trade-off for 4D dose accumulation [
The final 4D delivered dose is the sum of both accumulations.
Mathematically, the latter could be written as follows:
where pi corresponds to the tumor presence probability in phase i. Given the 4D CT reconstruction in 10 temporally equivalent phases, all pi are uniformly equal to 0.1.
Loss of dosimetric coverage induced by tumor motion was calculated by the relative differences between the Dx% (minimum dose received by x % of the volume of interest) planned to the PTV and that delivered to the extended GTV. These differences were reported for the D99%, D95% and D2% and were noted ΔD99%, ΔD95% and ΔD2%, respectively. The absolute difference of ΔDx% between MidP and directional MidP method corresponded to the final metric quantifying the difference in tumor dosimetric coverage between the two approaches.
2.7 Statical analysis
Statistical analyses were performed on IBM SPSS Statistics 26 software. The means of continuous variables were compared using two-tailed and paired t-tests when data sample were normally distributed or Wilcoxon tests in case of non-Gaussian distributions. The significance threshold was set at 5%.
3.1 Motion features
As shown in Table 1, major motion axis amplitude ranges from 10.26 mm to 35.74 mm. Lesions located in the upper or middle lobe show the smallest amplitudes. Two tumors (n° 5 and 10) feature a significant hysteresis of more than 8 mm.
Angles reported in Table 1 correspond to the deviation between conventional axis related to the largest motion amplitude and major motion axis direction defined by the eigenvector associated with the highest eigenvalue. There seems to be an inverse relationship between motion angulation and its amplitude. Linear coefficient refers to the root mean square of the linear correlation coefficients (LCC) of the projected tumor trajectory in the 3 spatial planes xy, xz and yz. A coefficient equal to 1 would mean a purely linear tumor motion morphology without any hysteresis. As can be observed in Table 1, hysteresis amplitude decreases with the linear coefficient.
3.2 Margin and PTV volume
As seen in Table 2, the directional MidP margin associated with the first principal component is systematically higher than the largest MidP margin. In contrast, the margin related to the third principal component is always smaller than the smallest MidP margin. Graphs in Fig 3 show a clear proportional relationship between margin of the first principal component and major motion axis amplitude. The same trend is observed between the hysteresis amplitude and the margin associated with the second principal component.
Table 2PTV random margins and volumes for MidP and directional MidP methods. Volume ratio refers to PTVMidPdir volume divided by PTVMidP volume.
From a volumetric point of view, the average PTVMidP dir volume did not significantly differ from the mean PTVMidP volume: 31351 mm3 IC95% [17242 – 45459] vs. 31003 mm3 IC95% [ 17347–44659], p = 0,477. However, as shown in Fig 4a, their spatial distribution appears not identical: PTVMidP dir is wider than PTVMidP along the main motion axis but narrower along perpendicular directions.
3.3 GTV dosimetric coverage
Plans’ conformity indices for both approaches were not significantly different. DSI index was 0.962 CI95% [0.955–0.968] vs. 0.957 CI95% [0.951–0.963], p = 0.067 for MidP and MidPdir respectively. Paddick's conformity index was 0.925 CI95% [0.912–0.938] vs. 0.916 CI95% [0.905–0.928], p = 0.063 for MidP and MidPdir respectively.
Directional MidP significantly improved dose coverage for D99% (mean difference of 1.087% CI95% [0.346–1.858], p = 0.011) and D95% (average difference of 0.739% CI95% [0.297–1.401], p = 0.032). The improvement of D99% and D95%, as shown in Fig. 5, is at first approximation driven by the major motion axis amplitude. This correlation is however subject to significant variability. By contrast, the D2% does not differ between the two techniques (average difference of 0.104 [-0.145–0.316], p = 0.429). The results for ΔD99%, D95% and D2% are compiled in Table 3.
Table 3ΔD99%, ΔD95% and ΔD2% for MidP and directional MidP methods.
Compared to PTVMidP, PTVMidPdir shows a slight but statistically significant trend towards improved D99% and D95% without any increase of the irradiation volume. D2% is not differently affected between both approaches. Higher tumor dosimetric coverage is particularly desirable as it allows potentially for a better local control [
]. Translation of the study results into clinical significance is thus not well established.
Beyond statistical analyses, it worth mentioning that the benefit of the directional MidP method is highly fluctuating from one lesion to another. Several determinants could theoretically explain this variability. As illustrated in Fig. 5, it seems that dosimetric gain of PTVMidP dir is proportional to the motion amplitude. Tumor motion being the only error treated differently by the proposed method, this may be due to its larger weight with regard to the other geometrical uncertainties. It could thus be conjectured that the dosimetric advantage of directional MidP margins becomes also apparent for motions of lower amplitudes in case of larger coefficient β e.g., normofractionated treatment of lung tumors [
]. Nevertheless, tumors with similar amplitudes can benefit very differently from PTVMidP dir e.g., lesions n° 8 and 5 or lesions n° 3 and 11. A first element could come from motion deviation from conventional axes. It is obvious to conclude that if principal motion axes are identical to the classical orthonormal axes, no difference will be observed between PTVMidP and PTVMidP dir. Lesions n°8 and 11 demonstrate indeed a larger motion obliquity compared to lesions n°5 and 3. A second determinant could be related to tumor trajectory shape: van Herk's model assumes a null covariance between the different spatial components of geometric uncertainties. It means that LCC of tumor positions between conventional axes are also equal to zero. As illustrated in Fig. 6B, this approximation is quite acceptable in case of large hysteresis where tumor trajectory projected in the 3 spatial planes exhibits scatterplots with circular-like distributions characterized by a low LCC. However, for motions with linear morphology, this assumption is no longer realistic as represented in Fig. 6A. Motion amplitude being the same, it is therefore logical to observe higher gain of PTVMidP dir for lesions with mean LCC close to 1 (n°8 and 11) in comparison with tumors having lower mean LCC (n°3 and 5). The influence of other motion pattern parameters on the absolute level of tumor dosimetric coverage highlighted by Wanet et al [
] such as motion amplitude or GTV size remain applicable. When motion amplitude increases or when random errors become non-negligible with respect to the GTV size, the model applicability is questionned. Some extreme cases can thus lead to inconsistent results regarding above described factors underlying PTVMidP dir benefit as for lesion 2, whose very small volume is of the order of 0.8 cc. Finally, 3 lesions with motion amplitude of less than 15 mm display a lower dosimetric coverage with directional margins (maximum difference of 0.49% for D95%). Apart from the fact that the clinical impact would be negligible, we attribute this small difference to the limit of similarity between PTVMidP and PTVMidP dir plan conformity for lesions where both approaches are most probably equivalent.
The patient preparation could affect the observed tumor motion, as well as its directionality. The impact of AC on tumor motion amplitude is well documented. It involves a significant reduction of motion amplitude for most of the liver and lung tumors [
] demonstrated that motion directionality was stable even when AC was applied. Indeed, AC was part of the treatment protocol for 5 of the 11 patients in their cohort. Unfortunately, the small sample size did not allow further conclusions to be drawn regarding its influence on the degree of motion obliquity. In our cohort, only 2 patients underwent AC. The same interpretative precautions as with the study of Worm et al [
Margins generated by the directional MidP strategy consistently incorporate major motion axis as well as hysteresis amplitudes. This results in a wider PTV along the main motion direction but narrower perpendicularly. From a dosimetric point of view, it involves, as shown in Fig 4b, an overdosage at the extreme inspiratory and expiratory positions of breathing cycle corresponding to those where the tumor spends the most time. This explains the improved dose coverage achieved with PTVMidP dir. Directional margin validity will therefore depend on the stability of motion directionality and hysteresis. Fortunately, these two parameters are quite constant over time [
]. These metrics calculated on a 4D CT can thus be reliably extrapolated to the entire treatment. Nevertheless, the proposed method does not solve the important methodological weakness of margin strategy which consists in neglecting variability of tumor motion amplitude [
] has been indeed the first to propose aligning margins to the directions of tumor motion. Our method is similar to his as it is based on the decomposition of the covariance matrix of the geometric errors into eigenvalues. In contrast, motion uncertainty is treated differently. After the CTV expansion for systematic errors, a directional dilation involving only tumor motion is performed. This operation is identical to the proposition of Worm et al [
]. They advocated that the alignment of motion’s margin should be made along the principal components of tumor motion. As part of of the directional MidP strategy, we have shown that the alignment of the directional margins only involves random uncertainties. It must also be performed along the principal components of the covariance matrix of the combined random errors and not those of the covariance matrix of tumor motion alone. The orthogonal basis defined by their eigenvectors is not the same indeed. Worm et al [
] also conjectured that the directional PTV volume would be smaller than the classical volume. In our study, the volumes obtained were similar between the investigated PTV definition methods.
In this work, directionality of baseline shift was not considered. Directional MidP method could be easily enhanced by assuming a non-null covariances between dimensions of this geometrical uncertainty. Specification of its covariances would require a population-based quantification trough intra-fraction tumor positions monitoring. This is however out of the scope of this study.
Several limitations of this study need to be pointed out. The first concerns the use of a convolution collapsed cone algorithm for dosimetric analysis. Although an excellent agreement with Monte Carlo calculations is achieved in most cases, a slight overestimation of the dose at the target volume can be observed in some extreme conditions (small lung tumor) [
]. Since comparisons between MidP and directional MidP methods were performed on the same CTs for the same patients, it is reasonable to assume that this uncertainty has been propagated symmetrically between both approaches and that its net impact was negligible. Second, the 4D dose accumulation method used in this work did not consider the interplay effect. Residual incertitude related to tumor motion should be very small as it has been shown to be clinically insignificant [
]. In the case of a large baseline shift, the study results should thus be extrapolated with caution. Lastly, rigid accumulation of the dose at the extended GTV did not take into account the dosimetric effects of tumor rotation.
The directional MidP method allows tumor motion to be better taken into account as a geometrical uncertainty without increasing the irradiation volume. This results in a higher dose delivered to the GTV than with the MidP approach, especially for highly mobile tumors whose trajectory has a linear shape and deviates from classical axes.
Loïc Vander Veken is funded trough a FNRS grant FC N°33411.
Kevin Souris is funded by the Walloon region (MECATECH/BIOWIN, grant number 8090).
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Appendix n. 1
Projected penumbra variances along axes of a new orthogonal basis (x’,y’,z’) can be written as follows :