Dose mapping sensitivity to deformable registration uncertainties in fractionated radiotherapy – applied to prostate proton treatments

Background Calculation of accumulated dose in fractionated radiotherapy based on spatial mapping of the dose points generally requires deformable image registration (DIR). The accuracy of the accumulated dose thus depends heavily on the DIR quality. This motivates investigations of how the registration uncertainty influences dose planning objectives and treatment outcome predictions. A framework was developed where the dose mapping can be associated with a variable known uncertainty to simulate the DIR uncertainties in a clinical workflow. The framework enabled us to study the dependence of dose planning metrics, and the predicted treatment outcome, on the DIR uncertainty. The additional planning margin needed to compensate for the dose mapping uncertainties can also be determined. We applied the simulation framework to a hypofractionated proton treatment of the prostate using two different scanning beam spot sizes to also study the dose mapping sensitivity to penumbra widths. Results The planning parameter most sensitive to the DIR uncertainty was found to be the target D95. We found that the registration mean absolute error needs to be ≤0.20 cm to obtain an uncertainty better than 3% of the calculated D95 for intermediate sized penumbras. Use of larger margins in constructing PTV from CTV relaxed the registration uncertainty requirements to the cost of increased dose burdens to the surrounding organs at risk. Conclusions The DIR uncertainty requirements should be considered in an adaptive radiotherapy workflow since this uncertainty can have significant impact on the accumulated dose. The simulation framework enabled quantification of the accuracy requirement for DIR algorithms to provide satisfactory clinical accuracy in the accumulated dose.


Background
The patient geometry can change significantly between fractions during radiotherapy treatments. This motivates the use of image sets to monitor the treatment progress and to serve as basis for optional re-planning to better fulfil the treatment objectives.
The ideal monitoring scenario would be to score dose for each individual tissue part, or cell, throughout all delivered treatment fractions. The closest realisation of such a scenario is to calculate dose on images taken for each treatment occasion, and use deformable image registration (DIR) to map and accumulate the dose contributions onto a representative, single image dataset [1]. Such dose distributions can be compared with the intended, planned dose distribution to provide a basis for corrective interventions in adaptive radiotherapy. The accuracy of the image registration will govern the accuracy of the mapped dose, and hence the relevance of the corrections. The validity of the accumulated dose, due to the image registration uncertainty, was recently debated by Schultheiss and Tomé [2]. It is therefore of fundamental interest to determine the accuracy requirements of the image registration process, separate from other uncertainties, in particular with regards to the accuracy of the quantities derived from the finally mapped dose.
The input to the registration also contains uncertainties, e.g. image noise and uncertain segmentations. The registration uncertainty can be considered as uncorrelated to other treatment uncertainties such as organ motion effects or dose calculation uncertainties, and can thus be investigated separately. There exists a plethora of DIR algorithms in the literature, see e.g. [3,4] for a review of available methods. Studies of the accuracy of different DIR algorithms applied to different body sites exist [5,6] with reported average errors in the range of 1-5 mm. The accuracy of dose mapping has also been studied explicitly by several authors; Rosu et al. [7] investigated dose grid resolution effects, Salguero et al. [8] utilised the image registration inverse inconsistency and found a maximum dose mapping uncertainty of >30% of the prescription dose for a lung patient, Yan et al. [9] related the lack of mass conservation in the registration to the dose uncertainty, Hub et al. [10] estimated the dose uncertainty from the registration parameter uncertainty and Murphy et al. [11] developed a method to sample image registration errors and demonstrated their effect on the mapped dose.
The aim of this study is to develop a methodology for determining the accuracy requirements for the use of deformable image registration as a basis for dose accumulation in fractionated radiotherapy. For this purpose we developed a simulation framework for dose accumulation where the uncertainty of the calculated treatment outcome (represented as dose, or as dose response determined with a biological models) as a function of the registration uncertainty could be investigated. The dosimetric impact of the registration uncertainty on the entire fractionated treatment has been studied by Risholm et al. [12] who presented the uncertainty of the total physical dose for a single estimated registration uncertainty, while we study how the final dose and the related response would vary with the registration uncertainty.
We apply the framework to a hypofractionated treatment of the prostate with spot scanned protons where a "plan of the day" is assumed to be tailored to the CTV for each individual fraction. The clinical importance of dose accumulation for prostate treatments [13] was investigated in a recent publication by Wen et al. We hypothesise that the sensitivity to dose mapping uncertainty increases with steeper dose gradients and therefore we use dose distributions from two different spot sizes, one based on measurements at the local clinic and one with smaller spot sizes producing a sharper and less forgiving penumbra. We also investigate the accuracy requirement versus the required size of the planning margin used for construction of the planning target volume (PTV) from the clinical target volume (CTV).

