Probabilistic optimization of dose coverage in radiotherapy.

BACKGROUND AND
PURPOSE: Probabilistic optimization is an alternative to margins for handling geometrical uncertainties in treatment planning of radiotherapy where uncertainties are explicitly incorporated in the optimization. We present a novel probabilistic method based on the same statistical measures as those behind conventional margin based planning.
MATERIAL AND METHODS: Percentile Dosage (PD) was defined as the dose coverage that a treatment plan meet or exceed to a given probability. For optimization, we used the convex measure Expected Percentile Dosage (EPD) defined as the average dose coverage below a given PD. An iterative method gradually adjusted the constraint tolerance associated with the EPD until the desired target PD was met. It was applied to planning of cervical cancer patients focusing on systematic uncertainty caused by organ deformation. The resulting plans were compared to margin based plans using target and organ at risk PDs.
RESULTS: The EPD tolerance converged in less than ten iterations to produce a PD within 0.1 Gy of the requested. The PD was on average within 0.5% of the requested PD when validated versus independent scenarios. The rectum volume, extracted from the PDs, receiving 90% of the intended target dose was decreased with 16% for the same target PD in comparison to margin based plans.
CONCLUSIONS: The proposed probabilistic optimization method enabled prescription of a dose volume histogram metric to a chosen confidence. The probabilistic plans showed improved target dose homogeneity and decreased rectum dose for the same target dose coverage compared to margin based plans.

Introduction

Geometrical uncertainties in radiotherapy, such as patient setup errors and organ motion, are commonly handled by expanding the clinical target volume (CTV) by a margin to create a planning target volume (PTV) to which the dose intended for the CTV is prescribed [1], [2]. A widely used recipe [3] for setting the margin was designed such that the entire CTV should be covered to receive 95% of the intended dose for 90% of the treatments. However, such margin recipes are intrinsically population based and do not explore individual patient geometries to reduce collateral dose burdens to surrounding normal tissues. Probabilistic or robust optimization strategies that explicitly include the uncertainties have been proposed to improve upon the margin based planning [4]. However, probabilistic strategies in general are not based on the same statistical basis as the margin recipe. Hence, introduction of these alter the principles of dose prescription as compared to margin based plans. Gordon et al. [5] formulated a probabilistic optimization strategy based on the statistics of the margin recipe which was extended by Mescher et al. [6] to include dose coverage as constraints. However, none of these formulations are convex thus rendering implementation of efficient optimization algorithms difficult.

The dose coverage is defined as the dose that the partial volume of a region of interest (ROI) at least receives, normally extracted by means of a cumulative dose volume histogram (DVH). The th percentile dosage (PD) (strictly really the th percentile) is defined such that the probability for at least a dose coverage is . We can use as a target minimum dose measure which is consistent with the common margin recipe. A mathematically equivalent measure to PD is the Value at Risk (VaR) that has been used in stock portfolio management to estimate the maximum loss of the th percentile worst outcome of the portfolio [7]. However, VaR and PD yields non-convex optimization problems [8]. Luckily, the Conditional Value at Risk (CVaR), defined as the expectation value of the distribution within the th percentile, gives a convex problem suitable for optimization. The concept was derived for investment management [9] but has also been applied to radiotherapy treatment planning [10], [11]. To suit the later, we introduce the term Expected Percentile Dosage (EPD), defined as the expectation value of the dose coverage for all scenarios with worse dose coverage than the PD (i.e. mathematically equivalent to CVaR).

In this work, we take advantage of the close statistical relationships between PD, EPD and the margin concept, creating a novel method for probabilistic and hence robust treatment planning. The algorithm satisfies a target PD constraint while minimizing the expectation value of the organ at risk (OAR) objective functions. To test the feasibility of this methodology, we applied it to radiotherapy planning of cervical cancer patients for which systematic geometrical uncertainties were modelled by a deformation model [12]. The resulting probabilistic treatment plans were evaluated versus margin based plans satisfying the same target PD. For simplicity, the effect of fractionation was not considered, although our method is general enough for future inclusion. To our knowledge, probabilistic planning taking organ deformation into account has previously never been explored for cervical cancer patients.

