A concept for anisotropic PTV margins including rotational setup uncertainties and its impact on the tumor control probability in canine brain tumors

2.1.1. Patient inclusion criteria

2.1.2. Positioning and verification

2.1.3. Definition of ROSMA

2.1.4. 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 a bone intensity range from 200-1700 Hounsfield units and limiting the region of interest to the skull. Each registration corresponds to a transformation matrix, which contains the rotations (roll ϕ, pitch θ, yaw ψ) and translations (x_t , y_t , z_t ). The desired roll and pitch rotation angles were read out from the transformation matrix as illustrated for the roll angle in figure 2. The registrations were done for a total of 44 dogs and 220 CBCTs.

Zoom In Zoom Out Reset image size

From the obtained set of roll and pitch variation angles, a Gaussian probability distribution was calculated for both rotation axes. The probability distributions were then used for Monte Carlo sampling in the PTV generation and the TCP estimation procedures, as will be explained in section 2.2.

For the selected patient, the original CTV of the irradiated brain tumor was used. With the intent to show that the generated margin is especially useful for constellations where the CTV has a cylindrical shape and the associated ROSMA lies on or close to its longitudinal axis, a second artificial CTV with similar volume, but with a cylindrical shape was drawn in the same region as the original CTV (see figure 3). In a first step, the CTV structure was voxelized. This is described in detail in section C.1 of the appendix. Figure 4 illustrates how rotation angles are sampled from the probability distributions derived from the CBCT registrations in a Monte Carlo procedure. Rotations (using the sampled rotation angles) are then applied to the voxelized CTV structure subsequently generating a three dimensional probability map. The PTV is then created by subtracting a threshold from the probability map and extracting corresponding contours. Specifically, 10'000 rotation angle sets were sampled and as many rotations were applied to the voxelized CTV. In order to enable the generation of a margin with 99.95% occurrence threshold using our method, 10'000 rotations are required. The procedure is explained in detail in section C.2 of appendix C. An additional application of the generation procedure for humans is presented in appendix A.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For RT planning the Varian Eclipse Treatment Planning System (TPS) was used. Varian provides an C#/.NET based API called ESAPI (Eclipse Scripting Application Programmer Inferface). The API provides methods for data access and manipulation and allows the automatization of various aspects of treatment planning. The computationally intensive procedures were implemented in C++ and OpenMP was used to parallelize parts of the code. To interface the C++ code with the C# code a wrapping layer was created using Microsoft C++/CLR. Google Protobufs [14] was used to create the data interface for the three dimensional dose distribution and the region of interest (ROI) structures.

A nine field equispaced coplanar intensity modulated radiation therapy (IMRT) plan was devised. A minimal set of optimization criteria was chosen such that a homogeneous coverage of the designated target (PTV) and optimal conformity were ensured. Utilizing ESAPI, a script was created which automatically adjusts the devised plan template to each designated target (PTV) and subsequently runs the plan optimization, dose calculation and normalization. For the plan optimization, the Varian Photon Optimizer Version 15.06.04 was used. The calculation model used was the Varian Anisotropic Analytical algorithm Version 15.06.04. The dose calculation grid precision was set to the maximal precision available which was one millimetre. For plan normalization a helper structure in the center of the CTV was defined by subtracting a uniform margin ${r}_{m}=\tfrac{1}{2}\cdot \sqrt[3]{\tfrac{3}{4\pi }{V}_{\mathrm{CTV}}}$ from the CTV, for V_CTV being the CTV volume. The plan was then normalized, such that the median dose inside the devised helper structure equalled the prescribed dose.

During the procedure, we noticed that repeated optimization runs of the same plan with all the same optimization criteria to a just minimally different target can lead to different spatial dose distributions. To ensure comparability of the different plans, the plan optimizations of the different targets were done iteratively, as is shown figure 5). In a first step, the plan was optimized and calculated for the smallest designated target T₀. Next, the plan was copied and adjusted to the next larger target T₀ → T₁. As a next step, the optimization was run 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 uniform margin target plans, with the most similar volume to the designated anisotropic target, was used as a starting point. For example, in our case T₀ was the 0.5 mm uniform margin target, which was used as starting point for plan optimization of the 1 mm uniform margin target.

Zoom In Zoom Out Reset image size

Fractionated radiotherapy was simulated in a Monte-Carlo procedure. The TCP as modelled by Nahum and Tait [15] is given by

$$\begin{eqnarray}&&{TCP}={e}^{-{N}_{S}}\end{eqnarray} \tag{ 1 }$$

where N_S is the number of surviving clonogenic cells. Webb and Nahum [12] derived a model for the TCP with non-uniform clonogenic cell density and non-uniform dose. We extended the model to incorporate fractions and the temporal variations of the dose distributions, which are thereby implicated. The extended model is used in the Monte-Carlo procedure.

