Fast and fully-automated multi-criterial treatment planning for adaptive HDR brachytherapy for locally advanced cervical cancer

Purpose: To develop and evaluate a fast, automated multi-criterial treatment planning approach for adaptive high-dose-rate (HDR) intracavitary + interstitial brachytherapy (BT) for locally advanced cervical cancer. Methods and materials: Twenty-two previously delivered single fraction MRI-based HDR treatment plans (SF clin ) were used to guide training of our in-house system for multi-criterial autoplanning, aiming for an autoplan quality superior to the training plans, while respecting the clinically desired ‘‘pear-shaped” dose distribution. Next, the conﬁgured algorithm was used to automatically generate treatment plans for 63 other fractions (SF auto ). The SF auto plans were compared to the corresponding SF clin plans in blind pairwise comparisons by an expert clinician. Then, the effect of adaptive autoplanning on total treatment (TT) plans (external beam + 3 BT fractions) was evaluated for 16 patients by simulating the clinically applied adaptive strategy to generate TT auto plans and compare them with the corresponding clinical treatments (TT clin ). Results: In the blind comparisons, all SF auto plans were considered clinically acceptable. In 62/63 comparisons, SF auto plans were considered at least as good as, or better than the corresponding SF clin . The average optimization time for autoplanning was 20.5 ± 19.2 s

In recent years, several approaches for high-dose-rate (HDR) brachytherapy (BT) automated treatment planning (autoplanning) have been proposed. For prostate cancer HDR-BT, Maree et al. and Cui et al. [13,14] approximated the Pareto-front for plan selection, allowing the user to make the final trade-off. A user-independent approach was presented by Breedveld et al. [15], who developed a system for automated multi-criterial treatment planning for prostate HDR-BT, based on their in-house, clinically applied autoplanning system for EBRT, and demonstrated superiority of automated treatment planning over conventional HDR-BT treatment planning.
A unique objective in locally advanced cervical cancer HDR-BT planning is to not only achieve adequate target coverage and acceptable organ-at-risk (OAR) doses, but to also realize this with a pear-shaped dose distribution. For cervical cancer HDR-BT autoplanning, Lessard et al. [16] proposed a dose-based objective function in which the weighted sum over the high-risk clinical target volume (CTV HR ), a manually contoured pear shape, and the OARs was optimized. The relative weights of the structures needed to be tuned per patient in order to achieve clinically favorable treatment plans. Hanania et al. [17] used tuned settings for the inverse planning algorithm from the TPS, but also required fine-tuning per patient. Guthier et al. [18] focused on target coverage maximization (weighted sum of CTV HR and intermediate risk CTV) under dose-volume based constraints on the OARs. They did not explicitly address other aims such as generating a pear-shaped dose distribution, and minimizing OAR dose. A way to avoid manual tuning for cervical cancer HDR-BT is presented by Shen et al. [19], who used deep reinforcement learning to mimic actions of a manual planner by predicting suitable weight adjustments. The optimization problem consisted of constraints on the CTV HR and an artificial pearshape, while the dose to the OARs was minimized by using a weighted-sum function. There was no explicit steering on tradeoffs between the dose in the target and OARs and there were no constraints on the dose to the OARs.
In this study, we follow the work by Breedveld et al. [15] to develop an approach for fully automated treatment planning for cervical cancer HDR-BT. Hereto, the autoplanning system was adapted for combined intracavitary + interstitial (IC + IS) HDR-BT for locally advanced cervical cancer. Two studies for autoplanning validation were conducted: (1) blind clinician comparisons of single fraction autoplans (SF auto ) with their corresponding clinical plans (SF clin ), and (2) comparisons of cumulative, total treatment (TT; EBRT + 3 BT fractions) dose distributions generated with adaptive autoplanning of the three BT fractions (TT auto ), with corresponding clinically delivered adaptive plans (TT clin ).