Methods
The impact of the registration uncertainty on the accumulated dose, and its estimated response, can be studied by adding known uncertainties to the registrations used in mapping the deformations. The simulations mimic a clinical workflow where a "plan of the day" is tailored to the actual position of the CTV for the individual fractions. The total accumulated effect, i.e. the estimated response for the dose mapped from all fractions to a fixed reference geometry is scored for various degree of registration uncertainties. As we are specifically interested in the effects of the registration uncertainty we assume that the patient is imaged and setup without errors at each fraction. According to the van Herk margin scheme [14] this means that we simulate a situation yielding only the random contribution from DIR based dose mapping to the CTV to PTV margin. In a real clinical situation also other sources of uncertainty must be considered. In section A we describe the simulation framework and the modelling of the image registration uncertainty, and in section B we apply the uncertainty model to a hypofractionated prostate treatment with spot scanned protons.
A. Simulation of a fractionated treatment with dose mapping uncertainties A "plan of the day" is prepared for each of the N treatment fractions, and the fraction dose is mapped to a reference (fixed) image set I F for evaluation of the cumulative radiation effect, see Figure 1. For each frac- The dose distribution d i ð Þ F r ð Þ will vary from fraction to fraction because of organ motion and the "plan of the day" approach. We adopt a view where we with a given uncertainty can track the position of a tissue element at all treatments and accumulate the radiobiological effect fraction by fraction. The D F (r) then becomes a stochastic quantity.
The resulting radiobiological effect was modelled based on the LQ-model [15]. The 4D dose distribution is thereby converted into a 3D equivalent fraction dose d eq F r ð Þ [16], calculated as the fraction dose giving the equivalent total radiobiological effect, S F (r), by solving for d eq F r ð Þ from The total equivalent dose, D eq F r ð Þ ¼ Nd eq F r ð Þ was used for analysis of the treatment DVHs. The tumour control probability, TCP, was calculated directly from S F (r), c.f. [17,18]. The normal tissue complication probability, NTCP, was calculated using the Lyman model [19] with  Figure 1 Simulated clinical workflow. The flow chart to the right shows how the clinical workflow, listed to the left, is simulated. The fraction dose is calculated on a randomly deformed virtual patient model for which the exact image registration to the fixed geometry I F is known. A registration uncertainty can be added (the red part of the flow chart), before the dose is mapped to I F where the radiation effect is accumulated. parameters from Burman et al. [20] and Emami et al. [21]. Also the equivalent uniform dose, EUD [22], was calculated from D eq F .