Rotation angles are sampled from probability distributions obtained from CBCT registrations, as described in section 2.1.4. As illustrated in figure 6, the voxelized CTV (see section C.1 in the appendix) is rotated around the center of mass of the delineated ROSMA inside the dose distribution of the RT plan. Each rotation corresponds to a treatment fraction. For each voxel of the rotated CTV the dose exposure for the current fraction is recorded. Therewith the characteristic of the fractionated radiotherapy is modelled, that due to geometric positioning uncertainties, a particular CTV voxel gets a different dose at each fraction.

Zoom In Zoom Out Reset image size

The recorded dose values are used to calculate the cell survival of a particular CTV voxel. The cell survival at fraction f of voxel i is given by

$$\begin{eqnarray}&&{S}_{i}^{f}=S({d}_{i}^{f},\alpha ,\beta ),\end{eqnarray} \tag{ 2 }$$

where S is the survival model function and d_i ^f is the dose in units of Gy, which the voxel i is exposed to at fraction f. In case of the LQ Model [16], the cell survival function is given by

$$\begin{eqnarray}&&S({d}_{i}^{f},\alpha ,\beta )={e}^{-\alpha {d}_{i}^{f}-\beta {\left({d}_{i}^{f}\right)}^{2}}\end{eqnarray} \tag{ 3 }$$

where α and β characterize the intrinsic radio-sensitivity of the cells. In the work of Dale [17] the LQ-model was combined with an assumed tumor repopulation factor. Further, Qi et al [18] characterized the dependence of the survival rate on the elapse time by an exponential decrease. The cell survival in voxel i after M fractions is

$$\begin{eqnarray}&&{S}_{i}=\prod _{f=0}^{M}{S}_{i}^{f},\end{eqnarray} \tag{ 4 }$$

where S_i ^f is the survival after a single fraction f. The patient survival after a follow-up period τ is then described by

$$\begin{eqnarray}&&{TCP}=\prod _{i=0}^{N}{e}^{-{\rho }_{i}{V}_{i}{S}_{i}{e}^{\gamma T}{e}^{a\tau }}\end{eqnarray} \tag{ 5 }$$

where V_i is the volume and ρ_i is the cell density of a voxel i. For this investigation a homogeneous cell density ρ = 10⁷ cells/cm³ was assumed inside the CTV, thus ∀ i: ρ_i = ρ. The term e^{γ T} accounts for the effective tumor-cell repopulation rate, where $\gamma =\mathrm{ln}(2)/{T}_{D}$, T_D being the tumor doubling time and T being the treatment time. For the treatment time T in days, seven-fifth times the number of fractions M was used if M > 5, else it was assumed to be equal to the number of fractions. e^{a τ} characterizes the exponential dependence of the survival rate on the follow-up time τ. The radio-biological parameters for canine brain tumors, used for the TCP calculation, obtained from Radonic et al [13], were α = 0.36 Gy⁻¹, α/β = 8 Gy, a = 0.9 yr⁻¹, T_D = 5.0 d. Clinical canine patient data such as survival depend on the follow-up time. The TCP model used in this work was fitted to clinical canine patient survival data, thus it incorporates the follow-up as well [13]. In this work the TCPs were calculated for a follow-up time of τ = 2 years.

In a Monte Carlo procedure, the calculation is repeated K times, which simulates a fictitious patient population undergoing the exact same radiation treatment protocol for an identical tumor. For each simulation run the resulting TCP value was recorded. This yields a TCP frequency distribution. We used K = 10000 for simulating single fraction RT schedules and K = 1000 for simulating RT schedules with 10 fractions.

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 chos to present the TCP_95%, which means that in 95% of the simulated treatments, the realised TCP was at least equal to the specified TCP_95% value.

As plan conformity might impact the TCP, in order to get a sense of the conformity of the resulting plans, the conformity index (CI) as defined in Feuvret et al [19]

$$\begin{eqnarray}&&{CI}=\displaystyle \frac{{V}_{{RI}}}{{TV}}\end{eqnarray} \tag{ 6 }$$

was calculated for all generated plans on the anisotropic and uniform margin PTVs for the original CTV. Where V_RI is the volume of the reference isodose, here we have specifically chosen 95% (of the prescribed dose). TV is the target volume.

We investigated the impact of rotational uncertainties on the TCP depending on the margin size, for uniform margins, as well as for anisotropic margins. For this purpose two CTVs with different geometrical shapes were selected. The first CTV was the original brain tumor CTV, which was actually irradiated. The original CTV resembles an ellipsoid with some distortions. With the intent of investigating the influence of the geometrical shape, a second CTV with a similar volume was artificially drawn in the same region as the original CTV, resembling a more cylindrical longitudinal shape. The CTVs and the ROSMA are depicted in figure 3.