Patients
The patients in this study were treated between 2015 and 2018 for stage 2.B-4.A cervical cancer. Treatment consisted of EBRT (23 fractions of 2 Gy or 25-27 fractions of 1.8 Gy) and BT (three or four fractions spread over three weeks). Patients were treated using the Utrecht applicator (Elekta AB, Stockholm, Sweden) with up to 10 interstitial needles. In each BT fraction, after applicator implantation an MRI or CT was obtained for adaptive planning.
For this study, only MRI-based fractions were considered because they were fully contoured (CTV HR + OARs), resulting in a database with MRIs and corresponding clinical dose distributions of 85 fractions belonging to 34 patients. The average CTV HR volume was 29.4 ± 12.5 cm 3 (range 8.5-92.0 cm 3 ). In 56 fractions interstitial needles were implanted (average 4.4 ± 2.4 needles). A subset of 48 fractions belonged to 16 patients with an MRI in all three BT sessions available for this study. The remaining 37 plans belonged to 18 patients for whom one or more fractions were not selected because they were CT-based, or problems occurred when retrieving the data.

Clinical planning
State-of-the-art clinical treatment planning [20] was performed manually in our clinical treatment planning system Oncentra-Brachy (version 4.5.1, Elekta AB, Stockholm, Sweden), using delineated scans with a reconstructed applicator, while considering already delivered EBRT and BT dose.
An overview of the clinical planning aims as used between 2015 and 2018 at our institution is provided in Table 1. Limits and goals for the dosimetric parameters are defined for total treatments (TT: EBRT + BT) and presented as equivalent doses of 2 Gy per fraction (EQD 2 ). The dose received during EBRT was assumed to be uniform. No deformable registration was applied before summing EBRT and BT dosimetric plan parameters [21].
The general goal at the start of BT treatment planning was to evenly distribute remaining OAR tolerance doses and the required additional CTV HR dose amongst three or four fractions. However, in the first fractions higher OAR doses could be accepted if it was expected to be necessary to obtain the minimum CTV HR dose of 85 Gy in the total treatment. The maximum OAR tolerances in the total treatments were however always respected.
Clinical BT planning started from a standard normalized treatment plan, in which the dwell positions in the intracavitary part of the applicator were evenly active. Dwell positions in the needles were then activated and dwell times in the standard plan were adjusted graphically by dragging isodose lines in the TPS. This process took 15-30 min per plan.

Automated planning
Erasmus-iCycle was used as a basis for automated plan generation [10,11,15]. The system uses a so-called 'wish-list' to steer the multi-criterial optimization. The wish-list contains hard planning constraints and prioritized objectives, and defines the lexicographic plan generation. A well-tuned wish-list results in clinically acceptable treatment plans with favorable trade-offs between the treatment objectives. The treatment site specific wish-lists are constructed in an iterative tuning process in close collaboration with the expert clinicians, using repetitive autoplanning, plan evaluation and wish-list adjustments for a small set of training patients (Electronic appendix of [11]).

