SISS-MCO: large scale sparsity-induced spot selection for fast and fully-automated robust multi-criteria optimisation of proton plans

Objective. Intensity modulated proton therapy (IMPT) is an emerging treatment modality for cancer. However, treatment planning for IMPT is labour-intensive and time-consuming. We have developed a novel approach for multi-criteria optimisation (MCO) of robust IMPT plans (SISS-MCO) that is fully automated and fast, and we compare it for head and neck, cervix, and prostate tumours to a previously published method for automated robust MCO (IPBR-MCO, van de Water 2013). Approach. In both auto-planning approaches, the applied automated MCO of spot weights was performed with wish-list driven prioritised optimisation (Breedveld 2012). In SISS-MCO, spot weight MCO was applied once for every patient after sparsity-induced spot selection (SISS) for pre-selection of the most relevant spots from a large input set of candidate spots. IPBR-MCO had several iterations of spot re-sampling, each followed by MCO of the weights of the current spots. Main results. Compared to the published IPBR-MCO, the novel SISS-MCO resulted in similar or slightly superior plan quality. Optimisation times were reduced by a factor of 6 i.e. from 287 to 47 min. Numbers of spots and energy layers in the final plans were similar. Significance. The novel SISS-MCO automatically generated high-quality robust IMPT plans. Compared to a published algorithm for automated robust IMPT planning, optimisation times were reduced on average by a factor of 6. Moreover, SISS-MCO is a large scale approach; this enables optimisation of more complex wish-lists, and novel research opportunities in proton therapy.