Methods

Dose coverage probability and the expected percentile dosage

In our simulations, we defined a scenario as one possible realization of a treatment with data for all relevant uncertainties sampled accordingly. We focused on patient setup errors as well as position and shape of the target and surrounding risk organs as they exhibit large variations for cervical cancer patients. Uncertainties in segmentation could have been included but were left out for simplicity. Also, single fraction scenarios were used implying that only systematic errors were modelled but the concepts of EPD is general enough to include also fraction specific random errors for fractionated treatments.

The shape and position of the patient anatomy was described by a statistical shape model (SSM) where the parameters are sampled from the probability density function of a random vector [12]. The SSM was derived through principal component analysis (PCA) of deformable image registration results of images acquired at different fractions during the radiotherapy. The sampled displacement transform from the SSM is applied to the spatial position of the reference image to generate a randomized geometrical state of the patient anatomy.

The dose distribution without positional error was calculated using the fixed reference geometry. Then the dose distribution for the anatomy in the reference geometry for a sampled single fraction scenario could be given bywhere represent the beamlet intensities (and hence the optimization variables) defining the irradiation. Note that Eq. (1) assumes that the dose distribution in the treatment machine coordinate system do not change due to patient or organ motion. This is a reasonable approximation for male pelvic patients [13] and a reasonable approximation also for females. Since is a random vector, it follows that for a fixed also the dose coverage is a random variable. Fig. 1 illustrates how the probability density function for can be determined from repeated treatment scenario simulations where the input elements of are sampled from their respective probability density functions for each scenario (see e.g. Tilly et al. [12] for a detailed description). Following the recommendations of the ICRU we used to represent the minimum target dose [14]. From we could calculate the probability for to be at least some number as

Fig. 1

In this schematic example the probability density function (right) is estimated from extracting the values determined from the scenario DVHs along the dashed arrow (left figure). The percentile dosage and expected percentile dosage can then be determined from the constructed . The is the that 90% of all scenarios meet or exceed (i.e. dashed line at the border of the red/white area). The is the average of all scenarios with worse dose coverage than (i.e. the average of the red area). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

We defined as the dose that exceeds with the probability

We could then express a probabilistic prescription using a PD criteriawhere is the prescribed PD with its probability explicitly stated. Margin based plans often use 95% of the intended homogenous target dose as the [3]. Similar criteria can be defined for OAR by requiring that a relevant does not exceed a tolerated dose.

PD is not convex in and therefore we instead used the EPD that is convex and defined as

For the implementation of the optimization we used surrogate functions closely related to rather than a direct implementing of Eq. (5).

The for any dose-volume combination and fixed , can be determined from simulation results by first creating dose volume coverage maps [5] describing the probability that a volume will receive at least a certain dose (i.e. the probability that a DVH will pass above/right of a dose-volume combination). We chose to use for targets and for OARs where was obtained from the dose coverage maps by extracting their 90% iso-probability line and in analogous way. The probabilistic and margin based plans were compared using the relevant for targets and for OARs. Analogous to DVHs for OARs the can be expressed in terms of the volume receiving at least a dose to the probability .

Probabilistic optimization using an EPD constraint

The quantities and are non-convex with respect to [15] and cannot easily be used to solve for with standard optimization methods. The planning paradigm chosen in this work was to treat the PD criterion of Eq. (4) as a constraint while minimizing the sum of objective functions describing all other aspects of the treatment plan (e.g. limiting the dose to the OARs). We therefore employed convex objective functions, routinely used in radiotherapy optimization, of the typefor a minimum dose coverage criterion, and similar for a maximum dose coverage criterion. is the volume of :th region . The function penalizes voxels with dose below whereas the analogous max dose function penalize voxels above . With a suitable choice of the behaviour of will be a good surrogate for a function directly measuring the deviation of from . We empirically found that setting 5% higher than yielded a very good proxy for the criterion.

