The ROSMA-PTV concept and the impact of rotational setup uncertainties on the tumor control probability in canine brain tumors

Purpose: In this modelling study we pursued two main goals. The first was to establish a new CTV-to-PTV expansion which considers the closest and most critical organ at risk (OAR). The second goal was to investigate the impact of the planning target volume (PTV) margin size on the tumor control probability (TCP) dependent on the geometrical setup uncertainties. The aim was to achieve a smaller margin expansion close to the OAR while allowing a moderately larger expansion in less critical areas further away from the OAR and most importantly while maintaining the TCP. Methods and Materials: Imaging data of radiation therapy (RT) plans from pet dogs which had undergone RT for brain tumor were used to estimate the clinic specific rotational setup uncertainties. A Monte-Carlo methodology used to quantify the implications of those uncertainties on the TCP. A computational CTV-to-PTV expansion method based on probability density was established. All required software modules were developed and integrated into a software package which directly interacts with the Varian Eclipse treatment planning system. Results: The rotational setup uncertainties were obtained from 44 dog patients. Their frequency distributions were analysed, and the mean values and standard deviations were determined. To investigate the impact on the TCP, several uniform and non-isotropic PTV targets were created. Standardized RT plans with equal optimization constraints were defined and automatically applied and calculated on these targets. The resulting TCPs were then evaluated and compared. Conclusion: The non-isotropic margins were found to result in larger TCPs with smaller margin excess volume. Further we have presented an additional application of the newly established CTV-to-PTV expansion method for radiation therapy of the spinal axis of human patients.


INTRODUCTION
The success and efficacy of radiotherapy are in part determined by the dose distribution inflicted upon the tumor target. In classical fractionated radiotherapy (without radiosensitizer adjuvants, hypothermia, etc.) the spatio-temporal dose distribution is the only controllable parameter. The precision of delivery of the planned and calculated dose to a patient is heavily dependent on the tumor location and the equipment and techniques used. Those determine the precision limiting factors such as setup uncertainties, intra-fractional movement and inter-fractional positioning discrepancies. To ensure adequate dose coverage of the region to be treated, the clinical target volume (CTV) which encompasses the macroscopic and microscopic tumor burden, is further expanded to a planning target volume (PTV) to compensate the aforementioned uncertainties. The resulting PTV often lies within close proximity or even overlaps with organs at risk (OAR), thus research has concentrated on the reduction of the margins in order to reduce the risk of toxicity while maintaining tumor control.
In principle all positional uncertainties could be corrected for in an ideal radiation therapy setting, using a linear accelerator and imager with perfect precision and a robotic couch. In such a scenario adding a margin to the CTV would be unnecessary. However, in realistic setups, especially without robotic couch an adequate margin remains crucial. The conventional approach to PTV definition is by uniform expansion of the CTV by a constant margin as proposed by ICRU [1]. Numerous studies investigated the impact of geometric uncertainties on the tumor control and proposed margin recipes and the direct incorporation thereof into the treatment planning optimization instead of relying on the uniform PTV margin concept [2][3][4][5][6][7][8][9][10][11]. In particular Selvaraj et al. [11] have investigated the impact of systematic and random geometric translational uncertainties on the tumor control depending on the margin size and fractionation. This was done for an artificial spherical CTV target with artificial spherical and four field brick dose distributions. 3 In this in-silico modelling study, our goals were to establish a new CTV-to-PTV expansion which considers the closest and most critical OAR, as well as to investigate as close to clinical reality as possible, the impact of the PTV margin size on the tumor control probability (TCP) dependent on the the geometrical setup uncertainties. The aim was to achieve a smaller margin expansion close to the OAR while allowing a moderately larger expansion in less critical areas further away from the OAR and most importantly while maintaining the TCP. A region of smallest margin (ROSMA) close to the most critical OAR was defined in an area where the PTV was desired to be smallest. A Monte-Carlo based non-isotropic PTV margin generation method, which is similar to the methodology proposed by Stroom et al. [2] was devised.
For quantification of the tumor control probability of the inhomogeneous dose distributions of the treatment plans we used a Monte-Carlo TCP calculation methodology which is reliant on the voxel-based TCP model proposed by Webb and Nahum [12], Radonic et al. [13]. The model was extended to incorporate inter-fractional temporal variations of the dose distribution. An IMRT plan template was applied, optimized and calculated for the different uniform and non-isotropic PTV targets.
Subsequently, the resulting TCPs were evaluated and compared. The non-isotropic margins were found to result in larger TCPs with smaller margin excess volume.