The roll and pitch uncertainties investigated were assumed to follow normal distributions ${ \mathcal N }(\mu ,{\sigma }^{2})$ centred around μ = 0. We considered uncertainty magnitudes with σ_θ = σ_ϕ = 0.1 rad (5.73°), as well as σ_θ = σ_ϕ = 0.02 rad (1.15°). The latter value is close to the observed uncertainties for our specific setup as will be shown in the corresponding section 3.1 in the results.

In a first step, for both CTVs a set of PTVs with uniform margins of different margin sizes was automatically generated using a functionality provided by the TPS. The margin sizes considered were 0.5 mm, 1 mm, 2 mm, 3 mm, 4 mm, 5 mm and 6 mm. Next, sets of anisotropic PTV were generated for different occurrence thresholds, namely 50%, 90%, 95%, 99%, 99.5%, 99.8%, 99.9% and 99.95% for both CTVs and both uncertainty distributions were considered. For clarification: an occurrence threshold of 90% means that from the probability map generated with 10 000 rotations, a threshold of 1000 is subtracted. The anisotropic margin generation procedure is described in more detail in section C.2 in the appendix. For each computed target, a RT plan was generated as described in section 2.4., The TCP was simulated in a Monte Carlo procedure for the resulting dose distribution for each plan, as described in section 2.5. This was done for a 10 × 4.0 Gy fractionation schedule, which corresponds to original treatment plan of the patient, 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.

The single fraction treatment schedule was deliberately chosen with the intent to show how the fractionation influences the impact of setup uncertainties on the TCP.

In figure 7 the roll and pitch angles obtained from registrations of positioning CBCTs are plotted in histograms. The means of the pitch and roll angles were found to be μ_θ = −0.003 rad (0.17°) and μ_ϕ = −0.001 rad (0.06°) while the standard deviations were σ_θ = 0.011 rad (0.63°) and σ_ϕ = 0.017 rad (0.97°).

Zoom In Zoom Out Reset image size

In figure 8, the TCP_95% is plotted for both CTVs as a function of the excess volume for σ_θ = σ_ϕ = 0.1 rad (5.73°). The excess volume V_excess stands for the additional volume which is added to the CTV volume by the particular PTV.

$$\begin{eqnarray}&&{V}_{\mathrm{excess}}={V}_{{PTV}}-{V}_{{CTV}},\end{eqnarray} \tag{ 7 }$$

where V_PTV is the volume of the PTV and V_CTV is the volume of the CTV. Figure 9 shows the same for σ_θ = σ_ϕ = 0.02 rad (1.15°). The theoretical maximum TCP value is the TCP which would result if each voxel of the CTV would always be exposed to the prescribed dose.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In tables 1 and 2, the conformity index as well as the ratio of the excess volume to the volume of the CTV values are provided for the plans computed on the generated uniform and anisotropic PTVs associated with the orginal CTV respectively. For the plans on all PTVs, the conformity index was found to be between 1 and 1.33, which according to Feuvret et al [19] can be interpreted as compliant with the treatment plan.

Table 1. Uniform margin PTVs of the original CTV for different margin sizes: Conformity index (CI) of resulting plans for reference isodose (RI) = 95%; Ratio of excess to CTV volume: $\tfrac{{V}_{{PTV}}-{VCTV}}{{V}_{{CTV}}}$.

Uniform margin [mm]	CI RI = 95%
		$\tfrac{{V}_{\mathrm{excess}}}{{V}_{\mathrm{CTV}}}$
0	1.33	0
0.5	1.23	0.18
1	1.20	0.32
2	1.15	0.71
3	1.10	1.14
4	1.09	1.61
5	1.05	2.15

Table 2. Anisotropic margin PTVs of the original CTV for different occurrence thresholds: Conformity index (CI) of resulting plans for reference isodose (RI) = 95%; Ratio of excess to CTV volume: $\tfrac{{V}_{{PTV}}-{VCTV}}{{V}_{{CTV}}}$.

Threshold [%]	CI RI = 95%
		$\tfrac{{V}_{\mathrm{excess}}}{{V}_{\mathrm{CTV}}}$
50.00	1.29	0.05
90.00	1.22	0.35
95.00	1.21	0.45
99.00	1.22	0.66
99.50	1.20	0.75
99.80	1.20	0.85
99.90	1.19	0.93
99.95	1.18	0.99

The contoured structures such as the CTV, PTV and OARs in DICOM format, as retrieved via the Varian ESAPI, are sets of slices. Each slice holds an ordered sequence of points (x,y,z—triplets) defining a contour. 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