Modelling the image registration uncertainty
The transformation T (i) (r) between the fixed image I F and the moving image I i ð Þ M , used by the dose mapping procedure in equation (1), was modelled through where T are in general known for a clinical situation and thus T (i) have to be treated as a stochastic quantity for which we assume that T i ð Þ e r ð Þ ¼ 0: Registration algorithms apply various regularisation techniques to produce a well behaved and physically realistic transformation. Clearly, there is a local correlation of the error for points close together whereas for points further apart the error will be more uncorrelated. In modelling of T i ð Þ e this is mimicked by sampling uncorrelated errors at a sparse grid of control points which by means of a 3D cubic B-Spline interpolation is applied to the denser dose grid. The registration deviation amplitudes is sampled at the control points, independently for each coordinate according to a multinormal distribution with zero mean and standard deviation σ e . A low resolution of the control points will produce a slowly varying vector field and thus mimic a high degree of regularisation and vice versa for a high resolution of the control points.
Registration algorithms are often validated using landmarks where the distances between known displacements of the landmarks and those calculated by the algorithm, which in our case is equivalent to the absolute registration error T i

.
A single patient instance is defined as a patient with a unique geometry per fraction. In the simulation we used ten different patient instances for evaluation of each parameter combination. A male pelvic virtual patient was for this purpose constructed from average data for 15 prostate patients treated in the supine position, see Figure 2. The prostate was modelled as a sphere with radius r CTV = 2.5 cm located at 20 cm depth and the bladder as an ellipsoid with radii 3.0 cm, 4.0 cm and 2.5 cm in the lateral, superior-inferior and axial directions, respectively, with a concave intrusion from the presence of the prostate. The rectum was modelled by a curved cylinder with outer radius 1.4 cm and length 10.0 cm with wall thickness 0.4 cm [23]. Other anatomical details, such as the femoral heads, were not included since they are in principle uncorrelated with the image registration dose mapping effects for the prostate, bladder and rectum regions.
The prostate is assumed to be incompressible and change location (but not shape) as the rectum and bladder change filling and shape. In prostate radiotherapy treatment it is common to also irradiate the lower parts of the seminal vesicles to the same dose as the prostate. The combined volume, prostate and lower parts of the vesicles, will constitute a rather convex volume which we for simplicity model as a sphere shaped CTV. We sample the displacement of the prostate centre for each fraction, v i ð Þ CTV ; according to a multinormal distribution with the standard deviations 0.4, 0.1 and 0.4 cm in the anterior-posterior, lateral and superior-inferior directions, respectively, consistent with the literature [24]. The fixed geometry I F , into which the radiation effect is accumulated, is simply chosen to be the one with zero displacement of the prostate, cf. Figure 2a.
A simple tissue deformation model was used to construct the true moving geometry I i ð Þ M from I F where the tissue displacement outside the prostate in I F is exponentially relaxed with the squared distance from the prostate edge according to where T 0 is the local displacement at position r. The relaxation parameter k was set to 0.1 to give reasonable volume differences in the rectum. The tissue displacement according to equation (6) for a prostate displace- = 0.5 cm and k = 0.1 is shown in Figure 3.
The virtual patient is shown for two prostate displacements in Figure 2b and c where the rectum and bladder are deformed according to the deformation model in equation (6). The model does not assume anything about the reasons why the prostate has moved, i.e. the motion can be seen as a consequence of the filling in the rectum or the bladder. The dose grid resolution was set to 0.1 cm while two different grid sizes, 1.0 cm and 3.0 cm, was used for the control points in the registration error. Analysis of all T (i) (r), i.e. the patient model and the registration error, showed values of the Jacobian determinant, J(T (i) (r), r), in the interval [0.4, 2.0] for the lower resolution indicating an inverse consistent transform without folds or tears since J > 0. The T (i) (r) produced using the higher control point resolution included small regions with negative values of the J for σ e ≥ 0.5 cm indicating a transform that is not well regularised. It is desirable that the accumulated dose is independent of the choice of reference image I F and this requires that J > 0, which is not in general guaranteed by registration algorithms. The resulting dose mapping uncertainty was very similar for both resolutions but slightly more sensitive with respect to the image registration uncertainty for the 3.0 cm   Figure 3 Tissue displacement model. The absolute value of the tissue displacement, |T 0 |, as a function of the distance |r| from the prostate centre, when a prostate with radius 2.5 cm is assumed to be displaced 0.5 cm. Inside the prostate, i.e. |r| < r CTV , the displacement is constant. Outside the CTV, i.e. |r| > r CTV , the displacement decreases exponentially and vanishes far from the prostate.
resolution. We have therefore concentrated on results for the 3.0 cm resolution.