METHODS AND MATERIALS
In fractionated radiotherapy before each fraction the patient has to be positioned appropriately. For each fraction the positioning slightly varies. The variation magnitudes and the ability to correct them depends on the specific tumor site treated and on the equipment and procedures which are used. To account for the setup uncertainties, it is an established procedure to expand the CTV by a safety margin, which yields the PTV. The following section describes our methodology of examining the magnitude of the setup uncertainties occurring in our specific setup, its impact on 4 the treatment outcome and the establishment of a CTV-to-PTV generation method.

Patient inclusion criteria
In this retrospective modelling study imaging data sets from pet dogs, which had undergone radiation therapy for brain tumor at the Vetsuisse Faculty of the University of Zurich in the period from 2007 to 2019 were used. Inclusion criteria were a simulation CT dataset as reference image and at least four cone-beam computed tomography (CBCT) datasets used for position verification during the course of radiation therapy.

Positioning and verification
All dogs underwent a short general anesthesia with endotracheal intubation using routine anesthetic protocols. They were positioned in a rigid positioning system consisting of a custom-made maxillary dental mold bite block (President The Original, Putty Soft, Coltène, Whaledent AG, Altstaetten, Switzerland) on a a polycarbonate tray that supported the maxilla and a vacuum cushion (BlueBag BodyFix, Elekta AB, Stockholm, Sweden) that supported the thorax and front legs. Position was verified with on-board imaging with kilovolt-(kV) CBCT using bone match (Clinac iX with OBI, Varian, Palo Alto, California, USA). In our setup, which uses a 4-degree of freedom couch, translational variations as well as the rotational yaw variation can be and are corrected for. However, roll and pitch variations can not be corrected for with the setup used.

Definition of volumes
To achieve a smaller margin expansion close to the OAR while allowing a moderately larger expansion in less critical areas further away from the OAR, a new structure called region of smallest margin (ROSMA) was created in the contouring workspace. The ROSMA is essentially a singular point with respect to which the daily CBCT to planning CT registrations is performed most carefully. For practical purposes however, a finite size structure is necessary. The ROSMA was delineated depending on the most critical organ at risk (OAR) nearby. It can be the OAR itself or a part thereof. In case of the optic chiasm as most critical OAR, the optic chiasm structure was duplicated and used as ROSMA. In case of the brainstem, brain or eye as most critical OAR, the 3D brush tool with a 0.5cm diameter was used to draw a sphere. The latter was placed in the region where the PTV margin was desired to be smallest; i.e. between the clinical target volume (CTV) and the brainstem, brain or ocular bulb, respectively. In general the ROSMA has to be manually defined by the radiation oncologist as performed in the present study.

Registration of CBCTs to simulation CT
Bone registration of the skull was performed in the image registration workspace of the External Beam Planning system (Eclipse™ Planning system, Varian Oncology Systems, Palo Alto, California, USA). Each CBCT was retrospectively registered to the planning CT using the auto matching tool with bone intensity range from 200-1700 Hounsfield units and limiting the region of interest to the skull. Each registration corresponds to a transformation matrix T (1), which contains the rotations 6 (roll φ, pitch θ, yaw ψ) and translations (x t , y t , z t ).
cos ψ cos θ cos ψ sin θ sin φ − sin ψ cos φ cos ψ sin θ cos φ + sin ψ sin φ x t sin ψ cos θ sin ψ sin θ sin φ + cos ψ cos φ sin ψ sin θ cos φ − cos ψ sin φ y t −sinθ cos θ sin φ cos θ cos φ z t The desired roll and pitch rotation angles were read out from the matrix. This is illustrated for the roll angle in Figure 1. The registrations were done for a total of 44 dogs and 220 CBCTs. Google Protobufs was used to create the data interface for the three dimensional dose distribution and the region of interest (ROI) structures.