For a given , is a random variable with a probability density function, . Percentile quantities of (i.e. and ) were defined in analogy with PD and EPD and used as surrogates for the for the latter. Rockafellar and Uryasev [9] introduced a help function that in our context (and discretized over scenarios) could be expressed as:where is the number of scenarios, are scenario specific samples from the probability density functions of and is the CTV region. Moreover, they showed that can be minimized by minimizing with respect to and , which we utilized in this work. We limited the dose to the OARs and prevented overdosage of the CTV by minimizing the weighted sums of functions of type with suitable values per ROI. The probabilistic optimization problem can then be formulated aswhere is the importance weight for , is number of ROIs, expresses the EPD constraint, and denotes the positivity requirement on the beamlet weights. The value of , and indirectly the value of , is then limited by the constraint tolerance that was iteratively adjusted until was fulfilled. Note that since and are convex, the whole problem is convex and hence standard methods for non-linear constrained optimization could be used. By controlling and via we also controlled the PD since is a surrogate for . In summary, a solution to problem (8) yields a treatment plan with a PD that is determined by while trying to limit the dose to OARs and also limit the overdosage to the CTV. We implemented the method introduced by Baum et al. [16] to efficiently calculate the weighted sums in the objective function by first determine the probability that occupies voxel and thus inter-changing the summation order over voxels and scenarios.

In the following we will outline a bilevel optimization procedure where the inner loop solves problem (8) and the outer loop determines the tolerance that yields the requested PD.

Optimization towards a percentile dosage

The planning paradigm used in this project was to optimize towards a given target dose coverage constraint . We hypothesized that given a good surrogate function for deviation of from the requested there is a value for which problem (8) yields . Finding by trial and error is possible but may be inefficient. Instead we repeatedly solved the inner problem (8) and used the difference between and to gradually tighten in an outer loop until the desired criterion was met ( sufficiently close to ). The inner loop problem (8) was implemented in C++ using the Ipopt [17] as a solver, see Supplement A for more details. Details of the outer loop are given in Supplement B.

Several hundred simulated scenarios are in general needed to provide sub-percent accuracy in the calculated . Due to high memory load and long calculation times we used  = 100 in the inner optimization problem and the influence of this was quantified in the results.

Probabilistic optimization applied to cervical cancer

Patient data

We used the same data set as in a previous work where we developed a probabilistic tool for treatment plan evaluation using a SSM to sample the large deformations of cervical cancer patients anatomy for five patients (with informed consent) treated at the Academic Medical Center, Amsterdam [12]. These patients were treated with single-arc rotational therapy 6MV and was prescribed . The primary CTV encompassed the GTV, cervix, corpus-uterus, and the upper 20 mm of the vagina with the left and right lymph nodes included into the total CTV.

Scenario sampling

The uncertainties considered were systematic setup-errors and organ motion for single fraction treatments. The deformations were sampled using the previously developed SSM. Segmentation uncertainty was not explicitly included but could to some degree be considered as part of the model since its learning data included manually drawn contours. The deformations were characterized by eigenmodes (seven was sufficient) resulting from a PCA with the first eigenmode representing the most dominating deformation pattern. As in our previous work [12] we made the assumption that a representative deformation scenario can be sampled as a weighted sum of eigenvectors. The setup-error was assumed to be normally distributed with zero mean and standard deviations of (0.3, 1.3, 1.9) mm in the left-right, superior-inferior and anterior-posterior directions [18].

Implementation of EPD constrained optimization

We used the CTV constraint to ensure the minimum dose requirement. All the OARs (bladder, rectum and the rest of the normal tissue) as well as the CTV were assigned a max dose objective per scenario (see details in Table 1). A conformance objective enforced a dose fall-off gradient from the edge of each scenario CTV as well as from the edge of the union of all scenario CTVs.

Table 1

Treatment planning criteria (implemented in Eq. (8)). The dose is given as percent of intended dose 46 Gy.