Wish-list
The wish-list used for autoplanning in this study ( Table 2) was configured with an expert clinician (Dr. Jan-Willem Mens, JWM), in line with clinical planning. Twenty-two of the 85 available treatment fraction MRIs with corresponding clinical dose distributions were used for training. Training fractions were selected to cover the full range of CTV HR volumes. Because of the adaptive planning, the wish-list contains fraction-specific dosimetric parameters. All activated dwell positions (IC + IS) were available during automatic optimizations.
In line with the clinical protocol (Table 1), the final wish-list contained hard constraints for the D 2cc of the OARs (bladder, rectum, sigmoid and small bowel). The DVH-based cost-functions in the wish-list were implemented similar to the approach presented in [15]. To constrain the dwell time modulation and enforce smooth dose distributions, a quadratic cost-function which penalizes the second derivative of the dwell times of adjacent dwell positions in the intracavitary applicator, was used as constraint, similar to fluence map smoothing in EBRT [22]. The dose in the CTV HR was optimized in three steps (priorities 1, 5 and 7 in Table 2). A pear-shaped dose distribution was created by optimizing on two artificially created structures ('pear' and 'pear-inside') which follow the dwell positions in the tandem and ovoids of the applicator Table 1 Overview of the clinical planning aims. Presented EQD 2 are for total treatments (TT), assuming a/b = 10 Gy for the CTV HR and a/b = 3 Gy for OARs.  Table 2 Wish-list for per-fraction autoplanning for HDR-BT for cervical cancer, and total cumulative treatment aims for EBRT + 3BT (last column). Per-fraction parameters (L i j dose limit for structure i in fraction j, and G i;p j dose goal for structure i in priority p and fraction j) in the wish-list were patient-specific and related to the performed adaptive planning (Section ''Clinical planning"). N equals the number of needles. All doses are presented in EQD 2 . The arrows indicate whether the objective was minimized (down) or maximized (up).

Constraints per fraction
Structure i Constraint function Dose limit in fraction j,L i j Cumulative dose limit    at distances of 9 mm and 5 mm respectively, excluding overlap areas with the target or OARs (priorities 2 and 3). The 4th priority in the wish-list aims at reducing the relative dose contribution from the needles using s needles as defined by: where t total represents the total accumulated dwell time for all dwell positions (IC + IS) and t needles the total dwell time for only the interstitial part of the implant. The goal value was related to the number of available needles N ( Table 2). The priority 6 objective aims at reducing OAR dose, where dose reduction in the small bowel was considered more relevant than other OARs and was thus assigned a relative weighting factor of 4.

Blind comparisons of single fraction autoplans (SF auto ) and clinical plans (SF clin )
To enable fair comparisons between SF auto and SF clin plans, the fraction-specific parameters L i j and G i;p j in Table 2 were based on the constraints and goals used for the corresponding SF clin . In this way, plans based on the same treatment goals could be compared, the only difference being the way of plan generation, manually or automated.

Comparison of total treatment adaptive autoplans (TT auto ) with clinical plans (TT clin )
For the TT comparisons, an automated adaptive strategy was simulated, in line with clinical practice (Section ''Clinical planning"). To compute L i j and G i;p j , the total dose received up to this fraction (EBRT and 0, 1 or 2 automatically generated BT fractions) was taken into account, while aiming at an equal spread of required or allowed dose to the CTV HR and OARs over the remaining fractions. However, if the absolute minimum cumulative dose goal of 85 Gy (EQD 2 ) to the CTV HR seemed to be infeasible, up to 5% more dose than L i j was allowed to the OARs in the first fraction and up to 3% in the second fraction. In the third and final fraction, the maximum tolerances (Table 1) were always respected.

Optimization details
The dose optimization points (voxels) were sampled with a density of 300 voxels/cc for all structures. To speed up computations, only the parts of the OARs within a 35 mm radius from the dwell positions were taken into account during optimizations. Beyond this distance, the maximum expected dose is anyway much less than the constrained values, which was verified by visual inspection for the training patients and had no impact on the resulting plan. Both the dose-volume cost-functions and the relative needle contribution cost-function (Eq. (1)) are nonconvex. We used our in-house developed interior-point solver with extended functionality for non-convex optimization for these functions, see [23,15]. The optimizations were performed on an Intel Core i7-7700 with 4 cores running at 3.6 GHz. After optimization in Erasmus-iCycle, the dwell-times were imported in the clinical TPS for final dose calculation and comparisons with clinical plans.

Validation of automated planning
First, blind comparisons of single fraction autoplans (SF auto ) with corresponding clinical plans (SF clin ) were performed for the 63 SF clin plans that were not used in wish-list tuning. The pairwise plan comparisons were performed in the clinical TPS by an expert clinician (JWM), using the 3D dose distributions and DVH parameters of these plans, while considering also DVH parameters of previous BT fractions and the EBRT course, that were also provided. There were no pre-defined criteria for the comparisons; the clinician followed his routine workflow for plan evaluation. The clinician first assessed the clinical acceptability of both plans, followed by assessments of differences in (i) overall quality, (ii) CTV HR dose and (iii) OAR dose, using a visual analog scale (Appendix A).
Second, total treatment adaptive autoplans (TT auto ) were compared with the clinical plans (TT clin ) for the 16 patients with an MRI in each BT fraction. Pairwise differences in the dosimetric plan parameters listed in Table 1 were assessed in the clinical TPS. After importing the dwell-times of the autoplans in the clinical TPS, the DVH parameters did not exactly match the Erasmus-iCycle optimized ones due to small differences in the implementation of volume and dose-point definition between the two systems. Therefore, D 90% CTV HR in the last TT auto fraction could be rescaled up to 93.8 Gy if possible within OAR limits, or down to 95 Gy in case the D 90% of TT auto exceeded 95 Gy. Statistical significance of differences in plan parameters was assessed with the paired Wilcoxon signed-rank test. Fig. 1 shows an example of an SF auto plan, compared to the corresponding SF clin plan, clearly showing the desired pear-shape in both dose distributions. All 63 SF auto plans were considered clinically acceptable by the clinician, while 1 clinical plan was (in retrospect) not. The clinician's scores, presented in Fig. 2, demonstrate superiority for the automatically generated plans in terms of overall plan quality. In 60/63 cases SF auto was preferred over SF clin , in 2/63 cases the quality of the plans was considered equal (score = 0) and for 1/63 cases SF clin was preferred over SF auto because of a more favorable small bowel dose. The scoring for CTV HR shows a similar trend with in 62/63 cases SF auto similar or better than the clinical plan. For the OARs, the differences were less pronounced. Still, in 29/63 cases SF auto was considered better, in 21/63 cases both plans were similar and for 13/63 cases SF clin was preferred.

Results
The average optimization time for the 63 SF auto plans was 20. 5 ± 19.2 s (range 4.4-106.4 s).
For the total treatment adaptive autoplans, the last fraction could be rescaled for 10/16 patients due to differences in TPSes to improve CTV HR dose, with an average absolute difference in D 90% of 0.26 ± 0.21 Gy (range 0.10-0.80 Gy). Total treatment dosimetric parameters are compared in Fig. 3. The required minimum CTV HR D 90% of 85 Gy (Table 1) was always obtained for TT auto , whereas for two TT clin plans it was not. Consequently, these two clinically delivered treatments did not strictly adhere to the requirements for clinical acceptability. This was due to an unfavorable anatomy of these patients, making manual treatment planning challenging. The first goal for the CTV HR was to achieve a minimum dose of 90 Gy. This was obtained in 14/16 TT auto plans, compared to 9/16 TT clin plans. For all 16 patients the CTV HR D 90% was highest in the TT auto plan. This improvement was statistically significant (p = 0.0004) and clinically relevant with a mean D 90% for the TT auto plans of 93.0 ± 2.0 Gy compared to 89.4 ± 3.2 Gy for TT clin , with differences ranging from +1.4 to +6.0 Gy.
For the D 2cc of the OARs for the total treatment adaptive plans, the prescribed limits (Table 1) were never exceeded, neither for TT auto , nor for TT clin (Fig. 3). The dose in the bladder was significantly reduced in the TT auto plans compared to TT clin (p = 0.05) with a mean reduction of 0.87 Gy (range À1.8 to 6.6 Gy). The rectum was also significantly more spared (p = 0.04), with a mean reduction in D 2cc of 1.4 Gy (range À4.9 to 6.3 Gy). There were no significant differences in D 2cc of the sigmoid (p = 0.3) and small bowel (p = 0.3).

Discussion
Optimization of the dose distribution for locally advanced cervical cancer HDR-BT is a complex problem, due to the special requirements of the dose distribution (pear-shaped), delivery restrictions (contributions by needles) and the adaptive treatment schedule. To the best of our knowledge, we have presented the first automated multi-criterial treatment planning solution considering each of the complexities, and capable of generating clinically favorable dose distributions. In the blind comparisons, all 63 single fraction autoplans were considered clinically acceptable by the expert clinician, and 62 of the 63 autoplans were scored similar or better than the clinically delivered plans.
In this study, we also proposed an adaptive treatment approach with a quantitative recipe for setting all fraction-specific treatment goals, depending on previously delivered doses. This adaptive Fig. 3. Scatter plots showing CTV HR D 90% , and bladder, rectum, sigmoid and small bowel D 2cc for the automatically generated and clinical total treatment plans (TT auto and TT clin ). Each patient has her own colored symbol. Improvements in CTV HR D 90% for TT auto plans, and bladder and rectum D 2cc were statistically significant. There were no significant differences for the sigmoid and small bowel.
autoplanning approach was in line with the (somewhat less rigidly defined) clinical adaptive planning. Adaptive autoplanning was compared with clinical planning based on cumulative total treatment doses, as there exists no clinical scenario for comparisons of single fraction doses in adaptive treatment. Nonetheless, all single fraction autoplans fulfilled the clinical requirements (Table 1). The proposed procedure with adaptive single fraction autoplanning resulted in overall higher quality of total treatment dose distributions compared to clinical total dose. For all 16 patients, the CTV HR D 90% was highest with automated planning. The number of patients who reached the desired minimum CTV HR dose increased from 9/16 in the clinical treatment to 14/16 when using automated treatment planning.
With an average optimization time of 20.5 s per fraction, this approach shortens treatment planning times considerably, which are currently between 15 and 30 min in clinical practice.
The wish-list in this study was tailored to the treatment planning procedure at the Erasmus MC at the time of the treatment of the patients. The Erasmus MC has recently made the transition to the Embrace II protocol 1 . As most challenging requirements in this protocol are similar to those for the current HDR-BT cervix plans, it is expected that the current wish-list can be adapted to generate clinically preferable dose distributions for Embrace II as well. In addition to limits and goals for the D 90% of the HR-CTV and D 2cc of the OARs, the EMBRACE II protocol includes more goals, for example for the D 98% of the CTV HR , intermediate-risk CTV and gross tumor volume. Because the Embrace II protocol contains more objectives compared to the protocol in this study, manual treatment planning becomes even more challenging and an automatic approach could become even more valuable. However, tuning the desired tradeoffs between the different objectives will be more time-consuming.
In the future, our HDR-BT autoplanning approach could be extended with optimization of the interstitial needle configuration (positions and insertion depth) of the applicator, by simulating the 10 possible needle positions and including an objective to mini-mize the number of used needles and their insertion depths in the wish-list. Needles could then be implanted based on an individualized and automatically generated treatment plan that is optimized both for the dwell times and number and position of needles.
In conclusion, this study demonstrated that fast, fullyautomated multi-criterial treatment planning for locally advanced cervical cancer HDR-BT was feasible, based on a well-tuned, clinically relevant wish-list. Blind pairwise clinician comparisons of single fraction manual-and autoplans pointed at a strong preference for the autoplans. Cumulative total doses resulting from adaptive treatment were also favorable when generated with perfraction autoplanning. With autoplanning, planning times reduced from the current 15-30 min to 20.5 s on average.

Disclosure
The Erasmus MC Cancer Institute has research collaborations with Elekta AB, Stockholm, Sweden, and Accuray Inc, Sunnyvale, USA. The Delft University of Technology-Radiation Science & Technology has research collaborations with Varian Medical Systems, Palo Alto, USA. The presented work was not sponsored by any of the above named companies.