Generation of treatment plans for simulation
The simulated treatments were assumed to be delivered with two opposed, scanned proton beams, indicated by the arrows in Figure 2a, aiming for a homogenous target dose. The PTV was constructed by adding a margin m isotropically around the CTV.
Clinically relevant planning objectives were chosen from literature values for intermediate risk patients for radiobiological evaluation of the prostate [25], see Table 1, to optimise the spot weights and thus shape the proton dose distribution. Additionally, an artificial OAR was created as a spherical shell around the PTV to suppress the normal tissue dose. The target dose prescription was chosen to give a TCP of 80% since the TCP curve is steep at that dose level and thus the TCP should be susceptible to dose mapping uncertainties. However, because of the suggested low α/β there is a trend towards applying hypofractionated treatments, see e.g. Ritter [26]. Therefore, an aggressive hypofractionation scheme of 6.7 Gy × 5fx still aiming for a TCP of 80% were chosen with the OAR planning objectives from Isacsson et al. [27] scaled with the new prescription dose. Zavgorodni [28] noted that it might be important to take the variable fraction dose into account, as in equation (3), when accumulating the dose for normal tissue and tumours with low α/β, such as the prostate, and for hypofractionated treatments.
The "plan of the day" treatment scenario described above requires a tailored dose distribution for each target position for all patient instances, i.e. d i ð Þ M in equation (1). To save computation time and make the result less dependent of the intrinsic details of the optimiser we displace the dose distribution according to the fraction specific CTV displacement. This results in one optimisation per margin size and spot size regimen.
The scanned proton dose distribution was calculated using an in-house pencil beam algorithm with the depth dose and lateral scattering calculated according to Bortfeld [29] and Russell et al. [30], respectively. Our normal beam model (NBM) parameters are based on data from Kimstrand et al. [31]. Dose mapping in the presence of a sharp penumbra will be sensitive to the registration uncertainty, especially for targets where the loss of dose coverage can greatly affect the outcome. Dose distributions were therefore also generated according to a sharp penumbra beam model (SBM) whose in-air spot sizes are smaller and with zero divergence thus producing a sharper penumbra. The spot sizes of SBM correspond approximately to what modern commercial proton machines can deliver, e.g. see [32]. The penumbra width is often reported as the distance that the dose falls from 80% to 20% of the target dose and here they were 1.16 and 0.88 cm for NBM and SBM respectively in the directions perpendicular to the beam axis. The spot size increase with treatment depth, due to multiple scatter of the protons, and note that a shallower target will have a sharper penumbra and thus more sensitive to dose mapping uncertainties.

Results
The simulations were performed for ten different instances of the virtual patient where each instant had its own set of fraction specific CTV positions with corresponding deformation fields. A single treatment simulation, according to Figure 1, modelled the treatment of a patient instance with the prescribed 6.7 Gy × 5fx. To ensure adequate statistics, each treatment simulation for the same patient instance was repeated 40 times with different sampled image registration T (i) according to equation (5)