ROI	Prob.	Margin	D_98%,90% [%]	Weight	Min dose [%]	Max dose [%]
PTV		✓		Constraint	95
PTV		✓		1000		101
CTV	✓		95	Constraint
CTV	✓			1000		101
Bladder	✓	✓		1		90
Rectum	✓	✓		10		90
Normal Tissue		✓		100		98%→50% in 2 cm from PTV edge
Normal Tissue	✓			100		98%→50% in 2 cm from edge of union of scenario CTVs
Normal Tissue	✓			10		98%→50% in 2 cm from edge of scenario CTVs

Note that the dose was calculated only in the fixed geometry, see Eq. (1). The beamlets for use in dose calculation were generated with a Monte Carlo dose engine [19], implemented in a research version of Monaco treatment planning system (Elekta AB, Stockholm, Sweden), to a 3% statistical uncertainty per beamlet. The 360° arcs for VMAT planning were approximated by 36 beam directions with 10° intervals. The resolution of the fluence map bixels were 10 × 10 mm2 and the dose grid voxel resolution was 3 × 3 × 3 mm3.

For verification of the generated probabilistic treatment plan, a set of 1000 scenarios (different from those used in the optimization) was generated for each patient. The PDs based on the verification calculations were evaluated versus the targeted probabilistic dose coverage criterion .

Comparison with conventional margin based plans

As far as possible the treatment planning parameters were kept the same for the conventional margin based plans as the probabilistic plans to facilitate a fair comparison. The clinically used PTV max dose criteria at our department enforces that  < 2% which required a high PTV max dose weight, see Table 1. Similar to the CTV in the probabilistic plans, a conformance objective enforced a dose gradient from the edge of the PTV. The conventional plans were designed to meet the same criteria as the probabilistic plans by stepwise increasing the margin size until the prescribed was fulfilled. The was determined by treatment scenario simulation as described in Fig. 1. Here, we found that a 20 mm CTV to PTV margin resulted in a equal to the . The same beamlets, objective functions and , objective weights and the same solver was used as in the probabilistic optimization.

Results

The iterative method of gradually finding the constraint tolerance that produces the requested percentile dosage converged as intended for all patients such that the were within 0.1 Gy (0.2%) of after 10 iterations, see Fig. 2. The from the verification calculations were on average 0.2 Gy (0.5%) lower than . The EPD, or , were near identical 42.6 Gy for the two planning paradigms.

Fig. 2

The evolution per iteration of the . The rightmost points (without connecting line) of the left panel shows the determined from the validation simulations using 1000 scenarios.

A representative dose distribution for a sagittal slice through the center of the CTV is shown in Fig. 3. The margin based plan exhibited the designed homogenous dose in the PTV and a good conformance. The probabilistic plan redistributed dose from areas with little deformation (superior, inferior and posterior) to the area with large deformation (anterior). Moreover, for the probabilistic plans the dose inside the CTV of the reference geometry was higher whereas it was lower in the rectum in the reference geometry.

Fig. 3

Dose distribution comparison between the margin based plan (left) and the probabilistic plan (right) for a sagittal plane through the middle of patient #4. The PTV is shown as a reference in the probabilistic dose although not used in the optimization. The dose is displayed relative the intended target dose 46 Gy. The CTV to PTV margin was 20 mm for the margin based plans.

The probabilistic evaluation of the dose distribution from Fig. 3 is shown as the CTV and the OAR in Fig. 4. The probabilistic planning consistently resulted in increased dose homogeneity in the CTV with near equivalent , and increased median dose closer to the intended dose, and lower max dose , see Table 2. Furthermore, the rectum volume receiving high dose was decreased. The decrease in the high dose volume for the bladder was however not consistent across all patients.

Fig. 4

A comparison between the probabilistic and margin based plan for patient #4 using the for the CTV and for the bladder and rectum.

Table 2

The treatment evaluation statistics. All doses are given in Gy. The prescribed was 43.7 Gy. The () is the value for the probabilistic plans minus the corresponding value for the margin based pans. The is the difference in relative volume for the given relative dose.

Patient #	CTV D_98%,90%		CTVΔD_50%,50%	CTVΔD_2%,90%	Bladder ΔV_90%,10%	Rectum ΔV_90%,10%
	Prob.	Margin