$$\begin{eqnarray}&&m\,:(x,y,z)\to (i,j,k)\end{eqnarray} \tag{ C1 }$$

$$\begin{eqnarray}&&{m}^{-1}\,:(i,j,k)\to (x,y,z)\end{eqnarray} \tag{ C2 }$$

between the points coordinates ${\bf{r}}={(x,y,z)}^{T}\in {\mathbb{R}}$ and grid indices ${\bf{IDX}}={\left(i,j,k\right)}^{T}\in {{\mathbb{Z}}}^{+}$. The mapping m is implied by:

$$\begin{eqnarray}&&{\bf{IDX}}=m(x,y,z)=\lfloor ({\bf{r}}+{\bf{S}})/{\bf{RES}}\rfloor ,\end{eqnarray} \tag{ C3 }$$

$$\begin{eqnarray}&&{\bf{r}}={m}^{-1}(i,j,k)={\bf{RES}}\odot {\bf{IDX}}-{\bf{S}}+\displaystyle \frac{{\bf{RES}}}{2},\end{eqnarray} \tag{ C4 }$$

where ${\bf{RES}}={({x}_{{res}},{y}_{{res}},{z}_{{res}})}^{T}$ is the grid resolution and S is the coordinate transformation which translates all coordinates into the positive domain and is given by

$$\begin{eqnarray}&&{\bf{S}}=-\min \{x,y,z\in \mathrm{CTV}\}.\end{eqnarray} \tag{ C5 }$$

For all slices, iterating through the ordered sequence of contour points and using the Bresenham's line algorithm [21] the values lying on the contour polygon are set to one. Subsequently using the even–odd rule [22], the value of each voxel which lies inside the boundaries of the contour polygons, is set to one. With that, the voxelization process is completed, all voxels belonging to the structure now have the integer value one, all other voxels have the value zero.

First a stationary global three-dimensional voxel grid (integer array) is initialized. As illustrated in figure 4 in a Monte Carlo procedure the rotation angles are sampled from the probability distributions derived from the CBCT registrations. Using the voxelized CTV structure obtained in the voxelization procedure, we iterate over all voxels. Using the sampled rotation angles (θ, ψ, ϕ), a rotational transformation is applied to the midpoint coordinates of each voxel W_IDX , rotating it around the center of mass of the delineated ROSMA r_ROSMA

$$\begin{eqnarray}&&{\bf{r}}^{\prime} ={{\bf{T}}}^{-1}({\bf{r}},{{\bf{r}}}_{{ROSMA}}){\bf{R}}(\theta ,\psi ,\phi ){\bf{T}}({\bf{r}},{{\bf{r}}}_{{ROSMA}})\cdot {\bf{r}}\end{eqnarray} \tag{ C6 }$$

where T(r, r_ROSMA ) is the translation to the center of mass of the ROSMA, and R(θ, ψ, ϕ) is the rotation matrix. The translation only pertains to the clinical matching process. The resulting coordinates ${\bf{r}}^{\prime}$ are mapped to the indices of the global voxel grid IDX_GLOB as described in equation (C3). We add the value of the rotated voxel to the corresponding voxel in the global grid.

$$\begin{eqnarray}&&{W}_{{{\bf{IDX}}}_{\mathrm{GLOB}}}^{\mathrm{GLOB}}\,+={W}_{{\bf{IDX}}}\end{eqnarray} \tag{ C7 }$$

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 ${TR}\lt N$ can then be defined and subtracted from all voxels of the global grid.

$$\begin{eqnarray}&&\forall \,i,j,k:{W}_{i,j,k}\,-={TR}\end{eqnarray} \tag{ C8 }$$

Now, the voxels with a value of zero or above define our newly generated PTV.

To write the generated PTV structure back to the TPS, the corresponding contour points have to be extracted from the voxel representation. The outermost voxels where W_i,j ≥ 0 are found, by iterating through all voxels W_i,j for each slice k. The midpoints of those voxels define a contour, as can be seen in figure C1. However the points are not yet ordered in a sequence as required by the TPS. To create an ordered polygon from the unordered set of points a concave hulling algorithm [23, 24] is used. The obtained polygons are written back into the TPS as structure ${{PTV}}_{\mathrm{GEN}}^{{\prime} }$. Subsequently a boolean OR operation between the segment volumes of the generated structure and the CTV is performed such, that it is guaranteed for the CTV to be fully included in PTV_GEN

$$\begin{eqnarray*}&&{CTV}\subset {{PTV}}_{\mathrm{GEN}}.\end{eqnarray*}$$

This is done to account for the possibility that for a very small threshold some CTV voxels are outside of the generated PTV margin. That concludes the generation procedure of the PTVs.

Zoom In Zoom Out Reset image size