Voxelization
The contoured structures such as the CTV, PTV and OARs in DICOM format, as retrieved via the Varian ESAPI, are basically sets of slices. Each slice holds an ordered sequence of points (x,y,z -triplets) defining a contour (cite DICOM standard). For the required processing the CTV structures need to be voxelized.
The grid resolution (voxel size) is chosen such that it matches the CT resolution, with which the required grid size is determined by the contour points with the minimal and maximal coordinates present in the structure. Accordingly, a mapping 8 between the points coordinates r = (x, y, z) T ∈ R and grid indices IDX = (i, j, k) T ∈ Z + and vice versa (3) is implied by (4) and (5).
where RES = (x res , y res , z res ) T is the grid resolution and S is the coordinate transformation which translates all coordinates into positive domain and is given by For all slices, iterating through the ordered sequence of contour points and using the Bresenham's line algorithm [14] the values lying on the contour polygon are set to one. Subsequently using the even-odd rule [15], the value of each voxel which lies inside the boundaries of the contour polygons, is set to one. Which-with the voxelization process is completed, all voxels belonging to the structure now have the integer value one, all other voxels have the value zero.

Generation of the dynamic margin
First a stationary global three-dimensional voxel grid (integer array is initialized.
As illustrated in Fig where T(r, r ROSMA ) is the translation to the center of mass of the ROSMA, and R(θ, ψ, φ) is the rotation matrix. The resulting coordinates r are mapped to the indices of the global voxel grid IDX GLOB as in Eq. (4). We add the value of the rotated voxel to the corresponding voxel in the global grid.
Thus after repeating the procedure for N times, the global voxel grid contains an occupancy probability distribution. The voxel values indicate how many times the particular voxel was inside the rotated CTV. An occurrence threshold can then be defined and subtracted from all voxels of the global grid. Which-after the voxels with a value of zero or above define our newly generated PTV.  During the procedure we have noticed that repeated optimization runs of the same plan with all the same optimization criteria to a just minimally different target leads to quite different spatial dose distributions. To ensure comparability of the different plans, we have done the plan optimizations of the different targets iteratively. In a first step the plan was optimized and calculated for the smallest designated target T 0 . Next the plan was copied and adjusted to the next larger target T 0 → T 1 .
Then the optimization was run but not from scratch but rather starting from the immediate dose already calculated. The iteration T n → T n+1 was continued up until the largest target. For the anisotropic margin targets, the plan of the uniform margin target, with the most similar volume to the designated anisotropic target, was used as the starting point.

Estimation of tumor control probability
Fractionated radiotherapy is simulated in a Monte-Carlo procedure. The TCP as modelled by Nahum and Tait [18] is given by where N S is the number of surviving clonogenic cells. Webb and Nahum [12] have de- where S is the survival model function and D f i is the dose [Gy] which the voxel i is exposed to at fraction f . In case of LQ Model [19] the cell survival function is given by The patient survival (TCP) after a follow-up period τ is then given by where e γT accounts for the effective tumor-cell repopulation rate and e aτ characterizes exponential dependence of the survival rate on the elapse time. The radio-biological parameters for animal brain tumors, used for the TCP calculation, In a Monte Carlo procedure the calculation is repeated K times, this simulates a fictitious patient population undergoing the exact same radiation treatment protocol for the exact same tumor. For each simulation run the resulting TCP value is recorded. This yields a TCP frequency distribution.
We integrated the frequency distribution and created a cumulative histogram.
It represents the percentage of the fictitious simulated patient population having a TCP greater than or equal to the value in the corresponding TCP bin. We have chosen to present the T CP 95% , which means that in 95% of the simulated treatments, the realised TCP was at least equal to the specified T CP 95% value.

Setup uncertainties
In Figure   The rotational roll and pitch uncertainties investigated were assumed to follow normal distributions N (µ, σ 2 ) around µ = 0. We have considered uncertainty magnitudes with σ θ = σ φ = 0.1, as well as σ θ = σ φ = 0.02 rad. This is close to the observed uncertainties for our specific setup as shown in the previous section.
In a first step, for each of the CTVs a set of PTVs with uniform margins of different margin size was automatically generated using TPS provided functionality. The margin sizes considered were specifically 0.05 cm, 0.  Figure 7: The plots show the comparison of the T CP 95% plotted versus the excess volume for the different uniform and anisotropic margin targets for the simulated single fraction and ten fractions treatments with assumed normally distributed rotational shifts with µ θ = µ φ = 0 rad and σ θ = σ φ = 0.  Figure 8: The plots show the comparison of the T CP 95% plotted versus the excess volume for the different uniform and anisotropic margin targets for the simulated single fraction and ten fractions treatments with assumed normally distributed rotational shifts with µ θ = µ φ = 0 rad and σ θ = σ φ = 0.02 rad for the artificial CTV section 2.2. For each plan, respectively for its resulting dose distribution the TCP was simulated in a Monte Carlo procedure as described in section 2.3. This was done for a 10×4.0 Gy fractionation schedule, which is how the patient was originally treated, as well as for a single fraction of 1 × 17.023 Gy. The single fraction dose was calculated to radio-biologically match the ten fractions treatment. In Fig. 7, for both CTVs the T CP 95% is plotted versus the excess volume for σ θ = σ φ = 0.1 rad. The excess volume V excess stands for the additional volume which is added to the CTV volume by the particular PTV target. Figure 8 is the same depiction for σ θ = σ φ = 0.02 rad. In the plots also the theo-19 retical maximum TCP value is shown as a horizontal line. That is the TCP which would result if each voxel of the CTV would always be exposed to the prescribed dose.

DVH and NTCP
To demonstrate the benefit of smaller target excess volume on the dose and NTCP of the OAR, which in our case is the brainstem, we plotted its cumulative dose volume histograms (DVH) (in Fig. 9) as well as computed and compared its NTCP. This was done for a uniform margin target and an anisotropic margin target which have approximately the same T CP 95% . Looking at the plots in Fig. 7 it can be observed that for the original CTV, such a case is given for the uniform 2mm margin target and the anisotropic target with 95% threshold, which result in a brainstem NTCP of 12.4% and 10.2% respectively. For the artificial CTV, this is the case for the 5 mm uniform and 99% threshold anisotropic targets, yielding brainstem NTCPs of 24.5% and 18.2% as well as for the 4 mm uniform and 95% anisotropic targets, which yield 22.7% and 17.2%. For the NTCP calculation parameters from [22] were used. The plan template used was created with the intent to ensure homogeneous dose coverage of the target and optimal conformity. Sparing of normal tissue or critical organs was not an objective, thus the absolute NTCP values are not relevant.
The uniform and anisotropic margins were systematically compared for their impact on the resulting target volume and TCP. Specifically analysed was the correlation of the excess volume, with the TCP which is at least achieved in 95% of the simulated treatments. The excess volume refers to the volume which a target adds to the primary CTV. We presume that the excess volume can be, within certain limits, regarded as proportional to the damage inflicted to normal tissue. Fig. 9 shows that the irradiated volume of the OAR is reduced and therewith the complication probability is smaller. Thus establishing a margin with a small excess volume which achieves the maximum possible TCP is desirable.
The absolute values of the rotational setup uncertainties occurring in our particular setup were found to be very small. Analysis of the histograms show that roll and pitch angles conform to a skew normal distribution around a mean angle which is very close to zero. This supports our hypothesis to assume the rotational uncertainties to be normally distributed around a mean equal zero. The currently used PTV generation procedure which consists of adding uniform margins of 1 - In setups without anaesthesia, which do not allow for such rigid positioning, such as in human patients, noticeably larger rotational uncertainties are to be expected.
In that case the anisotropic PTV was shown to be clearly superior to the uniform margin PTV recipe. It outperforms the uniform margin PTV by providing a larger TCP with a smaller excess volume.
As expected, the geometrical shape of the target and its position relatively to the position of the ROSMA was shown to have a major influence on the resulting TCP -margin dependence. For the artificial target with cylindrical shape the rotational uncertainties have a significantly larger impact on the TCP. Furthermore, the benefit of the anisotropic margin targets over the uniform margin targets is greater.
Another crucial aspect of this investigation is the fractionation. For the single fraction case, stochastically occurring setup uncertainties become systematic uncertainties. Fractionated treatment on the other hand can to some extent compensate for the stochastic uncertainties. This was also observed in the present study.
The stochastic CTV to PTV generation method sometimes leads to artifacts in the generated target, which require manual postprocessing or repetition of the generation process. The improvement and refinement of the procedure could be the subject of future work.
A shortcoming of the treatment simulation is that in reality the rotation would apply to the entire body of the patient, while in the simulation only the CTV is rotated inside the static dose distribution. Thus, the dose distribution in reality might be minimally different from the static dose distribution calculated by the TPS due to the slightly shifted passage of the beams. We believe it is reasonable to assume that this discrepancy is negligible for the rotational uncertainties of magnitudes investigated in this study.
The methodology and the tools established for conducting this study can be 23 easily extended to also investigate translational uncertainties and/or to simulate other relevant aspects such as inhomogeneous clonogenic cell density distribution within the CTV. This could be a potential subject of future research work. For the PTV generation we assumed that the rotation point lies at the cranial end of the spinal coord (around z = 0). We generated the PTVs for rotational variations around the yaw axis with σ φ =0.01 rad and σ φ = 0.02 rad. We assume that these magnitudes approximately match the variations encountered in clinical reality. Figure A1 shows the beam eye view of the whole spinal axis. The purple contour delineates the CTV. The dark blue and the yellow contours depict the PTVs for σ φ =0.01 rad for thresholds of 95% and 98% respectively. The red and the black contours are the PTV for σ φ =0.02 rad for thresholds of 95% and 98% respectively. Figure A2 shows the transversal view at the bottom end of the spinal cord (far away from the rotation point) with its latero-lateral margin expansion.
In Figure A3 we have plotted the lateral margin dependent on the distance from the cranial isocenter z. It can be observed that the margin is a linear function of the z. This allows a simple procedure for dynamic margin definition for spinal Figure A1: Beam eye view of the spinal cord CTV and generated PTVs at the bottom end of the spinal cord. The dark blue and the yellow contours are the PTV for σ φ =0.01 rad for thresholds of 95% and 98% respectively. The red and the black contours are the PTV for σ φ =0.02 rad for thresholds of 95% and 98% respectively. Figure A2: Transverse view of the spinal cord CTV and generated PTVs at the bottom end of the spinal cord. The dark blue and the yellow contours are the PTV for σ φ =0.01 rad for thresholds of 95% and 98% respectively. The red and the black contours are the PTV for σ φ =0.02 rad for thresholds of 95% and 98%, respectively.
axis irradiation in clinical routine. In the plane of the most cranial isocenter a (small) margin which accounts only for the non-rotational uncertainties is drawn.
At the most caudal plane of the spinal cord the lateral margin required is calculated and drawn using the slope. Then a simple interpolation is performed in between.
This functionality is provided by the TPS. As an example: the lateral expansion margin needed at a distance of 40 cm from the cranial isocenter is 1.5 cm (assuming Spinal cord irradiation margin = 0.01 rad 98% threshold = 0.01 rad 95% threshold = 0.02 rad 98% threshold = 0.02 rad 95% threshold slope: m=-0.021 slope: m=-0.017 slope: m=-0.041 slope: m=-0.033 Figure A3: Lateral margin between generated PTVs and the spinal cord CTV dependent on the distance from cranial isocenter z 28