1	43.7	43.8	0.61	−0.74	−18.2	−0.72
2	43.6	43.9	0.65	−0.77	−13.2	−5.9
3	43.4	43.7	0.79	−1.59	14.4	−19.9
4	43.4	43.3	0.91	−1.04	12.9	−33.5
5	43.4	43.8	0.78	−1.40	−10.5	−18.4
Average	43.5	43.7	0.75	−1.11	−2.9	−15.7

Another way to compare the plans is to calculate the treated volume, defined as the total volume receiving at least 95% of the intended dose 46 Gy. The 95% isodose mark the ideal anisotropic margin to produce the prescribed . The probabilistic plans showed an average treated volume reduction of 18.8% (ranging from 7.9% to 28.8%) compared to the margin based plans.

Each treatment optimization, i.e. solving the inner problem (8), with took ∼2 h on an Intel Core i7-6820HQ CPU @ 2.7 GHz computer with 32 GB RAM (single threaded).

Discussion

In this work we developed a novel probabilistic optimization framework based on the percentile dosage (PD) which is the same statistical basis as the van Herk margin recipe. This also demonstrates for the first time probabilistic optimization of cervical cancer radiotherapy where organ deformation is taken into account. The probabilistic plans showed an improved target conformance while at the same time reducing the rectum volume receiving high dose when compared to margin based plans. The method enables a dose planner to directly plan for a target percentile dosage.

The target PD criterion was met by gradually tightening with the proposed iterative algorithm. The convex reformulation of the optimization problem, in contrast to Gordon et al. [5] and Mescher et al. [6], enabled the use of standard optimization methods to find the global optima.

The sampled scenarios are from a population based shape model. Even so, the decrease in treated volume (18.8%) can be set in relation to individual margin creation. It has been shown that the margin (systematic and random) could be decreased to 7–10 mm with repeat imaging of the patients [20]. For the patients in this work, a ten mm PTV margin would be an average volume reduction of 25.1% compared to a 20 mm PTV margin. An alternative PTV might be created by adjusting the margin to the 95% isodose of the probabilistic plans. This would likely improve the margin based plans but this evaluation was out of scope for this work.

The transformation from an optimization problem using DP to instead use EDP enabled a convex optimization where the lower tail of is controlled to yield a requested . There could potentially be a solution to Eq. (8) that yields a highly irregular for which the iterative method of tightening the constraint tolerance might fail. Note that margin based plans do not consider the shape of the at all and the resulting differences in is worth further study. The near identical EPD in this work does not indicate that one or the other method produce more conservative plans.

The choice of threshold dose in the penalty function is non-trivial. The solution will be close to the requested PD given a close relationship between and . The naive choice of in would be . However, this would imply a very small to meet the dose coverage constraint. We found that the optimizer had problems finding a solution to problem (8) for very low . Increasing to 5% above and thereby increasing the for the same improved the performance of the optimizer. Note that the increased does not imply an increased target dose coverage since will be relaxed to ensure the prescribed . Even though the same is achieved using different surrogate functions the associated distribution may be different.

The beamlet weights as optimization variables must for photons be converted into machine parameters to enable delivery of the generated plan. This can be done using the resulting dose distribution as an input to a subsequent optimization to recreate the same dose distribution but with machine parameters as decision variables, see e.g. [21]. The current implementation was not optimized for performance, e.g. only single threaded, such that 10 iterations with progressively tighter constraint tolerance may take ∼20 h.

In conclusion, we have demonstrated a novel probabilistic optimization framework that enables radiotherapy treatment planning towards a requested dose coverage probability. The expected percentile dosage provided a convex measure that facilitates finding the global optima during optimization. An iterative scheme that gradually tightened the constraint tolerance of the expected percentile dosage provided convergence towards the percentile dosage planning goal. Verification calculations for the generated plans using 1000 new scenarios independent of the optimizations showed that the requested dose coverage was met within 1.2%. The resulting probabilistic plans outperformed conventional margin when evaluated using percentile dosages by showing improved target dose homogeneity and decreased rectum dose.

Conflict of interest statement

We declare no conflict of interest.