The DVH dependency on varying registration uncertainty
The spread of the calculated DVH curves for different registration uncertainties is shown in Figure 4 including the nominal DVH with no registration uncertainty as a reference. As expected, the spread of the DVH curves increases when the registration uncertainty 〈|T e |〉 increases, c.f. the upper left panel with the rest in Figure 4. For the CTV, use of a planning margin (m > 0 cm) reduces the sensitivity to the registration uncertainty, and the spread of the DVH curves decreases compared to planning without margin (m = 0 cm). Hence, an increased target margin reduces the difference in dose between the target and surrounding normal tissue and therefore reduces the effect of the registration uncertainty. A sharper penumbra, i.e. created by the SBM instead of the NBM, gives a slightly increased sensitivity to registration uncertainty.  The spread of the DVHs for the OAR is less dependent of the planning margin and the penumbra width, although the overall dose level increases when any of these two parameters increases. The dose tend to be overestimated in the low dose region but underestimated in high dose regions which is clearly seen for the rectum and bladder in Figure 4. The reason for the underestimation for the high dose is similar as for the target dose, a dose mapping uncertainty is more likely to lower the dose (and vice versa for the low dose region).
The dose spread for the 95% of the CTV volume, D 95 , and the volume of the OARs that receives 80% of the prescription dose, V 80 , is chosen to further illustrate the spread of the DVHs. The frequency distributions of ΔD 95 and ΔV 80 expressed as the change in D 95 and V 80 versus the nominal values are shown in Figure 5 for all sampled patient instances. These results demonstrate that the registration uncertainty causes a decreased value of the calculated D 95 (ΔD 95 ≤ 0%) for the target. It is evident that an increased planning margin will decrease the impact of the dose mapping uncertainty. There is also a trend that V 80 gets underestimated as the dose mapping uncertainty increases.
The absolute values of ΔD 95 and ΔV 80 that include 95% of all the DVH curves, i.e. |ΔD 95 | 95 and |ΔV 80 | 95 are shown in Figure 6 for dose distributions from beam models NMB and SBM. As expected, the |ΔD 95 | 95 is decreased when the planning margin is increased, or when a wider penumbra (NBM) is used.
To obtain the DIR accuracy such that |ΔD 95 | 95 is approximately 3% of the target dose, we found that the registration uncertainty 〈|T e |〉 must be less than 0.25 cm for NBM and SBM with SBM slightly more sensitive to 〈|T e |〉. If a planning margin of 0.1 cm is used, the 〈|T e |〉 registration uncertainty can be relaxed to 0.35 cm for NBM and 0.30 cm for SBM.
The value of |ΔV 80 | 95 was found to increase with increasing 〈|T e |〉 and weakly with increasing margin, as shown in the lower plots in Figure 6. The |ΔV 80 | 95 value was found to increase linearly as a function of 〈|T e |〉, and was approximately 1% and 4% of the ROI volume for 〈|T e |〉 = 0.1 cm and 0.35 cm, respectively. The |ΔV 80 | 95 increases weakly with increased margin and the difference between NBM and SBM is small. The results for the bladder |ΔV 80 | 95 were similar as the data found for the rectum.

The uncertainty in simulated radiation response
The sensitivity of the calculated radiobiological outcome quantities TCP and EUD for the CTV, and NTCP and EUD for the OARs, were determined for each treatment simulation. The chosen technique yields very low nominal NTCP values, less than 2% for the rectum and much smaller for the bladder. Accordingly, the EUD uncertainties for the risk organs were also very small, < 0.5 Gy. The decrease in TCP versus the nominal without dose mapping uncertainty, ΔTCP CTV , was found to increase with increasing registration uncertainty as shown in Figure 7. These results are analogous to the results of the DVH analysis in that the sensitivity to the registration uncertainty is reduced with increasing planning margin, and/or a shallower penumbra. The ΔTCP CTV was below 1% for both beam models NBM and SBM as long as the registration uncertainty 〈|T e |〉 was less than 0.25 cm. The standard deviation of the ΔTCP CTV (not shown) was small (<1%). A similar analysis of the ΔEUD CTV yielded quantitatively very similar results as for the ΔTCP CTV .