Introduction
Multiple solutions have been proposed for automated treatment plan generation for external beam radiation therapy (EBRT) with photon beams (Cotrutz and Xing 2003, Craft et al 2005, Monz et al 2008, Breedveld et al 2012, Rossi et al 2012, Yuan et al 2012, Xhaferllari et al 2013, Breedveld et al 2019, Bijman et al 2021, McIntosh et al 2021, Schipaanboord et al 2022).Often automation has led to better or comparable plan quality with vast reductions in planning time, planning workload, and improved consistency (Voet et al 2013, Fogliata et al 2014, Voet et al 2014, Sharfo et al 2016, Hansen et al 2017, Heijmen et al 2018, Hussein et al 2018, Marrazzo et al 2019, Bijman et al 2020, Fiandra et al 2020, Fjellanger et al 2021, Rossi et al 2021).Erasmus-iCycle has been developed in our centre for wish-list driven automated multi-criterial optimisation (MCO) of treatment plans (Breedveld et al 2012).
EBRT with proton beams is an emerging treatment modality with potential benefits in reducing OAR dose through the Bragg peak.However, plan quality may suffer due to anatomical variations or setup errors (Engelsman et al 2013).Several automated intensity modulated proton therapy (IMPT) planning methods have been proposed, including data-driven approaches like RapidPlan PT (van de Water et al 2013, Tol et (Stockholm, Sweden).RP PT uses a library of historical treatment plans to predict achievable dose volume histograms (DVHs) for OARs, enabling automated treatment planning in under 10 min (Delaney et al 2018b).RayStation's methods include photon-dose distribution mimicking (Kierkels et al 2019), an atlas regression forest-based model (van Bruggen et al 2022), and deep learning models with U-net architectures (Eriksson and Zhang 2022).
RP PT was validated for prostate cancer patients (PC), with Delaney et al using non-robust optimisation for 5 patients (Delaney et al 2018a) and Xu et al using a 13-scenario robust setting for 15 patients (Xu et al 2021).Delaney et al reported comparable plan quality to manual plans but noted that OAR dose degradation in RP PT did not always correlate with geometric-outlier warnings, making manual plan checks important (Delaney et al 2018a).For head and neck cancer (HNC), Delaney et al (2018b) validated RP PT in a multi-centre setting without robustness.A model was constructed from 50 manual HNC IMPT plans from a single centre.Three external centres each provided seven manual benchmark IMPT plans.Delaney et al achieved comparable or better results than benchmark plans.However, plan quality for test patients was compromised if some OARs or targets were not included in the model or were delineated differently (Delaney et al 2018b).Celik et al (2021) validated RP PT for gastro-oesophageal cancer with a training set of 45 and a test set of 15 patients in a robust setting.Although RP PT demonstrated added value, some discrepancies were observed, and the authors suggested interactive optimisation of dose-volume constraints for certain patients to reduce discrepancies (Celik et al 2021).Kierkels et al (2019) proposed a data-driven approach, using VMAT or IMRT plans for dose-mimicking and OAR dose reduction in IMPT scenarios.Evaluated for 40 HNC cases, the algorithm produced robust IMPT treatment plans within 2-4 h and was validated for proton therapy patient selection.In RayStation 10B, an atlas regression forest-based model was combined with a robust voxel and dose volume histogram-based dose mimicking optimisation algorithm.It was trained and tested on 88 and 25 oropharyngeal cancer patients respectively.Adequate robust target coverage was achieved in 24/25 manual plans and in 23/25 predicted plans with comparable NTCP (van Bruggen et al 2022).
In RayStation 11A, a pipeline was developed using deep learning models with U-net architecture for sequential prediction of nominal and scenario dose from input geometry using deep learning models with U-net architecture.Eriksson and Zhang (2022) trained the models on data from 37 prostate cancer cases and evaluated them on 5.The proposed method achieved robust plans using dose prediction and mimicking.Nominal dose was predicted from nominal geometry, and scenario dose was predicted by deforming the nominal geometry in accordance with the change to scenario geometry.Therefore, scenario dose was predicted for all robustness scenario dose, and plans were created using scenario dose mimicking to ensure a robust deliverable plan (Eriksson and Zhang 2022).
MCO is an approach to find the most suitable treatment plan based on multiple objectives.Often, these objectives conflict with each other.It is crucial to find a balance between them to ensure an optimal treatment outcome.MCO automates the process of finding the optimal treatment plan by considering multiple criteria simultaneously.These criteria typically include target coverage, OAR sparing, dose homogeneity and conformity and criteria related to treatment delivery (Breedveld et al 2019).
Automated MCO solutions for IMPT have been proposed by the work by Taasti et al (2020aTaasti et al ( , 2020b) ) and van de Water et al (2013), van de Water et al (2015).Taasti et al introduced an optimisation method called Expedited Constrained Hierarchical Optimisation (Zarepisheh et al 2019).This method maximises target coverage under dose constraints and then further improves the plan by minimising doses to Organs at Risk (OARs) while still respecting all dose constraints (Zarepisheh et al 2019).Taasti et al validated their approach on three HNC patients in a 13-scenario robust setting (Taasti et al 2020a) and, for four HNC patients, they included beam-angle optimisation without robustness considerations (Taasti et al 2020b).The robust optimisation times were approximately 22 and 36 h for two patients with three beam directions and 36 min for the patient with two beam directions.Van de Water et al proposed an iterative approach called iterative pencil beam resampling (van de Water et al 2013, van De Water et al 2015) (IPBR-MCO) for automated MCO of robust IMPT plans.This approach involves multiple iterations of pencil beam resampling, each followed by an optimisation process.In the last two iterations, they used a wish-list driven MCO of spot intensities as implemented in Erasmus-iCycle.The reported optimisation times for IPBR-MCO were on the order of magnitude of hours.This method has been implemented in Erasmus-iCycle and has been widely used in various studies on robust IMPT optimisation  (Breedveld et al 2012) was used for automated MCO of spot weights, applying for each of the three investigated tumour sites the same optimisation protocol ('wish list') for both approaches.In SISS-MCO, MCO of spot weights was applied once for every patient after a sparsity-induced spot selection (SISS) procedure for pre-selection of the most relevant spots from a large input set of candidate spots.IPBR-MCO comprised several iterations with spot re-sampling, each followed by MCO of the weights of the current spots.Automatically generated plans were not manually fine-tuned.

Methods and materials
2.1.Wish-list based MCO of spot MUs using Erasmus-iCycle In Erasmus-iCycle (Breedveld et al 2012) automated MCO, a single plan is generated for each patient that is Pareto-optimal while obeying all imposed hard constraints.Treatment site-specific wish-lists are used to define the MCO approach for each patient.Besides hard constraints, the wish-list has prioritised cost functions with goal values to define the planning objectives of goals.The final patient dose distribution is produced step-wise by sequential minimisation of objective cost functions from high to low priority, while always obeying all hard constraints.Each step results in an updated plan.The obtained objective value with slight relaxation is added to the problem as a constraint to ensure that following minimisations of lower priority cost functions will not jeopardise already obtained cost function values for the higher priorities.To ensure that final Pareto-optimal plans are also clinically favourable, a wish-list tuning procedure was performed for establishing applied cost functions, priorities, goal values, etc., see the appendix of Heijmen et al (2018).
For optimisation of IMPT plans, wish-lists generally have a constraint to adhere to the clinically applied maximum MU per spot constraint of our institution.An L1 norm was always included as the last objective to again maximally sparsify the solution, i.e. reduce the number of active spots.The wish-lists used in this study also had a minimum target dose constraint, in line with clinical practice.Wish-lists for IMPT plan generation also specify energy layer spacing, lateral spot spacing, and the robustness scenarios considered during optimisation.For the wish-lists used in this study, these settings were chosen in agreement with clinical practice at our institution or with pre-clinical research setting.Beam directions were also specified in the wish-list.

Overview of SISS-MCO
Figure 1 shows a flow diagram describing the generation of plans with SISS-MCO.After dose computations for a large number of candidate pencil beams, plan optimisation comprised three phases: (1) large scale SISS, (2) MCO of MUs for 20% of spots with highest MUs in phase 1, and (3) removal of spots with MUs below the clinically minimum MU (MU min ), and generation of the final Pareto-optimal plan with the Reference Point Method.Sections 2.2.2 to 2.2.4 provide detailed descriptions of these phases.SISS-MCO was embedded in the Erasmus-iCycle environment for automated treatment planning (Breedveld et al 2012).It uses mini-max (Fredriksson et al 2011) scenario-based robust optimisation (section 2.4).

Phase 1: sparsity-induced spot selection
The aim of SISS was to pre-select the most relevant spots from a large set of candidate spots for patients individually, which were then used in the following MCO phase.For each patient, SISS employed a surrogate large-scale bound-constrained optimisation problem where MU was bounded by (0, MU max ).The problem was addressed using a weighted sum cost function that approximated dosimetric aims based on the wish-list and was further expanded with a sparsity-inducing term.By minimising this cost function, MUs of all candidate spots were established and the spots with the highest MUs were considered the most relevant.
Prior to SISS, dose calculations were performed for 20 000 candidate pencil beams (spots) that originated from the beams listed in the wish-list.All spots were positioned in a volume consisting of the composite of the CTVs in all robustness scenarios, enlarged by a treatment site specific margin (section 2.4).Inside this volume, for each beam, energy layers were automatically selected considering the volume and the energy layer spacing as defined in the wish-list.Spots were laterally sampled uniformly in all energy layers, where the maximum lateral spacing within an energy layer was defined in the wish-list.
Next, the wish-list was used to derive the bound-constrained large-scale optimisation problem in equations (1) and (2) that comprises quadratic functions Q i (x) for all target underdose and all OAR hard constraints and for all target overdose objectives, weighted linear function Λ j (x) for all OAR objectives and a sparsity inducement term (L1 norm) λ||x|| 1 .The wish-lists for IMPT did not prescribe other types of cost functions, hence no other penalty functions were used. where indices of quadratic objectives indices of linear objectives number of quadratic objectives number of linear objectives spot intensity vector Q quadratic objectives equation 3 weighting coefficients for the linear objectives linear objective equation 4 coefficient of the vector norm.
The quadratic penalty is defined in equation (3), where the reference dose D ref was set to the value of the constraint or objective in the wish-list and where the value for M was chosen heuristically.An example of quadratic underdose and overdose penalty is given in figure 2. where number of voxels in structure reference dose value indices of voxels dose deposition matrix for quadratic objectives row of dose deposition matrix sgn 1, underdose penalty sgn 1, overdose penalty 10000, for constraints 50, for target dose objectives Equation (4) defines functions Λ j (x) that were used to minimise OAR objectives.In accordance to the wishlist, the mean dose in all scenarios was minimised, or the maximum dose in the worst performing scenario of each iteration.All Λ j (x) were weighted to their priority p j as defined in the wish-list, where the highest priority is 1, using the heuristic in equation (5).
, minimise maximum objective mean , minimise mean objective 4 where A p dose deposition matrix for linear objectives wish list priorities of linear objectives.
The parameter λ, which encourages sparsity, was determined based on the number of linear and quadratic penalties.With increasing number of quadratic penalties, their weights tends to dominate the overall costfunction.Conversely, when there are many linear penalties, λ might diminish their impact.Hence, in such cases, λ is chosen differently.Equation (6) illustrates this heuristic.
The SISS problems were solved using the L-BFGS-B (Byrd et al 1995, Zhu et al 1997, Morales and Nocedal 2011) solver with bound constraints on the MU (equation ( 2)) where MU max was equal to the clinically acceptable maximum MU per spot.The applied number of candidate spots (20 000) and parameters used in SISS-MCO resulted from an extensive grid search on 5 HNC test patients, where the increase in plan quality was weighted against the increased memory (RAM) usage with enhanced numbers of candidate spots.With 20 000 spots, plan quality increase had largely levelled off.The solver was run for either 5000 iterations or until machine floating-point precision convergence, driving the majority of the spot weights to zero.The top 20% highest weighted spots were selected from all spots with a weight greater than 0.1 MU for phase 2.

Phase 2: multi-criteria optimisation (MCO)
In phase 2 of SISS-MCO, the in-house CPU-based Erasmus-iCycle (Breedveld et al 2009(Breedveld et al , 2012) ) algorithm was used for fully automated robust MCO of the weights of the spots selected in phase 1, using a wish-list (see section 2.1) that also includes a constraint on MU max .

Phase 3: reference point method optimisation (RPM)
Phase 2 solutions could include spots with MU MU min < , the clinically minimum allowed MU.In phase 3, these spots were first removed.To satisfy constraints fully, a fast re-optimisation of MUs of the remaining spots was applied using the reference point method (RPM) (Wierzbicki et al 1986, van Haveren et al 2017), where a constraint was added such that MU MU min > .In this approach, the objective values resulting from Phase 2 were used as a reference point, which was projected onto the new Pareto front (due to spot removal).This procedure ensured a Pareto-optimal solution close to the one of Phase 2.

IPBR-MCO
The principle of the applied IPBR-MCO (van de Water et al 2013, van de Water et al 2015) for comparison with the novel SISS-MCO is dimensionality reduction, where a sequence of smaller problems is solved to find a good solution to the original large problem.The approach starts with a relatively small set of randomly chosen candidate pencil beams.After solving the reduced MCO problem, non-contributing spots are removed from the problem, and replaced by a new selection of candidate pencil beams.This is repeated until the solution reaches convergence, i.e. replacing non-contributing spots does not further improve the solution or after a predefined number of iterations (K).Since optimisation times typically grow polynomially with the number of pencil beams, IPBR-MCO reduces optimisation time compared to a single MCO for the full set of candidate pencil beams.
IPBR-MCO has been proposed by van de Water et al (2013), van de Water et al (2015).In this work, IPBR-MCO differs from van de Water et alʼs approach in iterations 1 to K − 2. Here, wish-list objectives, except priority 1, were jointly optimised with priority 2 to expedite the process.The same CPU-based MCO was used as for SISS-MCO.
In this study, 8 iterations with 3000 pencil beams were used for HNC and PC.For CC, the first iteration used 9000 pencil beams and the remaining 8 used 3000 pencil beams, as a lower number of pencil beams for the first iteration often led to an infeasible problem.In total, 24 000 spots were added over all iterations for HNC and PC and, 30 000 for CC.

Robustness and optimisation details
The wish-lists used in this study for both SISS-MCO and IPBR-MCO were previously developed and applied for patient selection for proton therapy, or for research on novel clinical IMPT applications.No modifications were made.The beam angles defined in the wish-list for both SISS-MCO and IPBR-MCO are identical and align with beam angles that have been in clinical use or utilised in pre-clinical research.For serial OARs, there were generally high priority Dmax objectives, while for more parallel structures there was a focus on Dmean.
In all phases of SISS-MCO and IPBR-MCO, plans were generated using scenario-based mini-max optimisation (Fredriksson et al 2011) for all targets and all structures that were flagged in the wish-list as 'robust' (see section 2.5).Nineteen scenarios were used including the nominal scenario, setup errors of 3 or 5 mm depending on the treatment site in positive and negative directions along 3 axes (6 scenarios), and a proton undershoot and overshoot of 3%, see table A1 of the supplementary material and (Korevaar et al 2019).
The available proton energies ranged from 70 to 244 MeV, with corresponding spot widths ranging from 5.9 to 3.3 mm (1σ in air at isocenter), respectively (ProBeam 4.0, Varian a Siemens Healthineers Company, Palo Alto, United States).A range shifter of 57 mm water-equivalent thickness could be inserted during the delivery of a field.Spots and energy layers were sampled from the same grid for fair comparison between SISS-MCO and IPBR-MCO.All dose computations were performed using the Astroid dose-engine (Kooy et al 2010).
All CPU-based computations for SISS-MCO and IPBR-MCO were done on an Intel Xeon Gold 6248 @ 2.5 GHz with 386 GB of available memory.Computations of the bound-constrained weighted-sum cost function used in SISS (section 2.2.2) were performed on a NVIDIA Quadro RTX 6000 with 4608 CUDA cores and 24 GB of available memory.

Patients
Clinical patient data sets were employed in this study.The investigated group of 40 HNC patients consisted of 8 nasopharynx cases, 19 oropharynx, 3 hypopharynx, 3 larynx, 4 tonsil, and 3 tongue base carcinoma patients.The dose was prescribed using a simultaneous integrated boost (SIB) scheme, prescribing 70 Gy RBE to the primary tumour and positive neck levels and 54 Gy RBE to elective neck levels.Robust optimisation was performed for all targets, the brain stem, brain, spinal cord, cochlea, optic nerves, optic chias, eyes and mandible bone, while for the parotids, submandibular glands, oral cavity, oesophagus, larynx, pharyngeal constrictor muscles (superior, middle and inferior) and cricopharyngeus only included the nominal scenario.
The 21 patients with locally advanced PC received a dose of 74 Gy RBE to the prostate and 55 Gy RBE to the seminal vesicles (SV), delivered with a SIB technique.Treatment plans were robustly optimised for all targets and OARs.
For the 19 CC patients, the low-dose ITV had a prescribed dose of 45 Gy RBE and 13/19 patients received SIB treatment with 55 or 57 Gy RBE delivered to positive nodes, according to the EMBRACE II (Pötter et al 2018) protocol.Generated plans had an emphasis on bone marrow sparing (Corbeau et al 2021, Zhou et al 2021).The targets were optimised using 19 scenarios and the OARs only for the nominal scenario.

Plan evaluation and comparison
For each patient, the SISS-MCO plan was compared to the IPBR-MCO plan in terms of plan quality, optimisation time, number of spots and number of energy layers.SISS-MCO optimisation times were measured separately for each of the three phases (section 2.2).For IPBR-MCO, only the total optimisation time was measured.
For all robustly optimised target and OAR dose parameters, plan evaluations and comparisons were performed accordingly.For targets, near-minimum doses, D99%, of the voxel-wise minimum (VWmin) dose distributions, and near-maximum doses, D1%, of voxel-wise maximum (VWmax) doses were used.For robustly optimised OARs, VWmax D1% and/or worst-case scenario mean dose (WC Dmean) were applied, in agreement with the applied wish-lists.For non-robustly optimised OARs, nominal scenario mean dose (N Dmean) or nominal scenario D1% (N D1%) were determined.
For HNC patients, normal tissue complication (NTCP) for grade II and III xerostomia and grade II and III dysphagia were determined according to the Dutch National Protocol for Model-Based Selection for Proton Therapy in HNC (Langendijk et al 2021).For all patients zero baseline toxicity was assumed, as this was unknown due to anonymisation of patient data.
Differences in dosimetric metrics and NTCPs were tested for significance (p < 0.05) using two-sided Wilcoxon signed-rank tests.

Comparisons of optimisation times, number of spots and numbers of energy layers
Average optimisation times, numbers of spots, and numbers of energy layers for the three investigated tumour sites are tabulated in table 1 with their corresponding standard deviations.Overall, total optimisation times were on average a factor 6 shorter with SISS-MCO compared to IPBR-MCO, with a maximum reduction by a factor of 13 for cervix cancer (from 210 to 16 min).For each of the three tumour sites, numbers of energy layers and numbers of spots were similar.
The optimisation times reported in table 1 do not include calculation times of dose deposition matrices for candidate spots as this was not the core of this research.As described in section 2.2.2, in SISS-MCO, 20 000 candidate spots were used for all patients for the three tumour sites.Average dose computation times for 20 000 candidate spots were 118 ± 6, 137 ± 3 and 97 ± 7 min for HNC, PC and CC, respectively.For IPBR more candidate spots were computed (section 2.3).As explained in the Discussion section, there is ample room to significantly reduce dose calculation times.

Dosimetric plan comparison
As demonstrated in figures 3-6, the novel SISS-MCO resulted in dose distributions of similar, or slightly higher, dosimetric quality than the previously published IPBR-MCO.
Figure 3 shows for each of the three investigated tumour sites, nominal SISS-MCO and IPBR-MCO dose distributions for a selected patient.The selected patient had a site-dependent metric difference between SISS-MCO and IPBR-MCO closest to the population mean.The metrics used were NTCP for HNC (figure 3(a)), rectum D99% for the PC case (figure 3(b)) and the mean dose of the whole bone for CC (figure 3(c)).For important structures, population-mean DVHs for doses in nominal scenarios are presented in figure 4.
For HNC, differences between SISS-MCO and IPBR-MCO in plan parameter values and calculated NTCPs are presented by box-and-whisker plots in figures 5(a) and (b), respectively.The positive valued outliers Table 1.Optimisation times, number of energy layers, and number of spots for SISS-MCO and IPBR-MCO.For SISS-MCO, apart from the total optimisation time, the optimisation times for the three consecutive phases (section 2.2 are also reported.Average results are presented for 40 head and neck, 19 cervix, and 21 prostate cancer cases.For cervix cancer, RPM solution was not accepted for both IPBR-MCO and SISS-MCO.observed in figure 5(a) of the VWmax D1% mandible bone, optical nerve L and spinal cord had SISS-MCO values that were far below the clinically acceptable levels in our institution.
Similarly, for PC and CC, differences in plan parameter values are presented in figures 6(a), and 6(b), respectively.In HNC, most differences in OAR plan parameters were in favour of SISS-MCO, and statistically significant (figure 5(a)).This is reflected by statistically significantly lower NTCPs (figure 5(b)).For PC, dose differences were small, with most of the IQR of the dose differences within [−1, 1] Gy (figure 6(a)).The PC targets received on average more dose and the OARs received on average a similar dose.For CC, dose to the targets was very similar in the SISS-MCO and IPBR-MCO plans.The VWmax D1% of the boosted lymph nodes (denoted CTVn) was slightly lower with SISS-MCO, which could be attributed to steeper DVH gradients in SISS-MCO.OAR doses were similar for both methods, with a slightly lower whole bone Dmean for SISS-MCO.

Discussion
Our aim was to develop a computational workflow for fully-automated MCO of high-quality, robust IMPT plans with short optimisation times and applicability to all tumour sites treated with IMPT.The dosimetric quality of the plans generated with the novel SISS-MCO was similar or slightly better than for plans generated with the previously proposed and validated IPBR-MCO for head and neck, prostate and cervix tumours, while the optimisation times were greatly reduced with reduction factors up to 13, depending on the tumour site.For investigated treatment sites, the SISS-MCO optimisation problem was fully defined by the treatment site specific wish-list, without any patient-specific tailoring.The SISS algorithm could be applied for all three tumour sites, without modifications (apart from using a different site-dependent wish-list).
Large improvements in optimisation times compared to IPBR-MCO were obtained notwithstanding the changes in IPBR-MCO compared to the published version, aiming to reduce the optimisation times (section 2.3).As explained in section 2.2.2, selection of relevant, patient-specific spots in SISS-MCO was performed using a bound-constrained surrogate optimisation problem that allowed for fast minimisation with the L-BFGS-B solver, and for each patient, constrained MCO was used only once.On the contrary, in IPBR- MCO, the final spots were iteratively established in multiple rounds of constrained MCO, followed by a last round of constrained MCO to generate the final plan (section 2.3).The observed large differences in optimisation times are largely attributed to these differences in spot selection and application of MCO.
In addition to faster optimisations, another advantage of SISS-MCO over IPBR-MCO in computation time was the lower total number of pencil beams in the optimisations, requiring fewer dose computations, with 20% reductions for HNC and PC and a 50% difference for CC.With the current hardware and implementation, observed mean optimisation/dose calculation times were 43/118, 77/137 and 16/97 min for HNC, PC and CC, respectively, meaning that final plans were generated in 2-4 h.In our clinical setting, this would be totally acceptable, given also the clinical acceptability and high quality of the generated plans.However, there is still a lot of room to further decrease these times as our dose-engine is not optimised for performance.The main limitation is that it is implementation allows usage of only 2 computing threads, resulting in long computation times for large number of spots.In comparison, dose computations for 20 000 spots only take some minutes using our commercial clinical TPS.
In this study, we have systematically compared our proposed SISS-MCO algorithm against the existing IPBR-MCO that used the same wish-lists and sampled from the same spot distributions.Additionally, in another study on HNC, SISS-MCO plans have been dosimetrically compared to clinically delivered plans generated with the RayStation TPS (RaySearch Laboratories, Stockholm, Sweden) (Huiskes et al 2023).For similar target coverage, all OAR doses were lower in the SISS-MCO plans, with sometimes large advantages for SISS-MCO.
While IPBR-MCO had several consecutive rounds of optimisation, each with a relatively small set of pencil beams, SISS-MCO started with an initial pre-selection of relevant pencil beams from large set of candidate pencil beams, followed by only one MCO run.A potential advantage of starting with pre-selection of relevant pencil beams from a large set of candidate pencil beams in SISS-MCO, compared to only working with smaller sets as in IPBR-MCO is that it might prevent spots, that are as a group important for high-quality plan generation, from ending up in different iterations and thereby lose their synergy for generation of the high-quality plans.Often when wish-lists contain many hard constraints, a large set of candidate pencil beams is needed (van de Sande et al 2016).
SISS-MCO has some advantages over data-driven optimisation approaches, as the quality of generated plans does not directly depend on availability of (a large set of) manually generated training plans.Instead, plan quality is fully determined by the optimisation protocol defined by the applied wish-list.As described in appendix of Heijmen et al (2018), wish-list tuning typically uses five example plans when available, but the procedure could also be performed without example plans.If example plans are used, there is no direct relation with the wish-list.For patients used in wish-list tuning, often the plan generated with the final wish-list has higher quality than the original patient plan (Heijmen et al 2018).The wish-list driven approach allows fast adaptation of automated plan generation to changes in clinical protocols; treatment modifications originating from the latest clinical findings can be directly applied to the wish-list to generate automated plans that suit the new treatment insights.The applied wish-list driven approach also has the advantage that final plans are Pareto-optimal.
SISS-MCO was able to automatically generate high-quality plans without the need for patient-specific finetuning afterwards.This may be different for planning approaches that use a weighted sum cost function for final plan generation, if the ambition is to have highest final plan quality for every patient.Although in SISS-MCO a single wish-list was used for all patients with a certain tumour site, this wish-list only described the global optimisation protocol; during each plan generation, patient-specific parameters were automatically established to adapt the plan generation to the specific patient anatomy (Breedveld et al 2012).Breedveld et al (2009) demonstrated that for each patient of a certain tumour site, a weighted sum cost function could result in the same plan as generated with a wish-list, but the weights in the weighted sum cost function had to be patientspecific (they could be derived from Lagrange multipliers obtained in the wish-list based optimisation).In the IMPT plan generation method by Matter et al (2019), importance of OAR constraints were controlled with relative weights.However, it remains unknown whether the weights were patient-specific, as only four patients were considered and no information on the weights was provided.
Sparsity inducement via the L1 norm was previously proposed.Examples include the work by Gu et al (2018Gu et al ( , 2019)), where the L1-norm was used to sparsify selected spots.It is unclear whether these plans guaranteed satisfaction of all constraints exactly, as the dose-fidelity term included non-differentiable penalties.Target dose and maximum dose constraints of the OARs were optimised using Moreau-Yosida regularisation to smooth the non-differentiable penalty.SISS-MCO plans always satisfied constraints, as the MCO phase and the final RPM phase enforce constraints on the solution.Moreover, Gu et al did not apply MU bound-constraints, which could potentially result in clinically unacceptable spot intensities.
This study has some limitations.During the development of SISS-MCO we found that the achievable plan quality increases with the resolution of the candidate spot distribution.More conformal dose and lower OAR burden could be achieved by decreasing the lateral spacing and energy layer spacing by including more candidate spots and allowing for more beam energies.The finally selected number of candidate spots (20 000) was determined from a grid search for HNC, which was a trade-off between the increase in plan quality and the increased memory (RAM) usage for an enhanced number of candidate spots.The number of energy layers was chosen to be consistent with clinical practice.For PC and CC, 20 000 candidate spots were chosen, without investigating optimality of this number for these sites.
The number of iterations and the number of candidate spots per iteration, used for IPBR-MCO, were provided together with the wish-lists.These numbers were previously selected by researchers and medical physicists that created the wish-lists.The numbers were chosen considering target coverage and OAR doses in the robust setting.
Our aim was to develop an automated planning approach that would work for all tumour sites relevant for IMPT.In this study, this was successfully tested for three sites.Possibly, adaptations could be needed for other tumour sites.
Optimality of final solutions after pre-selection of spots in phase 1 of SISS-MCO cannot be guaranteed, but for two arbitrarily selected HNC patients, we applied MCO of phase 2 on the full candidate spot set, and compared the resulting plan, denoted 'Total space-MCO', with the plan generated with SISS-MCO i.e. with upfront SISS spot selection (see figure A1 in the supplementary material).For these patients, quality of the SISS-MCO plan approached the quality of the Full-MCO plan.However, Full-MCO plans required over 30 h of optimisation time using the same CPUs as used for SISS-MCO.We believe that this is an additional indication for adequate spot selection with the (hyper)parameters used in SISS.
All objectives used in this study were linear in the spot weights, in accordance with the provided wish-lists.Quadratic penalties were only used to penalise constraint violations in phase 1 of SISS-MCO.Though other cost functions in IMPT have been tested at our institution, most consistent results were obtained using linear cost functions.However, future research may suggest other cost functions, which could affect the presented results.

Conclusion
A novel approach for fully-automated MCO of robust IMPT plans, designated SISS-MCO, was developed and validated for head and neck, prostate and cervix tumours.Configuring SISS-MCO did not require the availability of large sets of high-quality training plans.Instead, tumour site specific optimisation protocols called wish-lists were used, that can be derived with few of even without training plans.For each patient, plan generation began with an automatic SISS from a large candidate spot set to select the most relevant spots to be used in the final optimisation.
Compared to a previously published and validated automated IMPT MCO approach, plan quality was similar or slightly better for the investigated cases, depending on the tumour site, while optimisation times were greatly reduced, with an average reduction by a factor of 6 (from 287 to 47 min).Including time for dose calculations, final SISS-MCO plans were generated in 2-4 h.

Figure 1 .
Figure 1.Flow diagram showing the three phases of SISS-MCO.RPM = reference point method (see section 2.2.4).

Figure 2 .
Figure 2. Illustration of equation (3) for an OAR with a maximum 40 Gy RBE constraint and a target with dose constraints between 54 and 57 Gy RBE .

Figure 3 .
Figure 3. Side-by-side comparisons of SISS-MCO and IPBR-MCO single slice doses for selected patients (see text for patient selection).(a) Head and neck cancer: the CTV is delineated in red and the elective neck levels in yellow.(b) Prostate cancer: the CTVprostate is delineated in red.(c) Cervical cancer: the ITV is delineated in blue and the positive node is delineated in yellow.

Figure 6 .
Figure 6.The dosimetric difference between SISS-MCO and IPBR-MCO for (a) PC and (b) CC cases.Boxes show interquartile ranges (IQR) of dosimetric differences.In each box, the horizontal bar shows the population median value.Whiskers include all differences, except for outliers, defined as: outside Q1-1.5•IQR and Q3 + 1.5 • IQR.Significant dose differences (p < 0.05) are indicated by asterisks.WC Dmean=Dmean in worst-case scenario, N Dmean = Dmean in nominal scenario, FH = femoral head.
al 2015, van de Water et al 2015, Delaney et al 2018b, 2018a, Kierkels et al 2019, Taasti et al 2020a, 2020b, Celik et al 2021, Xu et al 2021, Eriksson and Zhang 2022, van Bruggen et al 2022) (RP PT) by Varian a Siemens Healthineers company (Palo Alto, CA, USA) (Tol et al 2015, Delaney et al 2018b, 2018a, Celik et al 2021, Xu et al 2021), and methods in RayStation by RaySearch Laboratories (Kierkels et al 2019, Eriksson and Zhang 2022, van Bruggen et al 2022) van de Water et al 2018, Kraan et al 2013, van de Water et al 2016, Arts et al 2017), intra-fraction plan adaptation to accommodate patient anatomy variations (Jagt et al 2017, 2018, Oud et al 2022), and patient selection for IMPT (Kouwenberg et al 2021).Several other published IMPT approaches rely on weighted sum optimisation (Albertini et al 2010, Liu et al 2012, Gu et al 2018, Ungun et al 2019) or weighted least squares optimisation (Fu et al 2022) and therefore require tuning of the cost function weights for each patient for optimal plan quality.In this study, we propose the novel 'SISS-MCO' method for fully automated and fast optimisation of robust IMPT plans, and we compare it with IPBR-MCO (van de Water et al 2013, van de Water et al 2015) for 40 HNC, 21 PC and 19 cervical cancer (CC).Scenario-based robust optimisation of spot weights (monitor units (MU)), and adherence to clinically imposed minimum and maximum MU constraints were applied in both planning approaches.Furthermore, in both SISS-MCO and IPBR-MCO, Erasmus-iCycle