The distribution of spatial dosimetric uncertainty with varying registration uncertainty
In order to analyse the spatial effects of dose mapping uncertainties we investigated how large volumes, V(ε D ) that got local dose mapping uncertainty (relative to the target prescription dose) larger than ε D . The V(3%) is shown in Figure 8 for the CTV, rectum and bladder for the beam models NBM and SBM. Without a planning margin (m = 0.0 cm), the V(3%) for the CTV vary between <1% and up to 10% as the 〈|T e |〉 vary between 0.2 cm and 0.4 cm. If a margin is applied the V(3%) decreases. The V(3%) vary in a logarithmic fashion for both rectum and bladder with the registration uncertainty. The use of a planning margin, or the larger spot sizes of NBM, result in larger V(3%) for the OARs, because of increased dose burden, but with a similar logarithmic behaviour with the registration uncertainty.    Figure 7 The uncertainty in TCP due to dose mapping errors. Isolevels of the decrease in the calculated TCP for the CTV, ΔTCP CTV , due to dose mapping errors caused by the image registration uncertainty 〈|T e |〉, shown for different planning margins m, and beam models NBM (left) and SBM (right).

Discussion
We have in this work developed a general framework for studying the influence image registration uncertainties have on the accuracy of the accumulated dose in fractionated radiotherapy. This framework was used to estimate the impact of the dose mapping uncertainties on the radiobiological quantities such as the TCP as well as the DVH quantities D 95 and V 80 for spot scanned proton therapy of the prostate. The results should be considered together to yield a complete view of the  The results indicate that the limiting accuracy requirement is that for the D 95 . A relaxed penumbra, or an increased target margin, decreases the uncertainty in the CTV while it, due to the increased dose burden, increases the uncertainties in the organs at risk. Jaffray et al. [33] highlighted the need for accurate dose accumulation for normal tissue and in that case the large V(3%) might be of concern, especially when accumulating dose to serial organs sensitive to local dose effects. The result of our particular simulations would to some degree depend on the choices of treatment site, patient geometry and deformation method. A shallower target would sharpen the penumbra of the proton treatment and increase the registration accuracy requirement. We assume a registration algorithm without systematic errors and stochastic, uncorrelated random deviations with a spatial frequency resulting from regularisation such that the registration error could be represented as a B-spline vector field with normal distributed control points. We consider these conditions for the registration uncertainty to be rather general for any regularised algorithm with random errors, although its applicability for certain clinically used DIR routines must be checked explicitly. The dosimetric parameters are scored locally ROI by ROI making the conclusions less sensitive to the fact that the DIR uncertainty in reality is not evenly distributed in space. Studies show [5,6,34] that registration algorithms can in some cases, e.g. lung and liver CT-CT registration, have a mean absolute error close to 0.1 cm and thus fulfilling the 3% above stated requirement for D 95 .
The uncertainty in the calculation of TCP obviously depends on the slope of the TCP sigmoid shape response curve, i.e. the γ 50 value. The parameters in this work were taken from the analysis of Cheung et al. [25] where γ 50 = 2.2 (95% CI 1.1, 3.2). A higher γ 50 would increase the sensitivity to an image registration uncertainty and thus increase ΔTCP. Roughly, using γ 50 = 3.2 instead of 2.2 would raise ΔTCP from 3% to approximately 4.4%, assuming the equivalent relative drop in dose when starting from TCP = 80%.
Future research using the presented framework could include other deformation models and other delivery techniques such as rotational therapy with photons.

Conclusions
We developed a simulation framework that described the total accumulated dose and the related response as a function of the deformable image registration uncertainty.
When applied to a spot scanned prostate treatment, the accuracy of D 95 on the accumulated dose is a limiting requirement on the deformable image registration when performing dose mapping. A required accuracy of 3% in D 95 would require a mean absolute image registration error uncertainty of 0.20 cm. By increasing the CTV to PTV planning margin with 0.1 cm, the mean absolute error can be relaxed to 0.3 cm with respect to the D 95 requirement. A few algorithms have reported accuracy better than 0.2 cm uncertainty for CT-CT registrations and thus meet the above stated [5,6,34]. However, for more challenging registration problems, such as multi-modality registrations, more development might be required to improve the registration accuracy, or to establish its accuracy as to define the dose mapping contribution to margin requirements in adaptive radiotherapy.